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链的搜索过程 万方数据