基于盒維數(shù)的雅魯藏布江流域水系分形維數(shù)及影響因素研究

費(fèi)俊源，吳鵬飛，劉金濤，2

（1.河海大學(xué)水文水資源學(xué)院，南京210098；2.河海大學(xué)水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，南京210098）

0 引 言

地貌水文過程研究的一大難點(diǎn)是定量化地描述水系結(jié)構(gòu)特征［1］，而水系結(jié)構(gòu)特征會(huì)對(duì)水資源空間分布［2］、聚居點(diǎn)規(guī)劃［3］、地貌演化［4］產(chǎn)生影響，這又使得此特征受到了廣泛關(guān)注。由于水系的局部與整體間具有自相似性，目前主要使用針對(duì)自相似結(jié)構(gòu)的分形理論描述水系結(jié)構(gòu)特征［5，6］。

在分形理論中，分形維數(shù)被用于量化描述非規(guī)則客體［7］。目前的水系分形特征研究主要依托數(shù)字高程模型（Digital Elevation Model，DEM），使用Horton 定律或計(jì)盒法獲取水系分維數(shù)［7，8］。除了針對(duì)中小流域的研究［5，9］，以較大流域?yàn)閷?duì)象的水系分形研究大多將流域劃分為少量子區(qū)域或子流域，單獨(dú)分析各子區(qū)域的分形特征。例如，竇明等［2］依照水資源分區(qū)將淮河流域劃分為13個(gè)區(qū)域，發(fā)現(xiàn)各區(qū)域的水系盒維數(shù)能夠反映區(qū)域的人類活動(dòng)、土地利用等情況；鄭楠炯等［8］將韓江流域劃分為9個(gè)子流域，發(fā)現(xiàn)子流域與整個(gè)流域間的分維數(shù)接近。但是，上述研究使用的分區(qū)較大，且分維數(shù)樣本數(shù)量較少，由于分維數(shù)與氣象［7］、地形地質(zhì)［10］間聯(lián)系緊密，少量的分區(qū)不足以充分反映大型流域內(nèi)氣象和地形條件強(qiáng)異質(zhì)性帶來(lái)的水系分維數(shù)的差異。

本研究以雅魯藏布江（以下簡(jiǎn)稱雅江）流域?yàn)檠芯繀^(qū)，采用自動(dòng)化提取方法得到流域內(nèi)大量不同面積尺度的子流域并通過計(jì)盒法獲取子流域的盒維數(shù)。基于此盒維數(shù)數(shù)據(jù)集分析雅江流域內(nèi)不同尺度的子流域水系分維數(shù)與降水、坡度、起伏度等氣象、地形特征的關(guān)系。

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

雅江流域位于西藏自治區(qū)南部，涉及拉薩市、山南市、日喀則市、林芝市4 個(gè)地級(jí)市部分地區(qū)，流域面積24.8 萬(wàn)km2，絕大多數(shù)區(qū)域海拔＞3 000 m，是世界上平均海拔最高的流域［11］。流域隨雅江自西向東延伸，南北方向相對(duì)狹窄。流域內(nèi)降水的空間分布差異大，同時(shí)存在干旱、半干旱半濕潤(rùn)、濕潤(rùn)地區(qū)［11］。

圖1 雅魯藏布江流域概況Fig.1 Overview of Yarlung Zangbo river basin

本研究使用1″（約30 m）分辨率的AW3D30 DEM，數(shù)據(jù)下載自日本宇宙航空研究開發(fā)機(jī)構(gòu)網(wǎng)站（https：//www.eorc.jaxa.jp/ALOS/en/aw3d30）。雅江流域范圍由更大范圍的DEM，經(jīng)填洼和流向識(shí)別［12，13］，依照雅江出口點(diǎn)提取。研究使用1 km 分辨率青藏高原地區(qū)（2000-2015）降水?dāng)?shù)據(jù)，下載自國(guó)家青藏高原科學(xué)數(shù)據(jù)中心（http：//data.tpdc.ac.cn），并以此計(jì)算出流域各位置的平均年降水量數(shù)據(jù)。

2 研究方法

2.1 水系及子流域提取

本研究首先對(duì)DEM 進(jìn)行填洼預(yù)處理，移除DEM 內(nèi)因采樣誤差導(dǎo)致的洼地和平地，再計(jì)算每個(gè)柵格單元的流向。然后，依據(jù)流向累積得到每個(gè)單元的上游匯水面積。通過設(shè)置臨界源面積（Critical Source Area，CSA），提取匯水面積大于CSA的單元作為雅江的主要水系。最后，以水系的全部源點(diǎn)和部分主要干支流交匯節(jié)點(diǎn)作為子流域窗口，提取其完整的上游區(qū)域作為不同尺度的子流域以供分析。具體方案如下：

首先，自行編程實(shí)現(xiàn)由Garbrecht 和Martz［12］提出的算法進(jìn)行填洼，然后使用Wu 等［13］提出的iFAD8 單流向算法獲取單元流向并進(jìn)一步得到上游匯水面積。

為了綜合考慮不同尺度流域間的分形維數(shù)差異，以250 km2作為CSA 閾值，提取出雅江流域的主要水系。將所有水系源點(diǎn)對(duì)應(yīng)的上游區(qū)域作為最小子流域。再選取2個(gè)上游入流匯水面積差異較小（較小入流的匯水面積不得少于較大入流的1/3）的單元作為主要干支流交匯點(diǎn)，提取對(duì)應(yīng)的完整上游流域作為較大面積尺度子流域。值得注意的是，不同尺度子流域間存在嵌套現(xiàn)象，即較小子流域可能被某一較大子流域覆蓋。這一方面是為了保證用于分析的是完整流域，另一方面則是為了保證充足的數(shù)據(jù)量。對(duì)此，本研究將所有子流域劃分至不同面積區(qū)間，主要針對(duì)尺度近似的子流域進(jìn)行分析。盡管依然存在少量處于同面積區(qū)間但相互嵌套的子流域，但由于后續(xù)使用的都是面平均參數(shù)，可以認(rèn)為這些子流域依然具有代表性。

對(duì)于每個(gè)子流域，通過設(shè)置新的CSA 閾值得到其內(nèi)部的水系，而本研究使用拐點(diǎn)法確定該閾值［7，14］。拐點(diǎn)法通過檢查不同匯水面積閾值下盒維數(shù)變化趨勢(shì)，選擇盒維數(shù)變化趨勢(shì)的拐點(diǎn)作為最終閾值。由于本研究使用的子流域數(shù)目較多，難以逐個(gè)確定每個(gè)子流域的最佳閾值，因此使用不同面積閾值下所有子流域的平均盒維數(shù)的變化趨勢(shì)確定拐點(diǎn)，將確定的拐點(diǎn)閾值作為所有子流域的最佳匯水面積閾值。

2.2 盒維數(shù)計(jì)算方法

首先將柵格DEM 區(qū)域劃分成大量邊長(zhǎng)為r的正方形窗口，然后檢索出存在水系的窗口數(shù)量N（r）。r、N（r）與盒維數(shù)D間存在如下關(guān)系［15］：

根據(jù)上述關(guān)系，使用一組r和N（r）數(shù)據(jù)點(diǎn)繪ln（r）-ln［N（r）］關(guān)系圖，由最小二乘法擬合的直線斜率的絕對(duì)值即為所求盒維數(shù)值。由于當(dāng)r的擴(kuò)大倍數(shù)超過2 時(shí)分維數(shù)會(huì)出現(xiàn)波動(dòng)［16］，因此研究將r的擴(kuò)大倍數(shù)固定為2，以5 倍柵格單元邊長(zhǎng)（約150 m）作為初始的正方形窗口邊長(zhǎng)r。

3 結(jié)果與分析

3.1 子流域劃分及匯水面積閾值確定

本研究從雅江流域中提取出了356 個(gè)子流域，其中依據(jù)259 個(gè)水系源點(diǎn)提取得到259 個(gè)最小子流域，并依據(jù)3.1 節(jié)的方法選取97 個(gè)主要水系節(jié)點(diǎn)提取了97 個(gè)較大尺度的子流域。然后，以250 個(gè)柵格為間距，共計(jì)算了9 個(gè)CSA 閾值（500～2 500 之間）對(duì)應(yīng)的平均盒維數(shù)，結(jié)果如圖2（a）所示?？梢园l(fā)現(xiàn)，平均盒維數(shù)隨著CSA 的增大而減小。當(dāng)以1 500 個(gè)柵格（灰色點(diǎn)）為分界點(diǎn)時(shí)，匯流柵格數(shù)目＜1 500 個(gè)（黑色點(diǎn)）和＞1 500 個(gè)（空心點(diǎn)）對(duì)應(yīng)的平均盒維數(shù)呈現(xiàn)出2 種不同的線性減小趨勢(shì)，且二者的擬合效果極佳。此外，如圖2（b）所示，各CSA 對(duì)應(yīng)的所有子流域盒維數(shù)間平均絕對(duì)誤差和均方根誤差都較小，整體與平均值接近，可見采用平均盒維數(shù)確定閾值具備代表性。因此，將拐點(diǎn)對(duì)應(yīng)的1 500 個(gè)柵格面積（約1.35 km2）作為匯水面積閾值提取水系。

圖2 平均盒維數(shù)隨上游匯流柵格數(shù)目的變化關(guān)系和所有子流域盒維數(shù)的平均絕對(duì)誤差和均方根誤差隨匯流柵格數(shù)目的變化關(guān)系Fig.2 The relationship between the average box dimension and the number of upstream grids and The relationship between the mean absolute error and root mean square error of all sub-basins’box dimension and the number of upstream grids

3.2 分形維數(shù)提取結(jié)果統(tǒng)計(jì)

根據(jù)流域的面積大小，研究以500、1 000、2 000、4 000 km2為面積尺度分區(qū)節(jié)點(diǎn)，分5 個(gè)面積區(qū)間對(duì)356 個(gè)子流域的氣候、地形特征以及分形維數(shù)分別進(jìn)行了統(tǒng)計(jì)（表1）。其中起伏度的定義為流域最大高程和最小高程之差?？梢钥闯觯崛〕龅拿娣e250～500 km2的子流域數(shù)目最多，為259 個(gè)，也就是有259 個(gè)由水系源點(diǎn)提取的最小子流域。值得注意的是，由于DEM單元匯水面積增長(zhǎng)的躍進(jìn)特征，由CSA 得到的水系源點(diǎn)的匯水面積不一定等于CSA 對(duì)應(yīng)的250 km2。其他4 個(gè)面積區(qū)間內(nèi)的子流域較少，但也均超過10個(gè)。本研究得到的分區(qū)盒維數(shù)均值最小為1.04，最大為1.17。這一數(shù)值范圍與鄰近流域的研究成果接近［5，16］，屬于較小的盒維數(shù)范圍。水系發(fā)育程度越高的區(qū)域盒維數(shù)越大［16］，由此體現(xiàn)出雅江流域水系發(fā)育程度較低。雅江流域的盒維數(shù)D≤1.6，表明雅江流域的河流地貌處于侵蝕發(fā)育的幼年期［17］。根據(jù)表1 中的數(shù)據(jù)可以發(fā)現(xiàn)，盒維數(shù)均值存在隨面積區(qū)間增大而增大的情況，但增幅并不明顯。這是由于較大流域包含更多下游河谷區(qū)域，該區(qū)域水系發(fā)育程度優(yōu)于上游山丘區(qū)。此外，不同面積子流域的平均坡度、平均海拔相近，起伏度呈現(xiàn)出隨流域面積增大而增加的趨勢(shì)。

表1 不同面積范圍的流域降水、地形特征及分維數(shù)統(tǒng)計(jì)Tab.1 Precipitation，topographic features and fractal dimension statistics of catchments with different area scales

3.3 分形維數(shù)和地形及氣象特征的關(guān)系

本研究分別分析了5 個(gè)面積區(qū)間內(nèi)盒維數(shù)隨平均坡度、起伏度以及年平均降水的變化關(guān)系，并用線性關(guān)系、指數(shù)型關(guān)系、對(duì)數(shù)型關(guān)系擬合了關(guān)系函數(shù)。對(duì)于所有數(shù)據(jù)組，3 種擬合關(guān)系的確定系數(shù)差值均小于0.04，因此下文均以最簡(jiǎn)單的線性關(guān)系進(jìn)行描述分析。

圖3 展示了盒維數(shù)隨平均坡度的變化。在5 種面積區(qū)間內(nèi)，盒維數(shù)均呈現(xiàn)隨平均坡度增大而減小的趨勢(shì)。這體現(xiàn)出雅江流域中高地形起伏的上游區(qū)域河流處于發(fā)育更早期的階段，越平坦的下游區(qū)域河流發(fā)育情況越好。此外，隨著面積區(qū)間的增大，盒維數(shù)與平均坡度間擬合的確定系數(shù)同樣增加。這是由于小流域的水系發(fā)育更容易受各類局地因素的影響，大流域更能體現(xiàn)盒維數(shù)（或河網(wǎng)密度）與平均坡度間的相關(guān)性。

圖4展示了盒維數(shù)與起伏度之間的相關(guān)關(guān)系。大致可以判斷出雅江流域內(nèi)子流域存在盒維數(shù)隨著起伏度的增加而減小的趨勢(shì)。但整體看來(lái)，盒維數(shù)與起伏度間的相關(guān)性不佳。由于高起伏度地區(qū)的平均坡度傾向于更大，這里的變化趨勢(shì)可以視作圖3中盒維數(shù)與流域平均坡度相關(guān)性的體現(xiàn)。

圖3 不同面積（s）范圍下分形盒維數(shù)與流域平均坡度的關(guān)系Fig.3 Relationship between box dimension and average slope of catchments with different area scales（s）

圖4 不同面積范圍（s）下分形盒維數(shù)與起伏度的關(guān)系Fig.4 Relationship between box dimension and average topographic relief of catchments with different area scales（s）

圖5 展示了不同條件下水系盒維數(shù)與年平均降水量的關(guān)系。雖然不同面積區(qū)間內(nèi)盒維數(shù)與平均年降水量的確定性系數(shù)存在波動(dòng)，但可以發(fā)現(xiàn)整體上不同面積區(qū)間的流域水系盒維數(shù)均呈現(xiàn)出隨年平均降水量增大而減小的趨勢(shì)。研究顯示的雅江流域內(nèi)部的盒維數(shù)隨年平均降水量的變化趨勢(shì)與王博等［16］統(tǒng)計(jì)全國(guó)多個(gè)流域得到的濕潤(rùn)半濕潤(rùn)地區(qū)分形維數(shù)普遍高于干旱半干旱地區(qū)的結(jié)果不同。結(jié)合分形維數(shù)與河網(wǎng)密度的對(duì)應(yīng)關(guān)系，此差異也許可以類比河網(wǎng)密度與降水量的U 型關(guān)系。Abrahams 等［18］證明了河網(wǎng)密度隨干旱程度呈U 型關(guān)系，即以某一干旱情境為節(jié)點(diǎn)，在相對(duì)此情境越干旱或越濕潤(rùn)的情況下，河網(wǎng)密度都會(huì)增加。因此，可以推斷分形維數(shù)與降水量可能也呈類似U 型關(guān)系，本研究和已有研究發(fā)現(xiàn)的特征規(guī)律大致分屬U 型的兩端。在本研究使用的子流域中，當(dāng)年平均降水量較大（＞600 mm）時(shí)，部分子流域并未遵從隨年平均降水量增大而遞減的趨勢(shì)，因此可以觀察到大量盒維數(shù)增大的子流域［如圖5（a）、圖5（b）、圖5（c）］，正是這些子流域?qū)е铝讼嚓P(guān)區(qū)間內(nèi)線性擬合效果較差，也間接佐證盒維數(shù)與平均年降水可能存在U型關(guān)系。隨著面積尺度的增加，受限于流域數(shù)目及降水量分布，難以觀測(cè)到水系分形維數(shù)與降水的U型關(guān)系。

圖5 不同面積范圍（s）下盒維數(shù)與年平均降水的關(guān)系Fig.5 Relationship between box dimension and average precipitation of catchments with different area scales（s）

4 結(jié)果與討論

針對(duì)已有水系分維研究中樣本數(shù)量過少的問題，本研究使用自動(dòng)化處理方法從雅江流域選取了合適的節(jié)點(diǎn)，提取出了356個(gè)不同面積的子流域，經(jīng)拐點(diǎn)法確定1.35 km2為最適合雅江流域分形研究的匯流面積閾值。然后使用計(jì)盒法計(jì)算了子流域的水系分形維數(shù)，發(fā)現(xiàn)雅江流域的河流處于地貌發(fā)育的幼年期。研究還劃分了5 個(gè)流域面積區(qū)間，分別分析了不同面積尺度下盒維數(shù)與流域平均坡度、起伏度以及平均年降水量間的關(guān)系。結(jié)果顯示，使用盒維數(shù)與平均坡度和平均年降水量間相關(guān)性較好，大致呈線性下降的關(guān)系，不同面積區(qū)間下擬合函數(shù)的確定系數(shù)則呈現(xiàn)出隨流域面積增大而增大的趨勢(shì)。而流域的分形維數(shù)與起伏度間的關(guān)系則較差。此外，根據(jù)本研究發(fā)現(xiàn)流域的盒維數(shù)與平均年降水量間關(guān)系呈現(xiàn)出非單調(diào)趨勢(shì)，推斷盒維數(shù)應(yīng)與河網(wǎng)密度一致，隨流域的干旱程度呈U 型關(guān)系。其中，較小面積尺度流域的盒維數(shù)與平均降水一定程度上呈U 型關(guān)系，而此關(guān)系在較大面積尺度的流域則不明顯，未來(lái)仍需要增加不同降水量的大面積尺度流域樣本數(shù)量做進(jìn)一步分析驗(yàn)證。

本研究主要關(guān)注雅江流域水系分形結(jié)構(gòu)與地形、氣象要素的關(guān)系。但是部分已有研究表明，水系的分形結(jié)構(gòu)與植被、土地利用等多種其他要素同樣存在關(guān)聯(lián)［2，19］，且這些要素與社會(huì)、經(jīng)濟(jì)關(guān)系密切。因此，后續(xù)研究需要關(guān)注相關(guān)領(lǐng)域，將水系分形特征在生產(chǎn)生活中的推廣應(yīng)用。 □