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

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

0 引言

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

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

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

綜上所述,本文提出一種基于線性特征檢測(cè)淤地壩并修正其所在高程的方法。首先綜合淤地壩和河道的地形特征進(jìn)行區(qū)域生長和形態(tài)學(xué)操作提取河道中心線,并結(jié)合改進(jìn)LSD算法檢測(cè)結(jié)果保留淤地壩輪廓直線,最后構(gòu)建十字模型檢測(cè)淤地壩,修正所在高程并填充剩余的洼地以提取河網(wǎng)。

1 材料和方法

1.1 研究區(qū)域

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

圖1 研究樣區(qū)

1.2 實(shí)驗(yàn)方法

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

圖2 技術(shù)路線圖

1.2.1DEM平滑處理

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

圖3 濾波前后坡度對(duì)比

1.2.2河道中心線提取

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

(1)

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

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

為每個(gè)柵格設(shè)置生長狀態(tài),將未標(biāo)記為河道且大于流量累積閾值rflow的流量柵格作為種子點(diǎn),以避免流量柵格的重復(fù)生長,提高算法效率。區(qū)域生長的生長準(zhǔn)則為

(2)

ha>hb

(3)

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

柵格a向相鄰柵格b生長,如果滿足上述兩條件,則將其標(biāo)記為河道柵格,直到確定河道邊界位置。式(2)能確保淤地壩壩頂以及河道等相對(duì)平緩地區(qū)的生長,而由于淤地壩兩側(cè)坡度變化較大,利用式(3)可保證位于壩頂?shù)纳L點(diǎn)向淤地壩上下游進(jìn)行區(qū)域生長。不同坡度閾值的生長結(jié)果如圖4所示,當(dāng)坡度閾值設(shè)置為3°時(shí),區(qū)域生長無法覆蓋大型淤地壩,導(dǎo)致細(xì)化后的河道中心線發(fā)生偏移,進(jìn)而影響后續(xù)角度濾波結(jié)果的準(zhǔn)確性。當(dāng)將坡度閾值設(shè)置為4°時(shí),可以保證河道的連續(xù)性和完整性。然而,當(dāng)閾值設(shè)置為5°時(shí),溝頭區(qū)域的生長結(jié)果可能不準(zhǔn)確,因?yàn)檩^大閾值可能使生長擴(kuò)展到坡度較緩的溝頭區(qū)域。此時(shí),溝頭區(qū)域的高程較高,而式(3)會(huì)將生長過程擴(kuò)展到所有低于溝頭柵格高程的區(qū)域,進(jìn)而導(dǎo)致生長結(jié)果的不準(zhǔn)確。因此,為確保河網(wǎng)連續(xù)性和準(zhǔn)確性,取閾值E為4°。

圖4 不同坡度閾值的生長結(jié)果

對(duì)分辨率10 m DEM的河道生長結(jié)果進(jìn)行形態(tài)學(xué)閉運(yùn)算以填補(bǔ)河道中未生長的小區(qū)域,并平滑河道邊界。然而,在分辨率10 m下,淤地壩的特征無法有效表示,因此需用高分辨率DEM來檢測(cè)淤地壩。為與高分辨率DEM嵌套來定位淤地壩,將河道生長結(jié)果上采樣成高分辨率河道數(shù)據(jù),最后利用細(xì)化算法生成單柵格大小的河道中心線。

1.2.3淤地壩特征提取

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

圖5 淤地壩區(qū)域示意圖

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

圖6 不同分辨率LSD檢測(cè)結(jié)果示意圖

1.2.4有效輪廓直線提取

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

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

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

(4)

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

1.2.5十字模型檢測(cè)

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

(5)

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

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

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

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

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

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

(4)根據(jù)垂直淤地壩橫截面直線所對(duì)應(yīng)的高程數(shù)據(jù),確定高程最高的柵格坐標(biāo)M(x,y)為淤地壩壩頂。

(6)根據(jù)垂直淤地壩橫截面直線上的每個(gè)柵格,分別向橫截面的對(duì)應(yīng)方向延長lenl和lenr來構(gòu)造淤地壩所對(duì)應(yīng)的矩形。

(7)利用形態(tài)學(xué)閉運(yùn)算來填補(bǔ)遺漏的柵格。

1.2.6淤地壩高程修正

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

1.3 結(jié)果對(duì)比

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

為了評(píng)估淤地壩的檢測(cè)效果,將本方法應(yīng)用于兩處典型的實(shí)驗(yàn)區(qū),并與人工標(biāo)記的淤地壩進(jìn)行比較。采用精確率(P)、召回率(R)和F1值來驗(yàn)證淤地壩檢測(cè)的有效性。此外,在相同河網(wǎng)閾值的條件下,將本文與填充法、裂開法提取的河網(wǎng)結(jié)果進(jìn)行比較。

2 結(jié)果與討論

2.1 淤地壩檢測(cè)結(jié)果

在2個(gè)實(shí)驗(yàn)區(qū)上進(jìn)行淤地壩檢測(cè),實(shí)驗(yàn)結(jié)果如表1所示。

表1 實(shí)驗(yàn)區(qū)1和實(shí)驗(yàn)區(qū)2的實(shí)驗(yàn)結(jié)果

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

圖9 淤地壩檢測(cè)結(jié)果

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

2.2 河網(wǎng)提取結(jié)果對(duì)比

僅使用填充法提取的河網(wǎng)結(jié)果誤差最大,尤其是位于淤地壩上游區(qū)域的河網(wǎng)會(huì)出現(xiàn)大角度的偏移。這是因?yàn)橛俚貕蝺蓚?cè)的道路柵格高程較低,可能被誤判為洼地的出口柵格,導(dǎo)致河網(wǎng)沿著兩側(cè)的道路穿過淤地壩。在到達(dá)出口柵格之前,河網(wǎng)流向與淤地壩的橫截面方向基本保持平行。裂開法提取的河網(wǎng)在地形平緩的實(shí)驗(yàn)區(qū)1中能較為準(zhǔn)確地穿過淤地壩,但相較于填充法的河網(wǎng)提取結(jié)果更為曲折,此外,在復(fù)雜地形中仍存在河網(wǎng)偏移問題,這表明裂開法的效果并不理想。

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

圖10 河網(wǎng)檢測(cè)結(jié)果對(duì)比

2.3 淤地壩檢測(cè)誤差分析

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

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

3 結(jié)束語

提出了一種通過對(duì)淤地壩進(jìn)行檢測(cè)并修正其高程的自動(dòng)化方法來實(shí)現(xiàn)河網(wǎng)的有效提取。在2個(gè)典型的實(shí)驗(yàn)區(qū)上進(jìn)行驗(yàn)證,結(jié)果表明,該方法能較為準(zhǔn)確地檢測(cè)不同規(guī)模的淤地壩,并能有效解決河道偏移問題。然而,該方法仍有一定的局限性,例如潰壩或淤地壩上下游植被的影響可能會(huì)導(dǎo)致漏檢和誤檢,使提取的河網(wǎng)結(jié)果發(fā)生輕微偏移。與目前的河網(wǎng)提取方法相比,本文方法提取的河網(wǎng)結(jié)果更符合真實(shí)地形,為黃土高原淤地壩地區(qū)的水土保持規(guī)劃以及流域數(shù)字地形分析提供了良好的補(bǔ)充與支持。