CN108256177A - 一种河流水沙模型的参数优化方法和系统 - Google Patents

Abstract

本发明实施例提供了一种河流水沙模型的参数优化方法和系统，该方法和系统具体为生成模型计算网格；采集河道的待预测区域的初始场和边界信息，形成待预测区域的信息文件；根据待预测区域的信息文件，离散求解水沙模型；采用拉丁超立方采样方法LHS对水沙模型的多个参数在可行域内进行分层采样，得到多个参数组；基于模型参数组合，求解河流水沙模型，获得每个参数组对应的似然度值；设定似然度阈值，将低于预设阈值的似然度值所对应的参数组予以删除；将剩余的参数组所对应的似然度值进行排序，确定不确定性区间；利用GLUE累积对似然度值进行区域敏感度分析。能够解决了该模型因参数不正确表达而导致其预测精度较差的问题。

Description

一种河流水沙模型的参数优化方法和系统

技术领域

本发明涉及水利工程技术领域，特别是涉及一种河流水沙模型的参数优化方法和系统。

背景技术

冲积河流过程水沙模型作为模拟预测水沙运动变化和有效分析水沙运动规律的重要手段之一，可以有效服务于河道整治、防洪调度、生态保护、水利工程建设等生产实践活动，一直是国内外学者重点关注和研究的对象。

随着全球气候变化、流域下垫面变化以及大坝修建等自然变异和人类活动改变了天然的来水来沙情势和水沙输移边界条件，由此可能会引发溃坝等极端事件的发生，这进一步会导致活跃的水沙运动和快速的河床变形，从而给科学认知冲积河流过程水沙相互作用机理带来了困难，给模拟预测水沙运动变化带来了更大的挑战。

具体来说，冲积河流水沙模型的发展水平日益成熟，但受模型结构、输入条件及模型参数误差等不确定性因素的影响，导致模型无法准确预测真实水沙的运动状态，这些不确定性因素的存在成为制约河流水沙模型精度提高的主要瓶颈。特别对于河流水沙模型来说，即使有精确的初始数据和输入数据，模型参数的不准确表达也会导致误差的增长，影响河流水沙模型准确预测真实水沙运动状态的能力。

发明内容

有鉴于此，本发明提供了一种河流水沙模型的参数优化方法和系统，以解决河流水沙模型因模型参数不正确表达而导致其预测精度较差的问题。

为了解决上述问题，本发明公开了一种河流水沙模型的参数优化方法，包括如下步骤：

基于卫星遥感影像数据和河道的地形高程数据生成计算网格；

采集所述河道的待预测区域的初始场和边界条件信息，并进行排列得到所述待预测区域的信息文件；

根据所述待预测区域的信息文件求解水沙模型；

采用拉丁超立方体采样对所述水沙模型的多个模型参数在可行域内进行分层采样，得到多个参数组；

利用所述模型参数组合，求解河流水沙模型得到每个所述参数组所对应的似然度值；

将所述多个参数组中，低于设定阈值的似然度值所对应的参数组予以删除；

将所述多个参数组中剩余的参数组所对应的似然度值进行排序，确定模型预测的不确定性区间；

利用GLUE累积似然度对模型参数进行区域敏感度分析。

可选的，所述基于卫星遥感影像数据和河道的地形高程数据生成计算网格，包括：

基于卫星遥感影像数据，提取水体边界，生成河道平面三角形网格；

利用河道地形高程数据对所述河道平面三角形网格进行插值运算，得到所述河道平面三角形网格的节点高程。

可选的，所述区域初始场和边界条件信息包括所述河道高程、河宽、上游边界的流量和含沙量、下游边界的水位、多个断面的初始流量场、初始水位场、初始含沙量场、悬沙级配、床沙级配中的部分或全部。

可选的，所述水沙模型为水动力学控制方程、悬移质不平衡输移方程、推移质不平衡输移方程、床面变形方程和挟沙力计算公式中的部分或全部。

可选的，所述模型参数包括水动力参数和泥沙参数，其中：

所述水动力参数包括糙率系数和/或紊动粘滞系数；

所述泥沙参数包括紊动扩散系数、悬移质泥沙恢复饱和系数、推移质泥沙恢复饱和系数、悬移质挟沙力经验参数和推移质挟沙力经验参数中的部分或全部。

相应的，为了保证上述方法的实施，本发明还提供了一种河流水沙模型的参数优化系统，包括：

网格生成模块，用于基于卫星遥感影像数据和河道的地形高程数据生成计算网格；

信息采集模块，用于采集所述河道的待预测区域的初始场和边界条件信息，并进行排列得到所述待预测区域的信息文件；

模型计算模块，用于根据所述待预测区域的信息文件求解水沙模型；

分层采样模块，用于采用拉丁超立方体采样对所述水沙模型的多个模型参数在可行域内进行分层采样，得到多个参数组；

似然度计算模块，用于利用所述模型参数组合，求解河流水沙模型得到每个所述参数组所对应的似然度值；

参数组筛选模块，用于设定似然度阈值，将所述多个参数组中，低于该阈值的似然度值所对应的参数组予以删除；

不确定性分析模块，用于将所述多个参数组中剩余的参数组所对应的似然度值进行排序，确定模型预测的不确定性区间；

敏感度分析模块，用于利用GLUE累积所述似然度值对模型参数进行区域敏感度分析。

可选的，所述网格生成模块包括：

平面网格生成单元，用于基于卫星遥感影像数据，提取水体边界，生成河道平面三角形网格；

插值计算单元，用于利用河道地形高程数据对所述河道平面三角形网格进行插值运算，得到所述河道平面三角形网格的节点高程。

可选的，所述区域初始场和边界条件信息包括所述河道高程、河宽、上游边界的流量和含沙量、下游边界的水位、多个断面的初始流量场、初始水位场、初始含沙量场、悬沙级配、床沙级配中的部分或全部。

可选的，所述水沙模型为水动力学控制方程、悬移质不平衡输移方程、推移质不平衡输移方程、床面变形方程和挟沙力计算公式中的部分或全部。

可选的，所述模型参数包括水动力参数和泥沙参数，其中：

所述水动力参数包括糙率系数和/或紊动粘滞系数；

所述泥沙参数包括紊动扩散系数、悬移质泥沙恢复饱和系数、推移质泥沙恢复饱和系数、悬移质挟沙力经验参数和推移质挟沙力经验参数中的部分或全部。

从上述技术方案可以看出，本发明实施例提供了一种河流水沙模型的参数优化系统，具体为基于卫星遥感影像数据和河道的地形高程数据生成河道计算网格；采集河道待预测区域的初始场和边界信息，形成待预测区域的信息文件；根据待预测区域的信息文件，离散求解水沙控制模型；采用拉丁超立方采样方法LHS对水沙模型的多个参数在可行域内进行分层采样，得到多个参数组；基于模型参数组合，求解河流水沙模型，获得每个参数组对应的似然度值；设定似然度阈值，将低于预设阈值的似然度值所对应的参数组予以删除；将剩余的参数组所对应的似然度值进行排序，确定不确定性区间；利用GLUE累积对似然度值进行区域敏感度分析。通过上述处理，可以给出河流水沙模型中模型参数的最优值。从而解决了该模型因参数不正确表达而导致其预测精度较差的问题。

另外，将LHS和GLUE相结合应用于河流水沙模型不确定性分析，更加精确反映了输入概率函数中值的分布，可以显著提高抽样效率。且利用实时观测数据的输入，能够实时更新模型参数，提高模型预测精度，同时实时给出不确定性区间，为风险管理提供了基础。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的一种河流水沙模型的参数优化方法的步骤流程图；

图2为本发明实施例提供的一种河流水沙模型的参数优化系统的结构框图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例一

图1为本发明实施例提供的一种河流水沙模型的参数优化方法的步骤流程图；

参照图1所示，本实施例提供的参数优化方法用于对河流水沙模型进行参数优化，具体步骤包括：

S101：基于卫星遥感影像数据和地形高程数据生成计算网格。

具体为，首先，基于卫星遥感影像数据，提取水体边界，在此基础上自动生成河道三角形网格；然后利用获取的河道的地形高程数据，插值生成河道三角形网格节点的河道高程。

S102：采集河道的待预测区域的区域信息。

该河道的区域信息包括前述得到的河道高程E，还包括河宽B、上游边界的流量Q_up和含沙量S_up、下游边界的水位Z_down，各个断面的初始流量场Q_初始、初始水位场Z_初始、初始含沙量场S_初始。悬沙级配、床沙级配等信息。并将上述信息按照采样时间的先后顺序分别进行排列，得到一张河道的待预测区域的信息文件。

S103：根据河道的待预测区域的信息文件离散求解水沙控制方程。

水沙控制方程包括水动力学控制方程、悬移质不平衡输移方程、推移质不平衡输移方程、床面变形方程、挟沙力计算公式等。

其中，动床条件下水动力学控制方程由二维浑水连续方程及x、y方向的运动方程组成，如下：

式中：t为时间；h为水深；u和v分别为x和y方向的水流流速；g为重力加速度；v_t为水流的紊动粘滞系数；Δρ＝ρ_s-ρ_f,ρ_s为泥沙密度，ρ_f为清水密度；ρ_m为浑水密度；S为总含沙量；ρ₀、ρ′分别为床沙的饱和湿密度及干密度；S_bx、S_by为床面底坡项，S_fx、S_fy为床面摩阻项，可分别表示为S_fx＝n²u(u²+v²)^1/2/h^4/3及S_fy＝n²v(u²+v²)^1/2/h^4/3，n为曼宁糙率系数；Z_b为河底高程。

二维悬移质泥沙的悬移质不平衡输移方程：

式中：ε_s为泥沙的紊动扩散系数；S_k、S_*k、ω_sk及α_sk分别为第k粒径组泥沙的含沙量、挟沙力、浑水沉速及恢复饱和系数。

二维推移质泥沙的推移质不平衡输移方程：

床面变形方程为由悬移质及推移质泥沙不平衡输移引起的床面变形方程，可分别表示为：

因此总的河床冲淤厚度可表示为：

式中，N为非均匀沙分组数；N_s为悬移质泥沙的分组数；ΔZ_sk、ΔZ_bk分别为悬移质及推移质泥沙在一个时间步长内引起的分组冲淤厚度；ΔZ_t为相应时段内总的床面冲淤厚度。因此时段末单元中心的河床底部高程可以表示为：(Z_b)^l+1＝(Z_b)^l+ΔZ_t，上标l指某一时间层。

挟沙力计算公式包括张瑞瑾公式、张红武公式、王兴奎公式、吴保生公式等。以吴保生公式为例：

式中，U为垂线平均流速，且U＝(u²+v²)^1/2；γ_s、γ_m分别为泥沙及浑水容重；ω_m为群体沉速。

推移质输沙力公式包括爱因斯坦公式、窦国仁公式等。以窦国仁提出的公式为例：

式中，q_b0k为第k粒径组推移质泥沙的单宽输沙力；K_b为一系数；U_ck为启动流速；C₀为无量纲的谢才系数，且C₀＝C/g^1/2，C＝h^1/6/n；ΔP_bk为床沙级配。故单位水体中推移质分组挟沙力可表示为q_b*k＝q_b0k/(hU)。

采用无结构三角网格下的有限体积法离散求解控制方程，同时采用单元中心的方式，即所有守恒变量值存储于单元中心，相邻两个单元的公共面为计算界面。无结构网格下的控制单元，如图2所示。

离散求解有限体积算法离散控制方程，

式中，U为守恒向量；E和G分别为x和y方向的对流通量；和分别为x和y方向的扩散通量；S为源项向量，包括床面坡降项、摩阻项及泥沙输移与床面变形引起的附加项。式中各项可以进一步表示为：

上式仅考虑悬移质输移方程。如忽略泥沙扩散项，推移质输移方程很容易考虑到上式中。将方程沿控制体A_i积分可得：

式中，假设单元平均值U存储在单元中心，则上式中的面积分可用线性积分进一步表示为：

其中，Γ为控制体A_i的边界；F_n(U)＝F·n＝En_x+Gn_y；n为边界Γ外法线方向的单位向量，n_x和n_y分别为n在x和y方向分量。式(13)中的线积分在三角单元中还可近似表示为：

式中，Δl_ij为第i单元第j条边的长度；F_ij、T_ij分别为通过界面Γ的对流通量与扩散通量。其中扩散通量一般可表示为式(13)表示为：

通过界面的水流及泥沙的对流通量F_ij采用Roe-MUSCL方法计算：

利用二维浅水方程的旋转不变性，可以将界面通量计算转换为求解局部的一维Riemann问题。通过单元A_i与A_i界面Γ的法向通量可以近似黎曼解求出，即：

F_ij＝F^*[(U_L)_ij，(U_R)_ij] (16)

式中，F^*为近似黎曼解；(U_L)_ij、(U_R)_ij分别为状态变量在界面两侧的值。采用不同的近似Riemann解及状态变量的插值方法，可以得到不同的计算格式。基于Roe近似Riemann解与MUSCL方法结合，即Roe-MUSCL方法的计算精度高。本发明实施例采用该方法计算水流对流通量，则上式进一步表示为：

式中，F(U_L)_ij、F(U_R)_ij分为为界面Γ两侧的法向通量；是基于Roe平均下的Jacobian矩阵。本发明实施例采用Van Leer提出的MUSCL方法重构界面两侧的状态变量值(U_R与U_L)，同时采用Roe和Baines提出的minmod函数对变量坡度进行限制。这类限制器可确保解的正性，且已广泛应用于具有TVD特性的计算格式中。因此采用Roe-MUSCL方法可使计算格式在空间方向具有二阶计算精度。

式(15)中的泥沙对流通量不采用近似Riemann求解，而是采用迎风格式计算泥沙对流通量，即可由泥沙浓度与通过界面的水流通量直接相乘求得。在图中，假设S_R和S_L分别界面两侧的浓度，P和Q分别为通过界面的单宽流量在x和y方向的分量。故通过界面的单宽流量F_f等于Pn_x+Qn_y，这样界面两侧的泥沙通量分别为F_f·S_L及F_f·S_R。采用迎风格式计算通过界面的泥沙对流通量，由下式表示：

水动力学控制方程中的源项处理。其中的源项主要包括床面底坡项、摩阻项、及由泥沙输移与床面冲淤引起的附加项。在三角形网格中，由于其3个顶点位于同一平面上，故床面底坡项较容易处理。根据已知3个顶点的坐标及高程，可求出该网格上的平面函数及相应底坡项。对于摩阻项，一般显式处理会导致计算不稳定，尤其在小水深时。故在本发明实施例中采用半隐式的计算格式处理该项。

通过直接求解单元内总含沙量变化梯度以及已知上一时刻的床面冲淤速率，可求出由泥沙输移与床面冲淤引起的附加项。

为获得高阶的时间离散格式，可以通过多种计算格式。本发明实施例采用Runge-Kutta方法获得二阶精度的时间离散格式，即常见的预估-校正格式：

U^l+1＝U^l+ΔtL(U^l+1/2)

式中，U^l+1/2＝(U^l+U^*)/2，U^*＝U_n+ΔtLU^l。由于该格式为显式，故时间步长受CFL条件限制。

在进行非恒定浅水问题模拟时，由于水位变化使得实际计算区域不断变化，为准确模拟这种动边界问题，通常需要引入界面干湿处理方法。在二维浅水控制方程的有限体积算法中，已有很多界面干湿处理的方法。如Sleigh、Zhao等提出了类似的界面干湿处理方法，他们在计算中将单元分为三类，即湿单元、干单元及半干单元。对半干单元的界面，仅考虑水流质量的输移，而不考虑动量传递。大部分模型在实际计算中，通常引入一最小水深来判断单元的干湿条件。在本发明实施例中，也引入最小水深来判断单元干湿，同时借鉴并改进了Falconer与Chen提出的规则计算网格中的干湿处理方法，使其能适用于无结构三角网格中的界面干湿处理。该方法的具体描述，详见相关文献。

实际的浅水流动问题一般包括两类边界：闭边界(陆地边界)及开边界。在陆地边界，一般采用滑移边界条件，即设边界处法向流速为零、切向流速不为零、且水深及含沙量等变量在边界上的法向梯度为零。对开边界，一般给定水位过程、或流量过程、或水位流量关系，而相应其他变量则采用Riemann不变问题求解。自由出流的开边界条件，一般设定各变量在边界上的法向梯度为零。

床面变形计算采用显格式离散床面变形方程式(6)，第l时间层到第l+1时间层内的ΔZ_sk及ΔZ_bk分别可由下式计算：

为模拟床面冲淤过程中的床沙粗化或细化现象，模型将床沙分为两大层，最上层的床沙层活动层及其该层以下的分层记忆层。床沙活动层厚度为H_m，相应级配为ΔP_bk。分层记忆层可根据实际情况共分n层，各层厚度及相应级配分别为ΔH_n、ΔP_nk。在计算中，当各粒径组泥沙发生淤积时，则记忆层数相应增加，且该层级配位t时刻的床沙活动层级配当各粒径组泥沙均发生冲刷时，根据冲刷量的大小，记忆层数响应减少若干层，且级配作用相应的调整。

S104：对水沙模型的多个参数进行分层采样，得到多个参数组。

具体为采用拉丁超立方采样LHS对水沙方程的模型参数中的水动力参数n、v_t，泥沙模型参数k、m、ε_s、α_sk、α_bk、k_b在可行域内进行分层采样，获取参数组。

利用拉丁超立方分层抽样方法分别将para₁,para₂,…,para_kk个参数空间范围划分为n个区间，在每个区间内均匀采样m次，每个参数共获取m×n个样本，将各个参数获取的m×n个样本进行随机组合，得到m×n个组合样本，每个样本代表一个参数组。

S105：利用水沙模型求解每个参数组对应的似然度值。

基于采样形成的水动力参数和泥沙参数构成的参数组，利用水沙模型进行求解，得到每个参数组对应的水位、流速、泥沙含量、河道高程等似然度值。

似然度函数是建立水位、流量、泥沙含量、河道高程等模拟值与实测值的定量关系，构成多准则的似然函数。

本发明实施例采用纳什效率系数作为似然度函数。

利用第i个参数组构建水动力模型，根据观测点处的水位、流量、含沙量和河床高程等观测值和观测点处的模拟值计算第i个参数组对应的似然度，计算公式如下：

式中，Z_it、Q_it、S_it和为第i个参数组所构建模型在t时刻的水位、流量、含沙量和河床高程模拟值；Z_ot、Q_ot、S_ot和为t时刻的水位、流量、含沙量观测值；和为水位观测值的均值，T为观测时间；L_z(i)为第z个观测点处水位观测值和模拟值之间的纳什效率系数，L_q(i)为第q个观测点处流量观测值和模拟值之间的纳什效率系数，L_s(i)为第s个观测点处水位观测值和模拟值之间的纳什效率系数，L_b(i)为第b个观测点处河床高程观测值和模拟值之间的纳什效率系数；num1、num2、num3和num4分别为水流、流量、含沙量和河床高程的观测点个数；L¹(i)、L²(i)、L³(i)和L⁴(i)分别为i个参数组对应的水位、流量、含沙量和河床高程的似然度；L(i)为加权似然度值，α₁、α₂、α₃和α₄分别表示L¹(i)、L²(i)、L³(i)和L⁴(i)的权重。

S106：设定似然度值阈值，将低于该阈值的似然度值对应的参数组删除。

设定似然度值的阈值，将似然度值低于该阈值的参数组进行剔除，将高于该阈值的参数组的似然度进行归一化，统计参数的后验分布，包括：

归一化处理：统计高于阈值的参数组所对应似然度的最大值L_max和最小值L_min，按下面公式对高于阈值的参数组所对应似然度值进行归一化处理：

式中，L_nor(i)为第i个参数组经归一化处理后的似然度。

统计参数后验分布：

式中，P(para_j|A)为参数j的后验分布；P(para_j)为参数j的先验分布；P(A|para_j)为似然度；k为模型参数个数。

S107：将剩余参数组对应的似然度值排序，确定不确定性量化区间。

删除部分参数组后，对高于该阈值的参数组对应的水位、泥沙含量、河流高程的值模拟值进行排序，确定5％和95％分位数，得到模拟水位、泥沙含量、河床底部高程90％不确定性区间。

S108：根据似然度值对模型参数进行敏感性分析。

利用GLUE累积似然度对模型参数敏感度性分析。如果模型参数取值对似然度值没有显著的影响，那么该参数累积概率密度分布应接近于原始分布，即均匀分布；如果模型参数取值对似然度值影响较大，则该参数累积概率密度分布与原始分布相差较大。本发明实施例采用Kolmogorov-Smirnov(K-S)统计量衡量似然度累计分布与原始分布的差别，其通过比较2个分布的最大垂直距离来计算2个分布的差异。2个分布之间的K-S距离越大，该参数的区域敏感度就越高。

D＝max|F_n(x)-F_m(x)|

式中，D为2个分布的距离；m、n为样本的个数；F_n、F_m为样本的分布函数。

从上述技术方案可以看出，本发明实施例提供了一种河流水沙模型的参数优化方法，包括为基于卫星遥感影像数据和河道的地形高程数据生成河道计算网格；采集河道待预测区域的初始场和边界信息，形成待预测区域的信息文件；根据待预测区域的信息文件，离散求解水沙控制模型；采用拉丁超立方采样方法LHS对水沙模型的多个参数在可行域内进行分层采样，得到多个参数组；基于模型参数组合，求解河流水沙模型，获得每个参数组对应的似然度值；设定似然度阈值，将低于预设阈值的似然度值所对应的参数组予以删除；将剩余的参数组所对应的似然度值进行排序，确定不确定性区间；利用GLUE累积对似然度值进行区域敏感度分析。通过上述处理，可以给出河流水沙模型中模型参数的最优值。从而解决了该模型因参数不正确表达而导致其预测精度较差的问题。

另外，将LHS和GLUE相结合应用于河流水沙模型不确定性分析，更加精确反映了输入概率函数中值的分布，可以显著提高抽样效率。且利用实时观测数据的输入，能够实时更新模型参数，提高模型预测精度，同时实时给出不确定性区间，为风险管理提供了基础。

需要说明的是，对于方法实施例，为了简单描述，故将其都表述为一系列的动作组合，但是本领域技术人员应该知悉，本发明实施例并不受所描述的动作顺序的限制，因为依据本发明实施例，某些步骤可以采用其他顺序或者同时进行。其次，本领域技术人员也应该知悉，说明书中所描述的实施例均属于优选实施例，所涉及的动作并不一定是本发明实施例所必须的。

实施例二

图2为本发明实施例提供的一种河流水沙模型的参数优化系统的结构框图；

参照图2所示，本发明实施例提供的参数优化系统用于对河流水沙模型进行参数优化，具体包括网格生成模块10、信息采集模块20、模型计算模块30、分层采样模块40、似然度计算模块50、参数组筛选模块60、不确定性分析模块70、敏感性分析模块80。

网格生成模块用于基于卫星遥感影像数据和地形高程数据计算河道高程。

该模块包括平面网格生成单元和高程插值计算单元，网格生成单元用于基于卫星遥感影像数据，提取水体边界，在此基础上生成河道多分辨率三角形网格；高程插值计算单元用于利用获取的河道的地形高程数据，插值生成河道三角形网格节点的河道高程。

信息采集模块用于采集河道的待预测区域的区域信息。

该河道的区域信息包括前述得到的河道高程E，还包括河宽B、上游边界的流量Q_up和含沙量S_up、下游边界的水位Z_down，各个断面的初始流量场Q_初始、初始水位场Z_初始、初始含沙量场S_初始。悬沙级配、床沙级配等信息。并将上述信息按照采样时间的先后顺序分别进行排列，得到一张河道的待预测区域的信息文件。

模型计算模块用于根据河道的待预测区域的信息文件求解水沙模型。

水沙模型方程包括水动力学控制方程、悬移质不平衡输移方程、推移质不平衡输移方程、床面变形方程、挟沙力计算公式等。

其中，动床条件下水动力学控制方程由二维浑水连续方程及x、y方向的运动方程组成，如下：

式中：t为时间；h为水深；u和v分别为x和y方向的水流流速；g为重力加速度；v_t为水流的紊动粘滞系数；Δρ＝ρ_s-ρ_f,ρ_s为泥沙密度，ρ_f为清水密度；ρ_m为浑水密度；S为总含沙量；ρ₀、ρ′分别为床沙的饱和湿密度及干密度；S_bx、S_by为床面底坡项，S_fx、S_fy为床面摩阻项，可分别表示为S_fx＝n²u(u²+v²)^1/2/h^4/3及S_fy＝n²v(u²+v²)^1/2/h^4/3，n为曼宁糙率系数；Z_b为河底高程。

二维悬移质泥沙的悬移质不平衡输移方程：

式中：ε_s为泥沙的紊动扩散系数；S_k、S_*k、ω_sk及α_sk分别为第k粒径组泥沙的含沙量、挟沙力、浑水沉速及恢复饱和系数。

二维推移质泥沙的推移质不平衡输移方程：

床面变形方程为由悬移质及推移质泥沙不平衡输移引起的床面变形方程，可分别表示为：

因此总的河床冲淤厚度可表示为：

式中，N为非均匀沙分组数；N_s为悬移质泥沙的分组数；ΔZ_sk、ΔZ_bk分别为悬移质及推移质泥沙在一个时间步长内引起的分组冲淤厚度；ΔZ_t为相应时段内总的床面冲淤厚度。因此时段末单元中心的河床底部高程可以表示为：(Z_b)^l+1＝(Z_b)^l+ΔZ_t，上标l指某一时间层。

挟沙力计算公式包括张瑞瑾公式、张红武公式、王兴奎公式、吴保生公式等。本发明实施例以吴保生公式为例：

式中，U为垂线平均流速，且U＝(u²+v²)^1/2；γ_s、γ_m分别为泥沙及浑水容重；ω_m为群体沉速。

推移质输沙力公式包括爱因斯坦公式、窦国仁公式等。本发明实施例以窦国仁提出的公式为例：

式中，q_b0k为第k粒径组推移质泥沙的单宽输沙力；K_b为一系数；U_ck为启动流速；C₀为无量纲的谢才系数，且C₀＝C/g^1/2，C＝h^1/6/n；ΔP_bk为床沙级配。故单位水体中推移质分组挟沙力可表示为q_b*k＝q_b0k/(hU)。

采用无结构三角网格下的有限体积法离散求解控制方程，同时采用单元中心的方式，即所有守恒变量值存储于单元中心，相邻两个单元的公共面为计算界面。具体计算方法见发明实施例一。

分层采样模块用于对水沙模型的多个模型参数进行分层采样，得到多个参数组。

具体为采用拉丁超立方采样LHS对水沙方程的模型参数中的水动力参数n、v_t，泥沙模型参数k、m、ε_s、α_sk、α_bk、k_b在可行域内进行分层采样，获取参数组。

利用拉丁超立方分层抽样方法分别将para₁,para₂,…,para_kk个参数空间范围划分为n个区间，在每个区间内均匀采样m次，每个参数共获取m×n个样本，将各个参数获取的m×n个样本进行随机组合，得到m×n个组合样本，每个样本代表一个参数组。

似然度计算模块用于利用水沙模型求解每个参数组对应的似然度值。

基于采样形成的水动力参数和泥沙参数构成的参数组，利用水沙模型进行求解，得到每个参数组对应的水位、流速、泥沙含量、河道高程等似然度值。

似然度函数是建立水位、流量、泥沙含量、河道高程等模拟值与实测值的定量关系，构成多准则的似然函数。

本发明实施例采用纳什效率系数作为似然度函数。

利用第i个参数组构建水动力模型，根据观测点处的水位、流量、含沙量和河床高程等观测值和观测点处的模拟值计算第i个参数组对应的似然度，计算公式如下：

式中，Z_it、Q_it、S_it和为第i个参数组所构建模型在t时刻的水位、流量、含沙量和河床高程模拟值；Z_ot、Q_ot、S_ot和为t时刻的水位、流量、含沙量观测值；和为水位观测值的均值，T为观测时间；L_z(i)为第z个观测点处水位观测值和模拟值之间的纳什效率系数，L_q(i)为第q个观测点处流量观测值和模拟值之间的纳什效率系数，L_s(i)为第s个观测点处水位观测值和模拟值之间的纳什效率系数，L_b(i)为第b个观测点处河床高程观测值和模拟值之间的纳什效率系数；num1、num2、num3和num4分别为水流、流量、含沙量和河床高程的观测点个数；L¹(i)、L²(i)、L³(i)和L⁴(i)分别为i个参数组对应的水位、流量、含沙量和河床高程的似然度；L(i)为加权似然度值，α₁、α₂、α₃和α₄分别表示L¹(i)、L²(i)、L³(i)和L⁴(i)的权重。

参数组筛选模块用于设定似然度值阈值，将低于该阈值的似然度值对应的参数组删除。

设定似然度值的阈值，将似然度值低于该阈值的参数组进行剔除，将高于该阈值的参数组的似然度进行归一化，统计参数的后验分布，包括：

归一化处理：统计高于阈值的参数组所对应似然度的最大值L_max和最小值L_min，按下面公式对高于阈值的参数组所对应似然度值进行归一化处理：

式中，L_nor(i)为第i个参数组经归一化处理后的似然度。

统计参数后验分布：

式中，P(para_j|A)为参数j的后验分布；P(para_j)为参数j的先验分布；P(A|para_j)为似然度；k为模型参数个数。

不确定性分析模块用于将剩余参数组对应的似然度值排序，确定不确定性量化区间。

删除部分参数组后，对高于该阈值的参数组对应的水位、泥沙含量、河流高程的值模拟值进行排序，确定5％和95％分位数，得到模拟水位、泥沙含量、河床底部高程90％不确定性区间。

敏感性分析模块用于根据似然度值对模型参数进行敏感性分析。

利用GLUE累积对似然度值进行区域敏感度分析。如果模型参数取值对似然度值没有显著的影响，那么该参数累积概率密度分布应接近于原始分布，即均匀分布；如果模型参数取值对似然度值影响较大，则该参数累积概率密度分布与原始分布相差较大。本发明实施例采用Kolmogorov-Smirnov(K-S)统计量衡量似然度累计分布与原始分布的差别，其通过比较2个分布的最大垂直距离来计算2个分布的差异。2个分布之间的K-S距离越大，该参数的区域敏感度就越高。

D＝max|F_n(x)-F_m(x)|

式中，D为2个分布的距离；m、n为样本的个数；F_n、F_m为样本的分布函数。

从上述技术方案可以看出，本发明实施例提供了一种河流水沙模型的参数优化系统，具体为基于卫星遥感影像数据和河道的地形高程数据生成河道计算网格；采集河道待预测区域的初始场和边界信息，形成待预测区域的信息文件；根据待预测区域的信息文件，离散求解水沙控制模型；采用拉丁超立方采样方法LHS对水沙模型的多个参数在可行域内进行分层采样，得到多个参数组；基于模型参数组合，求解河流水沙模型，获得每个参数组对应的似然度值；设定似然度阈值，将低于预设阈值的似然度值所对应的参数组予以删除；将剩余的参数组所对应的似然度值进行排序，确定不确定性区间；利用GLUE累积对似然度值进行区域敏感度分析。通过上述处理，可以给出河流水沙模型中模型参数的最优值。从而解决了该模型因参数不正确表达而导致其预测精度较差的问题。

另外，将LHS和GLUE相结合应用于河流水沙模型不确定性分析，更加精确反映了输入概率函数中值的分布，可以显著提高抽样效率。且利用实时观测数据的输入，能够实时更新模型参数，提高模型预测精度，同时实时给出不确定性区间，为风险管理提供了基础。

对于装置实施例而言，由于其与方法实施例基本相似，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。

本说明书中的各个实施例均采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似的部分互相参见即可。

本领域内的技术人员应明白，本发明实施例的实施例可提供为方法、装置、或计算机程序产品。因此，本发明实施例可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明实施例可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明实施例是参照根据本发明实施例的方法、终端设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理终端设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理终端设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理终端设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理终端设备上，使得在计算机或其他可编程终端设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程终端设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

尽管已描述了本发明实施例的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例做出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明实施例范围的所有变更和修改。

最后，还需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者终端设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者终端设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者终端设备中还存在另外的相同要素。

以上对本发明所提供的技术方案进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种河流水沙模型的参数优化方法，其特征在于，包括如下步骤：

基于卫星遥感影像数据和河道的地形高程数据计算生成计算网格；

采集所述河道的待预测区域的初始场和边界条件信息，并进行排列得到所述待预测区域的信息文件；

根据所述待预测区域的信息文件求解水沙模型；

采用拉丁超立方体采样对所述水沙模型的多个模型参数在可行域内进行分层采样，得到多个参数组；

利用所述模型参数组合，求解河流水沙模型得到每个所述参数组所对应的似然度值；

将所述多个参数组中，低于设定阈值的似然度值所对应的参数组予以删除；

将所述多个参数组中剩余的参数组所对应的似然度值进行排序，确定不确定性区间；

利用GLUE累积似然度值对模型参数进行区域敏感度分析。

2.如权利要求1所述的参数优化方法，其特征在于，所述基于卫星遥感影像数据和河道的地形高程数据生成计算网格，包括：

基于卫星遥感影像数据，提取水体边界，生成河道平面三角形网格；

利用河道地形高程数据对所述河道三角形网格进行插值运算，得到所述河道平面三角形网格的节点高程。

3.如权利要求1所述的参数优化方法，其特征在于，所述区域初始场和边界条件信息包括所述河道高程、河宽、上游边界的流量和含沙量、下游边界的水位、多个断面的初始流量场、初始水位场、初始含沙量场、悬沙级配、床沙级配中的部分或全部。

4.如权利要求1所述的参数优化方法，其特征在于，所述水沙模型为水动力学控制方程、悬移质不平衡输移方程、推移质不平衡输移方程、床面变形方程和挟沙力计算公式中的部分或全部。

5.如权利要求1所述的参数优化方法，其特征在于，所述模型参数包括水动力参数和泥沙参数，其中：

所述水动力参数包括糙率系数和/或紊动粘滞系数；

所述泥沙参数包括紊动扩散系数、悬移质泥沙恢复饱和系数、推移质泥沙恢复饱和系数、悬移质挟沙力经验参数和推移质挟沙力经验参数中的部分或全部。

6.一种河流水沙模型的参数优化系统，其特征在于，包括：

网格生成模块，用于基于卫星遥感影像数据和河道的地形高程数据生成计算网格；

信息采集模块，用于采集所述河道的待预测区域的初始场信息和边界条件信息，并进行排列得到所述待预测区域的信息文件；

模型计算模块，用于根据所述待预测区域的信息文件求解水沙模型；

分层采样模块，用于采用拉丁超立方体采样对所述水沙模型的多个模型参数在可行域内进行分层采样，得到多个参数组；

似然度计算模块，用于利用所述模型参数组合，求解河流水沙模型得到每个所述参数组所对应的似然度值；

参数组筛选模块，用于设定似然度阈值，将所述多个参数组中，低于预设阈值的似然度值所对应的参数组予以删除；

不确定性分析模块，用于将所述多个参数组中剩余的参数组所对应的似然度值进行排序，确定不确定性量化区间；

敏感度分析模块，用于利用GLUE累积的所述似然度值对模型参数进行区域敏感度分析。

7.如权利要求6所述的参数优化系统，其特征在于，所述网格生成模块包括：

平面网格生成单元，用于基于卫星遥感影像数据，提取水体边界，生成河道平面三角形网格；

插值计算单元，用于利用河道地形高程数据对所述河道平面三角形网格进行插值运算，得到所述河道平面三角形网格的节点高程。

8.如权利要求6所述的参数优化系统，其特征在于，所述区域初始场信息和边界条件信息包括所述河道高程、河宽、上游边界的流量和含沙量、下游边界的水位、多个断面的初始流量场、初始水位场、初始含沙量场、悬沙级配中的部分或全部。

9.如权利要求6所述的参数优化系统，其特征在于，所述水沙模型为水动力学控制方程、悬移质不平衡输移方程、推移质不平衡输移方程、床面变形方程和挟沙力计算公式中的部分或全部。

10.如权利要求6所述的参数优化系统，其特征在于，所述模型参数包括水动力参数和泥沙参数，其中：

所述水动力参数包括糙率系数和/或紊动粘滞系数；

所述泥沙参数包括紊动扩散系数、悬移质泥沙恢复饱和系数、推移质泥沙恢复饱和系数、悬移质挟沙力经验参数和推移质挟沙力经验参数中的部分或全部。