水系提取的集水面積閾值確定方法比較研究

鄭佳瑋，陳興偉，2,3*，李孝成，阮偉芳

(1.福建師范大學(xué)地理科學(xué)學(xué)院，福建 福州 350007；2.福建省陸地災(zāi)害監(jiān)測(cè)評(píng)估工程技術(shù)研究中心，福建 福州 350007；3.濕潤(rùn)亞熱帶山地生態(tài)國(guó)家重點(diǎn)實(shí)驗(yàn)室培育基地，福建 福州 350007；4.福建省水利水電科學(xué)研究院，福建 福州 350001)

流域內(nèi)所有河流、湖泊等各種水體組成的水網(wǎng)系統(tǒng)，稱為自然界的水系，其中河流是水系最為重要的組成部分[1]。水系的提取通常是根據(jù)局部地形特征，借助一定算法，輔助獲取研究區(qū)域內(nèi)的水文要素信息[2]。數(shù)字高程模型(Digital Elevation Model，DEM)是一組有序數(shù)組，用以描述地面高程值空間分布的實(shí)體地面模型，能夠在一定空間分辨率下反映地區(qū)局部的地形特征。因此，成為目前用于水系提取分析的主要數(shù)據(jù)[3]?；贒EM數(shù)據(jù)提取的水系是對(duì)真實(shí)水系的概化表達(dá)，最大程度上反映了水系的真實(shí)特征。應(yīng)用DEM提取水系主要包括2個(gè)步驟，一是流向計(jì)算[4-6]，目前根據(jù)水流方向計(jì)算的模型可以分為單流向算法(SFD)和多流向算法[7](MFD)。單流向算法以D8算法最為典型[8]，在大尺度區(qū)域河網(wǎng)提取中，D8算法應(yīng)用最為廣泛[9-10]。

應(yīng)用DEM提取水系的第2個(gè)步驟是水系精細(xì)程度的確定[11]，除了王林等[12]采用最小河流長(zhǎng)度以確定提取水系精細(xì)程度的方法外，目前主要采用控制集水區(qū)域面積的大小來確定提取水系的精細(xì)程度。集水區(qū)域面積越小，提取的水系越密集；反之，提取的水系就越稀疏。因此，理論上存在一個(gè)集水面積閾值。集水面積閾值一般定義為支撐一條河道永久性存在所需的最小集水區(qū)域面積[13-14]。Giannoni等[15]提出根據(jù)坡地水流路徑的頻率分布法確定集水面積閾值，孔凡哲等[16]提出根據(jù)集水區(qū)域面積與河網(wǎng)密度關(guān)系曲線尋求最佳閾值；李麗[17]將實(shí)際水系的分形維數(shù)與不同集水區(qū)域面積下的水系分形維數(shù)進(jìn)行比較，得出閾值；吳泰兵等[18]提出基于改進(jìn)適度指數(shù)法確定集水面積閾值；Jefferson等[19]總結(jié)出形為AxSy的流域內(nèi)坡度、集水面積尺度化指數(shù)，即給定閾值，超過該閾值的柵格屬于河道；張鵬舉等[20]通過河道分支比確定集水面積閾值。上述方法普遍利用集水區(qū)域面積與水系特征量的變化曲線來確定集水面積閾值，通過目視判斷水系的集水面積閾值，存在一定的主觀性。楊邦等[21]采用1∶82萬紙質(zhì)地圖的水系；程中陽等[22]以1∶25萬的全國(guó)水系圖作為實(shí)際水系，選用河道平均坡降法、河網(wǎng)密度法、流域?qū)挾确植挤?、分形維數(shù)法等方法計(jì)算集水面積閾值，討論各方法的適用性。這種提取水系與實(shí)際水系對(duì)比方法的主要問題是何種比例尺水系適合作為判別標(biāo)準(zhǔn)及其與DEM數(shù)據(jù)分辨率的匹配問題，如何定量比較實(shí)際水系和提取水系之間差異大小的方法，也仍需進(jìn)一步研究。

陳冬平等[23]提出利用提取的水系河源與實(shí)際河源的距離誤差最小的方法確定集水面積閾值，以避免目視判斷集水面積閾值的主觀性。河源即河流發(fā)源地,指該河干流離河口處最遠(yuǎn)開始補(bǔ)給的地方[24]。理論上，提取水系的河源應(yīng)該與實(shí)測(cè)河源重合，而集水面積閾值的合理確定，可以準(zhǔn)確反映提取河源的位置，提取水系的河源能否合理反映實(shí)測(cè)河源的位置，可以作為一種判斷提取水系合理性的評(píng)價(jià)標(biāo)準(zhǔn)。

因此，本研究以河源的合理提取為評(píng)判標(biāo)準(zhǔn)，以東南沿海的晉江流域?yàn)檠芯繀^(qū)，選用基于DEM提取水系的河網(wǎng)密度法、分形維數(shù)法、河源最小誤差法和河道分支比法等4種主要方法，以比較分析各方法提取水系的差異，探討各方法提取水系的合理性。

1 數(shù)據(jù)來源與研究方法

1.1 研究區(qū)概況

位于福建省東南部的晉江，是該省的五大河流之一。晉江流域所處的地理范圍在117°44′～118°47′E，24°31′～25°32′N，位于亞熱帶海洋季風(fēng)氣候區(qū)，多年平均降水量為1 147 mm。流域氣候暖熱，冬短夏長(zhǎng)，熱量豐富，日照充足，熱量資源豐富。降水較多而時(shí)空分配不均，春夏多雨，秋冬少雨，降水量的年際變化也較大。臺(tái)風(fēng)、暴雨、低溫、干旱等災(zāi)害性天氣時(shí)有發(fā)生。晉江發(fā)源于福建省中部戴云山，河源位置位于泉州市安溪縣境內(nèi)桃舟鄉(xiāng)達(dá)新村梯仔嶺東南坡[25]，坐標(biāo)位于117°41′42″E，25°21′30″N，河流全長(zhǎng)302 km，總面積5 269 km2。晉江上游分為東溪和西溪兩大支流，沿途流經(jīng)永春、安溪、南安、晉江、鯉城、豐澤等縣市區(qū)(圖1)。

圖1 晉江流域位置

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

研究區(qū)數(shù)據(jù)是由美國(guó)航空航天局(NASA)和國(guó)防部國(guó)家測(cè)繪局(NIMA)聯(lián)合測(cè)量的STRM DEM數(shù)據(jù)(http://dwtkns.com/srtm30m/?tdsourcetag=s_pcqq_aiomsg)，分辨率為30 m，數(shù)據(jù)覆蓋全球陸地表面的80%以上，該測(cè)量數(shù)據(jù)覆蓋中國(guó)全境。

采用ArcGIS 10.2 Hydrology模塊的相關(guān)功能，通過以下步驟提取水系：①運(yùn)用Fill工具進(jìn)行DEM洼地填充；②采用基于D8算法的水流方向計(jì)算工具，利用Flow Direction工具完成水流方向的計(jì)算；③運(yùn)用Flow Accumulation工具計(jì)算集水區(qū)域面積；④通過設(shè)定不同的集水區(qū)域面積，得到不同密度的水系圖。

1.3 研究方法

1.3.1河網(wǎng)密度法

河網(wǎng)密度是指流域內(nèi)干支流的總河長(zhǎng)∑L(km)與流域面積A(km2)的比值，基本單位為(km/km2)。河網(wǎng)密度的計(jì)算公式為：

(1)

通過設(shè)定不同的集水區(qū)域面積，生成一系列模擬的河網(wǎng)，計(jì)算各模擬河網(wǎng)的河網(wǎng)密度，得到模擬河網(wǎng)密度隨集水區(qū)域面積的變化關(guān)系曲線。河網(wǎng)密度變化趨于平緩時(shí)對(duì)應(yīng)的集水區(qū)域面積即為集水面積閾值[16]。

1.3.2分形維數(shù)法

目前有2種方法廣泛用于水系分形維數(shù)的估算，一是基于分形定義的盒維數(shù)法，二是基于Horton定理。其中，基于盒維數(shù)法對(duì)水系的分形維數(shù)估算可以適應(yīng)不同類型的線性體，有助于分析計(jì)算復(fù)雜河網(wǎng)的形態(tài)[26-27]。

盒維數(shù)法以邊長(zhǎng)為r的柵格覆蓋分形體，將與分形體相交的格網(wǎng)計(jì)為非空，否則給出空值，并統(tǒng)計(jì)非空柵格數(shù)計(jì)為N(r)。當(dāng)r持續(xù)變化時(shí)，即對(duì)應(yīng)一組r1，r2，r3，……，得到一組N(r1)，N(r2)，N(r3)……。以點(diǎn)[lgr，lgN(r)]為坐標(biāo)作散點(diǎn)圖，若測(cè)度r與非空格網(wǎng)數(shù)N(r)在雙對(duì)數(shù)坐標(biāo)軸上呈現(xiàn)線性關(guān)系，那么用最小二乘法可以擬合一條直線:

lgN(r)=-DHlgr+C

(2)

式中DH——水系的分維值；r——柵格邊長(zhǎng)；N(r)——包含河流片段的非空柵格數(shù)；C——待定常數(shù)[28-30]。

繪制不同集水區(qū)域面積水系的分維值隨集水區(qū)域面積的變化關(guān)系圖，水系分維值變化趨于平緩的點(diǎn)所對(duì)應(yīng)的集水區(qū)域面積即為集水面積閾值[11]。

1.3.3河源最小誤差法

河源最小誤差法采用有關(guān)部門綜合考察確定的河源來評(píng)估提取水系的合理性。受地形、算法等其他因素的影響，設(shè)定不同集水區(qū)域面積提取的河源與實(shí)際河源會(huì)存在一定的距離誤差。運(yùn)用實(shí)際河源與提取水系河源兩點(diǎn)之間的距離誤差作為判斷水系提取合理性的標(biāo)準(zhǔn)，將河源距離誤差最小時(shí)的集水區(qū)域面積看作是集水面積閾值[23]。

1.3.4河道分支比法

理論上，真實(shí)水系河道分支比的均值趨于一個(gè)常數(shù)。雖然河道數(shù)目隨集水區(qū)域面積的變化而發(fā)生變化，但各級(jí)河道數(shù)量比幾乎保持不變，河道分支比的平均值在某一范圍內(nèi)趨于穩(wěn)定。據(jù)此，可以繪制平均河道分支比隨集水區(qū)域面積的變化關(guān)系圖，以河道分支比變化趨于穩(wěn)定時(shí)的集水區(qū)域面積，作為集水面積閾值[20]。

1.3.5集水面積閾值判定的一階導(dǎo)數(shù)法

河網(wǎng)密度法、分形維數(shù)法以及河道分支比法都是基于計(jì)算的水系特征量，繪制與集水區(qū)域面積的變化關(guān)系曲線，認(rèn)為關(guān)系曲線趨于平緩時(shí)所對(duì)應(yīng)的點(diǎn)即為集水面積閾值。由于難以用目視的方法判別關(guān)系曲線趨于平緩的點(diǎn)，可用關(guān)系曲線的一階導(dǎo)數(shù)判定集水面積閾值。分別對(duì)上述水系提取方法得到的河網(wǎng)密度、分形維數(shù)、河道分支比與集水區(qū)域面積的變化關(guān)系曲線進(jìn)行曲線擬合，進(jìn)而求得擬合曲線的一階導(dǎo)數(shù)，判定集水面積的閾值。

2 結(jié)果與分析

2.1 河網(wǎng)密度法

依次設(shè)定不同的集水區(qū)域面積，求得河網(wǎng)密度與集水區(qū)域面積的變化關(guān)系見圖2。當(dāng)集水區(qū)域面積(即柵格數(shù))增加，河網(wǎng)密度呈現(xiàn)下降的趨勢(shì)，下降的速率由急劇變?yōu)榫徛?，?dāng)柵格數(shù)從0增加到10 000時(shí)，河網(wǎng)密度下降較為劇烈，之后趨于平緩。對(duì)河網(wǎng)密度擬合函數(shù)求一階導(dǎo)數(shù)，在10 000時(shí)一階導(dǎo)數(shù)趨于平緩且接近0，此時(shí)的集水區(qū)域面積為9 km2，即為河網(wǎng)密度法確定的晉江流域集水面積閾值。

2.2 分形維數(shù)法

集水區(qū)域面積與水系分維的關(guān)系變化及其擬合函數(shù)曲線的一階導(dǎo)數(shù)見圖3。結(jié)果表明，當(dāng)柵格數(shù)達(dá)到18 000左右時(shí)，一階導(dǎo)數(shù)曲線趨于平緩且接近0，則分形維數(shù)法確定的晉江流域集水面積閾值為16.2 km2。

圖3 分形維數(shù)與集水區(qū)域面積變化關(guān)系

2.3 河源最小誤差法

在500～24 000柵格數(shù)范圍內(nèi)，選取系列柵格數(shù)生成不同密度的水系圖。在ArcGIS 10.2下加載晉江河源數(shù)據(jù)與提取的水系圖，使用measure工具量測(cè)提取水系河源與實(shí)際河源坐標(biāo)的距離誤差，從而得到距離誤差與柵格數(shù)的變化見圖4。從圖4看出，當(dāng)柵格數(shù)增加到6 000個(gè)左右時(shí)，距離誤差由大到??；而當(dāng)柵格數(shù)從6 000個(gè)左右繼續(xù)增加時(shí)，距離誤差又呈現(xiàn)由小到大的變化。為了進(jìn)一步確定最小誤差及其對(duì)應(yīng)的柵格數(shù)，通過函數(shù)擬合，得到曲線的方程為：

圖4 河源距離誤差與集水區(qū)域面積變化關(guān)系

y=1E-05x2-0.173x+1526.1

(3)

上式一階導(dǎo)數(shù)等于零時(shí)的柵格數(shù)為6 315，此時(shí)提取河源與實(shí)際的距離誤差最小。相應(yīng)地，河源最小誤差法確定的晉江流域集水面積閾值為5.68 km2。

2.4 河道分支比法

通過設(shè)定不同的集水區(qū)域面積得到不同密度下的晉江流域河道分支比，所得河道分支比與集水區(qū)域面積的變化關(guān)系及其擬合函數(shù)曲線的一階導(dǎo)數(shù)見圖5。圖5表明，當(dāng)柵格數(shù)達(dá)到12 000時(shí)，河道分支比擬合函數(shù)的一階導(dǎo)數(shù)趨于平緩且接近0，對(duì)應(yīng)的10.8 km2被認(rèn)為是提取水系的集水面積閾值。

圖5 河道分支比與集水區(qū)域面積變化關(guān)系

3 討論

3.1 集水面積閾值確定的不確定性

河網(wǎng)密度法、分形維數(shù)法以及河道分支比法皆是基于河流形態(tài)變化趨于穩(wěn)定來確定閾值，對(duì)閾值的判定主觀性較強(qiáng)，存在較大的不確定性。而河源最小誤差法通過實(shí)測(cè)河源與實(shí)際的距離誤差和集水區(qū)域面積的變化關(guān)系，擬合出的函數(shù)曲線是一元二次函數(shù)，具有先減小后增大的特點(diǎn)，存在最小值，即用距離誤差最小時(shí)所對(duì)應(yīng)的集水區(qū)域面積作為閾值，保證了閾值確定的合理性。

3.2 水系提取結(jié)果合理性分析

4種集水面積閾值確定方法得到的晉江流域水系分布見圖6。4種方法提取的水系支流的形態(tài)特征變化較為明顯，集水區(qū)域面積影響提取水系信息的詳細(xì)程度。基于分形維數(shù)法和河道分支比法確定的集水面積閾值較大，提取出的水系形態(tài)結(jié)構(gòu)較為簡(jiǎn)單，而基于河源最小誤差法以及河網(wǎng)密度法確定的集水面積閾值較小，提取出的水系結(jié)構(gòu)較為復(fù)雜，河道數(shù)目較多，結(jié)構(gòu)較為精細(xì)。

理論上，4種方法提取河源與真實(shí)河源應(yīng)該重合，實(shí)際上，提取河源與真實(shí)河源之間存在距離，稱為距離誤差。采用ArcGIS 10.2的Measure測(cè)量工具對(duì)提取河源與實(shí)測(cè)河源進(jìn)行距離測(cè)量，結(jié)果為距離誤差。4種方法提取的晉江流域河源的距離誤差見表1。由于算法、數(shù)據(jù)精度等因素的影響，4種方法提取的晉江流域河源與實(shí)際河源都有一定誤差，但河源最小誤差法的誤差最小，所提取的水系應(yīng)能更合理反映晉江流域的真實(shí)水系。

a)河網(wǎng)密度法

c)河源最小誤差法

表1 不同方法的集水面積閾值及其距離誤差

從原理上看，河網(wǎng)密度法、分形維數(shù)法和河道分支比法是從河流形態(tài)的角度探究水系的形狀結(jié)構(gòu)特征，以找到河流形態(tài)特征趨于穩(wěn)定的集水區(qū)域面積，但由于判斷穩(wěn)定與否的指標(biāo)不盡合理，閾值的確定方法有待進(jìn)一步改進(jìn)。而河源最小誤差法引入實(shí)際河源位置來控制水系的提取，一定程度上保證了提取水系的合理性。實(shí)際上，河源最小誤差法將提取水系與實(shí)際水系的對(duì)比，簡(jiǎn)化為提取河源與真實(shí)河源的對(duì)比。河源即是水系起點(diǎn)，由于水系自身形態(tài)受地形控制，合理確定了水系的起點(diǎn)，水系本身誤差應(yīng)該也比較小。此操作方法簡(jiǎn)便，在難以獲取實(shí)際水系數(shù)據(jù)的情況下具有較高的可行性。

4 結(jié)論

以晉江流域?yàn)槔?，選用河網(wǎng)密度法、分形維數(shù)法、河源最小誤差法和河道分支比法等4種方法，分別確定水系提取的集水面積閾值，提取晉江流域水系；以實(shí)際河源為基準(zhǔn)，比較不同方法提取河源與實(shí)際河源的距離誤差，分析4種水系提取方法的合理性，結(jié)果表明。

a)水系提取的河網(wǎng)密度法、分形維數(shù)法和河道分支比法，其集水面積閾值確定的目視判別，具有較大的不確定性；這些方法提取的河源，與實(shí)際河源的距離誤差都較大。

b)河源最小誤差法提取的河源與實(shí)際河源的距離誤差具有駐點(diǎn)，據(jù)此確定集水面積閾值所得到的距離誤差最小，以此確定的河流水系具有唯一性，從而保證了提取水系的合理性。因此，河源最小誤差法是晉江流域4種基于DEM真實(shí)水系提取方法中比較合理的方法。

本文工作僅以東南沿海的晉江流域?yàn)槔_展研究，所得結(jié)論的推廣還需進(jìn)一步考慮不同地貌特征流域的適用性、不同類型和分辨率DEM數(shù)據(jù)可能的影響等問題。