基于SRTM的沅江流域水系提取与分形特征研究

Study on River Systems Extraction and Fractal Characteristics of Yuan River Basin Based on SRTM

摘要: 为完善沅江流域治理的基础数字地形信息，同时为沅江流域地理信息系统建设提供一定的信息支持，基于沅江流域的SRTM数据，利用ArcGIS软件的水文分析模块和提取分析模块等，提取了沅江流域的河网水系，同时采用网格法和分叉比河长比2种方法，计算了沅江流域的水系分维数，对比分析了网格法和分叉比河长比两种方法，并判断了此流域的发育情况、构造运动程度和侵蚀旋回阶段。研究表明，两种方法计算的沅江流域水系的分形维数分别为1.9327和1.9412，满足分支比河长比计算方式大于网格法计算方式的规律，前者计算工作量大于后者，但可靠度较后者高。从水系分维值可以判断出研究区的流域发育阶段属于老年期，该水系的发育情况较为完善，构造运动强烈，属于侵蚀旋回四级阶段。

Abstract: To perfect Yuan River basin based digital terrain information, at the same time for Yuan River basin geography information system construction to provide certain information support, based on Yuan River basin SRTM DEM data, using the hydrological analysis module and extracting analysis module of ArcGIS software, etc., to extract the Yuan River drainage basin, at the same time, using the grid method and bifurcation is longer than the river than the two methods. The fractal dimension of the river system in the Yuanjiang River Basin was calculated, the grid method and the bifurcation to river length ratio were compared and analyzed, and the development of the basin, the degree of tectonic movement and the cycle stage of erosion were judged. Research shows that: the fractal dimensions of the Yuanjiang River drainage system calculated by the two methods were 1.9327 and 1.9412, respectively, which satisfied the rule that the branch-to-river ratio calculated by the two methods was larger than that calculated by the grid method. The calculation workload of the former method was greater than that of the latter, but the reliability of the former method was higher than that of the latter. According to the fractal dimension of the river system, the development stage of the basin in the study area belongs to the old age, the development of the river system is relatively perfect, the tectonic movement is strong, and it belongs to the four stages of erosion cycle.

1. 引言

传统的水系描述上，主要采用水系形状等定性描述的方法来描述，而水流方向、汇流累积量等定量描述作为辅助。由于水系河网的外表的不可微分性，无法用欧氏几何进行描述，至今为止还未知如何准确地定量描述水系河网和流域地形地貌的结构特征 [1]。20世纪60年代，分形学的奠基人Mandelbrot [2] 在1957年提出了一种基于“粗糙化”海岸线的公式，第一次提出了分维的概念。近几十年来，采用分形理论对水系河网结构特征来定量描述也在学术界得到了充分的认可 [3] [4] [5]。Gupta等 [6] 运用统计学方法证明了分形理论中的部分经验公式。汪富泉 [7] 通过分析能量和熵的关系，阐述了冲击河流平面形态分形蜿蜒的原理。何自立 [8] 利用盒维数法计算了秦岭典型流域水系分维值，分析其与水文特征和流域形态的相关性及空间分异性。费俊源 [9] 提取了雅鲁藏布江流域356个不同面积尺度的子流域，并使用拐点法确定了最佳汇水面积阈值，计算了所有子流域的盒维数。张艳如 [10] 利用网格法计算燕沟流域水系分形维数，利用高程面积曲线法刻画了该流域的地貌侵蚀发育阶段。此后诸多学者对不同区域水系分形特征进行了大量研究。

为了完善沅江流域治理的基础数字地形信息，本研究使用ArcGIS软件提取沅江流域的SRTM数据和水流的流向、汇流累积量、河网水系等流域特征数据，同时采用网格法和分叉比河长比计算法2种方法，对比了这两种计算方式，同时对流域水系的分维含义作出了解释，借以分析流域的发育阶段、构造运动强度以及侵蚀旋回阶段等，并对分形方法对于沅江流域地貌研究的重要作用进行了探讨。本文的技术路线图如图1所示。

2. 研究区概况

沅江，别名为沅水，是湖南省洞庭湖流域四大水系之一，属于长江支流。源头为云贵高原苗岭东南部的斗篷山，流经贵州省、湖南省、湖北省、重庆市、广西省，在常德市德山汇入西洞庭湖。

Figure 1. Technology roadmap

图1. 技术路线图

沅江流域面积为89647平方公里，流域全长1033千米，是洞庭湖四大水系中最长的河流。沅江发育良好，支流众多，其中13条支流的流域面积超过1000平方千米，长度超过5千米以上的支流数目为1491条，属于羽状水系。沅江流域形状为平行四边形，由西南向东北斜，南北长东西窄。地势上跨越了我国的第二、第三级阶梯。流经的地貌主要为丘陵，植被覆盖较差，在坡面耕地上，有显著或严重的侵蚀现象发生 [11]。研究区如图2所示。

3. 研究方法与流程

传统分维计算几乎依赖手工操作，这种方法既耗时耗力，又精度低且容易出错。正是因为传统的分维计算的种种缺点，本研究选择基于SRTM的分析方法，从数据的下载到最后的结果分析，均是以计算机自动化的形式进行，精度好、效率高、省时省力。通过地图影像或遥感影像等来获取地表水系，只是停留在平均意义上的特征分析，其描述即各种参数的量化是从外部出发，不能对水系结构的本质特性进行剖析。数字高程模型(DEM)数据通过有序数组形式表示地表高程，包含着了丰富的水文资源信息与地形地貌信息，由DEM可以派生出等坡向、流量、坡度变化率等大量地貌特性信息 [12]。

Figure 2. Schematic diagram of the study area

图2. 研究区示意图

3.1. 数据来源

3.1.1. SRTM数据

本文所使用数据是SRTM数据，是DEM的一种。SRTM的中文全称为航天飞机雷达地形测绘使命(Shuttle Radar Topography Mission)，是美国国家航空航天局(NASA)和国防部国家测绘局(NIMA)联合测量。SRTM是基于WGS84和EGM96的DEM数据，在90%置信水平下的垂直误差均小于16 m。SRTM数据每经纬度方格提供一个文件，精度有1 arc-second和3 arc-seconds两种，称作SRTM1和SRTM3，或者称作30 M和90 M数据，SRTM1的文件里面包含3601 × 3601个采样点的高度数据，SRTM3的文件里面包含1201 × 1201个采样点的高度数据。本文使用的是90 M数据，即SRTM3数据。

3.1.2. SRTM数据预处理

一副生产出来的SRTM数据，由于其坐标系可能和研究区域数据有一定的出入，需要进行坐标校正。在生产过程中，也许会不可避免的产生很多误差，表现为SRTM数据中的错误点，需要进行消除操作。另外，下载的SRTM数据不是精确研究范围，所以需要镶嵌或者裁剪。

ArcGIS是一个可伸缩的、全面的GIS平台，可用其来收集、组织、管理、分析、交流和发布地理信息。ArcGIS包含了全球范围内的底图、地图数据、应用程序，可用于创建 Web 地图、发布GIS服务、共享地图、数据和应用程序。本文使用ArcGIS实现从DEM数据中提取沅江流域的河网，并完成水系分形维数的计算。使用ArcGIS的Create the Mosaic Data Set工具进行沅江流域SRTM的镶嵌，然后使用Hydrology工具生成沅江流域边界，并使用沅江流域边界和Crop Tool裁剪沅江流域SRTM，最后使用Mountain Shadow Tool进行山体阴影提取，并渲染，叠加相应色带的SRTM。

3.2. 沅江流域河流网络提取

3.2.1. 无洼地SRTM生成

无洼地SRTM生成包括计算初始流向确定Z限制值，并计算洼地、填充洼地 [13]，洼地填充后可以得到无洼地的SRTM，并根据使用fill工具设置Z限制值，确定要被填平洼地的最深值，小于该值的数据被保留，大于该值的洼地将被填充。

3.2.2. 水流方向计算

DEM中某像元的水流方向即为水流离开此栅格单元时的指向 [14]。在ArcGIS中，采用D8法计算水流方向。D8法将中心像元的8个相邻像元进行编码，每一个像元只有8种可能流向，中心像元的水流流向则由其中的某一值来决定 [15]。

D8算法计算水流方向，见公式(1)。

$S=\frac{H_z}{D_c}$ (1)

式中：H_z，D_c为2个像元中心之间的距离。

ArcGIS中水流方向计算：添加填洼后的SRTM数据后，使用Flow Direction计算流向。

3.2.3. 汇流累积量计算

在模拟地表径流时，通过计算水流方向数据来计算汇流累积量 [16]。汇流累积量的基本原理是：假设DEM每一个像元处有一个单位的水量，并根据自然水流由上往下的自然规律和研究区地形的水流方向数据，计算出各像元点所流过的水流量，最终计算出汇流累积量 [17]。ArcGIS中汇流累积量计算：使用Flow Accumulation进行汇流累积量计算。

3.2.4. 沅江流域河网提取

采用地表径流漫流模型，确定出河网的阈值，将栅格的汇流累积量与阈值比较，大于阈值的栅格被赋值为1，否则为NODATA。ArcGIS中河流网络提取：使用RasterCalculator工具提取网络，选择适宜的阈值，最终生成的沅江流域河网水系图如图3所示：

Figure 3. River Network and Water System in the Yuanjiang River Basin

图3. 沅江流域河网水系图

3.3. 沅江流域水系分形维数的计算

分形是非线性科学的前沿理论，主要用于描述整体与部分是相似的 [18]。根据水系的分维数可以判断一个水系的发育程度 [19]。本研究的水系分维数采用网格法与分叉比河长比2种方法进行计算。

3.3.1. 网格法计算计盒维数

网格法计算分形的计盒维数可以抽象的理解成需要多少网格才能覆盖线体，通过对网格边长的逐渐减少，并对比网格数随网格边长的变化，最后计算出分维数。其计算公式(2)如下。

$N\left(r\right)\propto r^{-D}$ (2)

等式两边取对数公式变换如下

$\lg N\left(r\right)=-D\lg r+A$ (3)

式中：N(r)是网格数，r是网格的边长，A为待定常数，D为分维数，式(3)斜率的绝对值为其取值。

在ArcGIS中，将矢量的水系河网转换成栅格，用像元的大小来控制正方形的网格大小，得到不同的栅格像元数。最后通过要素的属性表来获取每一次输入不同像元r大小时覆盖水系网络所用的网格数N(r)，并进行对数运算，结果见表1。将 $\lg r$ 和 $\lg N\left(r\right)$ 两列的值绘制成线性方程(见图4)，得二者关系式为 $\lg N\left(r\right)=-1.9327\lg r+11.308$， $R^2=0.9725$ (R为复相关系数)，即 $\lg N\left(r\right)$ 与 $\lg r$ 之间十分线性相关。最终算得沅江流域网格法计算的计盒维数即水系分维数为D = 1.9327。

Figure 4. Curves of lgr and lgN(r)

图4. lgr与lgN(r)曲线图

Table 1. Calculation table of dimension of grid method box

表1. 网格法计盒维数计算表

3.3.2. 分叉比与河长比计算分维数

河网水系结构的节点与分支构成了河网水系的基本结构。Horton [20] [21] [22] 在流域方面的研究提出了分形理论，在开展流域侵蚀发育的定量形态研究过程时，提出Horton定律——即水系组成定律。Horton-Strahler的水道数目定律公式表示为

$N_k=R_b\left(S-K\right)$ (4)

其中N_k是k级河流的数目，R_b是河流的分枝比

由Horton-Strahler定理可得计算K级河流的河流长度公式(5)。

$L\left(K\right)=N_K\times L_K=L_1{R}_L^{K-1}{R}_B^{\gamma -K}$ (5)

式(5)求和得河网总长计算公式如下(6)，(7)。

${\displaystyle \sum L}=L_1{R}_L^{\gamma -1}{\left|\frac{R_B}{R_L}\right|}^{\gamma -1}=L_1{R}_B^{\gamma -1}$ (6)

$L_1=L_{\gamma }{\left(\frac{1}{R_1}\right)}^{\gamma -1}$ (7)

取对数并将式(5)并代入式( 4)得式(8)。

${\displaystyle \sum L}=L_1^{1-\left(\frac{L_g{R}_B}{L_g{R}_L}\right)}$ (8)

当 $\gamma \to 0$，有 $L_1\to 0$，视L₁为量测长度r，得式(9)。

$L=N\times r~r^{1-d}$ (9)

由式(8)、式(9)得式(10)。

$D=\frac{L_g{R}_B}{L_g{R}_L}$ (10)

因为D取值范围[1, 2]，便得式(11)。

$D=\min \left[\max \left(1,\frac{L_g{R}_B}{L_g{R}_L}\right),2\right]$ (11)

由上述公式可知，河流分维数与河长比和分叉比有关。分支比越大，河流的分支数目越多，河道越蜿蜒，河流发育的程度越好，其水系的分维数也就越大。在计算分维数时，河网的平均分支比与平均河长比的比值即为分维数。

由水系河网分维值的计算公式，分别计算出沅江流域水系的平均分叉比(见表2)和平均河长比(见表3)，并将两者分别取对数，得到河网分维值D，见公式(12)

$D=\frac{L_g{R}_B}{L_g{R}_L}=1.9412$ (12)

Table 2. Calculation table of bifurcation ratio

表2. 分叉比计算表

Table 3. River length ratio calculation table

表3. 河长比计算表

有关资料表明，基于网格法的分形计算方式所得出的盒维数应小于分叉比河长比方法所得出的分形维数，本次计算符合该规律，由此可见是合理科学的。

4. 研究的结果与讨论

4.1. 分形计算方法对比分析

网格法与分叉比河长比计算结果分别为为1.9327和1.9412。分叉比河长比的工作量大、计算量大，需要对河网进行分级，并分别计算水系的平均河长比和平均分叉比，但是分叉比河长比的算法在水系分形研究中普遍使用，所以其科学性经过了大量研究成果检验，较网格法科学。网格法计算工作量小，效率高，无须河流分级操作，且只需用矢量转栅格控制像元大小。但是由于矢量河网无法描述河宽变化，一般视为同级的河流只能设置一致的宽度，所以在转换为栅格的过程中会有误差，无法真实反映河面面积，与真实网格数有差异，所以精度较分叉比河长比的算法低。

4.2. 检验分形维数的合理性

4.2.1. 水系合理性

为了考量SRTM提取河网的科学性，使用SRTM提取河流并设置不同阈值，最终与Google earth上的沅江流域河流地表形态进行对比，当阈值取1000时提取的较为符合实际情况。查询关于沅江流域的水文资料，沅江流域的确处于老年发育阶段，河网密度大，弯曲程度大。

4.2.2. 计算合理性

研究表明，同时基于网格法的分形计算方式所得出的计盒维数应小于分叉比河长比方法所得出的分形维数，本次计算符合该规律，由此可见是合理科学的，表明利用SRTM提取的水系并计算分维数具有可行性的。

4.3. 河网水系分形维数分析

4.3.1 . 流域发育阶段分析

根据李后强的流域地貌发育阶段划分，当水系分维值在1~1.6区间时，流域地貌发育为幼年期，在1.6~1.89区间时为壮年期，在1.89~2.0区间时为老年期。用网格法与分叉比河长比方法2种方法计算得到的水系分维值为1.9327、1.9412，按以上划分流域地貌发育阶段的原则进行判断，研究区域水系处于流域地貌发育的老年期，水系发育情况良好，河网密度大，地面侵蚀严重。并与沅江现有水文资料进行对比，该结果与区域实际情况相符合。

4.3.2 . 构造运动程度分析

根据孔凡臣提出的构造运动程度与水系分维值成正比，水系分维值越大，构造运动越强烈。水系分维值的取值范围是在1~2之间，基于SRTM提取的河网水系分别用网格法与分叉比河长比方法2种方法计算得到的水系分维值为1.9327、1.9412，本次计算的分维值趋于最大值2，说明沅江流域的构造运动强烈。

4.3.3 . 侵蚀旋回阶段分析

根据侵蚀旋回阶段划分原理，当水系分维值在1.0~1.4区间之内，属于侵蚀旋回一级阶段；当水系分维值在1.4~1.6属于侵蚀旋回二级阶段；当水系分维值在1.6~1.8属于侵蚀旋回三级阶段；当水系分维值在1.8~2.0属于侵蚀旋回四级阶段；基于SRTM提取的河网水系分别用网格法与分叉比河长比方法2种方法计算得到的水系分维值为1.9327、1.9412，水系分维值在1.8~2.0区间，本研究计算的沅江流域的分形值，表明沅江流域处于侵蚀旋回四级阶段。

5. 总结

研究结果表明：

1) 两种方法计算的沅江流域水系的分形维数分别为1.9327和1.9412，满足分叉比河长比方法计算分维数大于网格法计算方式的规律，前者计算工作大于后者，但是前者是大家公认的方法，已得到验证，前者可靠度较后者高。

2) 基于SRTM数据提取的河网水系，经过对比分析，选取不同阈值，最后发现阈值取1000时，用SRTM提取的沅江流域水系的流域特征与地貌分布与原始水系基本吻合，表明利用SRTM提取的水系并计算分维数具有可行性的。

3) 两种方式计算的分维值均大于1.89，其流域地貌发育位于老年期，水系河网密度大，发育情况良好，构造运动强烈，处于侵蚀旋回四级阶段。