基于超启发式-同伦算法的地下水有机污染源反演方法

1.本发明涉及地下水数值模拟反问题技术领域，尤其涉及一种基于超启发式-同伦算法的地下水有机污染源反演方法。


背景技术：

2.有机污染物具有低水溶性、高毒性和高界面张力等特性，进入地下水系统后会聚集滞留在含水层的顶部或底部(聚集位置取决于污染物密度小于水或大于水)，在与水接触的过程中不断向水中溶解释放，造成严重且持久的污染。在有机污染的实际修复过程中，常常面临污染物去除率低、修复过程时耗长和修复费用昂贵的困难。因此，制定合理高效的修复方案就显得尤为重要。而合理高效修复方案的制定则要以对于含水层中污染源状况的识别和掌握为前提。
3.但地下水埋藏于地下岩土介质之中，地下水污染通常具有存在的隐蔽性和发现的滞后性特点，致使人们对于地下水污染源的状况都缺乏了解和掌握，而且有机污染物在地下水系统中迁移过程涉及的许多关键参数难以通过目前的测量手段直接获取。针对这一问题，利用实际观测数据，结合数据同化方法的反演计算是目前主要的解决方案。


技术实现要素：

4.本发明要解决的技术问题是提供一种基于超启发式-同伦算法的地下水有机污染源反演方法，可以用于识别地下水有机污染源特征和估计多相流运移模型中的关键参数，从而实现地下水有机污染多相流运移过程的准确模拟预测，为地下水污染修复方案的合理制定、污染风险评估和污染责任认定提供重要前提基础条件。
5.为解决上述技术问题，本发明的技术方案为：1.一种基于超启发式
‑ꢀ
同伦算法的地下水有机污染源反演方法，其特征在于：具体方法包括如下步骤：
6.s1：通过野外现场调查、动态监测手段，掌握污染场地的水文地质条件，获得地下水水质动态监测数据；
7.s2：对场地水文地质条件进行概化处理，初步建立地下水有机污染多相流数值模拟模型，用以描述地下水中有机污染物的运移机理；
8.s3：确定污染源反演识别问题中的待识别的污染源特征及含水层参数；根据模型中待识别变量的取值范围，随机采样若干组样本，逐一代入s2中所建立的多相流数值模型，得到每组样本对应的模型响应，形成由“模型输入-模型响应”样本对构成的训练样本集和检验样本集；
9.s4：根据s3获得的输入-输出训练样本样集，采用不同机器学习方法建立多相流数值模型的单一替代模型；
10.s5：应用群智能优化算法，根据s3获得的检验样本集优化确定s4建立的各单一替代模型的权重及核心参数，建立群智能集成学习替代模型；
11.s6：建立污染源特征及含水层参数反演识别的非线性规划优化模型，并将s5建立
的群智能集成学习替代模型作为约束条件之一嵌入到优化模型之中；
12.s7：基于同伦理论，改写s6建立的优化模型，获得一系列反演同伦优化模型；
13.s8：采用超启发式算法依次求解s7建立的同伦优化模型系列，最后一个优化模型的解即为地下水污染源的反演识别结果。
14.进一步的，所述s2中用utchem程序构建地下水有机污染多相流运移模型。
15.进一步的，所述s3中待识别的污染源特征包括：污染源的纵向坐标、污染源的横向坐标、污染物迁移转化时长、污染物泄漏量；待识别的含水层参数包括：孔隙度、渗透率、纵向水相弥散度、横向水相弥散度；随机采样通过拉丁超立方采样方法实现。
16.进一步的，所述s4中的单一替代模型建模方法包括：克里金法、支持向量回归法、小波核极限学习机法；其中，克里金建模方法如下：克里金法(kriging)的回归方程可以表示为：
[0017][0018]
式中，f(x)＝[f1(x),f2(x),k,fk(x)]
t
为已知的基函数，β＝(β1,β2,k,βk)
t
为基函数对应的回归参数，由训练样本估计获得；z(x)为局部偏差项，均值为0，方差为其协方差可表示为：
[0019][0020]
r(u,v)为n维向量u和v的相关函数：
[0021][0022]
式中，αi为待定参数，ui和vi为u和v的第i个元素；
[0023]
对于给定的样本p＝[p1,p2,k,pm]
t
和相应的输出响应q＝[q1,q2,k,qm]
t
，任意输入向量x的预测输出响应y(x)可以写成：
[0024][0025]
式中，r是x与样本p之间的相关向量：
[0026]
r(x)＝[r(x,p1),r(x,p2),k,r(x,pm)]
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0027]
f为样本p的响应列向量：
[0028][0029]
r为样本p各样本点之间的相关矩阵：
[0030][0031]
是由广义最小二乘法估计得到β的估计值：
[0032][0033]
支持向量回归(svr)建模方法如下：
[0034]
svr是通过非线性映射函数把输入数据映射到高维空间进行线性回归，在训练拟合精度和预测精度之间进行权衡；对于给定的训练输入 x＝[x1,x2,k,xm]
t
每个要素代表一个n元输入：
[0035]
xi＝(x
i,1
,x
i,2
,k,x
i,n
),i＝1,2km和输出y＝(y1,y2,k,ym)
t
；非线性回归方程可表示为：
[0036]
f(x)＝《w,φ(x)》+b
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0037]
式中，w＝(w1,w2,k,wn)为权值向量，b为拟合误差，《w,x》表示w和φ(x) 的内积；φ(x)为非线性映射函数，可将输入向量从输入空间映射到高维空间；该问题可表示为如下优化问题：
[0038]
最小化
[0039]
约束条件
[0040]
式中，||w||2是w的范数，常数c及变量ξi和为惩罚因子和松弛变量，ε为可接受偏差；
[0041]
构造lagrange函数：
[0042][0043]
式中，l为拉格朗日算子，ηi、αi和为不小于0的拉格朗日乘子；对于最优解，和均为0：
[0044]
[0045][0046][0047][0048]
将以上性质与式11结合，可将式10中的优化问题改写为对偶形式：
[0049]
最小化
[0050]
约束条件
[0051]
由式16可知：
[0052][0053]
因此，可得到回归函数：
[0054][0055]
应用核函数代替内积《φ(x),φ(x')》，最为常用的为高斯函数：
[0056][0057]
式16可改写为：
[0058]
最小化
[0059]
约束条件
[0060]
最后得到：
[0061][0062]
拟合误差b可以利用karush-kuhn-tucker(kkt)条件计算；
[0063]
小波核极限学习机(wkelm)建模方法如下：
[0064]
对于训练样本(xj,tj),j＝1,k,n，极限学习机的输出yj可表示为：
[0065]
[0066]
q(xj)＝[p(ω1xj+b1),p(ω2xj+b2),k,p(ω
l
xj+b
l
)]
t
ꢀꢀ
(23)
[0067]
式中，p(
·
)为激励函数，输入节点通过权重向量ωi与第i个隐含神经元连接，i＝1,k,l，bi为第i个隐含神经元的阈值，p(ωixj+bi)为第 i个隐含神经元的输出函数，βi为第i个隐含神经元与输出神经元连接的权重向量；
[0068]
式22可表示为如下形式：
[0069]
qβ＝y
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(24)
[0070]
式中，β＝[β1,k,β
l
]
t
，y＝[y1,k,yn]
t
，q为elm的隐含层输出矩阵：
[0071]
如果具有l个隐含节点的极限学习机模型能够无偏的学习n个训练样本，那么存在下式：
[0072][0073]
式中，tj表示目标值；
[0074]
式26可以简写为：
[0075]
qβ＝t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(27)
[0076]
最小二乘法解决式27
[0077]
β＝q
+
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(28)
[0078]
式中，q
+
为q的moore
–
penrose广义逆；q
+
由下式计算
[0079]q+
＝q
t
(qq
t
)-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(29)
[0080]
对于训练样本(xj,tj),j＝1,k,n，kelm的原始优化问题表达为
[0081][0082]
s.t.q(xj)
t
·
β＝t
j-ξjꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(30)
[0083]
式中，c代表正则化参数，能够平衡训练误差和算法复杂性，ξj代表误差；
[0084]
优化问题可转化为拉格朗日对偶形式求解：
[0085][0086]
式中，θi代表拉格朗日算子；该对偶问题可运用karush-kuhn-tucker (kkt)条件求解：
[0087][0088]
[0089][0090]
根据式32可计算权重向量β的最小二乘解：
[0091][0092]
依据核函数理论，构建小波隐式映射函数k(xi,xj)代替随机映射函数 q(xj)：
[0093]kelm
＝qq
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(34)
[0094]kelm(i,j)
＝q(xi)
t
·
q(xj)＝k(xi,xj)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(35)
[0095]
训练后的kelm输出函数表达如下：
[0096][0097]
式中，t＝[t1,k,tn]
t
；
[0098]
小波核函数的表达式为：
[0099][0100]
式中αw,γ
1,w
,γ
2,w
为可调节参数。
[0101]
进一步的，所述s5中的集成学习替代模型是s4中的单一替代模型加权线性叠加构成，其输出的表达式如下：
[0102][0103]
式中，为集成替代模型的预测值，分别为克里金法模型，支持向量回归模型与小波核极限学习机模型的预测值，w1，w2，w3为权重值，权重之和为1；
[0104]
构造集成学习替代模型的一个关键环节是集成权重值和各单一替代模型核心参数的确定，替代模型的预测精度越高，其权重应越大；替代模型的精度越差，其权重应越小；采取以下方式作为集成权重和模型参数的确定方法；
[0105]
首先，建立如下优化模型，优化模型的决策变量为集成学习替代模型中的权重值和各单一替代模型的可变参数值，优化模型的最优解使得运用检验样本获得的集成学习替代模型输出与模拟模型输出之间的均方根误差最小；
[0106]
[0107][0108]
式中，yi(xk,pi)为第i个替代模型以参数向量pi和第k个检验样本输入xk为输入对应的输出响应，i＝1,2,k,n，k＝1,2,k,m，y
actual
(xk)为第 k个检验样本输入xk对应的模拟模型输出响应，m为检验样本的个数；应用群智能优化算法求解上述优化模型，获得集成权重和各单一替代模型可变参数的最优值，构建群智能集成学习替代模型。
[0109]
进一步的，所述s6中的非线性优化模型以监测井的模拟计算浓度值与实测浓度值的绝对偏差的平方和达到最小为目标函数；以地下水污染源反演识别问题中的待识别变量为决策变量，包括：污染源的纵向坐标、污染源的横向坐标、污染物迁移转化时长、污染物泄漏量、孔隙度、渗透率、纵向水相弥散度、横向水相弥散度；以群智能集成学习智能替代模型为污染物迁移转化规律的等式约束，以污染源特征和含水层参数取值范围为不等式约束条件；其表达式如下：
[0110][0111][0112]ck
为第k口监测井中污染物浓度的实际监测值，k＝1,2,k,n，为对应的模拟计算值；为污染物浓度的模拟计算值向量，即群智能集成学习机智能替代模型的输出，s为污染源相关变量取值向量，p为含水层参数取值向量；约束条件分别表示污染源纵向坐标(xi)、污染源横向坐标(yi)、污染物迁移转化时长(ti)、污染物泄漏量(vi)、孔隙度(θ)、渗透率(k)、纵向油相弥散度(d
oil,l
)、横向油相弥散度(d
oil,t
) 变化范围的上界及下界。
[0113]
进一步的，所述步骤7中的同伦优化模型系列构建首先需要借助同伦算法思想，引入同伦参数λ，构造同伦函数h(s,p,λ),使得当λ＝0 时，方程组h(s,p,λ)＝0的解为假设的地下水污染源及含水层参数，而当λ＝1时，方程组h(s,p,λ)＝0的解为待求的真实的地下水污染源及含水层参数；从任意假设的污染源信息及含水层参数值出发，通过路径跟踪，依次对一系列同伦方程组所对应的优化问题进行求解，逐步趋近并最终获得待求的真实的污染源信息及含水层参数值；
[0114]
同伦函数取线性同伦如下式：
[0115]
h(s,p,λ)＝λf(s,p)+(1-λ)g(s,p)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(40)
[0116]
其中，f(s,p)＝f(s,p)-c
obs
，g(s,p)＝f(s,p)-c0；
[0117]
式中，s为污染源特征相关变量取值向量，p为含水层参数取值向量；λ为同伦参数，取值范围[0,1]；f()表示多相流模拟模型或其替代模型；c
obs
为监测点实际监测浓度向量；c0表示将任意假设的地下水污染源特征代入模拟模型得到的监测点计算浓度向量；
[0118]
同伦函数h连续依赖于同伦参数λ，将λ的取值范围区间[0,1]进行划分λ0＝0＜λ1＜
…
＜λn＝1，得到方程组系列：
[0119]
h(s,p,λi)＝f(s,p)-(λi·cobs
+(1-λi)
·
c0)＝0,i＝1,2,k,n
ꢀꢀ
(41)
[0120]
如果λ
i+1-λi充分小，则相邻方程的解(s,p)
i+1
与(s,p)i是非常接近的；依据同伦算法思想，s6中的优化模型进行改写，可得到本文地下水污染源反演识别问题的同伦方程组系列所对应的优化模型系列，如式 (42)所示；表示将任意假设的地下水污染源特征及含水层参数值代入模拟模型得到的第k口监测井中的污染物计算浓度；
[0121][0122][0123]
进一步的，所述s7中的超启发式算法采用随机选择作为高层策略，即每次求解从给定的集合中随机选择低层启发式算法，集合中包括遗传算法、粒子群优化算法和模拟退火算法；
[0124]
运用超启发式算法(hha)依次对同伦优化模型系列进行求解，在求解过程，由于路径跟踪过程是逐渐演变的，因此前后相邻的两个优化问题的解也是十分接近的：
[0125][0126]
式中，为第i个同伦优化模型的初始种群，(s
i,opt
,p
i,opt
)为应用超启发式算法求得的第i个同伦优化模型的最优解，可作为生成下一优化模型初始种群的依据：
[0127]
[0128][0129][0130]
式中，v
s,
l
ow
，v
s,up
，v
p,
l
ow
，v
p,up
，表示超启发式算法中污染源特征与含水层参数对应的调整位移向量的上限和下限；
[0131]
应用第i个优化问题的解作为求解第i+1个优化问题初始值的生成依据，从而保证了寻优过程是逐渐演变的，避免因初始值选择距最优解较远时而产生的早熟收敛问题；最后一个优化模型的解即为地下水污染源的反演识别结果。
[0132]
本发明的优点在于：
[0133]
1、本发明中本发明建立了地下水有机污染多相流数值模型的群智能集成学习替代模型，提升了反演的计算效率。
[0134]
2、本发明构建了反演优化求解的超启发式-同伦算法，降低了反演的计算误差。
[0135]
3、本发明解决了地下水有机污染源特征的识别及多相流运移模型中关键参数的校正问题，从而实现地下水有机污染时空分布的准确模拟预测，为地下水污染修复方案的合理制定、污染风险评估和污染责任认定提供重要前提基础条件。
附图说明
[0136]
下面结合附图和具体实施方式对本发明作进一步详细的说明。
[0137]
图1为s7、s8中超启发式-同伦算法求解优化模型的流程架构。
[0138]
图2为实施例场地概化及水质监测井相对位置图。
[0139]
图3为s4、s5建立的各单一替代模型与集成学习智能替代模型检验输出与数值模拟模型输出的拟合散点图。
[0140]
图4为超启发式-同伦算法优化求解过程中不同变量识别值收敛曲线。
[0141]
图5为实施例实际污染羽与识别污染羽对比图：(a)实际污染羽；(b)识别污染羽。
具体实施方式
[0142]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
[0143]
因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0144]
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。
[0145]
在本发明的描述中，需要说明的是，术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，或者是该
发明产品使用时惯常摆放的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
[0146]
此外，术语“水平”、“竖直”等术语并不表示要求部件绝对水平或悬垂，而是可以稍微倾斜。如“水平”仅仅是指其方向相对“竖直”而言更加水平，并不是表示该结构一定要完全水平，而是可以稍微倾斜。
[0147]
在本发明的描述中，还需要说明的是，除非另有明确的规定和限定，术语“设置”、“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
[0148]
实施例1：一个有机污染潜水含水层，可概化为均质、各向同性的三维多相流模型；污染场地附近无天然边界，将边界划定在污染物迁移影响可以忽略不计的位置；其中，将东北边界、西南边界概化为一类边界；东南边界、西北边界由流面组成，概化为零通量边界；计算模拟区的下部为隔水层，可概化为零通量边界，上部为潜水面，是水量交换边界，因含水层厚度沿地下水流方向变化缓慢，故概化为等厚含水层；水和有机污染物氯苯的物理化学参数详见表1；
[0149]
表1水和有机污染物的物理化学参数
[0150][0151]
场地内的有机污染是由点状污染源造成的，污染物在短时间内进入含水层，而之后的大部分时间都是污染物在含水层中的自然溶解扩散阶段；因此，模拟模型中待识别的变量(污染源反演识别问题中待识别的变量)可最终确定为：1、污染源的横向、纵向坐标(污染源位置)；污染源默认位于含水层顶部，垂向坐标无需识别；2、污染物迁移转化时长(模拟期)；3、污染物泄漏量(污染物迁移转化时长及污染物泄漏量为该污染源反演识别问题中的
污染源释放历史)；4、含水层参数(包括孔隙度、渗透率、纵向水相弥散度、横向水相弥散度)；利用场地内的5口水位水质监测井获取污染物浓度监测数据作为反演求解的已知信息；场地概化情况及监测井位置见图2。
[0152]
实施例中给定的污染源信息及含水层参数的实际值见表2；
[0153]
表2实施例中污染源信息及含水层参数
[0154][0155]
基于上述信息，用utchem进行模型构建，建立关于实施例的地下水dnapls污染多相流模拟模型，并运用所建立的模拟模型做正演预报计算，获得渗流场内污染物浓度的时空分布状况，即监测井处的污染物浓度作为实际监测数据，末时刻五口监测井位置含水层底部的污染物浓度见表3；
[0156]
表3实施例正演模型计算输出(实际监测数据)
[0157][0158]
使用本发明的方法，根据表3中浓度监测数据对各个污染源特征及含水层参数进行反演计算；
[0159]
在本实施例中：
[0160]
步骤3中各参数先验区间见表4；
[0161]
表4参数先验分布特征
[0162][0163]
s3中的训练样本采样数为100，另采20检验样本；
[0164]
s4、s骤5中建立的各单一替代模型和集成学习替代模型检验输出与数值模型输出拟合情况见图3；
[0165]
s7中应用超启发式-同伦算法求解优化模型时，同伦参数的间隔取0.1，构造出包含10个同伦方程的方程系列，依据同伦方程系列对优化模型进行改写，可获得相应的10个优化模型。
[0166]
s8中对同伦优化模型系列求解时，前九个优化模型求解的最大进化代数为60，最后一个优化模型的大进化代数为120；得到反演识别结果及相对误差见表5；
[0167]
表5识别结果及相对误差
[0168]
[0169]
根据表5的结果，反演值与其真实值都非常接近，对误差均小于 5％；
[0170]
超启发式-同伦算法的循序渐进优化能力可以由不同污染源特征和含水层参数识别收敛曲线清晰地展示，见图4；在迭代优化的过程中，伴随着同伦参数的变化不同变量的识别值逐渐地向真实值趋近，特别是污染源的纵向坐标和横向坐标，以及污染物泄漏量。初始的大范围搜索问题被转化为若干个局部搜索的阶段，每一个阶段的最优值是容易求解的；尽管会出现陷于早熟收敛的不理想求解阶段，基于同伦算法的渐进求解机制可以减小和分散风险，后续的求解阶段可以再次将收敛趋势带向全局最优值。
[0171]
再将识别值代入s2中建立的多相流数值模型，计算末时刻含水层中污染物的分布情况，并与实施例末时刻污染物的实际分布进行对比，见图5；污染羽的形状十分接近，应用识别结果建立的模拟模型，能够较为准确地计算当前或预报未来含水层中污染物的分布情况，从而为地下水污染修复方案设计、风险评估提供可靠依据。
[0172]
如果用传统的模拟-优化方法来解决实施例中的污染源反演识别问题，需要调用多相流模拟模型39600(超启发式法的初始种群数量为60，最大进化代数为660次)次；该例中的污染含水层多相流模拟模型在cpu为intel core i5 3.0ghz，内存为8gb的计算机上运行一次平均需要大约600秒，整个优化求解过程将花费约275天时间。而集成学习替代模型运行一次仅需要1.2秒，如果以集成学习替代模型代替模拟模型，优化求解过程将缩短至13.2个小时。
[0173]
由此可见，发明能有效解决参数识别及模型校正问题，既提升了反演效率，同时降低了反演误差。
[0174]
本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

技术特征：
1.一种基于超启发式-同伦算法的地下水有机污染源反演方法，其特征在于：具体方法包括如下步骤：s1：通过野外现场调查、动态监测手段，掌握污染场地的水文地质条件，获得地下水水质动态监测数据；s2：对场地水文地质条件进行概化处理，初步建立地下水有机污染多相流数值模拟模型，用以描述地下水中有机污染物的运移机理；s3：确定污染源反演识别问题中的待识别的污染源特征及含水层参数；根据模型中待识别变量的取值范围，随机采样若干组样本，逐一代入s2中所建立的多相流数值模型，得到每组样本对应的模型响应，形成由“模型输入-模型响应”样本对构成的训练样本集和检验样本集；s4：根据s3获得的输入-输出训练样本样集，采用不同机器学习方法建立多相流数值模型的单一替代模型；s5：应用群智能优化算法，根据s3获得的检验样本集优化确定s4建立的各单一替代模型的权重及核心参数，建立群智能集成学习替代模型；s6：建立污染源特征及含水层参数反演识别的非线性规划优化模型，并将s5建立的群智能集成学习替代模型作为约束条件之一嵌入到优化模型之中；s7：基于同伦理论，改写s6建立的优化模型，获得一系列反演同伦优化模型；s8：采用超启发式算法依次求解s7建立的同伦优化模型系列，最后一个优化模型的解即为地下水污染源的反演识别结果。2.根据权利要求1所述的基于超启发式-同伦算法的地下水有机污染源反演方法，其特征在于：所述s2中用utchem程序构建地下水有机污染多相流运移模型。3.根据权利要求1所述的基于超启发式-同伦算法的地下水有机污染源反演方法，其特征在于：所述s3中待识别的污染源特征包括：污染源的纵向坐标、污染源的横向坐标、污染物迁移转化时长、污染物泄漏量；待识别的含水层参数包括：孔隙度、渗透率、纵向水相弥散度、横向水相弥散度；随机采样通过拉丁超立方采样方法实现。4.根据权利要求1所述的一种基于超启发式-同伦算法的地下水有机污染源反演方法，特征在于：所述s4中的单一替代模型建模方法包括：克里金法、支持向量回归法、小波核极限学习机法；其中，克里金建模方法如下：克里金法(kriging)的回归方程可以表示为：式中，f(x)＝[f1(x),f2(x),k,f
k
(x)]
t
为已知的基函数，β＝(β1,β2,k,β
k
)
t
为基函数对应的回归参数，由训练样本估计获得；z(x)为局部偏差项，均值为0，方差为其协方差可表示为：r(u,v)为n维向量u和v的相关函数：
式中，α
i
为待定参数，u
i
和v
i
为u和v的第i个元素；对于给定的样本p＝[p1,p2,k,p
m
]
t
和相应的输出响应q＝[q1,q2,k,q
m
]
t
，任意输入向量x的预测输出响应y(x)可以写成：式中，r是x与样本p之间的相关向量：r(x)＝[r(x,p1),r(x,p2),k,r(x,p
m
)]
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)f为样本p的响应列向量：r为样本p各样本点之间的相关矩阵：r为样本p各样本点之间的相关矩阵：是由广义最小二乘法估计得到β的估计值：支持向量回归(svr)建模方法如下：svr是通过非线性映射函数把输入数据映射到高维空间进行线性回归，在训练拟合精度和预测精度之间进行权衡；对于给定的训练输入x＝[x1,x2,k,x
m
]
t
每个要素代表一个n元输入：x
i
＝(x
i,1
,x
i,2
,k,x
i,n
),i＝1,2km和输出y＝(y1,y2,k,y
m
)
t
；非线性回归方程可表示为：f(x)＝<w,φ(x)>+b
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)式中，w＝(w1,w2,k,w
n
)为权值向量，b为拟合误差，<w,x>表示w和φ(x)的内积；φ(x)为非线性映射函数，可将输入向量从输入空间映射到高维空间；该问题可表示为如下优化问题：最小化
约束条件式中，||w||2是w的范数，常数c及变量ξ
i
和为惩罚因子和松弛变量，ε为可接受偏差；构造lagrange函数：式中，l为拉格朗日算子，η
i
、α
i
和为不小于0的拉格朗日乘子；对于最优解，和均为0：均为0：均为0：均为0：将以上性质与式11结合，可将式10中的优化问题改写为对偶形式：最小化约束条件且由式16可知：因此，可得到回归函数：应用核函数代替内积<φ(x),φ(x')>，最为常用的为高斯函数：
式16可改写为：最小化约束条件且最后得到：拟合误差b可以利用karush-kuhn-tucker(kkt)条件计算；小波核极限学习机(wkelm)建模方法如下：对于训练样本(x
j
,t
j
),j＝1,k,n，极限学习机的输出y
j
可表示为：q(x
j
)＝[p(ω1x
j
+b1),p(ω2x
j
+b2),k,p(ω
l
x
j
+b
l
)]
t
ꢀꢀꢀ
(23)式中，p(
·
)为激励函数，输入节点通过权重向量ω
i
与第i个隐含神经元连接，i＝1,k,l，b
i
为第i个隐含神经元的阈值，p(ω
i
x
j
+b
i
)为第i个隐含神经元的输出函数，β
i
为第i个隐含神经元与输出神经元连接的权重向量；式22可表示为如下形式：qβ＝y
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(24)式中，β＝[β1,k,β
l
]
t
，y＝[y1,k,y
n
]
t
，q为elm的隐含层输出矩阵：如果具有l个隐含节点的极限学习机模型能够无偏的学习n个训练样本，那么存在下式：式中，t
j
表示目标值；式26可以简写为：qβ＝t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(27)最小二乘法解决式27β＝q
+
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(28)式中，q
+
为q的moore
–
penrose广义逆；q
+
由下式计算q
+
＝q
t
(qq
t
)-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(29)对于训练样本(x
j
,t
j
),j＝1,k,n，kelm的原始优化问题表达为
s.t.q(x
j
)
t
·
β＝t
j-ξ
j
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(30)式中，c代表正则化参数，能够平衡训练误差和算法复杂性，ξ
j
代表误差；优化问题可转化为拉格朗日对偶形式求解：式中，θ
i
代表拉格朗日算子；该对偶问题可运用karush-kuhn-tucker(kkt)条件求解：tucker(kkt)条件求解：tucker(kkt)条件求解：根据式32可计算权重向量β的最小二乘解：依据核函数理论，构建小波隐式映射函数k(x
i,
x
j
)代替随机映射函数q(x
j
)：k
elm
＝qq
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(34)k
elm(i,j)
＝q(x
i
)
t
·
q(x
j
)＝k(x
i
,x
j
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(35)训练后的kelm输出函数表达如下：式中，t＝[t1,k,t
n
]
t
；小波核函数的表达式为：式中α
w
,γ
1,w
,γ
2,w
为可调节参数。5.根据权利要求1所述的一种基于超启发式-同伦算法的地下水有机污染源反演方法，特征在于：所述s5中的集成学习替代模型是s4中的单一替代模型加权线性叠加构成，其输出的表达式如下：
式中，为集成替代模型的预测值，分别为克里金法模型，支持向量回归模型与小波核极限学习机模型的预测值，w1，w2，w3为权重值，权重之和为1；构造集成学习替代模型的一个关键环节是集成权重值和各单一替代模型核心参数的确定，替代模型的预测精度越高，其权重应越大；替代模型的精度越差，其权重应越小；采取以下方式作为集成权重和模型参数的确定方法；首先，建立如下优化模型，优化模型的决策变量为集成学习替代模型中的权重值和各单一替代模型的可变参数值，优化模型的最优解使得运用检验样本获得的集成学习替代模型输出与模拟模型输出之间的均方根误差最小；型输出与模拟模型输出之间的均方根误差最小；式中，y
i
(x
k
,p
i
)为第i个替代模型以参数向量p
i
和第k个检验样本输入x
k
为输入对应的输出响应，i＝1,2,k,n，k＝1,2,k,m，y
actual
(x
k
)为第k个检验样本输入x
k
对应的模拟模型输出响应，m为检验样本的个数；应用群智能优化算法求解上述优化模型，获得集成权重和各单一替代模型可变参数的最优值，构建群智能集成学习替代模型。6.根据权利要求1所述的一种基于超启发式-同伦算法的地下水有机污染源反演方法，特征在于：所述s6中的非线性优化模型以监测井的模拟计算浓度值与实测浓度值的绝对偏差的平方和达到最小为目标函数；以地下水污染源反演识别问题中的待识别变量为决策变量，包括：污染源的纵向坐标、污染源的横向坐标、污染物迁移转化时长、污染物泄漏量、孔隙度、渗透率、纵向水相弥散度、横向水相弥散度；以群智能集成学习智能替代模型为污染物迁移转化规律的等式约束，以污染源特征和含水层参数取值范围为不等式约束条件；其表达式如下：表达式如下：
c
k
为第k口监测井中污染物浓度的实际监测值，k＝1,2,k,n，为对应的模拟计算值；为污染物浓度的模拟计算值向量，即群智能集成学习机智能替代模型的输出，s为污染源相关变量取值向量，p为含水层参数取值向量；约束条件分别表示污染源纵向坐标(x
i
)、污染源横向坐标(y
i
)、污染物迁移转化时长(t
i
)、污染物泄漏量(v
i
)、孔隙度(θ)、渗透率(k)、纵向油相弥散度(d
oil,l
)、横向油相弥散度(d
oil,t
)变化范围的上界及下界。7.根据权利要求1所述的一种基于超启发式-同伦算法的地下水有机污染源反演方法，特征在于：所述步骤7中的同伦优化模型系列构建首先需要借助同伦算法思想，引入同伦参数λ，构造同伦函数h(s,p,λ),使得当λ＝0时，方程组h(s,p,λ)＝0的解为假设的地下水污染源及含水层参数，而当λ＝1时，方程组h(s,p,λ)＝0的解为待求的真实的地下水污染源及含水层参数；从任意假设的污染源信息及含水层参数值出发，通过路径跟踪，依次对一系列同伦方程组所对应的优化问题进行求解，逐步趋近并最终获得待求的真实的污染源信息及含水层参数值；同伦函数取线性同伦如下式：h(s,p,λ)＝λf(s,p)+(1-λ)g(s,p)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(40)其中，f(s,p)＝f(s,p)-c
obs
，g(s,p)＝f(s,p)-c0；式中，s为污染源特征相关变量取值向量，p为含水层参数取值向量；λ为同伦参数，取值范围[0,1]；f()表示多相流模拟模型或其替代模型；c
obs
为监测点实际监测浓度向量；c0表示将任意假设的地下水污染源特征代入模拟模型得到的监测点计算浓度向量；同伦函数h连续依赖于同伦参数λ，将λ的取值范围区间[0,1]进行划分λ0＝0＜λ1＜
…
＜λ
n
＝1，得到方程组系列：h(s,p,λ
i
)＝f(s,p)-(λ
i
·
c
obs
+(1-λ
i
)
·
c0)＝0,i＝1,2,k,n
ꢀꢀꢀ
(41)如果λ
i+1-λ
i
充分小，则相邻方程的解(s,p)
i+1
与(s,p)
i
是非常接近的；依据同伦算法思想，s6中的优化模型进行改写，可得到本文地下水污染源反演识别问题的同伦方程组系列所对应的优化模型系列，如式(42)所示；表示将任意假设的地下水污染源特征及含水层参数值代入模拟模型得到的第k口监测井中的污染物计算浓度；及含水层参数值代入模拟模型得到的第k口监测井中的污染物计算浓度；8.根据权利要求1所述的一种基于超启发式-同伦算法的地下水有机污染源反演方法，特征在于：所述s7中的超启发式算法采用随机选择作为高层策略，即每次求解从给定的集合中随机选择低层启发式算法，集合中包括遗传算法、粒子群优化算法和模拟退火算法；
运用超启发式算法(hha)依次对同伦优化模型系列进行求解，在求解过程，由于路径跟踪过程是逐渐演变的，因此前后相邻的两个优化问题的解也是十分接近的：式中，为第i个同伦优化模型的初始种群，(s
i,opt
,p
i,opt
)为应用超启发式算法求得的第i个同伦优化模型的最优解，可作为生成下一优化模型初始种群的依据：初始种群的依据：初始种群的依据：式中，v
s,low
，v
s,up
，v
p,low
，v
p,up
，表示超启发式算法中污染源特征与含水层参数对应的调整位移向量的上限和下限；应用第i个优化问题的解作为求解第i+1个优化问题初始值的生成依据，从而保证了寻优过程是逐渐演变的，避免因初始值选择距最优解较远时而产生的早熟收敛问题；最后一个优化模型的解即为地下水污染源的反演识别结果。

技术总结
本发明涉及一种基于超启发式-同伦算法的地下水有机污染源反演方法，该方法包括以下步骤：根据观测资料，建立地下水有机污染多相流运移数值模型，确定模型中待反演的污染源特征、含水层参数以及各变量的取值范围；准备训练样本集；训练多相流数值模型的单一机器学习替代模型；建立组合智能替代模型；建立污染源特征及含水层参数反演识别的非线性规划优化模型，将集成学习智能替代模型耦合其中；以超启发式-同伦算法求解优化模型，反演识别待求的污染源特征值及含水层参数值；本发明可以解决地下水有机污染源及多相流运移参数反演问题，为地下水中有机污染物时空分布的准确模拟预测提供技术保障。预测提供技术保障。预测提供技术保障。


技术研发人员：王宇 卞建民 孙晓庆
受保护的技术使用者：吉林大学
技术研发日：2021.12.24
技术公布日：2022/7/4