重慶市梁灘河匯流閾值分析與地貌發(fā)育研究

胡煜佳，鄭莉萍，張森林，邵景安，2

(1.重慶師范大學地理與旅游學院，重慶 401331；2.三峽庫區(qū)地表過程與環(huán)境遙感重慶市重點實驗室，重慶 401331)

0 引 言

隨著GIS技術(shù)的不斷發(fā)展以及在水文模擬中的廣泛應(yīng)用，對流域河網(wǎng)的提取大體分為兩種思路：一種是基于已有的遙感影像進行自動提?。涣硪环N則是利用DEM數(shù)據(jù)進行提取[1]。以DEM和ArcGIS為載體進行河網(wǎng)數(shù)字化、自動化提取是眾多研究者采取的方法之一；而河網(wǎng)提取以及水文特征信息獲取的關(guān)鍵在于每個格網(wǎng)單元的匯流累積量。匯流累積量取決于每個格網(wǎng)單元的流向以及流量分配[2]，采取不同計算流向的方法以及匯流閾值的擬定其自動提取的河網(wǎng)和流域水文特征也不一。單流向算法和多流向算法是研究中常見的兩種流向算法[3]，單流向算法如D8算法[4-5]、Lea算法[6]、Rho8算法[7]等；而多流向算法則認為水流具有分散性，分散流入其他下游格網(wǎng)單元[8]，例如FD8算法[9]、Dinf算法[5]等。針對匯流閾值的擬定，學者也采取了不同方法進行研究。左穎等[10]采取多閾值約束法、單一t值確定法、河網(wǎng)密度法綜合確定合理閾值；唐杰等[11]采用河網(wǎng)和河源密度法確定集水面積閾值，函數(shù)曲線的拐點則為最佳集水面積閾值；張曉嬌等[12]采用河網(wǎng)密度法和水系分形法綜合確定閾值；許斌等[13]結(jié)合分形理論分析了不同匯流累積閾值下河網(wǎng)分形維數(shù)和密度變化率特征。綜上可看出，大部分學者確定匯流閾值采用頻率較高的是河網(wǎng)密度法，并且都用多種方法綜合確定。由于基于DEM數(shù)字高程模型提取河網(wǎng)以及流域水文信息可行性高、易于操作。因此，本文選取了重慶市污染情況最復(fù)雜的次級河流之一——梁灘河流域，基于研究區(qū)12.5 m分辨率的DEM，以ArcGIS軟件為平臺提取流域水文信息，運用數(shù)理方法確定最佳閾值，最終提取的河網(wǎng)與TM影像的歸一化水體指數(shù)提取的河網(wǎng)與三調(diào)數(shù)據(jù)成果提取的河網(wǎng)進行對比分析，并運用Horton算法與起伏比法對流域地貌發(fā)育階段進行識別判定。

1 研究區(qū)域概況

梁灘河流域位于重慶主城區(qū)縉云山與中梁山之間的北碚丘陵谷地及中梁山的狹長槽谷地帶，地貌分區(qū)為川東平行嶺谷區(qū)，地勢北低南高，東西高中間凹。該河流發(fā)源于九龍坡區(qū)走馬鎮(zhèn)廖家溝水庫，流經(jīng)九龍坡、沙坪壩區(qū)、北碚3區(qū)，最后匯入嘉陵江，是嘉陵江下游右岸的一條主要支流、主城區(qū)流域面積第三大的次級支流，整條流域共有55條支流，干流長約88.7 km，流域面積516.18 km2。流域內(nèi)常住人口超過20萬，并且由于梁灘河流域內(nèi)曾經(jīng)發(fā)展的養(yǎng)殖業(yè)和小型加工業(yè)興盛，其水質(zhì)受農(nóng)業(yè)活動、工業(yè)、城鎮(zhèn)污水的影響較大，是全市污染情況最復(fù)雜、治理難度最大的次級河流之一。因此，選取該流域為研究對象，提取河網(wǎng)等河流水文信息以期為研究該流域范圍內(nèi)的農(nóng)業(yè)面源污染或水質(zhì)評價奠定基礎(chǔ)。其DEM影像見圖1。

2 數(shù)據(jù)來源與預(yù)處理

本文所用的DEM數(shù)據(jù)來源于重慶市生態(tài)環(huán)境大數(shù)據(jù)應(yīng)用中心，其分辨率為12.5 m，根據(jù)梁灘河流域邊界對其進行裁剪處理。由于原始DEM存在凹陷單元格以及在平坦地區(qū)的水流方向存在不確定性；因此需對原始DEM進行預(yù)處理。即填洼操作，削峰填低，創(chuàng)建無凹陷DEM。這是進行流域水文分析的基礎(chǔ)。在填洼過程中，限制值Z為指定凹陷點深度和傾瀉點間的最大允許差值并確定要填充的凹陷點和保持不變的凹陷點，由于Z值無法精確限定，在本次操作過程中未對Z限制作出任何指定，因此將移除所有峰值。本文所用第三次國土調(diào)查(以下簡稱“三調(diào)”)數(shù)據(jù)來源于重慶市國土局。TM影像來源于地理空間數(shù)據(jù)云，其分辨率為30 m，利用ENVI5.3對影像進行拼接、大氣校正、輻射定標、裁剪等預(yù)處理。

3 流域水文特征信息提取

3.1 流域水文特征信息提取

3.1.1 水流方向的確定

在對原始的DEM進行填洼操作后，需要確定每個柵格單元的水流方向，而流向判斷主要分為單流向算法和多流向算法[14]。單流向算法目前應(yīng)用最廣泛的則是D8算法。該算法將水流方向簡化為8個方向，即東北、東、東南、南、西南、西、西北和北被定義為有效的水流方向，分別用128、1、2、4、8、16、32、64這8個有效特征碼表示；再根據(jù)高程值從而決定其每個柵格的水流方向[4]。多流向算法如FD8算法[9]、Dinf算法[5]等。本文利用ArcGIS10.2中D8算法計算水流方向，但此時要檢驗填洼是否完全；若不完全，則需反復(fù)填洼。

3.1.2 匯流累積量的計算

每個柵格都會產(chǎn)生一定的水流方向，某一柵格的匯流累積量就意味著從上游流至該柵格的流量累計；因此，通過計算匯流累積量從而判斷每個柵格的匯流能力。匯流累積量為柵格單元數(shù)目與柵格單元面積乘積之和[1]。

3.1.3 河網(wǎng)提取與分級

河網(wǎng)提取是在計算匯流累積量的基礎(chǔ)上進行的。即，將匯流累積量大于某一設(shè)定閾值的柵格提取出來的便是河網(wǎng)；因此，閾值擬定對最終提取河網(wǎng)尤為關(guān)鍵。隨后再利用水文分析工具集下的Strahler河網(wǎng)分級法進行分級。

3.1.4 柵格河網(wǎng)矢量化與水文特征信息提取

利用水文分析工具集下的柵格河網(wǎng)矢量化工具，將表示線狀網(wǎng)絡(luò)的柵格轉(zhuǎn)換為表示線狀網(wǎng)絡(luò)的要素；再使用河流鏈接、分水嶺工具生成集水區(qū)；最后，利用轉(zhuǎn)換工具集下的柵格轉(zhuǎn)面工具將柵格數(shù)據(jù)集轉(zhuǎn)換為面要素，從而計算不同閾值下的流域面積、河流長度、河道個數(shù)等水文信息。

3.2 TM影像的流域水文信息提取

利用遙感影像提取流域水體信息的方法主要有目視解譯和計算機自動解譯2種[15]。本文利用常用的水體信息提取方法——歸一化水體指數(shù)(NDWI)進行自動提取。該方法由Mcfeeters[16]在NDWI的基礎(chǔ)上提出，利用不同波段對地物反射的特征——近紅外波段對水體的強吸收而植被強反射的特點，采用綠波段和近紅外波段的比值可以較大程度上抑制植被信息，突出水體，從而較好地提取流域水體(見圖2)。這里歸一化水體指數(shù)

NDWI=(Green-NIR)/(Green+NIR)

(1)

式中，Green為綠光波段反射率；NIR為近紅外波段反射率。

從圖2可看出，利用ENVI對TM影像的波段識別最終得出水體歸一化處理效果不是很理想，對河流以及小面積水域提取效果不佳或存在漏提或斷流現(xiàn)象，導(dǎo)致諸類現(xiàn)象的原因可能是該遙感影像分辨率不高，對于地物的識別較模糊，有待采用更高分辨率的遙感影像進行提取與研究。

4 結(jié)果與分析

4.1 匯流閾值對流域水文特征的影響

閾值直接影響河網(wǎng)的生成以及流域水文信息，不同匯流閾值條件下其提取出來的河流長度、流域面積、河道個數(shù)等有所差異。圖3顯示，隨著匯流閾值的增大河網(wǎng)越稀疏，支流數(shù)量越少。

通過屬性表計算出不同匯流閾值下的河流長度、流域面積、河網(wǎng)密度以及河道個數(shù)(見表1)。

由表1可知，當匯流閾值由10 000增加至30 000時，河道個數(shù)由200個減至60個，減少了70%，河流長度、流域面積、河網(wǎng)密度分別降低了46.72%、3.42%、44.83%；當匯流閾值由30 000增加至50 000時，河道個數(shù)、河流長度、流域面積、河網(wǎng)密度分別降低了48.33%、25.87%、1.53%、24.72%；當匯流閾值由50 000增加至80 000時，河道個數(shù)由31個減至20個，減少了35.48%，其余分別降低了18.29%、1.85%、16.75%；匯流閾值由80 000增加至100 000時，河道個數(shù)由20個減至13個，減少了35%，其余分別降低了9.76%、2.7%、7.25%。由以上分析可知，隨著匯流閾值的增加流域各水文特征值都在不斷降低，但是變幅最小的為流域面積，其他3類特征值波動明顯，并且可發(fā)現(xiàn)匯流閾值由10 000增加至30 000 時流域各水文特征值變化最為明顯，各特征值急劇減??；閾值大于30 000后，各水文特征值變化幅度明顯減小。

由于流域面積受閾值波動影響遠小于其他特征值，因此，本文分別用對數(shù)、冪函數(shù)、指數(shù)函數(shù)、多項式函數(shù)等函數(shù)對匯流累積量與河流長度、河網(wǎng)密度、河道個數(shù)進行趨勢分析，尋找最佳函數(shù)類型以及最佳閾值，經(jīng)過多次實驗發(fā)現(xiàn)冪函數(shù)擬合性最優(yōu)。匯流閾值與河流長度、河網(wǎng)密度、河道個數(shù)的冪函數(shù)關(guān)系如下：

Y1=55 755.948X-0.562(R2=0.997 4)

(2)

Y2=74.479X-0.527(R2=0.995 1)

(3)

Y3=6.909X-1.134(R2=0.998 8)

(4)

式中，Y1為河流長度，km；Y2為河網(wǎng)密度，km/km2；Y3為河道個數(shù)；X為匯流閾值；R2為擬合度。

對匯流閾值與河流長度、河網(wǎng)密度、河道個數(shù)的冪函數(shù)關(guān)系求二階導(dǎo)數(shù)，并做非線性擬合與趨勢分析(見圖4)。

由圖4可看出，河流長度、河網(wǎng)密度、河道個數(shù)隨閾值變化的二階導(dǎo)數(shù)趨勢一致，曲線的變化趨勢由急劇下降至緩慢減少再至趨于平緩，都是從10 000 增至20 000時，流域各水文特征值急劇下降，20 000增至30 000時，各特征值降幅減小但降幅也較其他值段明顯。因此，閾值30 000左右是該流域明顯的一個特殊點。匯流閾值與河流長度、河網(wǎng)密度、河道個數(shù)的二階導(dǎo)數(shù)關(guān)系曲線趨于平緩時所對應(yīng)的點為合理的匯流閾值。從圖4還可看出，當閾值在35 000時為曲線趨于平緩時的變點，其二階導(dǎo)數(shù)值趨向0，各特征值趨于穩(wěn)定，后與實際河網(wǎng)進行對比驗證。即

(5)

(6)

(7)

4.2 提取河網(wǎng)結(jié)果的分析

通過最佳匯流閾值提取的河網(wǎng)需要檢驗與實際河網(wǎng)的差距。本文以三調(diào)數(shù)據(jù)成果提取的河流作驗證分析(見圖5)，發(fā)現(xiàn)通過數(shù)理方法確定的最佳匯流閾值35 000提取的河網(wǎng)(見圖6)與實際河網(wǎng)較為接近，但仍有一定出入；從DEM中提取的河網(wǎng)很好地保留了河流主干道信息，但存在有一小部分支流不能很好地體現(xiàn)或存在“偽支流”現(xiàn)象，尤其是流域中西部存在部分偽河道。通過其他閾值提取的河網(wǎng)發(fā)現(xiàn)，當閾值低于35 000時，河網(wǎng)越密集，提取的河網(wǎng)發(fā)育較好，設(shè)定的閾值較小時生成的河網(wǎng)與實際河網(wǎng)對比發(fā)現(xiàn)很多都是“偽支流”，在實際中并非真正的河網(wǎng)。當閾值高于35 000時，河流長度變短，雖然都較好地保留了主干道信息，但絕大部分支流無法體現(xiàn)，支流發(fā)育程度較差，與實際情況也不符合。綜上所述，該研究區(qū)最佳匯流閾值為35 000，提取的河網(wǎng)在保留河流主干道信息之外其支流也能較好地與實際河網(wǎng)吻合。

4.3 地貌發(fā)育階段識別

4.3.1 Horton算法

Horton算法是依據(jù)地形圖或數(shù)字高程數(shù)據(jù)利用ArcGIS提取水系，再用Strahler河網(wǎng)分級法對水系進行分級[17]。即

RB=NW-1/NW

(8)

RL=LW/LW-1

(9)

(10)

式中，RB為水系分叉比；RL為河長比；w為水系級別號；NW為第w級河道數(shù)目；LW為第w級河道的平均長度；Db為Horton-Strahler水系分維值。

由上述表達式得出梁灘河水系分維值為1.072 4，根據(jù)王玉成等[18]和任娟等[19]的相關(guān)研究，提出當水系分維值Db≤1.6、1.6

4.3.2 起伏比法

判斷地貌發(fā)育階段的方法除了Horton法，還有經(jīng)典的起伏比法，該方法由Pike和Wilson通過數(shù)學方法推導(dǎo)出的估算面積高程積分值的典型方法，其計算過程方便快捷、精度也較高[17]。計算式為

HI=(Hmean-Hmin)/(Hmax-Hmin)

(11)

式中，HI為流域的起伏比；Hmean為流域的平均高程，m；Hmax為流域的最大高程，m；Hmin為流域的最小高程，m。

依據(jù)地貌發(fā)育理論與Strahler對地貌發(fā)育階段的定量劃分(見表2)。梁灘河流域HI值為0.514 9，處于地貌發(fā)育階段的壯年(偏幼)期。由以上2種方法可推斷出，梁灘河流域的地貌發(fā)育階段可能為幼齡期向壯年(偏幼)期發(fā)展。

表2 地貌發(fā)育階段劃分標準

5 結(jié)論與討論

目前，利用DEM和ArcGIS快速提取流域河網(wǎng)的方法廣泛適用于各個流域的水文分析，是眾多研究者采取的主要方法。本文以重慶市梁灘河為例，利用ArcGIS10.2水文分析模塊的相關(guān)工具提取河網(wǎng)等水文信息，并通過數(shù)理方法確定最佳匯流閾值與流域地貌發(fā)育階段判定，得出了以下結(jié)論。

(1)隨著匯流閾值的增大，河流長度、流域面積、河網(wǎng)密度、河道個數(shù)等流域水文特征值不斷減小，提取的河網(wǎng)越稀疏，支流數(shù)量越少。

(2)根據(jù)匯流閾值與河流長度、河網(wǎng)密度、河道個數(shù)的二階導(dǎo)數(shù)關(guān)系曲線確定其最佳匯流閾值為35 000。

(3)通過匯流閾值35 000提取的河網(wǎng)與TM影像的水體歸一化處理結(jié)果、三調(diào)數(shù)據(jù)成果提取的河網(wǎng)進行對比發(fā)現(xiàn)該結(jié)果具有可靠性，但仍有一定出入。筆者認為原因可能如下：①在使用軟件計算與提取的過程中，由于算法的局限性以及某些參數(shù)的設(shè)定不同，會對運算結(jié)果產(chǎn)生很大的影響，如在填洼中Z值的限定、計算流量中匯流閾值的擬定都會對最后提取的結(jié)果產(chǎn)生很大的影響，眾多學者尤其對于Z值的限定無法做出較好的回答，這便成為日后研究中所需關(guān)注的一大重點；②利用TM影像進行水體歸一化處理效果并不太理想，存在漏提或斷流現(xiàn)象。這可能是由于TM遙感影像分辨率不高，因此有待采用更高分辨率的遙感影像進行提取與研究。

(4)針對梁灘河流域地貌發(fā)育階段的識別，本文通過Horton算法與起伏比法，最終推斷出該流域地貌發(fā)育階段可能為幼齡期向壯年(偏幼)期發(fā)展。

就總體而言，在流域提取或水文分析中越來越廣泛運用GIS技術(shù)，相關(guān)理論研究和提取方法也日益成熟。但由于DEM數(shù)據(jù)精度的差異、處理方法不同和人為因素影響等，流域水文特征提取的精度會受到很大影響，如何不斷進行技術(shù)創(chuàng)新、提高精度是日后研究的重點。