基于Arcgis-python斷面匯水面積批量提取方法研究

洪明海，汪仕偉，曾文治，黃介生，楊榮芳，蘇海鵬，李析男

(1. 貴州省水利水電勘測設計研究院，貴陽 550002；2. 武漢大學水資源與水電工程科學國家重點實驗室，武漢 430072)

0 引 言

匯水面積是由分水線包圍的集水面積，在閉合的匯水面積內(nèi)，構(gòu)成一定形狀的水路網(wǎng)系統(tǒng)[1]，是設計橋梁涵洞的孔徑大小、水庫大壩的壩高以及水庫來水量等的重要參數(shù)。因此，準確地勾繪出設計斷面的匯水面積對于水庫工程設計、生態(tài)環(huán)境治理以及水文統(tǒng)計分析提供數(shù)據(jù)尤為重要[2, 3]。

在地理信息系統(tǒng)軟件出現(xiàn)之前，匯水邊界主要以地形圖上的分水線進行劃分[4]，在地形圖上依據(jù)等高線分布繪制流域的分界線，然后量算出該分界線所包圍的面積，效率和精度依賴于人員操作熟練度；同時，對于地形復雜，山脊線和分水嶺等分界線難以區(qū)分，人工勾繪極易出錯進而嚴重影響匯水面積的成果，盧洋暘等[5]的研究也表明人工勾繪匯水面積的方法在處理復雜山地地形和范圍較大地形時存在較大的局限性；此外，人工勾繪的方法對于河道斷面位置無相應的地形圖或者已有地形圖的精度難以滿足工程計算的需求時顯得力不從心，李曉等[6]也認為需要一種新的方法應用于無地形圖地區(qū)的工程位置處流域面積的提取和計算。而現(xiàn)在，計算機程序可以從DEM中提取地形要素特征來生成匯水邊界，應用計算機技術，只需要傳統(tǒng)方法的一小部分時間就可以得到匯水區(qū)域的邊界[7]。目前，Arcgis、SWAT、MIKE、Mapgis、Civil 3D等眾多軟件都可以實現(xiàn)基于DEM的水文研究工作，比如匯水面積的勾繪和計算、提取河網(wǎng)信息以及流域特征分析等[8-12]?；贒EM分析計算的方法不依賴地形圖，不僅節(jié)約時間成本，而且受人為因素影響較小，且在地形復雜地區(qū)和大范圍地區(qū)均具有較好的適用性，但目前該方法多是針對單個斷面的一系列復雜操作，自動化程度還不高，對于短時間內(nèi)需要同時處里大量的斷面匯水面積并獲得數(shù)據(jù)結(jié)果仍是一個不小的挑戰(zhàn)。因此，有必要對該方法做進一步的拓展和研究以滿足工程實踐的需要。

本文擬采用ESRI公司的Arcgis10.4軟件來處理芙蓉江流域內(nèi)河道斷面的匯水面積，同時通過python調(diào)用Arcigs中相應的水文工具形成匯水面積處理模型，批量獲取河道斷面的匯水面積，將計算的結(jié)果和已有的數(shù)據(jù)進行比較分析，同時將該方法應用到洪渡河流域、桐梓河流域、綦江流域、赤水河流域、大渡口以上區(qū)域、大渡口以下干流區(qū)域以及六盤水區(qū)域等7個研究區(qū)進行驗證，為河道斷面匯水面積的勾繪提供一種新的借鑒方法。

1 材料與方法

1.1 研究區(qū)概況和數(shù)據(jù)來源

1.1.1 研究區(qū)概況

本次計算的斷面均位于貴州省遵義市的芙蓉江流域，根據(jù)《中國河湖大典·長江卷》[13]和《貴州河湖》記載，芙蓉江(107°02′E～107°52′E，28°08′N～29°20′N)發(fā)源于貴州省遵義市綏陽縣枧壩鎮(zhèn)，由西南向東北流經(jīng)貴州省遵義市正安縣、道真縣和重慶市武隆區(qū)、彭水縣，于武隆區(qū)江口鎮(zhèn)匯入烏江，全長234 km，流域面積7 406 km2，西鄰綦江水系的桐梓河、松坎河，北靠烏江，南鄰湄江，東鄰洪渡河；南北長約160 km，東西寬約195 km，流域地勢起伏大，東低西高，上游地形較開闊平坦，沿河發(fā)育有河谷階地和河漫灘，中下游河谷深切，河道狹窄，灘多流急，巖溶發(fā)育，多洞穴、井泉和伏流暗河段，屬中山峽谷及巖溶地貌特征，平均高程約850 m(黃?；?，下同)，總落差1 124 m，河流平均比降5.76‰，流域面積大于1 000 km2的支流有清溪河、三江和梅江，流域面積大于100 km2的支流有新民河、關壩等15條河流。流域?qū)賮啛釒駶櫦撅L氣候區(qū)，年平均氣溫16.4 ℃，年無霜期281.8 d，年日照時數(shù)為1 057 h，多年平均降水量1 106 mm，降水多集中在4-10月，占年降水量的73%，年水面蒸發(fā)量為603.2 mm。河口多年平均流量為169 m3/s，多年平均懸移質(zhì)輸沙量為277.3 萬t。流域多年平均水資源總量53.3 億m3，干流水力資源理論蘊藏量47.3 萬kW。流域內(nèi)無大型工礦企業(yè)，工業(yè)污水排放量小，污染物主要來自農(nóng)田的農(nóng)藥、化肥、養(yǎng)殖業(yè)的污染負荷和生活污水。

1.1.2 數(shù)據(jù)來源

本文使用的dem數(shù)據(jù)為Aster傳感器采集的30 m分辨率的柵格數(shù)據(jù)，WGS84橢球，數(shù)據(jù)下載于地理空間數(shù)據(jù)云(http:∥www.gscloud.cn/search)，由中國科學院計算機網(wǎng)絡信息中心提供。河流矢量數(shù)據(jù)來源于各水務局搜集資料整理，已經(jīng)項目三級校審制度審核。

根據(jù)已有的資料整理，本文分析選取的河道斷面均為貴州省遵義市境內(nèi)的芙蓉江流域上50 km2以上支流的斷面，共有38個，分三類，其中，一類是支流匯口斷面，如梅江匯口斷面、羅家河匯口斷面，共有29個，占比76.32%；第二類是已建成水庫斷面，以壩址為控制斷面，如沙壩水庫斷面、慶埡電站水庫斷面，共有6個，占總數(shù)的15.79%；第三類是河流上的水文站控制斷面，有3個，占總數(shù)的7.89%，它們是旺草水文站斷面、長壩水文站斷面和五家院子水文站斷面，從圖1可以看到本次分析選取的三種類型共38個河道斷面的分布情況，其中水文站類型的3個斷面有在芙蓉江干流上游的(旺草水文站斷面)，有在芙蓉江干流下游的(長壩水文站斷面)，還有位于芙蓉江一級支流梅江上的(五家院子水文站斷面)；6個水庫型斷面則分布在貴州省境內(nèi)的芙蓉江最大的三條一級支流清溪河、三江河梅江上；支流匯口型斷面則廣泛分布于芙蓉江流域的支流上，因此，本次計算分析選取的河道斷面具有較強的代表性。

研究區(qū)域的河流水系情況和本次分析選取的38個河道斷面分布情況如圖 1所示。

38個斷面的匯水面積數(shù)據(jù)均來源于水利普查資料，如表 2中“已有資料面積”欄所示，其匯水面積最大值為長壩水文站，為5 454 km2，最小值為鐵釘巖水庫，匯水面積為7 km2。本文根據(jù)中華人民共和國國家標準 GB/T 4883-2008(數(shù)據(jù)的統(tǒng)計處理和解釋·正態(tài)樣本離群值的判斷和處理)[14]，采用標準推薦的在未知總體標準差的情形下的離群值處理方法—格拉布斯(Grubbs)檢驗法對已有資料的流域面積進行檢驗，使用雙側(cè)檢驗，其計算公式如下所示：

圖1 研究區(qū)域河流水系及河道斷面分布情況Fig.1 Distribution of river systems and river sections of study area

(1)

(2)

(3)

(4)

表1 格拉布斯雙側(cè)檢驗計算結(jié)果表Tab.1 Table of calculation results of grubbs bilateral test

此外，本文提出的方法還在洪渡河流域、桐梓河流域、綦江流域、赤水河流域、大渡口以上區(qū)域、大渡口以下干流區(qū)域以及六盤水區(qū)域等7個研究區(qū)進行驗證，該7個研究區(qū)的匯水面積數(shù)據(jù)也包含水文站、水庫和匯口等3種類型的斷面數(shù)據(jù)，且所有7個研究區(qū)的數(shù)據(jù)均采用格拉布斯雙側(cè)檢驗法進行檢驗，保證對比資料在統(tǒng)計學上真實有效。

1.2 計算原理

基于DEM流域劃分和河流水系的提取一般都會對DEM柵格數(shù)據(jù)進行填洼處理以及計算流向，此外，當獲取的外部矢量數(shù)據(jù)(比如水文站等斷面數(shù)據(jù))和DEM不是來自同一份數(shù)據(jù)的話，矢量數(shù)據(jù)和DEM之間就有產(chǎn)生一定的誤差，這時就需要使用捕捉傾瀉點功能或者編輯功能對其適當修正。

洼地(或突起)是由于數(shù)據(jù)的分辨率或者將高程四舍五入到中心像元最近的整數(shù)值而產(chǎn)生的常見錯誤，應對洼地進行填充以確保盆地和河流得以正確劃界，如果未對洼地進行填充，則生成的河流水系網(wǎng)絡可能會呈現(xiàn)不連續(xù)性。填洼工具使用與焦點流、流向、洼地、分水嶺和區(qū)域填充等工具等效的功能來定位和填充洼地，該工具的執(zhí)行過程會進行迭代，直到指定 z限制內(nèi)的所有洼地均填充完畢，在填充洼地的同時，可能會在填充區(qū)域的邊界處創(chuàng)建其他洼地，這些洼地將在下個迭代中移除[15]，其原理示意圖如圖2。

圖2 填洼原理示意圖Fig.2 Schematic diagram of fill principle

圖3 流向原理示意圖Fig.3 Schematic diagram of flow direction

1.3 研究技術路線及步驟

本文的研究技術路線如圖 4所示，其主要操作步驟如下：

(1)對研究區(qū)DEM進行裁剪投影處理。

(2)根據(jù)DEM提取研究區(qū)的河流水系，其中填洼處理和計算流向以及設置流量閾值生成河流是比較關鍵的步驟，填洼影響流向的形成，計算流向決定了流量的匯集方向，而流量閾值的設置直接決定著河流水系的密度。

(3)對37個河道斷面數(shù)據(jù)進行處理分析，一般導入Arcgis的矢量河道斷面數(shù)據(jù)和所使用的DEM數(shù)據(jù)或多或少都會有一定的偏差，這時需要使用捕捉傾瀉點工具或者編輯工具并參照已有的矢量河流來適當調(diào)整以滿足軟件的計算需要。

(4)以上3步都處理好之后就可以使用分水嶺工具開始計算斷面的匯水區(qū)域了，注意斷面匯水區(qū)域之間的重疊情況，這里使用python批量分割出每個斷面單獨計算之后再合并到一起，避免了上游斷面的匯水區(qū)域未包含在下游的匯水區(qū)域中的特殊情況，此外計算得到的斷面匯水區(qū)域是柵格，需要將其轉(zhuǎn)為矢量格式。

(5)計算矢量匯水區(qū)域的各種屬性值，這里主要計算38個斷面的匯水面積值，需要在投影坐標系下方可計算。

圖4 本文研究技術路線圖Fig.4 Technology roadmap

1.4 分析指標

評價指標選擇絕對誤差Δ、相對誤差m、決定系數(shù)R2、均方根誤差RMSE、相對均方根誤差RRMSE、相對分析誤差RPD、偏差Bias：

Δ=Xgis-Xdata

(5)

(6)

(7)

R2=R×R

(8)

(9)

(10)

(11)

(12)

2 結(jié)果分析

2.1 匯水面積分析

表2給出的是python調(diào)用Arcgis工具批量計算的芙蓉江流域內(nèi)37個河道斷面匯水面積(下稱計算面積)和已有資料的匯水面積(下稱已有面積)結(jié)果，其中，按照公式(5)和公式(6)分別計算了計算面積和已有面積的絕對誤差和相對誤差。表2中，旺草水文站(8號)位于芙蓉江干流上，但其位于上游，從圖1中看到其上只有羅家河一條支流匯入，其匯水面積小得多(計算面積為401.83 km2，已有面積為413 km2)(表 2)，相對誤差為2.7%；此外，我們在芙蓉江支流梅江上也選取了一個水文站控制斷面—五家院子水文站斷面(圖 1)，其計算面積為1 212.14 km2，已有面積為1 230 km2，相對誤差為1.45%(表 2)；2個水文站控制斷面中，絕對誤差最大的是五家院子水文站控制斷面，為17.86 km2，相對誤差最大的則是旺草水文站控制斷面，為2.7%(表2)。

而6個水庫型斷面均位于芙蓉江的支流上，貴州省境內(nèi)芙蓉江最大的三條一級支流清溪河(清溪水庫斷面)、三江(慶埡電站水庫斷面)和梅江(洋渡電站水庫斷面和大河壩水庫斷面)均有分布(圖 1)，其中鐵釘巖水庫斷面位于巴漁河上(伏流型河流)。6個水庫型斷面中，匯水面積最大的是清溪水庫斷面，計算面積為1 237.44 km2，已有面積為1 261 km2(表 2)；絕對誤差最小的巴漁河上的鐵釘巖水庫斷面，為0.54 km2(表2)；相對誤差最小的為洋渡電站水庫斷面，為0.79%，最大的為鐵釘巖水庫斷面，為7.66%。

本次計算分析中最多的是支流型匯口斷面，共有29個，占比76.32%，廣泛分布在芙蓉江的各支流上。其中匯水面積最小的是芙蓉江一級支流茨竹河匯口斷面(圖 1)，計算面積為48.64 km2，已有面積為50 km2(表 2)；而最大的則是清溪河匯口斷面，計算面積為1 493.35 km2，已有面積為1 504 km2(表 2)。29個支流型匯口斷面中，絕對誤差最小的是池武溪匯口斷面，為0.08 km2，相對誤差最小的是梅江匯口斷面，為0.09%(表2)。王崢等[18]人采用DEM對涇河流域的1 000 km2以上部分支流的匯水面積進行提取并做了比較分析，絕對誤差也達到17 km2。其余河道斷面對應的匯水面積以及誤差等數(shù)據(jù)請參見表 2。

2.2 計算方法的驗證和誤差分析

37個河道斷面的匯水面積數(shù)據(jù)的相對誤差中，最小的為梅江匯口斷面，為0.09%，最大的則是鐵釘巖水庫斷面，為7.66%。圖5給出的是已有資料的匯水面積和python批量計算的匯水面積相關性曲線圖，從圖中可以看到基本上所有的數(shù)據(jù)點均落在1∶1線上，計算表明其決定系數(shù)R2為0.99，而決定系數(shù)是用來衡量理論值與計算值的離散程度，其值范圍0～1，越接近于1，效果越好，因此通過決定系數(shù)來看通過python調(diào)用Arcgis工具批量處理的方法精度很高。此外，本文也計算了37個數(shù)據(jù)點的均方根誤差RRMSE、相對分析誤差RPD以及偏差Bias，相對均方根誤差是用來衡量殘差的相對大小，該值越小，說明計算結(jié)果越好[19]，而數(shù)據(jù)點計算得到的相對均方根誤差僅為2.74%，表明計算結(jié)果很好；相對分析誤差大于2，說明計算結(jié)果很好，相對分析誤差小于1.4，說明計算較差，相對分析誤差介于1.4～2，說明計算較可靠，還可以進一步調(diào)整[20]，37個數(shù)據(jù)點的相對分析誤差為49.01，遠大于2，表明計算結(jié)果很好；37個河道斷面的匯水面積數(shù)據(jù)的偏差為2.88，偏差為理論值與計算值的偏差程度，大的偏差值表明計算值沒有很好地逼近理論值，計算效果差。

表2 已有資料和python批量計算的斷面匯水面積誤差表Tab.2 Table of error of available data and python-calculation

注：類型欄1代表水文站，2代表水庫，3代表匯口。

圖5 已有資料和python批量計算的匯水面積對比Fig.5 Comparison of catchment area of available data and python-calculation

同時，本文還將該方法應用其他研究區(qū)域以進一步驗證通過python調(diào)用Arcgis批量計算匯水面積的可靠性。除本文建立方法使用的芙蓉江流域以外，還在洪渡河流域、桐梓河流域、綦江流域、赤水河流域、大渡口以上區(qū)域、大渡口以下干流區(qū)域以及六盤水區(qū)域等7個研究區(qū)廣泛選取河道控制斷面來驗證該方法，選取的斷面類型包含匯口、水庫和水文站等3類。同樣，本文也采用格拉布斯雙側(cè)檢驗法對新選取的研究區(qū)已有的匯水面積資料進行檢驗，排除統(tǒng)計離群值的存在，保證用于比較分析的匯水面積在統(tǒng)計學上真實有效。之后，根據(jù)選取的斷面數(shù)據(jù)采用本文提出的方法批量計算了匯水面積數(shù)據(jù)，然后計算了平均相對誤差(m)、決定系數(shù)(R2)、相對均方根誤差(RRMSE)、對分析誤差(RPD)以及偏差(Bias)列于表3。從表3中可以看出，相對誤差平均值均較小，最大的僅為7.76%，而決定系數(shù)均在0.95以上，表明計算得到的匯水面積數(shù)據(jù)和已有的匯水面積數(shù)據(jù)具有很高的相關性，同時相對均方根誤差也很小，最小的為桐梓河流域2.1%，7個驗證研究區(qū)的相對分析誤差均大于2，說明該方法計算結(jié)果好，計算方法合理，可用于對斷面匯水面積的計算。其余統(tǒng)計數(shù)據(jù)請參見表3。

盡管本次選取的37個河道斷面的匯水面積和已有資料的匯水面積數(shù)據(jù)對比結(jié)果很好，相對誤差小，決定系數(shù)高，但在模型搭建和運行處理過程中也發(fā)現(xiàn)了許多需要進一步完善和改進的地方，誤差產(chǎn)生來源主要有以下4點，列于表4：第一，所使用的DEM產(chǎn)生的誤差。本文采用的DEM分辨率為30 m×30 m，和真實的現(xiàn)實地形相比，在某個像元處可能其地形會被概化為一個值，會產(chǎn)生不小的誤差，如果條件允許，可以考慮10 m×10 m、1 m×1 m甚至更高分辨率的DEM數(shù)據(jù)，但這樣也會增加模型的運行時間，同時還會對計算機等硬件性能有更高的需求。第二，對比值產(chǎn)生的誤差。本文所分析的兩組河道斷面的匯水面積數(shù)據(jù)，一組來源于我們通過python調(diào)用Arcgis計算得出，而另一組對應河道斷面的匯水面積數(shù)據(jù)則來源于水利普查等相關資料，其源頭是人工在地形圖上勾繪計算得到，這種就會存在人為誤差，在小的匯水區(qū)域可能誤差較小，但是在匯水面積很大，比如超過1 000 km2的斷面，河流水系錯綜復雜，地形圖上的山脊線和山谷線不易分辨，這樣就可能會產(chǎn)生較大的誤差。第三，Arcgis本身計算原理的誤差甚至可能是錯誤。在提取河道斷面的匯水面積過程中，其中有兩個步驟是計算流向和填洼，關于其計算原理在前面有所介紹。填洼是如果某中心像元值均低于其周圍的像元值，那么軟件會認為這個中心像元值是一個錯誤值，從而會對該像元值進行填洼，使其像元值等于周圍值的最小值，但這有時候并不符合現(xiàn)實情況，比如說自然或者人工形成的洼地，這是客觀存在的，但是軟件會將其“補平”；另外一個步驟是計算流向，Arcgis雖然提供了三種不同的流向計算方法，但是默認的(或者當您打開計算流向工具時該方法即為D8，不可更改)算法是D8算法，該算法認為任意一個像元值僅存在一個流向，從而會計算出某一點的唯一流向，但是現(xiàn)實中，某點的水流可能會有兩個流向，只是每個流向分配的流量不一樣而已，因此此種算法在某些情況下可能會造成很大的錯誤。第四，工具參數(shù)引起的誤差。在使用Arcgis自帶的工具“捕捉傾瀉點”時，對于搜索距離的控制也是相當?shù)年P鍵，必要的時候還需要打開相應的矢量河流數(shù)據(jù)人工校對軟件捕捉的結(jié)果是否準確，尤其是在河流水系復雜的支流匯口斷面，張瀟戈[21]等的研究也表明用DEM提取的河網(wǎng)匯口和真實河網(wǎng)匯口存在一定的偏差。

表3 計算方法在其他研究區(qū)的驗證Tab.3 Validation of python-calculation in other study area

表4 研究方法誤差來源分析Tab.4 Analysis of error sources of python-calculation

3 結(jié) 語

本文通過python調(diào)用Arcgis相應的工具來批量計算芙蓉江流域50 km2以上河流的河道斷面的匯水面積，河道斷面類型包含水庫、支流匯口以及水文站等3種控制斷面，分布廣泛，代表性強，其主要結(jié)論如下。

(1)計算結(jié)果準確性高。計算的河道斷面的匯水面積絕對誤差Δ中最小的僅為0.08 km2，相對誤差m最小的僅為0.09%，決定系數(shù)R2高達0.99，相對均方根誤差RRMSE僅有2.74%，相對分析誤差RPD為49.01，遠大于2，偏差Bias為2.88，表明本方法不管匯水面積大小都能獲得很好的計算結(jié)果。

(2)計算速度快，可極大地節(jié)省生產(chǎn)時間。本文計算的河道斷面均位于貴州省遵義市境內(nèi)的芙蓉江流域，根據(jù)已有的數(shù)據(jù)資料，為了方便計算結(jié)果的驗證和比較，本文選取的河道斷面共有37個，通過計算機計算匯水面積每分鐘可以計算7個，大約需要6 min，而通過人工在地形圖上勾繪，一個斷面勾繪匯水面積大約需要25 min，38個斷面共需要950 min，可見，通過python批量處理的時間僅為人工的0.63%，這在生產(chǎn)實踐中極大地節(jié)約工作時間，提高工作效率，贏得時間成本。

(3)為河流控制斷面的匯水面積勾繪提供一種新的計算方法。采用python調(diào)用Arcgis批量計算匯水面積，不僅可以極大地提高工作效率，還可以方便快速地查看和計算其他特征參數(shù)，為控制斷面流量的計算和壩址庫容曲線的計算分析奠定基礎。