董吳欣，施召云，李 銘，鄭永勝，王 銳，張海龍，虞美秀，陳煉鋼

（1. 河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098； 2. 雅礱江流域水電開發(fā)有限公司，四川 成都 610051； 3. 中國(guó)電建集團(tuán)成都勘測(cè)設(shè)計(jì)研究院有限公司，四川 成都 610072； 4. 南京水利科學(xué)研究院 水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029）

0 引 言

數(shù)字高程模型 （Digital elevation model，DEM）［1］中包含流域的高程、坡度等信息，根據(jù)這些信息可以提取流域內(nèi)的水系，DEM 提取水系的方法有ArcGIS 水文模塊、Arc SWAT、Arc Hydro Tools 等。隨著地理信息系統(tǒng)（GIS）的快速發(fā)展，多位學(xué)者［2-4］提出參數(shù)優(yōu)化分析、影像與DEM 融合等方法以期提高流域水系提取精度。DEM 的分辨率是影響水系提取的重要因素［5，6］，不同分辨率的DEM 有不盡相同的地形特征信息會(huì)影響水系的提取。此外，集水面積閾值的選取影響水系提取質(zhì)量［7，8］，研究表明可采用河網(wǎng)密度法、均值變點(diǎn)分析法、水系分維法等確定最佳集水面積閾值。迄今為止，已有學(xué)者基于DEM數(shù)據(jù)對(duì)流域水系提取進(jìn)行研究，但多是針對(duì)特定區(qū)域特定的DEM數(shù)據(jù)；隨著DEM數(shù)據(jù)源的增多，針對(duì)何種DEM數(shù)據(jù)提取的水系精度更高且適用于不同的流域，少有學(xué)者開展相關(guān)的對(duì)比研究。

鑒此，本文針對(duì)山區(qū)流域，選擇長(zhǎng)江上游山區(qū)及珠江上游山區(qū)作為研究對(duì)象，基于ALOS 12.5 m DEM、ASTER 30 m GDEM 與SRTM 90 m DEM 采用Arcpy針對(duì)不同流域選用不同集水面積閾值集合提取河網(wǎng)并評(píng)價(jià)精度，而后通過(guò)Arc Hydro Tools 結(jié)合藍(lán)線河網(wǎng)修正DEM，從而遴選出高適性的DEM 數(shù)據(jù)開展山區(qū)流域河網(wǎng)的提取。

1 研究區(qū)域概況與數(shù)據(jù)來(lái)源

1.1 研究區(qū)域概況

長(zhǎng)江上游山區(qū)流域選擇了雅礱江力丘河流域、岷江雜谷腦河流域與金沙江黑水河流域。力丘河流域?qū)賮啛釒Ъ撅L(fēng)氣侯，發(fā)源于四川省甘孜藏族自治州康定市，平均比降10.18‰，平均海拔4 125 m，河長(zhǎng)194 km，流域面積5 892 km2。雜谷腦河流域?qū)賮啛釒Ъ撅L(fēng)氣侯，發(fā)源于四川省阿壩藏族羌族自治州鷓鴣山，平均比降18‰，平均海拔3 657 m，河長(zhǎng)168 km，流域面積4 174 km2。黑水河流域?qū)俳鹕辰蔁岷庸鹊貛?，發(fā)源于四川省涼山彝族自治州，平均比降13.5‰，平均海拔2 327 m，河長(zhǎng)192 km，流域面積3 618 km2。

珠江上游山區(qū)流域選擇漓江流域的大溶江水文站以上的河源山區(qū)及支流甘棠江。漓江流域?qū)賮啛釒Ъ撅L(fēng)氣候區(qū)，大溶江水文站以上的漓江上游建有斧子口水庫(kù)與川江水庫(kù)，平均海拔565 m，總流域面積728 km2。甘棠江為漓江支流，河長(zhǎng)60 km，平均海拔428 m，流域面積784 km2，甘棠江上游建有青獅潭水庫(kù)，為廣西最大的水庫(kù)。

1.2 數(shù)據(jù)來(lái)源

本文收集并分析了3 種DEM：ALOS PALSAR RTC HIGH RES 12.5 m DEM、 ASTER 30 m GDEM 和SRTM 90 m DEM數(shù)據(jù)。

ALOS （Advanced Land Observation Satellite） 衛(wèi)星發(fā)射的L波段微波可穿透植被表層，獲取全天候高精度的高程數(shù)據(jù)［9］。ALOS PALSAR RTC HIGH RES DEM 空間分辨率為12.5 m，坐標(biāo)系為UTM 投影，來(lái)源于美國(guó)國(guó)家航空航天局（https：//earthdata.nasa.gov/），其獲取時(shí)間為2006-2011年。

ASTER GDEM（Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model）是通過(guò)光學(xué)立體攝影測(cè)量獲取的數(shù)據(jù)產(chǎn)品［10］，空間分辨率為30 m，坐標(biāo)系為WGS84，來(lái)源于地理空間數(shù)據(jù)云（http：//www.gscloud.cn/），原始數(shù)據(jù)獲取時(shí)間為2000-2009年。

SRTM（Shuttle Radar Topography Mission）DEM 不受天氣現(xiàn)象的干擾，但對(duì)地表覆蓋物較為敏感［11］。SRTM DEM 空間分辨率為90 m，坐標(biāo)系為WGS84，由2000年航天飛機(jī)拍攝獲取，研究采用的數(shù)據(jù)來(lái)源于地理空間數(shù)據(jù)云（http：//www.gscloud.cn/）。

2 研究方法

2.1 河網(wǎng)提取

研究采用Arcpy批量處理研究區(qū)共15個(gè)DEM 柵格文件，執(zhí)行步驟見(jiàn)圖2，其原理如下：首先將裁剪過(guò)的DEM 進(jìn)行填洼和削峰，目的是移除數(shù)據(jù)中的異常點(diǎn)，避免出現(xiàn)河網(wǎng)斷裂的情況；而后使用D8算法，在3×3的網(wǎng)格上以最陡坡度確定水流方向并形成流向圖；流量則為既定水流方向匯入柵格的累積量。以最佳集水面積閾值提取水系信息并矢量化，劃分子流域，最終形成流域圖。

表1 研究數(shù)據(jù)基本信息Tab.1 Basic information of research data

圖1 研究區(qū)域概化圖Fig.1 Outline of the study area

圖2 基于Arcpy的河網(wǎng)提取Fig.2 Arcpy-based river network extraction

2.2 提取結(jié)果精度評(píng)價(jià)

河網(wǎng)套合差是數(shù)字河網(wǎng)與藍(lán)線河網(wǎng)疊加產(chǎn)生的細(xì)碎多邊形面積與流域總面積的比值。選用河網(wǎng)套合差評(píng)價(jià)河網(wǎng)提取精度，河網(wǎng)套合差（D）小，則代表提取的水系與流域真實(shí)水系吻合度高。

式中：D為河網(wǎng)套合差，%；Ai為細(xì)碎多邊形面積，km2；S為流域總面積，km2。

2.3 最佳集水面積閾值的確定

集水面積閾值影響數(shù)字河網(wǎng)提取的疏密程度，為減小主觀誤差，采用均值變點(diǎn)分析法確定最佳集水面積閾值。其分析步驟如下：針對(duì)n（n≥2）序列，變點(diǎn)Xi（i=1，2，…，n）將序列分為兩段，計(jì)算兩段序列的離差平方和并疊加得到Si，總樣本離差平方和S與Si的差值的最大值點(diǎn)為最佳集水面積。

2.4 基于Arc Hydro Tools的DEM修正

Arc Hydro Tools 模型基于GIS 平臺(tái)研發(fā)兼具地形處理與流域處理兩大功能，研究采用地形處理中的DEM 修正模塊（DEM Reconditioning）修正不同類型DEM。DEM 修正模塊采用AGREE 算法，其基本原理如下：①將輸入的藍(lán)線河網(wǎng)（.shp 文件）柵格化并進(jìn)行緩沖區(qū)分析；②對(duì)緩沖區(qū)柵格、臨近河網(wǎng)柵格與緩沖區(qū)外柵格三者按距離比例線性插值；③對(duì)河網(wǎng)所在柵格進(jìn)行降低高程處理，使河網(wǎng)所在柵格高程低于非河網(wǎng)柵格；④輸出修正后DEM。

3 結(jié)果與分析

3.1 多源DEM對(duì)河網(wǎng)提取的影響

3.1.1 河網(wǎng)精度分析

DEM 的分辨率是影響數(shù)字河網(wǎng)提取的一大重要因素［12］。研究選用河網(wǎng)套合差評(píng)價(jià)不同分辨率DEM提取效果，對(duì)分辨率為12.5、30 與90 m 的DEM 提取河網(wǎng)并與藍(lán)線河網(wǎng)疊加，疊加得出的細(xì)碎面積依據(jù)公式1 計(jì)算河網(wǎng)套合差，結(jié)果見(jiàn)表2?？梢?jiàn)，ALOS DEM 的細(xì)碎面積在本研究流域范圍內(nèi)最小，ASTER GDEM 最大，兩者在研究區(qū)的面積差分別為10.255、8.991、1.632、13.487、7.449 km2。河網(wǎng)套合差能更好地衡量提取的精度，研究區(qū)域ALOS DEM 的河網(wǎng)套合差比SRTM DEM 小0.121%、0.092%、0.106%、0.354%、0.070%，而SRTM DEM 的河網(wǎng)套合差又比ASTER GDEM 小0.053%、0.124%、0.009%、1.500%、0.881%，表明ALOS DEM 河網(wǎng)提取效果最好，SRTM DEM次之，ASTER GDEM 在研究區(qū)提取效果最差。

表2 不同分辨率DEM河網(wǎng)提取對(duì)比Tab.2 Comparison of DEM river network extraction at different resolutions

3.1.2 DEM修正后河網(wǎng)提取分析

不同分辨率DEM在地形起伏變化時(shí)存在一定偏差，故依據(jù)藍(lán)線河網(wǎng)修正DEM有一定的必要性。研究采用基于AGREE算法的Arc Hydro Tools 模型優(yōu)化提取精度，以河網(wǎng)矢量圖層為基線調(diào)整高程值并使用線性插值的方法生成高程緩沖。對(duì)修正后DEM 計(jì)算河網(wǎng)套合差（表3），可見(jiàn)，修正后的細(xì)碎多邊形面積與河網(wǎng)套合差明顯減小，表明修正后3 種DEM 提取的水系更加精確。研究區(qū)ALOS DEM 的河網(wǎng)套合差減小85.82%、85.50%、83.88%、85.37%、88.37%，其減小幅度最大且河網(wǎng)套合差值最小，修正后ALOS DEM 河網(wǎng)套合差為0.054%、0.040%、0.050%、0.128%與0.160%，表征其精度最高。ASTER GDEM 與SRTM DEM 在不同流域的修正效果亦不盡相同，力丘河流域ASTER GDEM 修正后僅減小26.57%，修正效果在該流域最差，但ASTER GDEM 在其他研究區(qū)比SRTM DEM 減小幅度大17.31%、37.01%、18.39%、10.89%，修正效果更優(yōu)。由此導(dǎo)致DEM 修正后ASTER GDEM 與SRTM DEM 提取河網(wǎng)的精度在不同的流域表現(xiàn)不一，但總體上ASTER GDEM 表現(xiàn)更佳。因此DEM 修正后，ALOS DEM 河網(wǎng)提取效果依然最好，ASTER GDEM次之，而SRTM DEM在研究區(qū)提取效果變?yōu)樽畈睢?/p>

表3 修正后不同精度DEM提取效果對(duì)比Tab.3 Comparison of DEM extraction effects with different precisions after correction

3.2 最佳集水閾值的確定

通過(guò)對(duì)數(shù)據(jù)散點(diǎn)的擬合，發(fā)現(xiàn)集水面積閾值與河網(wǎng)密度的關(guān)系曲線起初隨著閾值的增大河網(wǎng)密度急劇下降而后趨于平緩，較為符合冪函數(shù)特征，并且已有學(xué)者研究表明冪函數(shù)擬合的效果最佳［8］，故通過(guò)設(shè)置多個(gè)閾值擬合y=kxa（式中：x為集水面積閾值，y為數(shù)字河網(wǎng)密度），擬合公式與決定系數(shù)R2見(jiàn)圖3，其中漓江上游（大溶江水文站以上的漓江流域）的SRTM DEM 的決定系數(shù)R2為0.930 2，雜谷腦河流域ALOS DEM 的R2為0.985 4，其余工況的R2均大于0.99，表征其相關(guān)性較好。

圖3 河網(wǎng)密度與集水面積閾值擬合對(duì)比Fig.3 Comparison of river network density and catchment area threshold

考慮流域面積的差異性，長(zhǎng)江上游山區(qū)流域同種DEM選用相同閾值序列，珠江上游山區(qū)則采用另一序列，具體閾值見(jiàn)表4。最佳集水面積閾值的選取通常采用冪函數(shù)的拐點(diǎn)來(lái)確定，但擬合曲線拐點(diǎn)較難尋找且受主觀因素影響較大，為減小誤差，研究采用均值變點(diǎn)分析法選取流域的最佳集水面積閾值，結(jié)果見(jiàn)圖4。可見(jiàn)，總樣本的離差平方和S與變點(diǎn)兩側(cè)離差平方和之和Si的差值隨著閾值的增加呈先漲后落的趨勢(shì)，故各曲線均存在一個(gè)峰值點(diǎn)，即最佳集水面積閾值點(diǎn)。圖4（a）中ALOS DEM 在長(zhǎng)江上游（力丘河流域、雜谷腦河流域、黑水河流域）的峰值為第6 個(gè)點(diǎn)，珠江上游（漓江上游流域、甘棠江流域）的峰值為第8 個(gè)點(diǎn)；圖4（b）中研究區(qū)ASTER GDEM 的峰值點(diǎn)序列號(hào)均為3；圖4（c）中研究區(qū)SRTM DEM 的峰值點(diǎn)序列號(hào)均為5。研究區(qū)不同DEM 起伏有明顯差異，但峰值點(diǎn)的集水面積閾值一致，即不同流域同種類型的DEM有相同的最佳集水面積閾值。將均值變點(diǎn)分析所獲取的最佳集水面積閾值與圖3數(shù)據(jù)對(duì)比，得出三類DEM 最佳集水面積閾值對(duì)應(yīng)的河網(wǎng)密度，發(fā)現(xiàn)所得點(diǎn)數(shù)據(jù)位于河網(wǎng)密度的陡降與平緩的過(guò)渡帶。采用ALOS DEM 提取河網(wǎng)的最佳集水面積閾值為20 000，ASTER GDEM 與SRTM DEM 提取河網(wǎng)的最佳集水面積閾值則為500。

表4 各類DEM在不同流域的積水面積閾值序列Tab.4 Threshold sequences of various types of DEMs in different watersheds

表5 DEM修正前后不同地形的面積響應(yīng)Tab.5 DEM corrects the area response of different terrains before and after

圖4 S與Si差值變化Fig.4 Changes in the difference between S and Si

3.3 地形特征對(duì)數(shù)字河網(wǎng)提取的影響

地形地貌對(duì)數(shù)字河網(wǎng)的提取亦有影響，根據(jù)梁麗芳等［13］基于DEM 的地形分類研究，將山區(qū)按坡度分為小起伏山地（2～6°）、中起伏山地（6～25°）與大起伏山地（＞25°）。按不同地貌分別對(duì)原始DEM 與修正后DEM 所提取的面積進(jìn)行對(duì)比，發(fā)現(xiàn)修正后大起伏山地面積均增加，而小起伏山地面積均減小。力丘河流域與甘棠江流域的ALOS DEM 修正后中起伏山地面積增大，其余工況中起伏山地的面積均減小。大起伏山地由于坡度大且河網(wǎng)兩側(cè)的山體坡度也較大，因此在匯流累積量的計(jì)算過(guò)程中，DEM 各格網(wǎng)累積量數(shù)值不易與周圍的格網(wǎng)重復(fù)，提取的水系更加精確［14］。從不同地形面積方面考慮，DEM 修正將部分小起伏山地轉(zhuǎn)化為大起伏山地，修正后DEM的數(shù)字河網(wǎng)提取精度提高。ALOS DEM 在大起伏山地提取時(shí)相對(duì)誤差最小，分別為1.247%、0.869%、1.249%、2.808%、2.334%，故認(rèn)為其精度最高。ASTER GDEM 在力丘河流域修正后面積差較大，對(duì)比SRTM DEM 相對(duì)誤差大3.919%，在其他研究區(qū)的表現(xiàn)則更優(yōu)，兩者相對(duì)誤差之差為0.941%、5.741%、26.946%、22.632%。

為精確化分析地形地貌對(duì)數(shù)字河網(wǎng)提取的影響，采用套合差分析研究區(qū)各類地形的河網(wǎng)提取情況。根據(jù)表6 套合差值，可以得出在同一流域內(nèi)大起伏山地的套合差均小于中起伏山地與小起伏山地，即大起伏山地河網(wǎng)提取效果最佳。中起伏山地提取效果也優(yōu)于小起伏山地，故印證地形起伏越大越有利于河網(wǎng)的提取。

表6 地形特征對(duì)多源DEM數(shù)字河網(wǎng)提取對(duì)比Tab.6 Comparison of topographic features on multi-source DEM digital river network

4 結(jié) 論

本文基于Arcpy 與Arc Hydro Tools 對(duì)山區(qū)流域多源DEM 提取河網(wǎng)并修正分析。修正前后ALOS DEM 提取效果均最優(yōu)且適用于不同流域，未修正前的SRTM 90m DEM 優(yōu)于ASTER 30 m GDEM，修正后力丘河流域與漓江上游流域SRTM 90 m DEM仍優(yōu)于ASTER 30 m GDEM，而雜谷腦河流域、黑水河流域與甘棠江流域采用修正后ASTER 30 m GDEM 提取精度較高。通過(guò)河網(wǎng)密度與均值變點(diǎn)分析法確定最佳集水面積閾值，綜合考慮河網(wǎng)套合差與地形面積變化并得出大起伏山地利于河網(wǎng)的提取。研究采用的均值變點(diǎn)分析法、Arc Hydro Tools 修正、河網(wǎng)提取精度評(píng)價(jià)方法、地形特征影響分析可應(yīng)用于相似流域DEM的選取。