CN111241758A - 基于可溶性污染物在水环境中输移扩散模型的评估方法 - Google Patents

Abstract

本发明涉及基于可溶性污染物在水环境中输移扩散模型的评估方法，以圣维南方程和一维对流扩散方程为基础，设置汊点连接，建立一维河网汊点模型，通过Preissmann四点隐式格式离散转化为线性方程组，基于三级联解法的思想，应用追赶法，求得污染物浓度时空变化的递推关系，带入初始、边界条件，将未知数集中到汊点上，运用超松弛迭代法求解汊点未知数，然后带入各单一直流河段分别求解。针对概化断面数据支撑不足的情况，本发明创新性地采用Google Earth遥感与图像处理技术相结合，设计合理可行的数据提取方法，得到基础数据，并通过水力学公式，计算河网模型必要的水力参数；关于参数率定方面，采用改进的Bayesian‑MCMC方法对模型的参数进行率定。

Description

基于可溶性污染物在水环境中输移扩散模型的评估方法

技术领域

本发明涉及水环境评估领域，尤其是一种一维河网水动力水质模型数据的获取方法及模型针对具体河流区域的实际应用。

背景技术

水库是我国重要的淡水资源，随着社会经济的发展和人口数量的增多，给我国水库资源安全带来了巨大压力，水污染层出不穷。水库上游地区的河流污染直接影响着库区的水质安全，因而研究库区上游区域河网的污染物迁移和转化规律，对水库的日常治理及水污染突发事件的决策处理有着十分重要的意义。随着CFD的兴起与高速发展，国内外对水动力水质模型研究的学者颇多，模型的理论研究及数值计算已趋于成熟，对于污染物在特定区域的河网水系统中的迁移和转换规律的模拟，多建立一维河网水动力水质模型，模型的求解需将整个河网系统划分为大量的概化断面，因而需要庞大的基础数据支撑。

目前主流的方法是通过测量几个典型断面的实测数据，然后采用插值法获得全部断面的数据，而对于河道地形复杂的研究区域，建立监测站的成本太高，所以有实测数据的典型断面数量很少，结果导致插值得到数据与实际数据相差较大，使得模型的准确性不高。为解决这一难题，卢敏等人提出一种利用Google Earth遥感影像提取河流多期水面宽度，结合实测最大水深获取河流概化断面参数的方法，并分别用概化断面和实测断面在MIKE11下构建少资料河流一维水动力水质模型，对两种断面模型的计算结果进行对比分析。但只应用Google Earth遥感影像，提取的信息有限。

发明内容

本发明目的在于提供一种设计合理、计算精确、应用广泛的基于可溶性污染物在水环境中输移扩散模型的评估方法。

为实现上述目的，采用了以下技术方案：本发明所述方法包括以下步骤：

S1，数据收集，基于公共部门提供、查阅文献资料、实地勘测得到的部分实测数据，采用Google Earth遥感技术和图像处理技术，设计基础数据提取方法；

计算得到河道宽度、断面间距、河床海拔、河道形状、糙率、水深等基础数据，通过水力学相关的计算公式，分别计算湿周、流量、流量模数、断面面积、水位、水力半径、水力坡度等必要水力参数，为一维河网水动力水质模型提供数据支撑；

S2，应用系统工程的设计思想，将整个流域概化为一个由单一河段和汊点组成的系统，把上游河流化分为主干流和支流，并通过汊点连接；

本发明以桃林口水库上游河流流域为案例分析，进行河网系统的设计和模拟计算。桃林口位于温带大陆性季风区，降水主要集中在7、8两个月，上游河流以青龙河水系为主，除主干河流(青龙河)和几大支流外，其他河流的经流量均小于10m³/s，而且河流的径流量集中于汛期7、8月，其他时间段内流量较小，甚至出现断流的情况。

为简化问题，研究主要矛盾，本发明根据上游流域特点，忽略其他细小河流的影响，只选择汛期的青龙河主干道及4条径流量较大的支流作为主要的河流研究对象，将流域中的河流概化为9个河段和4个汊点，然后分别建立河段结构体和汊点结构体，通过汊点连接机制连接，进行河网系统的构建。并在每条河段选取多个典型断面(如图4)，在此基础上，运用分段样条插值，细化河段的划分。为满足伯努利方程的微分性质，本发明取500m的断面间距。

S3，运用欧拉法，将河段划分为多个研究断面，建立一维河网水动力模型对上游河流的水动力条件进行模拟，进而建立与水动力特性相对应的污染输移物扩散(水质)模型，计算污染输移物扩散的过程；

S4，一维河网水动力水质模型的求解，采用Preissmann四点隐式格式将河网中每条单一河段的水动力水质模型的控制方程离散为差分方程，并整理成线性方程组，通过汊点条件进行衔接；

基于基础数据与河流的自然条件，确定各河段的初始和边界条件。基于三级联解法的思想，运用追赶法，求得各断面污染物浓度的时空递推关系，将各河段的未知数集中在汊点上，根据汊点连接条件即汊点边界条件形成封闭的汊点方程组，应用超松弛迭代法求解此方程组，然后回代到单一河段中，最终得到各断面的河流流量、水位和污染物浓度等结果；

S5，模型参数率定，采用改进的Bayesian-MCMC方法，将模型参数的率定问题视为贝叶斯估计问题，根据有限差分方法和贝叶斯推理得到参数的后验概率密度函数，通过改进的Metropolis-Hastings抽样方法得到合理的参数值，提高模型的针对性和准确性。

进一步的，在步骤S1中，利用Google Earth遥感技术，获取河道两个相邻断面间的周边卫星遥感地图，并基于图像处理技术，对遥感地图进行灰度化、去燥、滤波等操作后，得到图像的二值化矩阵；(见图3)

其中，f(x_i,x_j)为坐标为(x_i,y_j)的像素点的灰度值；

首先，进行阈值调试和连通性检验，保留最大连通分支Ψ_max(本发明案例为最大连通分支，可根据具体的实际情况进行合理的调整)的同时，合理设置阈值m₀和M₀，使得当f(x_i,x_j)∈[m₀,M₀]∩f(x_i,x_j)∈Ψ_max时，像素点(x_i,x_j)在河流集合H＝{f₁,f₂,f₃,...,f_s,}中；

然后将提取的河流集合H与原图的河流进行相似度判断，如果满足相似度判断，输出集合H，如果不满足相似度判断，继续调整阈值m₀和M₀或对二值图像进行开闭运算等形态学操作，扩大、减小集合H，重复上述步骤，多次迭代，直到满足相似度判断，最终得到准确的河流集合H；

(1)河道宽度

提取在河流的像素点集合H中研究断面的河宽坐标，计算该河段第k个断面附近20组河宽的像素距离D_x ^k(20)，通过数据齐全的断面建立像素距离和实际距离之间的比例尺，计算该20组河宽的实际距离，用其均值来代表断面处的水面宽度；

(2)断面间距

假设两个断面间河流的河岸长度近似为断面间距，基于DFS算法思想，设置优先搜索方向和末断面坐标(搜索路径从首断面河岸点开始，当横坐标到达末断面时，搜索结束)，按照特定的搜索方向提取河岸的像素点集P＝{p₁,p₂,...,p_l}，通过比例尺，计算断面的间距；

(3)河床海拔

在Google Earth7.3版本中，有记录地形海拔的功能，在断面附近提取20组河道的海拔高度，用均值代替断面处的河床海拔；

(4)河道形状

将河道概化为抛物线与等腰梯形的连接，利用Google Earth遥感影像，提取的枯水期和丰水期的水面宽度并结合结合枯水期和丰水期的最大水深，通过插值拟合，即可得到各断面形状概化方程；

(5)糙率

糙率是影响水体动态特征的重要指标，天然山区河流河道糙率一般介于0.025～0.035之间，建立专家评价方法，以m个在糙率测量和水力计算等方面有相当经验的专家组成评价小组组，根据Google Earth遥感影像为断面处的糙率进行打分，得分记为n_i ^k,i＝1,2,...,m，并根据专家的权威性，得到专家的权值向量P＝(p₁,p₂,...,p_m)，由计算式

得到专家评价下的糙率

然后通过实测断面的糙率插值得到的第k个断面的糙率为

最后通过计算式

计算的第k个断面的糙率，其中θ为专家的置信系数；

(6)水深

以第k个断面为例，在断面附近的河道中心区域、两个岸边分别提取20组河床的海拔高度，记为h_0i ^k、h_li ^k、h_ri ^k，i＝1,2,...,20，并由公式

近似计算第k个断面的水深

进一步的，在步骤S3中，一维河网水动力水质模型的建立：

针对河段部分，基于河流水体在河道流动过程中质量守恒和动量守恒与污染物的质量守恒，分别建立一维圣维南方程组(水动力)和一维对流扩散方程(水质)；针对汊点部分，基于汊点处河流水体的质量守恒和能量守恒及污染物的质量守恒，分别建立水动力模型和水质模型的汊点连接条件；为方便计算机求解，采用Preissmann四点加权隐式差分格式，将连续方程离散化；

一、河网水动力模型

(1)河网水动力控制方程

描述明渠非恒定流的一维圣维南方程组为：

式中，x和t分别为空间和时间坐标，z为水位，Q为过水流量，B为过水宽度，A为过水断面面积，K为流量模数，g为重力加速度，q₁为旁侧入流流量；

(2)河网水动力方程离散

利用Preissmann四点加权隐式差分格式，将上述圣维南方程组离散得单一河段的差分方程组：

a1_jΔQ_j+b1_jΔz_j+c1_jΔQ_j+1+d1_jΔz_j+1＝e1_j (7)

a2_jΔQ_j+b2_jΔz_j+c2_jΔQ_j+1+d2_jΔz_j+1＝e2_j (8)

式中，a1_j、b1_j、c1_j、d1_j、e1_j、a2_j、b2_j、c2_j、d2_j、e2_j为时间步长Δt内河段断面j的差分方程的系数，Δz_j、Δz_j+1分别为第j，j+1断面在Δt时间内的水位增量，ΔQ_j、ΔQ_j+1分别为第j，j+1断面在Δt时间内的流量增量；

(3)河网水动力节点连接条件

流量守恒条件：进出某一汊点水量与该汊点实际水量增减平衡，表示如下：

式中，m为与汊点i相连的河段数，A_i为汊点i处的面积，z_i为汊点的水位，

为连接i汊点n时刻的流量之和；

能量守恒条件：连接汊点的各河段的水位增量与汊点的水位增量相同，表示如下：

Δz_i＝Δz_j，(i,j＝1,2,......,m) (10)

二、河网水质模型

(1)河网水质控制方程

式中，C为水流输送的水质变量浓度，E_x为污染物纵向离散系数，K₁为污染物衰减系数，Q为流量，A为断面面积。

(2)河网水质方程的离散

对方程(11)采用前差分离散时间项，隐式迎风格式离散对流项，中心差分离散扩散项；可得到三对角方程：

a_jC_j-1+b_jC_j+c_jC_j+1＝z_j，(j＝2,......,L2-1) (12)

(3)河网水质汊点连接条件

汊点处可给出质量平衡方程：

作为河网水质的连接条件。

进一步的，在步骤S4中，河网模型的求解：采用三级联解法，将河网系统拆分成单一河段及连接各个河段的汊点，在各河段上进行断面划分，在断面上将圣维南方程组离散化，利用追赶法得河段方程，辅以汊点连接条件形成以汊点水位为待求变量的汊点方程组，应用超松弛迭代法求解得到各汊点水位，然后将各汊点水位回代至各单一河段方程，最终求得各断面水位及流量；(模型具体的求解过程见图5)

一、水动力模型的求解：

根据河段是否连接外边界，将河段分成内河段和外河段。

对于内河段，子河段方程(7)、(8)经递推运算得如下形式方程：

ΔQ_j＝θ_j+η_jΔz_j+γ_jΔz₁，(j＝2,3，...,L2) (15)

其中，L2为该河段划分的断面数目，α_j，β_j，

θ_j，η_j，γ_j为追赶系数，可逐步递推求得；内河段总有两个汊点与之相连，首末断面水位各有一追赶方程：

ΔQ_L2＝θ_L2+η_L2Δz_L2+γ_L2Δz₁ (17)

对于外河段，首末断面关系有如下的线性方程组：

ΔQ_j＝F_jΔz_j+G_j (18)

Δz_j＝H_jΔQ_j+1+I_jΔz_j+1+J_j (19)

外河道用追赶法求解时，在追的过程中求得追赶系数H_j、I_j、J_j、F_j和G_j，而后在赶的过程中求出

和

同时给出边界条件，确定F₁和G₁初始值，进行求对单个汊点，建立其连接的内、外河段的边界方程，代入式(16)并与式(17)联立，得封闭的以汊点水位为未知量的方程组，代入边界条件，同理对其它汊点建立对应的汊点方程组，最终得到河网汊点水位方程组；应用超松弛迭代法求解该汊点水位方程组，再根据追赶法并结合初始条件逐步回代可求得各河段每个断面的水位和流量；

二、水质方程的求解：

三对角方程(12)的离散系数为：

各系数在顺流、逆流等不同流动类型时系数作相应变化；

汊点处可给出质量平衡方程：

式中，Ω是汊点处的水面面积，j是节点编号，i是与节点j相连的河段编号，NL是与节点j相连的河段总数；

与汊点相连的断面流向为流出汊点时，设该断面浓度等于汊点浓度，若该断面流向为流入汊点，则根据该断面所在河段的递推方程组获得该断面的浓度表达式；根据质量平衡方程，建立汊点方程，由上述方程代入各断面浓度的递推关系式，辅以边界条件，可得到包含整个河网汊点浓度的代数方程；通过超松弛迭代法可求得各汊点的水质浓度，根据流动方向的不同，选择不同的递推公式，计算得河段各个断面的水质浓度值。

进一步的，在步骤S5中，模型参数的率定：

在河网水动力水质模型建立求解过程中，出于抽象问题、简化计算、增强模型准确性和适用性等原因，引入了一些参数。基于未知参数的不确定性分布信息，在一定程度上为避免因“最优”参数失真而带来的决策风险，本发明将参数的获取作为贝叶斯估计问题，采用基于贝叶斯推理采用改进的马尔科夫链蒙特卡罗(MCMC)方法，通过构造合适的马尔科夫链进行抽样而使用蒙特卡罗方法进行积分计算，求得待求参数的后验概率分布及其统计特征值，进一步获取待求参数点估计来率定相关参数。相比较一般的优化方法，该方法能很好解决由观测数据噪声带来的非唯一解的问题，稳定性和准确性更高。

1)根据变量个数N及其部分先验信息，确定未知参数的样本空间和先验概率密度函数p(θ)；

2)在其先验范围内随机生成N个初始值X＝{x_i(1),x_i(2),x_i(3),....,x_i(n)}，并设定i＝1；

3)设定Proposal分布U(x_i(s)-step,x_i(s)+step)，并生成x'(s)，其中U表示均匀分布，step为随机游走的步长；

4)分别计算出x(s)和x'(s)对应的污染物浓度值Y和Y₀，即：B＝∑|Y-Y₀|；

5)如果B>0.6，则接受该测试参数并设定为当前模型参数，即x_i(s)＝x'(s)；否则不接受该测试参数，x_i(s)＝x(s)；

6)利用分布U(x_i(s)-step,x_i(s)+step)生成X^*＝{x^*(1),x^*(2),.....x^*(N)}；

7)计算能够反映模型参数和观测数据之间关系的似然函数p(θ|y)；

8)计算未知参数的后验概率密度p(θ|y)；

9)计算Markov链从X⁽ⁱ⁾位置移动到X^(*)；

10)产生一个0～1间均匀分布的随机数R，如果R＜A(X⁽ⁱ⁾,X^(*))，则接受该测试参数并设定为当前模型参数，即X⁽ⁱ⁺¹⁾＝X^(*)；否则，不接受该测试参数，X⁽ⁱ⁺¹⁾＝X⁽ⁱ⁾；

11)重复步骤1)到步骤10)，直至达到预定迭代次数。

与现有技术相比，本发明方法具有如下优点：基于一维河网水动力水质模型常用的求解过程，将Google Earth遥感和图像处理等技术结合，得到河网系统特殊处理后的灰度图，设计数据提取算法，以计算任意断面更多必要的水力数据。并根据桃林口水库上游流域的具体水文、气象条件，在桃林口库区上游河网进行了模型的创新性应用。

附图说明

图1是本发明的技术路线图。

图2为小菜峪断面附近的遥感图。

图3为图2经处理后的灰度图。

图4为桃林口水库上游河网系统的概化图。

图5为一维河网水动力水质模型的求解流程图。

具体实施方式

下面结合附图对本发明做进一步说明：

如图1所示，本发明所述方法包括以下步骤：

S1，数据收集，基于公共部门提供、查阅文献资料、实地勘测得到的部分实测数据，采用Google Earth遥感技术和图像处理技术，设计基础数据提取方法；

利用Google Earth遥感技术，获取河道两个相邻断面间的周边卫星遥感地图，并基于图像处理技术，对遥感地图进行灰度化、去燥、滤波等操作后，得到图像的二值化矩阵；如图3所示。

其中，f(x_i,x_j)为坐标为(x_i,y_j)的像素点的灰度值；

首先，进行阈值调试和连通性检验，保留最大连通分支Ψ_max(本发明案例为最大连通分支，可根据具体的实际情况进行合理的调整)的同时，合理设置阈值m₀和M₀，使得当f(x_i,x_j)∈[m₀,M₀]∩f(x_i,x_j)∈Ψ_max时，像素点(x_i,x_j)在河流集合H＝{f₁,f₂,f₃,...,f_s,}中；

然后将提取的河流集合H与原图的河流进行相似度判断，如果满足相似度判断，输出集合H，如果不满足相似度判断，继续调整阈值m₀和M₀或对二值图像进行开闭运算等形态学操作，扩大、减小集合H，重复上述步骤，多次迭代，直到满足相似度判断，最终得到准确的河流集合H；

(6)河道宽度

提取在河流的像素点集合H中研究断面的河宽坐标，计算该河段第k个断面附近20组河宽的像素距离

通过数据齐全的断面建立像素距离和实际距离之间的比例尺，计算该20组河宽的实际距离，用其均值来代表断面处的水面宽度；

(7)断面间距

假设两个断面间河流的河岸长度近似为断面间距，基于DFS算法思想，设置优先搜索方向和末断面坐标(搜索路径从首断面河岸点开始，当横坐标到达末断面时，搜索结束)，按照搜索方向提取河岸的像素点集P＝{p₁,p₂,...,p_l}，通过比例尺，计算断面的间距；

(8)河床海拔

在Google Earth7.3版本中，有记录地形海拔的功能，在断面附近提取20组河道的海拔高度，用均值代替断面处的河床海拔；

(9)河道形状

将河道概化为抛物线与等腰梯形的连接，利用Google Earth遥感影像，提取的枯水期和丰水期的水面宽度并结合结合枯水期和丰水期的最大水深，通过插值拟合，即可得到各断面形状概化方程；

(10)糙率

糙率是影响水体动态特征的重要指标，天然山区河流河道糙率一般介于0.025～0.035之间，建立专家评价方法，以m个在糙率测量和水力计算等方面有相当经验的专家组成评价小组组，根据Google Earth遥感影像为断面处的糙率进行打分，得分记为n_i ^k,i＝1,2,...,m，并根据专家的权威性，得到专家的权值向量P＝(p₁,p₂,...,p_m)，由计算式

得到专家评价下的糙率

然后通过实测断面的糙率插值得到的第k个断面的糙率为

最后通过计算式

计算的第k个断面的糙率，其中θ为专家的置信系数；

(6)水深

以第k个断面为例，在断面附近的河道中心区域、两个岸边分别提取20组河床的海拔高度，记为h_0i ^k、h_li ^k、h_ri ^k，i＝1,2,...,20，并由公式

近似计算第k个断面的水深

计算得到河道宽度、断面间距、河床海拔、河道形状、糙率、水深等基础数据，通过水力学相关的计算公式，分别计算湿周、流量、流量模数、断面面积、水位、水力半径、水力坡度等必要水力参数，为一维河网水动力水质模型提供数据支撑；

S2，应用系统工程的设计思想，将整个流域作为一个由单一河段和汊点组成的系统，把上游河流化分为主干流和支流，并通过汊点连接；

本发明以桃林口水库上游河流流域为案例分析，进行河网系统的设计和模拟计算。桃林口位于温带大陆性季风区，降水主要集中在7、8两个月，上游河流以青龙河水系为主，除主干河流(青龙河)和几大支流外，其他河流的经流量均小于10m³/s，而且河流的径流量集中于汛期7、8月，其他时间段内流量较小，甚至出现断流的情况。

为简化问题，研究主要矛盾，本发明根据上游流域特点，忽略其他细小河流的影响，只选择汛期的青龙河主干道及4条径流量较大的支流作为主要的河流研究对象，将流域中的河流概化为9个河段和4个汊点，然后分别建立河段结构体和汊点结构体，通过汊点连接机制连接，进行河网系统的构建。并在每条河段选取多个典型断面，如图4所示。在此基础上，运用分段样条插值，细化河段的划分。为满足伯努利方程的微分性质，本发明取500m的断面间距。

S3，运用欧拉法，将河段划分为多个研究断面，建立一维河网水动力模型对上游河流的水动力条件进行模拟，进而建立与水动力特性相对应的污染输移物扩散(水质)模型，计算污染输移物扩散的过程；

一维河网水动力水质模型的建立：

针对河段部分，基于河流水体在河道流动过程中质量守恒和动量守恒与污染物的质量守恒，分别建立一维圣维南方程组(水动力)和一维对流扩散方程(水质)；针对汊点部分，基于汊点处河流水体的质量守恒和能量守恒及污染物的质量守恒，分别建立水动力模型和水质模型的汊点连接条件；为方便计算机求解，采用Preissmann四点加权隐式差分格式，将连续方程离散化；

一、河网水动力模型

(1)河网水动力控制方程

描述明渠非恒定流的一维圣维南方程组为：

式中，x和t分别为空间和时间坐标，z为水位，Q为过水流量，B为过水宽度，A为过水断面面积，K为流量模数，g为重力加速度，q₁为旁侧入流流量；

(2)河网水动力方程离散

利用Preissmann四点加权隐式差分格式，将上述圣维南方程组离散得单一河段的差分方程组：

a1_jΔQ_j+b1_jΔz_j+c1_jΔQ_j+1+d1_jΔz_j+1＝e1_j (28)

a2_jΔQ_j+b2_jΔz_j+c2_jΔQ_j+1+d2_jΔz_j+1＝e2_j (29)

式中，a1_j、b1_j、c1_j、d1_j、e1_j、a2_j、b2_j、c2_j、d2_j、e2_j为时间步长Δt内河段断面j的差分方程的系数，Δz_j、Δz_j+1分别为第j，j+1断面在Δt时间内的水位增量，ΔQ_j、ΔQ_j+1分别为第j，j+1断面在Δt时间内的流量增量；

(3)河网水动力节点连接条件

流量守恒条件：进出某一汊点水量与该汊点实际水量增减平衡，表示如下：

式中，m为与汊点i相连的河段数，A_i为汊点i处的面积，z_i为汊点的水位，

为连接i汊点n时刻的流量之和；

能量守恒条件：连接汊点的各河段的水位增量与汊点的水位增量相同，表示如下：

Δz_i＝Δz_j，(i,j＝1,2,......,m) (31)

二、河网水质模型

(1)河网水质控制方程

式中，C为水流输送的水质变量浓度，E_x为污染物纵向离散系数，K₁为污染物衰减系数，Q为流量，A为断面面积。

(2)河网水质方程的离散

对方程(13)采用前差分离散时间项，隐式迎风格式离散对流项，中心差分离散扩散项；可得到三对角方程：

a_jC_j-1+b_jC_j+c_jC_j+1＝z_j，(j＝2,......,L2-1) (33)

各系数在顺流、逆流等不同流动类型时系数作相应变化；

(3)河网水质汊点连接条件

汊点处可给出质量平衡方程：

作为河网水质的连接条件。

S4，一维河网水动力水质模型的求解，采用Preissmann四点隐式格式将河网中每条单一河段的水动力水质模型的控制方程离散为差分方程，并整理成线性方程组，通过汊点条件进行衔接；

基于基础数据与河流的自然条件，确定各河段的初始和边界条件。基于三级联解法的思想，运用追赶法，求得各断面污染物浓度的时空递推关系，将各河段的未知数集中在汊点上，根据汊点连接条件即汊点边界条件形成封闭的汊点方程组，应用超松弛迭代法求解此方程组，回代到单一河段中，最终得到各断面的河流流量、水位和污染物浓度等信息；

河网模型的求解：采用三级联解法，将河网系统拆分成单一河段及连接各个河段的汊点，在各河段上进行断面划分，在断面上将圣维南方程组离散化，利用追赶法得河段方程，辅以汊点连接条件形成以汊点水位为待求变量的汊点方程组，应用超松弛迭代法求解得到各汊点水位，然后将各汊点水位回代至各单一河段方程，最终求得各断面水位及流量；模型具体的求解过程如图5所示。

一、水动力模型的求解：

根据河段是否连接外边界，将河段分成内河段和外河段；

对于内河段，子河段方程(28)、(29)经递推运算得如下形式方程：

ΔQ_j＝θ_j+η_jΔz_j+γ_jΔz₁，(j＝2,3，...,L2) (36)

其中，L2为该河段划分的断面数目，α_j，β_j，

θ_j，η_j，γ_j为追赶系数，可逐步递推求得；内河段总有两个汊点与之相连，首末断面水位各有一追赶方程：

ΔQ_L2＝θ_L2+η_L2Δz_L2+γ_L2Δz₁ (38)

对于外河段，首末断面关系有如下的线性方程组：

ΔQ_j＝F_jΔz_j+G_j (39)

Δz_j＝H_jΔQ_j+1+I_jΔz_j+1+J_j (40)

外河道用追赶法求解时，在追的过程中求得追赶系数H_j、I_j、J_j、F_j和G_j，而后在赶的过程中求出

和

同时给出边界条件，确定F₁和G₁初始值，进行求解；

对单个汊点，建立其连接的内、外河段的边界方程，代入式(37)并与式(38)联立，得封闭的以汊点水位为未知量的方程组，代入边界条件，同理对其它汊点建立对应的汊点方程组，最终得到河网汊点水位方程组；应用超松弛迭代法求解该汊点水位方程组，再根据追赶法并结合初始条件逐步回代可求得各河段每个断面的水位和流量。

二、水质方程的求解：

对方程(32)采用前差分离散时间项，隐式迎风格式离散对流项，中心差分离散扩散项，可得到三对角方程：

a_jC_j-1+b_jC_j+c_jC_j+1＝z_j，(j＝2,......,L2-1) (41)

三对角方程的离散系数为：

各系数在顺流、逆流等不同流动类型时系数作相应变化；

汊点处可给出质量平衡方程：

式中，Ω是汊点处的水面面积，j是节点编号，i是与节点j相连的河段编号，NL是与节点j相连的河段总数；

与汊点相连的断面流向为流出汊点时，设该断面浓度等于汊点浓度，若该断面流向为流入汊点，则根据该断面所在河段的递推方程组获得该断面的浓度表达式；根据质量平衡方程，建立汊点方程，由上述方程代入各断面浓度的递推关系式，辅以边界条件，可得到包含整个河网汊点浓度的代数方程；通过超松弛迭代法可求得各汊点的水质浓度，根据流动方向的不同，选择不同的递推公式，计算得河段各个断面的水质浓度值。

S5，模型参数率定，采用改进的Bayesian-MCMC方法，将模型参数的率定问题视为贝叶斯估计问题，根据有限差分方法和贝叶斯推理得到参数的后验概率密度函数，通过改进的Metropolis-Hastings抽样方法得到合理的参数值，提高模型的针对性和准确性。模型参数的率定：

在河网水动力水质模型建立求解过程中，出于抽象问题、简化计算、增强模型准确性和适用性等原因，引入了一些参数。基于未知参数的不确定性分布信息，在一定程度上为避免因“最优”参数失真而带来的决策风险，本发明将参数的获取作为贝叶斯估计问题，采用基于贝叶斯推理采用改进的马尔科夫链蒙特卡罗(MCMC)方法，通过构造合适的马尔科夫链进行抽样而使用蒙特卡罗方法进行积分计算，求得待求参数的后验概率分布及其统计特征值，进一步获取待求参数点估计来率定相关参数。相比较一般的优化方法，该方法能很好解决由观测数据噪声带来的非唯一解的问题，稳定性和准确性更高。

1)根据变量个数N及其部分先验信息，确定未知参数的样本空间和先验概率密度函数p(θ)；

2)在其先验范围内随机生成N个初始值X＝{x_i(1),x_i(2),x_i(3),....,x_i(n)}，并设定i＝1；

3)设定Proposal分布U(x_i(s)-step,x_i(s)+step)，并生成x'(s)，其中U表示均匀分布，step为随机游走的步长；

4)分别计算出x(s)和x'(s)对应的污染物浓度值Y和Y₀，即：B＝∑|Y-Y₀|；

5)如果B>0.6，则接受该测试参数并设定为当前模型参数，即x_i(s)＝x'(s)；否则不接受该测试参数，x_i(s)＝x(s)；

6)利用分布U(x_i(s)-step,x_i(s)+step)生成X^*＝{x^*(1),x^*(2),.....x^*(N)}；

7)计算能够反映模型参数和观测数据之间关系的似然函数p(θ|y)；

8)计算未知参数的后验概率密度p(θ|y)；

9)计算Markov链从X⁽ⁱ⁾位置移动到X^(*)；

10)产生一个0～1间均匀分布的随机数R，如果R＜A(X⁽ⁱ⁾,X^(*))，则接受该测试参数并设定为当前模型参数，即X⁽ⁱ⁺¹⁾＝X^(*)；否则，不接受该测试参数，X⁽ⁱ⁺¹⁾＝X⁽ⁱ⁾；

11)重复步骤1)到步骤10)，直至达到预定迭代次数。

通过Google Earth遥感影像，得到丰水期和枯水期的研究区域河网的遥感图，利用图像处理技术，对河网遥感图进行目标提取，得到河网的像素点集合，然后建立准确的图像像素坐标系与实际坐标系的比例尺和对基础水力数据提取方法，实现对河道宽度、断面间距、河床海拔、河道形状、糙率、水深等数据的计算，进而为一维河网水动力水质模型提供数据支撑。

以上所述的实施例仅仅是对本发明的优选实施方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。

Claims (5)

1.一种基于可溶性污染物在水环境中输移扩散模型的评估方法，其特征在于，所述方法包括以下步骤：

S1，数据收集，基于公共部门提供、查阅文献资料、实地勘测得到的部分实测数据，采用Google Earth遥感技术和图像处理技术，设计基础数据提取方法；

计算得到各过水断面河道宽度、断面间距、河床海拔、河道形状、糙率、水深等基础数据，通过水力学相关计算公式，分别计算湿周、流量、流量模数、断面面积、水位、水力半径、水力坡度等必要水力参数，为一维河网水动力水质模型提供数据支撑；

S2，应用系统工程的设计思想，将整个流域概化为一个由单一河段和汊点组成的系统，把上游河流化分为主干流和支流，并通过汊点连接；

S3，运用欧拉法，将河段划分为多个研究断面，建立一维河网水动力模型对上游河流的水动力条件进行模拟，进而建立与水动力特性相对应的污染输移物扩散模型，计算污染输移物扩散的过程；

S4，一维河网水动力水质模型的求解，采用Preissmann四点隐式格式将河网中每条单一河段水动力水质模型的控制方程离散为差分方程，并整理成线性方程组，通过汊点条件进行衔接；

根据基础数据与河流的自然条件，确定各河段的初始和边界条件。基于三级联解法的思想，运用追赶法，求得各断面污染物浓度的时空递推关系，将各河段的未知数集中在汊点上，依据汊点连接条件即汊点边界条件形成封闭的汊点方程组，应用超松弛迭代法求解此方程组，然后回代到单一河段中，最终得到各断面的河流流量、水位和污染物浓度等结果；

S5，模型参数率定，采用改进的Bayesian-MCMC方法，将模型参数的率定问题视为贝叶斯估计问题，根据有限差分方法和贝叶斯推理得到参数的后验概率密度函数，通过改进的Metropolis-Hastings抽样方法得到合理的参数值，提高模型的针对性和准确性。

2.根据权利要求1所述的基于可溶性污染物在水环境中输移扩散模型的评估方法，其特征在于：步骤S1中，利用Google Earth遥感技术，获取河道两个相邻断面间的周边卫星遥感地图，并基于图像处理技术，对遥感地图进行灰度化、去燥、滤波操作后，得到图像的二值化矩阵；

其中，f(x_i,x_j)为坐标为(x_i,y_j)的像素点的灰度值；

首先，进行阈值调试和连通性检验，保留最大连通分支Ψ_max的同时，合理设置阈值m₀和M₀，使得当f(x_i,x_j)∈[m₀,M₀]∩f(x_i,x_j)∈Ψ_max时，像素点(x_i,x_j)在河流集合H＝{f₁,f₂,f₃,...,f_s,}中；

然后将提取的河流集合H与原图的河流进行相似度判断，如果满足相似度判断，输出集合H，如果不满足相似度判断，继续调整阈值m₀和M₀或对二值图像进行插值、形态学开闭运算等操作，扩大、减小集合H，重复上述步骤，多次迭代，直到满足相似度判断，最终得到准确的河流集合H；

(1)河道宽度

提取在河流的像素点集合H中研究断面的河宽坐标，计算该河段第k个断面附近20组河宽的像素距离

通过数据齐全的断面建立像素距离和实际距离之间的比例尺，计算该20组河宽的实际距离，用其均值来代表断面处的水面宽度；

(2)断面间距

假设两个断面间河流的河岸长度近似为断面间距，基于DFS算法思想，设置优先搜索方向和末断面坐标，按照特定的搜索方向提取河岸的像素点集P＝{p₁,p₂,...,p_l}，通过比例尺，计算断面的间距；

(3)河床海拔

在Google Earth7.3版本中，在断面附近提取20组河道的海拔高度，用均值代替断面处的河床海拔；

(4)河道形状

将河道概化为抛物线与等腰梯形的连接，利用Google Earth遥感影像，提取的枯水期和丰水期的水面宽度并结合结合枯水期和丰水期的最大水深，通过插值拟合，即可得到各断面形状概化方程；

(5)糙率

糙率是影响水体动态特征的重要指标，天然山区河流河道糙率一般介于0.025～0.035之间，建立专家评价方法，以m个在糙率测量和水力计算等方面有相当经验的专家组成评价小组组，根据Google Earth遥感影像为断面处的糙率进行打分，得分记为n_i ^k,i＝1,2,...,m，并根据专家的权威性，得到专家的权值向量P＝(p₁,p₂,...,p_m)，由计算式

得到专家评价下的糙率

然后通过实测断面的糙率插值得到的第k个断面的糙率为

最后通过计算式

计算的第k个断面的糙率，其中θ为专家的置信系数；

(6)水深

以第k个断面为例，在断面附近的河道中心区域、两个岸边分别提取20组河床的海拔高度，记为h_0i ^k、h_li ^k、h_ri ^k，i＝1,2,...,20，并由公式

近似计算第k个断面的水深

3.根据权利要求1所述的基于可溶性污染物在水环境中输移扩散模型的评估方法，其特征在于：在步骤S3中，一维河网水动力水质模型的建立：

针对河段部分，基于河流水体在河道流动过程中质量守恒和动量守恒与污染物的质量守恒，分别建立一维圣维南方程组和一维对流扩散方程；针对汊点部分，基于汊点处河流水体的质量守恒和能量守恒及污染物的质量守恒，分别建立水动力模型和水质模型的汊点连接条件；采用Preissmann四点加权隐式差分格式，将连续方程离散化；

一、河网水动力模型

(1)河网水动力控制方程

描述明渠非恒定流的一维圣维南方程组为：

式中，x和t分别为空间和时间坐标，z为水位，Q为过水流量，B为过水宽度，A为过水断面面积，K为流量模数，g为重力加速度，q₁为旁侧入流流量；

(2)河网水动力方程离散

利用Preissmann四点加权隐式差分格式，将上述圣维南方程组离散得单一河段的差分方程组：

a1_jΔQ_j+b1_jΔz_j+c1_jΔQ_j+1+d1_jΔz_j+1＝e1_j (7)

a2_jΔQ_j+b2_jΔz_j+c2_jΔQ_j+1+d2_jΔz_j+1＝e2_j (8)

式中，a1_j、b1_j、c1_j、d1_j、e1_j、a2_j、b2_j、c2_j、d2_j、e2_j为时间步长Δt内河段断面j的差分方程的系数，Δz_j、Δz_j+1分别为第j，j+1断面在Δt时间内的水位增量，ΔQ_j、ΔQ_j+1分别为第j，j+1断面在Δt时间内的流量增量；

(3)河网水动力节点连接条件

流量守恒条件：进出某一汊点水量与该汊点实际水量增减平衡，表示如下：

式中，m为与汊点i相连的河段数，A_i为汊点i处的面积，z_i为汊点的水位，

为连接i汊点n时刻的流量之和。

能量守恒条件：连接汊点的各河段的水位增量与汊点的水位增量相同，表示如下：

Δz_i＝Δz_j，(i,j＝1,2,......,m) (10)

二、河网水质模型

(1)河网水质控制方程

式中，C为水流输送的水质变量浓度，E_x为污染物纵向离散系数，K₁为污染物衰减系数；Q为流量，A为断面面积，I为水力坡降，B为河道宽度，h为断面平均水深；

(2)河网水质方程的离散

对方程(13)采用前差分离散时间项，隐式迎风格式离散对流项，中心差分离散扩散项；可得到三对角方程：

a_jC_j-1+b_jC_j+c_jC_j+1＝z_j，(j＝2,......,L2-1) (12)

(3)河网水质汊点连接条件

汊点处可给出质量平衡方程：

作为河网水质的连接条件。

4.根据权利要求1所述的基于可溶性污染物在水环境中输移扩散模型的评估方法，其特征在于，在步骤S4中，河网模型的求解：采用三级联解法，将河网系统拆分成单一河段及连接各个河段的汊点，在各河段上进行断面划分，在断面上将圣维南方程组离散化，利用追赶法得河段方程，辅以汊点连接条件形成以汊点水位为待求变量的汊点方程组，应用超松弛迭代法求解得到各汊点水位，然后将各汊点水位回代至各单一河段方程，最终求得各断面水位及流量；

一、水动力模型的求解：

根据河段是否连接外边界，将河段分成内河段和外河段。

对于内河段，子河段方程(7)、(8)经递推运算得如下形式方程：

ΔQ_j＝θ_j+η_jΔz_j+γ_jΔz₁，(j＝2,3，...,L2) (15)

其中，L2为该河段划分的断面数目，α_j，β_j，

θ_j，η_j，γ_j为追赶系数，可逐步递推求得；内河段总有两个汊点与之相连，首末断面水位各有一追赶方程：

ΔQ_L2＝θ_L2+η_L2Δz_L2+γ_L2Δz₁ (17)

对于外河段，首末断面关系有如下的线性方程组：

ΔQ_j＝F_jΔz_j+G_j (18)

Δz_j＝H_jΔQ_j+1+I_jΔz_j+1+J_j (19)

外河道用追赶法求解时，在追的过程中求得追赶系数H_j、I_j、J_j、F_j和G_j，而后在赶的过程中求出

和

同时给出边界条件，确定F₁和G₁初始值，进行求解；

对单个汊点，建立其连接的内、外河段的边界方程，代入式(16)并与式(17)联立，得封闭的以汊点水位为未知量的方程组，代入边界条件，同理对其它汊点建立对应的汊点方程组，最终得到河网汊点水位方程组；应用超松弛迭代法求解该汊点水位方程组，再根据追赶法并结合初始条件逐步回代可求得各河段每个断面的水位和流量；

二、水质方程的求解：

对方程(11)采用前差分离散时间项，隐式迎风格式离散对流项，中心差分离散扩散项，可得到三对角方程：

a_jC_j-1+b_jC_j+c_jC_j+1＝z_j，(j＝2,......,L2-1) (20)

三对角方程的离散系数为：

各系数在顺流、逆流等不同流动类型时系数作相应变化；

汊点处可给出质量平衡方程：

式中，Ω是汊点处的水面面积，j是节点编号，i是与节点j相连的河段编号，NL是与节点j相连的河段总数；

与汊点相连的断面流向为流出汊点时，设该断面浓度等于汊点浓度，若该断面流向为流入汊点，则根据该断面所在河段的递推方程组获得该断面的浓度表达式；根据质量平衡方程，建立汊点方程，由上述方程代入各断面浓度的递推关系式，辅以边界条件，可得到包含整个河网汊点浓度的代数方程；通过超松弛迭代法可求得各汊点的水质浓度，根据流动方向的不同，选择不同的递推公式，计算得河段各个断面的水质浓度值。

5.根据权利要求1所述的基于可溶性污染物在水环境中输移扩散模型的评估方法，其特征在于，在步骤S5中，模型参数的率定：

1)根据变量个数N及其部分先验信息,确定未知参数的样本空间和先验概率密度函数p(θ)；

2)在其先验范围内随机生成N个初始值X＝{x_i(1),x_i(2),x_i(3),....,x_i(n)}，并设定i＝1；

3)设定Proposal分布U(x_i(s)-step,x_i(s)+step),并生成x'(s),其中U表示均匀分布,step为随机游走的步长；

4)分别计算出x(s)和x'(s)对应的污染物浓度值Y和Y₀,即：B＝∑|Y-Y₀|；

5)如果B>0.6,则接受该测试参数并设定为当前模型参数,即x_i(s)＝x'(s)；否则不接受该测试参数,x_i(s)＝x(s)；

6)利用分布U(x_i(s)-step,x_i(s)+step)生成X^*＝{x^*(1),x^*(2),.....x^*(N)}；

7)计算能够反映模型参数和观测数据之间关系的似然函数p(θ|y)；

8)计算未知参数的后验概率密度p(θ|y)；

9)计算Markov链从X⁽ⁱ⁾位置移动到X^(*)；

10)产生一个0～1间均匀分布的随机数R,如果R＜A(X⁽ⁱ⁾,X^(*))，则接受该测试参数并设定为当前模型参数,即X⁽ⁱ⁺¹⁾＝X^(*)；否则,不接受该测试参数,X⁽ⁱ⁺¹⁾＝X⁽ⁱ⁾；

11)重复步骤1)到步骤10)，直至达到预定迭代次数。

Images