A辑第25卷第2期 水动力学研究与进展 Vol.25,No.2 2010年3月 CHINESE JOURNAL OF HYDRODYNAMICS Mar.,2010 D0L:10.3969j.issn.10004874.2010.02.001 对流一扩散方程源项识别 反问题的MCMC方法 曹小群，宋君强，张卫民，张理论 (国防科学技术大学计算机学院，长沙410073， Email:qunxiaocao2000@yahoo.com.cn) 痛要：给出了利用马尔科夫链蒙特卡罗(MCMC)方法求解对流-扩散方程源项识别反问题的一种新方法。该方法 把源项识别反问题视为贝叶斯估计问题，然后用MCMC法求解。首先利用贝叶斯公式，导出了多点源中源强和位置等 未知参数分布规律的后验概率密度函数：接着以后验概率分布为目标分布采用自适应Metropolis算法构造Markov链：然 后截取收敛的链序列，.利用后验均值方法估计出源项中的未知参数。数值试验结果表明，该方法具有精度高，收敛速 度快且易于计算机实现等优点。 关德阔：对流-扩散方程：源项：反问题：马尔科夫链蒙特卡罗方法 中图分类号：TV14 文献标识码：A MCMC method on an inverse problem of source term identification for convection-diffusion equation CAO Xiao-qun,SONG Jun-qiang, ZHANG Wei-min,ZHANG Li-lun (College of Computer,National Univ.of Defense Technology,Changsha 410073,China) Abstract:A new approach based on Markov Chain Monte Carlo(MCMC)method for source term identification of convection-diffusion equation was proposed.It views the inverse problem of source term identification as the problem of Bayesian estimation resolved by MCMC algorithm.Firstly,the posterior probability density function for unknown parameters of multiple point sources was deduced with the Bayesian formula.Secondly,taking the posterior probability as the target 幸收精日期：2009-03-192009-11-30修致稿) 基金项目：国家自然科学基金项目(40775064)资助 作者简介：曹小群(1980一)，男，湖南益阳人，助理研究员，博士. 万方数据
A辑第25卷第2期 20lO年3月 水动力学研究与进展 CHINESE JOURNALOF HYDRODYNAMICS V01．25．№．2 舰Br．．2010 DOI：lO．396州．issn．1000—4874．20lO．02．00l 对流一扩散方程源项识别 反问题的MCMC方法木 曹小群，宋君强，张卫民，张理论 (国防科学技术大学计算机学院，长沙410073， Emil：qun)【iaoca02000@yahoo．com．cn) 囊要：给出了利用马尔科夫链蒙特卡罗(MCMc)方法求解对流．扩散方程源项识别反问题的一种新方法。该方法 把源项识别反问题视为贝叶斯估计问题，然后用MCMC法求解。首先利用贝叶斯公式，导出了多点源中源强和位置等 未知参数分布规律的后验概率密度函数；接着以后验概率分布为目标分布采用自适应Me昀polis算法构造Markov链；然 后截取收敛的链序列，利用后验均值方法估计出源项中的未知参数。数值试验结果表明，该方法具有精度高，收敛速 度快且易于计算机实现等优点。 关■词：对流．扩散方程；源项；反问题；马尔科夫链蒙特卡罗方法 中圈分类号：Tvl4 文献标识码；A MCMC method on an inVerse problem of source term ●■ J●一 J● 一 J● 1●一一 ● J● ldentlIlCanon 10r C0nVeCnon-dnIUSlOn eqUatlOn CAO Xiao—qun， SONG Jun-qiang， ZH舢NG Wei—min， ZHj～NG Li—lun (College of Computer，National UniV．of Defense Technology，Challgsha 4 l 0073，China) Abstract：A new approach b勰ed 0n M砒刚Ch抽Monte Callo(MCMC)method for so眦e tem id即tification of conVection-di丘hsi∞equation w鹤pmpo辩d．It V洒Vs the inVer∞problem of s伽rce te珊idemification鸹the problem of Bayesi姐estin扭tion rcsolved by MCMc alg鲥thrn．F戤ly，Ⅱle posterior probabil时dcns时舢1ction f．0f蚰lmawn par锄et哪 ofmultiple p0砬∞llrces w髓dcduccd winl the Bayesi锄f砷咖la．Secondly，taking me postedor probabil时 船the target ·收稿日期l 2009．03．19(2009．1l一30修改稿) 基金项目：国家自然科学基金项目(407750“)资助 作者简介t曹小群(1980一)，男，湖南益阳人，助理研究员，博士． 万方数据
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=丁时的 分布已知，那么源项识别反问题就是根据这些已知 万方数据
曹小群，等：对流-扩散方程源项识别反问题的MCMC方法 129 的浓度分布观测来确定源项∑sδ(x一x),即确定 2e(s,-Pp..(x》 多个点污染源的位置和排放强度。 反问题在求解过程中通常需要将未知参数向 将(4)和(5)式代入(3)式中，则可求得T(t)的解为 量的估计值映射成观测空间的值，这就需要获得正 问题的解算子。本文采用的是Fourier方法来求解对 流-扩散方程，闵涛等1和吴自库等21分别采用相 同的方法对系统(1)进行了求解。首先引入下面的函 Z0-ri4归++2euo-e】何 B 数变换 其中 C(x,t)=V(x,t)e(E-4) (2) 2 B= sem.(x) 将其代入(1)中，则(1)可转化为等价的定解问题： L间 将(6)式代入(4)式中，可以求出方程组(3)的解，然 at ou-rn) 后将其代入(2)式中，最后求得方程组(1)的解为 V(0,)=0 V(L,t)=0 2e.() c=e品 V(x,0)=0 后u(4E)+2n+k (3) 由于方程组(3)的特征值和特征函数分别为 e0-e4R因 (7) n=气m/L)2 为了下面表示方便，()式可以简化成函数映射关系 和 C(x,t)=M(x,t,,x2,…,xg,52,…,5g）= pn(x)=sin(x/L),(n=1,2,) M(x,t,m) (⑧) 因此可利用特征函数展开法求解方程(3)。令 其中M表示对流-扩散方程(1)的解算子，m表示 由多个点污染源的位置x,和排放强度S,, rx,）=∑T.0p. (4) (i=1,2,…9)构成的需要识别的未知向量。在实际 数值计算中通常要对解算子M进行截断，截断阶数 用N表示。 将(3)中第一式的右端源项也利用特征函数进行展 开： 3MCMC方法 5) 在贝叶斯统计理论中，将观测数据采集前所有 其中 关于未知参数向量m的先验信息概率表述为先验 分布p(m)。获取观测后，根据对观测概率分布规 律的了解，使用贝叶斯公式可将未知参数的先验分 布改进为后验分布p(md),即 p(md)=p(d m)p(m)/p(d) (9) 万方数据
曹小群，等：对流一扩散方程源项识别反问题的McMc方法 的浓度分布观测来确定源项>^6@一石i)，即确定 j．1 多个点污染源的位置和排放强度。 反问题在求解过程中通常需要将未知参数向 量的估计值映射成观测空间的值，这就需要获得正 问题的解算子。本文采用的是Fo删er方法来求解对 流一扩散方程，闵涛等㈣和吴自库掣211分别采用相 同的方法对系统(1)进行了求解。首先引入下面的函 数变换 兰e笔7(喜s，e～c峨Ⅸ2E)，垆i(‘)) 129 将(4)和(5)式代入(3)式中，则可求得乙(f)的解为 【e“2f『(4扪一e—fI厶+。】】 (6) c(x，f)=y(x，彬枷胁’懈’ (2) ，'口 将其代入(1)中，则(1)可转化为等价的定解问题： 召2芝善s，e一嘶肥即识‘t’ 肾E警坩“矗Ⅲ刮2鄹扣H， {矿(O，，)=0 l矿(￡，，)=o 【矿(x，o)=o 矿(x，f)=∑乙(f)织(x) 刀=l 将(3)中第一式的右端源项也利用特征函数进行展 开： e舻r／(4E卜埘，(2鳓圭s，6(x—t)=主z纯(工) (5) j=l 行=l 其中 以=三re似2f，(4。w，(2∞’喜蚋@一t概(x)血= 将(6)式代入(4)式中，可以求出方程组(3)的解，然 后将其代入(2)式中，最后求得方程组(1)的解为 备e-(》酏) z，2／(4E)+以+后 [e(笔)一e叫¨1】纯(x) 为了下面表示方便，(7)式可以简化成函数映射关系 C(墨f)=膨(x，，，五，屯，…，％，毛，s2，…，sg)= M(石，f，肌) (8) 其中朋表示对流．扩散方程(1)的解算子，胁表示 由多个点污染源的位置薯和排放强度墨， (f=1，2，…g)构成的需要识别的未知向量。在实际 数值计算中通常要对解算子M进行截断，截断阶数 用Ⅳ表示。 3 MCMC方法 在贝叶斯统计理论中，将观测数据采集前所有 关于未知参数向量肼的先验信息概率表述为先验 分布p(朋)。获取观测后，根据对观测概率分布规 律的了解，使用贝叶斯公式可将未知参数的先验分 布改进为后验分布p(彤Id)，即 p(J"ld)=p(d I朋)p(朋)／p(口) 蒜再 。∑扁 卉石 一 = 旦蛆 C／L 而 D 万方数据
130 水动力学研究与进展 A辑2009年第2期 其中p(d|m)表示观测的条件概率密度。d是长 协方差矩阵在每一次迭代后都要自适应地调整。协 度为M的观测向量，本文中可表示为 方差的计算公式如(11)式所示，在初始。次迭代中， 协方差矩阵C,取固定值C。,之后进行自适应更新。 dr=(Ca,C2,…,CoM) 应，是未知模式参数向量m中某个元素在第i次迭 代的估计值，即 它包含了污染物在T时刻不同位置的浓度观测。因 为观测数据已经给出，所以p(d)是一个与m无关 [Co isio 的常数，于是(9)式可写成 C,= s Cov(mo,m)+sel i2i p(md)∝p(dm)p(m) (10) (11) (10)式是进行贝叶斯推理的基础，通过它理论上可 其中E是一个非常小的数，引入它是为了确保C 以获得参数的任何统计矩，如：每个参数的均值、 不成为奇异矩阵：S。是比例因子，它依赖于未知随 方差和其它高阶统计量。但实际应用中会遇到巨大 机向量空间的维数d,引入它是为了保证接受率在 困难：一方面除了非常简单的情况，后验概率密度 一个合适的范围内。本文分别取￡=106和 都不存在明确的数学表达式：另一方面，采用通常 sa=(2.4)2/d。La为d维单位矩阵。当进行第i+1 的数值积分方法（如：Monte Carlo方法）时，计算量 次迭代时，由公式(12)可导出协方差的计算公式， 将随未知向量维数的增加而呈指数增长。因此贝叶 斯方法几乎不能直接解决实际问题。但是近期发展 - 的马尔可夫链蒙特卡洛(MCMC)方法使得这种情况 Ci=- C,+(m-(+)mm+ i 得到改善。 MCMC算法可以对定义在高维随机向量空间 上无明确数学表达式的概率分布p进行抽样，其基 mm +el) (12) 本思想是产生大量服从分布P的随机向量序列 {m1,m2,L,m1},其中I为抽样数2四。如果向量序 其中m和m分别是前面i-1和i次迭代的参数 均值。自适应Metropolis算法的计算步骤如下： 列满足马尔可夫性质：向量m,的产生仅依赖于前 (1)设定i=0,对不同变量进行初始化。 一个向量m,而与i-1,i-2,…,1步骤的状态向量 (2)随机量的生成和接受，构造Markov链。 m,m-2,,m1都无关，则该向量序列称为马尔 ①利用公式(12)计算协方差矩阵C,: 可夫链。马尔可夫性质的另一种描述是：若抽样算 ②产生服从正态分布的推荐参数值 法当前访问的是m,点，则下一步访问另一点m,的 m~N(mi,C): ③利用下式计算接受概率 概率只依赖于m,而与先前访问的点无关。马尔 科夫性质意味着抽样算法完全可由转移概率矩阵 p(d m)p(m') P描述，矩阵元素P表示算法在当前访问m &min p(d m,)p(m) 的条件下接着将要访问m,的条件概率。按照构造 Markov链所用转移概率矩阵的不同，MCMC方法 ④产生服从均匀分布的随机数4~U(0,1): 的主要抽样算法有：Gibbs抽样算法24、 Metropolis-Hastings算法2阿和自适应Metropolis ⑤若<a,则接受m=m，否则 算法2询 m+1=m:。 本文采用自适应Metropolis算法，以对流-扩散 (3)重复上面的步骤①⑤，直到产生预先指定 方程源项未知参数的后验分布为目标分布来构造 数量的样本为止。 Markov链。自适应Metropolis算法是Haario在20OI 自适应Metropolis算法的最大优点是推荐分布 年提出的一种改进MCMC抽样算法2。与传统的 随计算过程自动更新，不再需要预先指定。同时， Metropolis-Hastings算法相比，自适应Metropolis?算 相比Metropolis-Hastings算法2，因为参数不再需 法不再需要预先确定参数的推荐分布，而是由后验 要分组更新，所以计算量大大减少。 参数的协方差矩阵来估算参数分布2。后验参数的 万方数据
水动力学研究与进展 A辑2009年第2期 其中p(d I朋)表示观测的条件概率密度。d是长 度为M的观测向量，本文中可表示为 dr=(E缸，c气，…，c篡) 它包含了污染物在丁时刻不同位置的浓度观测。因 为观测数据已经给出，所以p似)是一个与朋无关 的常数，于是(9)式可写成 p(J"ld)oc p(dI胧)p(朋) (10) (10)式是进行贝叶斯推理的基础，通过它理论上可 以获得参数的任何统计矩，如：每个参数的均值、 方差和其它高阶统计量。但实际应用中会遇到巨大 困难：一方面除了非常简单的情况，后验概率密度 都不存在明确的数学表达式；另一方面，采用通常 的数值积分方法(如：Monte Carlo方法)时，计算量 将随未知向量维数的增加而呈指数增长。因此贝叶 斯方法几乎不能直接解决实际问题。但是近期发展 的马尔可夫链蒙特卡洛(MCMC)方法使得这种情况 得到改善。 McMc算法可以对定义在高维随机向量空间 上无明确数学表达式的概率分布p进行抽样，其基 本思想是产生大量服从分布p的随机向量序列 {朋”朋，，厶朋。)，其中，为抽样数【22】。如果向量序 列满足马尔可夫性质：向量朋Ⅲ的产生仅依赖于前 一个向量册，，而与f—l，f一2，…，1步骤的状态向量 朋H，朋m，…，肌。都无关，则该向量序列称为马尔 可夫链。马尔可夫性质的另一种描述是：若抽样算 法当前访问的是m，点，则下一步访问另一点m，的 概率只依赖于肌i，而与先前访问的点无关。马尔 科夫性质意味着抽样算法完全可由转移概率矩阵 P描述，矩阵元素p。表示算法在当前访问朋， 的条件下接着将要访问眠的条件概率。按照构造 Markov链所用转移概率矩阵的不同，MCMC方法 的主要抽样算法有：Gibbs抽样算法例、 Me仃opolis—H嬲tings算法例和自适应Me仃0polis 算法【261。 本文采用自适应Me仃opolis算法，以对流．扩散 方程源项未知参数的后验分布为目标分布来构造 Mad【ov链。自适应Me仃opolis算法是Ha撕。在2001 年提出的一种改进MCMC抽样算法【26J。与传统的 Me仃('polis．Hastings算法相比，自适应Me仃opolis算 法不再需要预先确定参数的推荐分布，而是由后验 参数的协方差矩阵来估算参数分布【26】。后验参数的 协方差矩阵在每一次迭代后都要自适应地调整。协 方差的计算公式如(11)式所示，在初始乇次迭代中， 协方差矩阵C取固定值G，之后进行自适应更新。 碗是未知模式参数向量肼中某个元素在第f次迭 代的估计值，即 e《墓y。驴九m“，篇 (11) 其中￡是一个非常小的数，引入它是为了确保C 不成为奇异矩阵；劫是比例因子，它依赖于未知随 机向量空间的维数d，引入它是为了保证接受率在 一个合适的范围内。 本文分别取￡=10_6和 s。=(2．4)2肛。L为d维单位矩阵。当进行第f+1 次迭代时，由公式(12)可导出协方差的计算公式 e+。：掣c+粤(厩一。历：。一(f+1)厩砑+ Z j 聊f矿+￡厶) (12) 其中历一和而，分别是前面f一1和f次迭代的参数 均值。自适应Me仃opolis算法的计算步骤如下： (1)设定f=0，对不同变量进行初始化。 (2)随机量的生成和接受，构造Markov链。 ①利用公式(12)计算协方差矩阵C； ②产生服从正态分布的推荐参数值 m+一Ⅳ(m}，C)： ⑨利用下式计算接受概率 ～{t，甓矧 ④产生服从均匀分布的随机数“～U(O，1)； ⑤若“<口，则接受聊f+1=聊‘，否则 mtn2 Inl o (3)重复上面的步骤①橱，直到产生预先指定 数量的样本为止。 自适应Me仃opolis算法的最大优点是推荐分布 随计算过程自动更新，不再需要预先指定。同时， 相比Me仃0polis．H嬲tings算法【川，因为参数不再需 要分组更新，所以计算量大大减少。 万方数据
曹小群，等：对流-扩散方程源项识别反问题的MCMC方法 131 样的目标区域，有利于提高参数的估计精度。一般 来说，可以通过经验知识和历史统计确定更复杂和 4数值模拟试验 准确的先验分布。自适应Metropolis算法的一个优点 是对于m的任何先验分布都能够收敛于目标分布。 由式(13)、(14)和(15)可知，在给定观测数据的条件 在数值试验中，设定x和s,可由(7)式计算 下对流扩散方程中点源位置的后验概率密度函数 出C(x,t)在t=T时的浓度分布，并将其作为观测 必 资料C(x,T)。试验中取T=1。河道长度取L=1, 观测位置取为取(02,0.3，…，0.9)。式(7)中的截断阶 p(m d)= 数取为N=4000。试验表明N分别取4000和 (2πo)p 5000时，C(x,t)的误差小于机器误差。 第一种情形：假定污染源的强度S己知，且只 exp[-ld-M(x,T.m/(202)]U(x)(16) 考虑四个点污染源，需要识别污染源的位置 m=(名1，x2,x3,x4)。如果将公式(9)中分母表示的 未知污染源位置x,(i=1,2,…,4)的先验分布的形 常数项正规化为一，则对流扩散方程未知参数的后 验概率函数可以表示为 式为 1/L p(md)=p(dx,…,4)p(x,…,x4) 0<x<L (13) U(x)= (17) 0 otherwise 式(13)中已经假设d=M(x,T,m)+0,其中 M(x,T,m)表示对流-扩散方程的正模式，0为包 得到后验概率密度函数后，利用自适应Metropolis 含观测误差的独立分布的随机噪声，并且假定随机 算法按照计算步骤①~⑤来构造未知参数向量 噪声0服从均值为零、标准偏差为σ。的正态 m=(,x2,x3,x4)的Markov链，迭代次数预先指 定为2000，未知参数初始值取其对应先验均匀分布 分布，即0~N(0,σ)。同时假设由解算子M引 中的随机值。构造Markov链的过程实际上是在由 进的误差远小于观测误差。通过上述假设，可以得 先验分布界定的空间内对未知参数进行随机最优 到观测似然函数 搜索以使后验概率密度函数值最大化的过程。抽样 过程中未知参数m的每次更新，都要利用对流扩 1 L(d m)=- 散方程的解算子M计算出T时刻观测位置的观测 (2πo2)w2 相当量，以便计算后验概率密度函数的大小。各个 点源位置的Markov链如图1、图2、图3和图4所 exp[-d-M(x,T,m)/(202)] 14) 示，从图中可以看出，马尔科夫链在经历大约800 次左右的迭代后，链值就十分接近参数的真值，而 且收敛后的马尔科夫链具有平稳和分段光滑的特 式中n表示观察资料数量，H川表示欧几里得范数。 征。为了减少统计误差，取各个参数在1000-2000 在贝叶斯推理中，认为未知向量m是随机向 步之间的链序列求数学期望，分别得到污染源位置 量，需要设定每个未知参数的先验分布。在本文中 的后验均值，结果如表1所示。通过对精确值和估 假设所有未知参数的先验分布都满足独立的均匀 计值的对比分析可以得出两个结论：第一，从两者 分布，即参数(x,…,x4)分别满足由U(x)、 的差值可以看出，在无观测噪声的情况下，利用 U(x2)、U(x)和U(x4)表示的独立均匀分布。总 MCMC方法识别对流-扩散方程中点源位置至少具 的先验分布表可示为 有O10)的精度：第二，不同污染源位置的识别 精度不相同，其中x2和x的精度最高，达到了 p(m)=ΠU(x) O104),x4次之，精度在O(103)以上，x相对较 15) 差，精度接近O(102)。闵涛等采用遗传算法对类 似的源项识别反问题（试验设置不同）进行了研究 均匀分布是一种最简单的先验分布，虽然它只能指 ,方法具有收敛速度快、易于计算机实现和精度 定未知参数变化的上下界，但可以缩小参数随机抽 高等特点。实验中对两个源项位置的估计精度达到 万方数据
曹小群，等：对流一扩散方程源项识别反问题的McMc方法 4数值模拟试验 在数值试验中，设定五和s，，可由(7)式计算 出CO，f)在f=丁时的浓度分布，并将其作为观测 资料eh(x，丁)。试验中取r=1。河道长度取三=1， 观测位置取为取(O．2，0．3，…，0．9)。式(7)中的截断阶 数取为Ⅳ=4 000。试验表明Ⅳ分别取4 000和 5 000时，C(z，f)的误差小于机器误差。 第一种情形：假定污染源的强度鼠已知，且只 考虑四个点污染源，需要识别污染源的位置 朋=(五，恐，而，五)。如果将公式(9)中分母表示的 常数项正规化为一，则对流．扩散方程未知参数的后 验概率函数可以表示为 p(，”ld)=p(dl五，…，五)p(五，…，_) 13l 样的目标区域，有利于提高参数的估计精度。一般 来说，可以通过经验知识和历史统计确定更复杂和 准确的先验分布。自适应Me仃opolis算法的一个优点 是对于m的任何先验分布都能够收敛于目标分布。 由式(13)、(14)和(15)可知，在给定观测数据的条件 下对流一扩散方程中点源位置的后验概率密度函数 为 咖∽2面匆 exp[一od一^f(x，丁，胁)112／(2刃)]．血u(誓)(16) 未知污染源位置薯(f-1，2，…，4)的先验分布的形 式为 ∞，吣，彤慧主e ∽， 式(13)中已经假设d=M(x，丁，柳)+彩，其中 膨(x，丁，柳)表示对流．扩散方程的正模式，缈为包 含观测误差的独立分布的随机噪声，并且假定随机 噪声国服从均值为零、标准偏差为以的正态 分布，即缈～Ⅳ(o，一)。同时假设由解算子M引 进的误差远小于观测误差。通过上述假设，可以得 到观测似然函数 圳∽2面扣 exp[一8d—M(x，丁，胁)112／(2《)】 (14) 式中，z表示观察资料数量，…表示欧几里得范数。 在贝叶斯推理中，认为未知向量柳是随机向 量，需要设定每个未知参数的先验分布。在本文中 假设所有未知参数的先验分布都满足独立的均匀 分布，即参数(而，…，矗)分别满足由U(五)、 U(而)、U(屯)和U(矗)表示的独立均匀分布。总 的先验分布表可示为 4 p(朋)=nu(t) f=l (15) 均匀分布是一种最简单的先验分布，虽然它只能指 定未知参数变化的上下界，但可以缩小参数随机抽 得到后验概率密度函数后，利用自适应Metropoli8 算法按照计算步骤①～⑤来构造未知参数向量 肌=(五，荔，五，兄)的Markov链，迭代次数预先指 定为2000，未知参数初始值取其对应先验均匀分布 中的随机值。构造Markov链的过程实际上是在由 先验分布界定的空间内对未知参数进行随机最优 搜索以使后验概率密度函数值最大化的过程。抽样 过程中未知参数朋的每次更新，都要利用对流．扩 散方程的解算子M计算出丁时刻观测位置的观测 相当量，以便计算后验概率密度函数的大小。各个 点源位置的Markov链如图1、图2、图3和图4所 示，从图中可以看出，马尔科夫链在经历大约800 次左右的迭代后，链值就十分接近参数的真值，而 且收敛后的马尔科夫链具有平稳和分段光滑的特 征。为了减少统计误差，取各个参数在1000—2000 步之间的链序列求数学期望，分别得到污染源位置 的后验均值，结果如表1所示。通过对精确值和估 计值的对比分析可以得出两个结论：第一，从两者 的差值可以看出，在无观测噪声的情况下，利用 MCMC方法识别对流．扩散方程中点源位置至少具 有D(10。3)的精度；第二，不同污染源位置的识别 精度不相同，其中而和黾的精度最高，达到了 D(10一)，无次之，精度在D(10。)以上，五相对较 差，精度接近D(10-2)。闵涛等采用遗传算法对类 似的源项识别反问题(试验设置不同)进行了研究 【l 91，方法具有收敛速度快、易于计算机实现和精度 高等特点。实验中对两个源项位置的估计精度达到 万方数据