128 水动力学研究与进展 A辑2009年第2期 distribution,the Adaptive Metropolis algorithm was used to construct the Markov Chains of unknown parameters.And the converged samples were used to estimate the unknown parameters of source term.The results of numerical experiments show that the method has many virtues,such as high accuracy,quick convergent speed and easy to program and implement with computer. Key words:convection-diffusion equation;source term;inverse problem;Markov Chain Monte Carlo method 题进行了数值研究。综上所述,由于对流-扩散方程 1引言 反问题的不适定性,所以它的求解一般要采用特殊 方法,如Tikhonov正则化方法、变分伴随方法和遗 传算法等等。本文在贝叶斯理论的基础上,提出采 对流-扩散方程是描述粘性流体运动的非线性 用马尔科夫链蒙特卡罗(Markov Chain Monte Carlo, Burgers方程的线性化模型,它可以刻画许多自然现 简称MCMC)方法222来识别对流扩散方程中多个 象,如:水体和大气中污染物的输移、扩散和降解, 点源中的未知参数。 海水盐度和温度的扩散,流体流动与传热和电化学 结合利用贝叶斯方法和MCMC算法求解反问 反应等。研究对流扩散模型具有重要的理论价值和 题,具有以下优点:1)能方便地将各种先验信息和 实际意义,它已经广泛应用于环境科学、能源开发、 误差信息高效地融合到问题求解过程中,减小问题 流体力学和电子科学等领域。总的来说,目前关于 的不确定性:2)和确定性算法不同,反问题的不适 对流-扩散方程的研究大致可以分为两个方面。一方 定性不再是MCMC算法要考虑的问题,且计算获得 面是在给定初边值条件下,通过不同的数值计算方。 的是全局最可能解,而通常的最优化算法可能陷入 法求解对流扩散方程,以模拟研究对象(例如:温 目标函数局部极小值:3)能对定义在高维空间且无 度、盐度和污染物等)在时间和空间上的发展演化, 明确数学表达式的概率分布密度函数进行数值计 这类问题可以统称为正问题。迄今为止已经有很多 算,而确定性方法无法解决此类问题;4)MCMC算 成熟方法求解对流-扩散方程,如有限差分方法 法通过构造Markov链来进行随机模拟,是一种动态 (FDM12)、有限体积方法(FVM56和有限元方法 Monte Carlo方法,计算速度高于一般的Monte Carlo (FEM07,89等。 方法和模拟退火算法,而且计算复杂度不依赖于计 另外一方面是关于对流-扩散方程反问题的研 算空间的维数。 究,即通过所研究对象的观测资料来估计和识别方 程中的参数、源项、边界和初始条件等。从某种意 义上讲,反问题的求解是对流扩散模型研究中一个 2反问题模型 更重要的问题,因为它的正确与否直接影响模型的 可靠性。由于偏微分方程反问题固有的非线性和不 适定性1,对流-扩散方程反问题的求解会存在巨 不失一般性,用对流扩散方程来模拟污染物在 大困难,通常的方法常常导致求解失败。近年来,国 河道中的扩散,考虑对流扩散方程的初边值问题 内外学者关于对流-扩散反问题开展了广泛研究。 92川, 公式如下: Andreas Kirsch对一维扩散方程逆过程反问题进行 了稳定性分析,并给出了误差估计公式。 acac 8 Yildiz1、刘继军等3-l61对相关问题采用Tikhonov +u- t --kC+s,(x-x). 正则化方法进行了深入研究。闵涛等以函数逼近 和Tikhonov正则化为基础,利用算子识别摄动法和 (x,t)∈(0,L)x(0,T)(1) 线性化技术,建立了河流水质纵向弥散系数反问题 C(0,t)=0,C(L,t)=0, tE(0,T) 的迭代算法,并进行了数值试验。闵涛等18利用有 C(x,0)=0, x∈(0,L) 限元法求解了二维稳态对流扩散方程,并利用迭代 法对二维稳态对流扩散方程参数反演进行了研究。 其中C为污染物的浓度,4为流速,E为扩散系数, 闵涛等例利用遗传算法就对流扩散方程的源项识 k为污染物的降解率,L表示河道长度。δ是狄拉 别反问题进行了研究。潘军峰等20对一维对流-扩 散方程的反问题利用Tikhonov正则化方法进行了研 克函数,x和s,(i=1,2,…g)分别表示多个点污 究。吴自库等2结合利用伴随同化方法和处理数学 染源的位置和排放强度。假定C(x,t)在1=T时的 物理反问题的技巧就对流-扩散方程逆过程的反问 分布已知,那么源项识别反问题就是根据这些已知 万方数据128 水动力学研究与进展 A辑2009年第2期 distribution,the AdaptiVe Me仃opolis algorithm w勰吣ed t0 cons仇Jct the MarkoV Chains of幽own p啪mete娼.And the converged s锄ples were llsed t0 estimate tlle unl【nown parameters of source tc衄.The rcsults of numerical experimems show that me method has m锄y Vin∞s,鲫ch鹊high accuracy,quick conVergem speed锄d e觞y t0 program加d iInpl锄em with computer. Key words:conVection—dim塔ion equation;source telm;inVerse problem;MarkoV Chain Monte Carlo method 1引言 对流.扩散方程是描述粘性流体运动的非线性 Burgers方程的线性化模型,它可以刻画许多自然现 象,如:水体和大气中污染物的输移、扩散和降解, 海水盐度和温度的扩散,流体流动与传热和电化学 反应等。研究对流.扩散模型具有重要的理论价值和 实际意义,它已经广泛应用于环境科学、能源开发、 流体力学和电子科学等领域。总的来说,目前关于 对流.扩散方程的研究大致可以分为两个方面。一方 面是在给定初边值条件下,通过不同的数值计算方 法求解对流一扩散方程,以模拟研究对象(例如:温 度、盐度和污染物等)在时间和空间上的发展演化, 这类问题可以统称为正问题。迄今为止已经有很多 成熟方法求解对流.扩散方程,如有限差分方法 (FDM)oⅥ’3J、有限体积方法(FvM)‘4'5,6J和有限元方法 (FEM)【7’8'91等。 另外一方面是关于对流.扩散方程反问题的研 究,即通过所研究对象的观测资料来估计和识别方 程中的参数、源项、边界和初始条件等。从某种意 义上讲,反问题的求解是对流一扩散模型研究中一个 更重要的问题,因为它的正确与否直接影响模型的 可靠性。由于偏微分方程反问题固有的非线性和不 适定性【l 01,对流.扩散方程反问题的求解会存在巨 大困难,通常的方法常常导致求解失败。近年来,国 内外学者关于对流.扩散反问题开展了广泛研究。 Andreas Kirsch对一维扩散方程逆过程反问题进行 了稳定性分析,并给出了误差估计公式¨¨。 Yildiz【12J、刘继军等【13-16J对相关问题采用Til(1lonov 正则化方法进行了深入研究。闵涛等【l 7J以函数逼近 和Til(1lonov正则化为基础,利用算子识别摄动法和 线性化技术,建立了河流水质纵向弥散系数反问题 的迭代算法,并进行了数值试验。闵涛等【l8j利用有 限元法求解了二维稳态对流.扩散方程,并利用迭代 法对二维稳态对流.扩散方程参数反演进行了研究。 闵涛等【l圳利用遗传算法就对流.扩散方程的源项识 别反问题进行了研究。潘军峰等【20J对一维对流.扩 散方程的反问题利用Ti姓onov正则化方法进行了研 究。吴自库掣21j结合利用伴随同化方法和处理数学 物理反问题的技巧就对流.扩散方程逆过程的反问 题进行了数值研究。综上所述,由于对流.扩散方程 反问题的不适定性,所以它的求解一般要采用特殊 方法,如Til【llonov正则化方法、变分伴随方法和遗 传算法等等。本文在贝叶斯理论的基础上,提出采 用马尔科夫链蒙特卡罗(Markov Chain Monte Carlo, 简称McMc)方法【2理到来识别对流一扩散方程中多个 点源中的未知参数。 结合利用贝叶斯方法和MCMC算法求解反问 题,具有以下优点:1)能方便地将各种先验信息和 误差信息高效地融合到问题求解过程中,减小问题 的不确定性;2)和确定性算法不同,反问题的不适 定性不再是MCMC算法要考虑的问题,且计算获得 的是全局最可能解,而通常的最优化算法可能陷入 目标函数局部极小值;3)能对定义在高维空间且无 明确数学表达式的概率分布密度函数进行数值计 算,而确定性方法无法解决此类问题;4)MCMC算 法通过构造Markov链来进行随机模拟,是一种动态 Monte Carlo方法,计算速度高于一般的Monte Carlo 方法和模拟退火算法,而且计算复杂度不依赖于计 算空间的维数。 2反问题模型 不失一般性,用对流一扩散方程来模拟污染物在 河道中的扩散,考虑对流一扩散方程的初边值问题 【19,211,公式如下: 等+“篆=E等一圮+和¨,, (x,f)∈(0,三)×(0,丁) (1) C(0,f)=O,C(厶f)=0, f∈(0,丁) C(x,0)=O, x∈(O,三) 其中C为污染物的浓度,甜为流速,E为扩散系数, 七为污染物的降解率,£表示河道长度。6是狄拉 克函数,‘和墨,(f=l,2,…g)分别表示多个点污 染源的位置和排放强度。假定CO,f)在f=丁时的 分布已知,那么源项识别反问题就是根据这些已知 万方数据