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

費俊源，吳鵬飛，劉金濤，2

（1.河海大學水文水資源學院，南京210098；2.河海大學水文水資源與水利工程科學國家重點實驗室，南京210098）

0 引 言

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

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

本研究以雅魯藏布江（以下簡稱雅江）流域為研究區，采用自動化提取方法得到流域內大量不同面積尺度的子流域并通過計盒法獲取子流域的盒維數。基于此盒維數數據集分析雅江流域內不同尺度的子流域水系分維數與降水、坡度、起伏度等氣象、地形特征的關系。

1 研究區及數據源

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

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

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

2 研究方法

2.1 水系及子流域提取

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

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

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

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

2.2 盒維數計算方法

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

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

3 結果與分析

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

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

圖2 平均盒維數隨上游匯流柵格數目的變化關系和所有子流域盒維數的平均絕對誤差和均方根誤差隨匯流柵格數目的變化關系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 分形維數提取結果統計

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

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

3.3 分形維數和地形及氣象特征的關系

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

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

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

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

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

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

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

4 結果與討論

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

本研究主要關注雅江流域水系分形結構與地形、氣象要素的關系。但是部分已有研究表明，水系的分形結構與植被、土地利用等多種其他要素同樣存在關聯［2，19］，且這些要素與社會、經濟關系密切。因此，后續研究需要關注相關領域，將水系分形特征在生產生活中的推廣應用。 □