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，方法具有收敛速度快、易于计算机实现和精度 高等特点。实验中对两个源项位置的估计精度达到 万方数据
132 水动力学研究与进展 A辑2009年第2期 了OI0)。本文的实验结果中对x2和x的估计精 0.4 度与其相当，但另外两个点的估计精度稍差，主要 Q.42 是因为MCMC方法具有一定的随机性。进一步的 实验结果表明，当Markov链迭代次数增加到4000 0.40 时，对源项位置的估计精度都能到达O104),见 0.38 型 表2。在计算效率方面，闵涛等采用的遗传算法9) 需要产生1000代包含50个个体的群体，即需要生 034 成50000个样本：而MCMC方法当前的估计量只 人 依赖于前一个估计向量，即使采用长度为4000的 0.30 Markov链，即对于单个未知参数也只需要产生4000 0.28 个样本。另外在每次迭代中MCMC方法只需要计 0.26 算和比较概率值的大小，而遗传算法需要在群体中 0 400 8001200 1600 2000 进行交叉、选择和变异等步骤。因此MCMC方法 取样次数 的CPU计算时间比遗传算法要少，计算效率占优。 图2污染源位置x2的Markov链 样本标准偏差衡量了统计样本相对于估计值的发 散程度。从表1和2中可以看出各个位置参数估计 0.62 量的标准偏差量级都等于或小于(Q0),这从一个 方面表明了MCMC方法对点源位置识别的收敛性。 0.60 另外，数值试验结果还表明：由于MCMC方法是 一种基于随机模拟的算法，因此估计结果具有随机 性，在同样条件下不同试验产生的结果会不相同。 0.56 0.22r 0.54 0.20F 0.52 0.18 0.16 0.50 0 400 800120016002000 取样次数 0.14 0.12 图3污染源位置x3的Markov链 0.10 0.90 0.08 0 00 800 12001600 2000 取样次数 0.85 图1污染源位置x的Markov链 举080 表1选代2000步时污染源位量识别的结果 0.75 为 X2 3 XA 0.70 精确值 0.2 0.4 0.6 0.8 0.65 400 800120016002000 估计值 0.20190.39990.59980.7992 取样次数 样本标准偏差0.00210.00500.00500.0021 图4污染源位置x4的Markov链 万方数据
132 水动力学研究与进展 A辑2009年第2期 了D(10一)。本文的实验结果中对艺和黾的估计精 度与其相当，但另外两个点的估计精度稍差，主要 是因为MCMc方法具有一定的随机性。进一步的 实验结果表明，当Markov链迭代次数增加到4000 时，对源项位置的估计精度都能到达D(10q)，见 表2。在计算效率方面，闵涛等采用的遗传算法【19J 需要产生1000代包含50个个体的群体，即需要生 成50000个样本；而MCMC方法当前的估计量只 依赖于前一个估计向量，即使采用长度为4000的 Markov链，即对于单个未知参数也只需要产生4000 个样本。另外在每次迭代中MCMC方法只需要计 算和比较概率值的大小，而遗传算法需要在群体中 进行交叉、选择和变异等步骤。因此MCMC方法 的CPU．计算时间比遗传算法要少，计算效率占优。 样本标准偏差衡量了统计样本相对于估计值的发 散程度。从表l和2中可以看出各个位置参数估计 量的标准偏差量级都等于或小于(10q)，这从一个 方面表明了McMc方法对点源位置识别的收敛性。 另外，数值试验结果还表明：由于MCMC方法是 一种基于随机模拟的算法，因此估计结果具有随机 性，在同样条件下不同试验产生的结果会不相同。 0．兹 o．加 o．18 嚣 舂o．16 轧h o．12 o．10 o．∞ O 400 咖 1200 l鲫枷0 取样次数 图l污染源位置五的Mad【0v链 裹1迭代2000步时污染源位I识别的结果 o-62 o．∞ o-弱 掰 扣e 辗 一0．54 0．52 0．50 图2污染源位置屯的M破oV链 O 枷 800 1200 l啪2鲫 取样次数 图3污染源位置屯的MarkoV链 0 枷 800 1200 1600 2咖 取样次数 图4污染源位置丘的Markov链 ∞ {8 鼢 {2 ∞ 曲 n n n n n n Ⅻ警量蛊．k 万方数据
曹小群，等：对流-扩散方程源项识别反问题的MCMC方法 133 表2选代4000步时污染源位置识别的结果 3 0.45 精确值 0.2 0.4 0.6 0.8 0.40 035 估计值 0.20010 0.39995 0.59997 0.80002 030 样本标准偏差 0.00013 0.00026 0.00012 0.00011 0.25 第二种情形：假定污染源项中的位置x,及强度 020 3,均未知，但只考虑两个点污染源，即需要利用观 0.1 400 800120016002000 测数据确定x及，i=1,2。此时把需要识别的未 取样次数 知向量表示为m=(:,x2,1,S2),则m的后验概率 图5污染源位置方的Markov链 函数可以表示为 7.0 p(md)=p(d,x2,1,S2)p(x1,x2,,S2)(18) 65 观测似然函数采用(14)式，参数m=(:,x2,S,S2) 分别满足U(s)、U(s2)、U(:)和U(x)表示的 ao 独立均匀分布。总的先验分布可表示为 p(m)=IIU(X(s) (19) 5.0 由式(14)、(18)和(19)可知，在给定观测数据的条件 00 8001200 16002000 下对流-扩散方程中两个点源未知参数的后验概率 取样次数 密度函数为 图6污染源强度S,的Markov链 m=(1,x2,S,S2)的Markov链，迭代次数预先指 p(mld)-2yexpt--M(.T.m)/ 定为2000，未知参数初始值取其对应先验均匀分布 中的随机值。图5、6、7和8所示的是(：1，x2,S1,S2) 的Markov链。从图中可以看出，马尔科夫链在经 (2U(U(s) (20) 历大约800次左右的迭代后达到收敛，且收敛后的 马尔科夫链具有平稳和光滑的特征：污染源强度不 到800次迭代就能收敛，且在真实值附近的振荡很 未知污染源位置x(i=1,2)的先验分布取(17)式， 小。为了减少统计误差，取1000-2000步的链序列 污染源强度5，(i=1,2)的先验分布取为 计算各个未知参数的数学期望，得到未知参数的后 验估计值，识别结果如表3所示。另外，还利用 1/10 Us)=0 0<3,<10 1000-2000样本序列计算了各个参数的标准偏差， (21) 量级为O103),这在一定程度上验证了MCMC方 otherwise 法同时识别点源位置和源强参数时的收敛性。从与 闵涛等采用遗传算法得到的估计结果比较来看， 得到后验概率密度函数后，利用自适应Metropolis 点源位置的估计精度比其稍差，点源强度的估计精 算法按照计算步骤①⑤来构造未知参数向量 度相当：但如果将Markov链的长度增加到4000， 不但点源位置估计精度能达到O10),而且源强 也能达到同样的估计精度（见表4）。因此，结合前面 万方数据
曹小群，等：对流一扩散方程源项识别反问题的McMc方法 五 恐 屯 矗 精确值 O．2 0．4 0．6 0．8 估计值 0．200 10 O．399 95 O．59997 0．80002 样本标准偏差 0．000 13 O．00026 O．000 12 0．ooOll 第二种情形：假定污染源项中的位置五及强度 量均未知，但只考虑两个点污染源，即需要利用观 测数据确定五及墨，f=1，2。此时把需要识别的未 知向量表示为朋=(五，屯，Jl，s2)，则朋的后验概率 函数可以表示为 。 p(肌ld)=p(dljcl，恐，sl，s2)p(_，jc2，sl，s2) (18) 观测似然函数采用(14)式，参数朋=(五，恐，■，s2) 分别满足【，(量)、U(s2)、U(五)和U(而)表示的 独立均匀分布。总的先验分布可表示为 (19) 由式(14)、(18)和(19)可知，在给定观测数据的条件 下对流．扩散方程中两个点源未知参数的后验概率 密度函数为 咖护面知eXp[一卜M(妒圳2／ 一 2 (2《)】．兀u(五)u(墨) (20) f=I 未知污染源位置薯(f=l，2)的先验分布取(17)式， 污染源强度s，(f=l，2)的先验分布取为 fl／10 叭墨户佰 0<s：<10 oth未ise (21) 得到后验概率密度函数后，利用自适应Me仃opolis 算法按照计算步骤∞来构造未知参数向量 0．的 o．45 o．20 o．15 7．0 矗5 掣＆o 墨 皇 g 5．5 5．0 生5 O 400 8∞ 1200 1600 2咖 取样次数 图5污染源位置五的Markov链 0 400 800 1200 1600 2000 取样次数 图6污染源强度sl的MarkoV链 133 朋=(五，五，s，岛)的Markov链，迭代次数预先指 定为2000，未知参数初始值取其对应先验均匀分布 中的随机值。图5、6、7和8所示的是(五，％，量，s，) 的Markov链。从图中可以看出，马尔科夫链在经 历大约800次左右的迭代后达到收敛，且收敛后的 马尔科夫链具有平稳和光滑的特征；污染源强度不 到800次迭代就能收敛，且在真实值附近的振荡很 小。为了减少统计误差，取1000~2000步的链序列 计算各个未知参数的数学期望，得到未知参数的后 验估计值，识别结果如表3所示。另外，还利用 1000乏000样本序列计算了各个参数的标准偏差， 量级为D(10一3)，这在一定程度上验证了McMC方 法同时识别点源位置和源强参数时的收敛性。从与 闵涛等采用遗传算法得到的估计结剁19J比较来看， 点源位置的估计精度比其稍差，点源强度的估计精 度相当；但如果将Markov链的长度增加到4000， 不但点源位置估计精度能达到D(10一)，而且源强 也能达到同样的估计精度(见表4)。因此，结合前面 ∞ ：穹 ∞ 龉 n n n n 端，卫皇‘宅．k U 薯 、， U0 、， ：n阿 lI p，L 肼 、， 万方数据
134 水动力学研究与进展 A辑2009年第2期 的分析可知，在对点源位置和强度进行同时识别的 表4选代4000步时对污染源位量和强度同时识别的结果 情况下，如果要达到同样的估计精度，MCMC方法 的计算效率要优于遗传算法。总的来说，第二种情 为 X2 2 形的试验结果与第一种情形类似，说明了利用 MCMC方法能同时识别对流-扩散方程中源项的位 精确值 0.2 0.4 5.6 2.8 置和强度，而且具有较高的估计精度。 估计值 0.2011 0.399985.599982.80003 0.5刷f 0.52 样本标准偏差 0.000130.000070.000110.00015 0.50 0阳 5结束语 0.42卡 本文在贝叶斯理论的框架下，提出了一种用 MCMC方法识别对流-扩散方程中源项的新方法。 0.38 首先利用贝叶斯公式，导出了多点源中未知参数分 0.36 400 800120016002000 布规律的后验概率密度函数。然后采用自适应 取样次数 Metropolis算法，以后验分布为极限不变分布构造了 图7污染源位置x,的Markov链 Markov链，并截取收敛后的样本序列进行分析，成 功识别了对流扩散方程中源项的未知参数。主要结 论归纳如下： 3.6 (I)利用MCMC方法能准确地识别对流-扩散 3.4 方程源项中的多个未知参数，未知参数的估计精度 3.2 较高，从而验证了用MCMC方法解决对流-扩散方 程源项识别反问题的有效性。 3.0 (2)利用自适应Metropolis算法构造未知参数 2.8 的Markov链能较快向参数真实值靠近。通过对样本 标准偏差的分析可知，Markov链具有收敛特性，且 2.6 收敛后的Markov链是平稳和光滑的。 2.4 (3)由于MCMC方法是一种基于随机模拟的算 法，因此估计结果具有随机性，在同样条件下不同 2.2 0 00 800 1200 1600 2000 试验产生的结果可能会不相同。 取样次数 总之，MCMC方法是一种通过构造Markov链进 图8污染源强度S,的Markov链 行随机模拟的方法。构造Markov链的过程实际上是 在由先验分布界定的空间内对未知参数进行随机 最优搜索的过程。本文的研究结果表明，MCMC方 表3选代2000步时对污染源位置和强度同时识别的结果 法能准确识别对流-扩散方程源项中的多个未知参 数，同时具有较高的估计精度。 为 S2 精确值 0.2 0.4 5.6 2.8 考文献： 估计值 0.20100.3984 0.5985 2.7951 [ 王彩华.一维对流扩散方程的一类新型高精度紧致差 样本标准偏差0.00040.0001 0.0082 0.0022 分格式[)。水动力学研究与进展，A辑，2004,19(5： 655-663. 万方数据
水动力学研究与进展 A辑2009年第2期 的分析可知，在对点源位置和强度进行同时识别的 情况下，如果要达到同样的估计精度，McMc方法 的计算效率要优于遗传算法。总的来说，第二种情 形的试验结果与第一种情形类似，说明了利用 MCMC方法能同时识别对流．扩散方程中源项的位 置和强度，而且具有较高的估计精度。 O 400 800 1200 16∞ 2000 取样次数 图7污染源位置而的Markov链 3．6 3．4 3．2 撼 呈3·o 塞2．。 铀 2．6 Z4 2．2 O 400 800 1200 1600 2咖 取样次数 图8污染源强度s2的M棚链 衰3迭代2000步时对污染源位置和置度同时识别的结果 裹4迭代4000步时对污染源位置和蟊度同时识别的结果 5结束语 本文在贝叶斯理论的框架下，提出了一种用 MCMc方法识别对流．扩散方程中源项的新方法。 首先利用贝叶斯公式，导出了多点源中未知参数分 布规律的后验概率密度函数。然后采用自适应 MenDpolis算法，以后验分布为极限不变分布构造了 Markov链，并截取收敛后的样本序列进行分析，成 功识别了对流一扩散方程中源项的未知参数。主要结 论归纳如下： (1)利用MCMC方法能准确地识别对流一扩散 方程源项中的多个未知参数，未知参数的估计精度 较高，从而验证了用MCMC方法解决对流．扩散方 程源项识别反问题的有效性。 (2)利用自适应Me仃opolis算法构造未知参数 的Markov链能较快向参数真实值靠近。通过对样本 标准偏差的分析可知，Markov链具有收敛特性，且 收敛后的Markov链是平稳和光滑的。 (3)由于McMC方法是一种基于随机模拟的算 法，因此估计结果具有随机性，在同样条件下不同 试验产生的结果可能会不相同。 总之，MCMC方法是一种通过构造MarkoV链进 行随机模拟的方法。构造Markov链的过程实际上是 在由先验分布界定的空间内对未知参数进行随机 最优搜索的过程。本文的研究结果表明，McMC方 法能准确识别对流一扩散方程源项中的多个未知参 数，同时具有较高的估计精度。 参考文献： 【l】 王彩华．一维对流扩散方程的一类新型高精度紧致差 分格式【J】．水动力学研究与进展，A辑，2004，19(5)： 655．663． 矾 匏 ∞ 鸺 稻 “ 眈 ∞ 嚣 {5； 仉 n n n n n n n n 仉 埘>罩毒雹HK 万方数据
曹小群，等：对流-扩散方程源项识别反问题的MCMC方法 135 WANG Cai-hua.A class of high-order compact for nonlinear convection-diffusion equations[J.Acta difference scheme for a convection-diffusion equation[. Scientiarum Naturalium Universitatis Nankaiensis,2006, Journal of Hydrodynamics,Ser.A,2004,19(5):655-663. 392:40-45. [2 葛水斌，田振夫，吴文权.含源项非定常对流扩散方 [10)黄思训，伍荣生.大气科学中的数学物理问题M, 程的高精度紧致隐式差分方).水动力学研究与进 北京：气象出版社，2001. 展，A辑.2006,21(5)：619-625. HUANG Si-xun,WU Rong-sheng.Mathematical GE Yong-bin,TIAN Zhen-fu,WU Wen-quan.A Physics Problems in Atmosphere Science[M].Beijing: high-order compact implicit difference method for the Meteorology Press,2001. unsteady convection-diffusion equation with source [11]KIRSCH A.An introduction to the mathematical theory term[J].Journal of Hydrodynamics,Ser.A,2006,21(5): of inverse problem[M].New York:Springer-Verlag, 619-625. 1996. 31 杨雪源。王彩华，齐海涛，等，对流扩散问题的一种 [12] YILDIZ B,YETISKIN H,SEVER A.A stability 紧致差分方法).水动力学研究与进展，A辑，2008， estimate on the regularized solution of the backward 23(4:426437. heat equation[].Applied Mathematics and Computation, YANG Xue-yuan,WANG Cai-hua,QI Hai-tao,et al.A 2003,135(2):561-567. compact finite difference method for convection- [13]LIU J.Numerical solution of forward and backward diffusion problem[J].Chinese Journal of problem for 2-D heat conduction equation [J].Journal of Hydrodynamics,2008.23(4):426437. Computational and Applied Mathematics,2002,145(2): [4 汪继文，窦红.求解对流扩散方程的一种高效的有限 459-482. 体积法).应用力学学报，2008,25(3：480-484. [14]QUN Chen,LIU Ji-jun.Solving an inverse parabolic WANG Ji-wen,DOU Hong.High-efficient finite problem by optimization from final measurement data[J]. volume method for advection-diffusion equations[J]. Journal of Computational and Applied Mathematics, Chinese Joumal of Applied Mechanics,2008,25(3): 2006,193(1:183-203. 480-484. [15]LIU Ji-jun.On uniqueness and linearization of an [5) WANG Tong-ke.A class of upwind finite volume inverse electromagnetic scattering problem[J].Applied schemes for one-dimensional steady convection- Mathematics and Computation,2005,171(1):406-419. dominated diffusion equations[.Mathematica [16]LIU Ji-jun,SEO J.On stability for a translated obstacle Applicata.2004,17(4):544-550. with impedance boundary condition[J].Nonlinear [6 高智，柏戚.对流扩散方程的摄动有限体积方法及讨 Analysis,2004,59S:731-744. 论0.力学学报，2004,36(1少88-93. [1刀闵涛，周孝德.河流水质纵向弥散系数反问题的迭代 GAO Zhi,BAI Wei.Perturbational finite volume 算法[J.水动力学研究与进展，A辑，2003,18(5)： method for convective diffusion equation and 547.552 dissussion[J.Acta Mechanica Sinica,2004,36(1): MIN Tao,ZHOU Xiao-de.An iteration method of the 88-93. inverse problem for the dispersion coefficient in water [7] 孙澈，秦树杰.非定常线性对流扩散问题的算子分裂 quality model[J].Journal of Hydrodynamics,Ser.A, 半显式有限元分析[).计算数学，2003,25(1片23-34. 2003,18(5):547-552. SUN Che,QIN Shu-jie.Analysis of operator splitting [18】闵涛，刘相国，张海燕，等.二维稳态对流-扩散方程 semi-explicit finite element method for time-dependent 参数反演的迭代算法).水动力学研究与进展，A辑， linear convection-diffusion problems[J.Mathematica 2007,22(6):744-752. Numerica Sinica,2003,25(1):23-34. MIN Tao,LIU Xiang-guo,ZHANG Hai-yan,et al. [8] 超志勇，胡健伟，孙琳.对流扩散方程迎风有限元的 Iterative algorithms for the inversion of parameter in a 自适应方法[.计算数学，2005,27(4：337-354. two-dimensional stable-state convection-diffusion ZHAO Zhi-yong,HU Jian-wei,SUN Lin.Adaptive equation[J].Journal of Hydrodynamics,Ser.A,2007, upwind finite element method for convection-diffusion 22(6):744-752. equations[J].Mathematica Numerica Sinica,2005,27(4): [19]闵涛，周孝德，张世梅，等。对流扩散方程源项识别 337-354. 反问题的遗传算法[).水动力学研究与进展，A辑， [9]赵志勇，胡健伟.非线性对流扩散方程迎风有限元的 2004,19(4):520-525. 自适应方法).南开大学学报，2006,39(2)：40-45. MIN Tao,ZHOU Xiao-de,ZHANG Shi-mei,et al. ZHAO Zhi-yong,HU Jian-wei.Adaptive upwind FEM Genetic algorithm to an inverse problem of source term 万方数据
曹小群，等：对流～扩散方程源项识别反问题的McMC方法 Ⅵ7ANG Cai—lm．’A class of Kgh．order compact dif断eIlce scheme 6吖a convectio小dimlsion equation阴 Jo哪al of HydrodynarIlics，S既A，2004，19(5)：655—663． 【2】 葛永斌，田振夫，吴文权．含源项非定常对流扩散方 程的高精度紧致隐式差分方法m．水动力学研究与进 展，A辑，2006，2l(5)：619·625． GE Yong-bin，1’队N ZheIl-fil，WU W曲—q呦．A high-ordcf c彻叩act implicit di伍珊mcc Imthod fbr the 吼steady convection-dif如sion eqmtion with∞urce tenn【J】．Jo呦l of Hydro由咖ics，Ser．A，2006，2l(5)： 619—625． 【3】 杨雪源，王彩华，齐海涛，等．对流扩散问题的一种 紧致差分方法【J】．水动力学研究与进展，A辑，2008， 23(4)：426-437． YANG Xue—”啪，WANG Cai—hua，QI Hai—tao，et a1．A compact f．mite difI．erence meIhod fbf ∞nvecti∞． dif如si∞pfoblem[J】． ChiIlese Joumal of Hydrody瑚mics，2008，23(4)：426-437． 【4】 汪继文，窦红．求解对流扩散方程的一种高效的有限 体积法【J】．应用力学学报，2008，25(3)：480-484． WANG Ji．w锄，DOU Hong． High-e伍cient fi出 vol啪e me删f．0r adve甜on-difmsion equations『J]． Chinese Jo哪l of Applied Mech柚ics，2008，25(3)： 480_484． [5】 wANG Tc’ng—ke．A class of upwind finite volu眦 schemes for o北．dimensio∞I 毗ady convection_ domina自ed djahsjon equati衄s叨．Ma岫石ca Applicata，2004，17(4)：544—550． [6】 高智，柏威．对流扩散方程的摄动有限体积方法及讨 论叨．力学学报，2004，36(1)：88．93． GAO Zhi． BAI Wei． Perturbational finite volume method for convective difmsion equation 锄d diss璐sion【J】．Acta Mechanica Sinica，2004，36(1)： 88．93． 。 【7】 孙澈，秦树杰．非定常线性对流扩散问题的算子分裂 半显式有限元分析叨．计算数学，2003，25(1)：23．34． SUN Che，Q矾Shu．jie．Analysis of opemtor spli仕吨 ∞Ini·expIicit fitl慨elemem metllod f．or time-d印endent l№ar convection面肺si蚰probl锄s阴．Math锄atica Num甜ca Sillica，2003，25(1)：23—34． 【8】 赵志勇，胡健伟，孙琳．对流扩散方程迎风有限元的 自适应方法[J】．计算数学，2005，27(4)：337—354． ZHAO Zhi-yon岛HU Ji柚-wei，SUN Lin．A峨ptivc upw枷蠡nite eleInent嘣hod f-0f convection．difmsi∞ equatio衄【J】．Mathenlatica Numerica SiIlica，2005，27(4)： 337-354． f9】 赵志勇，胡健伟．非线性对流扩散方程迎风有限元的 自适应方法叨．南开大学学报，2006，39(2)：40-45． ZHAo Zhi—yon＆HU Ji粕一wei．Ada嘶ve upw砌FEM 135 for nonlinear convecti∞．did’usion eqlla瞳io璐[J】．Acta Sciemianlm Natll鞠Ii啪Universitatis NaIlkaiensis。2006。 39(2)：40-45． [1 o】黄思训，伍荣生．大气科学中的数学物理问题[M】． 北京：气象出版社．2001． HUANG Si—x眦，、ⅣU Rong．sh％g． Mamer∞tical Physics Problems in Annosphe陀Sci钮cc[M】．Beijing： Meteorolog)r Press，2001． 【11】 ⅪRSCH A．An in响duction to t11e删lthematical ttleo叮 of研erse probl锄【M】．N哪Y破：s阿ng*Venag， 1996． 【12】 ⅥLDIz B，ⅥmSⅪN H，SEⅦR A．A鼢il时 estitnate on the regularized solution of the backward heat equation【J】．Applied Mathematics and C．)mputation' 2003，135(2)：56l一567． 【13】 LⅣJ．N啪erical soluti∞of fo刑ard加d backward problem for 2一D hcat conducti∞e叩眦ion叨．Jo明诅l of Comp眈Ltio∞l锄d Applied Mathematics，2002，145(2)： 459-482． 【14】QUN Ch％，LⅢ州吼．solv吨狮illverse p啪bolic problem by唯im妇i∞fbm f-mal me髂urcmem data【J】． Jo硼1al of C唧蝴io∞l龃d Applied M砒ematics， 2006，193(1)：183-203． 【15】 LIU Ji0啦．on lllliqu∞ess and linearization of姐 inVer∞electromagnetic scatte]riI唱probl锄【J】．Applied Math锄atics柚d Comp吡ltion，2005，17l(1)：406419． 【16】 Lm Ji-j咖，SE0 J．0n s油ilit)r for a柏璐la刚obstacle with impedance bollIIda巧 conmtion【J】． Nonline缸 Amalysis，2004，59(5)：731-744． 【17】闵涛，周孝德．河流水质纵向弥散系数反问题的迭代 算法[J】．水动力学研究与进展，A辑，2003，18(5)： 547．552． M烈T∞．ZHOU xi∞-de．An ite|吼ion method of me inver∞probl锄for the dispersion coefncient in w栅 q眦l畸model川．JouⅡml of Hydro蛳ics，ser．～ 2003，18(5)：547—552． 【18】闵涛，刘相国，张海燕，等．二维稳态对流．扩散方程 参数反演的迭代算法川．水动力学研究与进展，A辑， 2007，22(6)：7“-752． MIN T∞，LIU)【i锄g—gllo，ZHANG Hai-y蛐，et a1． I删ve algorithms for the inve璐ion of paramcter in a two—dimensional stable—sta_tc convection．di肺sion equation【J】．J0响扭l of HydrodynaⅡlics，Ser．～2007， 22(6)：744．752． 【19】闵涛，周孝德，张世梅，等．对流．扩散方程源项识别 反问题的遗传算法叨．水动力学研究与进展，A辑， 2004，19(4)：520-525． MIN T的，ZHOU)(i∞-de，ZHANG Shi．mei．et a1． (艳netic alg硎Ⅱlln t0矩inV蹦沁pmbIem of∞urcc telln 万方数据
136 水动力学研究与进展 A辑2009年第2期 identification for convection-diffusion equation(J]. Markov chain monte carlo in practice[M].London: Journal of Hydrodynamics,Ser.A,2004,19(4): Chapman Hall,1996. 520-525. [23]曹小群，张卫民，宋君强，等.海气振子系统中未知 [20潘军峰，闵涛，周孝德，等.对流-扩散方程逆过程反 参数的MCMC方法识别[U.物理学报，2009,58(9： 问题的稳定性及数值求解).武汉大学学报（工学 6050-6057. 版)，2005,38(1：10-13. CAO Xiao-qun,ZHANG Wei-min,SONG Jun-qiang,et PAN Jun-feng,MIN Tao,ZHOU Xiao-de,et al.Stability al.Estimating parameters of air-sea oscillator with analysis and numerical solution of inverse problem in MCMC method[J].Acta Physica Sinica,2009,58(9): reverse process of convection-diffusion equation[. 6050-6057. Engineering Journal of Wuhan University,2005,38(1): [24]SMITH A F M,ROBERT G O.Bayesian computation 10-13. via the gibbs sampler and related markov chain monte [21】吴自库，范海梅，陈秀荣，等.对流扩散方程逆过程 carlo methods[J].Journal of Royal Statistical Society, 反问题的伴随同化研究U.水动力学研究与进展，A Ser.B,1993,551:3-23. 辑，2008,23(2121-125. [25] TIERNEY L.Markov-chains for exploring posterior WU Zi-ku,FAN Hai-mei,CHEN Xiu-rong,et al.The distributions[J].Annals of Statistics,1994,22(4): numerical study of the inverse problem in reverse 1701-1762. process of convection-diffusion equation with adjoint [26 HAARIO H,SAKSMAN E,TAMMINEM J.An assimilation method[J].Chinese Journal of adaptive metropolis algorithm[J].Bernoulli,2001,7(2): Hydrodynamics,Ser.A,2008,23(2):121-125 223-242. [22]GILKS W R,RICHARDSON S,SPIEGELHALTER D J. 万方数据
id舶tification fbr conve商on-di航sion equation叨． Jo啪al of HydrodyIlalIlics，S盯．八2004，19(4)： 520．525． 【20】潘军峰，闵涛，周孝德，等．对流．扩散方程逆过程反 问题的稳定性及数值求解田．武汉大学学报(工学 版)，2005，38(1)：10一13． PANJ衄-feng，MIN TaIo，ZHOU Xiao·de，et a1．S切Ibility 鲫alysis锄d numerical solution of inver：辩probl锄in 化verse process of convection—difmsion eqmtion【J】． Enginecring Jo哪l of Wuh锄UniVersity，2005，38(1)： 10—13． 【2l】吴自库，范海梅，陈秀荣，等．对流-扩散方程逆过程 反问题的伴随同化研究田．水动力学研究与进展，A 辑，2008，23(2)：121-125． ⅥrU Zi—l[Il，FAN Hai—mei，CHEN Xiu-rong，et a1．The numerical study of the iIlverse problem in RVer∞ pmcess of convection-difmsion eq眦ion with adjoint 勰siIIlilation method川． Chi鹏sc JoIlrnal of HydrodymInics，Scr．A，2008，23(2)：12l-125． [22】 GILKS WR，lUCH■J5U)SON S，SPIEGELHALTER D J． Mad【ov cllain monte cano in pmctice[M】．London： Cllapn瑚&Hall，1996． [23】曹小群，张卫民，宋君强，等．海-气振子系统中未知 参数的McMc方法识别【刀。物理学报，2009，58(9)： 6050一6057． CAO Xiao．qun，Z地～NG Wei—m．m，SONG J岫一qi粕g，et a1．Estimating讲咖etcrs of air．∞a oscillator wim MCMC memod明．Acta Physica SiIlica，2009，58(9)： 6050．6057． 【24】 SMITH A F M，ROBERT G O．Bayesi锄computatlon via the gibbs s卸叩ler锄d rclated m诎ov chain monte carlo methods【J]．Jo啪al of Royal Statistical Society， Ser．B，1993，55(1)：3．23． 【25】 TIERNEY L．Markov—cllai璐f．or exploring posterior distributions[J]．Ammls of Statistics，1994，22(4)： 1701．1762． [26】 HAARIO H，SAKShL～N E，TA^伽INEM J．An adaptive me仃opolis algorithm田．B锄oulli，200l，7(2)： 223．242． 万方数据