基于IDW插值的測量數據等值線生成方法*

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

(長江航道測量中心，湖北 武漢 430010)

長江是貫穿我國東、中、西部地區的水路運輸大通道。《長江經濟帶發展規劃綱要》的審議通過，標志著長江航道需要突破當前瓶頸、承擔更多的經濟建設與上下游聯通工作，因此需要航道單位提供更好、更便捷的電子航道圖信息服務。電子航道圖的數據基礎為內業制圖成果，等值線繪制及檢查作為內業制圖最重要的一環，需要耗費5名內業人員4～5工作日，其效率難以滿足快速成圖、為電子航道圖及時提供數據的需求，亟待提升制作工藝。

目前，在航道測繪中，繪制等值線的基本流程為：1)將測量數據導入CASS軟件；2)框選數據生成三角網；3)利用三角網計算生成等值線；4)將該等值線導入清華山維軟件；5)在清華山維軟件內手動調整等值線；6)將調整的等值線導出為CAD格式。

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

大部分的等值線自動生成方法[1-3]基于不規則三角網(triangular irregular network,TIN)模型或其改進模型。王海濤[4]提出一種基于薄板樣條函數(TPS)模型構建的等深線生產方法，并通過實際數據驗證其有效性。王學潮[5]將測深點網格化后，設計了通過等值點追蹤生成等值線、并予以平滑的算法。

本文設計了基于IDW(inverse distance weighted)插值的等值線生成方法，并使用空間查詢方法對等值線進行線形優化和拓撲錯誤修正，最后基于生成和優化方法開發了等值線生成平臺。

1 基本方法

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

同時，進行一次測量可能會產生數萬個測量點，等值線的生成中需要對這些測量點進行存儲、空間查詢及計算，以提升等值線計算的效率。本文使用MySQL數據庫進行測量點和等值線數據的存儲、空間查詢與計算。

1.1 插值方法

空間插值方法是基于已知的地理特征點數據，推求同一區域內的未知點的數據。

目前，國內外常見的空間插值方法有IDW插值[6-8]、不規則三角網方法、趨勢面(trend)分析法、普通克里格(ordinary Kriging)插值、簡單克里格(simple Kriging)插值、泛克里格(universal Kriging)插值等。

IDW插值全稱為反距離加權插值，它以插值點與樣本點之間的距離為權重進行加權平均，離插值點越近的樣本點賦予的權重越大。

IDW插值的公式為：

(1)

(2)

因為在航道測繪后，測量點經過點位校正在地圖上分布均勻，密度也足以在計算中反映河床局部表面變化，本文選取IDW插值方法進行等值線生成的計算。

1.2 MySQL的空間數據庫應用

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

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

本次研究所用的空間拓撲函數的作用為：MBRContains(A,B),判斷幾何形狀A的最小邊界矩形是否包含幾何形狀B的最小邊界矩形；ST_Contains(A,B),判斷幾何形狀A是否完全包含幾何形狀B(當且僅當B中的所有點都位于A中，并且至少有一個B中的點位于A的內部，為完全包含)；ST_Intersection(A,B),計算幾何形狀A和B的相交或重疊區域；ST_Disjoint(A,B),判斷幾何形狀A與幾何形狀B是否有交集；ST_PointN(A,n),返回折線A的第n個點；ST_Area(A)：計算多邊形A的面積；ST_Buffer(A,x)：計算生成幾何形狀A、長度為x的多邊形緩沖區。

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

本文通過IDW插值算法生成測量區域的等值線圖，再從等值線圖中提取等值線，根據測量區域的邊界進行裁剪。等值線具體的生成流程見圖1。

圖1 等值線生成流程

將水深測量原始數據進行預處理后，對測深點進行IDW插值計算，得到測量區域的等值面圖，并按照高程間距1 m提取等值線。

隨后根據測深點計算測深區域邊界，計算流程為：1)對每個測深點都生成緩沖區；2)求所有緩沖區的并集；3)將求得的并集生成柵格圖層；4)生成柵格圖層的線狀外接多邊形；5)對線狀外界多邊形進行轉換，生成面狀外接多邊形；6)利用測深區域邊界與等值線的拓撲關系裁剪等值線，最終生成裁剪后的等值線。

根據上述等值線生成流程開發等值線生成平臺。在Visual Studio環境下，MySQL作為數據庫，以ArcGIS Engine工具進行空間數據顯示及處理，采用C#語言搭建框架進行開發。

本文的測量原始數據來源于武漢白沙洲橋區2017年9月實地測量數據，測區面積約為7.5萬m2,采用多波束測深儀進行河床掃測，采集約3.3萬個測深點，并經過數據預處理。等值線生成平臺和導入平臺的測量原始數據顯示見圖2，利用上述流程在平臺中計算生成的等值線見圖3。

圖2 等值線生成平臺及測量原始數據

圖3 生成的等值線

3 等值線優化方法

3.1 等值線優化需求

為確保等值線可以準確地反映河道地理信息，需要查驗等值線與測深點的對應關系，以及等值線自身與等值線之間的拓撲關系。因為水深原始數據的水深或高程值保留一位小數，所以會出現一批值為“x.0”(x=…,-1,0,1,…)的測深點。此類測深點會使生成的等值線與測深點對應關系出現錯誤，并且在拓撲關系、線形上出現錯誤，如圖4所示，細線為計算生成的等值線，粗線為等值線的理想走向，數字為原始數據的測深值注記。

圖4 等值線錯誤狀況示例

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

圖5 等值線線形優化需求

3.2 特征點與等值線對應關系

此處的特征點是值為“x.0”(x=…,-1,0,1,…)的測深點，為確保等值線與特征點對應關系正確，可以考慮對特征點的部分特征值進行修改，以對IDW插值結果進行修正。

因此，本文根據特征點與周圍8個點的數值對比判斷是否需要修改特征值，比對步驟為：第一步，利用MySQL的MBRContains()和ST_Buffer()函數生成該點的緩沖區，通過空間拓撲查詢獲取最鄰近的8個點；第二步，根據點的坐標進行排序，最終得到9條點位的記錄，特征點與臨近8個點的排列順序見圖6；第三步，根據設定規則判斷是否修改特征值，部分規則見圖7，左邊為調整之前的測深點與IDW方法生成的等值線，右邊為數值調整之后的測深點與等值線。可以看到，經過數值調整后的測深點生成的等值線可以與特征點有更好的對應關系。

圖6 特征點與鄰近點點位排序

圖7 修改特征值部分規則

3.3 等值線優化及平滑

對特征點的值進行修改并通過IDW生成等值線后，需要對等值線的線型進行修改，以使部分曲線準確地反映測深點的高程特征，并使等值線更加美觀。

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

曲線優化后，利用ArcGIS Engine的Smooth()函數，對曲線進行貝塞爾平滑處理，并且利用Densify()函數增加平滑后曲線的節點密度。等值線優化平滑效果見圖8，粗線為優化前的曲線，細線為優化后的曲線，數字為原始數據的測深值注記。

圖8 等值線優化平滑效果

3.4 拓撲錯誤修正

等值線進行平滑處理后，部分曲線的走向會發生改變，在等值線分布密集的區域可能出現相交、重疊等情況，等值線與測深點的對應關系可能發生改變。為了保證等值線圖的準確性，需要對等值線進行拓撲檢查與修正，檢查的圖層為經過優化平滑的等值線圖層與修改部分特征點特征值的測深點圖層。

拓撲檢查的內容為：

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

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

3)檢查線要素自身重疊或相交的情況，通過曲線的節點遍歷查詢判斷是否有重復節點，若有則說明自身有重疊或相交的情況。

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

5)檢查特征值是“x.0”(x=…,-1,0,1,…)的測深點上沒有線要素的情況，通過MySQL的ST_contains()函數判斷。

在這里，需要定義元組TopoError(geometryA,geometryB,isclosedA,isclosedB)，其中geometryA和geometryB用來記錄發生拓撲錯誤的形狀的id；isclosedA和isclosedB用來記錄形狀在等值線圖的范圍內是否閉合，0為閉合，1為不閉合，2表示形狀不是曲線。在進行拓撲錯誤的檢查后，將產生拓撲錯誤的形狀記錄在TopoError的數組中，并進行拓撲錯誤的修正，并生成拓撲修正后的等值線圖層，拓撲修正流程見圖9。

圖9 等值線拓撲修正流程

等值線拓撲修正的效果見圖10，細線為修正前的曲線，粗線為修正后的曲線，數字為原始數據的測深值注記。

圖10 等值線拓撲修正效果

3.5 線形渲染及輸出

等值線經過拓撲檢查修正，已與測深點的對應關系準確，走向合理，可以輸出為CAD圖層。輸出為CAD前，針對等值線的特征值，進行首曲線與計曲線的區分，并根據內業制圖規范賦予不同的顏色和線型。最終利用ArcGIS Engine的ConversionTools工具將圖層轉換為CAD文件，得到等值線生成的最終成果。經質檢人員檢查，最終成果不存在拓撲錯誤。

使用等值線生成平臺進行等值線的生成，其流程簡化為“導入測量數據——點擊‘等值線生成’按鈕計算”兩步。大幅降低了內業人員的操作難度與工作強度，提升了內業成圖效率，可以及時地為電子航道圖提供數據，滿足電子航道圖實時性需求。

4 結語

1)利用IDW插值法進行計算，得到自動生成的等值線；其次通過空間查詢方法對等值線進行線型優化與拓撲修正處理；最后根據等值線生成優化方法開發等值線生成平臺，驗證整個流程的高效性與準確性。本文提出的方法對等值線自動化生成研究有很大的借鑒意義。

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