基于数字地形信息的拉萨河流域河宽预测模型研究

倪飞宇，吴鹏飞，王 静，江玉吉，刘金涛,3,4，王万杰,4

(1. 河海大学 水文水资源学院，江苏 南京 210098； 2. 西藏自治区水文水资源勘测局，西藏 拉萨 850000； 3.河海大学 水文水资源与水利工程科学国家重点实验室，江苏 南京210098；4. 河海大学 水利学科专业实验教学中心，江苏 南京 210098)

1 研究背景

河道宽度作为河流的一项基本地形特征，被许多河流水动力学模型作为重要参数[1-2]。在这些水文模型中，河宽数据的精度对洪水流量的模拟结果具有显著影响[3-5]。受资金和人力资源限制，详细地手动测量河道各处的河宽分布在较大研究尺度下并不可行[6]。因此，当前的水文研究需要有能够自动获取流域内河宽分布的可靠手段[1]。

河宽的自动化获取手段主要分为两类，即图像识别方法[7-9]和统计拟合方法[4,6]。其中，图像识别方法可以从高精度卫星影像中直接判读出水体区域，然后测量出河流任意位置的水面宽度作为输出结果[8]。由于水面宽度存在季节性波动，使用旱季影像得到的水面宽度与真实的河道宽度间差异较大，因此测量河道宽度时还应包含以稳定河岸为边界的经常性过水的河床[6]。此外，地表覆盖、云雾等影响因素使得借助卫星影像判读出的河道往往并不连续，需要额外的手动修正[1]。统计拟合方法通过在河宽与上游集水面积、坡度等地形要素间构建经验性的统计关系，并在此基础上直接使用由数字高程模型(digital elevation model，DEM)计算得到的地形要素预测河宽[6,10-12]。由于图像识别方法主要应用于宽度达数千米的大型河流，在类似拉萨河流域的河道宽度较小的区域缺乏应用评价[12]，且高精度的青藏高原夏季卫星影像较难获取，因此只需要少量资料的统计拟合方法更适合于青藏高原的河宽获取研究。

现有的统计拟合方法普遍认可上游集水面积是河宽的重要影响因素[6,12]。此外，研究表明河宽同样受其他地形要素影响[11-13]。例如，石书缘等[11]利用Google Earth发现较大河道曲率下的河道宽度随点坝长度增大而增大；何蒙等[6]将坡度原点矩和集水面积结合成河宽因子，发现该因子与河宽间具有良好的相关性。蒋成伟等[12]使用Google Earth分别测量了拉萨河10个主要河段的平均河宽与5种地形要素平均值，以10组数据分析相关性。他们的研究发现河宽与集水面积、河道纵比降分别呈良好的线性、幂函数关系，与流域的平均坡度、河流凹度、河流弯曲度间相关性不佳。但是，蒋成伟等[12]使用的是自冬季枯水期影像中读取的水面宽度，其结果未必能准确反映真实河道宽度的分布规律。此外，该研究仅使用10组平均值构建的相关关系，存在一定偶然性，且并未对所构建的经验性模型进行精度评价。

在已有研究选取的地形要素中，河道纵比降与坡度原点矩相似，都被用于描述河道地形起伏变化[6]，而河流凹度的测量相对复杂，难以通过编程实现[12]。因此本文选取了上游集水面积、局部坡度、坡度原点矩、河道曲率4种地形要素，以拉萨河流域主要河段为研究对象，分析了河宽与各地形要素间的相关关系。然后选取其中相关性较好的地形要素，探讨建立河宽预测模型的可行性，并进行了模型验证。

2 数据来源与研究方法

2.1 研究区概况

拉萨河发源于念青唐古拉山，流经拉萨市后于曲水县汇入雅鲁藏布江，总面积32 588 km2，地理坐标为东经90°05′～93°21′，北纬29°15′～31°25′。当地属高原温带半干旱气候，年平均气温5.3℃，年平均降水量约500 mm[12]。流域内高山密布，河谷宽窄相间，河道上游为单一蜿蜒型河道，下游多呈分汊、辫状河道。拉萨河流域高程及主要水系分布如图1所示。

本研究以流域内集水面积大于700 km2的主要河段为研究对象，主要河段从1″(约30 m)分辨率的SRTM DEM使用D8方法并设置上游集水面积阈值提取得到。

图1 拉萨河流域高程及主要水系分布

2.2 数据获取

2.2.1 河宽数据读取 本文使用的河宽数据通过Google Earth提供的卫星影像测量得到。Google Earth作为一款应用广泛的地理信息软件，其提供的高精度卫星影像在目前的水文、地理研究中应用广泛[14-15]。本文依据拉萨河流域主要河段走向，大致每隔500 m量取一次，共量取河宽1 000个，其中800个用于寻找与地形要素之间的规律及建模，200个用于模型验证。所量取的河宽同时包含了水体部分和经常性的过水河床。为了保证有效河宽的准确性，量测过程中避开了部分大坝拦水区和辫状河道。由于部分河段尤其在下游辫状河道过多，量取的河宽点较少，因此读取得到的较大河宽数值较少。

2.2.2 地形要素计算 本文选取了上游集水面积、局部坡度、坡度原点矩、河道曲率4种地形要素进行分析，相关要素主要通过SRTM DEM提取得到。集水面积可以通过多种基于DEM的流向算法计算得到，但不同方法对山丘区河道栅格的集水面积影响不大[16]。例如使用ArcGIS内置的D8方法[17]和SAGA内置的D∞方法[18]计算得到的拉萨河主要河段栅格的集水面积间的平均绝对偏差约为0.2%，最大绝对偏差仅为1.4%。因此本文只使用了目前应用更为广泛的D8方法获取集水面积。

由于河道的DEM高程只是水面高程而非河底高程[19]，传统的DEM计算局部坡度的方法只能反映水面起伏程度而非河道断面所处位置的地形特征。针对这一问题，何蒙等[6]提出了坡度原点矩的定义，通过计算一定范围内所有DEM栅格局部坡度的平均值度量河道的地形起伏特征。本文同时计算了所有河道栅格的局部坡度和坡度原点矩用于分析，其中局部坡度由ArcGIS计算得到，坡度原点矩则根据何蒙等[6]的方法获取。

河道曲率指的是一段河流的实际长度与两端点间直线距离的比值。已有研究主要通过卫星影像直接测量曲率[11-12]。在Google Earth中，为确定河段的实际长度，本研究于正视状态下以“添加路径”工具大致沿河道中泓绘制一段大致以选定断面为中点的长度约200～300 m的曲折路径。以该路径长度作为河流实际长度，除以路径首尾端点间的直线距离即可得到实测的断面处真实河道曲率。由于该方法不适用于提取全流域河道所有局部位置的曲率信息，本文仿照常见的基于DEM的河道形态提取方法[20]，尝试借助提取的栅格河道测量曲率，探讨使用DEM提取河道曲率的可行性。该方法以需要计算曲率的河道栅格为中心，向前向后各追踪读取5个栅格，计算共11个栅格间的流径长度及两端的2个栅格中心间的直线距离，求取两者比值，进而得到计算的曲率。此外，本文在读取河宽的同时实测了相应位置的真实河道曲率用于验证。

2.3 研究方法

2.3.1 河宽与各地形要素间相关性分析 为了更好地发现不同河宽对应地形要素的整体规律，仿照McNamara等[21]的策略对相关要素值进行了整体平均处理。其原理是：从0 开始对每2 m范围内的所有河宽对应的要素值进行平均，获得的平均值作为该范围内中值河宽对应的地形要素值，用以更加清晰地表现河宽所对应的各地形要素的整体特征。使用线性函数、指数函数、对数函数、幂函数4种方法拟合河宽与各地形要素值的关系曲线，得到相应的确定性系数和回归方程，并取其中确定性系数最高的一种方法。

2.3.2 河宽预测模型构建 河宽预测模型采用2种方法构建：一是采用常见的方法，以河宽与单一地形要素的回归方程构建模型[12]；二是通过对何蒙等[6]构建河宽因子的方法进行优化，选取与河宽相关性较好的多个地形要素结合成新型河宽因子，以河宽与河宽因子的回归方程构建模型。

河宽因子的确定主要通过将地形要素划分为趋势性因素和局地性因素。趋势性因素反映河宽在流域上下游的变化趋势，主要考虑集水面积。局地性因素反映局部地形和河道形态变化对河宽的影响，具体影响因素从局部坡度、坡度原点矩、河道曲率3种地形要素中根据相关性选取。DEM中任一河道栅格的趋势性因素可以利用公式(1)进行表达，即：

(1)

式中：fa为趋势性因素；T为集水面积阈值，km2；A为当前栅格的集水面积，km2； max(A)为流域最大集水面积，km2。

考虑到河宽随局部地形特征的变化趋势，局地性因素的计算分为2种情况。

当河宽随第i个局地地形要素数值增大而增大时，局地性因素fsi用公式(2)进行表达：

(2)

式中：Meri为第i个地形要素的数值; max(Meri)为第i个地形要素的最大值。

当河宽随第i个局地地形要素数值增大而减小时，局地性因素fsi用公式(3)进行表达：

(3)

为反映多个地形要素的综合影响，将趋势性因素和所有n个局地性因素乘积作为河宽因子，fte即：

(4)

使用线性函数、指数函数、对数函数、幂函数4种方法拟合河宽与河宽因子间的回归方程，选取确定性系数最大的回归方程构建河宽预测模型。

3 结果分析

3.1 河宽与各地形要素间相关性分析

图2为实测河宽与各地形要素(集水面积、局部坡度、坡度原点矩和河道曲率)间的相关关系，图2中的拟合曲线均为4种拟合方式中确定性系数最高者。

图2(a)显示河宽与集水面积呈良好的幂函数关系，确定性系数达到了0.92。这是因为冲积河流的河宽与流量有良好的相关性，同时流量与集水面积之间的关系也十分密切，因此河宽与集水面积之间的关系较好。河宽与局部坡度的相关性较弱(图2(b))，确定性系数仅为0.41，这是因为借助DEM测得的局部坡度是水面比降而非河底坡度，河底坡度反映了河道地形起伏的特点，而水面比降不能准确反映河宽与地形起伏的关系。与只考虑可见水面的已有研究相比[12]，本研究判读的包含了经常性过水河道的河宽与集水面积间确定性系数更高，相关性更好。图2(c)表明，随着河宽的增大，坡度原点矩存在先增大后减小的趋势。经过反复试验，发现以30 m河宽为界，对数据点分别进行拟合，可以各自拟合出良好的函数关系。其中窄于30 m的河道与坡度原点矩呈线性关系，确定性系数为0.87，宽于30 m的河道与坡度原点矩呈对数关系，确定性系数为0.77，二者均优于河宽与局部坡度的相关性。此处以坡度原点矩反映河道起伏变化，得到河宽随地形起伏呈两段式变化规律，这与蒋成伟等[12]以河道比降反映河道起伏变化得到的河道起伏越缓则河道越宽的规律不一致。这是由于蒋成伟等[12]考虑的河段主要属于本研究选取河段的偏下游部分，不能反映上游偏窄河段的变化趋势。图2(d)分析了河宽与实测河道曲率间的幂函数关系，发现二者间存在较好的相关性，确定性系数为0.59。这是因为河流越曲折，河道所能容纳的水量就越大，而河宽就会相应减小。此处较好的相关性与已有研究任务的河宽与河道曲率无关的结论不一致，可能是由于已有研究只使用了经过平均化的10个点值确定规律，不能很好地体现规律。然而，河宽与使用DEM直接提取计算得到的河道曲率间相关性不佳(图2(e))，确定性系数仅为0.22，且计算的河道曲率与实测曲率间并不接近(图2(f))。这是因为基于DEM提取河道栅格的方法受限于栅格尺寸，栅格尺寸无法准确反映宽度不同于栅格尺寸的河道的形变趋势，进而无法准确描述河宽与地形因子之间的关系。

图2 实测河宽与各地形要素的相关关系

结果表明，河宽与集水面积的关系最为密切，且与坡度原点矩、实测河道曲率之间的相关性较好。但是，河道曲率难以准确地自动提取。因此，构建河宽模型主要考虑集水面积和坡度原点矩两个地形要素。鉴于实测河道曲率与河宽间较好的相关性，本研究同样采用了实测河道曲率参与建模计算，以探讨是否需要发展能够自动计算且准确反映真实河道曲率的方法用于河宽预测模型。此外，由于局部坡度在相关研究中应用广泛，本研究同样将其作为建模的一个可选参数。

3.2 河宽预测模型

3.2.1 模型构建 根据2.1节的相关性分析结果，分别选取与河宽相关性良好的集水面积和坡度原点矩构建河宽与单一地形要素的回归模型。

根据图2(a)的拟合结果，河宽W与集水面积A的回归方程为：

W=0.1319A0.6912

(5)

根据图2(c)的拟合结果，分两段确定河宽W与坡度原点矩s的回归方程。为了量化临界点特征，这里借助了30 m河宽在图2(a)中大致对应的2 100 km2集水面积。当集水面积A≥2 100 km2时，回归方程为：

W=-52.67 ln(s)+203.39

(6)

当集水面积A大于阈值700 km2，且A<2 100 km2时，回归方程为：

W=2.3926s0.8586

(7)

以河宽与河宽因子的回归方程构建模型时，除使用集水面积提供河宽因子的趋势性因素外，局地性因素由坡度原点矩、实测河道曲率、局部坡度中的部分项提供。根据局地性因素的不同，本研究共构建了4个河宽因子模型，其局地性因素分别包含坡度原点矩、坡度原点矩和局部坡度、坡度原点矩和实测河道曲率、所有3个地形指标(坡度原点矩、实测河道曲率、局部坡度)。

仅使用坡度原点矩提供局地性因素时，以2 100 km2为阈值，当集水面积A<2 100 km2时，局地性因素使用公式(2)计算，当集水面积A≥2 100 km2时，局地性因素使用公式(3)计算。再用公式(4)拟合趋势性因素和唯一的局地性因素得到河宽因子，对河宽W和河宽因子fte分段拟合的结果如图3(a)和3(b)所示。

当集水面积A<2 100 km2时的回归方程为：

W=28.78fte0.284

(8)

当集水面积A≥2 100 km2时的回归方程为：

W=153.66fte+42.758

(9)

使用坡度原点矩与局部坡度提供局地性因素时，分段拟合的结果如图3(c)和3(d)所示。

当集水面积A<2 100 km2时的回归方程为：

W=31.85fte0.2861

(10)

当集水面积A≥2 100 km2时的回归方程为：

W=147.08fte+47.078

(11)

使用坡度原点矩与实测曲率提供局地性因素时，分段拟合的结果如图3(e)和3(f)所示。

当集水面积A<2 100 km2时的回归方程为：

W=34.759fte0.2631

2.2 2012-2016年剖宫产指征及其构成比 2012-2016年本院剖宫产的主要指征为社会因素、瘢痕子宫、胎儿窘迫、胎儿异常、妊娠合并症及并发症，其中社会因素的构成比变化明显，由2012年的19.74%逐年下降至2016年的6.7%，差异有统计学意义(χ2=18.25，P<0.05)。见表2。

(12)

当集水面积A≥2 100 km2时的回归方程为：

W=311.24fte+43.489

(13)

使用坡度原点矩、局部坡度与实测曲率提供局地性因素，分段拟合的结果如图3(g)和3(h)所示。

当集水面积A<2 100 km2时的回归方程为：

W=38.302fte0.2655

(14)

当集水面积A≥2 100 km2时的回归方程为：

(15)

3.2.2 模型验证 使用200个实测点对集水面积模型、坡度原点矩模型及4种河宽因子模型效果进行验证，实测河宽和预测河宽的对比如图4所示。表1展示了对这6种模型预测的河宽与实测河宽进行误差统计的结果。由表1可知，仅以坡度原点矩为参数构建的坡度原点矩模型的预测结果最差，在集水面积A<2 100 km2和A≥2 100 km2的2个范围内都产生了最大的相对误差。以坡度原点矩和实测曲率为局地性因素构建的河宽因子回归模型的预测结果也较为不理想。相较坡度原点矩模型，使用集水面积为参数构建的回归模型具有更好的应用效果，这与图2展示的集水面积与河宽间较坡度原点矩与河宽间更好的相关性吻合，相关性越高，拟合的回归方程预测效果越好。以集水面积为趋势性因素，坡度原点矩与局部坡度或坡度原点矩、局部坡度与实测曲率为局地性因素构建的河宽因子模型效果较为理想，但预测效果均差于仅以坡度原点矩为局地性因素的河宽因子模型。以坡度原点矩为局地性因素的河宽因子预测模型在集水面积A<2 100 km2的上游区域拥有最好的预测效果，预测结果的相对误差小于集水面积模型。但在集水面积A≥2 100 km2的下游区域，该河宽因子预测模型的预测效果则差于集水面积模型，这同样与两个要素和河宽间的相关性强弱关系吻合。鉴于以坡度原点矩和实测曲率为参数构建的回归模型预测效果较差，可认为实测曲率不适合作为局地性因素，因此不进一步探讨自动获取准确河道曲率的方法。

根据预测模型的精度评价可以得知，集水面积是影响拉萨河流域河宽分布规律的主要因素，但在集水面积较小的上游河段，地形起伏同样具有相当明显的影响。最佳的拉萨河流域河宽预测模型可以以2 100 km2集水面积为阈值，当集水面积A<2 100 km2时使用基于坡度原点矩和集水面积的河宽因子模型，当集水面积A≥2 100 km2时使用集水面积模型。考虑到这种方法的处理过程较为繁杂，针对所有河段直接使用基于坡度原点矩和集水面积的河宽因子模型以避免复杂的情景划分步骤同样可取，其上下游相对误差均稳定且较小。

图3 以不同地形指标提供局地性因素时河宽因子与测得河宽相关关系

表1 6种模型平均预测误差统计结果

图4 模型预测河宽与实测河宽对比

4 结论与展望

本文依据Google Earth影像，对拉萨河流域主要河道的1 000个局部河段的河宽进行了测量。同时利用该流域的DEM数据，计算了实测点对应位置的地形要素(上游集水面积、局部坡度、坡度原点矩、河道曲率)，分析了河宽与地形要素间的相关关系，基于不同要素建立了6种河宽预测模型并进行了验证。研究得到以下结论:

(1)拉萨河流域河宽与上游集水面积呈良好的幂函数关系，确定性系数达到了0.92。虽然局部坡度不能很好地展现出河道地形的起伏变化，但坡度原点矩与河宽间明显存在两段不同的相关关系。以坡度原点矩表现的流域内主要河道的坡度从上游到下游先变陡后变缓，以30 m河宽为界拟合的两段相关关系的确定性系数均大于0.75。尽管实测河道曲率与河宽同样存在一定相关性，但是由于使用DEM难以准确提取河道曲率，因此不适用于河宽预测模型的建立研究。

(2)单独使用坡度原点矩建立的河宽预测模型预测效果不佳。使用综合了集水面积和坡度原点矩生成的河宽因子构建的预测模型整体较好。对比发现在集水面积A<2 100 km2的上游较窄河段使用以坡度原点矩为局地性因素的河宽因子模型预测效果最优，在集水面积A≥2 100 km2的下游河段使用集水面积模型预测效果最优。

(3)拉萨河流域精确河宽模型应同时采用以坡度原点矩为局地性因素的河宽因子模型和集水面积模型，在上、下游分别使用河宽因子模型和集水面积模型以达到最高精度。该精确综合模型的平均误差在上、下游均小于20%。此外，为追求简便性，全河道使用河宽因子模型同样不会造成太大误差。

考虑到拉萨河流域的Google Earth等影像全部为低水期的卫星影像和航拍影像，读取实测点时需要大量人眼判别经常性过水河道，人工判别的精度对结果易产生较大影响。因此，对青藏高原河宽的进一步研究有待无人机摄影测量等高精度雨季影像数据的推广。本研究发现的河道地形起伏自上游向下游呈先变陡后变缓的规律，与部分其他流域的已有研究成果一致[22]，表明该现象具有一定的普遍性。鉴于较多流域内集水面积、坡度原点矩(或河道比降)与河宽间的相关性最好[6,10-12]，本研究提出的基于河宽与坡度原点矩相关规律的拐点分段建模，根据集水面积模型、坡度原点矩模型及结合上述两个要素的河宽因子模型中选取最优者的河宽预测方法适合推广到其他流域。