曹小群，等：对流-扩散方程源项识别反问题的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，方法具有收敛速度快、易于计算机实现和精度 高等特点。实验中对两个源项位置的估计精度达到 万方数据