基于雙側(cè)均值變點法的數(shù)字河網(wǎng)閾值劃定
——以南苕溪流域為例

聶啟陽

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

1 引言

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

2 研究區(qū)域及數(shù)據(jù)

2.1 研究區(qū)域概況

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

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

圖1 研究區(qū)域地理地貌

2.2 數(shù)據(jù)來源及處理

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

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

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

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

3 研究方法

3.1 雙側(cè)均值變點法

雙側(cè)均值變點法為均值變點法的變式，均屬于數(shù)理統(tǒng)計方法。針對河網(wǎng)閾值問題其離散化模式為以下狀態(tài)：

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

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

河網(wǎng)密度計算使用：

(1)

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

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

計算各段內(nèi)的算術(shù)平均值：

(2)

(3)

(4)

進(jìn)而求得各段統(tǒng)計量：

(5)

(6)

(7)

求和各段統(tǒng)計量：

Smk=SⅠ+SⅡ+SⅢ

(8)

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

3.2 研究流程

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

4 研究結(jié)果

4.1 雙側(cè)均值變點結(jié)果

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

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

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

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

4.2 河網(wǎng)劃定結(jié)果

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

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

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

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

圖4 雙閾值劃定數(shù)字河網(wǎng)對比驗證

5 結(jié)論

(1)雙側(cè)均值變點法分割的閾值密度曲線3個區(qū)域內(nèi)部變化趨勢一致性最高，3個區(qū)域?qū)?yīng)河網(wǎng)密度隨集水面積閾值增大而呈現(xiàn)3個階段：急速下降階段、較慢的速度下降階段、平穩(wěn)下降階段。

(2)雙側(cè)均值變點法所得兩個集水面積閾值生成的數(shù)字河網(wǎng)均能達(dá)到準(zhǔn)確、合理的基礎(chǔ)要求。而高閾值所生成數(shù)字河網(wǎng)收縮程度更高，河道分叉程度更低，接近非汛期及長晴期階段流域河網(wǎng)條件；低閾值所生成數(shù)字河網(wǎng)對于小型溝道的分叉、匯流、形態(tài)描述更為細(xì)致準(zhǔn)確，接近汛期及降雨條件下的流域河網(wǎng)狀態(tài)。

(3)基于雙側(cè)均值變點法所劃定的數(shù)字河網(wǎng)可針對研究問題選擇高低閾值的數(shù)字河網(wǎng)，以達(dá)到更符合實際場景的表達(dá)。低閾值生成的河網(wǎng)適合做短時間序列洪水等災(zāi)害模擬，高閾值生成的河網(wǎng)適合做長時間序列徑流、生態(tài)環(huán)境模擬。