第41卷第5期四川大学学报(工程科学版)Vo.l41No.52009年9月JOURNALOFSICHUANUNIVERSITY(ENGINEERINGSCIENCEEDITION)Sept.2009
文章编号:1009-3087(2009)05-0030-06水动力–水质耦合模型污染源识别的贝叶斯方法
朱 嵩1,刘国华1*,王立忠1,毛根海1,程伟平1,黄跃飞2(1.浙江大学建筑工程学院,浙江杭州310027;2.清华大学水利水电工程系,北京100084)
摘 要:环境水力学系统存在诸多不确定性,如测量数据的不确定性等,这导致水体中污染源识别这一类反问题具有不适定性,尤其表现为反演结果的非唯一性。经典的正则化方法和最优化方法由于只能获得参数的/点估计0,因而在求解不确定性较强的问题时存在较大的困难。此外水质模型和流场控制方程(Navier-Stokes方程)耦合,使得正问题的解具有较强的非线性特征。为解决上述问题,针对水动力-水质耦合模型,建立了基于贝叶斯推理的污染物点源识别的数学模型,通过马尔科夫链蒙特卡罗(MarkovchainMonteCarlo,MCMC)后验抽样获得了污染源位置和强度的后验概率分布和估计量,较好地处理了模型的不确定性和非线性。算例结果表明,结合MCMC抽样的贝叶斯推理方法能很好地描述及求解水动力-水质耦合场条件下的污染源识别反问题。关键词:环境水力学;反问题;贝叶斯推理;污染源识别中图分类号:TV13;X192文献标识码:A
ABayesianApproachfortheIdentificationofPollutionSourceinWaterQualityModelCoupledwithHydrodynamics
ZHUSong1,LIUGuo-hua1,WANGLi-zhong1,MAOGen-hai1,CHENGWei-ping1,HUANGYue-fei2(1.CollegeofCivilEng.andArchitecture,ZhejiangUniv.,Hangzhou310027,China;2.Dept.ofHydraulicEng.,TsinghuaUniv.,Beijing100084,China)
Abstract:Muchuncertaintyliesintheenvironmentalhydraulicssystem,suchastheuncertaintyofthemeasurementdata,havecausedthepollutionsourceidentificationil-lposed,especiallythenon-unique.Inordertosettletheproblem,forthehydrodynamics-waterqualitycoupledmode,lapollutionpointsourceidentificationmodelwaspro-posedbasedonBayesianinference.MarkovchainMonteCarlosamplingmethodwasusedtogettheposteriorproba-bilitydistributionofthesource.spositionandintensity,thussolvingtheuncertaintyandthenonlinearitywel.lComputationresultshowedthattheBayesianinferencewithMCMCsamplingcandescribeandsolvethepollutionsourceidentificationinverseproblemforthehydrodynamics-waterqualitycoupledmodelbetter.Keywords:environmentalhydraulics;inverseproblem;Bayesianinference;pollutionsourceidentification
环境水力学是当前水力学及河流动力学研究的收稿日期:2008-05-03基金项目:浙江省自然科学基金资助项目(Y506138);国家自然科学基金项目(50609024);973课题(2005CB724202)作者简介:朱 嵩(1981-),男,博士后.研究方向:环境水力学反问题.E-mai:lmigao@zju.edu.cn
一个主要方向,其主要研究内容是水体中的污染物迁移、扩散和转化机理以及相应的求解方法[1]。环境水力学研究自上个世纪70年代末以来得到了迅猛的发展。目前关于污染物在水体中的迁移转化规律及相关预测数学模型的研究已经较为成熟,如对流扩散方程的各种数值解法等。正问题(预测问题)研究的成熟为环境水力学反问题研究奠定了基础。水体中污染源识别反问题就是环境水力学反问题中一类重要的问题,其解决对于水质管理、水环境控制等具有较大的实际意义。反演的主要困难来源于测量数据有限且带有噪声,这样使得污染源识别是一个不适定的问题。针对不适定问题,目前求解采用的主要方法有脉冲谱法、离散优化法、控制论方法和正则化方法等[2-3]。水体中污染源识别反问题是指根据实测的水环境观测数据(如污染物浓度等)来反演污染源的强度、位置等信息。闵涛采用优化算法(遗传算法)研究了对流-扩散方程的源项识别反问题[3]。Vo-lkanAk®elik等采用Tikhonov正则化和总变差正则化的方法求解了对流扩散输运系统的源项反演问题[4]。此外,LeevanLing等从数学理论的角度研究了不考虑测量误差的二维扩散方程的源项位置识别问题,获得了关于反演唯一性的结论[5]。由于实际水体浓度观测总是带有误差,因而研究带有测量误差条件下的污染物源识别是需要解决的重要课题。测量噪声的考虑将使得污染源反演可能存在多解,然而正则化方法和遗传算法等优化算法只能获得反演变量的/点估计0,即只能获得最优解,因而这类方法在处理考虑测量噪声的反问题中较为困难[6-7]。与此相反,贝叶斯方法在处理环境水力学反问题上有较大的优点。贝叶斯方法把不适定的反问题认为是在一个扩展随机空间上的适定问题。它不是仅仅计算/点估计0,而是计算反演变量的后验概率分布,因而它能很好解决由观测数据噪声带来的非唯一解问题。Snodgrass和Kitanidis较早地将贝叶斯方法引入到污染源的识别的研究中[8]。此后,JingboWang等研究了多孔介质流动中的污染物瞬时源的识别问题[9]。针对由Navier-Stokes方程和含源对流扩散方程耦合的水动力-水质模型的点源的位置与强度的联合识别问题,建立了基于贝叶斯推理的反演数学模型,采用马尔科夫链蒙特卡罗(MarkovchainMonteCarlo,MCMC)抽样方法对后验空间进行抽样获得了污染源的后验概率分布规律以及相应的估计。1 正问题数学模型1.1 控制方程以二维水动力-水质耦合模型作为正问题数学模型。流场采用稳态Navier-Stokes方程描述,污染物浓度场采用保守物质的非稳态含源对流扩散方程描述,水动力-水质耦合模型方程如下:Ñ#u=0(1)Q(u#Ñ)u-GÑ2u+Ñp=0(2)
5c5t+u#Ñc=Ñ#(DÑc)+R(3)
其中,u为流速矢量,Q为流体的密度,G为流体的动力粘度,p为压强,c为污染物的浓度,D为污染物综合扩散系数,R为污染物源项。当仅有点源存在时,R
=Eq-1i=0MiD(x-xi,y-yi),其中,q为污染源的个数,Mi为点源的强度,(xi,yi)为点源的位置。此外,水的密度及动力粘度均取20e时的值。1.2 求解方法有限单元法(finiteelementmethod,FEM)在处理不可压Navier-Stokes方程和对流扩散方程已经较为成熟,文中使用FEM作为水动力-水质耦合模型的求解方法,限于篇幅,此部分将不做详细讨论。由于在污染物点源附近浓度梯度较大,在点源位置附近进行网格自动加密,以保证计算精度。
2 污染物源项识别反问题由上节可知,污染物的强度、位置等信息是控制污染物浓度分布的重要参数。但往往这些信息是模糊的或未知的,因而需要通过获得的实测的水动力/水质观测数据来对污染源的一些信息进行反演计算。2.1 环境水力学反问题的不适定性正问题和反问题的一个根本性的区别是反问题是不适定的(il-lposed),而正问题往往是适定的(wel-lposed)。关于适定和不适定的概念是Had-amard为了描述数学物理问题与定解条件的合理搭配,于20世纪初引入的[10]。根据Hadamard,一个物理问题的数学模型必须满足如下三条性质,可认为是适定的。1)该问题存在一个解(存在性);2)该问题至多存在一个解(唯一性);3)解连续依赖数据(稳定性)。一般地,对于一个实际环境水力学问题,若定解条件提法得当,其正问题总是适定的。而于此相反,环境水力学反问题一般而言总是不适定的,其主要原因是环境水力学系统的不确定性[11]。
31 第5期朱 嵩,等:水动力–水质耦合模型污染源识别的贝叶斯方法 对于环境流动反问题,不确定性主要来源于测量数据、数学模型、模型参数、定解条件以及污染物源项。流体系统本身的非线性和不确定性给环境流动反问题求解带来了很大的困难。2.2 结合蒙特卡罗方法的贝叶斯推理反演方法如何处理由不确定性给环境水力学反问题带来的不适定性(尤其是反演结果的不唯一性)是环境水力学反问题研究首先需要解决的课题。贝叶斯推理属于统计反演方法,由于它建立在贝叶斯统计学的基础上,因而能以概率语言来描述和解决工程反问题。2.2.1 贝叶斯定理贝叶斯推理的基础是贝叶斯定理,它可以表述如下,
p(H|y)=p(H)p(y|H)p(y)Wp(H)p(y|H)(4)其中,H是模型参数,y是观测数据,p(H)是参数的先验概率密度函数,p(y|H)是似然函数,p(H|y)是参数的后验概率密度函数。在贝叶斯推理中,p(H)表示在未获得测量数据之前,对模型参数分布的认识,主要来源于以往数据、经验和主观判断等。p(y|H)代表模型参数拟合测量数据的程度,越大表示拟和效果较好,反之则差。p(H|y)表示了获得测量数据以后模型参数的分布规律,即为在统计反演意义下的反问题的解。总体来说,贝叶斯推理可以分为如下3个步骤:1)基于未知参数的所有先验信息,确定一个先验概率密度函数;2)找到能够反映模型参数和测量数据之间关系的一个似然函数;3)对后验概率密度函数抽样,进而获得参数的估计值。一般情况下,对于参数Hi,如果仅知道它分布的区间[ai,bi],那么先验概率密度函数可以写成,