匯水盆地算法的研究與實(shí)現(xiàn)

楊龍淵，王星捷，黃 平

（1.成都理工大學(xué)工程技術(shù)學(xué)院，四川 樂(lè)山 614007；2.四川經(jīng)準(zhǔn)檢驗(yàn)檢測(cè)集團(tuán)股份有限公司，四川廣安 638500）

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

在國(guó)內(nèi)對(duì)于匯水盆地的生成大多直接利用地理信息軟件平臺(tái)來(lái)識(shí)別和建立匯水盆地模型，如：洪明海、汪仕偉等在文獻(xiàn)[2]中利用地理信息軟件對(duì)斷面進(jìn)行了匯水盆地及面積提?。粡埾檩x、張昌民等在文獻(xiàn)[3]中利用地理信息軟件平臺(tái)生成了蘇干湖盆地匯水盆地；匡增桂等在文獻(xiàn)[4]中用匯水盆地對(duì)水合物進(jìn)行了研究。

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

1 算法設(shè)計(jì)

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

圖1 匯水盆地算法設(shè)計(jì)

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

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

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

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

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

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

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

2.3 三角網(wǎng)轉(zhuǎn)DEM

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

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

當(dāng)T1、T2、T3 同時(shí)滿足大于或等于0 時(shí)，證明P點(diǎn)位于三角形內(nèi)部。即C、P點(diǎn)位于向量左側(cè)；B、P位于向量右側(cè)；A、P位于向量向量左側(cè)時(shí)，P點(diǎn)位于三角形內(nèi)部。

其次，當(dāng)計(jì)算該點(diǎn)在DEM 中的插值時(shí)，利用了三點(diǎn)式平面方程：

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

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

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

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

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

2.4 填 洼

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

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

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

3 流向流量模塊

3.1 流向生成

流向的生成采用了D8 算法，該算法假設(shè)每個(gè)單元格的水流方向有8 個(gè)，即流入與之相鄰的8 個(gè)單元格，它通過(guò)最陡坡度法來(lái)確定水流方向。在3×3的DEM 網(wǎng)格中，計(jì)算中心網(wǎng)格與八方向網(wǎng)格間的高程落差，取高程落差最大的單元格方向?yàn)橹行膯卧竦牧鞒龇较騕1]，如圖3 所示。

圖3 八方向及其對(duì)應(yīng)像元值

具體步驟：

1）遍歷每一個(gè)高程像元，進(jìn)入步驟2）。

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

3）取8 個(gè)高程差數(shù)據(jù)遍歷，取最大值方向?qū)?yīng)數(shù)，寫(xiě)入對(duì)應(yīng)位置的流向數(shù)據(jù)像元，返回步驟1）。

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

每個(gè)像元值的生成公式如下，選擇最大值的方向?yàn)閿?shù)值輸入像元。

3.2 流量生成

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

具體步驟：

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

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

代碼設(shè)計(jì)如下：

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

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

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

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

具體步驟：

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

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

3）對(duì)應(yīng)位置的像元賦新值，清空棧，返回步驟1）。將中心像元移動(dòng)至上游像元；將流入像元的位置存入棧，返回步驟2）。

代碼設(shè)計(jì)如下：

5 匯水盆地模塊

5.1 一級(jí)匯水盆地

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

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

具體步驟：

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

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

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

5.2 多級(jí)匯水盆地

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

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

具體步驟：

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

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

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

代碼設(shè)計(jì)：

6 應(yīng)用案例展示

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

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

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

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

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

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

7 結(jié)論

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