楊龍淵，王星捷，黃 平

（1.成都理工大學工程技術學院，四川 樂山 614007；2.四川經(jīng)準檢驗檢測集團股份有限公司，四川廣安 638500）

匯水盆地是地球表面以分水嶺為邊界的自然降水匯集區(qū)域[1]。匯水盆地表示的是在分水嶺作用下，自然降水所形成的地表徑流匯集范圍。自然水系流域的形成嚴格受匯水盆地范圍的制約。匯水盆地模型以數(shù)字模型方式對匯水盆地進行描述，是礦產(chǎn)資源調查、預測與評價工作不可缺少的基礎數(shù)據(jù)之一。

在國內對于匯水盆地的生成大多直接利用地理信息軟件平臺來識別和建立匯水盆地模型，如：洪明海、汪仕偉等在文獻[2]中利用地理信息軟件對斷面進行了匯水盆地及面積提??；張祥輝、張昌民等在文獻[3]中利用地理信息軟件平臺生成了蘇干湖盆地匯水盆地；匡增桂等在文獻[4]中用匯水盆地對水合物進行了研究。

文中采用VC++語言為開發(fā)語言引用開源柵格空間數(shù)據(jù)轉換庫（Geospatial Data Abstraction Library，GDAL），實現(xiàn)了匯水盆地模型的生成，并對主要算法進行了研究。整個算法分為數(shù)據(jù)輸入模塊、DEM 生成模塊、流向流量模塊、河網(wǎng)生成模塊以及數(shù)字盆地模塊。

1 算法設計

文中算法先將高程數(shù)據(jù)生成TIN 三角網(wǎng)，然后生成DEM 圖，通過DEM 圖得出研究區(qū)域的水流方向以及水流量，再通過水流方向初步得出匯水盆地模型，同時將水流量按閾值提取出柵格水系網(wǎng)，最后由水系網(wǎng)以及水流方向綜合生成低級別匯水盆地模型（多級匯水盆地模型），具體算法設計如圖1 所示。

圖1 匯水盆地算法設計

2 數(shù)據(jù)輸入模塊及DEM生成模塊

2.1 數(shù)據(jù)輸入處理

數(shù)據(jù)輸入包括點數(shù)據(jù)與線數(shù)據(jù)，點數(shù)據(jù)與線數(shù)據(jù)的加載利用gdal 庫來讀取矢量高程數(shù)據(jù)，然后將實際地理坐標轉化為屏幕坐標顯示到電腦屏幕上。世界坐標向屏幕坐標轉化公式如下：

其中，x、y是屏幕坐標，X、Y是世界坐標，minX、maxY是最小、最大世界坐標。

2.2 TIN三角網(wǎng)生成

傳統(tǒng)的算法是通過凸包的Delaunay-TIN 生成算法生成TIN 三角網(wǎng)，該算法方便實現(xiàn)，但是速度較慢。文中利用改進的Delaunay 三角剖面算法對TIN三角網(wǎng)算法進行優(yōu)化[5-6]。

基本的思路：生成凸包、檢測凸包、向凸包中生成Delaunay 三角形，再檢測凸包和三角形。其中，在第三個步驟又將點集劃分為多個區(qū)域來構成三角網(wǎng)，然后采用重心九格檢測法檢測三角形是否為Delaunay 三角形。

2.3 三角網(wǎng)轉DEM

假設有三角網(wǎng)中的一個三角形ABC，三角網(wǎng)高程關系如圖2 所示，三角網(wǎng)轉DEM 需要首先判斷插值點P是否在三角形內[7]，根據(jù)混合積計算如下：

圖2 三角網(wǎng)高程關系

當T1、T2、T3 同時滿足大于或等于0 時，證明P點位于三角形內部。即C、P點位于向量左側；B、P位于向量右側；A、P位于向量向量左側時，P點位于三角形內部。

其次，當計算該點在DEM 中的插值時，利用了三點式平面方程：

給定XX、YY值可確定一個ZZ值，從而得出該點的高程值。具體步驟如下：

1）以高程點的最大X、Y值作為DEM 的寬和高，DEM 初始化高程為0。

2）遍歷DEM，判斷每個DEM 像元位置是否在該三角網(wǎng)的三角形中，是則進入步驟3），否則繼續(xù)步驟2）。

3）通過插值得出該位置的像元值，返回步驟2）。

通過以上算法實現(xiàn)三角網(wǎng)轉DEM。

2.4 填 洼

洼地在自然界的地形中比較少見，但在生成DEM高程的過程中，由于數(shù)據(jù)的輸入錯誤，高程內插值精度以及網(wǎng)格單元內高程信息取平均往往會出現(xiàn)偽洼地[7-8]。填洼對處理偽洼地、使流域水系順暢連通方面具有重要作用。

文中采用M&V 算法,該算法于1993 年由Moarn和Vezina 提出，2001 年由Olivier Planehon 和Frederic DarbouX 實現(xiàn)。王建平等在《一種新的DEM 填洼處理算法》中詳細論述了實現(xiàn)填洼的兩種算法，即J&D算法以及M&V算法，并比較了兩種算法之間的差異。

M&V 算法思路：除邊界外初始化一個極大高程，將原始高程覆蓋，然后通過迭代新的高程進行逐步收斂，迭代結束后新高程必須滿足大于或等于原始高程，且任意一點高程都有一條遞減的路徑通向邊界[8-9]。

3 流向流量模塊

3.1 流向生成

流向的生成采用了D8 算法，該算法假設每個單元格的水流方向有8 個，即流入與之相鄰的8 個單元格，它通過最陡坡度法來確定水流方向。在3×3的DEM 網(wǎng)格中，計算中心網(wǎng)格與八方向網(wǎng)格間的高程落差，取高程落差最大的單元格方向為中心單元格的流出方向[1]，如圖3 所示。

圖3 八方向及其對應像元值

具體步驟：

1）遍歷每一個高程像元，進入步驟2）。

2）以步驟1）中像元為中心像元，遍歷八方向上的像元值，如圖3 所示，如果像元不是無意義像元，則計算中心像元與該方向上像元的高程差并保存，進入步驟3）。

3）取8 個高程差數(shù)據(jù)遍歷，取最大值方向對應數(shù)，寫入對應位置的流向數(shù)據(jù)像元，返回步驟1）。

在步驟3）中，需要考慮一個情況，如果碰到計算得出的流向是流向中心像元，數(shù)據(jù)表現(xiàn)為小于零的情況，如果按最大值寫入流向數(shù)據(jù)，則會導致結果與實際不符，此時將該像元值設為-1，并在后續(xù)進行特定處理。

每個像元值的生成公式如下，選擇最大值的方向為數(shù)值輸入像元。

3.2 流量生成

流量數(shù)據(jù)利用流向數(shù)據(jù)追蹤水流方向，在水流方向上將每個DEM 單元格進行累加。

具體步驟：

1）遍歷每一個流向數(shù)據(jù)，如果流向數(shù)據(jù)不為無，則進入步驟2）。

2）如果像元值為-1 和無值，則返回步驟1）；如果是其他值則對指向像元位置對應的流量像元值加1，將中心像元移動至上游像元，繼續(xù)進行步驟2）。

代碼設計如下：

4 河網(wǎng)生成模塊

河流鏈接是形成多級匯水盆地模型的關鍵數(shù)據(jù)，因為一級匯水盆地區(qū)域是一個比較大的流域，在很多水文分析中，需要將大的流域盆地劃分為更小的流域單元進行分析，因此，需要河流鏈接按大的流域進行劃分，隨后進行流域的分割。

流域的分割是確定河流出水口即匯水點的位置。如果沒有出水點的柵格或矢量數(shù)據(jù),可以采用河流連接數(shù)據(jù)作為匯水區(qū)的出水點數(shù)據(jù)。

河流鏈接的生成需要按規(guī)定閾值將流量數(shù)據(jù)分割出河網(wǎng)，再沿著河網(wǎng)向上游搜索，遇到分支則標記為新河流。

具體步驟：

1）遍歷流量數(shù)據(jù)，設置像元值閾值大于210，且未賦值的河流鏈接標識的數(shù)據(jù)像元值，進入步驟2）。

2）遍歷該像元八方向的像元值，如果流入中心像元數(shù)為1，則進入步驟3）。如果大于1，則對棧中進行賦值河流鏈接標識操作。

3）對應位置的像元賦新值，清空棧，返回步驟1）。將中心像元移動至上游像元；將流入像元的位置存入棧，返回步驟2）。

代碼設計如下：

5 匯水盆地模塊

5.1 一級匯水盆地

一級匯水盆地模型的生成難度與多級匯水盆地模型較相似，但在實際應用中一級匯水盆地模型使用較少，且一級匯水盆地模型生成僅需要流向數(shù)據(jù)[10-13]。

盆地模型生成時根據(jù)流向數(shù)據(jù)的八領域不斷遞歸向上游搜索，直到?jīng)]有流入像元時停止搜索。

具體步驟：

1）遍歷流向數(shù)據(jù)，如果盆地數(shù)據(jù)未賦值、流向數(shù)據(jù)不為無，則進入步驟2）。

2）遍歷該像元八方向的像元值，如果沒有流入像元則進入步驟3）；其他情況則將所有流入中心像元的上游像元位置以及像元值存入棧，將中心像元移動至上游像元，繼續(xù)進行步驟2）。

3）尋找棧中最大的像元值，并且對匯水盆地中對應位置賦該值，清空棧后返回步驟1）。

5.2 多級匯水盆地

多級匯水盆地模型的實現(xiàn)是建立在河流鏈接的基礎上的，該文生成的河流鏈接數(shù)據(jù)中隱含著河網(wǎng)中每一條河流段的鏈接信息，包括河流段的起點和終點，就是該匯水區(qū)域出水口所在的位置[14-16]。

多級匯水盆地模型仍然是向上不斷遞歸搜索的。與一級匯水盆地模型不同的是需要在匯水點位置終止。

具體步驟：

1）遍歷河流鏈接，如果是在河網(wǎng)上，且在多級匯水盆地區(qū)域未賦值，則進入步驟2）。

2）遍歷該像元八方向的像元值，如果有流入像元，則進入步驟3）；否則將棧中與多級匯水盆地對應位置的像元賦該段河流鏈接值，清空棧，返回步驟1）。

3）流入像元未賦值或值等于該段河流鏈接像元值，則將流入像元值以及位置存入棧中，將中心像元移動至上游像元，返回步驟2）。

代碼設計：

6 應用案例展示

文中以米拉山矢量高程點為基礎數(shù)據(jù)，經(jīng)過TIN 三角網(wǎng)轉化生成柵格DEM 數(shù)據(jù)，用于后續(xù)的匯水盆地生成。實驗設備為普通i5 四核的PC 機，其在執(zhí)行過程中運行速度較快，結果準確性高。DEM 數(shù)據(jù)為249 行×207 列的點陣式數(shù)字高程數(shù)據(jù)，其中無值區(qū)域用像元值0 表示。生成的高程點、TIN 三角網(wǎng)、DEM 高程柵格數(shù)據(jù)如圖4 所示。

圖4 高程點、三角網(wǎng)、DEM混合圖

河網(wǎng)以流量圖為基礎，在流量圖的基礎上設定閾值為210 后，得出柵格河網(wǎng)圖，一級與多級匯水盆地皆與柵格河網(wǎng)的結構與分布一致，即一塊主河流存在于一塊一級匯水盆地之中，每條直流存在于一塊多級匯水盆地之中。通過文中算法實現(xiàn)的一級匯水盆地模型及河網(wǎng)如圖5 所示，多級匯水盆地模型及河網(wǎng)效果如圖6 所示。

圖5 一級匯水盆地模型及河網(wǎng)效果

圖6 多級匯水盆地模型及河網(wǎng)效果

與ArcGIS 軟件中工箱處理的效果對比，該文算法實現(xiàn)的一級匯水盆地模型和多級匯水盆地模型的效果一致，其中河網(wǎng)數(shù)量、河網(wǎng)分級和盆地分級的數(shù)量保持一致，充分證明了文中匯水盆地生成算法的準確性。

7 結論

文中主要研究了匯水盆地模型生成算法，生成DEM、填洼處理、流向處理、流量分析、河網(wǎng)、一級匯水盆地模型、多級匯水盆地模型。通過多個算法的融合，最終實現(xiàn)了匯水盆地數(shù)字模擬。對每個部分的算法都與ArcGIS 軟件處理的結果進行了對比，保證了結果的準確性和一致性。在一級、多級匯水盆地模型和河流鏈接算法的研究中，通過對基礎理論的研究，設計了簡潔而實用的算法，并進行了實驗驗證，為該方向算法的研究提供了全新的想法。文中實現(xiàn)的匯水盆地算法具有自主創(chuàng)新性，也具有一定的研究參考價值，為匯水盆地模型生成算法的研究和改進提供了一種新的思路。