基于雙側均值變點法的數字河網閾值劃定
——以南苕溪流域為例

2022-06-23 10:11:50聶啟陽

綠色科技 2022年10期

聶啟陽

(上?？睖y設計研究院有限公司，上海 200335)

1 引言

近年來，隨著衛星遙感、無人機傾斜攝影等測量手段的逐漸豐富，數字高程模型(Digital Elevation Model，DEM)為水文、生態環保等學科研究提供了極其廣闊的運用場景[1]。DEM數據在使用各種流向算法后可生成表面流向柵格數據，再使用匯流運算可獲得表征每個網格匯流數量的匯流累積柵格數據，目前已有HydroSHEDS、GDBD、JapanDir等分析產品[2～4]。在匯流累積柵格數據的基礎上，往往使用一個集水面積閾值(CSA)來劃定流域內的河道區域與非河道區域，因此閾值選取的大小將直接影像河網提取的細致程度和后期研究的適用程度[5]。目前閾值選擇方法的主要基于2個方向：①根據已有藍線數據進行對比，尋求與藍線形態擬合最佳的閾值；②根據流域閾值-河網密度曲線，選取變化趨于平緩的點對應閾值最為最佳閾值[6～8]。第一種方式受制于藍線數據的精度及獲取難度，在小區域或無資料地區難以運用；后者應用簡單，且不受制于其他資料，自身基于地形演化原理具有較高的理論支撐，已有眾多研究使用均值變點法檢測這一平穩點的位置，用于改進其主觀性。然而不同研究對象的實際運用研究中，汛期和非汛期河網具有較大差異，已有研究表明流域河網存在多個閾值，主要可分為：①曲線由急劇驟減到逐漸減少的“閾值點”，表征坡面上河網鏈消失的閾值分界點；②曲線由逐漸減少到逐漸平穩的點，用于表征小型溝道上河網鏈消失的閾值分界點[9]。本次研究基于雙側均值變點法，在南苕溪流域開展數字河網提取驗證，尋求一種客觀有效的河網雙閾值劃定方式，以期為后續水文、生態環境模型研究的數據前期處理提供參考。

2 研究區域及數據

2.1 研究區域概況

南苕溪流域位于浙江省杭州市，流域面積約720 km2，全長約63 km。流域地處中亞熱帶季風氣候區南緣，年雨量充沛多年平均降雨量為1460 mm，降水集中在4～9月份，多年平均氣溫15.8 ℃[10]。

南苕溪自東天目山水竹塢發源，并由浪口溪、南溪和錦溪3條主要支流組成，浪口溪南流經里畈水庫至浪口匯入南溪后稱南苕溪，東流穿過臨安城區于青龍口匯入錦溪后入青山水庫。浪口溪和南溪河段地處浙西山丘區，河道坡降較大，土地利用類型也存在差異，浪口溪以天然林覆被為主，南溪沿岸雷竹林廣泛分布。錦溪河段地勢相對平坦，依次流經上游農業區、玲瓏工業區和臨安市區[11]。南苕溪在匯入東苕溪后，最終匯入太湖(圖1)。

圖1 研究區域地理地貌

2.2 數據來源及處理

本研究使用的DEM數據來源為ASTER GDEM V3，此數據集由日本經濟產業省(METI)和美國國家航空航天局(NASA)在2019年8月聯合發布。數據垂直分辨率為1 m，水平分辨率為30 m。具有較高的準確度，已得到國內外的廣泛運用[12～15]。(https://asterweb.jpl.nasa.gov/gdem.asp)。

DEM后續處理采用一種優先洪水算法改進的D8流向算法[16]。該算法無須對DEM進行前期填洼，因此不會改變DEM信息，可解決過度填洼產生的平行河網問題，對于本研究區下游較為平坦的城市區域較為合適[17]。

該算法以優先洪水算法的漫水思想為基礎，對DEM虛擬出由外至內的淹水過程，從外側最低網格淹沒開始至DEM最高網格被淹沒為止，反向記錄此流程就形成了整個地形內水流的整體流向趨勢。二維示意圖詳見圖2。圖2a反應的是漫水過程，外圍漫水過程從右端f開始至e，隨后由a至b，b至c，最終e和c至d完成整個柵格的淹沒過程。反向記錄整個流程即可得到流向柵格數據，同時因為是單流向算法，因此結合D8的最陡下降原則判定具有多流向趨勢的柵格，其中d具有流向c和e的趨勢，根據最陡下降原則確定d流向e。因此整個流向方向為圖2b：d流向e再流向f，c流向b再流向a，最終流出整個流域。整體來看對于外流型流域，該流向算法識別性強，有較強的物理過程支撐，無需填洼，能有效地解決填洼導致的平行河網問題。

DEM經過處理后得到地表流向柵格數據，再通過匯流累積算法獲得匯流累積柵格數據。常用的匯流累積算法采用一套遞歸：首先初始化匯流累積柵格數據均為1，對每一柵格進行判定，若其匯流累積值為1則根據地表流向柵格數據判定其四周8個網格流向它的柵格，求和滿足要求的柵格匯流累積值，在加上該柵格的初始值1。遍歷所有柵格，獲得整個區域內的匯流累積柵格數據，因為分水嶺僅流出，周圍無網格流向它，因此其匯流累積值為0，因此再對整個區域匯流累積柵格數據-1，得到最終的匯流累積柵格數據。

3 研究方法

3.1 雙側均值變點法

雙側均值變點法為均值變點法的變式，均屬于數理統計方法。針對河網閾值問題其離散化模式為以下狀態：

閾值序列為a，長度為n，其值為a1、a2、…、an-1、an，對應計算的河網密度序列為b，長度為n，其值為b1、b2、…、bn-1、bn。

圖2 水流整體流向趨勢二維示意

河網密度計算使用：

(1)

式(1)中，Ln為閾值序列對應的河網長度，A為流域面積。

令m=2、3、…、n-4、n-3，k=m+2、m+3、…、n-2、n-1，m和k將序列河網密度n分割為3段。

計算各段內的算術平均值：

(2)

(3)

(4)

進而求得各段統計量：

(5)

(6)

(7)

求和各段統計量：

Smk=SⅠ+SⅡ+SⅢ

(8)

式(8)中，Smk為m、k在河網密度序列上的統計量，Smk的最小值所的m和k所對應的am與ak為曲線由急劇驟減到逐漸減少的“閾值點”和曲線由逐漸減少到逐漸平穩的“閾值點”。

3.2 研究流程

對研究區域內DEM處理后得到的匯流累積柵格數據進行集水面積閾值采樣獲得閾值序列，劃定各采樣集水面積閾值下的河網柵格數據，進而求得對應的河網密度序列。將兩序列使用雙側均值變點法進行運算，獲得對應的兩個“變點”。得到兩個集水面積閾值參數，提取對應的數字河網數據在衛星圖上疊加對比，分析提取準確度及特點。

4 研究結果

4.1 雙側均值變點結果

本研究中，南苕溪流域DEM經上述處理后取樣序列長度為100，得到閾值密度曲線成果詳見圖3a，雙側均值變點檢測成果詳見圖3b。

由圖3a可知，隨著集水面積閾值的增加，河網密度的下降有顯著的3個階段：①首先是河網密度隨集水面積閾值增大的急速下降階段，表征各微小坡面匯流的驟縮，由于流域中大量的坡面網格匯流累積數較低，因此該階段，集水面積閾值的變化將導致河網長度的急劇減小，從而引起河網密度的急劇陡降，直至達到表征坡面上河網鏈消失的閾值分界點，曲線才開始由急劇驟減到逐漸減少；②緊接著是以一個較慢的速度下降階段，表征小型溝道部分的逐步排除過程，由于流域中存在不少的坡面網格匯流累積數為中低等級，因此該階段，集水面積閾值的變化導致河網長度的逐漸減小，從而引起河網密度的逐漸下降，直至達到表征坡面上小型溝道上河網鏈消失的閾值分界點，曲線才開始由逐漸減少到平穩下降階段；③最后是一個平穩下降階段，流域中微小坡面匯流與小型溝道已完全排除，留下顯著的中大型溝道，由于流域中坡面網格匯流累積數較大的網格存在數量較少，因此該階段，集水面積閾值的變化導致河網長度的減小程度微弱，從而引起河網密度的下降程度也顯著平穩。

由圖3b可知，整體統計量Smk變化較為連續，整體起伏變化不復雜，呈現中間低，四周高的形態。其中，m=1，k=2；m=99，k=100；m=1，k=100這三類情況下統計量Smk最高，分割的閾值密度曲線3個區域內部變化趨勢一致性最低。而當m=21，k=45時統計量最低，表明此狀態下分割的閾值密度曲3個區域的內部變化趨勢一致性最高。對應的河網密度隨集水面積閾值增大而呈現的3個階段：急速下降階段、較慢的速度下降階段、平穩下降階段對應統計尺度上最合理。因此可得到該狀態下的對應集水面積閾值分別為：a21=2517和a45=8936，分別提取上述集水面積閾值下的流域數字河網疊加至衛星影像進行對比驗證。

圖3 (a)河網密度-閾值曲線，(b)雙側均值變點檢測

4.2 河網劃定結果

集水面積閾值為a21=2517和a45=8936所生成數字河網詳見圖4，由圖4可見，數字河網與地形吻合程度均較高，河網在上游山丘區山谷發育較為集中，下游平緩城市地帶未發現平行河網問題。集水面積閾值為2517時相較于集水面積閾值為8936時，河網往上游延伸的程度顯著增加，河網長度與河網密度也更高。從河網分布的狀況來看，南苕溪流域內小型溝道部分集中在中游及上游，下游多為坡面的微小坡面匯流。

局部衛星影像對比上更明顯看出數字河網在山谷地形間所體現的較高準確度，圖4(a)、(b)均為西北上游部分，該部分局部對比可看出a45=8936所生成數字河網在a21=2517所生成數字河網基礎上往干流部分收縮程度顯著提高，對于山谷上游部分匯流面積較小的部分均排除，僅保留匯流規模有所保證的主干河道與溝槽，因此對于長時間序列徑流、生態環境模擬，可選擇a45=8936所生成數字河網。關注程度集中在顯著的河道溝槽中，排除上游匯流面積較小干擾。

局部衛星影像對比圖4(c)、(d)位于西側和南側的河流末端分叉口，由對比圖可見，a45=8936所生成數字河網在a21=2517所生成數字河網基礎上河道分叉顯著降低，往上游的小型溝道未包含入數字河網范圍，而后者對于小型溝道的分叉、匯流、形態描述更為細致準確，因此對于短時間序列洪水等災害模擬，可選擇a21=2517所生成數字河網，暴雨情況之下，河流整體水位抬升，流域整體土壤含水量較高，平時水流較少的小型溝道內也將形成顯著的徑流過程，因此對于該部分描述更為詳盡的數字河網更適合雨洪情景下的地理表達。

總的來看，雙側均值變點法所得集水面積閾值為a21=2517和a45=8936生成的南苕溪流域數字河網均能達到準確、合理的基礎要求。而a45=8936所生成數字河網收縮程度更高，河道分叉程度更低，接近非汛期及長晴期階段流域河網條件；a21=2517所生成數字河網對于小型溝道的分叉、匯流、形態描述更為細致準確，接近汛期及降雨條件下的流域河網狀態。后期研究中可針對研究問題選擇不同特征的數字河網，以達到更符合實際場景的表達。