基于數字高程模型的淤地壩地區河網提取方法

張宏鳴 樊世豪 陳茹雪 魚文虎 張 琦 袁琳琳

(西北農林科技大學信息工程學院,陜西楊凌 712100)

0 引言

從數字高程模型(Digital elevation model,DEM)中提取河網是數字流域[1]、分布式水文模擬[2-3]、土壤侵蝕模型建立[4]等的關鍵技術環節,可為進一步劃分流域和計算匯流面積等其他水文參數提供基礎[5]。但在黃土高原土流失嚴重的地區,淤地壩作為重要而獨特的滯洪、攔淤、蓄水工程被廣泛修筑于開闊的溝谷區域,在該地區利用DEM進行流域水系分析時,直接提取河網會受到干擾[6]。原始DEM中存在許多中心柵格高程比周圍柵格低的洼地[7]。在這些洼地內,水流不能確定正確的流向,導致水流路徑不連續甚至出現逆流現象[7],進而影響河網的正常提取。因此,國內外研究人員提出了許多DEM預處理方法以改善河網的流動路徑,從而實現河網提取[8-11]。

DEM數據的預處理是通過修正洼地,使得每個柵格都有明確的水流方向。目前,洼地的處理方法主要分為填充法和裂開法[7-8]。填充法通過抬升洼地中的柵格來修改地形,直到匹配最低出口柵格的高程[8,11-13]。而裂開法則是通過降低洼地邊緣處阻擋物的柵格高程,來模擬水流下切裂開洼地[8,14]。雖然這些方法能保證河網的連續性,但在淤地壩地區的河網提取結果存在河網偏移問題。由于淤地壩兩側道路高程較低,河網將其作為洼地出水口穿過淤地壩,而到達出水口位置之前,河網流向與淤地壩的橫截面方向基本保持平行。因此針對淤地壩流域,需開發一種新的DEM預處理方法提取河網。

CARLSON等[15]提出一種橋梁識別并修正其柵格高程的河網提取方法,與填充法相比,這種先識別后處理的方法能建立更真實的河網。然而,位于山地的梯田具有相似的地形特征,對淤地壩檢測及河網提取造成負面影響[16]。ZHU等[17]提出利用區域生長和形態學處理提取黃土高原溝道,雖然該方法并沒有考慮黃土高原地區的淤地壩,但能將淤地壩的檢測范圍縮小至河道區域,從而避免梯田的影響。此外,SOFIA等[18]提出使用直線段檢測(Line segment detector,LSD)算法[19]提取具有線性特征的梯田,為檢測具有類似線性特征的淤地壩提供了一個新的思路。

綜上所述,本文提出一種基于線性特征檢測淤地壩并修正其所在高程的方法。首先綜合淤地壩和河道的地形特征進行區域生長和形態學操作提取河道中心線,并結合改進LSD算法檢測結果保留淤地壩輪廓直線,最后構建十字模型檢測淤地壩,修正所在高程并填充剩余的洼地以提取河網。

1 材料和方法

1.1 研究區域

如圖1所示,本文選擇的兩塊實驗區均位于黃土高原地區,且具有不同的地形特征。實驗區1為王茂溝流域,位于陜西省綏德縣韭園溝中游左岸,是我國最早的治理實驗小流域之一[20]。該流域地貌復雜,以梁、峁為主,溝壑縱橫,坡陡溝深,屬典型的黃土丘陵溝壑地貌。多年平均降水量513.1 mm,汛期7—9月降水量占年降水量的73.1%。經過近70年治理,壩系布局較為完整,存在淤地壩淤滿的情況,坡耕地大幅減少,梯田面積和植被覆蓋度提高。實驗區2為周屯溝流域的一條支流,位于陜西省延安市安塞區。流域內地貌復雜多樣,境內川道狹長、梁峁遍布,土壤以黃綿土為主,物理風化作用強,土壤疏松,多粉沙,易受水力侵蝕。多年平均降水量500 mm左右,汛期7—9月降水量占年降水量的74%。該流域淤地壩和梯田相對較少,淺溝區域治理以種植草灌為主要措施。圖1中兩流域對于淤地壩區域的河網提取和分析,具有典型的代表意義。

圖1 研究樣區

1.2 實驗方法

技術路線如圖2所示,淤地壩地區的河網提取主要為:將無人機航測分辨率0.15 m DEM采用最鄰近插值法生成高分辨率DEM和低分辨率DEM。將低分辨率DEM提取河道中心線,高分辨率DEM導出山體陰影圖并提取輪廓直線。依據河道中心線進行角度濾波保留淤地壩輪廓直線,構建十字模型定位淤地壩。修正淤地壩柵格高程并填充剩余洼地來提取河網。

圖2 技術路線圖

1.2.1DEM平滑處理

高分辨率DEM提供詳細的地形信息,但其中存在的大量噪聲不僅會增加計算復雜性,還會對河道提取產生負面影響。本文主要根據河道兩側高、中間低的地形特征提取河道,不過多關注河道內的地形細節。有研究表明,分辨率10 m DEM的坡度統計百分比與分辨率5 m DEM相比變化不大,足以滿足地形和水文建模應用的需求[21-23]。分辨率更低的DEM會導致更多地形信息的丟失,使得區域生長結果無法有效反映河道的實際情況,因此本方法使用分辨率10 m DEM數據進行區域生長。分辨率10 m DEM仍存在少部分噪聲,為了去除這些噪聲并保持河道和山谷的邊緣信息,本文使用中值濾波窗口(3×3)進行濾波。如圖3所示,濾波能有效去除河道內部噪聲,且能降低部分淤地壩所在的柵格高程,從而有利于后續區域生長提取連續的河道。

圖3 濾波前后坡度對比

1.2.2河道中心線提取

對低分辨率DEM的洼地進行填充,并用最陡坡降計算流向和流量。經過中值濾波和填充處理后,河道的整體地形相對平坦。利用坡度和高程之間的關系進行區域生長得到河道,坡度S定義為

(1)

式中H——相鄰柵格的高程差

L——相鄰柵格之間的水平距離

為每個柵格設置生長狀態,將未標記為河道且大于流量累積閾值rflow的流量柵格作為種子點,以避免流量柵格的重復生長,提高算法效率。區域生長的生長準則為

(2)

ha>hb

(3)

式中h——柵格高程E——坡度閾值

柵格a向相鄰柵格b生長,如果滿足上述兩條件,則將其標記為河道柵格,直到確定河道邊界位置。式(2)能確保淤地壩壩頂以及河道等相對平緩地區的生長,而由于淤地壩兩側坡度變化較大,利用式(3)可保證位于壩頂的生長點向淤地壩上下游進行區域生長。不同坡度閾值的生長結果如圖4所示,當坡度閾值設置為3°時,區域生長無法覆蓋大型淤地壩,導致細化后的河道中心線發生偏移,進而影響后續角度濾波結果的準確性。當將坡度閾值設置為4°時,可以保證河道的連續性和完整性。然而,當閾值設置為5°時,溝頭區域的生長結果可能不準確,因為較大閾值可能使生長擴展到坡度較緩的溝頭區域。此時,溝頭區域的高程較高,而式(3)會將生長過程擴展到所有低于溝頭柵格高程的區域,進而導致生長結果的不準確。因此,為確保河網連續性和準確性,取閾值E為4°。

圖4 不同坡度閾值的生長結果

對分辨率10 m DEM的河道生長結果進行形態學閉運算以填補河道中未生長的小區域,并平滑河道邊界。然而,在分辨率10 m下,淤地壩的特征無法有效表示,因此需用高分辨率DEM來檢測淤地壩。為與高分辨率DEM嵌套來定位淤地壩,將河道生長結果上采樣成高分辨率河道數據,最后利用細化算法生成單柵格大小的河道中心線。

1.2.3淤地壩特征提取

如圖5所示,為提高淤地壩的提取精度,本文總結淤地壩的4個特征:①淤地壩修建在截面兩側具有陡峭山地的河道上。②淤地壩與河道的流向近乎垂直。③淤地壩為橫向線性結構壩工建筑物。④淤地壩的壩頂、上游和下游地勢相對平坦。

圖5 淤地壩區域示意圖

由于淤地壩整體呈線性,本文使用LSD算法對分辨率1、2、5、10 m的山體陰影圖進行輪廓直線提取。如圖6所示,分辨率10 m山體陰影圖只能檢測到河道兩側的輪廓直線,分辨率5 m雖能檢測到大型淤地壩的輪廓直線,但對于中小型淤地壩的檢測并不理想。在分辨率2 m下,能很好地檢測各種規模淤地壩的輪廓直線。而在分辨率1 m下,檢測到的淤地壩輪廓直線與分辨率2 m相比并無較大區別,但河道內的噪聲輪廓直線數量卻大大增加,使得提取淤地壩輪廓直線的方法更為復雜,檢測效果也會受到影響。因此,本文選擇分辨率2 m DEM進行淤地壩檢測,以獲得較好的檢測效果。

圖6 不同分辨率LSD檢測結果示意圖

1.2.4有效輪廓直線提取

利用淤地壩與河道近似垂直的位置關系,對輪廓直線進行角度濾波,以保留與河道中心線相交且近乎垂直的輪廓直線。如圖7所示,淤地壩橫跨河道而建,淤地壩的輪廓直線與河道的輪廓直線可能相鄰。為了正確計算淤地壩輪廓直線的傾斜角,本文改進LSD算法分離傾斜角不同的直線。改進方法如圖7a所示,由于檢測到的輪廓直線按先后順序輸出,若差值在0°～5°范圍內,則認為這兩條直線屬于同一直線,保留柵格F并接著檢測下一柵格是否滿足該條件;相反則將該柵格丟棄,以分離不同傾斜角的輪廓直線。

分離輪廓直線后,利用最小二乘法擬合輪廓直線l與河道中心線j的傾斜角,分別為θl和θj。如果這兩傾斜角的差值滿足

|θl-θj|∈[90°-α,90°+α]

(4)

則認為直線l屬于淤地壩輪廓直線;反之則認為直線l為噪聲輪廓直線并丟棄。由于淤地壩和河道中心線并非完全垂直,且基于柵格的傾斜角計算存在誤差,本文將α設置為30°,盡可能避免淤地壩輪廓直線的丟失及噪聲輪廓直線的影響。

1.2.5十字模型檢測

角度濾波后仍存在少部分的噪聲直線,需對每條輪廓直線構建十字模型以確定最后淤地壩的位置。由于多條輪廓直線表示同一個淤地壩,將其合并以避免十字模型的重復構建。合并的判斷依據:若距離較近的兩直線傾斜角差值小于5°,則合并兩直線。如圖8所示,確定合并之后的輪廓直線中點U,則有

(5)

k——垂直淤地壩橫截面的直線斜率

根據淤地壩壩頂和上下游區域地勢平坦的特征,設置平緩值SA以確定垂直淤地壩橫截面直線的延伸長度,平緩值SA只根據相鄰兩柵格整型高程是否相等來變化。

建立十字模型檢測淤地壩具體步驟為:

(1)將輪廓直線中點U作為垂直淤地壩橫截面直線的起點,并存儲在鏈表的第1個節點。利用式(5)計算垂直淤地壩橫截面的直線斜率k,其中標準高程ch初始化為中心點U的高程,平緩值SA初始化為3,可根據分辨率進行調整。

(2)根據Bresenham算法沿斜率k確定下一點的坐標u1。判斷u1與ch的整形高程值是否相等,若相等且SA=0,則停止向該方向的延伸并執行步驟(3),否則將u1存儲至鏈表,并把ch更新為u1的高程;相反若兩柵格高程相等且SA>0,則將SA減1,并重復步驟(2),計算下一個點u2。

(3)根據反向斜率-k,重復步驟(1)～(2)。

(4)根據垂直淤地壩橫截面直線所對應的高程數據,確定高程最高的柵格坐標M(x,y)為淤地壩壩頂。

(6)根據垂直淤地壩橫截面直線上的每個柵格,分別向橫截面的對應方向延長lenl和lenr來構造淤地壩所對應的矩形。

(7)利用形態學閉運算來填補遺漏的柵格。

1.2.6淤地壩高程修正

確定淤地壩位置后,使用侵蝕操作符修正淤地壩高程。侵蝕操作符是一個9×9的窗口,以淤地壩所在柵格為窗口中心,窗口內的高程最小值替代為窗口中心的高程。根據經驗,連續5個侵蝕操作符能完全處理淤地壩高程,這些參數可根據分辨率進行調整。除淤地壩外,一些因其他原因形成的洼地仍然存在,使用填充法處理剩余的洼地,最后采用D8算法提取河網。

1.3 結果對比

本研究將所述方法與LINDSAY[24]提出改進的填充法(填充法)和最小代價裂開法(裂開法)進行對比,其中裂開法的最大搜索距離設為100 m。填充法對填充洼地后的平坦地區進行微小抬升來避免平行流向問題,而裂開法提取的河網更自然,能夠在合適的位置穿過洼地[24]。

為了評估淤地壩的檢測效果,將本方法應用于兩處典型的實驗區,并與人工標記的淤地壩進行比較。采用精確率(P)、召回率(R)和F1值來驗證淤地壩檢測的有效性。此外,在相同河網閾值的條件下,將本文與填充法、裂開法提取的河網結果進行比較。

2 結果與討論

2.1 淤地壩檢測結果

在2個實驗區上進行淤地壩檢測,實驗結果如表1所示。

表1 實驗區1和實驗區2的實驗結果

該方法在2個實驗區中檢測淤地壩的F1值分別為86.67%和86.95%,結果表明該方法能較為準確地檢測淤地壩。淤地壩的檢測結果如圖9所示,所述方法中改進LSD算法,確保在進行角度濾波時能正確計算淤地壩輪廓直線的傾斜角,且角度濾波操作能有效避免位于山地中梯田輪廓直線的干擾。此外,淤地壩存在淤滿和未淤滿的兩種地形特征,如圖8所示上游虛線部分表示淤地壩已蓄滿,其上游高程近似等于壩頂高程。本文構建的十字模型能有效識別這兩種情況的淤地壩,因為十字模型中設置的平緩值SA僅用于判斷淤地壩上下游高程是否平緩,而不考慮其上下游和壩頂之間的高程關系,從而提高了淤地壩的整體檢測效果。本方法利用線性特征檢測黃土高原地區的淤地壩,但該方法也可拓展到具有相似線性特征的建筑物,如橋梁、道路等。

圖9 淤地壩檢測結果

由于淤地壩地區地形復雜,檢測過程中可能存在淤地壩誤檢和漏檢的情況。位于河道交匯處兩側的道路輪廓直線可能與河道中心線近似垂直,并且道路與淤滿淤地壩的特征相似,導致誤檢情況的發生。此外,潰壩或淤地壩周圍的植被可能會破壞淤地壩的地形特征,導致部分淤地壩漏檢。當淤地壩上下游存在植被時,植被高度會使所在區域的高程整體高于真實地面高程,垂直誤差高達10 m以上,這些高程更高的植被柵格會代替原本的壩頂,導致在十字模型建模時會因橫截面高程差過大而將被舍棄。

2.2 河網提取結果對比

僅使用填充法提取的河網結果誤差最大,尤其是位于淤地壩上游區域的河網會出現大角度的偏移。這是因為淤地壩兩側的道路柵格高程較低,可能被誤判為洼地的出口柵格,導致河網沿著兩側的道路穿過淤地壩。在到達出口柵格之前,河網流向與淤地壩的橫截面方向基本保持平行。裂開法提取的河網在地形平緩的實驗區1中能較為準確地穿過淤地壩,但相較于填充法的河網提取結果更為曲折,此外,在復雜地形中仍存在河網偏移問題,這表明裂開法的效果并不理想。

如圖10所示,本文通過先修正淤地壩所在高程,再使用填充法提取河網,能有效解決河網偏移問題。相較于僅使用填充法直接改變整個淤地壩上游高程,修正后的淤地壩上游不存在洼地,提取的河網結果能更好地反映淤地壩上游的真實地形。與裂開法的河網提取結果相比,本文方法的河網提取結果更為平滑,且能有效地穿過淤地壩。在平坦河道中,這3種方法都存在共有的問題,即當多條支流匯入主干河網時,支流可能會沿著河道穿過一段距離后才匯入主干河網,這主要由于平坦區域的河網流向分配不理想所導致的。

圖10 河網檢測結果對比

2.3 淤地壩檢測誤差分析

淤地壩存在誤檢和漏檢的情況,需要使用填充法來填補剩余的洼地。如圖11a所示,因潰壩而導致淤地壩漏檢時對河網的影響最小,因為破損淤地壩上游并未形成洼地,河網會沿破損處穿過淤地壩。如圖11b所示,當受淤地壩下游植被影響而導致漏檢時,河網提取結果仍存在偏移問題,而由于淤地壩檢測的精確率較高,整體的河網提取結果仍能得到有效的改善。圖11c中誤判為淤地壩區域的最小高程替代成整個區域的高程,改變了原始區域的地形,導致該區域的河網發生輕微偏移。將來可利用高程關系判斷淤地壩的上下游,在其兩側確定一條合理的單柵格路徑,并只修正路徑所在的柵格高程,這樣既能減少對DEM的修改,又能盡可能降低誤判對河網提取結果的影響。

圖11 淤地壩漏檢、誤檢的影響

3 結束語

提出了一種通過對淤地壩進行檢測并修正其高程的自動化方法來實現河網的有效提取。在2個典型的實驗區上進行驗證,結果表明,該方法能較為準確地檢測不同規模的淤地壩,并能有效解決河道偏移問題。然而,該方法仍有一定的局限性,例如潰壩或淤地壩上下游植被的影響可能會導致漏檢和誤檢,使提取的河網結果發生輕微偏移。與目前的河網提取方法相比,本文方法提取的河網結果更符合真實地形,為黃土高原淤地壩地區的水土保持規劃以及流域數字地形分析提供了良好的補充與支持。