基于DEM的格尔木河流域水系分维分析

原晓平，刘少峰，田贵中，陈 李，喻 静

(1.中国地质大学地球科学与资源学院，北京 100083;2.辽宁工程技术大学土木与交通学院，阜新 123000)

0 引言

分形理论既是非线性科学的前沿和重要分支，又是一门新兴的横断学科，在当今世界十分风靡和活跃。其概念由美籍数学家Mandelbrot首先提出并应用于水文学等相关领域[1]。作为一种新的世界观和方法论[2]，分形理论最显著的特点是:本来看起来十分复杂的事物，事实上大多数可用仅含有很少参数的简单公式来描述[3]。这样便可抛开复杂的表面现象，实现对其本质的定量刻画，揭示内在量化特性，因此分形理论为地貌形态描述与模拟、地貌发育过程等的非线性研究开辟了新思路[4]。

目前大多采用平均坡度、流域面积、河流长度、分支比、相对高差和切割深度等单一指标中的一个或若干个指标进行流域地貌形态描述[5－6]。这些指标是传统地貌形态量化方法的简单沿用，只能表达流域地貌某一方面(如长度、起伏度)的形态特征，不能准确地刻画三维立体复杂流域地貌形态的整体性和综合性[7－8]。分形理论为流域的复杂地貌形态研究提供了定量依据，分形维数作为该流域地貌形态特征的综合指标，可揭示其复杂程度和整体性，克服了平均坡度等单一指标的缺点[9]。对于流域水系分形的研究，国内外一些学者已经进行了一定程度的研究[10－18]，但总体而言研究还比较初步。本研究在这方面进行了更深一步地探索，以期对深化相关研究有所裨益。

本文将GIS与分形理论相结合，通过提取不同汇流累积量阈值下的水系河网，并用盒维数法计算相应的水系分维数来描述流域的地貌形态特征，预测流域地貌的发育演化趋势。

1 研究区概况

选取格尔木流域为实验样区，其位于东昆仑山脉与柴达木盆地边界，是研究青藏高原环境演化的关键地区之一[19－20]。发源于昆仑山北麓的格尔木河是柴达木盆地的第二大河，分东、西2支:西支叫昆仑河，又叫奈齐格勒河，发源自东昆仑博卡雷克塔克山的冰川;东支叫舒尔干河，也叫雪水河，发源自东昆仑玉珠峰的冰川，是格尔木河的正源。2支流相向流经东昆仑山脉，出山口在纳赤台附近汇流后称格尔木河，向北流经格尔木市，最终注入柴达木盆地达布逊湖。格尔木河全长468 km，流域面积2.188 ×104km2，属于内陆水系，水利资源丰富(图1)。

图1 格尔木流域分布Fig.1 Distribution of Golmud River basin

2 实验原理与方法

DEM能反映一定分辨率的局部地形特征，配合一定算法就可以自动提取给定地理空间范围内的自然水系。虽然基于DEM提取流域水系特征的精度会受到DEM质量和分辨率的影响，但已有研究表明[15]，用30 m空间分辨率的DEM数据分析水系分维能得到合理结果，且满足计算的精度要求，与研究目的和研究尺度相一致。

本文采用2011年10月公布的ASTER－GDEM数据，根据NASA新一代对地观测卫星Terra的观测结果制作完成，空间分辨率约30 m，以GIS为技术平台，利用现有的嵌入式水文分析模块自动提取格尔木流域盆地内部的地表水系网络，并在Excel中进行统计、模拟及分维数计算。

2.1 提取水系网络

采用地表漫流模型算法进行水系网络的提取。该模型根据DEM栅格单元和8个相邻栅格单元之间的最大坡度来确定水流方向，并沿水流方向计算每个栅格单元的上游汇流能力，然后选定一个汇流阈值，将高于所选阈值的单元标记为河网。

1)生成无洼地DEM。DEM是对光滑地形表面的模拟。由于内插失真等原因，DEM表面存在一些畸变凹陷，使得计算水流方向时产生错误。因此，在计算水流方向之前，应对原始DEM进行洼地填平，得到无洼地的DEM。

2)计算水流方向。水流方向是指水流离开每一个栅格单元时的指向。在ArcGIS软件中水流方向采用D8算法，即在3像元×3像元的DEM栅格上，计算中心栅格与各相邻栅格间的距离权落差，取距离权落差最大的栅格为中心栅格的流出栅格。将中心栅格的8个邻域栅格编码，则中心栅格的水流方向可由其中的某一值来确定。

3)生成汇流栅格。在地表径流模拟过程中，汇流累积量是基于水流方向数据计算的。对每一个栅格而言，汇流累积量的大小代表上游有多少个栅格的水流方向最终汇流经过该栅格。汇流累积值较大者视为河谷，汇流累计值较小者视为较高的地方。通过计算汇流累积量生成汇流栅格。

4)提取河网。在汇流栅格的基础上，设定一个阈值，当某个栅格点上的累积流量超过了该阈值，则可以认定该栅格点属于某个水系范围。对于给定的流域，集水面积阈值越小，河网则越稠密，水系分维值越大。通过这种方法最终形成水系网络。

2.2 计算水系分维数

本文采用简单而实用的盒维数法计算水系的分维数，基本思路是用不同长度的正方形网格去覆盖被测水系，即假设粗视化网格边长为ε，包含河流片段的网格数目为N(ε)，当 ε不断变化时，N(ε)与ε－D成正比关系，即

两边分别取以a(a＞0)为底的对数，即

对应粗视化网格边长的一组系列(ε1，ε2，…，εk)，得到一组(N(ε1)，N(ε2)，…，N(εk))，以点(logaε，logaN(ε))为坐标作双对数图，用最小二乘法可拟合出一条直线，其中M为直线的截距，D为斜率，即

在无标度区范围内所拟合直线的斜率D为流域地貌形态分形信息维数。利用这种方法得到的分形维数称作盒维数。

由于分形维数只有存在于一定的尺度范围内才能满足自相似性，因此这种方法的关键在于无标度区间的确定[21]。经过多次尝试，本文采用实际距离为500～10 000 m不等的网格边长来进行网格分析。计算不同汇流累积量阈值下的河网分形维数如表1所示，其中相关系数均在0.995 3之上，最高可达0.998 9，说明在此区间内，格尔木流域具有很好的分形特征，且该区间必然处于无标度区间之内。

表1 河网分形维数Tab.1 Fractal dimensions of river network

续表

具体步骤:将创建出的不同边长的网格图转为面文件，用其覆盖不同汇流累积量阈值下的水系矢量图，并统计非空网格的数目，将得到的若干点对作散点图，斜率即为河网分维数。

2.3 计算高程面积积分值

一般设全流域的面积为S，流域内某条等高线以上的面积为s，流域最高点与流域最低点的高差为H，该等高线与流域最低点的高差为h，记面积比重 X=s/S，高程比重 Y=h/H，显然 X，Y 均在[0，1]内取值。根据一系列(X，Y)值，以X为横坐标，Y为纵坐标画出高程面积曲线，那么曲线左下方与坐标轴之间的相对面积即为流域的高程面积积分值[22]。本文借助ArcGIS Desktop软件平台中ArcToolbox中的 Surface Volume模块进行计算[23]。

3 结果分析

3.1 汇流累积量阈值与分维数的关系

根据分形定义可知，一条河流所对应的水系分维是唯一的，但在基于GIS技术自动提取河网水系的过程中，由于汇流累积量阈值的不同，提取出的河网密度也不同，从而导致计算出的水系分维存在差异，即最后得到一个分维数范围。格尔木河流域的汇流累积量阈值和分维数对数拟合结果见图2。

图2 汇流累积量阈值与分维数拟合图Fig.2 Fitting chart of confluence accumulation thresholds and fractal dimensions

由图2可以看出，汇流累积量阈值在3 000～8 000的范围内递增时，格尔木河流域的水系分维值在1.777～1.592之间递减。当汇流累积量阈值从3 000增加到6 300时，水系分维由1.777递减为1.645，降幅明显;而当汇流累积量阈值由6 300增加到8 000时，水系分维仅由1.645减少为1.592，递减趋势较为平缓:随着阈值由小到大，其对应的分维数开始以较快的速度递减，继而速度变缓，但总体上还是呈现递减趋势，说明同一流域汇流累积量阈值与其分维之间存在着密切关系，也说明水系分维有效地反映了流域的形态特征。二者拟合曲线为

式中:3 000≤x≤8 000;拟合曲线的相关系数为0.994，具有良好的对数拟合效果。因此可以通过该方程直接计算格尔木河流域在不同汇流累积量阈值的设定下，其对应的水系分维值。

3.2 水系的地貌发育阶段

地表水系的分形维数大小可反映流域地貌的侵蚀发育程度。河网密度与分维数大小呈正相关，即分维数越大，河网密度越大，河流发育越成熟，因此水系分维可作为该流域地貌形态分形特征的具体量化指标。据此，何隆华等提出以分维数Dn值来划分流域地貌侵蚀发育阶段[24]:当 Dn＜1.6时，流域地貌处于侵蚀发育阶段的幼年期，水系尚未发育成熟，河网密度小，地面比较完整，河流深切侵蚀剧烈，河谷呈V形;Dn越接近1.6，流域地貌越趋于幼年晚期，河流下蚀作用逐渐减弱，侧蚀作用加强，地面分割得越来越破碎，谷坡的分水岭变成了锋锐的岭脊。此时地势起伏最大，地面最为破碎、崎岖;当Dn=1.6时，地貌趋于幼年期末、壮年期初;当1.6＜Dn≤1.9时，流域地貌处于侵蚀发育阶段的壮年期，在河流的侧蚀、重力作用和坡面冲刷下，尖锐的分水岭山脊不断变低，谷坡变得平缓，山脊变得浑圆，地面由原来的峭峰深谷变成低丘宽谷;当1.9＜Dn≤2.0时，流域地貌处于侵蚀发育阶段的老年期，河流作用主要为旁蚀和堆积，下蚀作用已很微弱，地势起伏微缓，形成宽广的谷底平原。这种方法为划分流域地貌侵蚀发育阶段提供了定量化指标。依据上述划分方案，计算得到的格尔木河流域水系分维值大致处于1.6～1.8之间，因此可以判断出其地貌发育阶段处于壮年期，并向壮年晚期发展。流域地势起伏大，地面切割得支离破碎、崎岖不平，河流以侧蚀为主，下切作用相对较弱。

为了验证该结论的正确性，本文采用流域的高程面积积分值进行判断。高程面积积分的分析方法由美国地貌学家Strahler于1952年提出[25]，目前在侵蚀地貌发育阶段定量研究中被广泛应用。利用流域的高程面积积分值可以表示流域地貌面受侵蚀的程度，并以此可以判断出流域的发育阶段。按照戴维斯的地貌旋回理论[26]，当流域高程面积积分值＞0.6时，流域侵蚀剧烈，是地貌发育的不均衡阶段，称为幼年期;当高程面积积分值≤0.6时，侵蚀过程变得非常缓和，流域地貌形态基本上趋向于稳定状态，是地貌发育的均衡阶段。该阶段又可划分为2个时期:当0.35＜高程面积积分值≤0.6时，地貌发育阶段为壮年期;当高程面积积分值≤0.35时，地貌发育阶段为老年期[27]。本文以2 502 m为高程下限，6 212 m为高程上限，相对高差3 710 m，流域面积21 881.86 km2，画出河流的Strahler曲线图(图3)，计算得到格尔木河流域高程面积积分值为0.505，由此判断该流域属于地貌发育均衡阶段的壮年期。这与上述通过分维数计算所得到的结论相一致，进一步证明了格尔木河流域处于地貌发育阶段壮年期的结论是可靠的。

图3 格尔木河流域Strahler曲线图Fig.3 Strahler curve of Golmud River basin

4 结论与讨论

4.1 结论

本文运用ArcGIS的水文分析工具提取了格尔木河流域在不同汇流累积量阈值下的水系河网，并用盒维数法计算相应的水系分维数，得到了在一定无标度区间范围内汇流累积量阈值和分维数的拟合方程。最后根据分维数大小判断出格尔木河流域处于地貌侵蚀发育壮年期。主要得到以下结论:

1)将分形理论应用于地貌学是一种新的研究方法，基于GIS技术的分形维数计算方法既提高了工作效率，也改善了研究精度，再结合常用的盒维数法计算能够有效促进分形理论在水文学领域的快速发展，为地貌学者解释复杂地形提供了可能。

2)利用所得到汇流累积量阈值与水系分维数的拟合方程，可以快速确定流域河网阈值及相应水系分维数。

3)通过分维值的大小来判断流域侵蚀发育阶段，并用高程面积积分值来验证该结论，使得到的结果可信度更高。

4.2 讨论

虽然本文深化了GIS在分形研究中的应用，但仍然存在许多需要改进的地方:

1)用地表漫流模型从DEM中自动提取河网的算法存在误差，该算法只适用于受到人类活动影响较小的流域，是一种理想状态下的模拟河网，这是由于只考虑到地形因素，而实际水系的形成是受多方面因素综合影响造成的。因此，要得到符合实际的高精度模拟河网水系，还应考虑当地人类活动对自然水流流向的影响，从而对DEM数据进行修正以及对算法进行相关改进等。

2)本文得到的分形维数是表征流域整体地貌形态特征，而实际地表上的流域系统十分复杂，应该用多重分形来刻画。多重分形可以看作是大量具有不同无标度区分形的集合。由于各个子级流域所处的地理环境不同，要详细了解流域内局部沟谷地貌的侵蚀发育情况，需将流域划分成若干个小沟谷水系，对局部分形维数与整体分形维数间的关系进行深入研究。

3)决定分维数的因素有很多，仅从计算结果去分析流域的水系情况容易产生错误的结论，因此还应结合流域地貌类型、构造活动、植被类型等自然因素和土地利用、人类活动等人为因素来进行综合分析。本文得到的分维数的物理意义还需进一步验证。

[1]Mandelbrot B B.How long is the coast of Britain?Statistical selfsimilarity and fractional dimension[J].Science，1967，156(3775):636－638.

[2]王 桥，毋河海.地图信息的分形描述与自动综合研究[M].武汉:测绘科技大学出版社，1998:28－29.Wang Q，Wu H H.Description of the map information’s fractal dimension and automatic comprehensive study[M].Wuhan:Technical University of Surveying and Mapping Press，1998:28－29.

[3]陈 颙，陈 凌.分形几何学[M].北京:地震出版社，1998:12－14.Chen Y，Chen L.Fractal geometry[M].Beijing:Seismology Press，1998:12－14.

[4]崔灵周，李占斌，肖学年.岔巴沟流域地貌形态分形特征量化研究[J].水土保持学报，2004，18(2):41－44.Cui L Z，Li Z B，Xiao X N.Study on quantifying topographical fractal character of Chabagou watershed[J].Journal of Soil and Water Conservation，2004，18(2):41－44.

[5]张会平，杨 农，张岳桥，等.岷江水系流域地貌特征及其构造指示意义[J].第四纪研究，2006，26(1):126－135.Zhang H P，Yang N，Zhang Y Q，et al.Geomorphology of the Minjiang drainage system(Sichuan，China)and its structural implications[J].Quaternary Sciences，2006，26(1):126－135.

[6]张 廷.基于DEM的洮河流域构造地貌分析[D].北京:中国地质大学，2010.Zhang T.Tectonic geomorphology analysis based on DEM in Tao River basin[D].Beijing:China University of Geosciences，2010.

[7]王协康，方 绎.流域地貌系统定量研究的新指标[J].山地研究，1998，16(1):8－12.Wang X K，Fang Y.A new index of quantitative study on the drainage geomorphic system[J].Journal of Mountain Research，1998，16(1):8－12.

[8]朱永清，李占斌，崔灵周.流域地貌形态特征量化研究进展[J].西北农林科技大学学报:自然科学版，2005，33(9):149－154.Zhu Y Q，Li Z B，Cui L Z.The quantification study of watershed topography characteristics[J].Journal of Northwest Sci－ Tech University of Agriculture and Forestry:Natural Suence Edition，2005，33(9):149－154.

[9]崔灵周，李占斌，郭彦彪.基于分形信息维数的流域地貌形态与侵蚀产沙关系[J].土壤学报，2007，44(2):197－203.Cui L Z，Li Z B，Guo Y B.Fractal－ information－ dimension－based relationship between sediment yield and topographic feature of watershed[J].Acta Pedologica Sinica，2007，44(2):197－203.

[10]André R，André G R.On the fractal interpretation of the mainstream length－ drainage area relationship[J].Water Resources Research，1990，26(5):839－842.

[11]Rosso R，Bacchi B，La Barbera P.Fractal relation of mainstream length to catchment area in river networks[J].Water Resources Research，1991，27(3):381－387.

[12]Tarboton D G.Fractal river networks，Horton’s laws and Tokunaga Cyclicity[J].Hydrology，1996，187(1/2):105－117.

[13]Daya Sagar B S，Chockalingam L.Fractal dimension of non－ network space of a catchment basin[J].Geophysical Research Letters，2004，31(12):1－4.

[14]张宏才，汤国安.基于GIS的河网分形研究[J].西北大学学报:自然科学版，2006，36(4):659－660.Zhang H C，Tang G A.The research on drainage networks fractal by GIS[J].Journal of Northwest University:Natural Science Edition，2006，36(4):659－660.

[15]王 林，陈兴伟.基于DEM的流域水系分维计算与结果分析[J].地球信息科学学报，2007，9(4):133－137.Wang L，Chen X W.Study on relationship between extracted river network and fractal dimension based on DEM[J].Geo－ Information Science，2007，9(4):133－137.

[16]杨锦玲.基于分形的数字水系集水面积阈值确定研究[J].测绘科学，2011，36(4):33－34.Yang J L.Digital determination of catchment area threshold based on fractal dimension[J].Science of Surveying and Mapping，2011，36(4):33－34.

[17]梁春玲，梁海清，张祖陆.祖厉河流域水系分维与地貌发育阶段浅析[J].水土保持研究，2006，13(3):187－191.Liang C L，Liang H Q，Zhang Z L.Fractal dimension of the river networks and developing stage of drainage geomorphic of Zuli River basin[J].Research of Soil and Water Conservation，2006，13(3):187－191.

[18]丰 满，张 征，朱 凌，等.基于DEM的滇池流域水系提取及分维值探讨[J].环境科学与技术，2010，33(s2):11－14.Feng M，Zhang Z，Zhu L，et al.River network information extraction and fractal dimension values discussion of Dianchi Basin based on DEM[J].Environmental Science and Technology，2010，33(s2):11－14.

[19]Lewis A O，Robert C F，Ma H Z，et al.Late quaternary landscape evolution in the Kunlun Mountains and Qaidam Basin，northern Tibet:A framework for examining the links between glaciation，lake level changes and alluvial fan formation[J].Quaternary International，2006，154(155):73－86.

[20]陈艺鑫，李英奎，张 跃，等.末次冰期以来格尔木河填充—切割及驱动机制初探[J].第四纪研究，2011，31(2):347－359.Chen Y X，Li Y K，Zhang Y，et al.Late quaternary deposition and incision sequences of the Golmud River and their environmental implication[J].Quaternary Sciences，2011，31(2):347－359.

[21]林 峰，陈兴伟，王 林.基于DEM的九龙江流域水系分维估算[J].水资源与水工程学报，2009，20(1):29－32.Lin F，Chen X W，Wang L.Calculation of the fractal dimension in Jiulongjiang River basin based on DEM[J].Journal of Water Resources and Water Engineering，，2009，20(1):29－32.

[22]Harlin J M.Watershed morphometry and time to hydrograph peak[J].Journal of Hydrology，1984，67(4):141－154.

[23]党安荣，贾海峰.ArcGIS地理信息系统应用指南[M].北京:清华大学出版社，2003:325－328.Dang A R，Jia H F.ArcGIS geographical information system application guide[M].Beijing:Tsinghua University Press，2003:325－328.

[24]何隆华，赵 宏.水系的分形维数及其含义[J].地理科学，1996，16(2):124－128.He L H，Zhao H.The fractal dimension of river networks and its interpretation[J].Scientia Geographica Sinica，1996，16(2):124－128.

[25]Strahler A N.Hypsometric(area－altitude)analysis of erosional topography[J].Geological Society of America Bulletin，1952，63(11):1117－1142.

[26]曹伯勋.地貌学及第四纪地质学[M].武汉:中国地质大学出版社，1995:17－24.Cao B X.Geomorphology and quaternary geology[M].Wuhan:China University of Geosciences Press，1995:17－24.

[27]孙希华，姚孝友，周 虹，等.基于DEM的山东沂沭泗河流域地貌演化与水土流失研究[J].水土保持通报，2005，25(4):24－28.Sun X H，Yao X Y，Zhou H，et al.Research on erosion landforms evolution and water and soil loss in Yishusi Valley based on DEM[J].Bulletin of Soil and Water Conservation，2005，25(4):24－28.