732 军事医学2012年10月第36卷第10期Mil Med Sci,Vol36,No10,0ct,2012 GECCCCCCCEECCCCCe 基于MCMC方法的生物气溶胶袭击 风险管 理 施放源项参数反演 许晴，祖正虎，张文斗，徐致靖，黄培堂，郑涛 [摘要]生物气溶胶施放源项参数反演是生物气溶胶袭击危害评估的反问题，对危害评估及应急响应具有重要指 导意义。本文基于贝叶斯推理方法，利用生物传感器检测数据和正向大气扩散模型，构造似然函数，采用结合 Metropolis-Hasting算法的马尔可夫链蒙特卡洛(Markov chain Monte Carlo,MCMC)抽样，对施放源位置、高度、施放 剂量进行反演。统计分析表明，反演结果和初始源项参数设置吻合非常好，证明了方法的有效性。 [关健词]源项参数反演；贝叶斯推理：马尔可夫链蒙特卡洛；似然函数；生物安全 [中图分类号]E919 [文献标志码]A [文章编号]1674-9960(2012)10073204 Source inversion of bioaerosol attack based on MCMC method XU Qing,ZU Zheng-hu,ZHANG Wen-dou,XU Zhi-jing,HUANG Pei-tang,ZHENG Tao" (Institute of Biotechnology,Academy of Military Medical Sciences,Beijing 100071,China) Corresponding author,E-mail:zt19721@hotmail.com [Abstract]The inversion of bioaerosol release source parameters is the inverse problem of hazard assessment of bioaero- soattacks,and is of great significance for hazard assessment and emergency response.Based on observations of biosensor and concentrations predicted by an atmospheric dispersion model,a likelihood function was assigned,with which the Markov chain Monte Carlo(MCMC)sampling based on Bayesian inference was used to invert the source parameters,including the source location,source height,and dispersion strength.statistic analysis shows that the inversion results fit the initial source parameters very well.The validity of the method is proved. [Key words]source parameter inversion;Bayesian inference;Markov chain Monte Carlo;likelihood function;biosecurity 在城市上空施放气溶胶化生物剂是具有极大危 接反算法、优化法和随机逼近法3种。直接反算法 害性的生物恐怖袭击方式。在该类事件发生后，迅 主要是将大气扩散方程反解找到其数值或解析解， 速掌握施放源项参数如施放地点、高度、施放剂量等 早期反问题的研究主要集中于构造模型的控制微分 参数对于划定污染区域、评估人员伤亡、科学预测事 方程的反问题并直接求解。Tikhonov)的研究工作 态发展以及制定合理有效的应急响应策略具有十分 为这一领域做出了重要贡献。优化法的基本思想是 重要的意义。然而在实际中，该类恐怖事件具有隐 利用传感器测量值和数值模型模拟的信息，构造目 蔽性和突发性的特点，往往无法预知施放源的空间 标函数，通过优化的手段，选择最优的泄漏源项参 位置、施放剂量等关键信息，遭受生物恐布袭击的信 数，使其产生的模拟结果最大程度地与测量信息相 息通常是通过生物气溶胶侦查设备的检测数据或者 匹配，从而达到参数反演的目的，具体算法主要包括 大量非常规病例的暴发等后继现象才能获知。因 遗传算法、模拟退火算法等。Wang2I、McKinney 此，如何通过这些可获得的后验信息反演出施放源 等[3)在水文、地表水污染研究中采用遗传算法进行 关键参数对后续应急响应具有重要的指导意义。 模型数识别。采用概率方法研究源项参数反演问题 生物气溶胶施放源项参数反演是生物气溶胶袭 主要基于贝叶斯推断理论。该方法假设采样结果与 击危害评估的反问题。源项参数反演方法主要有直 泄漏源项参数的先验分布一致，观测的浓度值结合 [基金项目】国家自然科学基金资助项日(90924019)：“艾滋病和 贝叶斯推断可以获得施放源项参数的后验概率分 病卷性肝炎等重大传染病防治”国家科技重大专项(2008ZX10403) 布，通过预测模型参数计算出在这些参数下传感器 [作者简介】许晴，男，博士研究生，研究方向：生物恐饰事件危害 观测值出现的似然值，似然值反过来改进模型参数 评估与应急处置，Tel:010-63833523,E-mail:xugingl4443@yahoo.com.cn [作者单位].军事医学科学院生物工程研究所，北京100071 的估计，由于直接计算后验分布需要在整个参数空 [遁讯作者】]郑涛，Td:010-66948809,E-mail:zt19721@hotmail..com 间进行多维积分，在数值上难以实现，故使用马尔可 万方数据
732 军事医学2012年10月第36卷第10期Mil Med Sci，Vol 36，No 10，Oct，2012 基于MCMC方法的生物气溶胶袭击 施放源项参数反演 许晴，祖正虎，张文斗，徐致靖，黄培堂，郑涛 [摘要] 生物气溶胶施放源项参数反演是生物气溶胶袭击危害评估的反问题，对危害评估及应急响应具有重要指 导意义。本文基于贝叶斯推理方法，利用生物传感器检测数据和正向大气扩散模型，构造似然函数，采用结合 Metropolis-Hasting算法的马尔可夫链蒙特卡洛(Markov chain Monte Carlo，MCMC)抽样，对施放源位置、高度、施放 剂量进行反演。统计分析袁明，反演结果和初始源项参数设置吻合非常好，证明了方法的有效性。 【关键词] 源项参数反演；贝叶斯推理；马尔可夫链蒙特卡洛；似然函数；生物安全 [中图分类号] E919 [文献标志码] A [文章编号] 1674-9960(2012)10-0732-04 Source inversion of bioaerosol attack based on MCMC method XU口i愕，ZU Zheng-hu，ZHANG Wen—dou，XU Zhi-jing，HUANG Pei-tang，ZHENG Tao+ (Institute of Biotechnology，Academy of Military Medical Sciences，Beijing 100071，China) }Corresponding author，E-mail：ztl9721@hotmml．com [Abstract]The inversion of bioaerosol release source parameters is the inverse problem of hazard assessment of bioaerosol aRacks，and is of great significance for hazard assessment and emergency response．Based on observations of biosensors and．concentrations predicted by an atmospheric dispersion model，a likelihood function懈assigned，with which the Markov chain Monte Carlo(MCMC)sampling based on Bayesian inference Was used to invert the source parameters，including the source location，source height，and dispersion strength．statistic analysis shows that the inversion results fit the initial source parameters very well．The validity of the method is proved． [Key words]source parameter inversion；Bayesian inference；Markov chain Monte Carlo；likelihood function；biosecurity 在城市上空施放气溶胶化生物剂是具有极大危 害性的生物恐怖袭击方式。在该类事件发生后，迅 速掌握施放源项参数如施放地点、高度、施放剂量等 参数对于划定污染区域、评估人员伤亡、科学预测事 态发展以及制定合理有效的应急响应策略具有十分 重要的意义。然而在实际中，该类恐怖事件具有隐 蔽性和突发性的特点，往往无法预知施放源的空间 位置、施放剂量等关键信息，遭受生物恐怖袭击的信 息通常是通过生物气溶胶侦查设备的检测数据或者 大量非常规病例的暴发等后继现象才能获知。因 此，如何通过这些可获得的后验信息反演出施放源 关键参数对后续应急响应具有重要的指导意义。 生物气溶胶施放源项参数反演是生物气溶胶袭 击危害评估的反问题。源项参数反演方法主要有直 [基金项目]国家自然科学基金资助项目(90924019)；“艾滋病和 病毒性肝炎等重大传染病防治”国家科技重大专项(2008ZXl0403) [作者简介]许晴，男，博士研究生，研究方向：生物恐怖事件危害 评估与应急处置，Tel：010-63833523，E—mail：xuqin91443@yahoo．COIn．cn 【作者单位】．军事医学科学院生物工程研究所，北京100071 【通讯作者]郑涛，Tel：010-66948809，E-marl：ztl9721@hotmail．com 接反算法、优化法和随机逼近法3种。直接反算法 主要是将大气扩散方程反解找到其数值或解析解， 早期反问题的研究主要集中于构造模型的控制微分 方程的反问题并直接求解。Tikhonov…的研究工作 为这一领域做出了重要贡献。优化法的基本思想是 利用传感器测量值和数值模型模拟的信息，构造目 标函数，通过优化的手段，选择最优的泄漏源项参 数，使其产生的模拟结果最大程度地与测量信息相 匹配，从而达到参数反演的目的，具体算法主要包括 遗传算法、模拟退火算法等。Wang【2】、McKinney 等日1在水文、地表水污染研究中采用遗传算法进行 模型数识别。采用概率方法研究源项参数反演问题 主要基于贝叶斯推断理论。该方法假设采样结果与 泄漏源项参数的先验分布一致，观测的浓度值结合 贝叶斯推断可以获得施放源项参数的后验概率分 布，通过预测模型参数计算出在这些参数下传感器 观测值出现的似然值，似然值反过来改进模型参数 的估计，由于直接计算后验分布需要在整个参数空 间进行多维积分，在数值上难以实现，故使用马尔可 万方数据
军事医学2012年10月第36卷第10期Mil Med Sci,Vol36,No10,0et,2012 733 夫链蒙特卡洛(Markov Chain Monte Carlo,MCMC) 在确定似然概率时，需要考虑： 抽样方法得到参数的估计。 (1)传感器测量值与理论值的误差，Y:=T,+c; 利用MCMC方法求解反问题能方便地将各种 (2)数值模型预测值与理论值的误差，F:(X)= 先验信息融人到问题的求解过程中，减小了不确定 T,.(X)+e,通常这两种误差近似任务均服从 性，能对定义在高位空间且无明确数学表达式的概 Gauss分布，e:~Gau(0,o.i)和e:~Gau(0,o.i)。 率密度函数进行数值计算。而且该方法通常获得的 从而得到 是全局最可能解，而通常的最优化算法可能陷入目 标函数局部极小值。由于以上特点，Bayes方法被 p(Y.I7)exp(] 2d3i (2) 广泛应用于统计、生物、环境科学、雷达电磁波回波 反演等各个领域的反问题研究4】。本文采用该方 p(TF.(X)Xexp(X-F:(X) 2o (3) 法研究生物气溶胶施放源项参数反演问题研究。 这两种误差相互独立，可得 1方法 p(Y.IX)=[p(Y.ITX)p(Tx)dT exp 1.1贝叶斯推理 [E.(X)-Y (4) 贝叶斯推理是一种将先验信息转化为后验信息 2(3+) 的方法，其数学基础是贝叶斯定理。假设X代表模 各个生物传感器检测到的浓度值数据相互独 型参数，Y代表测量值，则模型参数的后验概率密度 立，可得似然概率 函数可以由贝叶斯定理计算得到 p(x)p(x)ep(- (5) 2(:+) p(xIY)=B(YIX)p(X) (1) p(Y) 1.3 MCMC抽样 对于生物气溶胶施放源项参数反演问题，可以 将式(5)带人式(1)，即可得到泄漏源项参数的 利用上述公式，式中， 后验概率密度。其中 X={X1,X2,4Xn}为施放源项参数，是需要进 p(=(x)p(xdF(-yE 2(c+o) p(X)dX 行反演的随机变量，包括施放点坐标(x,y)、施放点 高度(h)、施放剂量(s)。 (6) Y={Y,Y2,AY}为各生物传感器检测到的生 通过贝叶斯定理可以得到施放源项参数X(x, 物气溶胶浓度信息。 y,h,s)的后验概率分布，即为反问题的解。通常情 p(X)为施放源项参数X的先验概率密度函数， 况下，反问题的解后验概率分布只有在极简单的情 代表未获得观测值Y前参数X的分布规律，一般来 况下（例如，只有一个最大值等）才能用解析的方法 源于专家经验、主观判断等。 分析。直接利用式(1)计算后验概率分布，需要计 P(YIX)称为似然函数，表示在给定施放源项参 算分母的积分项（式6）。由于施放源项参数X为高 数X(x,y,h,s)情况下，正向模型数值结果拟合生物 维向量，对式(6)的计算将会产生维数灾难，即计算 传感器检测到的生物剂浓度值的程度，值越大表示拟 量将随着维数的增加而指数增加，会超出大部分计 合程度越好，反之越差，该函数的确定通常需要综合 算机的计算能力，因此必须采取一定的抽样方法回 考虑生物传感器的相关特性和正向数值模型的误差。 避该问题。 p(XIY)为X的后验概率密度函数，代表在获得 目前，抽样方法已经有很多种，如舍选法(rejec- 了生物传感器检测数据Y后，泄漏源项参数X的概 tion sampling),重要性抽样法(importance sampling）) 率分布。 和MCMC法。MCMC的计算量随维数线性增加，它 1.2似然函数的构造 不是在给定的点计算概率密度，而是根据概率密度 似然函数表示模型参数拟合观测数据的程度， 自动产生满足目标分布的样本。 因而它的构造是贝叶斯推理中重要的一环，直接影 本文采用MCMC抽样方法，直接产生一组目标 响到反演结果的可靠性、稳定性和计算效率。 分布为后验概率的随机抽样点，这些抽样点构成一 似然函数p(YIX)反映了生物传感器检测值 条Markov链，通过合理的构造转移概率，其收敛后 Y={Y,Y2,AY}与理论T={T1,T2,AT}的偏差。 的静态分布即为所需的后验分布。满足上述要求的 其中，理论值T可通过求解正向数值模型（如高斯烟 转移概率的MCMC抽样方法主要有Gibbs算法、 团模型或其他CFD模型等)获得，将正向数值模型的 Metropolis算法、Metropolis-Hasting算法等。本文采 计算结果记为F(X)={F(X),F2(X),AF.(X)}。 用结合Metropolis-Hasting算法的MCMC抽样，该算 万方数据
军事医学2012年10月第36卷第10期Mil Med Sci，Vol 36，No 10，Oct，2012 733 夫链蒙特卡洛(Markov Chain Monte Carlo，MCMC) 抽样方法得到参数的估计。 利用MCMC方法求解反问题能方便地将各种 先验信息融入到问题的求解过程中，减小了不确定 性，能对定义在高位空间且无明确数学表达式的概 率密度函数进行数值计算。而且该方法通常获得的 是全局最可能解，而通常的最优化算法可能陷入目 标函数局部极小值。由于以上特点，Bayes方法被 广泛应用于统计、生物、环境科学、雷达电磁波回波 反演等各个领域的反问题研究H圳。本文采用该方 法研究生物气溶胶施放源项参数反演问题研究。 1 方法 1．1 贝叶斯推理 贝叶斯推理是一种将先验信息转化为后验信息 的方法，其数学基础是贝叶斯定理。假设x代表模 型参数，y代表测量值，则模型参数的后验概率密度 函数可以由贝叶斯定理计算得到 p(Xl’，)：越掣 (1) pt J J 对于生物气溶胶施放源项参数反演问题，可以 利用上述公式，式中， X={x。，x：，以x。}为施放源项参数，是需要进 行反演的随机变量，包括施放点坐标(x，Y)、施放点 高度(h)、施放剂量(S)。 Y={yl，y2，A K}为各生物传感器检测到的生 物气溶胶浓度信息。 P(X)为施放源项参数x的先验概率密度函数， 代表未获得观测值y前参数x的分布规律，一般来 源于专家经验、主观判断等。 P(YIX)称为似然函数，表示在给定施放源项参 数X(戈，)，，h，s)情况下，正向模型数值结果拟合生物 传感器检测到的生物剂浓度值的程度，值越大表示拟 合程度越好，反之越差，该函数的确定通常需要综合 考虑生物传感器的相关特性和正向数值模型的误差。 P(XI y)为x的后验概率密度函数，代表在获得 了生物传感器检测数据y后，泄漏源项参数盖的概 率分布。 1．2似然函数的构造 似然函数表示模型参数拟合观测数据的程度， 因而它的构造是贝叶斯推理中重要的一环，直接影 响到反演结果的可靠性、稳定性和计算效率。 似然函数P(Yl X)反映了生物传感器检测值 Y={y。，y2，A k}与理论T={r。，疋，A丁。}的偏差。 其中，理论值T可通过求解正向数值模型(如高斯烟 团模型或其他CFD模型等)获得，将正向数值模型的 计算结果记为F(X)={F。(x)，疋(x)，以n(x)}。 在确定似然概率时，需要考虑： (1)传感器测量值与理论值的误差，y：=t…+占。； (2)数值模型预测值与理论值的误差，E(X)= t…(x)+ei，通常这两种误差近似任务均服从 Gauss分布，si～Gau(0，盯2，，。)和ei～Gau(O，矿2“)。 从而得到 p(tIt．。。，石)。e。p{一上兰_=窖掣}(2) ‘Cry．i p(I，。。IFi(X)，x)∞expl一上至塑生坚2掣}(3) Z4Tf,i 这两种误差相互独立，可得 ．p(KIx)=￡p(EIt，。。，x)p(t．，。Ix)dt．。。z exp{一 [t(x)一一]2l ／¨ 2(巧ii。)‘ ”’ 各个生物传感器检测到的浓度值数据相互独 立，可得似然概率 p(yIx)=玉(K Ix)*exp{一圭端} (5) 1．3 MCMC抽样 将式(5)带入式(1)，即可得到泄漏源项参数的 后验概率密度。其中 JP(1I)=fp(YIX)p(X)dⅨ*『唧{一量嬲错}P(x)ax (6) 通过贝叶斯定理可以得到施放源项参数X(戈， Y，h，s)的后验概率分布，即为反问题的解。通常情 况下，反问题的解后验概率分布只有在极简单的情 况下(例如，只有一个最大值等)才能用解析的方法 分析。直接利用式(1)计算后验概率分布，需要计 算分母的积分项(式6)。由于施放源项参数x为高 维向量，对式(6)的计算将会产生维数灾难，即计算 量将随着维数的增加而指数增加，会超出大部分计 算机的计算能力，因此必须采取一定的抽样方法回 避该问题。 目前，抽样方法已经有很多种，如舍选法(rejection sampling)，重要性抽样法(importance sampling) 和MCMC法。MCMC的计算量随维数线性增加，它 不是在给定的点计算概率密度，而是根据概率密度 自动产生满足目标分布的样本。 本文采用MCMC抽样方法，直接产生一组目标 分布为后验概率的随机抽样点，这些抽样点构成一 条Markov链，通过合理的构造转移概率，其收敛后 的静态分布即为所需的后验分布。满足上述要求的 转移概率的MCMC抽样方法主要有Gibbs算法、 Metropolis算法、Metropolis—Hasting算法等。本文采 用结合Metropolis-Hasting算法的MCMC抽样，该算 万方数据
734 军事医学2012年10月第36卷第10期Mil Med Sci,.Val36,No10,0ct,2012 法的基本思想是构造一个以目标分布π(x)为不变 对炭疽大气扩散应急响应的研究等。作为验证本文 分布的Markov链，根据π(x)和提议分布(proposal 所提出的方法有效性的数值试验示例，下文采用该 distribution)g(x·Ix),可以从当前位置x获得候选 模型计算炭疽菌随大气扩散后的地面浓度分布。采 抽样值x·,Markov链从位置x移到x·的接受概率 用其他正向数值模型不影响本文方法的正常使用。 为Ax,x)=min1,[ex4x]l。该算法 根据文献[11]中已有的工作，在上述想定情况下， p(x)q(x"Ix) 扩散后达到稳态的浓度分布如图2所示。9个生物 伪代码如下： 侦查设备分别检测出相应的浓度值，本文的工作就是 (1)Initializa x(0)/产生初始点 利用检测出的这些浓度值反演出施放源项参数X(x, (2)for i=0 to N-1 y,h,s),为后续危害评估及应急响应工作提供参考。 Sample x·~g(x·Ix(i))/根据提议分布在当 4001000200030 100050008000 前参数x()状态下计算获得新的参数x· 300 Sample u~U[0,1]/产生[0,1]之间的随机数 200 100 ifu()=min( p(x)q(x'Ix(i)) -100 -200 x(i+1)=x" //接受测试参数 -300 else -400 x(i+1)=x(i) /拒绝测试参数 图2正向数值模型计算结果 end 令X表示源项参数，Y表示生物传感器测量的 2 结果 浓度信息，F表示前向模型获得的传感器位置浓度 按照上文所述的技术路线，以实时在线的生物 理论值，T表示传感器处浓度的真实检测值。综合 传感器所监测到的气溶胶浓度稳态值为后验信息， 以上，本文的技术路线可由图1表示。 结合正向大气扩散模型，构造似然函数，采用结合 证向数值计第 Metropolis-Hasting算法的MCMC抽样，对施放源位 正向数值 PTF.X) 1T.0 生物传感器 置(x,y)、高度(h)、施放剂量(s)进行反演。整个抽 模型误差 贝叶斯推理 测量值 样的Markov链长取为10000，bum-in即过渡阶段取 施放源参数 传感器误差 后验概率分布 1000步。施放源各参数(x,y,h,s)的先验概率分布 IAXY 取均匀分布，x~U[0,6000],y~U[-400-400], MCMC采样 h~U[0,300],s~U[10°，104]，可以看出这是一个 采样点统计分析 相当大的状态空间，状态空间中每一个点都以一定 施放源参数估计 概率满足施放源项参数，而本文采用的方法就是从 这个庞大的状态空间中找出满足后验信息的施放源 图1技术路线 项参数最大概率值。 1.4正向数值模拟 从该状态空间中随机选取3个点作为Markov 不失一般性，本文假定所研究的地理范围为图 链搜索起始点，经初始bum-in阶段后，3条Markov链 2所示的区域（长度单位：m)。在该区域内分布着 很快收敛到真实泄漏源(0,0)的位置附近（图3）。 编号1~9共9个生物气溶胶侦查设备。想定如下炭 400 10002000 30 005600 疽芽孢杆菌气溶胶袭击，施放地点位于坐标(0,0)处， 300 施放点高度H=100m,施放剂量s=102芽孢，生物 200 气溶胶侦查设备检测到浓度信息时风速u=5m/s, 100 风向沿X轴正向，正向数值模型采用高斯烟团模 型。该模型比较适用于均匀、稳定的大气条件，大多 -100 应用于地面开阔平坦区域的大气扩散模拟，在诸多 炭疽菌大气扩散相关研究中被采用，取得了良好的 效果。如Meselson等9)对1979年前苏联斯维尔德 400 洛夫斯克炭疽泄漏事件的调查研究以及Wein等1o] 图3三条独立Markov链的搜索过程 万方数据
734 军事医学2012年10月第36卷第lO期Mil Med Sci，Vol 36，No 10，Oct，2012 法的基本思想是构造一个以目标分布7r(戈)为不变 分布的Markov链，根据7r(菇)和提议分布(proposal distribution)q(戈+I戈)，可以从当前位置戈获得候选 抽样值石+，Markov链从位置茗移到并+的接受概率 为A(算，x‘)=min I 1,[赵P乏黼q ]}。 该算法 L戈J L并 I戈J 伪代码如下： (1)Initializa x(0)／／产生初始点 (2)for i=0 to N一1 Sample x+～q(x+I戈(i))／／根据提议分布在当 前参数石(i)状态下计算获得新的参数戈+ Sample“～U[0，1]／／产生[0，1]之间的随机数 if u<A(戈(i)，戈+)=min{1，[气；i揣]} z(i+1)=彳+ ／／接受测试参数 else x(i+1)=戈(i) ／／拒绝测试参数 end 令x表示源项参数，l，表示生物传感器测量的 浓度信息，F表示前向模型获得的传感器位置浓度 理论值，丁表示传感器处浓度的真实检测值。综合 以上，本文的技术路线可由图l表示。 正向数值计算 模型误差l 莲璺鍪篁卜j卫也l贝叶斯推理 ’J“…Ⅻm5 施放源参数 后验概率分布 』P(X"1Y) MCMC采样 墨壁皇堕生坌塑l 巫函 』卜生盏墨霍器 L坐生 传感器误差I 1．4正向数值模拟 不失一般性，本文假定所研究的地理范围为图 2所示的区域(长度单位：m)。在该区域内分布着 编号1-9共9个生物气溶胶侦查设备。想定如下炭 疽芽孢杆菌气溶胶袭击，施放地点位于坐标(0，0)处， 施放点高度H=100 m，施放剂量s=1012芽孢，生物 气溶胶侦查设备检测到浓度信息时风速／／,=5 m／s， 风向沿x轴正向，正向数值模型采用高斯烟团模 型。该模型比较适用于均匀、稳定的大气条件，大多 应用于地面开阔平坦区域的大气扩散模拟，在诸多 炭疽菌大气扩散相关研究中被采用，取得了良好的 效果。如Meselson等∽1对1979年前苏联斯维尔德 洛夫斯克炭疽泄漏事件的调查研究以及Wein等¨叫 对炭疽大气扩散应急响应的研究等。作为验证本文 所提出的方法有效性的数值试验示例，下文采用该 模型计算炭疽菌随大气扩散后的地面浓度分布。采 用其他正向数值模型不影响本文方法的正常使用。 根据文献[11]中已有的工作，在上述想定情况下， 扩散后达到稳态的浓度分布如图2所示。9个生物 侦查设备分别检测出相应的浓度值，本文的工作就是 利用检测出的这些浓度值反演出施放源项参数X(戈， Y，h，S)，为后续危害评估及应急响应工作提供参考。 图2正向数值模型计算结果 2结果 2284—4568 4568—6852 6852—9316 9136—11420 1420—13705 13 705—1 5989 按照上文所述的技术路线，以实时在线的生物 传感器所监测到的气溶胶浓度稳态值为后验信息， 结合正向大气扩散模型，构造似然函数，采用结合 Metropolis-Hasting算法的MCMC抽样，对施放源位 置(戈，Y)、高度(h)、施放剂量(s)进行反演。整个抽 样的Markov链长取为10 000，bum—in即过渡阶段取 1000步。施放源各参数(戈，Y，h，s)的先验概率分布 取均匀分布，算～叭0，6000]，Y～U[一400—400]， h～U[0，300]，s～U[1010，1014]，可以看出这是一个 相当大的状态空间，状态空间中每一个点都以一定 概率满足施放源项参数，而本文采用的方法就是从 这个庞大的状态空间中找出满足后验信息的施放源 项参数最大概率值。 从该状态空间中随机选取3个点作为Markov 链搜索起始点，经初始burn—in阶段后，3条Markov链 很快收敛到真实泄漏源(0，0)的位置附近(图3)。 图3三条独立Markov链的搜索过程 万方数据
军事医学2012年10月第36卷第10期Mi1 Med Sci,Vol36,No10,0ct,2012 735 图4A~D分别展示了某次抽样过程中施放源 表1施放源项参数反演统计结果 项参数(x,y,h,s)的迭代曲线，通过MCMC和似然 源项参数真实值 先验区间 平均值 函数检验的样本占总样本数量的36.8%，此时样本 0 [0,60001 0.6 被接受。 y 0 [-400.400] -26 h 100 [0,300] 98.05 表1所示为施放源项参数反演统计量，可以看 1012 [100,1014 9.8×101 到反演出的结果与真实值吻合的非常好。 120- 500- 10 104 110E+012 400 88 102 105E+012 200 c100 M1.00E+012 100 0 98 9.50E+011 -2 0 WL0WY 40 96 9.00E+011 010002000300040005000 010002000300040005000 01000200030004000 2000 4000 迭代次数 迭代次数 迭代次数 迭代次数 B C D 图4 某次抽样过程中施放源项参数(x,y,h,s)的迭代曲线 3讨论 要更深入的进一步研究。 对有可能发生的在空气中大规模施放生物气溶 【参考文献】 胶形式的生物恐怖袭击，快速准确地判定施放源地 [1]Tikhonov AN.Arsenin VY.Solution of Ill-posed Problem[M] 点、施放点高度、施放剂量等关键信息对后续的危害 New-York John Wiley and Sons.1977. 评估及应急救援等工作具有重要意义。本文提出了 [2]Wang QJ.The genetic algorithm and its application to ealibrating 基于MCMC方法的生物气溶胶施放源项参数反演方 conceptual rainfall-runoff models[J].Water Resour Res,1991. 27(9):2467-2471. 法，结合外界环境信息、传感器探测到的稳态生物剂 [3]MeKinney DC,Min DL.Genetic algorithm solution of ground wa- 浓度以及正向大气扩散模型，反演施放源项参数(x, ter management models J].Water Resour Res,1994.30(6): y,h,s),反演结果和数值试验想定初始值符合较好， 1897-1906. 证实了该算法的有效性。必须说明的是，本文的研 [4]Gilks WR,Richardson S.Spiegelhalter DJ.Markov Chain Monte 究基于一定的理想假设，如采用简单的高斯烟团模 Carlo in Practice M].London:Chapman Hall/CRC,1996. 型计算炭疽菌气溶胶扩散后的稳态浓度值，施放源 [5]郭少东，杨锐，翁文国.基于MCMC方法的城区有毒气体扩 散源反演[J].清华大学学报.2009.49(5)：11-16. 项参数先验概率分布及先验区间的设定等。生物剂 [6]曹小群，宋君强.张卫民，等.MCMC方法在生物逆问题求解 大气扩散是个复杂的过程，涉及到生物剂的存活、衰 中的应用[J],计算机工程与应用.2011,47(26)：219-221. 亡、扩散、沉积等多种现象，扩散后的地面浓度受生 [7]盛峥，黄思训.曾国栋.利用Bayesian-MCMC方法从雷达回 物剂自身理化特性，天气因素如日光、温度、湿度、风 波反演海洋波导[J月，物理学报，2009,58(6)：4335-4341. 速、风向，下垫面情况等诸多因素的影响，而目前尚 [8]毛星，刘茂，刘冬华，等。应用贝叶斯理论的河流污染源 重建探讨[J】.安全与环境学报，2009,9(1)：100-103 没有能够精确描述生物剂随大气扩散过程的数学模 [9]Meselson M.Guillemin J,Hugh-Jones M.et al.The Sverdlovsk an 型，构造一个合理的且计算可行的生物剂大气扩散 thrax outbreak of 1979[J].Science.1994.266 (5188)1202 数学模型是研究如何有效应急处置此类生物恐怖事 1208. 件的一个重要研究课题。另外，本文所采用的方法 [10]Wein LM,Craft DL,Kaplan EH.Emergency response to anan- 适用于生物气溶胶扩散达到稳态后的源项参数反演 thrax attack[J].Proe Natl Acad Sci USA.2003,100(7):4346- 4351. 问题，而在实际情况中，在线生物气溶胶传感器会实 [11】许晴，祖正虎，张文斗，等，生物恐怖事件计算实验支持平 时探测到生物剂浓度信息，如何充分利用不同时刻 台及研究实例[J.军事医学，2011,35(11)：814-818 的传感器数据以更精确地反演施放源项参数，还需 (钱全编辑2012-08-09收稿) 万方数据
军事医学2012年10月第36卷第10期Mil Med Sci，Vol 36．No 10．Oct．2012 735 图4 A～D分别展示了某次抽样过程中施放源 项参数(x，Y，h，s)的迭代曲线，通过MCMC和似然 函数检验的样本占总样本数量的36．8％，此时样本 被接受。 表l所示为施放源项参数反演统计量，可以看 到反演出的结果与真实值吻合的非常好。 120 100 80 60 毛40 X 20 0 —20 —40 表l 施放源项参数反演统计结果 0 1 000 20003000 4000 5000 0 1 000 2000 3000 4000 迭代次数 迭代次数 B C 图4 某次抽样过程中施放源项参数(』，Y，h，s)的迭代曲线 3讨论 对有可能发生的在空气中大规模施放生物气溶 胶形式的生物恐怖袭击，快速准确地判定施放源地 点、施放点高度、施放剂量等关键信息对后续的危害 评估及应急救援等工作具有重要意义。本文提出了 基于MCMC方法的生物气溶胶施放源项参数反演方 法，结合外界环境信息、传感器探测到的稳态生物剂 浓度以及正向大气扩散模型，反演施放源项参数(戈， Y，h，s)，反演结果和数值试验想定初始值符合较好， 证实了该算法的有效性。必须说明的是，本文的研 究基于一定的理想假设，如采用简单的高斯烟团模 型计算炭疽菌气溶胶扩散后的稳态浓度值，施放源 项参数先验概率分布及先验区间的设定等。生物剂 大气扩散是个复杂的过程，涉及到生物剂的存活、衰 亡、扩散、沉积等多种现象，扩散后的地面浓度受生 物剂自身理化特性，天气因素如日光、温度、湿度、风 速、风向，下垫面情况等诸多因素的影响，而目前尚 没有能够精确描述生物剂随大气扩散过程的数学模 型，构造一个合理的且计算可行的生物剂大气扩散 数学模型是研究如何有效应急处置此类生物恐怖事 件的一个重要研究课题。另外，本文所采用的方法 适用于生物气溶胶扩散达到稳态后的源项参数反演 问题，而在实际情况中，在线生物气溶胶传感器会实 时探测到生物剂浓度信息，如何充分利用不同时刻 的传感器数据以更精确地反演施放源项参数，还需 1 10E+012 1 0fiE+012 ∞1 OOE+012 9 50E+011 9 OOE+011 要更深入的进一步研究。 【参考文献】 [1]Tikhonov AN，Arsenin VY．Solution of Ill—posed Problem[M]． New—York：John Wiley and Sons．1977 [2] Wang QJ．The genetic algorithm and its application to calibrating conceptual rainfall—runoff models[J]．Water Recur Res，1 99 1， 27(9)：2467—2471． [3] McKinney DC，Min DL．Genetic algorithm solution of ground wa— ter management models[J]．Water Resour Res．1994，30(6)： 1897一1906． [4] Gilks WR，Richardson S，Spiegelhaher DJ．Markov Chain Monte Carlo in Practice[M]．London：Chapman＆HalI／CRC，1 996． [5] 郭少东，杨锐，翁文国．基于MCMC方法的城区有毒气体扩 散源反演[J]．清华大学学报，2009．49(5)：ll—16． [6] 曹小群，宋君强．张卫民．等．MCMC方法在生物逆问题求解 中的应用[J]．计算机工程与应用，201 1，47(26)：219—221． [7] 盛 峥，黄思训，曾国栋．利用Bayesian．MCMC方法从雷达回 波反演海洋波导[J]．物理学报．2009，58(6)：4335—4341． [8] 毛星，刘茂，刘冬华．等．应用贝叶斯理论的河流污染源 重建探讨[J]．安全与环境学报，2009．9(1)：100—103． [9] Meselson M，Guillenfin J，Hugh-Jones M，el a1．The Sverdlovsk an— thr&x outbreak of 1979[J]．Science，1994，266(5188)：1202一 1208 [10] Wein LM，Craft DL，Kaplan EH．Emergency response to an an— thra=x attack[J]．Proc Natl Acad Sci USA．2003．100(7)：4346— 4351． [11] 许晴，祖正虎．张文斗．等．生物恐怖事件计算实验支持平 台及研究实例[J．军事医学．201 1，35(1 1)：814—818． (钱全编辑 2012一08—09收稿) ∞∞∞∞∞o 5 4 3 2， 一广仁一x 万方数据