对流_扩散方程源项识别反问题的MCMC方法_曹小群

DOI：10.3969/j.issn.1000-4874.2010.02.001水动力学研究与进展A辑2009年第2期 1281 引言对流-扩散方程是描述粘性流体运动的非线性Burgers方程的线性化模型,它可以刻画许多自然现象,如：水体和大气中污染物的输移、扩散和降解，海水盐度和温度的扩散，流体流动与传热和电化学反应等。

研究对流-扩散模型具有重要的理论价值和实际意义,它已经广泛应用于环境科学、能源开发、流体力学和电子科学等领域。

总的来说，目前关于对流-扩散方程的研究大致可以分为两个方面。

一方面是在给定初边值条件下，通过不同的数值计算方法求解对流-扩散方程，以模拟研究对象(例如：温度、盐度和污染物等)在时间和空间上的发展演化，这类问题可以统称为正问题。

迄今为止已经有很多成熟方法求解对流-扩散方程,如有限差分方法(FDM)[1,2,3]、有限体积方法(FVM)[4,5,6]和有限元方法(FEM)[7,8,9]等。

另外一方面是关于对流-扩散方程反问题的研究，即通过所研究对象的观测资料来估计和识别方程中的参数、源项、边界和初始条件等。

从某种意义上讲，反问题的求解是对流-扩散模型研究中一个更重要的问题,因为它的正确与否直接影响模型的可靠性。

由于偏微分方程反问题固有的非线性和不适定性[10], 对流-扩散方程反问题的求解会存在巨大困难，通常的方法常常导致求解失败。

近年来, 国内外学者关于对流-扩散反问题开展了广泛研究。

Andreas Kirsch对一维扩散方程逆过程反问题进行了稳定性分析，并给出了误差估计公式[11]。

Yildiz[12]、刘继军等[13-16]对相关问题采用Tikhonov 正则化方法进行了深入研究。

闵涛等[17]以函数逼近和Tikhonov正则化为基础，利用算子识别摄动法和线性化技术，建立了河流水质纵向弥散系数反问题的迭代算法，并进行了数值试验。

闵涛等[18]利用有限元法求解了二维稳态对流-扩散方程,并利用迭代法对二维稳态对流-扩散方程参数反演进行了研究。

闵涛等[19]利用遗传算法就对流-扩散方程的源项识别反问题进行了研究。

潘军峰等[20]对一维对流-扩散方程的反问题利用Tikhonov正则化方法进行了研究。

吴自库等[21]结合利用伴随同化方法和处理数学物理反问题的技巧就对流-扩散方程逆过程的反问题进行了数值研究。

综上所述，由于对流-扩散方程反问题的不适定性，所以它的求解一般要采用特殊方法，如Tikhonov正则化方法、变分伴随方法和遗传算法等等。

本文在贝叶斯理论的基础上，提出采用马尔科夫链蒙特卡罗(Markov Chain Monte Carlo,简称MCMC)方法[22,23]来识别对流-扩散方程中多个点源中的未知参数。

结合利用贝叶斯方法和MCMC算法求解反问题，具有以下优点：1) 能方便地将各种先验信息和误差信息高效地融合到问题求解过程中，减小问题的不确定性；2) 和确定性算法不同，反问题的不适定性不再是MCMC算法要考虑的问题，且计算获得的是全局最可能解，而通常的最优化算法可能陷入目标函数局部极小值；3) 能对定义在高维空间且无明确数学表达式的概率分布密度函数进行数值计算，而确定性方法无法解决此类问题；4) MCMC算法通过构造Markov链来进行随机模拟，是一种动态Monte Carlo方法，计算速度高于一般的Monte Carlo 方法和模拟退火算法，而且计算复杂度不依赖于计算空间的维数。

2 反问题模型不失一般性，用对流-扩散方程来模拟污染物在河道中的扩散，考虑对流-扩散方程的初边值问题[19,21]，公式如下：221(),(,)(0,)(0,) (0,)0,(,)0,(0,)(,0)0,(0,)qi iiC C Cu E kC s x xt x xx t L TC t C L t t TC x x Lδ=⎧∂∂∂+=−+−⎪∂∂∂⎪⎪∈×⎨⎪==∈⎪=∈⎪⎩∑(1)其中C为污染物的浓度，u为流速，E为扩散系数，k为污染物的降解率，L表示河道长度。

δ是狄拉克函数，ix和is，(1,2,)i q= 分别表示多个点污染源的位置和排放强度。

假定(,)C x t在t T=时的分布已知，那么源项识别反问题就是根据这些已知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曹小群，等：对流-扩散方程源项识别反问题的MCMC 方法129的浓度分布观测来确定源项1()qiii s x x =δ−∑，即确定多个点污染源的位置和排放强度。

反问题在求解过程中通常需要将未知参数向量的估计值映射成观测空间的值，这就需要获得正问题的解算子。

本文采用的是Fourier 方法来求解对流-扩散方程，闵涛等[18]和吴自库等[21]分别采用相同的方法对系统(1)进行了求解。

首先引入下面的函数变换2(/2/4)(,)(,)e ux E u t E C x t V x t −= (2)将其代入(1)中，则(1)可转化为等价的定解问题：22(/4/2)21e()(0,)0(,)0(,0)0qu t E ux E i i i V V E kV s x x t x V t V L t V x −=⎧∂∂=−+δ−⎪∂∂⎪⎪=⎨⎪=⎪=⎪⎩∑ (3)由于方程组(3)的特征值和特征函数分别为2(π/)n n L λ= 和()sin(π/)n x n x L ϕ=，(1,2,)n =因此可利用特征函数展开法求解方程(3)。

令1(,)()()n n n V x t T t x ϕ∞==∑ (4)将(3)中第一式的右端源项也利用特征函数进行展开：2(/(4)/(2))11e ()()q u t E ux E i i n n i n s x xf x ϕ∞−==δ−=∑∑ (5)其中2(/(4)/(2))012e ()()d qL u t E ux E n i i n i f s x x x x L ϕ−==δ−=∑∫2(/(2))412e (e ())i u q t ux E E i n i i s x L ϕ−=∑将(4)和(5)式代入(3)式中，则可求得()n T t 的解为2[]/(4)2()[e e ]/(4)n t k u t E n n B T t u E kλλ−+=−++ (6) 其中/(2)12e ()i qux E i n i i B s x L ϕ−==∑将(6)式代入(4)式中，可以求出方程组(3)的解，然后将其代入(2)式中，最后求得方程组(1)的解为2()2()124212e ()(,)e/(4)i uxqE ux u t i n i i E En n s x L C x t u E kϕλ−∞−===++∑∑i2()[]4[ee ]()n u t t k En x λϕ−+− (7)为了下面表示方便，(7)式可以简化成函数映射关系1212(,)(,,,,,,,,,)q q C x t x t x x x s s s == M(,,)x t M m (8)其中M 表示对流-扩散方程(1)的解算子，m 表示由多个点污染源的位置i x 和排放强度i s ，(1,2,)i q = 构成的需要识别的未知向量。

在实际数值计算中通常要对解算子M 进行截断，截断阶数用N 表示。

3 MCMC 方法在贝叶斯统计理论中，将观测数据采集前所有关于未知参数向量m 的先验信息概率表述为先验分布()p m 。

获取观测后，根据对观测概率分布规律的了解，使用贝叶斯公式可将未知参数的先验分布改进为后验分布()p m |d ，即()()()()p p p p =m d d m m d (9)水 动 力 学 研 究 与 进 展 A 辑2009年第2期130其中(|)p d m 表示观测的条件概率密度。

d 是长度为M 的观测向量，本文中可表示为12(,,,)T Mobs obs obs C C C = d它包含了污染物在T 时刻不同位置的浓度观测。

因为观测数据已经给出，所以()p d 是一个与m 无关的常数，于是(9)式可写成()()()p p p ∝m d d m m (10)(10)式是进行贝叶斯推理的基础，通过它理论上可以获得参数的任何统计矩，如：每个参数的均值、方差和其它高阶统计量。

但实际应用中会遇到巨大困难：一方面除了非常简单的情况，后验概率密度都不存在明确的数学表达式；另一方面，采用通常的数值积分方法(如：Monte Carlo 方法)时，计算量将随未知向量维数的增加而呈指数增长。

因此贝叶斯方法几乎不能直接解决实际问题。

但是近期发展的马尔可夫链蒙特卡洛(MCMC)方法使得这种情况得到改善。

MCMC 算法可以对定义在高维随机向量空间上无明确数学表达式的概率分布p 进行抽样，其基本思想是产生大量服从分布p 的随机向量序列12I {,,,}L m m m ，其中I 为抽样数[22]。

如果向量序列满足马尔可夫性质：向量1i +m 的产生仅依赖于前一个向量i m ，而与1,2,,1i i −− 步骤的状态向量121,,,i i −− m m m 都无关，则该向量序列称为马尔可夫链。

马尔可夫性质的另一种描述是：若抽样算法当前访问的是j m 点，则下一步访问另一点i m 的概率只依赖于 j m ，而与先前访问的点无关。

马尔科夫性质意味着抽样算法完全可由转移概率矩阵P 描述，矩阵元素 ij p 表示算法在当前访问 j m 的条件下接着将要访问 i m 的条件概率。

按照构造Markov 链所用转移概率矩阵的不同，MCMC 方法的主要抽样算法有：Gibbs 抽样算法[24]、Metropolis-Hastings 算法[25]和自适应 Metropolis 算法[26]。