基于Burn-in和D8法的巖溶區(qū)流域邊界提取研究

123123123123123

(1.廣西大學(xué)土木建筑工程學(xué)院， 廣西南寧530004；2.工程防災(zāi)與結(jié)構(gòu)安全教育部重點(diǎn)實(shí)驗(yàn)室， 廣西南寧530004；3.廣西防災(zāi)減災(zāi)與工程安全重點(diǎn)實(shí)驗(yàn)室， 廣西南寧530004)

0 引言

流域邊界提取是區(qū)域水量平衡、水資源評價、水環(huán)境分析等研究的基礎(chǔ)。隨著計算機(jī)科學(xué)的發(fā)展，DEM(數(shù)字高程模型)已取代紙質(zhì)地形圖成為最主要的流域特征提取數(shù)據(jù)源[1-2]。為提高流域特征提取精度，不少學(xué)者進(jìn)行了相應(yīng)研究，例如TURCOTTE等[3]將Chaudière流域的數(shù)字河流信息與DEM進(jìn)行融合，對柵格的流向進(jìn)行修正，從而提高了該流域河網(wǎng)的提取精度；劉學(xué)軍等[4]以比匯水面積為研究對象，對D8、Rho8、Dinf、FMFD和DEMON這5種流向確定方法進(jìn)行對比，結(jié)果表明Dinf和DEMON的模擬結(jié)果要優(yōu)于其他算法；彭培等[5]利用Agree算法對湖北陽新縣地區(qū)的DEM進(jìn)行預(yù)處理，使得提取出的河網(wǎng)更接近真實(shí)情況；張維等[6]將數(shù)字化河道信息與DEM進(jìn)行融合，并利用D8法在平原區(qū)的秦淮河流域特征提取中得到了較好的效果；吳用等[7]利用Burn-in法對內(nèi)蒙古河套平原地區(qū)的DEM進(jìn)行預(yù)處理，最終提高了平原區(qū)流域邊界提取的準(zhǔn)確性。目前對基于DEM的流域特征提取研究，主要集中在河網(wǎng)的識別上，而對于流域邊界提取研究相對較少。特別在巖溶區(qū)流域，由于地下巖溶發(fā)育，明暗河流錯綜交替，流域邊界往往難以準(zhǔn)確劃定，從而對巖溶區(qū)的水量平衡、水資源評價、水環(huán)境分析等相關(guān)研究帶來較大影響。因此，對巖溶區(qū)流域邊界劃分研究具有重要意義。論文基于Burn-in法和D8法對巖溶區(qū)流域邊界提取進(jìn)行研究，以資為巖溶區(qū)流域水文水資源相關(guān)的基礎(chǔ)性研究提供服務(wù)。

1 研究區(qū)概況

圖1 澄碧河流域水系分布Fig.1 Chengbihe reservoir water distribution

澄碧河流域?yàn)閺V西典型巖溶區(qū)流域，流域內(nèi)主要河流為澄碧河。澄碧河屬西江水系3級支流，河流發(fā)源于凌云縣青龍山北麓，沿青龍山腳向北流，在弄桃附近潛入地下，反向伏流約20 km后于凌云縣城東側(cè)2 km的水源洞流出形成明河；再折往東南方向，流過約30 km后穿過八仙洞，流出后轉(zhuǎn)向西南，經(jīng)過15 km后進(jìn)入浩坤溶洞；河流從浩坤溶洞流出后，往下游再流經(jīng)坡帖、平塘、林河后到達(dá)澄碧河水庫壩首，壩下河流再流經(jīng)7 km最終匯入右江[8]，澄碧河流域水系分布如圖1所示。流域在浩坤—弄塘以上部分多見伏流河、落水洞和天窗，為典型的中國西南巖溶區(qū)地貌；弄塘以下地區(qū)植被覆蓋良好，林木密布，屬丘陵地貌。流域從2001年起布設(shè)了13個水情測報站點(diǎn)，其中壩上、平塘、浩坤和下甲為水文站，下塘、百練、林河、凌云、朝里、弄塘、介福和東和為雨量站，壩下為水位站[9-11]。

2 研究方法與數(shù)據(jù)

2.1 Burn-in法

Burn-in法由SAUNDERS[12]于1996年提出，該算法也被稱為Stream Burning法、河道燒錄法或凹陷化算法。其原理是通過強(qiáng)制降低河網(wǎng)所在位置的DEM柵格高程，從而使得河網(wǎng)部分的匯流能力提高(原理示意見圖2)。Burn-in法具體操作可分為4步：①輸入河網(wǎng)的矢量信息；②將河網(wǎng)矢量轉(zhuǎn)化為柵格；③降低河網(wǎng)柵格對應(yīng)DEM高程；④對被降低DEM周邊的柵格進(jìn)行平滑化處理。本次研究Burn-in法使用的巖溶區(qū)地下河網(wǎng)信息來源于2013年澄碧河水庫除險加固設(shè)計報告[13]。

2.2 D8法

D8法由O’CHALLAGHAN等[14]于1984年提出，該方法又稱為最大距離權(quán)落差法或坡面流累積法，屬于單一流向算法，通常運(yùn)用于流域河網(wǎng)提取或邊界劃定。D8法在劃定流域邊界時，以柵格DEM數(shù)據(jù)為基礎(chǔ)，假定每一個柵格周邊有8種可能流向，而實(shí)際流向是沿著坡度最陡的方向。通過確定每個柵格的流向，建立流向矩陣，在指定流域出口位置后，通過對流向矩陣進(jìn)行搜索，找出所有匯水至出口位置的柵格，從而確定流域邊界(原理示意見圖3)。D8法具有原理簡明、操作簡單、運(yùn)算效率高的特點(diǎn)，被集成到了ArcGIS水文分析工具箱、Arc-SWAT、Arc Hydro Tool等軟件當(dāng)中，運(yùn)用較為廣泛。本次研究使用的原始DEM來源于地理空間數(shù)據(jù)云網(wǎng)站(http://www.gscloud.cn/)的ASTER GDEM 30 m分辨率數(shù)據(jù)。

2.3 SWAT模型

為了檢證提取流邊界的合理性，論文采用SWAT模型來進(jìn)行徑流模擬檢驗(yàn)。SWAT模型于20世紀(jì)90年代由美國農(nóng)業(yè)部農(nóng)業(yè)研究中心(USDA-ARS)開發(fā)提出，模型可以運(yùn)用于流域水量平平衡、地表徑流模擬、輸沙量模擬、污染物模擬等方面。SWAT模型是基于水文響應(yīng)單元(HRU)的分布式水文模型[15]，進(jìn)行地表徑流模擬時需要輸入DEM數(shù)字高程信息、土地覆蓋類型數(shù)據(jù)、土壤類型數(shù)據(jù)以及降雨等氣象資料。

本文驅(qū)動SWAT模型的DEM數(shù)據(jù)來源于ASTER GDEM 30 m，土地覆蓋類型數(shù)據(jù)來源于黑河數(shù)據(jù)中心網(wǎng)站的1km精度IGBP DISCover數(shù)據(jù)集，土壤類型數(shù)據(jù)來源于聯(lián)合國糧農(nóng)組織FAO網(wǎng)站的世界土壤數(shù)據(jù)庫HWSD，降雨數(shù)據(jù)來源于澄碧河水庫12個雨量站2002～2014年實(shí)測降雨數(shù)據(jù)，氣溫、輻射、風(fēng)速等天氣數(shù)據(jù)來自CMADS(SWAT模型中國大氣同化驅(qū)動數(shù)據(jù)集)和SWAT-Weather(SWAT天氣發(fā)生器)的模擬數(shù)據(jù)，徑流數(shù)據(jù)來自壩首水文站2002～2014年的實(shí)測資料。

SWAT模型率定參數(shù)選用土壤容重(SOL_BD)、SCS徑流曲線系數(shù)(CN2)、基流消退系數(shù)(ALPHA_BF)、主河道曼寧系數(shù)(CH_N2)和土壤蒸發(fā)補(bǔ)償系數(shù)(ESCO)[16]；模型率定期選取為2003～2011年；模型驗(yàn)證期選取為2012年～2014年；模型迭代次數(shù)選取為30次；模擬效果利用R2(決定系數(shù))和Ef(納什效率系數(shù))檢驗(yàn)，計算公式如下：

(1)

(2)

3 結(jié)果與分析

3.1 流域邊界提取

研究首先使用D8法對原始DEM進(jìn)行流域提取，投影坐標(biāo)選擇WGS84-UTM-48N，在經(jīng)過填洼、指定流向、累積匯流、指定出口、流域提取等操作后，得到澄碧河流域劃分如圖4所示，將該流域劃分命名為流域A，對應(yīng)面積為1 558.14 km2。之后，采用Burn-in法對流域DEM進(jìn)行預(yù)處理，利用Arc-SWAT軟件的Burn-in工具將流域的地下河分布輸入模型并使地下河對應(yīng)位置的DEM高程下降，此時地下河被“顯化”成為地表河流，在此基礎(chǔ)上再使用D8法來確定地下河對應(yīng)的巖溶區(qū)流域邊界，補(bǔ)充后的澄碧河水庫流域邊界如圖5所示，將該域劃分命名為流域B，對應(yīng)流域面積為2 035.26 km2。Burn-in法處理后，提取出的流域面積擴(kuò)大了447.12 km2，增加的部分集中在凌云、東和、介福等地，土壤類型主要為鐵硅鋁土、黃棕壤土和裸露巖石。

圖4 原始DEM流域邊界
Fig.4 Watershed from Original DEM

圖5 DEM Burn-in處理后流域邊界
Fig.5 Watershed from Burn-in Treatment DEM

3.2 流取邊界合理性分析

基于SWAT模型的2種逐月徑流模擬結(jié)果見圖6；徑流模擬的R2和Ef見表1；國際常用的Ef評價標(biāo)準(zhǔn)[17]見表2。

圖6 徑流模擬結(jié)果Fig.6 Runoff simulation results

表1 徑流模擬系數(shù)對比Tab.1 Comparison of runoff simulation coefficient

表2 Ef評價標(biāo)準(zhǔn)Tab.2 Ef evaluation criterion

由圖5可以看到，流域A與流域B徑流模擬結(jié)果整體上與實(shí)測值較為吻合。在率定年的汛期，流域A的峰值模擬略微呈現(xiàn)低估，流域B略微呈現(xiàn)高估；在率定年的非汛期，流域A與流域B模擬結(jié)果均存在低估；在驗(yàn)證期，流域A的模擬結(jié)果整體呈現(xiàn)低估，流域B則擬合較好。

由表1的計算結(jié)果，從R2分析，流域A和流域B徑流模擬的R2均較高，兩者數(shù)值接近，擬合相關(guān)度較好。流域A在率定期的Ef為0.76，驗(yàn)證期為0.60，擬合效果出現(xiàn)較明顯下滑，按Ef評價標(biāo)準(zhǔn)，驗(yàn)證期的模擬結(jié)果為“不合格”；流域B在率定期的Ef為0.85，驗(yàn)證期為0.83，擬合效果較為穩(wěn)定，按Ef評價標(biāo)準(zhǔn)，驗(yàn)證期的模擬結(jié)果為“良”。

在基于SWAT模型的徑流模擬中，在相同的模擬次數(shù)和參數(shù)個數(shù)條件下，流域B檢證期的Ef達(dá)到了0.83，要顯著優(yōu)于流域A的Ef為0.60的模擬結(jié)果。結(jié)合研究區(qū)流域的實(shí)際地物分析，流域B的邊界，涵蓋了流域源頭的青龍山北麓，以及東和與介福雨量站；流域A的邊界則未包含上述地點(diǎn)。綜合上述分析，流域B的劃分優(yōu)于流域A，其結(jié)果更為接近研究區(qū)的真實(shí)情況。

4 問題和討論

① 對于包括巖溶區(qū)在內(nèi)的非閉合流域，流域邊界提取應(yīng)依據(jù)地下分水線。本次流域提取研究，只收集到了地表DEM數(shù)據(jù)，在使用D8法進(jìn)行流域提取時假定了地表分水線與地下分水線重合，因此提取出的流域邊界與真實(shí)情況存在一定誤差。

② 流域B面積為2 035.26 km2，流域提取基于水平方向30 m精度的DEM數(shù)據(jù)，流域邊界為30 m×30 m方格構(gòu)成，因?yàn)檫吔缟系姆礁裨谒秸`差范圍內(nèi)，這些方格可能存在也可能不存在，因此，根據(jù)邊界上方格的數(shù)目就可以計算出其對應(yīng)的面積誤差。經(jīng)統(tǒng)計，本次流域B的邊界共有8 007個方格，對應(yīng)的面積誤差為±7.21 km2，因此本次提取流域面積的真實(shí)值介于2 028.05 km2到2 042.47 km2之間。

5 結(jié)論

由于巖溶區(qū)地下暗河發(fā)育，基于原始DEM提取出的流域邊界會出現(xiàn)不同程度的缺失；把Burn-in法與D8法結(jié)合，對巖溶區(qū)的DEM進(jìn)行預(yù)處理，將地下河流“顯化”，結(jié)果使得流域邊界得到補(bǔ)充；將Burn-in法處理后的流域邊界利用SWAT模型進(jìn)行徑流模擬，模擬出率定期的Ef為0.85，驗(yàn)證期為0.83，其結(jié)果要顯著優(yōu)于Burn-in法處理前的流域邊界。研究結(jié)果表明，在巖溶區(qū)將Burn-in法與D8法結(jié)合提取出的流域邊界，比于單純利用D8法提取出的流域邊界更為接近流域的真實(shí)情況。