CN108986222B - 无汊河道数字地形生成方法 - Google Patents

Abstract

本发明属于河流数值模拟技术领域，提供了一种无汊河道数字地形生成方法，包括(1)获取基础数据，(2)目标河段网格剖分，(3)横断面概化和(4)生成数字地形这四个步骤。该方法基于特征纵向控制线与较少实测横断面地形数据即可实现无汊河道数字地形生成，可解决现有技术为了保证网格化插值结果真实合理，需要提供覆盖模拟区域且分布均匀的大量地形数据而造成构建河网模型成本高、周期长和难度大的问题，有效降低了河道数字地形的生成难度与河网模型的构建成本。

Description

无汊河道数字地形生成方法

技术领域

本发明属于河流数值模拟技术领域，涉及一种无汊河道数字地形生成方法。

背景技术

在河流数值模拟中，网格节点的地形值能否反映真实地形非常重要，它直接决定着模拟结果的可靠性。很多水动力计算软件，例如Mike、SMS等，在河道数字地形插值中为保证网格化插值结果真实合理，均要求提供覆盖数模计算区域且分布均匀的大量地形数据。由于地形测量周期长，测量难度大，投入经费大，对于一般河道而言，在多数情况下并不采用对整个河道布设均匀测点的方式来进行地形量测，通常测得的只是具有相当布设间距的横断面地形数据。如果仍完全依靠常规软件自带的插值功能来构建河网模型，其地形精度势必较低，难于精准反映实际河床的起伏变化，也就不能保证数值模拟结果准确可信，因此有必要探寻一种基于实测横断面地形数据生成河道数字地形的可行方法。

发明内容

针对现有技术为保证网格化插值结果真实合理，需要提供覆盖模拟区域且分布均匀的大量地形数据而造成构建河网模型成本高、周期长和难度大的不足，本发明提出了一种基于特征纵向控制线与较少实测横断面地形数据的无汊河道数字地形生成方法，以有效降低河道数字地形的生成难度与河网模型的构建成本。

本发明提供的无汊河道数字地形生成方法，步骤如下：

(1)获取基础数据

选定需要进行数字地形生成的无汊目标河段，获取目标河段的包括河道边界线在内的特征纵向控制线上若干控制点的平面坐标数据，将位于最左侧与最右侧特征纵向控制线之间的纵向控制线记作中间特征纵向控制线；通过现场勘测获取目标河段的若干横断面的地形数据，横断面的地形数据包括横断面左、右端点，以及左、右端点之间的各测点的平面坐标数据和高程数据；将位于目标河段两端的横断面记作初始横断面和终止横断面；

(2)目标河段网格剖分

①计算初始横断面、终止横断面与最左侧、最右侧特征纵向控制线的交点的平面坐标，将初始横断面、终止横断面与最左侧、最右侧特征纵向控制线的交点分别记作初始交点和终止交点，将初始交点、终止交点以及介于初始交点与终止交点之间的部分称为有效部分，提取最左侧、最右侧特征纵向控制线的有效部分的控制点坐标，将位于最左侧、最右侧特征纵向控制线的有效部分的控制点记作有效控制点，假设最左侧与最右侧特征纵向控制线的有效控制点个数分别为N_左+1、N_右+1，将最左侧、最右侧纵向控制线的各有效控制点的坐标分别记作(x_左(i),y_左(i)),i＝1,2,3,…,N_左+1，(x_右(i),y_右(i)),i＝1,2,3,…,N_右+1，i＝1代表初始交点，计算最左侧特征纵向控制线的各有效控制点相对于最左侧初始交点的沿线累加距离L_左(i)i＝1,2,3,…,N_左+1，以及最右侧特征纵向控制线的各有效控制点相对于最右侧初始交点的沿线累加距离L_右(i)i＝1,2,3,…,N_右+1；

②采用定数等分法或定距等分法对最左侧特征纵向控制线的有效部分进行沿线剖分生成剖分节点，具体如下：

定数等分法：剖分份数取n_左，则沿线剖分步长s_左＝L_左(N_左+1)/n_左，剖分后将生成n_左+1个剖分节点，剖分过程如下：依据沿线累加距离将最左侧特征纵向控制线的有效部分分为N_左个累加距离区间[L_左(1),L_左(2)]，[L_左(2),L_左(3)]，…，[L_左(N_左),L_左(N_左+1)]，然后以最左侧特征纵向控制线的初始交点为起点，沿着最左侧特征纵向控制线的有效部分取沿线步进距离j·s_左,j依次取1,2,…n_左-1，根据j·s_左的值判断各沿线步进距离的终点落于哪一个累加距离区间中，当L_左(k)≤j·s_左≤L_左(k+1)时，则沿线步进距离j·s_左的终点落于累加距离区间[L_左(k),L_左(k+1)]中，k＝1,2,3,…,N_左，由式(Ⅰ)～(Ⅱ)计算出最左侧特征纵向控制线的有效部分的剖分节点坐标(x_左,节(j+1),y_左,节(j+1))，并记最左侧初始交点与最左侧终止交点的平面坐标分别为(x_{左,初始交点},y_{左,初始交点})＝(x_左,节(1),y_左,节(1))，(x_{左,终止交点},y_{左,终止交点})＝(x_左,节(n_左+1),y_左,节(n_左+1))；

定距等分法：定距等分距离取L_0,左，则沿线剖分步长s_左＝L_0,左，剖分后生成的剖分节点个数与定距等分距离取值有关，当L_左(N_左+1)/L_0,左的余数等于0时，相当于采用定数等分法进行剖分，剖分份数n_左＝L_左(N_左+1)/L_0,左，剖分后将生成n_左+1个剖分节点，剖分过程及剖分节点坐标的计算方法同定数等分法；当L_左(N_左+1)/L_0,左的余数不等于0时，剖分份数n_左＝[L_左(N_左+1)/L_0,左]+1，[L_左(N_左+1)/L_0,左]表示取L_左(N_左+1)/L_0,左的整数部分，剖分后将生成n_左+1个剖分节点，剖分过程及剖分节点坐标的计算方法同定数等分法；

③采用定数等分法或定距等分法对最右侧特征纵向控制线的有效部分进行沿线剖分生成剖分节点，由于是基于四边形网格对目标河段进行剖分，所以最右侧与最左侧特征纵向控制线的有效部分的剖分份数应相等，即n_右＝n_左，具体如下：

定数等分法：剖分份数为n_右，则沿线剖分步长s_右＝L_右(N_右+1)/n_右，剖分后将生成n_右+1个剖分节点，剖分过程如下：依据沿线累加距离将最右侧特征纵向控制线的有效部分分为N_右个累加距离区间[L_右(1),L_右(2)]，[L_右(2),L_右(3)]，…，[L_右(N_右),L_右(N_右+1)]，然后以最右侧特征纵向控制线的初始交点为起点，沿着最右侧特征纵向控制线的有效部分取沿线步进距离j·s_右,j依次取1,2,…n_右-1，根据j·s_右的值判断各沿线步进距离的终点落于哪一个累加距离区间中，当L_右(k)≤j·s_右≤L_右(k+1)时，则沿线步进距离j·s_右的终点落于累加距离区间[L_右(k),L_右(k+1)]中，k＝1,2,3,…,N_右，由式(Ⅲ)～(Ⅳ)计算出最右侧特征纵向控制线的有效部分的剖分节点坐标(x_右,节(j+1),y_右,节(j+1))，并记最右侧初始交点与最右侧终止交点的平面坐标分别为(x_{右,初始交点},y_{右,初始交点})＝(x_右,节(1),y_右,节(1))，(x_{右,终止交点},y_{右,终止交点})＝(x_右,节(n_右+1),y_右,节(n_右+1))；

定距等分法：为保证n_右＝n_左，最右侧特征纵向控制线的定距等分距离L_0,右应在区间[L_右(N_右+1)/n_右，L_右(N_右+1)/(n_右-1))内取值，沿线剖分步长s_右＝L_0,右，剖分后将生成n_右+1个剖分节点，剖分过程及剖分节点坐标的计算方法与步骤③中的定数等分法相同；

④依次连接最左侧与最右侧特征纵向控制线上对应的剖分节点，计算所连线段与各中间特征纵向控制线的交点的平面坐标，以实现对各中间特征纵向控制线的剖分，至此即完成对目标河段的四边形网格初步剖分；

⑤对各纵向网格线之间的横向线段进行剖分以减小各纵向网格线之间的横向间距，具体方法是参照步骤②和③中的定数等分法对各纵向网格线之间的横向线段进行剖分，计算剖分后的节点坐标，然后纵向依次连接该步骤剖分生成的对应的剖分节点，即完成对目标河段四边形网格初步剖分结果的横向加密；

(3)横断面概化

计算各横断面与各纵网格线的交点的平面坐标，然后基于横断面上各测点的平面坐标数据和高程数据，采用距离加权法对各横断面与各纵向网格线的交点进行高程插值，实现横断面的地形概化；

(4)生成数字地形

利用步骤(2)完成网格剖分后得到的网格节点的平面坐标与步骤(3)完成横断面地形概化的横断面地形数据，采用距离加权法沿纵向网格线对网格节点进行高程插值，得到所有网格节点的高程数据，即完成无汊河道数字地形的生成。

上述无汊河道数字地形生成方法的技术方案中，所述的距离加权法是基于相近相似原理的插值方法，即两点离得越近，其性质就越相似。该方法利用待插值点与样本点间的距离构造插值权重，离待插值点越近的样本点被赋予的权重越大。如下所示，假设两样本点的高程分别为z₁、z₂，待插值点与两样本点间的距离分别为d₁、d₂，待插值点的高程z₀则由右式(Ⅴ)计算。

具体地，对于步骤(3)中横断面概化时的各横断面与各纵向网格线的交点进行高程插值而言，此处的距离d₁、d₂是指直线距离，对于步骤(4)中的网格节点进行高程插值而言，此处的距离d₁、d₂是指纵向沿线距离。

步骤(3)中基于横断面上各测点的平面坐标数据和高程数据，采用距离加权法进行高程插值时，的方法如下：

假设待进行高程插值的点为A点，与A点处于同一横断面上且与A点直线距离最近的两个测点为B点和C点，B点、C点的高程分别为z_B，z_C，A点到B点、C点的距离分别为d_AB、d_AC，则A点的高程

步骤(4)中采用距离加权法沿纵向网格线对网格节点进行高程插值的方法如下：

假设待进行高程插值的点为D点，与D点纵向沿线距离最接近的两个已知高程的点为E点和F点，E点和F点为概化后的横断面上的点，D点位于E点和F点所在的概化后的横断面之间，E点、F点的高程分别为z_E，z_F，D点到E点、F点的纵向沿线距离分别为d_DE、d_DF，则D点的高程

上述无汊河道数字地形生成方法的技术方案中，通过对卫星图片、遥感图像、河势图解译或现场勘测获取目标河段的特征纵向控制线上若干控制点的平面坐标数据。

上述无汊河道数字地形生成方法的技术方案中，特征纵向控制线还包括深泓线、滩槽分界线及水边线，在需要考虑控导工程与生产堤等工程边界对河道数字地形的影响的情况下，还需提供相应的工程边界位置坐标。

上述无汊河道数字地形生成方法的技术方案中，步骤(2)②中，剖分份数n_左的取值越大，步骤(2)④中完成对目标河段的四边形网格初步剖分得到的网格的纵向间距就越小，步骤(2)⑤中采用定数等分法对各纵向网格线之间的横向线段进行剖分的剖分份数越大，步骤(2)⑤中完成对目标河段四边形网格初步剖分结果的横向加密形成的网格的横向间距就越小，最终剖分得到的四边形网格就越密集。

上述无汊河道数字地形生成方法的技术方案中，通常能够获取的特征纵向控制线的数量是较为有限的，通过对卫星图片、遥感图像、河势图解译这些技术手段获取的特征纵向控制线通常只包括河道边界线、深泓线、滩槽分界线、水边线，由于特征纵向控制线数目较少，故步骤(2)④完成对目标河段的四边形网格初步剖分后得到的纵向网格线的横向间距较大。为了解决各纵向网格线横向间距过大的问题，本发明采用步骤(2)⑤对各纵向网格线之间的横向线段再进行剖分，为保证剖分网格横向疏密均匀，须采用定数等分法对各纵向网格线之间的横向线段进行剖分，不同相邻纵向网格线间的横向定数剖分份数可以不同，而同一相邻纵向网格线间的横向定数剖分份数必须相同，定数剖分份数的具体取值应根据各相邻纵向网格线的横向疏密程度确定。当两相邻纵向网格线之间的间距较大时，则横向线段的定数剖分份数取值应较大；当两相邻纵向网格线之间的间距较小时，则横向线段的定数剖分份数取值应较小。剖分后再沿纵向依次连接各相应的新增节点，对各相邻纵向网格线之间的区域进行补线操作，即完成对河道区域四边形网格初步剖分结果的横向加密。

上述无汊河道数字地形生成方法的技术方案中，横断面和特征纵向控制线数量越多、横断面上的测点布置越密集，步骤(2)②、⑤中采用的剖分份数越大，最终生成的数字地形的精确度就越高。但在实际用时获取的特征纵向控制线的数量是较为有限的，因此，为了提高生成的数字地形的精度，最好是提供尽量多的横断面地形数据，增加横断面地形数据，同时增加步骤(2)②、⑤中采用的剖分份数。在横断面和特征纵向控制线数量确定的情况下，步骤(2)②、⑤中的剖分份数取值越大，目标河段网格剖分后形成的网格就越精细，最终生成的数字地形的精度就越高。剖分份数的增加会造成计算量的增加，实际应用时，剖分份数根据实际工程对数字地形的精度要求、计算量以及计算设备计算能力的情况综合考虑确定。

上述无汊河道数字地形生成方法的技术方案中，步骤(2)①中，计算初始横断面、终止横断面与最左侧、最右侧纵向控制线的交点的平面坐标的方法是将与横断面相交处的最左侧或最右侧特征纵向控制线看作一小段线段，先计算该线段的直线方程，以及与该线段相交的横断面的直线方程，然后计算该线段与横断面的交点坐标即可。

上述无汊河道数字地形生成方法的技术方案中，步骤(3)横断面概化是否合理决定着生成的河道数字地形精度，最理想的情况是各横断面的概化形态与其实际地形完全一致，也就是说相同水位条件下概化后的横断面与相应的横断面实际地形具有相同的过流面积。而实际概化效果的优劣与网格剖分的横向疏密直接相关，当网格剖分的横向间距较小时，则基于距离加权法进行高程插值概化得到的河床横断面形态就与实际地形就更为接近。若需要对实测横断面的概化合理度作出评判，可分别计算实际横断面与概化横断面在一系列水位条件下的过流面积参数(如过流面积、湿周、水力半径、水面宽度与平均水深)进行比较。

借助计算机编程(如Fortran、Matlab等)，将本发明所述方法生成的数字地形数据按照Tecplot、Mike、SMS、南方Cass、Auto-CAD等软件的可读文本格式或批量绘图命令序列进行整合，可实现其数据交换通道的完美搭建，继而充分发挥所述软件的独特优势对所生成的河道数字地形进行多形式可视化与深层次再利用。

与现有技术相比，本发明具有以下有益效果：

1.本发明提供了一种无汊河道数字地形生成方法，该方法通过将河道边界线、深泓线、滩槽分界线及水边线等特征地形界线纳入河道区域的四边形网格剖分过程，使所剖分的网格能够较好地适应河道边界变化，管控地形纵向走势，反映横断面形态变率；基于距离加权法利用实测横断面地形数据对实测横断面与纵向网格线交点进行高程插值，实现对实测横断面地形的合理概化；利用概化横断面数据基于距离加权法沿纵向网格线进行高程内插以获得所有网格节点的高程数据。该方法物理概念清晰，执行简单高效，普遍适用于以实测横断面数据为基础的天然河道数字地形生成。

2.本发明提供的无汊河道数字地形生成方法是一种基于特征纵向控制线与实测横断面地形数据生成河道数字地形的可行方法，可解决现有技术为了保证网格化插值结果真实合理，需要提供覆盖模拟区域且分布均匀的大量地形数据而造成构建河网模型成本高、周期长和难度大的问题，有效降低河道数字地形的生成难度与河网模型的构建成本。

附图说明

图1是实施例1中目标河段的特征纵向控制线和横断面分布图。

图2是实施例1中目标河段的四边形网格初步剖分结果。

图3是实施例1中对目标河段四边形网格初步剖分结果进行横向加密后的结果。

图4是实施例1中目标河段的第

号～第

号横断面与各纵网格线的交点示意图。

图5是实施例1中

号横断面的实际地形与概化后的地形比较图。

图6是实施例1中生成的河道散点数字地形图。

图7是实施例2中Tecplot绘图文件示意图。

图8是实施例2中基于Tecplot的河床形态。

图9是实施例2中Mike四边形网格文件示意图。

图10是实施例2中Mike四边形网格文件节点与单元的编码规则图。

图11是实施例2中Mike四边形网格单元构成节点的编号计算示意图。

图12是实施例2中基于Mike的地形云图。

图13是实施例2中基于Mike的数模局部流场分布图。

图14是实施例2中SMS四边形网格文件示意图。

图15是实施例2中SMS四边形网格文件节点与单元的编码规则图。

图16是实施例2中SMS四边形网格单元构成节点的编号计算示意图。

图17是实施例2中基于SMS的着色河网。

图18是实施例2中基于SMS的数模流场分布图。

图19是实施例2中南方Cass可读文件示意图。

图20是实施例2中基于南方Cass的河道等高线图。

图21是实施例2中Auto-CAD绘图命令序列示意图。

图22是实施例2中基于Auto-CAD的三维网格图。

具体实施方式

下面结合附图通过实施例对本发明提供的无汊河道数字地形生成方法作进一步说明。有必要指出，以下实施例只用于对本发明作进一步说明，不能理解为对本发明保护范围的限制，所属领域技术人员根据上述发明内容，对本发明做出一些非本质的改进和调整进行具体实施，仍属于本发明的保护范围。

实施例1

本实施例中，详细说明无汊河道数字地形生成方法，步骤如下：

(1)获取基础数据

选定需要进行数字地形生成的无汊目标河段，目标河段的特征纵向控制线和横断面分布情况如图1所示。

通过卫星图片获取目标河段特征纵向控制线的地形数据，具体是将一张高清晰度的卫星图片导入CAD软件中，描出河道边界线和槽滩分界线，得到4条特征纵向控制线，以及这4条纵向控制线上若干控制点的平面坐标数据，将位于河道左边界线(最左侧特征纵向控制线)与河道右边界线(最右侧特征纵向控制线)之间的特征纵向控制线(也就是滩槽分界线)称为中间特征纵向控制线，左侧与右侧的判定基于面向水流的方向进行。

通过实测获取目标河段若干横断面地形数据，本实施例中共包括21个横断面，横断面地形数据包括横断面左、右端点，以及左、右端点之间各测点的平面坐标数据和高程数据；将位于目标河段最上游的横断面记作初始横断面、最下游的横断面记作终止横断面，按照从初始横断面到终止横断面的方向，依次对各横断面编号为①号～

号横断面。

(2)目标河段网格剖分

①计算初始横断面、终止横断面与最左侧、最右侧特征纵向控制线的交点的平面坐标，将初始横断面、终止横断面与最左侧、最右侧特征纵向控制线的交点分别记作初始交点和终止交点，将初始交点、终止交点以及介于初始交点与终止交点之间的部分称为有效部分，提取最左侧、最右侧特征纵向控制线的有效部分的控制点坐标，将位于最左侧、最右侧特征纵向控制线的有效部分的控制点记作有效控制点，假设最左侧与最右侧特征纵向控制线的有效控制点个数分别为N_左+1、N_右+1，将最左侧、最右侧纵向控制线的各有效控制点的坐标分别记作(x_左(i),y_左(i)),i＝1,2,3,…,N_左+1，(x_右(i),y_右(i)),i＝1,2,3,…,N_右+1，i＝1代表初始交点，即最左侧、最右侧特征纵向控制线的有效部分分别由N_左、N_右条顺次相连的线段组成。

计算最左侧特征纵向控制线的各有效控制点相对于最左侧初始交点的沿线累加距离L_左(i)i＝1,2,3,…,N_左+1，以及最右侧特征纵向控制线的各有效控制点相对于最右侧初始交点的沿线累加距离L_右(i)i＝1,2,3,…,N_右+1。

实际上最左/右侧特征纵向控制线的终止交点相对于相应初始交点的沿线累加距离就等于其有效部分总长度，而后一沿线累加距离与前一沿线累加距离的差值则为两对应控制点的沿线距离。例如，L_右(1)＝0，L_右(2)表示最右侧特征纵向控制线上的第二个有效控制点到最右侧初始交点的距离，L_右(3)表示最右侧特征纵向控制线上的第二个有效控制点到最右侧初始交点的距离与最右侧特征纵向控制线上的第三个有效控制点到第二个有效控制点的距离之和，L_右(3)-L_右(2)表示最右侧特征纵向控制线上的第三个有效控制点与第二个有效控制点之间的距离。

②采用定数等分法对最左侧特征纵向控制线的有效部分进行沿线剖分生成剖分节点，具体如下：

剖分份数取n_左，则沿线剖分步长s_左＝L_左(N_左+1)/n_左，剖分后将生成n_左+1个剖分节点，剖分过程如下：依据沿线累加距离将最左侧特征纵向控制线的有效部分分为N_左个累加距离区间[L_左(1),L_左(2)]，[L_左(2),L_左(3)]，…，[L_左(N_左),L_左(N_左+1)]，然后以最左侧特征纵向控制线的初始交点为起点，沿着最左侧特征纵向控制线的有效部分取沿线步进距离j·s_左,j依次取1,2,…n_左-1，根据j·s_左的值判断各沿线步进距离的终点落于哪一个累加距离区间中，当L_左(k)≤j·s_左≤L_左(k+1)时，则沿线步进距离j·s_左的终点落于累加距离区间[L_左(k),L_左(k+1)]中，k＝1,2,3,…,N_左，由式(Ⅰ)～(Ⅱ)计算出最左侧特征纵向控制线的有效部分的剖分节点坐标(x_左,节(j+1),y_左,节(j+1))，并记最左侧初始交点与最左侧终止交点的平面坐标分别为(x_{左,初始交点},y_{左,初始交点})＝(x_左,节(1),y_左,节(1))，(x_{左,终止交点},y_{左,终止交点})＝(x_左,节(n_左+1),y_左,节(n_左+1))；

③采用定数等分法对最右侧特征纵向控制线的有效部分进行沿线剖分生成剖分节点，由于是基于四边形网格对目标河段进行剖分，所以最右侧与最左侧特征纵向控制线的有效部分的剖分份数应相等，即n_右＝n_左，具体如下：

剖分份数为n_右，则沿线剖分步长s_右＝L_右(N_右+1)/n_右，剖分后将生成n_右+1个剖分节点，剖分过程如下：依据沿线累加距离将最右侧特征纵向控制线的有效部分分为N_右个累加距离区间[L_右(1),L_右(2)]，[L_右(2),L_右(3)]，…，[L_右(N_右),L_右(N_右+1)]，然后以最右侧特征纵向控制线的初始交点为起点，沿着最右侧特征纵向控制线的有效部分取沿线步进距离j·s_右,j依次取1,2,…n_右-1，根据j·s_右的值判断各沿线步进距离的终点落于哪一个累加距离区间中，当L_右(k)≤j·s_右≤L_右(k+1)时，则沿线步进距离j·s_右的终点落于累加距离区间[L_右(k),L_右(k+1)]中，k＝1,2,3,…,N_右，由式(Ⅲ)～(Ⅳ)计算出最右侧特征纵向控制线的有效部分的剖分节点坐标(x_右,节(j+1),y_右,节(j+1))，并记最右侧初始交点与最右侧终止交点的平面坐标分别为(x_{右,初始交点},y_{右,初始交点})＝(x_右,节(1),y_右,节(1))，(x_{右,终止交点},y_{右,终止交点})＝(x_右,节(n_右+1),y_右,节(n_右+1))；

④依次连接最左侧与最右侧特征纵向控制线上对应的剖分节点，计算所连线段与各中间特征纵向控制线的交点的平面坐标，以实现对各中间特征纵向控制线的剖分，至此即完成对目标河段的四边形网格初步剖分，结果如图2所示。

从图2可以发现，步骤(2)④完成对目标河段的四边形网格初步剖分后得到的纵向网格线的横向间距较大，这是由于特征纵向控制线的数量较少的原因造成的。为了解决各纵向网格线横向间距过大的问题，需要采用步骤⑤对各纵向网格线之间的横向线段再进行剖分。

⑤对各纵向网格线之间的横向线段进行剖分以减小各纵向网格线之间的横向间距，具体方法是参照步骤②和③中的定数等分法对各纵向网格线之间的横向线段进行剖分，定数剖分份数的具体取值应根据各相邻纵向网格线的横向疏密程度确定，当两相邻纵向网格线之间的间距较大时，则横向线段的定数剖分份数取值应较大，当两相邻纵向网格线之间的间距较小时，则横向线段的定数剖分份数取值应较小，计算剖分后的节点坐标，然后纵向依次连接该步骤剖分生成的对应的剖分节点，即完成对目标河段四边形网格初步剖分结果的横向加密，结果如图3所示。

(3)横断面概化

计算各横断面与各纵网格线的交点的平面坐标，目标河段的第

号横断面～第

号横断面与各纵网格线的交点如图4所示，然后基于横断面上各测点的平面坐标数据和高程数据，采用距离加权法对各横断面与各纵向网格线的交点进行高程插值，实现横断面的地形概化。

该步骤中基于横断面上各测点的平面坐标数据和高程数据，采用距离加权法对各横断面与各纵向网格线的交点进行高程插值的方法如下：

假设待进行高程插值的点为A点，与A点处于同一横断面上且与A点直线距离最近的两个测点为B点和C点，B点、C点的高程分别为z_B，z_C，A点到B点、C点的距离分别为d_AB、d_AC，则A点的高程

经过本步骤的横断面概化后，

号横断面的实际地形与概化后的地形比较图如图5所示，由图5可知，该横断面的概化形态与其实际地形高度吻合，说明该步骤横断面概化准确合理。

(4)生成数字地形

利用步骤(2)完成网格剖分后得到的网格节点的平面坐标与步骤(3)完成横断面地形概化的横断面地形数据，采用距离加权法沿纵向网格线对网格节点进行高程插值，得到所有网格节点的高程数据，即完成无汊河道数字地形的生成，结果如图6所示。

该步骤中，采用距离加权法沿纵向网格线对网格节点进行高程插值的方法如下：

假设待进行高程插值的点为D点，与D点纵向沿线距离最接近的两个已知高程的点为E点和F点，E点和F点为概化后的横断面上的点，D点位于E点和F点所在的概化后的横断面之间，E点、F点的高程分别为z_E，z_F，D点到E点、F点的纵向沿线距离分别为d_DE、d_DF，则D点的高程

实施例2

本实施例中，举例说明采用本方法生成的数字地形的延拓性应用。

借助计算机编程(如Fortran、Matlab等)将实施例1生成的目标河段的数字地形按照Tecplot、Mike、SMS、南方Cass、Auto-CAD等软件的可读文本格式或批量绘图命令序列进行整合，实现数据交换通道的完美搭建，充分发挥各软件独特优势对所生成的河道数字地形进行多形式可视化与深层次再利用。

(1)Tecplot

借助计算机编程将实施例1生成的河道数字地形保存为如图7所示扩展名为“plt”的有序数据文件，就可基于Tecplot强大的图形可视化功能对河道形态进行生动显示，如图8所示。在所述有序数据文件中，全部地形数据隶属一个区域(zone)；i、j、k为数据组维数下标(i代表纵向网格线数目，j代表横向网格线数目)，i、j大于1而k等于1表示地形数据构成一数据点个数为i·j的二维数组；地形数据间以空格分隔，依次为横坐标、纵坐标和高程，与第二行(variables＝"X","Y","Z")变量名的指定顺序一致；第一行(title＝”river bedform”)用于指定图名。

(2)Mike

Mike是由丹麦DHI公司研发的水流模拟组件，集合了目前广泛使用的Mike11和Mike21。Mike21适用于宏观上的流域控制性工程规模论证分析、流域洪水调度研究以及微观水流模拟等领域，其常用网格类型中就包括四边形网格。Mike21四边形网格文件的扩展名为“mesh”，其内部数据包括节点标题行、节点行、单元格标题行与单元格行，如图9所示。节点标题行又分为整数形式的条目类型、整数形式的条目单位、节点数量和投影类型字符串，条目类型为“高程”，其整数形式为“100079”；条目单位即高程单位，整数形式的“1000”代表保存在Z坐标的高程数值其单位为m；其后的整数“1846”为节点数量；最后的字符串“NON-UTM”为投影类型。每一节点行代表一个节点，节点行的行数与节点标题行中的节点数量相同，每行节点信息包括节点编号、X、Y、Z与边界代码，边界代码“0”代表内部节点，“1”代表水陆边界，“2”代表进口边界，“3”代表出口边界。单元格标题行的三个数字分别表示单元格数量、单个单元格的最大节点数量与单元格类型代码(“25”代表四边形单元格)。每一单元格行代表一个单元，单元格行的总行数与单元格标题行中定义的单元格数量相同，每行单元格信息包括单元格编号和构成单元格的节点编号。

为使实施例1生成的河道数字地形可用于Mike21数模计算，必须借助计算机编程将其以Mike的可读网格格式保存。在具体的网格转化过程中，节点与单元的编码可遵循图10所示规则，图中填充单元内的值为单元编号，填充单元四角处的值为节点编号，均从河道进口至出口沿纵向进行，横向输出次序则由河道左岸开始到河道右岸结束。

假设实施例1生成的河道数字地形的纵向网格线数目为m，横向网格线数目为n，纵向网格线编号顺序由左到右，横向网格线编号顺序从进口至出口，如图10所示。当将各纵向网格线分别视为一行，各横向网格线分别视为一列，行号以i标记，列号以j标记，单元以其四顶点的最小行号、最小列号组合(i,j)标记时，则Mike四边形网格文件中的单元编号及构成该单元的节点编号与横纵向网格线编号之间存在一一对应关系。如图11左侧所示，当已知Mike网格单元的行列号组合标记为(i,j)(其中i＝1,2,…,m-1；j＝1,2,…,n-1)时，便可计算出该单元的单元编号为(i-1)·(n-1)+j，如果再将构成该单元的各节点沿逆时针方向分别记为①、②、③、④，则各节点编号就可由图11右侧各式算得；若已知某单元的单元编号为N，亦可计算出该单元的行列号组合标记(i,j)，其中i取不小于N/(n-1)的最小整数，j等于N-(i-1)·(n-1)，当计算出i、j后，代入图11右侧各式即可求出相应的节点编号。

当借助编程语言基于上述规则生成Mike的可读网格文件后，即可将其导入Mike21水动力模块进行数模计算，其中网格地形云图与数模流场见图12、13。

(3)SMS

地表水模拟系统(Surface Water Modeling System)简称SMS，是美国陆军工程兵水利工程实验室(United States Army Corps of Engineers Hydraulics Laboratory)和扬·伯明翰大学(Brigham Young University)等合作开发的商业软件，其四边形网格文件扩展名为“2dm”，如图14所示，内部数据主要包括单元格行、节点行与开边界节点串(进出口边界)。单元格行以“E8Q”开头，其后依次为单元格编号、构成单元格的节点编号(SMS中的一个四边单元由八个节点构成-四个顶点与四边中点)与材质编号。节点行以“ND”开头，其后依次为节点编号、节点横纵坐标与高程。开边界节点串以“NS”开头，其后为组成节点串的各节点编号，编号顺序一般从右岸至左岸并以负编号标记结束。

为使实施例1生成的河道数字地形可用于SMS数模计算，必须借助计算机编程将其以SMS的可读网格格式保存。然由于SMS可读网格文件中的一个四边形单元包含八节点(四边形单元的四个顶点与四边中点)，而实施例1剖分得到的仅是四边形网格单元的顶点坐标，因此在进行网格转化前还需基于各四边形网格单元的顶点坐标计算出相应的各边中点坐标。在具体的网格转化过程中，SMS四边形网格文件中的节点与单元编码可遵循图15所示规则(图中填充单元内的值为单元编号，填充单元周围的值为节点编号)：先对各四边形网格单元的顶点进行编码，然后对横向边中点编码，最后对纵向边中点编码；四边形网格单元的顶点与横向边中点编码均从河道进口至出口沿着纵向网格线进行，横向输出次序则由河道左岸开始到河道右岸结束；四边形网格单元的纵向边中点编码从河道左岸至右岸沿着横向网格线进行，纵向输出次序则由河道进口开始到河道出口结束；相邻两纵向网格线间的网格单元编码从河道进口开始到河道出口沿结束，横向输出次序则沿横向网格线方向由左至右。

假设实施例1生成的河道数字地形的纵向网格线数目为m，横向网格线数目为n，纵向网格线编号顺序由左到右，横向网格线编号顺序从进口至出口(如图15所示)。当将各纵向网格线分别视为一行，各横向网格线分别视为一列，行号以i标记，列号以j标记，单元以其四顶点的最小行号、最小列号组合(i,j)标记时，则SMS四边形网格文件中的单元编号及构成该单元的节点编号与横纵向网格线编号之间存在一一对应关系。如图16左侧所示，当已知SMS网格单元的行列号组合标记为(i,j)，其中i＝1,2,…,m-1；j＝1,2,…,n-1时，便可计算出该单元的单元编号为(i-1)·(n-1)+j，如果再将构成该单元的各节点沿逆时针方向分别记为①、②、③、④、⑤、⑥、⑦、⑧，则各节点编号就可由图16右侧各式算得；若已知某单元的单元编号为N，亦可计算出该单元的行列号组合标记(i,j)，其中i取不小于N/(n-1)的最小整数，j等于N-(i-1)·(n-1)，当计算出i、j后，代入图16右侧各式即可求出相应的节点编号。

当借助编程语言基于上述规则生成SMS的可读网格文件后，即可将其导入SMS二维水沙输移模块进行数模计算，其河网着色图与数模流场见图17、18。

(4)南方Cass

借助计算机编程将生成的河道数字地形散点以如图19所示的文本文件格式保存，则可基于南方Cass生成河道的地形等高线，见图20。南方Cass可读文件的扩展名为“dat”，其内部数据分为四列，分别为散点编号、X、Y、Z，其中第一、二两列以两个逗号分隔，其余各列以一个逗号分隔。

(5)Auto-CAD

借助计算机编程将河网数字地形分别沿纵向与横向输出成如图21所示的直线批量绘图命令序列，则可通过Auto-CAD进行河道地形的三维网格显示。所述批量绘图命令序列由若干命令单元组成，各命令行均不含空格(空格行除外)，地形河网各节点的横纵坐标与高程之间用逗号分隔；一个命令单元起始于“line”而终止于空格行，若将一个直线绘图命令序列复制、粘贴至Auto-CAD命令行将绘制出一条三维纵向网格线或横向网格线；当将全部直线绘图命令序列复制、粘贴至Auto-CAD命令行，待Auto-CAD运行完毕后便可显示河道的三维河网(见图22)。

Claims (4)

1.一种无汊河道数字地形生成方法，其特征在于步骤如下：

(1)获取基础数据

选定需要进行数字地形生成的无汊目标河段，获取目标河段的包括河道边界线在内的特征纵向控制线上若干控制点的平面坐标数据，将位于最左侧与最右侧特征纵向控制线之间的纵向控制线记作中间特征纵向控制线；通过现场勘测获取目标河段的若干横断面的地形数据，横断面的地形数据包括横断面左、右端点，以及左、右端点之间的各测点的平面坐标数据和高程数据；将位于目标河段两端的横断面记作初始横断面和终止横断面；

(2)目标河段网格剖分

①计算初始横断面、终止横断面与最左侧、最右侧特征纵向控制线的交点的平面坐标，将初始横断面、终止横断面与最左侧、最右侧特征纵向控制线的交点分别记作初始交点和终止交点，将初始交点、终止交点以及介于初始交点与终止交点之间的部分称为有效部分，提取最左侧、最右侧特征纵向控制线的有效部分的控制点坐标，将位于最左侧、最右侧特征纵向控制线的有效部分的控制点记作有效控制点，假设最左侧与最右侧特征纵向控制线的有效控制点个数分别为N_左+1、N_右+1，将最左侧、最右侧纵向控制线的各有效控制点的坐标分别记作(x_左(i),y_左(i)),i＝1,2,3,…,N_左+1，(x_右(i),y_右(i)),i＝1,2,3,…,N_右+1，i＝1代表初始交点，计算最左侧特征纵向控制线的各有效控制点相对于最左侧初始交点的沿线累加距离L_左(i)i＝1,2,3,…,N_左+1，以及最右侧特征纵向控制线的各有效控制点相对于最右侧初始交点的沿线累加距离L_右(i)i＝1,2,3,…,N_右+1；

②采用定数等分法或定距等分法对最左侧特征纵向控制线的有效部分进行沿线剖分生成剖分节点，具体如下：

定数等分法：剖分份数取n_左，则沿线剖分步长s_左＝L_左(N_左+1)/n_左，剖分后将生成n_左+1个剖分节点，剖分过程如下：依据沿线累加距离将最左侧特征纵向控制线的有效部分分为N_左个累加距离区间[L_左(1),L_左(2)]，[L_左(2),L_左(3)]，…，[L_左(N_左),L_左(N_左+1)]，然后以最左侧特征纵向控制线的初始交点为起点，沿着最左侧特征纵向控制线的有效部分取沿线步进距离j·s_左,j依次取1,2,…n_左-1，根据j·s_左的值判断各沿线步进距离的终点落于哪一个累加距离区间中，当L_左(k)≤j·s_左≤L_左(k+1)时，则沿线步进距离j·s_左的终点落于累加距离区间[L_左(k),L_左(k+1)]中，k＝1,2,3,…,N_左，由式(Ⅰ)～(Ⅱ)计算出最左侧特征纵向控制线的有效部分的剖分节点坐标(x_左,节(j+1),y_左,节(j+1))，并记最左侧初始交点与最左侧终止交点的平面坐标分别为(x_{左,初始交点},y_{左,初始交点})＝(x_左,节(1),y_左,节(1))，(x_{左,终止交点},y_{左,终止交点})＝(x_左,节(n_左+1),y_左,节(n_左+1))；

定距等分法：定距等分距离取L_0,左，则沿线剖分步长s_左＝L_0,左，剖分后生成的剖分节点个数与定距等分距离取值有关，当L_左(N_左+1)/L_0,左的余数等于0时，相当于采用定数等分法进行剖分，剖分份数n_左＝L_左(N_左+1)/L_0,左，剖分后将生成n_左+1个剖分节点，剖分过程及剖分节点坐标的计算方法同定数等分法；当L_左(N_左+1)/L_0,左的余数不等于0时，剖分份数n_左＝[L_左(N_左+1)/L_0,左]+1，[L_左(N_左+1)/L_0,左]表示取L_左(N_左+1)/L_0,左的整数部分，剖分后将生成n_左+1个剖分节点，剖分过程及剖分节点坐标的计算方法同定数等分法；

③采用定数等分法或定距等分法对最右侧特征纵向控制线的有效部分进行沿线剖分生成剖分节点，由于是基于四边形网格对目标河段进行剖分，所以最右侧与最左侧特征纵向控制线的有效部分的剖分份数应相等，即n_右＝n_左，具体如下：

定数等分法：剖分份数为n_右，则沿线剖分步长s_右＝L_右(N_右+1)/n_右，剖分后将生成n_右+1个剖分节点，剖分过程如下：依据沿线累加距离将最右侧特征纵向控制线的有效部分分为N_右个累加距离区间[L_右(1),L_右(2)]，[L_右(2),L_右(3)]，…，[L_右(N_右),L_右(N_右+1)]，然后以最右侧特征纵向控制线的初始交点为起点，沿着最右侧特征纵向控制线的有效部分取沿线步进距离j·s_右,j依次取1,2,…n_右-1，根据j·s_右的值判断各沿线步进距离的终点落于哪一个累加距离区间中，当L_右(k)≤j·s_右≤L_右(k+1)时，则沿线步进距离j·s_右的终点落于累加距离区间[L_右(k),L_右(k+1)]中，k＝1,2,3,…,N_右，由式(Ⅲ)～(Ⅳ)计算出最右侧特征纵向控制线的有效部分的剖分节点坐标(x_右,节(j+1),y_右,节(j+1))，并记最右侧初始交点与最右侧终止交点的平面坐标分别为(x_{右,初始交点},y_{右,初始交点})＝(x_右,节(1),y_右,节(1))，(x_{右,终止交点},y_{右,终止交点})＝(x_右,节(n_右+1),y_右,节(n_右+1))；

定距等分法：为保证n_右＝n_左，最右侧特征纵向控制线的定距等分距离L_0,右应在区间[L_右(N_右+1)/n_右，L_右(N_右+1)/(n_右-1))内取值，沿线剖分步长s_右＝L_0,右，剖分后将生成n_右+1个剖分节点，剖分过程及剖分节点坐标的计算方法与步骤③中的定数等分法相同；

④依次连接最左侧与最右侧特征纵向控制线上对应的剖分节点，计算所连线段与各中间特征纵向控制线的交点的平面坐标，以实现对各中间特征纵向控制线的剖分，至此即完成对目标河段的四边形网格初步剖分；

⑤对各纵向网格线之间的横向线段进行剖分以减小各纵向网格线之间的横向间距，具体方法是参照步骤②和③中的定数等分法对各纵向网格线之间的横向线段进行剖分，计算剖分后的节点坐标，然后纵向依次连接该步骤剖分生成的对应的剖分节点，即完成对目标河段四边形网格初步剖分结果的横向加密；

(3)横断面概化

计算各横断面与各纵网格线的交点的平面坐标，然后基于横断面上各测点的平面坐标数据和高程数据，采用距离加权法对各横断面与各纵向网格线的交点进行高程插值，实现横断面的地形概化；

(4)生成数字地形

利用步骤(2)完成网格剖分后得到的网格节点的平面坐标与步骤(3)完成横断面地形概化的横断面地形数据，采用距离加权法沿纵向网格线对网格节点进行高程插值，得到所有网格节点的高程数据，即完成无汊河道数字地形的生成。

2.根据权利要求1所述无汊河道数字地形生成方法，其特征在于，步骤(3)中基于横断面上各测点的平面坐标数据和高程数据，采用距离加权法对各横断面与各纵向网格线的交点进行高程插值的方法如下：

假设待进行高程插值的点为A点，与A点处于同一横断面上且与A点直线距离最近的两个测点为B点和C点，B点、C点的高程分别为z_B，z_C，A点到B点、C点的距离分别为d_AB、d_AC，则A点的高程

3.根据权利要求1所述无汊河道数字地形生成方法，其特征在于，步骤(4)中采用距离加权法沿纵向网格线对网格节点进行高程插值的方法如下：

假设待进行高程插值的点为D点，与D点纵向沿线距离最接近的两个已知高程的点为E点和F点，E点和F点为概化后的横断面上的点，D点位于E点和F点所在的概化后的横断面之间，E点、F点的高程分别为z_E，z_F，D点到E点、F点的纵向沿线距离分别为d_DE、d_DF，则D点的高程

4.根据权利要求1至3中任一权利要求所述无汊河道数字地形生成方法，其特征在于，通过对卫星图片、遥感图像、河势图解译或现场勘测获取目标河段的特征纵向控制线上若干控制点的平面坐标数据。

Images