基于地貌視角的喀斯特流域水系分維估算方法適應(yīng)性分析

安 全，賀中華，趙翠薇，梁 虹，焦樹林，楊朝暉

(1.貴州師范大學(xué)地理與環(huán)境科學(xué)學(xué)院，貴陽 550001；2.貴州師范大學(xué)國(guó)家喀斯特石漠化防治工程技術(shù)研究中心，貴陽 550001；3.貴州省山地資源與環(huán)境遙感應(yīng)用重點(diǎn)實(shí)驗(yàn)室，貴陽 550001；4.貴州省應(yīng)急管理廳，貴陽 550001)

0 引言

自然界許多系統(tǒng)、要素的空間組合結(jié)構(gòu)具有某種意義上的分形特點(diǎn)。諸多河流發(fā)育的支流與干流在空間結(jié)構(gòu)上呈現(xiàn)出一種典型的分形特征。通過傳統(tǒng)計(jì)算水系分維方法工作量大、費(fèi)時(shí)耗力，利用地理信息系統(tǒng)(geographic information system，GIS)技術(shù)并結(jié)合1∶5萬地形圖水系與30 m空間分辨率ASTER-GDEM數(shù)字高程數(shù)據(jù)估算流域水系分維，探討流域地貌發(fā)育與水文水系特征的關(guān)系是當(dāng)前研究的熱點(diǎn)問題。分形是自然界諸多客體表現(xiàn)局部特征與整體的相似性客觀幾何規(guī)律，在現(xiàn)代數(shù)學(xué)的發(fā)展中逐步演化成一支新的具有獨(dú)特理論體系與認(rèn)識(shí)的世界觀與方法論。在地學(xué)領(lǐng)域，地貌發(fā)育與河流水系的演化規(guī)律中本質(zhì)上都潛藏著分形的自相似性特點(diǎn)。分形理論在地學(xué)界的快速發(fā)展始于20世紀(jì)90年代，早在1977年Mandelbrot首先將分形理論引入水文學(xué)，成為水系分維理論的創(chuàng)始人。進(jìn)入20世紀(jì)90年代后，國(guó)內(nèi)外對(duì)水系分維的研究迅速展開。何隆華等[1]利用計(jì)盒法對(duì)全國(guó)14個(gè)大流域與67個(gè)小流域進(jìn)行了計(jì)算，并得出水系分維值1.6是流域地貌發(fā)育階段的臨界值。進(jìn)入21世紀(jì)后，基于數(shù)字高程模型(digital elevation model,DEM)提取水系估算水系分維的研究得到快速發(fā)展[2-5]。往后隨著科學(xué)技術(shù)的快速發(fā)展，特別是GIS與DEM的發(fā)展與建立使得基于ArcGIS系列技術(shù)與DEM模型的流域水系分維研究進(jìn)入大發(fā)展時(shí)期[6-15]。如王倩[16]等利用GIS技術(shù)，對(duì)秦淮河流域水系分維展開了研究，王林等[17]基于ArcGIS8.3和ArcVIEW3.2進(jìn)行了DEM模型的流域水系分維估算分析。

目前，水系分維的研究主要是區(qū)域流域水系分維與地貌發(fā)育關(guān)系，人工渠系與自然水系分維對(duì)比，水系分維與洪澇、徑流過程，水系分維與巖性識(shí)別，水系分維與居名點(diǎn)分布特征，水系分維的生態(tài)意義等方面的研究[18-20]。而針對(duì)中國(guó)西南喀斯特筑壩流域地區(qū)與地貌發(fā)育對(duì)水文特征影響方面的研究報(bào)道還較少。鑒于此，本文以黔中筑壩工程區(qū)龍場(chǎng)橋流域?yàn)檠芯繉?duì)象，以1∶5萬地形圖水系、30 m空間分辨率的ASTER-GDEM和1∶10萬貴州省綜合地貌圖為數(shù)據(jù)源，利用基于ArcGIS10.2的Horton-Strahler理論、水系柵格法、漁網(wǎng)法估算黔中筑壩工程區(qū)龍場(chǎng)橋流域水系分維，探討筑壩區(qū)地貌發(fā)育對(duì)流域水文特征的影響，以期為喀斯特地區(qū)的生態(tài)建設(shè)，經(jīng)濟(jì)的可持續(xù)發(fā)展，筑壩帶的洪澇風(fēng)險(xiǎn)評(píng)估以及流域綜合管理等方面提供一定參考。

1 研究區(qū)概況及數(shù)據(jù)源

1.1 研究區(qū)概況

研究區(qū)是以龍場(chǎng)橋水文控制斷面提取的流域，位于黔中水利工程區(qū)上游地區(qū)(圖1)。黔中水利樞紐工程是貴州省首個(gè)大型跨地區(qū)、跨流域長(zhǎng)距離水利調(diào)水工程，也是黔中地區(qū)生存和發(fā)展的生命線工程[23]。工程區(qū)內(nèi)是長(zhǎng)江流域和珠江流域兩大流域交錯(cuò)地帶，研究區(qū)流域面積4 074.15 km2，地形西高東低，海拔在1 138～2 800 m之間，流域的西部地貌類型為典型的喀斯特巖溶高原[27]，中東部峽谷、峰叢洼地相間分布，氣候?qū)儆趤啛釒Ъ撅L(fēng)氣候，夏季高溫多雨，冬季溫和濕潤(rùn)，四季分明，多年平均降雨量在1 100～1 400 mm之間，最熱月平均氣溫22 ℃以上，最冷月0～15 ℃左右。植被為亞熱帶常綠闊葉林、灌木林、草地和裸巖相間分布。土壤以黃壤為主，兼有非地帶性石灰土分布。流域內(nèi)發(fā)育的烏江水系南源三岔河干流全長(zhǎng)325.6 km，流域集水面積7 264 km2，是黔中地區(qū)大型水利工程集中地。

(a)龍場(chǎng)橋流域位置 (b)龍場(chǎng)橋流域DEM

圖1 研究區(qū)概況圖

Fig.1Researchareaoverview

1.2 數(shù)據(jù)源及其預(yù)處理

本文數(shù)據(jù)源包括1∶5萬地形圖水系(以下簡(jiǎn)稱地形圖水系)和30 m空間分辨率的ASTER-GDEM等。在對(duì)數(shù)據(jù)源的坐標(biāo)進(jìn)行一致性處理的基礎(chǔ)上，使用30 m空間分辨率 ASTER-GDEM數(shù)據(jù)進(jìn)行水系提取，利用ArcGIS10.2軟件的Hydrology水文分析工具集對(duì)從中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心的地理空間數(shù)據(jù)云(http://www.gscloud.cn)獲取的ASTER-GDEM數(shù)字高程數(shù)據(jù)進(jìn)行填洼、流向、匯流累積量計(jì)算處理，利用ArcGIS10.2軟件的Con命令設(shè)定匯流閾值1 500提取研究區(qū)河網(wǎng)柵格水系(以下簡(jiǎn)稱ASTER-GDEM水系)。

2 研究方法

利用ASTER-GDEM水系、結(jié)合0.6 m空間分辨率Google Earth影像、30 m空間分辨率Landsat 8 OLI遙感影像，對(duì)提取的筑壩水庫(kù)區(qū)平行狀水系依據(jù)Strahler理論進(jìn)行了分級(jí)。依次將提取水系柵格的邊長(zhǎng)設(shè)定為500～10 000 m，步長(zhǎng)為500 m；然后利用ArcGIS軟件的SQL查詢功能分別統(tǒng)計(jì)不同邊長(zhǎng)水系柵格所對(duì)應(yīng)的柵格數(shù)目，再對(duì)統(tǒng)計(jì)的柵格邊長(zhǎng)、數(shù)目與Strahler分級(jí)的矢量水系導(dǎo)入相關(guān)統(tǒng)計(jì)軟件，分別求柵格法與Horton定理的水系分維[28-29]。

2.1 Horton定理

Horton-Strahler計(jì)算原理為利用ArcGIS 10.2 提取1∶5萬地形圖水系，應(yīng)用Strahler分級(jí)原理對(duì)提取的水系進(jìn)行河流分級(jí)；同理，以ASTER-GDEM為數(shù)據(jù)源，利用ArcGIS 10.2的con函數(shù)對(duì)30 m空間分辨率DEM分別設(shè)定匯流閾值100～10 000自動(dòng)提取水系，并與1∶5萬地形圖尺度下提取的水系不斷做疊加試驗(yàn)，通過匯流閾值與河網(wǎng)密度的擬合曲線變點(diǎn)分析，確定用以提取研究區(qū)河網(wǎng)水系的最佳匯流閾值為1 500，再根據(jù)Strahler河網(wǎng)分級(jí)理論對(duì)矢量河網(wǎng)進(jìn)行分級(jí)；最后統(tǒng)計(jì)地形圖水系和ASTER-GDEM通過最佳閾值1 500提取水系的每一級(jí)河流分級(jí)的河道數(shù)及其河道長(zhǎng)度，計(jì)算2個(gè)數(shù)據(jù)源提取水系的的分叉比與河長(zhǎng)比，對(duì)分叉比與河長(zhǎng)比取對(duì)數(shù)，比值即為所求分維。計(jì)算公式如下：

Rb=Ni-1/Ni，i=2,3,…,n，

(1)

(2)

(3)

式中：Rb為分叉比；RL為河長(zhǎng)比；Ni為各級(jí)河道數(shù)；i為河道等級(jí)；DH為Horton-Strahler水系分維值。

2.2 柵格法

柵格法又稱覆蓋法、網(wǎng)格法、計(jì)盒維數(shù)法，其原理是：取邊長(zhǎng)為r的正方形網(wǎng)格與水系圖求取交集，得到水系覆蓋的網(wǎng)格數(shù)目為N(r)，當(dāng)r不斷變化時(shí)，與之相對(duì)應(yīng)會(huì)得到一系列的N(r)值，二者關(guān)系為

N(r)∝r-D。

(4)

左右兩邊求取對(duì)數(shù)，以點(diǎn)(lgr,lgN(r))為坐標(biāo)作雙對(duì)數(shù)圖，采用最小二乘法可擬合出一條直線，即

lgN(r)=-Dlgr+b，

(5)

式中：r為正方形網(wǎng)格的邊長(zhǎng)；N(r)為對(duì)應(yīng)邊長(zhǎng)網(wǎng)格與水系圖求取交集所覆蓋的網(wǎng)格數(shù)目；b為待定系數(shù)；D為雙對(duì)數(shù)曲線的斜率值，即所求的水系分維[25]。

3 結(jié)果與分析

3.1 黔中筑壩流域水系分維提取分析

利用ArcGIS10.2軟件對(duì)提取水系進(jìn)行Strahler二次分級(jí)，通過Horton-Strahler法求得地形圖水系與研究區(qū)ASTER-GDEM水系的水系分維，再通過水系柵格法、漁網(wǎng)法對(duì)提取地形圖水系與ASTER-GDEM水系，分別依次將柵格邊長(zhǎng)設(shè)定為500～10 000 m，步長(zhǎng)為500 m，分別統(tǒng)計(jì)不同柵格邊長(zhǎng)的柵格個(gè)數(shù)，并求取柵格邊長(zhǎng)與不同邊長(zhǎng)柵格數(shù)的對(duì)數(shù)，利用Origing9.1統(tǒng)計(jì)軟件作柵格邊長(zhǎng)與不同邊長(zhǎng)柵格數(shù)雙對(duì)數(shù)擬合線，求其斜率即為分維值，統(tǒng)計(jì)結(jié)果如表1—4所示。

表1 基于地形圖水系的Horton-Strahler法參數(shù)Tab.1 Horton-Strahler method parameters based on water system of topographic map

表2 基于ASTER-GDEM水系的Horton-Strahler法參數(shù)Tab.2 Horton-Strahler method parameters based on water system of ASTER-GDEM

表1和表2的統(tǒng)計(jì)結(jié)果表明，河流數(shù)與河流平均長(zhǎng)度隨著河流級(jí)數(shù)的增大而減少。河流的分叉比隨著河流級(jí)數(shù)的增大而增大，與河流數(shù)、河流平均長(zhǎng)度成反比；河長(zhǎng)比隨河流級(jí)數(shù)越大而減少?？傮w而言，河流分叉比、河長(zhǎng)比除了受到河流數(shù)及其河流長(zhǎng)度的影響，與河流的分級(jí)存在一定的相關(guān)性。

表3和表4是地形圖水系與ASTER-GDEM水系柵格參數(shù)和漁網(wǎng)法參數(shù)。統(tǒng)計(jì)結(jié)果表明，柵格水系法與漁網(wǎng)法統(tǒng)計(jì)的柵格數(shù)目隨著試驗(yàn)設(shè)定的柵格邊長(zhǎng)的增加而減少，柵格邊長(zhǎng)與柵格數(shù)目之間存在明顯的線性關(guān)系。

表3 地形圖水系與ASTER-GDEM水系柵格參數(shù)Tab.3 Water system of topographic map and ASTER-GDEM grid parameters

表4 地形圖水系與ASTER-GDEM漁網(wǎng)法參數(shù)Tab.4 Water system of topographic map and ASTER-GDEM fishing net method parameters

表5是結(jié)合表1至表4統(tǒng)計(jì)數(shù)據(jù)，利用分維計(jì)算公式根據(jù)表1和表2分別求出1∶5萬地形圖提取水系和ASTER-GDEM提取水系的Horton-Strahler分維；利用Origing9.1統(tǒng)計(jì)軟件根據(jù)表3和表4分別作1∶5萬水系和ASTER-GDEM的柵格水系法及漁網(wǎng)水系法試驗(yàn)得到柵格邊長(zhǎng)與柵格數(shù)的雙對(duì)數(shù)擬合曲線。結(jié)果見表5。

表5 不同方法提取地形圖水系與ASTER-GDEM分維值對(duì)比表Tab.5 Comparison of fractal values from water system of topographic map and ASTER- GDEM by different methods

從表5中可以看出，同一地區(qū)不同方法、不同數(shù)據(jù)源估算的水系分維相差較大。其中Horton-Strahler法估算地形圖水系差別顯著，根據(jù)Horton-Strahler法估算地形圖水系與ASTER-GDEM水系的分維值分別為1.69和0.66，二者相差1.03。根據(jù)河數(shù)定律與河長(zhǎng)定律的分叉比、河長(zhǎng)比理論取值范圍可知，針對(duì)二元三維地表地下結(jié)構(gòu)，溶洞、裂隙發(fā)育的典型喀斯特地貌而言，由于諸多干流沿途遇到溶洞，大型溶蝕管道、裂隙時(shí)會(huì)轉(zhuǎn)入地下流一段距離等特殊現(xiàn)象，使得在地形圖上提取水系時(shí)無法判斷地下的水流狀況導(dǎo)致提取水系誤差增大，根據(jù)地形圖水系計(jì)算的河流的分叉比、河長(zhǎng)比的比值范圍未在理論值范圍，以至于所計(jì)算出的水系分維值與水系柵格法、漁網(wǎng)法估算的分維值大相徑庭，不能穩(wěn)定的判斷研究區(qū)地貌發(fā)育狀況。與此同時(shí)，表1和表2根據(jù)Horton-Strahler法提取的地形圖水系與ASTER-GDEM水系的河段數(shù)差別尤為顯著，總體上相差1 174段，河流分級(jí)的河段數(shù)差別主要體現(xiàn)在1～3級(jí)。這種顯著差異可能是由于復(fù)雜的喀斯特地貌組合結(jié)構(gòu)所致。

然而，柵格水系法與漁網(wǎng)法估算的水系分維無論是地形圖水系還是ASTER-GDEM水系，估算的分維值均比較接近。就表5數(shù)據(jù)表明，通過ASTER-GDEM水系利用2種方法估算的分維略比地形圖水系大。另外，從水系提取的R2值來看，1∶5萬地形圖尺度下柵格水系法與漁網(wǎng)法提取水系分維的R2值分別為0.996 6和0.996 4，而基于ASTER-GDEM采用柵格水系法與漁網(wǎng)法提取水系分維的R2值分別為0.994 1和0.993 4，這表明前者提取水系的無標(biāo)度區(qū)間在0.996 4～0.996 6之間，而后者為0.993 4～0.994 1之間。

3.2 黔中筑壩流域水系分維與地貌發(fā)育關(guān)系分析

圖2為研究區(qū)根據(jù)1∶10萬貴州省綜合地貌圖矢量化的研究區(qū)地貌發(fā)育現(xiàn)狀圖和30 m空間分辨率ASTER-GDEM數(shù)字高程數(shù)據(jù)利用ArcGIS10.2制作的研究區(qū)山體陰影圖。表6為研究區(qū)各類型地貌類型及其面積比例。

(a)研究區(qū)地貌類型 (b)研究區(qū)山體陰影

圖2 研究區(qū)各流域地貌類型與山體陰影圖

Fig.2Landformtypesandmountainshadowmapsofvariouswatershedsinthestudyarea

表6 不同流域地貌發(fā)育參數(shù)及其面積比例Tab.6 Geomorphological development parameters and area ratio of different watersheds

①字母F，K1，K2分別表示地貌發(fā)育的侵蝕-剝蝕類型、溶蝕為主型、溶蝕-侵蝕類型。

從圖2(a)中可以看出，整個(gè)研究區(qū)地貌發(fā)育類型以中山谷地、峰叢谷地、峰叢洼地為主，其中龍場(chǎng)橋流域、陽長(zhǎng)流域少部分地區(qū)出現(xiàn)深切中山型地貌。從圖2(b)山體陰影圖可以看出，研究區(qū)地貌發(fā)育破碎程度大，河谷多呈“V”型，越往研究區(qū)的東部即陽長(zhǎng)、龍場(chǎng)橋一帶，地勢(shì)起伏越大。地貌發(fā)育組合以深切中山、中山谷地、峰叢洼地復(fù)型組合為主。從表6中可以看出，向陽流域地貌發(fā)育組合以峰叢谷地、峰林溶原(盆地)、淺切中山為主，所占流域面積比例分別為23.31 %，20.83%和38.04%，研究區(qū)比例分別為4.76%，4.25%和7.77%；陽長(zhǎng)流域地貌發(fā)育組合以K化中山谷地、峰叢谷地、峰叢洼地為主，所占流域面積比例分別為26.75%，47.43%和11.62%，研究區(qū)比例分別為10.48%，18.59%和4.56%；龍場(chǎng)橋流域地貌發(fā)育以K化中山谷地、峰叢谷地、峰叢洼地為主，所占流域面積比例分別為43.23%，11.95%和41.14%，研究區(qū)比例分別為17.46%，4.83%和16.61%。

3.3 不同流域地貌發(fā)育成因分析

從表6中可以看出，向陽流域、陽長(zhǎng)流域和龍場(chǎng)橋流域3個(gè)流域的地貌發(fā)育成因主要為侵蝕-剝蝕類型，溶蝕為主類型、溶蝕-侵蝕類型，向陽流域由于侵蝕-剝蝕成因下發(fā)育的地貌類型深切中山面積占整個(gè)流域的38.04%，這說明整個(gè)流域地貌的發(fā)育主要以侵蝕-剝蝕為主，其次為溶蝕、溶蝕-侵蝕型。陽長(zhǎng)流域由于溶蝕、溶蝕-侵蝕成因下發(fā)育的地貌類型峰叢谷地面積占整個(gè)流域的47.43%，其地貌發(fā)育主要以溶蝕、溶蝕-侵蝕為主。龍場(chǎng)橋流域由于溶蝕、溶蝕-侵蝕成因下發(fā)育的地貌類型K化中山谷地、峰叢洼地面積占整個(gè)流域的34.07%，故其地貌發(fā)育主要以溶蝕、溶蝕-侵蝕為主。從地貌發(fā)育的角度看皆屬于地貌發(fā)育幼年期的中晚期階段。

綜合水系估算分維值來看，根據(jù)根據(jù)何隆華的計(jì)算的水系分維與地貌發(fā)育階段的關(guān)系可知，當(dāng)D或DH≤1.6時(shí)，流域地貌發(fā)育屬于侵蝕發(fā)育階段的幼年期；當(dāng)D>1.6或DH≤1.89時(shí)，流域地貌發(fā)育屬于侵蝕發(fā)育階段的壯年期；當(dāng)D或DH>1.89時(shí)，流域地貌發(fā)育屬于侵蝕發(fā)育階段的老年期。Horton-Strahler法、水系柵格法、漁網(wǎng)法提取除了Horton-Strahler法估算水系分維與研究區(qū)實(shí)際出入較大外，水系柵格法、漁網(wǎng)法估算的水系分維與研究區(qū)地貌發(fā)育較為吻合。結(jié)合表5和表6數(shù)據(jù)來看，綜合分析典型喀斯特復(fù)雜的地貌結(jié)構(gòu)與水系分維關(guān)系，漁網(wǎng)法估算的水系分維與研究區(qū)實(shí)際地貌現(xiàn)狀最為吻合。根據(jù)漁網(wǎng)估算的水系分維可知，研究區(qū)利用漁網(wǎng)法估算的1∶5萬地形圖提取水系分維值為1.54，通過ASTER-GDEM提取水系估算的分維值約為1.60，這說明研究區(qū)正處于地貌發(fā)育階段的幼年晚期、壯年期早期。此結(jié)果與研究區(qū)實(shí)際發(fā)育地貌吻合。

4 結(jié)論與展望

本文以黔中筑壩工程區(qū)龍場(chǎng)橋流域?yàn)檠芯繉?duì)象，以1∶5萬地形圖水系、30 m空間分辨率的ASTER-GDEM和1∶10萬貴州省綜合地貌圖為數(shù)據(jù)源，利用基于ArcGIS10.2的Horton-Strahler理論、水系柵格法、漁網(wǎng)法估算流域水系分維，探討筑壩區(qū)地貌發(fā)育對(duì)流域水文特征的影響，結(jié)果表明：

1)喀斯特地區(qū)復(fù)雜地貌組合結(jié)構(gòu)下不同方法、不同數(shù)據(jù)源估算的水系分維相差較大。Horton-Strahler法、水系柵格法、漁網(wǎng)法估算地形圖水系分維值分別為1.69，1,53和1.54；估算ASTER-GDEM水系的分維值分別為0.66，1.59和1.60。其中Horton-Strahler法估算分維值差別顯著，差值達(dá)到1.03。

2)綜合分析Horton-Strahler理論、水系柵格法、漁網(wǎng)法估算喀斯特筑壩區(qū)不同數(shù)據(jù)源水系分維與實(shí)際地貌發(fā)育的關(guān)系可知，漁網(wǎng)法估算的水系分維與研究區(qū)實(shí)際地貌現(xiàn)狀最為吻合。根據(jù)漁網(wǎng)估算的水系分維可知，研究區(qū)利用漁網(wǎng)法估算的地形圖水系分維值為1.54，ASTER-GDEM水系估算的分維值約為1.60，這說明研究區(qū)正處于地貌發(fā)育階段的幼年晚期、壯年期早期。此結(jié)果與研究區(qū)實(shí)際發(fā)育地貌吻合。此外，3種方法估算喀斯特筑壩流域的水系分維精度排序?yàn)椋簼O網(wǎng)法>水系柵格法>Horton-Strahler法。