中国环境科学  2021,41(2):713~719 China  Environmental  Science 基于U-D分解卡尔曼滤波地下水污染源溯源辨识
贾顺卿1,2,卢文喜1,2*,李久辉1,2,白玉堃1,2(1.吉林大学地下水与资源环境教育部重点实验室,吉林长春 130012;2.吉林大学新能源与环境学院,吉林长春 130012)
摘要:采用基于U-D分解的卡尔曼滤波与非线性规划优化模型相结合,溯源辨识出地下水污染源的个数、位置与释放强度.基于一个假想例子,建立地下水污染质数值模拟模型,运用灵敏度分析筛选出对模型影响较大的参数作为模型中的随机变量.然后,应用基于U-D分解的卡尔曼滤波辨识出污染源的个数与位置.在此基础上建立辨识污染源释放强度的优化模型,应用克里格插值法建立地下水污染质运移数值模拟模型的替代模型,代替模拟模型,作为约束条件嵌入优化模型中,运用遗传算法求解优化模型辨识出地下水污染源源强.结果表明:采用基于U-D分解的卡尔曼滤波方法能够保证滤波的稳定性,有效识别出污染源的个数和位置;非线性规划优化模型,可以辨识出污染源释放强度.在优化模型的求解过程中,应用克里格方法建立模拟模型的替代模型嵌入优化模型,能在保证一定精度的情况下,大幅度减少计算负荷和计算时间.
关键词:地下水;污染源溯源辨识;U-D分解卡尔曼滤波;替代模型;遗传算法
中图分类号:X523      文献标识码:A 文章编号:1000-6923(2021)02-0713-07
Inversion identification of groundwater contamination source based on U-D factorization Kalman filter. JIA Shun-qing1,2, LU Wen-xi1,2*, LI Jiu-hui1,2, BAI Yu-kun1,2 (1.Key Laboratory of Groundwater Resources and Environment Ministry of Education, Jilin University, Changchun 130012, China;2.College of Environment and Resources, Jilin University, Changchun 130012, China). China Environmental Science, 2021,41(2):713~719
Abstract:This paper uses U-D factorization Kalman filtering and a nonlinear programming optimization model to identify the number, location, and release intensity of groundwater pollution sources. Based on a hypothetical example, a numerical simulation model of groundwater contamination was established, and the parameters having large influences on the model were selected as random variables in the model by using sensitivity analysis. Then, Kalman filtering based on U-D factorization was used to identify the number and location of pollution sources. On the basis of these processes, an optimization model for identifying the release intensity of contamination sources was established, and a Kriging interpolation method was used to establish an Surrogate model for the numerical simulation model of groundwater pollution transport, as an alternative of the simulation model, which was embedded in the optimization model as a constraint condition, and the genetic algorithm was applied to solve the optimization model. Finally, the source strength of ground
water pollution was identified. The results show that the Kalman filter method based on U-D factorization could ensure the stability of the filter and effectively identify the number and location of pollution sources; the nonlinear programming optimization model could identify the release intensity of pollution sources. In the process of solving the optimization model, a substitute model of the simulation model embedded the optimization model was built with the Kriging method, which could greatly reduce the calculation load and time under the condition of a certain accuracy.
Key words:groundwater;contamination source identification;Kalman filter with UD factorization;Surrogate model;genetic algorithm
地下水作为人类生存与社会发展必要的自然资源,近年来遭受的污染问题愈发严重.地下水污染难以监测,不易治理.地下水污染溯源识别是地下水污染修复与治理过程中的重要一环,能提供污染源空间分布、污染物释放历史等信息,为其治理提供参考依据.
卡尔曼滤波是求解地下水污染反问题的一种有效方法[1].卡尔曼滤波作为一种最优状态估计方法,可以应用于受随机干扰的动态系统[2-3]. Herrera[4], Schmidt等[5]应用卡尔曼滤波对地下水水质监测网进行时空优化;Dokou[6]尝试使用卡尔曼滤波创建一个最优的搜索策略,使用最少的水质样本来识别DNAPL来源[6],Jiang等[7],江思珉等[8-9],顾文龙等[10]应用基于卡尔曼滤波技术和模糊集理论识别污染羽的
方法.常规卡尔曼滤波可能出现矩阵协方差非对称或者负定的情况[11-12],崔尚进将卡尔曼滤波协方差矩阵进行U-D分解,提高了卡尔曼滤波的稳定
收稿日期:2020-06-10
基金项目:国家自然科学基金资助项目(41972252)
* 责任作者, 教授,***************
714 中国环境科学 41卷
性[13].Chen等[14], Xu等[15-16]使用集合卡尔曼滤波对二维确定性含水层中的污染源进行了溯源辨识[14-16];场地信息存在不确定性[17-18],白玉堃等[19]采用灵敏度分析筛选出灵敏度最高的参数作为随机变量进行污染源反演.
结合前人研究,本文运用U-D分解卡尔曼滤波识别污染源的个数与位置,在此基础上建立优化模型,采用克里格插值法建立替代模型并嵌入优化模型,运用遗传算法求解优化模型得到污染源源强.
1研究方法
1.1技术路线
图1  污染源反演技术路线
Fig.1 Flow chart of contamination source identification
首先应用U-D分解卡尔曼滤波反演污染源的个数与位置.在现场调查与动态监测的基础上,结合专家意见初步估计污染源潜在的个数、位置以及初始权重,构建地下水溶质运移模拟模型.考虑场地信息存在
不确定性,对场地参数进行灵敏度分析,筛选出灵敏度最高的参数作为随机变量,其余参数作为确定型变量.给出随机变量一个取值范围,由拉丁超立方抽样生成参数随机场,采用蒙特卡罗法将参数随机场输入模拟模型生成溶质浓度场库,结合各潜在污染源初始权重计算得到初始综合浓度场与误差协方差矩阵.对误差协方差矩阵进行U-D分解,结合采样点数据运用U-D分解卡尔曼滤波对污染浓
度场和潜在污染源权重进行更新,权重稳定后判断
污染源的个数与位置.
在识别出污染源个数与位置的基础上建立优
化模型,进行污染源释放强度反演.以污染物监测浓
度与模拟计算浓度拟合误差极小化为目标函数,污
染源释放强度为待求变量,替代模型作为优化模型
等式约束条件,同时考虑源强上下限等约束条件,建
立非线性规划优化模型,运用遗传算法求解优化模
water pollution型反演出污染源释放强度.为了减小求解优化模型
时调用模拟模型产生的运算负荷,建立溶质运移模
拟模型的克里格替代模型.
1.2灵敏度分析与浓度场库建立
灵敏度分析可反应模型输出对参数变化的敏
感程度,本文采用Morris全局灵敏度分析方法评价
参数变化对模型输出的影响,选出影响程度最大的
参数进行研究.Morris设计采用“一次只改变一个参
数”的抽样取值方法,轮流计算各参数的目标函数值,
从而得到各个参数全局灵敏度及参数间相关性和
非线性的定性描述[20-21]:
11
EE[(,...,,,,...,)()]/
i i i i i i k i
y x x x x x y x
−+
=+Δ−Δ (1)
式中:EE i为第i个参数变化引起的变化效应, {}
12
,,,k
X x x x
=…为模型()
y x的k个参数.
采用文献[20]中的修正化方法进行评价,用影响
值S i判断参数变化对模型输出值的影响
1
000
()
i i
i i
i
y y y x x
S
x
−−
=Δ=
Δ
(2)
多次抽样取均值
1
00
11
()
11
n n
i
i
i i
i
y y y
S S
n n
==
==
Δ
∑∑ (3) 把研究区用网格剖分,将筛选出的参数进行拉
丁超立方抽样得到n组随机变量参数场,只考虑一
个潜在污染源时,输入模拟模型得到对应的n组溶
质浓度场,m个潜在污染源共m n×组.计算每个网
格上参数场取不同值时浓度值的均值,就得到了该
污染源的均值浓度场.在同一随机变量参数场下,m
个污染源的溶质浓度场按照初始权重进行加权叠
加得到叠加浓度场.将m个污染源的均值浓度场按
初始权重进行加权叠加得到初始综合浓度场.
1.3U-D分解卡尔曼滤波
常规卡尔曼滤波存在数值稳定性较差,计算复
2期 贾顺卿等:基于U -D 分解卡尔曼滤波地下水污染源溯源辨识 715
杂等问题,针对这一点可引入U -D 分解改进.
常规卡尔曼滤波更新方程[8-
9]:
1
[][][]T T r C C Z C I −−−+−−+−
=×+=+−=−K P H HP H K H P KH P  (4)
式中:K 为卡尔曼增益矩阵,P 为n 阶误差协方差矩阵,C 为n 维状态向量,表示浓度估计值;-、+号表示先
验估计与后验估计;Z 为采样值,r 是采样误差的方差,H 是1×n 阶量测矩阵,采样点处矩阵元素为1,其余位置为0,[00,1,00]H =K K ;I 是n 阶单位向量.
P 矩阵更新时矩阵相减易导致矩阵稳定性降低,甚至引发滤波发散,将P 矩阵分解为P =UDU T 的形式,方程表示如下:
1
[]g a C C K Z A B
−−+−−+−+=××=+−=×=K U HC U U D  (5)
其中
1
T
T
T
T
a g r
g f f
H
a g g ABA
−−−
−−=××+=×=×−××=H U D U D  (6) 式中:U 为单位上三角矩阵;D 为正定对角矩阵. 1.4  权重更新
每个潜在污染源均值浓度场对应的污染羽为单个污染羽,综合浓度场对应的污染羽为综合污染羽.确定污染源的个数与位置需要通过污染源权重判断,将综合污染羽与单个污染羽标准化,用模糊集表示(所有浓度值除以最大浓度值),每个模糊集合的元素隶属度都大于等于给定的αi ,本文取值为0.2、0.4、0.6、0.8.
将以模糊集形式表示的综合污染羽与各单个污染羽进行比较,记录二者的公共面积S i ,计算出全局相
似度并用各污染源全局相似度除以各污染源全局相似度中的最大者,即可算得各污染源的权重.全局相似度为
[7-9]
:
1n
i i i g S α==∑ (7) 式中:g 为全局相似度;αi 为模糊集标准值,
1,2,i =…;S i 为综合污染羽与单个污染羽在同一αi
下的交叉面积.
图2  污染羽对比示意
Fig.2  Comparison of pollution plumes
1.5  优化模型建立
以污染物监测浓度与模拟计算浓度拟合误差极小化为目标函数,污染源释放强度为待求变量,替代模型作为优化模型等式约束条件,同时考虑源强上下限等约束条件,建立非线性规划优化模型
目标函数:
21ˆmin ()(-)n
k k k Z q c c
==∑ (8) 约束条件:
min max q q q ≤≤ (9)
()k c f q = (10)
式中, q 为污染源释放强度,n 为采样点个数,c k 为替
代模型输出的模拟值,ˆk c 为实测值,f 为替代模型. 1.6  替代模型
应用遗传算法求解优化模型的过程中会大量调用溶质运移模拟模型,使用克里格方法建立溶质运移模拟模型的替代模型能有效减少运算负荷.替代模型回归方程为
1ˆ()()()()()n
T i i y
x f x Z x f x Z x ββ=+=+∑ (11) 式中:ˆ()y
x 为()y x 的估计值,()y x 为要求解的污染物浓度值;()T f x β是线性回归部分,()Z x 是随机部分;12()[(),()()]n f x f x f x f x =K 为模型的基函数,12[,]n ββββ=K 为基函数系数,可应用训练数据求出,()Z x 满足以下条件:
22
(())0
(())cov[(),()](,)i j i j E Z x D Z x Z x Z x R x x σ
σ⎧=⎪⎪=⎨⎪=⎪⎩ (12) (,)i j R x x 为任意两点,i j x x 之间的关联函数,本文采用高斯型关联函数表示:
716 中  国  环  境  科  学 41卷
1
(,)exp(||)n
i
j i j k k k R x x x x θ=−−∑ (13)
式中:k θ为待定参数,i
k x 为第i 个样本第k 维坐标.
克里格模型中,点x 的响应值()y x 的估计值为
1ˆ()()()T T y x f r x R y f ββ−=+− (14) 式中:()r x 为点x 与n 个采样点12(,,,)n x x x K 的相关向量, 12()[(,),(,),,(,)]n r x R x x R x x R x x =K ;y 为n 个采样点对应的响应值,为1n ×阶向量;β为待定参数,可以通过最优线性无偏估计求得:
11()T T T f R f f R y β−−= (15)
相关矩阵R 可表示为:
1111(,)(,)(,)(,)n n n n R x x R x x R R x x R x x ⎛⎞⎜⎟
=⎜⎟⎜⎟
⎝⎠L M O M L  (16) 方差2
σ可表示为:
21()()/T y f R y f n σββ−=−− (17)
替代模型通过求解上面的非线性无约束优化问题来实现,待定参数k θ可通过一个无约束优化问题求得:
2min{ln ln ||}n R σ+ (18)
2  假想算例 2.1  问题概述
图3  污染场地
Fig.3  The plan of contaminated site
假定研究区大小为1000m ×1000m,潜水含水层厚度为20m,介质为粗砂,含水层非均质各向同性.忽略非饱和带的影响.南北边界为已知水头边界,以含
水层底板为基准面,北边界水位为16m,南边界水位为14m;东西边界为隔水边界.模型模拟时间长度为500d,分为5个时段,每个时段为100d.研究区接受降雨入渗补给,5个时段降水量分别为100,300, 100,50, 100mm.忽略研究区蒸发蒸腾作用的影响.结合场地
信息与专家经验,估计潜在污染源个数为4个(图3),并给出潜在污染源的初始权重(权重表示潜在污染源为真实污染源的可能性大小,取值0~1之间),Ⅰ、Ⅱ、Ⅲ、Ⅳ号污染源分别为0.6、0.7、0.7、0.6.设定Ⅲ号为真实污染源,污染源泄漏流量为2000mg/d. 2.2  问题分析
2.2.1  位置个数反演  将研究区剖分成20×20的
网格,应用GMS 软件的MODFLOW 和MT3DMS 模块建立溶质运移模拟模型,考虑场地信息的不确定
性,对模型中的参数进行Morris 全局灵敏度分析,筛选出渗透系数K 、纵向弥散系数L d 、给水度S y 、降雨入渗补给系数α灵敏度最高的参数作为随机变量.共有15种组合方式,具体方法参照文献[21].
表1  参数组合方式
Table 1  Co mbinatio n -patterns of parameters
组号 1 2 3 4 5 参数 K L d
S y  α K,L d
组号 6 7 8 9 10
参数 K,S y  K,α L d ,S y
L d ,α S y ,α
组号 11 12 13 14 15
参数
K,L d ,S y  K,L d ,α K,S y ,α L d ,S y ,α K,L d ,S y ,α
0.20.40.60.811.21.41234567 8 9 10 11 12 131415
组号
灵敏度
图4  参数灵敏度分析
Fig.4  Sensitivity analysis of parameters
从图4可以看出,渗透系数灵敏度最高,本研究将渗透系数作为随机变量,其余参数作为确定性变量,应用拉丁超立方抽样生成110组渗透系数随机场并进行相关性排序,场地存在4个潜在污染源,输
2期 贾顺卿等:基于U -D 分解卡尔曼滤波地下水污染源溯源辨识 717
入模型得到4×110组溶质浓度场,计算得出初始综合浓度场和误差协方差矩阵,Ⅲ号污染源的均值浓度场作为真实浓度场.
x (m)
y (m )
图5  真实污染羽 Fig.5  True po llutio n plume
将采样点数据带入U -D 分解卡尔曼滤波更新方程对综合浓度场和污染源权重进行更新,图为浓度场和权重更新结果.
图5为污染源真实污染羽,图6为初始综合污染羽.从图7可以看出,经过4次采样更新后,综合污染羽逐
渐收敛于Ⅲ号污染源,与真实污染羽相似.4个污染源权重从0.6,0.7,0.7,0.6变为0,0.2,1,0.05,可以判断出Ⅲ号污染源为真实污染源,更新完成后污染羽中浓度最高处就是污染源位置.
y (m )
x (m)
图6  初始综合污染羽
Fig.6  Initial composite pollution plume
x (m)
y (m )
x (m)
y (m )
a.1次更新
b. 2次更新
y (m )
x (m)
x (m)
y (m )
c. 3次更新
d. 4次更新