一種基于網(wǎng)格延展法的等值線生成算法

陳 勝 ，黃詩峰 ，臧文斌

（1. 中國水利水電科學(xué)研究院，北京 100038；2. 水利部防洪抗旱減災(zāi)工程技術(shù)研究中心，北京 100038）

0 引言

等值線在水利信息化系統(tǒng)和科學(xué)研究中有著廣泛的應(yīng)用，是數(shù)據(jù)可視化的重要手段之一。通過雨量等值線可以表達(dá)降雨量的空間分布特性[1]；通過地下水埋深等值面可以實時掌握地下水動態(tài)，為水資源開發(fā)、水資源量管理提供依據(jù)[2-3]；通過旱情等值線分布，可研究和監(jiān)測旱情的時空分布與演變[4]。等值線是將有限的觀測點數(shù)據(jù)，通過專門的等值線算法生成圖形數(shù)據(jù)，用于表達(dá)數(shù)據(jù)在空間的分布特性。

等值線生成算法從生成方法上分主要有基于三角網(wǎng)[5]和規(guī)則網(wǎng)格[6]2 種生成算法?；谌蔷W(wǎng)的算法計算效率高、精度高，但算法實現(xiàn)過程比較復(fù)雜，當(dāng)數(shù)據(jù)點分布不規(guī)律時容易出現(xiàn)異常，如出現(xiàn)等值線交叉、分叉等現(xiàn)象；基于規(guī)則網(wǎng)格的算法計算方法簡單，但是精度不夠高。已有的研究中二者在等值線追蹤方面都需要考慮等值線交叉、連通域等問題。在生產(chǎn)應(yīng)用中基于三角網(wǎng)進(jìn)行等值線計算的算法，因為算法比較復(fù)雜，往往需要借助第三方軟件和組件實現(xiàn)，如 ArcGIS，MapX，Surface等[7-9]工具軟件中包含等值線組件，處理過程不但需要考慮如何在算法組件和應(yīng)用系統(tǒng)之間進(jìn)行集成，而且在實際應(yīng)用時常常出現(xiàn)因為數(shù)據(jù)點分布原因表現(xiàn)出的算法健壯性不夠的情況，比如出現(xiàn)等值線交叉，等值面局部區(qū)域缺失等異常現(xiàn)象。等值線的生成算法近幾年研究較少，等值線的生產(chǎn)應(yīng)用及相關(guān)問題的研究解決近幾年相對較多。

本研究立足于水利信息化生產(chǎn)應(yīng)用，考慮降雨和水資源分布 2 類等值線，主要用于表達(dá)和分析降雨量、水資源的宏觀分布，并不需要很高的精度，同時數(shù)據(jù)點空間分布也不規(guī)律。因此研究一種基于網(wǎng)格延展法的等值線算法，目標(biāo)是使等值線算法簡單、健壯，易于實現(xiàn)，并能夠滿足生產(chǎn)應(yīng)用的需求。

1 算法原理

網(wǎng)格延展法是一種基于規(guī)則網(wǎng)格的等值線生成新方法，該算法沒有采用已有三角網(wǎng)或規(guī)則網(wǎng)格等值線生成研究中的基于線的追蹤方法，而是利用研究區(qū)域內(nèi)同一等值面在網(wǎng)格空間上是連續(xù)分布這一特性進(jìn)行研究的，從研究區(qū)域內(nèi)部的一個網(wǎng)格點向四周進(jìn)行遞歸搜索，可以延展到研究區(qū)域內(nèi)的其他所有連通且取值范圍相同的網(wǎng)格點，從而形成一個等值面。采用網(wǎng)格延展法生成等值線的基本流程如圖1 所示。

圖1 基于網(wǎng)格延展法的等值線生成過程

1.1 研究區(qū)域網(wǎng)格構(gòu)建

基于網(wǎng)格延展法生成等值線，首先構(gòu)建用于等值線計算的規(guī)則網(wǎng)格，網(wǎng)格的范圍應(yīng)覆蓋整個研究區(qū)域，研究區(qū)域往往是某一區(qū)域的行政邊界或某一流域的邊界。整體網(wǎng)格的外包絡(luò)矩形記錄有地理坐標(biāo)信息，便于最終生成的等值線從網(wǎng)格坐標(biāo)轉(zhuǎn)換為地理坐標(biāo)。網(wǎng)格的精度根據(jù)實際應(yīng)用需要確定，網(wǎng)格數(shù)目越大生成的等值線越細(xì)致，但需要耗費的計算量也越大。實際應(yīng)用研究表明，對于雨量等值線來說，一般生成規(guī)模為 1 000×1 000 的網(wǎng)格即可滿足生產(chǎn)應(yīng)用需求。研究區(qū)域網(wǎng)格在算法程序中表示為一個二維數(shù)組，為提高算法效率可以創(chuàng)建整數(shù)類型的二維數(shù)組存儲網(wǎng)格信息，這樣可避免浮點類型的計算，減少運行內(nèi)存的開銷，等值線的計算效率會有明顯提升?；诰W(wǎng)格延展法的等值線算法程序中需要構(gòu)建多個同等大小的輔助網(wǎng)格，用來保存邊界標(biāo)記信息、中間計算成果等。

1.2 等值線邊界確定

等值線邊界確定是為了將等值線的計算限定在一定的邊界范圍內(nèi)進(jìn)行，不同于已有研究中先生成等值線再根據(jù)邊界裁剪的確定方法。

多數(shù)已有的等值線生成算法研究中，一般是先生成 1 個較大區(qū)域（如等值線外包絡(luò)矩形）的等值線，然后根據(jù)某一區(qū)域的行政邊界或某一流域的邊界進(jìn)行裁剪，得到指定邊界范圍內(nèi)的等值線?；诰W(wǎng)格延展法進(jìn)行等值線計算時，將等值線生成計算限定在研究區(qū)域邊界范圍內(nèi)進(jìn)行，這樣可避免復(fù)雜的等值線裁剪計算過程，也省去無效區(qū)域的計算工作量。邊界數(shù)據(jù)一般為 1 組矢量數(shù)據(jù)，如 ArcGIS 里的面狀圖層，是由折線（折線由點序列構(gòu)成）構(gòu)成的不規(guī)則多邊形，截取的其中一部分如圖2 中的線段所示。此時需要將這些折線對應(yīng)到網(wǎng)格上，并形成連續(xù)的封閉區(qū)域作為等值線計算的范圍。

圖2 等值線邊界柵格化

本研究通過對研究區(qū)域邊界數(shù)據(jù)進(jìn)行柵格化處理和標(biāo)記，從而在構(gòu)建的網(wǎng)格上形成連續(xù)的封閉區(qū)域。算法步驟如下：

1）步驟 1。在研究區(qū)域邊界點序列上任取一點作為初始點，將初始點的坐標(biāo)轉(zhuǎn)換為網(wǎng)格坐標(biāo)并標(biāo)記為一個起始的邊界網(wǎng)格點。因為網(wǎng)格坐標(biāo)是整數(shù)，所以邊界點地理坐標(biāo)轉(zhuǎn)為網(wǎng)格坐標(biāo)時，因為要進(jìn)行取整計算，故會產(chǎn)生一定的誤差。

2）步驟 2。在研究區(qū)域邊界點序列上取下一個點坐標(biāo)并轉(zhuǎn)為網(wǎng)格坐標(biāo)，此時需要判斷 2 個網(wǎng)格點之間應(yīng)該是哪些柵格點需要標(biāo)記為邊界網(wǎng)格點。本算法計算時首先將 2 個點的網(wǎng)格坐標(biāo)連成一線，如圖2 中的 P1和 P22 個點（P1，P2為構(gòu)成失量邊界點序列中相鄰的 2 個點），得到 P1到 P2網(wǎng)格 2 個方向的梯度變化 Δ x 和 Δ y，以 P1點為起點向 8 個方向移動 1 個網(wǎng)格，分別計算各個網(wǎng)格點到 P1的梯度變化 Δ xn和 Δ yn，將梯度變化最接近 P1到 P2的網(wǎng)格標(biāo)記為邊界網(wǎng)格，即將 (Δ y/Δ x - Δ yn/Δ xn)min的網(wǎng)格點標(biāo)記為邊界網(wǎng)格點，記為 Pn。再以 Pn為起點向 8 個方向移動 1 個網(wǎng)格，分別計算各個網(wǎng)格到P1的梯度變化 Δ xn和 Δ yn，得到下一個邊界網(wǎng)格點，直至 Pn到達(dá) P2點，從而形成 P1到 P2所有邊界網(wǎng)格點的標(biāo)記。

3）步驟 3。將 P2點作為 P1，從研究區(qū)域邊界點序列中取下一個相鄰的點作為 P2，重復(fù)步驟 2，直至矢量邊界的點序列處理完成。

完成這幾個步驟后所有標(biāo)記為邊界的網(wǎng)格點形成一個閉合的區(qū)域，該區(qū)域等值線生成計算的范圍。

1.3 網(wǎng)格插值計算

插值計算中網(wǎng)格范圍的確定、等值面的提取、重分類等等值線生成重要計算步驟都應(yīng)用了網(wǎng)格延展法，體現(xiàn)了算法過程簡單易行、健壯性等特點。

網(wǎng)格插值計算并不是對所有的網(wǎng)格點都進(jìn)行插值計算，所以首先要確定待插值網(wǎng)格。網(wǎng)格延展法利用需要插值的網(wǎng)格區(qū)域是連續(xù)分布的特性，從確定邊界區(qū)域內(nèi)任選一點，并遞歸遍歷及插值計算所有待插值的網(wǎng)格點。待插值網(wǎng)格點確定的遞歸搜索算法步驟如下：

1）步驟 1。從邊界標(biāo)記區(qū)域內(nèi) 1 個網(wǎng)格點開始向該柵格點的 4 個方向搜索，如果存在相鄰網(wǎng)格點且不是邊界標(biāo)記網(wǎng)格點（同樣未被標(biāo)記為待插值網(wǎng)格）則標(biāo)記為待插值網(wǎng)格點，并加入到遞歸搜索的網(wǎng)格點隊列中。

2）步驟 2。從遞歸搜索網(wǎng)格點隊列中選取 1 個網(wǎng)格點為起點向 4 個方向查找，如果存在相鄰網(wǎng)格點且不是邊界標(biāo)記點（且是未被標(biāo)記為待插值網(wǎng)格）則標(biāo)記為待插值點，同時加入到遞歸搜索網(wǎng)格點隊列中。如果當(dāng)前搜索的網(wǎng)格點 4 個方向都沒有找到新的待插值網(wǎng)格點，則將該網(wǎng)格點從遞歸搜索網(wǎng)格點隊列中刪除。

3）步驟 3。如果遞歸搜索網(wǎng)格點隊列已經(jīng)為空，則搜索過程結(jié)束，否則轉(zhuǎn)向步驟 2。

網(wǎng)格點的遞歸搜索步驟和結(jié)果如圖3 所示。圖3 a 中編號為 1 的點表示遞歸搜索的起始網(wǎng)格點，數(shù)字 2 表示遞歸遍歷的第 2 步，以此類推，最終形成圖3 b 所示的遞歸搜索結(jié)果，每個數(shù)字表示遞歸搜索的次序。對所有待插值的網(wǎng)格點進(jìn)行空間插值計算，以雨量等值線生成為例，網(wǎng)格插值計算是根據(jù)雨量觀則站點的數(shù)據(jù)對研究區(qū)域內(nèi)的所有網(wǎng)格進(jìn)行插值計算。空間插值算法有克里金（Kriging）、反距離加權(quán)平均（IDW）、樣條函數(shù)等算法[10]，本研究選用反距離加權(quán)平均插值算法作為插值方法，實際應(yīng)用時可根據(jù)實際需要替換為其他插值算法。

圖3 網(wǎng)格點的遞歸搜索

反距離加權(quán)平均算法基本公式如下：

1.4 插值結(jié)果重分類

重分類是根據(jù)空間梯度的展示需求對網(wǎng)格的數(shù)值進(jìn)行重新劃分，1 個分類值代表1 個取值范圍，提取相鄰的且分類相同的 1 組網(wǎng)格就是 1 個等值面，這不同于已有研究中一般通過等值線追蹤生成等值線的方法。

重分類算法是將插值后的數(shù)據(jù)按照一定條件進(jìn)行歸類，每一類代表該等值面的取值范圍。重分類的每一個類別的取值范圍可以根據(jù)空間梯度的應(yīng)用需求進(jìn)行定義。比如雨量等值面中 0～10 mm 雨量的網(wǎng)格點歸為一類，10～20 mm 雨量的網(wǎng)格點歸為一類，以此類推。當(dāng)研究區(qū)域內(nèi)某些位置的空間梯度很大導(dǎo)致等值線分布較密集時，由于網(wǎng)格化坐標(biāo)帶來的誤差可能會略過一些中間取值范圍的等值線，無法表現(xiàn)等值線的細(xì)部分布特點，但這種方法生成的等值線不影響宏觀信息分布的表達(dá)，也不會有常規(guī)等值線生成算法追蹤過程出現(xiàn)的等值線交叉現(xiàn)象。重分類算法比較簡單，計算時將所有插值計算網(wǎng)格點遍歷 1 次，按照指定的重分類條件對網(wǎng)格進(jìn)行重新賦值。

圖4 插值結(jié)果重分類

圖4 a 是插值后的網(wǎng)格數(shù)值，網(wǎng)格中的每個數(shù)字代表每個網(wǎng)格的雨量值；圖4 b 是重分類后的網(wǎng)格數(shù)值，網(wǎng)格中的數(shù)字代表該網(wǎng)格所屬于的分類，屬于同一分類的代表取值范圍相同的網(wǎng)格點。從圖4 b 可以看出同一分類所有網(wǎng)格構(gòu)成的連續(xù)區(qū)域即為一個等值面。

1.5 等值線邊界追蹤

等值線邊界追蹤即通過算法得到同一分類且連續(xù)的所有網(wǎng)格點所構(gòu)成的邊界。然而等值線的邊界往往很不規(guī)則，需要通過一定的算法追蹤確定。本研究同樣通過網(wǎng)格延展法將每個等值面的網(wǎng)格點先提取出來，然后對這些網(wǎng)格的外包絡(luò)邊界進(jìn)行追蹤?；诰W(wǎng)格延展法追蹤等值線的算法步驟如下：

1）步驟 1。使用網(wǎng)格延展法從重分類結(jié)果網(wǎng)格中提取分類值相同且連續(xù)分布的所有網(wǎng)格點，提取時從研究區(qū)域范圍內(nèi)任意一個點開始，以該點為起始點向四周遞歸搜索所有分類值相同的網(wǎng)格點，得到分類值相同且連續(xù)的所有網(wǎng)格點。同時將這些網(wǎng)格點在重分類結(jié)果中標(biāo)記為已空網(wǎng)格，從而避免此后被再次提取。提取的網(wǎng)格點放入一個新的臨時網(wǎng)格中（與原始構(gòu)建的網(wǎng)格大小一致，所有其他網(wǎng)格標(biāo)志為空網(wǎng)格）進(jìn)行處理。

2）步驟 2。從步驟 1 的臨時網(wǎng)格中選取 1 個邊緣網(wǎng)格點的邊作為起始邊。識別起始邊時，可以從網(wǎng)格中的 1 個空白網(wǎng)格向有數(shù)值的網(wǎng)格區(qū)域進(jìn)行遍歷查找，碰到的第 1 個邊即為該等值面的起始邊（起始邊一邊的網(wǎng)格有數(shù)值，另外一邊的網(wǎng)格標(biāo)志為空），并將該起始邊加入到邊界結(jié)果隊列中。以起始邊為起點向 4 個方向搜索下一個邊，不在邊界結(jié)果隊列中且為邊界邊的網(wǎng)格邊即為下一個符合條件的邊界邊，同樣加入到邊界結(jié)果隊列中。以此方式繼續(xù)搜索，直至所有的邊界遍歷完成（即遍歷到達(dá)起始邊），邊界結(jié)果隊列中所有的邊界數(shù)據(jù)構(gòu)成一個等值面的邊界。某一等值線追蹤計算的結(jié)果如圖5所示，粗線的所有網(wǎng)格邊構(gòu)成一個等值面邊界追蹤的結(jié)果。

3）步驟 3。重復(fù)步驟 1～2 直至所有的重分類結(jié)果被提取完，則完成等值線邊界的追蹤計算。

1.6 等值線平滑

圖5 等值線邊界追蹤結(jié)果

由網(wǎng)格邊界追蹤直接形成的等值線在視覺上不夠平滑，所以需要對等值線追蹤結(jié)果進(jìn)行平滑計算。等值線的平滑算法既要達(dá)到平滑效果，還要保證平滑后的等值線之間不能發(fā)生交叉現(xiàn)象。本研究使用一種簡單的平滑算法，可以消除梯形邊界，也能保證等值線的邊界平滑后不相交，滿足大多數(shù)應(yīng)用場合對等值線平滑的要求。

簡單的平滑算法是逐段將等值線邊界追蹤結(jié)果的每個邊界邊的中點依次連接起來，形成一個新的邊界邊，新生成的邊界從視覺上顯得更加平滑。研究區(qū)域內(nèi)所有的等值線使用同樣的平滑算法，這樣平滑后的邊依然是相切的，不會出現(xiàn)等值線相交現(xiàn)象。圖5 所示的等值線追蹤后的原始網(wǎng)格邊界構(gòu)成的梯形比較明顯，經(jīng)過一次平滑處理后得到的邊界結(jié)果如圖6 所示，如果需要更平滑的效果可以連續(xù)進(jìn)行多次平滑。

圖6 等值線的平滑處理結(jié)果

2 算法應(yīng)用

基于網(wǎng)格延展法生成等值線算法主要的操作是對規(guī)則網(wǎng)格的操作，即對二維數(shù)組的操作，所以可以使用 C#，Java，JavaScript，Python 等多種語言實現(xiàn)。本研究使用 Java 語言實現(xiàn)等值線算法，形成一個可供水利信息化應(yīng)用系統(tǒng)直接調(diào)用的算法組件。該算法組件在不需要安裝 GIS 軟件或等值線組件的情況下，可以根據(jù)雨量站和研究區(qū)域邊界等數(shù)據(jù)生成等值線，并能方便和應(yīng)用系統(tǒng)進(jìn)行集成。

以北京市懷柔區(qū)雨量等值面為例，其網(wǎng)格規(guī)模設(shè)定為 1 000×1 000，雨量觀測站點為 26 個，雨量站點信息及監(jiān)測數(shù)據(jù)由北京市懷柔區(qū)水務(wù)局提供。等值面邊界為北京市懷柔區(qū)的行政區(qū)邊界，提取自全國 1∶250 000 水利電子地圖，邊界為 Shp 格式的矢量面狀圖層。等值線算法組件依據(jù)邊界點序列，得到研究區(qū)域網(wǎng)格的坐標(biāo)范圍為東經(jīng) 116°17′～116°63′，北緯 40°41′～41°4′，以該坐標(biāo)范圍構(gòu)建矩形網(wǎng)格區(qū)域。然后根據(jù)邊界點序列數(shù)據(jù)在網(wǎng)格上進(jìn)行邊界標(biāo)記，再根據(jù) 26 個站點的雨量數(shù)據(jù)進(jìn)行空間插值，得到邊界范圍內(nèi)每個網(wǎng)格的值。按照雨量 0～100 mm范圍內(nèi)每隔 10 mm 為 1 個分類，100 mm 以上單獨作為 1 個分類，對邊界范圍內(nèi)所有網(wǎng)格點進(jìn)行重分類計算。將等值線追蹤結(jié)果按照本研究的平滑算法處理 3 次，最終生成的北京市懷柔區(qū)雨量等值線、面如圖7 所示。

圖7 等值線生成結(jié)果

網(wǎng)格規(guī)模為 1 000×1 000，站點數(shù)量為 26 個的等值線生成計算表明，等值線生成計算所需時間約為 3 s，在計算效率上能夠滿足實際業(yè)務(wù)的需要。任意改變站點的數(shù)量和雨量值，該等值線算法組件均能正常生成雨量等值線，表現(xiàn)出較好的算法健壯性。該算法組件在全國 200 多個縣級山洪災(zāi)害監(jiān)測預(yù)警系統(tǒng)中，以及多個省、市級水資源監(jiān)測管理和中小河流監(jiān)測預(yù)報預(yù)警等系統(tǒng)中得到應(yīng)用，表現(xiàn)出良好的站點數(shù)據(jù)適應(yīng)能力，并能夠較好地表達(dá)降雨、水資源的空間分布信息，滿足水利信息化相關(guān)應(yīng)用需求。

3 結(jié)語

基于網(wǎng)格延展法生成等值線算法利用研究區(qū)域內(nèi)網(wǎng)格分布的特性，在計算范圍的確定、等值面的提取上提供了一種新的方法，避免基于等值線追蹤時出現(xiàn)交叉、分叉等復(fù)雜問題，避免了等值線生成結(jié)果的裁剪計算。網(wǎng)格延展法生成等值線在重分類時，可以根據(jù)生產(chǎn)應(yīng)用需要，生成不同空間梯度的等值線。實際應(yīng)用表明，網(wǎng)格延展法生成等值線算法對不同的研究區(qū)域邊界、測站的數(shù)據(jù)具有較好的適應(yīng)性，算法表現(xiàn)為良好的健壯性。基于網(wǎng)格延展法生成等值線算法簡單，可使用多種語言實現(xiàn)，能與多數(shù)的水利信息化應(yīng)用系統(tǒng)進(jìn)行集成。但本研究生成的雨量等值線圖（如圖7 b）只能表達(dá)不同雨強(qiáng)分布區(qū)域的分界線，不能表達(dá)相同數(shù)值點構(gòu)成的等值線，在這一方面基于網(wǎng)格延展法生成等值線尚需進(jìn)一步研究和改進(jìn)。