基于IDW插值的測(cè)量數(shù)據(jù)等值線生成方法*

胡添毅，朱江彥，向晉祥，解 遜，周冠男

(長(zhǎng)江航道測(cè)量中心，湖北 武漢 430010)

長(zhǎng)江是貫穿我國(guó)東、中、西部地區(qū)的水路運(yùn)輸大通道。《長(zhǎng)江經(jīng)濟(jì)帶發(fā)展規(guī)劃綱要》的審議通過(guò)，標(biāo)志著長(zhǎng)江航道需要突破當(dāng)前瓶頸、承擔(dān)更多的經(jīng)濟(jì)建設(shè)與上下游聯(lián)通工作，因此需要航道單位提供更好、更便捷的電子航道圖信息服務(wù)。電子航道圖的數(shù)據(jù)基礎(chǔ)為內(nèi)業(yè)制圖成果，等值線繪制及檢查作為內(nèi)業(yè)制圖最重要的一環(huán)，需要耗費(fèi)5名內(nèi)業(yè)人員4～5工作日，其效率難以滿足快速成圖、為電子航道圖及時(shí)提供數(shù)據(jù)的需求，亟待提升制作工藝。

目前，在航道測(cè)繪中，繪制等值線的基本流程為：1)將測(cè)量數(shù)據(jù)導(dǎo)入CASS軟件；2)框選數(shù)據(jù)生成三角網(wǎng)；3)利用三角網(wǎng)計(jì)算生成等值線；4)將該等值線導(dǎo)入清華山維軟件；5)在清華山維軟件內(nèi)手動(dòng)調(diào)整等值線；6)將調(diào)整的等值線導(dǎo)出為CAD格式。

該繪制流程存在兩點(diǎn)不足：1)現(xiàn)有等值線生成流程步驟較多，需要制圖人員對(duì)CASS、清華山維軟件有充分了解，并掌握各項(xiàng)繪圖參數(shù)，人員培訓(xùn)成本較大；2)手動(dòng)調(diào)整等值線工作量過(guò)大，并且在調(diào)整過(guò)程中可能出現(xiàn)錯(cuò)繪、漏繪等情況。為了提升等值線繪制的效率與準(zhǔn)確性，等值線自動(dòng)生成方法被廣泛研究。

大部分的等值線自動(dòng)生成方法[1-3]基于不規(guī)則三角網(wǎng)(triangular irregular network,TIN)模型或其改進(jìn)模型。王海濤[4]提出一種基于薄板樣條函數(shù)(TPS)模型構(gòu)建的等深線生產(chǎn)方法，并通過(guò)實(shí)際數(shù)據(jù)驗(yàn)證其有效性。王學(xué)潮[5]將測(cè)深點(diǎn)網(wǎng)格化后，設(shè)計(jì)了通過(guò)等值點(diǎn)追蹤生成等值線、并予以平滑的算法。

本文設(shè)計(jì)了基于IDW(inverse distance weighted)插值的等值線生成方法，并使用空間查詢方法對(duì)等值線進(jìn)行線形優(yōu)化和拓?fù)溴e(cuò)誤修正，最后基于生成和優(yōu)化方法開(kāi)發(fā)了等值線生成平臺(tái)。

1 基本方法

在航道測(cè)繪作業(yè)中，測(cè)量得出的水深值或高程值均是在測(cè)量區(qū)域內(nèi)處處都有定義的地理特征值，即空間連續(xù)數(shù)據(jù)(spatial continuous data)。為表示完整的河床地形地貌，在源數(shù)據(jù)中只記錄有限個(gè)樣點(diǎn)值，樣點(diǎn)以外各點(diǎn)的值通過(guò)插值方法計(jì)算。

同時(shí)，進(jìn)行一次測(cè)量可能會(huì)產(chǎn)生數(shù)萬(wàn)個(gè)測(cè)量點(diǎn)，等值線的生成中需要對(duì)這些測(cè)量點(diǎn)進(jìn)行存儲(chǔ)、空間查詢及計(jì)算，以提升等值線計(jì)算的效率。本文使用MySQL數(shù)據(jù)庫(kù)進(jìn)行測(cè)量點(diǎn)和等值線數(shù)據(jù)的存儲(chǔ)、空間查詢與計(jì)算。

1.1 插值方法

空間插值方法是基于已知的地理特征點(diǎn)數(shù)據(jù)，推求同一區(qū)域內(nèi)的未知點(diǎn)的數(shù)據(jù)。

目前，國(guó)內(nèi)外常見(jiàn)的空間插值方法有IDW插值[6-8]、不規(guī)則三角網(wǎng)方法、趨勢(shì)面(trend)分析法、普通克里格(ordinary Kriging)插值、簡(jiǎn)單克里格(simple Kriging)插值、泛克里格(universal Kriging)插值等。

IDW插值全稱為反距離加權(quán)插值，它以插值點(diǎn)與樣本點(diǎn)之間的距離為權(quán)重進(jìn)行加權(quán)平均，離插值點(diǎn)越近的樣本點(diǎn)賦予的權(quán)重越大。

IDW插值的公式為：

(1)

(2)

因?yàn)樵诤降罍y(cè)繪后，測(cè)量點(diǎn)經(jīng)過(guò)點(diǎn)位校正在地圖上分布均勻，密度也足以在計(jì)算中反映河床局部表面變化，本文選取IDW插值方法進(jìn)行等值線生成的計(jì)算。

1.2 MySQL的空間數(shù)據(jù)庫(kù)應(yīng)用

MySQL是一種關(guān)系型數(shù)據(jù)庫(kù)管理系統(tǒng)[9-10]，MySQL5.7版本遵循OpenGIS聯(lián)盟(open geospatial consortium，OGC)的幾何數(shù)據(jù)模型規(guī)則，對(duì)空間數(shù)據(jù)庫(kù)進(jìn)行支持，并支持空間索引[11]，可以提高對(duì)數(shù)據(jù)集空間檢索的效率。

空間查詢利用空間索引機(jī)制，從數(shù)據(jù)庫(kù)中找出符合該條件的空間數(shù)據(jù)，MySQL支持MBRContains()、ST_Contains()、ST_Disjoint()等空間拓?fù)浜瘮?shù)，并可使用ST_PointN(A,n)、ST_Area()、ST_Buffer()等空間數(shù)據(jù)相關(guān)方法。采用MySQL數(shù)據(jù)庫(kù)可以有效加快等值線生成平臺(tái)的計(jì)算速度，進(jìn)而提升內(nèi)業(yè)制圖人員的工作效率。

本次研究所用的空間拓?fù)浜瘮?shù)的作用為：MBRContains(A,B),判斷幾何形狀A(yù)的最小邊界矩形是否包含幾何形狀B的最小邊界矩形；ST_Contains(A,B),判斷幾何形狀A(yù)是否完全包含幾何形狀B(當(dāng)且僅當(dāng)B中的所有點(diǎn)都位于A中，并且至少有一個(gè)B中的點(diǎn)位于A的內(nèi)部，為完全包含)；ST_Intersection(A,B),計(jì)算幾何形狀A(yù)和B的相交或重疊區(qū)域；ST_Disjoint(A,B),判斷幾何形狀A(yù)與幾何形狀B是否有交集；ST_PointN(A,n),返回折線A的第n個(gè)點(diǎn)；ST_Area(A)：計(jì)算多邊形A的面積；ST_Buffer(A,x)：計(jì)算生成幾何形狀A(yù)、長(zhǎng)度為x的多邊形緩沖區(qū)。

2 基于ArcGIS的等值線生成流程

本文通過(guò)IDW插值算法生成測(cè)量區(qū)域的等值線圖，再?gòu)牡戎稻€圖中提取等值線，根據(jù)測(cè)量區(qū)域的邊界進(jìn)行裁剪。等值線具體的生成流程見(jiàn)圖1。

圖1 等值線生成流程

將水深測(cè)量原始數(shù)據(jù)進(jìn)行預(yù)處理后，對(duì)測(cè)深點(diǎn)進(jìn)行IDW插值計(jì)算，得到測(cè)量區(qū)域的等值面圖，并按照高程間距1 m提取等值線。

隨后根據(jù)測(cè)深點(diǎn)計(jì)算測(cè)深區(qū)域邊界，計(jì)算流程為：1)對(duì)每個(gè)測(cè)深點(diǎn)都生成緩沖區(qū)；2)求所有緩沖區(qū)的并集；3)將求得的并集生成柵格圖層；4)生成柵格圖層的線狀外接多邊形；5)對(duì)線狀外界多邊形進(jìn)行轉(zhuǎn)換，生成面狀外接多邊形；6)利用測(cè)深區(qū)域邊界與等值線的拓?fù)潢P(guān)系裁剪等值線，最終生成裁剪后的等值線。

根據(jù)上述等值線生成流程開(kāi)發(fā)等值線生成平臺(tái)。在Visual Studio環(huán)境下，MySQL作為數(shù)據(jù)庫(kù)，以ArcGIS Engine工具進(jìn)行空間數(shù)據(jù)顯示及處理，采用C#語(yǔ)言搭建框架進(jìn)行開(kāi)發(fā)。

本文的測(cè)量原始數(shù)據(jù)來(lái)源于武漢白沙洲橋區(qū)2017年9月實(shí)地測(cè)量數(shù)據(jù)，測(cè)區(qū)面積約為7.5萬(wàn)m2,采用多波束測(cè)深儀進(jìn)行河床掃測(cè)，采集約3.3萬(wàn)個(gè)測(cè)深點(diǎn)，并經(jīng)過(guò)數(shù)據(jù)預(yù)處理。等值線生成平臺(tái)和導(dǎo)入平臺(tái)的測(cè)量原始數(shù)據(jù)顯示見(jiàn)圖2，利用上述流程在平臺(tái)中計(jì)算生成的等值線見(jiàn)圖3。

圖2 等值線生成平臺(tái)及測(cè)量原始數(shù)據(jù)

圖3 生成的等值線

3 等值線優(yōu)化方法

3.1 等值線優(yōu)化需求

為確保等值線可以準(zhǔn)確地反映河道地理信息，需要查驗(yàn)等值線與測(cè)深點(diǎn)的對(duì)應(yīng)關(guān)系，以及等值線自身與等值線之間的拓?fù)潢P(guān)系。因?yàn)樗钤紨?shù)據(jù)的水深或高程值保留一位小數(shù)，所以會(huì)出現(xiàn)一批值為“x.0”(x=…,-1,0,1,…)的測(cè)深點(diǎn)。此類測(cè)深點(diǎn)會(huì)使生成的等值線與測(cè)深點(diǎn)對(duì)應(yīng)關(guān)系出現(xiàn)錯(cuò)誤，并且在拓?fù)潢P(guān)系、線形上出現(xiàn)錯(cuò)誤，如圖4所示，細(xì)線為計(jì)算生成的等值線，粗線為等值線的理想走向，數(shù)字為原始數(shù)據(jù)的測(cè)深值注記。

圖4 等值線錯(cuò)誤狀況示例

同時(shí)，生成的等值線在線形美觀方面相比手工繪制的曲線有較大的差距。如圖5所示，需要對(duì)曲線部分節(jié)點(diǎn)的點(diǎn)位進(jìn)行校正，并對(duì)曲線進(jìn)行平滑。

圖5 等值線線形優(yōu)化需求

3.2 特征點(diǎn)與等值線對(duì)應(yīng)關(guān)系

此處的特征點(diǎn)是值為“x.0”(x=…,-1,0,1,…)的測(cè)深點(diǎn)，為確保等值線與特征點(diǎn)對(duì)應(yīng)關(guān)系正確，可以考慮對(duì)特征點(diǎn)的部分特征值進(jìn)行修改，以對(duì)IDW插值結(jié)果進(jìn)行修正。

因此，本文根據(jù)特征點(diǎn)與周圍8個(gè)點(diǎn)的數(shù)值對(duì)比判斷是否需要修改特征值，比對(duì)步驟為：第一步，利用MySQL的MBRContains()和ST_Buffer()函數(shù)生成該點(diǎn)的緩沖區(qū)，通過(guò)空間拓?fù)洳樵儷@取最鄰近的8個(gè)點(diǎn)；第二步，根據(jù)點(diǎn)的坐標(biāo)進(jìn)行排序，最終得到9條點(diǎn)位的記錄，特征點(diǎn)與臨近8個(gè)點(diǎn)的排列順序見(jiàn)圖6；第三步，根據(jù)設(shè)定規(guī)則判斷是否修改特征值，部分規(guī)則見(jiàn)圖7，左邊為調(diào)整之前的測(cè)深點(diǎn)與IDW方法生成的等值線，右邊為數(shù)值調(diào)整之后的測(cè)深點(diǎn)與等值線?？梢钥吹?，經(jīng)過(guò)數(shù)值調(diào)整后的測(cè)深點(diǎn)生成的等值線可以與特征點(diǎn)有更好的對(duì)應(yīng)關(guān)系。

圖6 特征點(diǎn)與鄰近點(diǎn)點(diǎn)位排序

圖7 修改特征值部分規(guī)則

3.3 等值線優(yōu)化及平滑

對(duì)特征點(diǎn)的值進(jìn)行修改并通過(guò)IDW生成等值線后，需要對(duì)等值線的線型進(jìn)行修改，以使部分曲線準(zhǔn)確地反映測(cè)深點(diǎn)的高程特征，并使等值線更加美觀。

對(duì)于節(jié)點(diǎn)數(shù)超過(guò)2的等值線，優(yōu)化流程為：1)遍歷等值線中不是起點(diǎn)和終點(diǎn)的節(jié)點(diǎn)；2)空間查詢判斷節(jié)點(diǎn)是否落于測(cè)深點(diǎn)上；3)如果否，則查找該節(jié)點(diǎn)的上一節(jié)點(diǎn)和下一節(jié)點(diǎn)；4)計(jì)算該節(jié)點(diǎn)與上下節(jié)點(diǎn)連線的長(zhǎng)度；5)根據(jù)余弦定理判斷該節(jié)點(diǎn)與上下節(jié)點(diǎn)連線的夾角；6)根據(jù)矢量叉積判斷該節(jié)點(diǎn)是否為拐點(diǎn);7)如果該節(jié)點(diǎn)非拐點(diǎn)，上下節(jié)點(diǎn)連線的長(zhǎng)度不超過(guò)閾值,并且角度小于60°或大于160°，刪除該節(jié)點(diǎn)。

曲線優(yōu)化后，利用ArcGIS Engine的Smooth()函數(shù)，對(duì)曲線進(jìn)行貝塞爾平滑處理，并且利用Densify()函數(shù)增加平滑后曲線的節(jié)點(diǎn)密度。等值線優(yōu)化平滑效果見(jiàn)圖8，粗線為優(yōu)化前的曲線，細(xì)線為優(yōu)化后的曲線，數(shù)字為原始數(shù)據(jù)的測(cè)深值注記。

圖8 等值線優(yōu)化平滑效果

3.4 拓?fù)溴e(cuò)誤修正

等值線進(jìn)行平滑處理后，部分曲線的走向會(huì)發(fā)生改變，在等值線分布密集的區(qū)域可能出現(xiàn)相交、重疊等情況，等值線與測(cè)深點(diǎn)的對(duì)應(yīng)關(guān)系可能發(fā)生改變。為了保證等值線圖的準(zhǔn)確性，需要對(duì)等值線進(jìn)行拓?fù)錂z查與修正，檢查的圖層為經(jīng)過(guò)優(yōu)化平滑的等值線圖層與修改部分特征點(diǎn)特征值的測(cè)深點(diǎn)圖層。

拓?fù)錂z查的內(nèi)容為：

1)檢查線與線之間相互重疊情況，通過(guò)MySQL的ST_disjoint()函數(shù)判斷。

2)檢查線與線之間相交情況，通過(guò)MySQL的ST_disjoint()函數(shù)判斷。

3)檢查線要素自身重疊或相交的情況，通過(guò)曲線的節(jié)點(diǎn)遍歷查詢判斷是否有重復(fù)節(jié)點(diǎn)，若有則說(shuō)明自身有重疊或相交的情況。

4)檢查線要素被特征值不是“x.0”(x=…,-1,0,1,…)的測(cè)深點(diǎn)覆蓋的情況，通過(guò)MySQL的ST_contains()函數(shù)判斷。

5)檢查特征值是“x.0”(x=…,-1,0,1,…)的測(cè)深點(diǎn)上沒(méi)有線要素的情況，通過(guò)MySQL的ST_contains()函數(shù)判斷。

在這里，需要定義元組TopoError(geometryA,geometryB,isclosedA,isclosedB)，其中g(shù)eometryA和geometryB用來(lái)記錄發(fā)生拓?fù)溴e(cuò)誤的形狀的id；isclosedA和isclosedB用來(lái)記錄形狀在等值線圖的范圍內(nèi)是否閉合，0為閉合，1為不閉合，2表示形狀不是曲線。在進(jìn)行拓?fù)溴e(cuò)誤的檢查后，將產(chǎn)生拓?fù)溴e(cuò)誤的形狀記錄在TopoError的數(shù)組中，并進(jìn)行拓?fù)溴e(cuò)誤的修正，并生成拓?fù)湫拚蟮牡戎稻€圖層，拓?fù)湫拚鞒桃?jiàn)圖9。

圖9 等值線拓?fù)湫拚鞒?/p>

等值線拓?fù)湫拚男Ч?jiàn)圖10，細(xì)線為修正前的曲線，粗線為修正后的曲線，數(shù)字為原始數(shù)據(jù)的測(cè)深值注記。

圖10 等值線拓?fù)湫拚Ч?/p>

3.5 線形渲染及輸出

等值線經(jīng)過(guò)拓?fù)錂z查修正，已與測(cè)深點(diǎn)的對(duì)應(yīng)關(guān)系準(zhǔn)確，走向合理，可以輸出為CAD圖層。輸出為CAD前，針對(duì)等值線的特征值，進(jìn)行首曲線與計(jì)曲線的區(qū)分，并根據(jù)內(nèi)業(yè)制圖規(guī)范賦予不同的顏色和線型。最終利用ArcGIS Engine的ConversionTools工具將圖層轉(zhuǎn)換為CAD文件，得到等值線生成的最終成果。經(jīng)質(zhì)檢人員檢查，最終成果不存在拓?fù)溴e(cuò)誤。

使用等值線生成平臺(tái)進(jìn)行等值線的生成，其流程簡(jiǎn)化為“導(dǎo)入測(cè)量數(shù)據(jù)——點(diǎn)擊‘等值線生成’按鈕計(jì)算”兩步。大幅降低了內(nèi)業(yè)人員的操作難度與工作強(qiáng)度，提升了內(nèi)業(yè)成圖效率，可以及時(shí)地為電子航道圖提供數(shù)據(jù)，滿足電子航道圖實(shí)時(shí)性需求。

4 結(jié)語(yǔ)

1)利用IDW插值法進(jìn)行計(jì)算，得到自動(dòng)生成的等值線；其次通過(guò)空間查詢方法對(duì)等值線進(jìn)行線型優(yōu)化與拓?fù)湫拚幚?；最后根?jù)等值線生成優(yōu)化方法開(kāi)發(fā)等值線生成平臺(tái)，驗(yàn)證整個(gè)流程的高效性與準(zhǔn)確性。本文提出的方法對(duì)等值線自動(dòng)化生成研究有很大的借鑒意義。

2)本文提出的方法存在不足：自動(dòng)生成的等值線在美觀方面有所欠缺、沒(méi)考慮分布離散的陸域地形測(cè)量點(diǎn)等問(wèn)題。這些問(wèn)題有待日后深入研究。