CN114757036A - 一种基于逐次逼近法的河道断面自动插补方法 - Google Patents

Abstract

本发明公开了属于河道水动力数值模拟技术领域的一种基于逐次逼近法的河道断面自动插补方法。包括以下步骤：步骤1：以观测断面水位和流量水文要素的纳什系数和最大为目标函数，以水位或流量误差和断面间距为约束条件，建立基于一维非恒定流数学模型的河道断面自动插补优化模型；步骤2：确定河道断面自动插补优化模型的参数取值范围；步骤3：采用以形定位模式或以位定形模式的逐次逼近法对河道断面自动插补优化模型进行求解。本发明实现了对缺少地形资料的山区河道插补断面的个数、位置、形状和河道糙率等参数的自动优化；可避免人工调试的主观任意性和不确定性，提高模拟精度。

Description

一种基于逐次逼近法的河道断面自动插补方法

技术领域

本发明涉及河道水动力数值模拟技术领域，尤其涉及一种基于逐次逼近法的河道断面自动插补方法。

背景技术

河道地形资料是河道水动力建模的重要依据，但常存在河道地形资料不足或资料粗糙等现象。为此，可充分利用有限的地形和水文资料，适当插补断面以更好地模拟河道的沿程水位和流量变化情况。

目前现有技术多是采用人工试算法，首先给定插补断面形状和位置的初始值，然后代入一维非恒定流数学模型进行计算，比较观测断面水位流量计算值与实测值之间的差距，再人工调整插补断面形状和位置参数的取值重新进行计算，重复这一过程直至观测断面水位流量计算值与实测值之间的差距满足精度要求，从而选出最优的插补断面形状和位置参数。采用人工试算法进行调参不仅费时费力、工作量大，且具有较大的经验性和任意性，精度难以保证。尤其是插补断面时还要考虑插补断面形状和位置之间的相互作用，计算难度较大，人工试算法不一定能够找出最优的插补断面形状和位置参数，使得观测断面水位流量计算值与实测值拟合误差较小。所以，有必要研究一种成熟有效的河道断面自动插补方法。

针对人工试算法在河道断面插补中的问题，基于逐次逼近算法，需要一种河道断面自动插补方法，在缺少地形资料的情况下，实现对河道断面的自动插补，得到最优的插补断面形状和位置参数。与人工调参相比，利用计算机智能调参不仅更省时省力，还可避免人工调试的主观任意性和不确定性，提高模拟精度。

发明内容

本发明的目的是提出一种基于逐次逼近法的河道断面自动插补方法，其特征在于，包括以下步骤：

步骤1：以观测断面水位与流量水文要素的纳什系数之和最大为目标函数，以水位、流量误差和断面间距为约束条件，建立基于一维非恒定流数学模型的河道断面自动插补优化模型；

步骤2：根据河道的长度、宽度和比降，确定河道断面自动插补优化模型的参数取值范围；所述参数包括河道糙率、插补断面个数、插补断面位置、插补断面形状参数，其中插补断面形状参数包括断面底宽、河底高程和边坡系数；

步骤3：根据基于逐次逼近法的并行求解方法，采用以形定位计算模式或以位定形计算模式对河道断面自动插补优化模型进行求解。

所述步骤1中河道断面自动插补优化模型为：

目标函数：

式中，i为观测断面的序号；N为观测断面总数；NSE(i)为观测断面i的纳什系数值，具体为：

式中，T为计算周期，y_t(i)为i断面的水位或流量实测值；

为i断面水位或流量实测序列的平均值；

为i断面水位或流量计算值；

约束条件：

式中：abs为绝对值函数，ε_i为第i个断面的最大水位或流量误差；

L_min≤L(k)≤L_max (4)

式中，L(k)为插补断面k与前一断面的断面间距；L_max为给定的最大断面间距；L_min为给定的最小断面间距，取为0。

所述步骤3中的基于逐次逼近法的并行求解方法具体包括以下步骤：

步骤S1：对参数进行离散，设插补断面个数为k₁，插补断面位置的个数为k₂，插补断面形状中断面底宽的个数为k₃、河底高程的个数为k₄、边坡系数的个数为k₅，河道糙率的个数为k₆，则共有

个参数组合；

步骤S2：固定插补断面个数、插补断面位置和河道糙率的取值，仅将各插补断面形状中的

个参数组合代入基于一维非恒定流数学模型的河道断面自动插补优化模型进行逐次寻优；

步骤S3：结合针对单主机多CPU并行计算而设计的OpenMP编译处理技术对基于一维非恒定流数学模型的河道断面自动插补优化模型进行并行化处理，将步骤S2中的组合数

分配给不同的线程同时计算。

所述步骤3中的以形定位计算模式是先优化插补断面形状，再确定插补断面位置：首先假定插补断面等距分布，然后按照目标函数对各插补断面的形状进行优化计算，得到断面底宽、河底高程和边坡系数，最后对各插补断面的断面间距进行寻优计算，得到各插补断面最优的间距；具体包括以下步骤：

步骤A1：对参数进行初始化，确定各参数的取值范围和各参数的个数k_i，i＝1,2,…,6；

步骤A2：确定河道糙率取值；

步骤A3：确定插补断面个数k₁，等间距布置插补断面，并对各插补断面按照从下游到上游的顺序进行编号；

步骤A4：令断面形状寻优迭代次数L＝1，插补断面j＝1；

步骤A5：将插补断面j的

个形状参数组合方案代入模型，结合河道断面自动插补优化模型进行并行寻优计算，再转到步骤A6；

步骤A6：固定已优化的插补断面形状参数，再转到步骤A7；

步骤A7：判断插补断面j的寻优编号是否等于插补断面个数k₁，若否，令j＝j+1，再转到步骤A5；若是，则得到各插补断面形状优化结果和本次形状迭代目标函数值F_L，再转到步骤A8；

步骤A8：判断断面形状寻优迭代次数L是否大于1，若否，令L＝L+1，再转到步骤A5；若是，则转到步骤A9；

步骤A9：判断F_L和上次形状迭代目标函数值F_L-1两次迭代结果是否满足收敛条件，若否，令L＝L+1，再转到步骤A5；若是，则得到当前河道糙率和插补断面个数下的各插补断面形状的最优方案；

步骤A10：令断面位置寻优迭代次数M＝1，插补断面j＝1；

步骤A11：固定各插补断面形状，对插补断面j的k₂个位置参数方案代入模型，结合河道断面自动插补优化模型进行并行寻优计算；再转到步骤A12；

步骤A12：固定已优化的插补断面位置参数；再转到步骤A13；

步骤A13：判断插补断面j是否等于k₂，若否，令j＝j+1，再转到步骤A11；若是，则得到各插补断面位置优化结果和本次位置迭代目标函数值F_M，再转到步骤A14；

步骤A14：判断断面位置寻优迭代次数M是否大于1，若否，令M＝M+1，再转到步骤A11；若是，则转到步骤A15；

步骤A15：判断F_M和上次位置迭代目标函数值F_M-1两次迭代结果是否满足收敛条件，若否，令M＝M+1，再转到步骤A11；若是，则得到当前河道糙率和插补断面个数下的各插补断面形状和位置最优方案，再转到步骤A16；

步骤A16：依次改变插补断面数量和河道糙率取值，重复上述步骤A1～步骤A15，得到不同插补断面数量和河道糙率取值方案下的各插补断面形状和位置最优方案。

所述步骤3中的以位定形计算模式是先确定插补断面位置，再优化插补断面形状：首先设定数量足够多且形状相同的插补断面，并将断面形状设定为矩形，然后仅改变各插补断面的河底高程进行寻优计算，对河底高程相近的断面进行合并，以此选出具有代表性的位置，最后再对各插补断面形状进行寻优，得到各插补断面的最优形状参数；具体包括以下步骤：

步骤B1：对参数进行初始化，确定各参数的取值范围和各参数的个数k_i，i＝1,2,…,6，的取值；

步骤B2：确定河道糙率取值；

步骤B3：确定矩形插补断面个数k_矩，等间距布置插补断面，并对各插补断面按照从下游到上游的顺序进行编号；

步骤B4：令断面位置寻优迭代次数M＝1，插补断面j＝1；

步骤B5：将插补断面j的河底高程计算方案代入模型，结合河道断面自动插补优化模型进行并行寻优计算，再转到步骤B6；

步骤B6：固定已优化的插补断面河底高程，再转到步骤B7；

步骤B7：判断插补断面j的寻优编号是否等于矩形插补断面个数k_矩，若否，令j＝j+1，再转到步骤B5；若是，则得到各插补断面的河底高程优化结果和本次位置迭代目标函数值F_M，再转到步骤B8；

步骤B8：判断断面位置寻优迭代次数M是否大于1，若否，令M＝M+1，再转到步骤B5；若是，则转到步骤B9；

步骤B9：判断F_M和上次位置迭代目标函数值F_M-1两次迭代结果是否满足收敛条件，若否，令M＝M+1，再转到步骤A5；若是，则得到当前河道糙率和插补断面个数下的各插补断面的河底高程变化情况，合并河底高程相近的断面，选出具有代表性的位置；

步骤B10：令断面形状寻优迭代次数L＝1，插补断面j＝1；

步骤B11：确定插补断面个数k₁，将插补断面j的形状参数方案代入模型，结合河道断面自动插补优化模型进行并行寻优计算，再转到步骤B12；

步骤B12：固定已优化的插补断面形状参数，再转到步骤B13；

步骤B13：判断插补断面j是否等于k₁，若否，令j＝j+1，再转到步骤B11；若是，则得到各插补断面位置优化结果和本次形状迭代目标函数值F_L，再转到步骤B14；

步骤B14：判断断面形状寻优迭代次数L是否大于1，若否，令L＝L+1，再转到步骤B11；若是，则转到步骤B15；

步骤B15：判断F_L和上次形状迭代目标函数值F_L-1两次迭代结果是否满足收敛条件，若否，令L＝L+1，再转到步骤B11；若是，则得到当前河道糙率和插补断面个数下的各插补断面形状和位置最优方案；

步骤B16：依次改变插补断面数量和河道糙率取值，重复上述步骤B1～步骤B15，得到不同插补断面数量和河道糙率取值方案下的各插补断面形状和位置最优方案。

本发明的有益效果在于：

1、本发明实现了对缺少地形资料的山区河道插补断面的个数、位置、形状和河道糙率等参数的自动优化；

2、本发明基于逐次逼近法的并行求解方法，减少了计算量，降低了计算难度；

3、本发明所建立的模型，对于断面形状规则的概化断面(如矩形、梯形、三角形等)，具有一定的通用性，仅需修改形状参数取值，即可得到最优的插补方案；与手动调参相比，不仅更省时省力，还可避免人工调试的主观任意性和不确定性，提高模拟精度；

4、所建立的河道插补断面自动优化模型具备较好的并行计算条件，通过OpenMP编译处理技术对计算程序进行并行化处理，可大幅提高计算效率，节省计算时间。

附图说明

图1为本发明基于逐次逼近法的河道断面自动插补方法的流程图；

图2为“以形定位”模式的计算流程图；

图3为“以位定形”模式的计算流程图；

图4为插补断面(D0₁～D0₅)位置示意图；

图5为D0断面的流量过程曲线图；

图6为D12断面的水位过程曲线图；

图7为断面形状变量示意图；

图8为“以形定位”模式D0断面水位实测值和计算值对比图；

图9为“以形定位”模式D2断面流量实测值和计算值对比图；

图10为“以形定位”模式D4断面水位实测值和计算值对比图；

图11为117个插补断面的河底高程参数变化曲线图；

图12为“以位定形”模式D0断面水位实测值和计算值对比图；

图13为“以位定形”模式D2断面流量实测值和计算值对比图；

图14为“以位定形”模式D4断面水位实测值和计算值对比图。

具体实施方式

本发明提出一种基于逐次逼近法的河道断面自动插补方法，下面结合附图和具体实施例对本发明做进一步说明。

图1为本发明基于逐次逼近法的河道断面自动插补方法的流程图；具体包括：

1模型建立

将观测断面水位和流量等水文要素的纳什系数和最大作为目标函数，水位、流量误差和断面间距作为约束条件，建立基于一维非恒定流数学模型的河道插补断面自动优化模型。

1.1一维非恒定流数学模型

采用Saint-Venant方程组描述河道一维水流运动，其水流连续方程和运动方程分别为

式中：t为时间；x为流程；Q为流量；Z为水位；A为过水断面面积；B为河宽；R为水力半径；n为糙率；V为断面平均流速；q_l和u_l为河段单位长度侧向入流量和侧向入流在x方向的分量；α₁为动量修正系数，α₁＝(∫_Au²dA)/(Q²/A)；g为重力加速度。

采用线性化的Preissmann四点隐格式对方程(1)和(2)进行离散，并采用追赶法求解离散后的方程组。

1.2河道插补断面自动优化模型

1.2.1目标函数

通过一维非恒定流数学模型计算得到观测断面的水文要素(水位和流量)的计算值，采用纳什系数反映观测断面的水文要素实测值与计算值的吻合程度，以观测断面的水文要素纳什系数和最大为目标函数，可表示为：

式中，i为观测断面的序号；N为观测断面总数；NSE(i)为观测断面i的纳什系数值，具体可表示为：

式中，T为计算周期，y_t(i)为i断面的水位(流量)实测值；

为i断面水位(流量)实测序列的平均值；

为i断面水位(流量)计算值。

1.2.2约束条件

(1)最大水位(流量)误差约束。河道插补断面自动优化模型的目标函数是要获得使在计算周期内的水位(流量)纳什系数和最大的河道断面插补个数、位置、形状和河道糙率，是对河道断面插补个数、位置、形状和河道糙率取值的综合考虑，不排除个别时段某个断面水位、流量误差过大的情况，因此应给出最大水位、流量误差约束，即

式中：abs为绝对值函数，ε_i为第i个断面的最大水位(流量)误差。

(2)断面间距约束。插补断面间距不得超过给定的最大断面间距，并不得小于给定的最小断面间距，即

L_min≤L(k)≤L_max (6)

式中，L(k)为插补断面k与前一断面的断面间距，模型中k＝2,3,4,5；L_max为给定的最大断面间距；L_min为给定的最小断面间距，一般取为0。

2模型的主要参数

由于插补断面的个数、位置、形状和河道糙率均会在不同程度上影响河道的沿程水位和流量变化，所以河道插补断面自动优化模型的主要参数为河道糙率、插补断面个数、位置和形状。又由于插补断面一般为概化断面(如矩形、梯形、三角形等)，所以断面形状参数主要包括底宽、河底高程和边坡系数等。根据河道长度、宽度和比降等资料，可大致确定主要参数的取值范围。

3模型求解方法

本实施例所建模型不仅涉及的主要参数较多，且每次寻优都要求解一维非恒定流数学模型，而一维非恒定流数学模型的圣维南方程组是一个非线性方程组，需通过不断叠代进行求解。又由于要考虑插补断面形状和位置之间的相互作用，所以河道插补断面自动优化模型构成一个非常复杂的高维、多变量、非线性优化问题，计算难度较大。

因此，为减少计算量，降低计算难度，本实施例以河道断面插补个数、位置、形状(包括断面底宽、河底高程、边坡系数等)和河道糙率作为状态变量，以观测断面的实测水位过程以及流量过程作为决策变量，提出了基于逐次逼近法的并行求解方法，计算时首先对主要参数进行离散，设插补断面的个数为k₁，位置变量个数为k₂，形状变量包括断面底宽、河底高程、边坡系数等变量个数分别为k₃、k₄、k₅，河道糙率变量个数为k₆，则共有

(即k₁×k₂×k₃×k₄×k₅×k₆)个变量组合；然后固定插补断面个数、位置和河道糙率的取值，仅将各插补断面形状变量组合(变量组合数为

)代入模型进行逐次寻优，这样便将复杂的多断面寻优问题分解为多个简单的单断面寻优问题，大大降低河道插补断面优化计算的求解难度。但变量组合数即使减少为

计算量仍然很大，故本发明结合针对单主机多核/多CPU并行计算而设计的OpenMP编译处理技术对模型进行并行化处理，将变量组合数

分配给不同的线程同时进行计算，进一步提高计算效率，节省计算时间。

3.1“以形定位”模式

图2为“以形定位”模式的计算流程图；“以形定位”模式是先优化插补断面形状，再确定插补断面位置。首先假定插补断面等距分布；然后按照目标函数对各插补断面的形状进行优化计算，得到插补断面的形状参数(底宽、河底高程、边坡系数)；最后在一定范围内，对各插补断面的断面间距进行寻优计算，得到各插补断面最优的间距。

“以形定位”模式的求解主要分为两个阶段，第一个阶段主要计算各插补断面的最优形状，第二个阶段主要计算各插补断面的最优位置，具体计算步骤如下：

3.1.1形状优化阶段

步骤A1：对参数进行初始化，确定各参数的取值范围和各参数的个数k_i，i＝1,2,…,6；

步骤A2：确定河道糙率取值；

步骤A3：确定插补断面个数k₁，等间距布置插补断面，并对各插补断面按照从下游到上游的顺序进行编号；

步骤A4：令断面形状寻优迭代次数L＝1，插补断面j＝1；

步骤A5：将插补断面j的

个形状参数组合方案代入模型，结合河道断面自动插补优化模型进行并行寻优计算，再转到步骤A6；

步骤A6：固定已优化的插补断面形状参数，再转到步骤A7；

步骤A7：判断插补断面j的寻优编号是否等于插补断面个数k₁，若否，令j＝j+1，再转到步骤A5；若是，则得到各插补断面形状优化结果和本次形状迭代目标函数值F_L，再转到步骤A8；

步骤A8：判断断面形状寻优迭代次数L是否大于1，若否，令L＝L+1，再转到步骤A5；若是，则转到步骤A9；

步骤A9：判断F_L和上次形状迭代目标函数值F_L-1两次迭代结果是否满足收敛条件，若否，令L＝L+1，再转到步骤A5；若是，则得到当前河道糙率和插补断面个数下的各插补断面形状的最优方案；

3.1.2位置优化阶段

步骤A10：令断面位置寻优迭代次数M＝1，插补断面j＝1；

步骤A11：固定各插补断面形状，对插补断面j的k₂个位置参数方案代入模型，结合河道断面自动插补优化模型进行并行寻优计算；再转到步骤A12；

步骤A12：固定已优化的插补断面位置参数；再转到步骤A13；

步骤A13：判断插补断面j是否等于k₂，若否，令j＝j+1，再转到步骤A11；若是，则得到各插补断面位置优化结果和本次位置迭代目标函数值F_M，再转到步骤A14；

步骤A14：判断断面位置寻优迭代次数M是否大于1，若否，令M＝M+1，再转到步骤A11；若是，则转到步骤A15；

步骤A15：判断F_M和上次位置迭代目标函数值F_M-1两次迭代结果是否满足收敛条件，若否，令M＝M+1，再转到步骤A11；若是，则得到当前河道糙率和插补断面个数下的各插补断面形状和位置最优方案，再转到步骤A16；

步骤A16：依次改变插补断面数量和河道糙率取值，重复上述步骤A1～步骤A15，得到不同插补断面数量和河道糙率取值方案下的各插补断面形状和位置最优方案。

3.2“以位定形”模式

“以位定形”模式是先确定插补断面位置，再优化插补断面形状。首先设定数量足够多且形状相同的插补断面；然后仅改变各插补断面的河底高程进行寻优计算，根据计算结果确定具有代表性的插补断面位置；最后按照目标函数再对各插补断面形状进行寻优，得到各插补断面的最优形状参数。

“以位定形”模式的求解主要分为两个阶段，第一个阶段主要计算插补断面的代表性位置，在河道中插补足够多的断面，并假定每个断面的形状相同，一般将断面形状设定为矩形。然后仅改变各断面的河底高程进行寻优计算，对河底高程相近的断面进行合并，以此选出具有代表性的位置；第二个阶段主要计算各插补断面的最优形状，具体计算步骤如下：

3.2.1位置优化阶段

步骤B1：对参数进行初始化，确定各参数的取值范围和各参数的个数k_i，i＝1,2,…,6，的取值；

步骤B2：确定河道糙率取值；

步骤B3：确定矩形插补断面个数k_矩，等间距布置插补断面，并对各插补断面按照从下游到上游的顺序进行编号；

步骤B4：令断面位置寻优迭代次数M＝1，插补断面j＝1；

步骤B5：将插补断面j的河底高程计算方案代入模型，结合河道断面自动插补优化模型进行并行寻优计算，再转到步骤B6；

步骤B6：固定已优化的插补断面河底高程，再转到步骤B7；

步骤B7：判断插补断面j的寻优编号是否等于矩形插补断面个数k_矩，若否，令j＝j+1，再转到步骤B5；若是，则得到各插补断面的河底高程优化结果和本次位置迭代目标函数值F_M，再转到步骤B8；

步骤B8：判断断面位置寻优迭代次数M是否大于1，若否，令M＝M+1，再转到步骤B5；若是，则转到步骤B9；

步骤B9：判断F_M和上次位置迭代目标函数值F_M-1两次迭代结果是否满足收敛条件，若否，令M＝M+1，再转到步骤A5；若是，则得到当前河道糙率和插补断面个数下的各插补断面的河底高程变化情况，合并河底高程相近的断面，选出具有代表性的位置；

3.2.2形状优化阶段

计算步骤同“以形定位”模式的步骤1～9。

“以位定形”模式的计算流程图如图3所示。

将所建模型和方法应用于全长为25.5km的河道中，其间有D0～D12共13个断面，其中上游D0断面～D1断面(含D0断面)缺少河道大断面资料，需要插补断面并进行优化，插补断面(D01～D05)位置如图4所示。

D0断面为河道进口断面，D12断面为河道出口断面，D0断面～D12断面距D0断面里程如表1所示，D0断面的流量过程和D12断面的水位过程如图5和图6所示。

表1断面里程

4河道插补断面自动优化计算方案如下：

(1)断面个数

断面个数取值分为为2～5，共4种方案。

(2)河道糙率

假定河道所有断面的糙率相同，取值范围为0.034～0.043，共10种方案。

(3)断面形状

假定插补断面的形状为等腰梯形，通过河底高程Z_p、底宽b和边坡系数m三个变量可计算插补断面相应的过水断面面积A、水面宽B、水面宽随水位的变化率

再结合一维非恒定流数学模型，即可计算得到河道沿程水位和流量过程。具体关系式如式(7)～(10)所示和断面形状变量示意图如图7所示：

A＝(b+m(Z-Z_p))(Z-Z_p) (7)

B＝b+m(Z-Z_p) (8)

m＝cotθ (10)

给定插补断面形状河底高程Z_p、底宽b和边坡系数m的取值范围，便可确定断面形状优化的计算方案。本发明中坡角θ取值范围为2°～85°，将其离散为20个值，即可得到对应的20个边坡系数m的离散值；底宽取值范围b为50m～500m，将其离散为10个值；河底高程Z_p取值范围为329m～348m，将其离散为20个值，则共计4000种插补断面形状优化计算方案。具体参数取值如表2所示。

表2插补断面底宽b、河底高程Z_p和边坡系数m的离散值

5“以形定位”模式的计算结果

采用“以形定位”模式得到插补2～5个断面不同糙率下的优化计算结果如表3所示。

表3“以形定位”模式优化计算结果

由表3可知，“以形定位”模式下，插补断面为2～5时目标函数值均随糙率的增加而增大，在糙率为0.043时，目标函数值最大，分别为2.778、2.821、2.860和2.861。且在糙率为0.043时，插补断面个数为2～5时的D0断面水位、D2断面流量和D4断面水位纳什系数均大于0.8，表明通过在“以形定位”模式下可以得到精度较高的结果。其中插补断面个数为5时的目标函数值最优，D0断面水位、D2断面流量和D4断面水位的纳什系数分别为0.993、0.958和0.910，具体计算结果如表4所示。

表4“以形定位”模式最优计算结果

插补5个断面的最优计算结果中D0断面水位、D2断面流量、D4断面水位计算值和实测值分别如图8～图10所示。从图中可以看出计算值和实测值模拟效果良好，其中D0断面水位、D2断面流量、D4断面水位的平均相对误差分别为0.03％、6.74％、0.05％。

6“以位定形”模式的计算结果

在D0断面～D1断面间等间距插补117个(含D0断面)矩形断面，对117个插补断面的河底高程进行优化计算，得到各插补断面的河底高程参数变化情况如图11所示。

由图11可以看出，在距D0断面里程0m～700m、2250m～3050m、3450m～3950m、4050m～5200m和5350m～5850m的区域内，插补断面的河底高程有明显的上下调整趋势，表明此区域对目标函数值的影响较大。因此，本实施例根据插补断面的河底高程调整趋势选定了6个具有代表性的断面位置，分别为距D0断面里程0m、700m、3050m、3950m、5200m和5850m的位置，其中前5个为5个插补断面对应的位置，第6个为D1断面所在的位置。

将5个插补断面的代表性位置进行有序组合，即当插补断面个数为2时，插补断面位置组合数为4；插补断面个数为3时，插补断面位置组合数为6；插补断面个数为4时，插补断面位置组合数为4。

根据“以形定位”模式的糙率最优结果，将糙率取值为0.043，对不同插补断面个数下的断面位置组合依次进行寻优，得到不同插补断面个数下的最优断面位置组合，结果如表5所示。

表5不同插补断面个数下的最优断面位置

由表5可知，除了进口断面D0₅(D0)断面位置不变外，插补断面个数从2个增加至5个，断面插补顺序依次为D0₁、D0₂、D0₄和D0₃。

给定插补2～5个断面的断面位置后，不同糙率下的优化计算结果，如表6所示。

表6“以位定形”模式优化计算结果

由表6可知，“以位定形”模式下，插补断面为2～5时目标函数值均随糙率的增加而增加，在糙率为0.043时，目标函数值最大，分别为2.827、2.849、2.860和2.863。且在糙率为0.043时，插补断面个数为2～5时的D0断面水位、D2断面流量和D4断面水位纳什系数均大于0.8，表明在“以位定形”模式下可以得到精度较高的结果。其中插补断面个数为5时的目标函数值最优，D0断面水位、D2断面流量和D4断面水位的纳什系数分别为0.996、0.959和0.908，具体计算结果如表7所示。

表7“以位定形”模式最优计算结果

“以位定形”最优计算结果中，D0断面水位、D2断面流量、D4断面水位的计算值和实测值分别图12～图14所示。从图中可以看出计算值和实测值模拟效果良好，其中D0断面水位、D2断面流量、D4断面水位的平均相对误差分别为0.02％、6.68％、0.04％。

7计算耗时

本模型采用Fortran编程语言在Microsoft Visual Studio 2019开发平台上实现，利用OpenMP编译处理技术中对程序进行并行化处理。计算程序运行环境为：①四十核处理器Intel(R)Xeon(R)Gold 5218R二颗，主频2.1GHz；②内存128GB；③操作系统Windows10专业版；④配置管理器为x64，Release模式。“以形定位”模式每次寻优时的变量组合数为4000，在并行计算时将并行计算线程数设为75，“以形定位”模式下并行与串行计算时间如表8所示。

表8“以形定位”模式下不同插补断面个数下并行与串行计算耗时(线程数75)

从表8中可知，“以形定位”模式下当河道插补断面个数从2个增到5个时，计算量随之增加，串行计算耗时从3771.56s增加到42870.83s，增加了11.37倍。采用并行计算时，其计算耗时相比于串行计算将明显减少，如当插补断面个数为5时，串行计算耗时42870.83s(约12h)，并行计算耗时587.14s(约10min)，并行计算耗时仅为串行耗时的1.37％，加速比达73.02，效率为0.974。同时随着插补断面个数的增加，并行计算加速比和效率也在逐渐增大，计算性能显著增强。

本发明采用基于逐次逼近法的并行求解方法对河道插补断面自动优化模型进行求解，计算结果表明：(1)“以形定位”模式和“以位定形”模式的计算值和实测值模拟效果均比较良好，“以形定位”模式下D0断面水位、D2断面流量和D4断面水位的纳什系数分别为0.993、0.958和0.910，“以位定形”模式下纳什系数分别为0.996、0.959和0.908，模拟精度较高；(2)“以形定位”模式和“以位定形”模式下的最优目标函数值分别为2.861和2.863，“以位定形”模式模拟精度要高于“以形定位”模式，先优化断面位置后优化断面形状的“以位定形”计算模式能更好地反映河道的变化特征；(3)本实施例所建立的河道插补断面自动优化模型具备较好的并行计算条件，通过OpenMP编译处理技术对计算程序进行并行化处理，可大幅提高计算效率，节省计算时间，在“以形定位”模式下并行计算耗时不超过串行耗时的2％，最大加速比为73.02，效率为0.974。

此实施例仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (5)

1.一种基于逐次逼近法的河道断面自动插补方法，其特征在于，包括以下步骤：

步骤1：以观测断面水位与流量水文要素的纳什系数之和最大为目标函数，以水位、流量误差和断面间距为约束条件，建立基于一维非恒定流数学模型的河道断面自动插补优化模型；

步骤2：根据河道的长度、宽度和比降，确定河道断面自动插补优化模型的参数取值范围；所述参数包括河道糙率、插补断面个数、插补断面位置、插补断面形状参数，其中插补断面形状参数包括断面底宽、河底高程和边坡系数；

步骤3：根据基于逐次逼近法的并行求解方法，采用以形定位计算模式或以位定形计算模式对河道断面自动插补优化模型进行求解。

2.根据权利要求1所述基于逐次逼近法的河道断面自动插补方法，其特征在于，所述步骤1中河道断面自动插补优化模型为：

目标函数：

式中，i为观测断面的序号；N为观测断面总数；NSE(i)为观测断面i的纳什系数值，具体为：

式中，T为计算周期，y_t(i)为i断面的水位或流量实测值；

为i断面水位或流量实测序列的平均值；

为i断面水位或流量计算值；

约束条件：

式中：abs为绝对值函数，ε_i为第i个断面的最大水位或流量误差；

L_min≤L(k)≤L_max (4)

式中，L(k)为插补断面k与前一断面的断面间距；L_max为给定的最大断面间距；L_min为给定的最小断面间距，取为0。

3.根据权利要求1所述基于逐次逼近法的河道断面自动插补方法，其特征在于，所述步骤3中的基于逐次逼近法的并行求解方法具体包括以下步骤：

步骤S1：对参数进行离散，设插补断面个数为k₁，插补断面位置的个数为k₂，插补断面形状中断面底宽的个数为k₃、河底高程的个数为k₄、边坡系数的个数为k₅，河道糙率的个数为k₆，则共有

个参数组合；

步骤S2：固定插补断面个数、插补断面位置和河道糙率的取值，仅将各插补断面形状中的

个参数组合代入基于一维非恒定流数学模型的河道断面自动插补优化模型进行逐次寻优；

步骤S3：结合针对单主机多CPU并行计算而设计的OpenMP编译处理技术对基于一维非恒定流数学模型的河道断面自动插补优化模型进行并行化处理，将步骤S2中的组合数

分配给不同的线程同时计算。

4.根据权利要求1所述基于逐次逼近法的河道断面自动插补方法，其特征在于，所述步骤3中的以形定位计算模式是先优化插补断面形状，再确定插补断面位置：首先假定插补断面等距分布，然后按照目标函数对各插补断面的形状进行优化计算，得到断面底宽、河底高程和边坡系数，最后对各插补断面的断面间距进行寻优计算，得到各插补断面最优的间距；具体包括以下步骤：

步骤A1：对参数进行初始化，确定各参数的取值范围和各参数的个数k_i，i＝1,2,…,6；

步骤A2：确定河道糙率取值；

步骤A3：确定插补断面个数k₁，等间距布置插补断面，并对各插补断面按照从下游到上游的顺序进行编号；

步骤A4：令断面形状寻优迭代次数L＝1，插补断面j＝1；

步骤A5：将插补断面j的

个形状参数组合方案代入模型，结合河道断面自动插补优化模型进行并行寻优计算，再转到步骤A6；

步骤A6：固定已优化的插补断面形状参数，再转到步骤A7；

步骤A7：判断插补断面j的寻优编号是否等于插补断面个数k₁，若否，令j＝j+1，再转到步骤A5；若是，则得到各插补断面形状优化结果和本次形状迭代目标函数值F_L，再转到步骤A8；

步骤A8：判断断面形状寻优迭代次数L是否大于1，若否，令L＝L+1，再转到步骤A5；若是，则转到步骤A9；

步骤A9：判断F_L和上次形状迭代目标函数值F_L-1两次迭代结果是否满足收敛条件，若否，令L＝L+1，再转到步骤A5；若是，则得到当前河道糙率和插补断面个数下的各插补断面形状的最优方案；

步骤A10：令断面位置寻优迭代次数M＝1，插补断面j＝1；

步骤A11：固定各插补断面形状，对插补断面j的k₂个位置参数方案代入模型，结合河道断面自动插补优化模型进行并行寻优计算；再转到步骤A12；

步骤A12：固定已优化的插补断面位置参数；再转到步骤A13；

步骤A13：判断插补断面j是否等于k₂，若否，令j＝j+1，再转到步骤A11；若是，则得到各插补断面位置优化结果和本次位置迭代目标函数值F_M，再转到步骤A14；

步骤A14：判断断面位置寻优迭代次数M是否大于1，若否，令M＝M+1，再转到步骤A11；若是，则转到步骤A15；

步骤A15：判断F_M和上次位置迭代目标函数值F_M-1两次迭代结果是否满足收敛条件，若否，令M＝M+1，再转到步骤A11；若是，则得到当前河道糙率和插补断面个数下的各插补断面形状和位置最优方案，再转到步骤A16；

步骤A16：依次改变插补断面数量和河道糙率取值，重复上述步骤A1～步骤A15，得到不同插补断面数量和河道糙率取值方案下的各插补断面形状和位置最优方案。

5.根据权利要求1所述基于逐次逼近法的河道断面自动插补方法，其特征在于，所述步骤3中的以位定形计算模式是先确定插补断面位置，再优化插补断面形状：首先设定数量足够多且形状相同的插补断面，并将断面形状设定为矩形，然后仅改变各插补断面的河底高程进行寻优计算，对河底高程相近的断面进行合并，以此选出具有代表性的位置，最后再对各插补断面形状进行寻优，得到各插补断面的最优形状参数；具体包括以下步骤：

步骤B1：对参数进行初始化，确定各参数的取值范围和各参数的个数k_i，i＝1,2,…,6，的取值；

步骤B2：确定河道糙率取值；

步骤B3：确定矩形插补断面个数k_矩，等间距布置插补断面，并对各插补断面按照从下游到上游的顺序进行编号；

步骤B4：令断面位置寻优迭代次数M＝1，插补断面j＝1；

步骤B5：将插补断面j的河底高程计算方案代入模型，结合河道断面自动插补优化模型进行并行寻优计算，再转到步骤B6；

步骤B6：固定已优化的插补断面河底高程，再转到步骤B7；

步骤B7：判断插补断面j的寻优编号是否等于矩形插补断面个数k_矩，若否，令j＝j+1，再转到步骤B5；若是，则得到各插补断面的河底高程优化结果和本次位置迭代目标函数值F_M，再转到步骤B8；

步骤B8：判断断面位置寻优迭代次数M是否大于1，若否，令M＝M+1，再转到步骤B5；若是，则转到步骤B9；

步骤B9：判断F_M和上次位置迭代目标函数值F_M-1两次迭代结果是否满足收敛条件，若否，令M＝M+1，再转到步骤A5；若是，则得到当前河道糙率和插补断面个数下的各插补断面的河底高程变化情况，合并河底高程相近的断面，选出具有代表性的位置；

步骤B10：令断面形状寻优迭代次数L＝1，插补断面j＝1；

步骤B11：确定插补断面个数k₁，将插补断面j的形状参数方案代入模型，结合河道断面自动插补优化模型进行并行寻优计算，再转到步骤B12；

步骤B12：固定已优化的插补断面形状参数，再转到步骤B13；

步骤B13：判断插补断面j是否等于k₁，若否，令j＝j+1，再转到步骤B11；若是，则得到各插补断面位置优化结果和本次形状迭代目标函数值F_L，再转到步骤B14；

步骤B14：判断断面形状寻优迭代次数L是否大于1，若否，令L＝L+1，再转到步骤B11；若是，则转到步骤B15；

步骤B15：判断F_L和上次形状迭代目标函数值F_L-1两次迭代结果是否满足收敛条件，若否，令L＝L+1，再转到步骤B11；若是，则得到当前河道糙率和插补断面个数下的各插补断面形状和位置最优方案；

步骤B16：依次改变插补断面数量和河道糙率取值，重复上述步骤B1～步骤B15，得到不同插补断面数量和河道糙率取值方案下的各插补断面形状和位置最优方案。

Images