482 仪器仪表学报 第33卷 法进行了研究。而在电容成像领域，目前只有少数学者 到介电常数的分布。 对基于统计方法的电容成像图像重构方法进行了研究， 如Watzenig及Fox02009年在对基于统计方法的ET图 3 Markov chain Monte Carlo重构算法 像重构方法进行综述的基上对其在ECT图像重构中的 应用进行了初步研究。 3.1MCMC模型 本文重点对利用MCMC方法来进行ECT图像重构 在利用Bayesian推理对参数进行估计的时候，如果 的相关问题进行了研究，重点讨论了ECT图像重构的 后验分布是复杂的、高维的分布，则以利用马尔科夫链 MCMC实现方法，为ECT图像重构提供了一种新的解决 蒙特卡罗(MCMC)的方法对后验概*进行采样。MCMC 思路。 总能得到一条或几条收敛的马尔科夫链，该马尔科夫链 的极限分布就是所需的后验分布。 2 ECT传感器 假设在获得电容测量数据的过程中引入了一个与其 独立的噪声，该噪声具有高斯分布，均值为0，标准偏 ECT传感器如图1所示，12个电极均匀地排列在 差为d。 管道外壁处，电极外部有屏蔽罩，可以防止外界的电磁 入=Sg+p (5) 干扰。测量过程中，E1首先作为激励电极，其他的电极 式中：p是高斯随机变量，p~N(0,d)。本文中，取d= 接地，测其E1和其他电极之的电容。然后将E2作 0.001。 为激励电极，同样其他的电极接地，测量E2和除E1外 介电常数分布g的似然函数与exp(L(g))成正比。 其他的电极间的电容。重复这个步骤，直至E11与E12 其中，L(g)是g的对数似然函数，可以表示为： 间的电容测量完毕，则测量得到的总的独立电容个 L(g)=-(I Sg -Al+ulg)/2d (6) 数为： M=N(N-1)/2 式中：“是正则化参数，引入此参数，可以提高重建图像 (1) 的质其。在本文中，u初设为0.001。 式中：N是电极的数目，在此例中N=12,M=66。 参考Nicholls和Fox的T作，以Markov随机场的 屏蔽罩 管道外壁 .10E9 方式给出介电常数分布g的先验概率。此种方式给出的 E11 先验概率，认为每个像素点的取值只和与其直接相邻的 12 E7 4个格子的像素值有关，具体可以由式(7)来描述： F6 Pr(g)cexp(∑Hn(g) (7) 3 E4 电板 式中：.(g）=二8表示对于标号为m的格子，与 图112电极ECT传感器示意图 其相邻的、且与其介电常数相同的格子的总数。6。.。表示 Fig.1 A 12-electrode ECT sensor 当a=b时，其取值为1。 在此，定义一个随机变量为X,的Markov链{X,}。, 图！所示的传感器，电极的张角与两个电极之间夹 当进行t步选代后，Pr(X=g,IX。=go)趋近于当t+ 角的比值为9：1，管壁的相对介电常数为2.6。 o时的后验概卒分布P,(gIA)。Markov链的转移核利 在ECT系统中，测敏得到的电容值可以看成管道中 用Metropolis-Hastings方法得到，如(8)所示： 介电常数的函数： Pr (X=giI X:=g:)=qp(X=gX, C=(8) (2) g）·a,(X+1=g1IX=g) (8) 对于一个介电常数的微小变化而言，相应的电容变 式中：9。是产生此步转移的概率，是一个潜在的转移核， 化可以用式(3)的离散形式来描述： a,为接受此步转移的概率，也就是说，对于一个新的状 AC=d1 (3) 态，以a,的概率接受转移到X1,以1-，的概率拒绝 式中：J是Jacobi矩阵（灵敏度矩阵），L是待重构图像的 转移到X1。参考Nicholls和Fox1的工作，在本文中使 像素个数。与式(3)对应的归一化形式可表示为： 用到的状态转移包括以下4种形式。 A=Sg (4) Movel:一个像素发生改变。在所有的像素中，以均 式中：入是归一化电容向其，S是归一化灵敏度矩阵，g 匀分布随机选择一个像素，并为其指定一个新的像索值。 是归一化介电常数向量。重构算法的任务就是利用测址 9(g→g)=1/M(C-1) 得到的电容数据和预先计算的敏感矩阵求解式(4)来得 a,(g-+g')=minl,e2".-.》×e)-{ 万方数据仪器仪表学报 第3 3卷 法进行了研究。而在电容成像领域，目前只有少数学者 对基于统计方法的电容成像图像鼋构方法进行了研究， 如Watzenig及Fox瑚12009年在对基于统计方法的EIT图 像重构方法进行综述的基上对其在ECT图像重构中的 应用进行了初步研究。 本文重点对利用MCMC方法来进行ECT图像蘑构 的相关问题进行了研究，重点讨论r ECT图像重构的 MCMC实现方法，为ECT图像莺构提供r一种新的解决 思路。 2 ECT传感器 删。◆ 图l所示的传感器，电极的张角与两个电极之问央 角的比值为9：1，管壁的相对介电常数为2．6。 在ECT系统中，测肇得到的电容值叮以看成管道中 介电常数的函数： C=孝(占) (2) 对于一个介电常数的微小变化而言，相应的电容变 化可以用式(3)的离散形式来描述： AC=J Ae (3) W¨ 村x￡￡×I 式中：J是Jaeobi矩阵(灵敏度矩阵)，￡是待重构图像的 像素个数。与式(3)对应的归一化形式可表示为： A=sg (4) 式中：A是归一化电容向量，S是归一化灵敏度矩阵，g 是归一化介电常数向量。重构算法的任务就是利用测量 得到的电容数据和预先计算的敏感矩阵求解式(4)来得 到介电常数的分布。 3 Markov chain Monte Carlo重构算法 3．1 MCMC模型 在利用Bayesian推理对参数进行估计的时候，如果 后验分布是复杂的、高维的分布，则町以利用马尔科夫链 蒙特膏罗(MCMC)的方法对后验概率进行采样。MCMC 总能得到一条或几条收敛的马尔科夫链，该马尔科夫链 的极限分布就是所需的后验分布。 假设在获得电容测量数据的过程中引入了一个与其 独讧的噪声p，该噪声具有高斯分布，均值为0，标准偏 差为d。 A=sg+p (5) 式中：9是高斯随机变量，垆一N(o，d2)。本文中，取d= 0．001。 介电常数分布g的似然函数与exp(￡(暑))成正比。 其中，Lig)是g的对数似然函数，可以表示为： ￡(g)=一(I Sg—A l+芦}g 1)／2d2 (6) 式中：弘是正则化参数，引入此参数，可以提高重建图像 的质麓。在本文中，肛初设为0．001。 参考Nieholls和Fox‘is]的：【作，以Markov随机场的 方式给出介电常数分布g的先验概率。此种方式给出的 先验概率，认为每个像素点的取值只和与其直接相邻的 4个格子的像素值有关，具体可以由式(7)来描述： 肼 pr(g)oc exp(∑以(g)) (7) 式中：以(g)=∑以．．。，表示对于标号为m的格子，与 其相邻的、且与其介电常数相同的格子的总数。艿“表示 当口=b时，其取值为1。 在此，定义一个随机变量为x。的Markov链{x。}；， 当进行t步迭代后，Pr(置=g。f x0=go)趋近于当t一 ∞时的后验概率分布Pr(g A)。Markov链的转移核利 用Metropolis—Hastings方法得到，如(8)所示： Pr”’(X川=g川I X。=g。)=qp(X¨I=g川I X。= g。)‘瑾。(Xl+l=g川l X。=g，) (8) 式中：q。是产生此步转移的概率，是一个潜在的转移核， a。为接受此步转移的概率，也就是说，对于一个新的状 态，以％的概率接受转移到．)if川，以1一叱的概率拒绝 转移到x川。参考Nicholls和Fox¨引的工作，在本文中使 用到的状态转移包括以下4种形式。 Morel：一个像素发生改变。在所有的像素中，以均 匀分布随机选择一个像素，并为其指定一个新的像素值。 ql(g—g’)=1／M(C一1) al(g—+譬’)=min{I，e"‘辩-‘r’一H·‘童’’x eL(Z7)一“D i 万方数据