（应用数学专业论文）基于贝叶斯理论的电容层析成像算法研究.pdf

丫节 o r e c o n s t r u c t i o na l g o r i t h mi ne l e c t r i c a lc a p a c i t a n c e t o m o g r a p h yb yb a y e s i a ni n f e r e n c e at h e s i ss u b m i t t e dt o d a l i a nm a r i t i m eu n i v e r s i t y i np a r t i a lf u l f i l l m e n to ft h er e q u i r e m e n t sf o rt h ed e g r e eo f m a s t e ro fs c i e n c e b y z h a n g z h i l i n ( a p p l i e dm a t h e m a t i c s ) j t h e s i ss u p e r v i s o r ：a s s o c i a t ep r o f e s s o rf uh o n g s u n j u n e2 0 1 1 、，r 了心 丫一 b q 气 _ 譬 懒18 9 5 9 9 9 大连海事大学学位论文原创性声明和使用授权说明 原创性声明 本人郑重声明：本论文是在导师的指导下，独立进行研究工作所取得的成果， 撰写成硕士学位论文= ! 基王巫吐堑堡途的鱼空屋堑盛堡篡选班究：。除论文中已经 注明引用的内容外，对论文的研究做出重要贡献的个人和集体，均已在文中以明确 方式标明。本论文中不包含任何未加明确注明的其他个人或集体已经公开发表或 未公开发表的成果。本声明的法律责任由本人承担。 学位论文作者签名：自盘查拯 学位论文版权使用授权书 本学位论文作者及指导教师完全了解大连海事大学有关保留、使用研究生学 位论文的规定，即：大连海事大学有权保留并向国家有关部门或机构送交学位论文 的复印件和电子版，允许论文被查阅和借阅。本人授权大连海事大学可以将本学位 论文的全部或部分内容编入有关数据库进行检索，也可采用影印、缩印或扫描等复 制手段保存和汇编学位论文。同意将本学位论文收录到中国优秀博硕士学位论 文全文数据库( 中国学术期刊( 光盘版) 电子杂志社) 、中国学位论文全文数据 库( 中国科学技术信息研究所) 等数据库中，并以电子出版物形式出版发行和提供 信息服务。保密的论文在解密后遵守此规定。 本学位论文属于：保密口在年解密后适用本授权书。 不保密( 请在以上方框内打“ ) 珠琳名：锄彳 日期：讼j 1 年多月必日 y 心 丫咚 、 t q 气 飞， 中文摘要 摘要 电容层析成像( e c t ) 是一种可获取封闭区域内介质分布图像的过程成像技 术，它具有非侵入性、适用范围广、价格低廉等优点，因此在多种工业生产过程 中得到广泛运用。 传统的确定性算法获得的是介质分布的近似最优解，而基于贝叶斯的抽样方 法可以得到场域内介质分布的全面描述。 本文从贝叶斯理论出发，构造了一种新的e c t 图像重建算法。首先，在正问 题部分，通过自己编写的有限元仿真软件，进行有限元的计算，得到场域内介质 电位分布情况。其次，在反问题部分，以标准f 态分布作为先验信息，根据贝叶 斯理论，获得了e c t 介电常数的后验分布模型。最后，引入线性反投影的结果作 为马尔可夫链的初状态，应用m e t r o p o l i s h a s t i n g s 抽样算法对后验分布进行抽样， 根据抽样获得的样本期望值可重建图像，并获得了不同建议分布方差的成像结果。 关键词：电容层析成像；后验分布；m e t r o p o l i s - h a s t i n g s 算法；有限元法；正反 问题 电 l y 吣 气 英文摘要 a b s t r a c t e l e c t r i c a lc a p a c i t a n c et o m o g r a p h y ( e c t ) i sap r o m i s i n gt e c h n i q u eo fp r o c e s s t o m o g r a p h y ，w h i c hc a no b t a i nt h ed i s t r i b u t i o no ft h em e d i ai m a g ei na ne n c l o s e da r e a i t h a st h ea d v a n t a g eo fb e i n gn o n - i n t r u s i v e ，b r o a da p p l i c a t i o n ，a n dl o wp r i c e ，e t c ，s oi th a s b e e nw i d e l yu s e di nav a r i e t yo fi n d u s t r i a lp r o c e s s e s t r a d i t i o n a le c ti m a g er e c o n s t r u c t i o na l g o r i t h mc a no n l yo b t a i nt h ea p p r o x i m a t e o p t i m a ls o l u t i o no ft h ed i e l e c t r i cc o n s t a n t ，b u tt h es a m p l i n gm e t h o db a s e d o nb a y e s i a n g a m e sc a nb ed i s t r i b u t e dw i t h i nac o m p r e h e n s i v ed e s c r i p t i o no f t h ed i e l e c t r i cc o n s t a n t i nt h i sp a p e r , b a s e do nb a y e s i a nt h e o r y , w ec o n s t r u c tan e wi m a g er e c o n s t r u c t i o n a l g o r i t h mf o re c t f i r s t l y , i nt h ef o r w o r dp r o b l e m ，w eo b t a i nt h ed i s t r i b u t i o no ft h e p o t e n t i a lf i e l dw i t h i nt h em e d i ab yf i n i t ee l e m e n tc a l c u l a t i o n s s e c o n d l y , i nt h ei n v e r s e p r o b l e m ，b a s e d o n b a y e s i a nt h e o r y a n d p r i o r i i n f o r m a t i o no fs t a n d a r dn o r m a l d i s t r i b u t i o n ，w eo b t a i nt h ep o s t e r i o rd i s t r i b u t i o nm o d e lo ft h ed i e l e c t r i cc o n s t a n to f e c t f i n a l l y , w ei n t r o d u c et h el i n e a rb a c k - p r o j e c t i o na st h ei n i t i a l s t a t eo fm a r k o v c h a i n ，s a m p l et h ep o s t e r i o rd i s t r i b u t i o nw i t h t h ea p p l i c a t i o no fm e t r o p o l i s - h a s t i n g s s a m p l i n ga l g o r i t h m ，a n dr e c o n s t r u c tt h ei m a g eb ye x p e c t a t i o n sb a s e do ns a m p l e s t h e n w eo b t a i nv a r i a n c ei m a g i n gr e s u l t so fd i f f e r e n tp r o p o s a ld i s t r i b u t i o n k e yw o r d s ：e l e c t r i c a lc a p a c i t a n c et o m o g r a p h y ；p o s t e r i o r d i s t r i b u t i o n ； m e t r o p o l i s - h a s t i n g sa l g o r i t h m ；f i n i t ee l e m e n tm e t h o d ；t h ef o r w a r dp r o b l e m s a n d i n v e r s ep r o b l e m 、， 妒 f - v，“、 t k 1 1 i q i 伶 目录 目录 第1 章绪论1 1 1e c t 概述1 1 2 电容层析成像算法研究进展2 1 3 本文的主要结构和工作4 1 4 本文的创新点5 第2 章e c t 系统及其工作原理6 2 1e c t 成像系统结构和测量原理6 2 2e c t 数学模型7 2 3e c t 正问题9 2 3 1 正问题求解方法9 2 3 2 有限元法1 0 2 4e c t 逆问题1 8 2 4 1e c t 基本方程1 8 2 4 2 灵敏度系数矩阵的建立1 9 2 4 3e c t 反问题的不适定性2 0 2 4 4 图像重建算法2 1 2 4 5 线性反投影( l b p ) 算法2 1 2 5 本章小结2 2 第3 章贝叶斯理论和e c t 反问题2 3 3 1 贝叶斯理论2 3 3 1 1 贝叶斯定理2 3 3 1 2 先验信息的选择2 4 3 1 3 似然函数的选择2 5 3 1 4 后验概率密度函数的抽样方法2 6 3 2e c t 后验分布模型2 6 3 3 本章小结2 7 第4 章m h 抽样在e c t 图像重建中的应用2 8 4 1 m c m c 统计抽样2 8 4 1 1m o n t ec a r l o 概述2 8 4 1 2 舍选法抽样2 9 4 1 3重要性抽样2 9 4 1 4m o n t ec a r l o 算法过程3 0 4 2m h 抽样3 0 4 2 1m h 抽样的基本原理3 0 4 2 2m h 抽样的过程3 l 4 2 3 m h 抽样的数值例子及结果分析3 2 4 3基于m h 抽样的e c t 图像重建3 3 4 3 1介电常数边际分布3 3 4 3 2 m e t r o p o l i s h a s t i n g s 抽样3 4 4 3 3 样本的获取与处理3 5 4 3 4 仿真结果对比分析取l b p 投影结果为初状态3 5 第5 章本文的主要工作与展望3 8 5 1 主要工作与结论分析3 8 5 2 展望3 8 参考文献4 0 致谢4 4 j 一 p t ， k 协 嵋 基于贝叶斯理论的电容层析成像算法研究 第1 章绪论 1 1e c t 概述 由于计算机、电子技术、近代物理技术与信号处理技术的快速发展，以微波、 红外、激光与声发射技为特征的新型检测技术得到了长足进步，使传统以工程量 ( 如温度、压力、流量、流速、液压、料位等) 为主的检测技术，逐渐过渡到热 工量、电学量、力学量和光学量等多种物理量的近代物理检测技术。 进入2 1 世纪，随着信息时代的来临，科学技术的发展更为迅速，测试技术作 为人们认识客观世界的一个重要手段也显得更加现代化。不但工业过程要依靠各 种先进的测量方法来实现自动控制，各种工业产品资料要通过先进的在线监测来 保证，工业生产设备的安全运行要通过在线检测来实现；而且各种更为复杂的科 学实验要通过测量提供可靠的数据。 近代物理检测技术的发展呈现以下特点： 1 、传统的测温、测速方法都是接触式的，而近代的激光测速、激光测温都是 非接触式的，这种非接触式的测量方法，避免了传感器对被侧物理量场的干扰， 代表了当今测量技术的发展方向。 2 、在测量的时间域上，有物理量的静态测量发展为物理量的动态测量。 3 、在数据处理上，由被测数据的手工采集或仪表记录打着到计算机采集、存 储与处理。 2 0 世纪8 0 年代，由英国曼彻斯特大学提出了一种类似于医学的c t 技术的电 容层析成像技术( e l e c t r i c a lc a p a c i t a n c et o m o g r a p h y ，简记e c t ) ，这种方法根据 被测场域内不同介质介电常数的不同，当内部介质分布发生变化时，由边界电极 对间测量的电容值，运用相关的图像重建算法，就可重建出被测场域内的介质分 布图。 电容层析成像( e c t ) 是一种可获取封闭区域内介质分布图像的过程成像技 术，它具有非侵入性、适用范围广等优点，因此在多种工业生产过程中得到广泛 运用。并且它的操作简单和价格低廉比较适合我国的国情。但现有的相关检测技 第1 章绪论 术和设备还不能满足工业飞速发展的需求，因此研究这项检测技术对我国相关工 业的发展具有重大意义。 e c t 包括两个计算过程：即正演过程和反演过程。其中正演过程由已知场域 内介电常数分布情况求解边界电极对间的电容值；而相应的反演过程由已测得的 边界电极对间的电容值估计被测场域内介质分布情况。 1 2 电容层析成像算法研究进展 h o u s e f i e l d 在1 9 7 2 年研制出头颅x 射线断层成像装置，即c t 机，并获得了 1 9 7 9 年诺贝尔医学奖，其原理是：在不损伤物体本身内部结构的情况下，发射可 通过物体的讯号( 如射线、波、粒子、电磁场等) ，然后通过从体外接收到的信号， 利用数学方法和计算机处理后，反演出物体内部结构的分布信息，亦即图像反演 重建。e c t 的图像重建原理，基本类似于医学c t 的思想，但是c t 射线成像的性 质与e c t 电学成像不同，现有的医学c t 机器主要是通过发射x 射线，而且操作复 杂，成本昂贵。 e c t 技术作为一种非侵入的对封闭区域可视化的成像技术，其操作简单性为 工业生产的实时监测提供了可能，但目前此项技术在国内研究的并不多。从上个 世纪8 0 年代开始，很多科研机构和国家都对e c t 技术进行了卓有成效的研究。 表1 1 国内外相关研究机构及算法 t a b 1 1d o m e s t i ca n df o r e i g nr e s e a r c hi n s t i t u t i o n sa n da 培o d t h m s 国家研究机构 算法 中国科学院线性反投影算法( l b p ) ，l a n d w e b e r 迭代算法 中国 清华大学 水平集法，最速下降法，贝叶斯迭代方法 浙江大学 同步迭代法( s i r t ) ，共扼梯度算法 德国汉诺威人学 l b p ， s i r t ，迭代的代数法( a r t ) 墨西哥墨西哥石油研究所 l b p 挪威卑尔根大学基于模型的迭代法 英国曼切斯特大学理下学院l b p ，l a n d w e b e r 迭代算法，正则化 吣 v ， k 协 口 b 基于贝叶斯理论的电容层析成像算法研究 基于e c t 发展现状，一般将e c t 反演算法分为：确定性方法【1 7 】和概率统计方 法【件1 8 】。其中确定性方法又分为迭代算法和非迭代单步算法。 单步非迭代算法主要包括：t i k h o n o v 正则化算法、线性反投影算法( l b p ) 、 基于正则化和多元线性回归的图像重建算法( m l r r ) 、基于奇异值分解直接法 ( s v d ) 等单步算法，单步算法成像速度快，而且计算量小，但成像的精度较差。 迭代算法主要包括：l a n d w e b e r - f r i d m a n 迭代法、最速下降法、截断共扼梯度算 法、n e w t o n c g 算法、迭代的代数法( a r t ) 、水平集方法、同步迭代法( s i r t ) 等，迭代算法成像速度慢，但是成像质量相对较好。 概率统计方法主要包括基于贝叶斯理论的相关蒙特卡洛抽样方法，包括 g i b b s ，模拟退火方法、遗传算法、神经网络方法等。 国内方面，1 9 9 6 年，天津商学院的马军海等用最速下降法解决e i t 反问题【1 1 。 1 9 9 9 年，清华大学的宫莲等用广义线性增量函数理论和奇异值分解法解反问题1 2 】。 清华大学江鹏等( 2 0 0 7 ) 采用基于g a u s s i a n 窗函数的正则化方法【3 1 。天津大学的 王超等【4 】( 2 0 0 7 ) 将贝叶斯理论引入e c t 反演问题的中，运用g i b b s 统计抽样方法， 得到了不错的反演效果，但文中选择的后验边际分布相对简单。清华大学江鹏等 ( 2 0 0 8 ) 5 1 。采用基于贝叶斯理论的迭代算法，但是只针对了噪声和介质分布先验 概率相同的情况。西安交通大学胡红利等1 7 1 ( 2 0 1 0 ) ，利用r b f 神经网络方法进行 反演成像，对环流和核心流都得到了不错的效果。 随着计算机技术的不断发展，国内在传感器系统研究方面也得到了大量成果。 第四军医大学郑万松等( 2 0 0 6 ) 7 1 改进了腹部电极测量系统，提高了电极测量精度。 重庆大学罗辞勇等( 2 0 0 9 ) 【8 】采用电极交叉测量模式，提高了系统的抗噪声性能。 天津大学徐苓安教授和清华大学张宝芬教授成功研制了电阻层析成像系统原理性 样机【9 1 ，并在实验室中获得了不错的成像效果。 国际方面，x i e ( 1 9 9 2 ) 1 1 0 】等将反投影算法运用到e i t 图像重建中；i s a k s e n ( 1 9 9 6 ) 1 1 】运用基于优化理论模型的重建算法，对环状流和层状流有较好的反演 效果；n o o r a l a h i y a n ( 1 9 9 7 ) 1 1 2 j 等提出了神经网络方法在e c t 反演中的应用；曼彻 斯特理工学院的y a n g 等( 1 9 9 9 ) 1 1 3 提出在e c t 问题中使用l a n d w e b e r 迭代算法。 第1 章绪论 进入2 1 世纪，随着计算机技术的飞速发展，需要进行大量统计计算的抽样的方法 越来越多的应用到工程中，而将统计算法应用到e c t 反演问题方面，比较有代表 性的是奥地利格雷兹大学的d a n i e lw a t z e n i n g 、东芬兰大学的j a r ipk a i p i o 、新西兰 奥塔哥大学的c o l i nf o x 。 由于应用图像重建概率统计方法可以获得场域内介质分布的大量统计信息， 对于e c t 的后续分析问题有很大帮助，因此本文将基于贝叶斯理论的 m e t r o p o l i s - h a s t i n g s 算法运用到e c t 反演成像中，并对其进行了改进，获得了不错 的反演结果。 1 3 本文的主要结构和工作 本篇论文共分五章： 第一章绪论，简要介绍了层析成像技术的原理【1 2 1 4 l ，电容层析成像的发展历 程及研究意义； 第二章e c t 系统及工作原理，介绍了e c t 的系统构成及测量原理，分析了 e c t 的正问题和e c t 的逆问题的相关研究方法。并建立e c t 的正向模型，通过自 己编写的有限元程序进行仿真计算，将线性反投影算法引入到反演过程中，并将 其反演结果作为马尔可夫链的初值，避免落入局部极值的陷阱。 第三章贝叶斯理论和e c t 反问题，介绍了贝叶斯的相关原理【1 5 - 1 7 1 ，针对e c t 问题选取合适的j 下态分布先验概率，和高斯白噪声的测量误差，建立了介电常数 的后验分布模型。 第四章介绍了蒙特卡洛的相关抽样方法，分析了马尔可夫蒙特卡洛( m c m c ) 抽样的算法过程1 1 8 锄】，利用m c m c 算法中最经典的m h 算法进行反演成像，并 对结果进行了相关的分析和处理。提出了利用线性反投影的结果作为马尔可夫链 的初状态的方法，并可获得相对较优的重建效果。 第五章结论分析与建议。总结本文基于贝叶斯方法的e c t 图像重建算法的研 究成果，提出了一些建议与展望。 1 ： 一 、| 、， 。 k k h 嘣 ( 基于贝叶斯理论的电容层析成像算法研究 1 4 本文的创新点 作者依据贝叶斯理论以及有限元方法，构造了以标准的j 下态分布作为先验信 息的e c t 后验模型。 基于贝叶斯理论的m h 抽样方法处理e c t 反演问题，对不同方差给反演结果 的影响进行了研究，并获得了相对较优的反演成像结果，从数值仿真结果表明了 该方法的有效性。 利用以l b p 成像结果为初状态的m h 抽样算法实现了e c t 图像重建，避免了 马尔可夫链的局部极值，加快了收敛，获得了介电常数的大量后验分布统计信息， 为进一步的研究e c t 图像重建统计算法奠定了基础。 第2 章e c t 系统及其工作原理 第2 章e c t 系统及其工作原理 2 1e c t 成像系统结构和测量原理 图2 - 1 表示的是一个十二电极的e c t 系统结构乜卜矧，硬件部分主要由阵列电 极、图像重建与分析的计算机和数据采集、主控制系统三部分构成。首先传感器 将管壁内的多相流介质电容分布转换成输出电信号；通过数据采集与信号处理单 元将电信号转化为数字量，最后计算机则根据图像重建算法，重建出被测场域内 的介质分布情况并在相应的计算机上显示。 图2 1e c t 系统结构 f i g 2 1e c ts y s t e ms t r u c t u r e 以8 电极系统为例，从电极1 开始施加激励电压，电极2 到8 为测量电极， 可测得7 个测量数据；再对电极2 施加激励电压，电极l 接地，电极3 8 作为测量 电极，可测得6 个测量数据。依此方法下去，直到电极7 被激励，电极8 作为测 量电极，这个过程共测得2 8 个独立的电容值【3 4 1 。对于e c t 系统的n 个电极，设 其测量次数为m ，则有下式成立 m ：n ( n - 1 ) ( 2 1 ) 2 由于电容测量数据的好坏会直接影响到图像重建的质量和效率，因此，如何 提高正演精度，这也是本文研究的一个重点内容。 基于贝叶斯理论的电容层析成像算法研究 2 2e c t 数学模型 世间的事物或现象之间往往存在着一定的自然顺序，如时间顺序、空间顺序、 因果顺序，等等。所谓正问题( f o r w a r dp r o b l e m ) ，一般是按着一个物理系统被完 全描述的情况下来预测它的某些观测结果。反问题( i n v e r s ep r o b l e m ) 则是根据事 物的演化结果，由可观测的现象和一些测量结果来探求事物的内部规律或所受的 外部影响，由表及里，索隐探秘，起着倒果求因的作用【5 7 】。 图2 2 正问题和反问题 f i g 2 2f o r w a r dp r o b l e ma n di n v e r s ep r o b l e m ( d ) 图2 - 2 ( a ) 描述了正问题的发生过程；图2 - 2 ( b ) ( c ) ( d ) 描述了反问题的发生过 程 在研究正问题之前，根据被测场域的特性建立相关数学模型。对于e c t 问题 而言，我们假定被测场域的电磁场是静电场，根据电磁场相关理论，任何电磁场 都可用如下的麦可斯韦方程组描述【2 5 2 6 】 第2 章e c t 系统及其工作原理 v e ：一a b 0 t v b _ i - ， ( 2 2 ) v d = p 、7 v 日= 以+ 詈 其中，v 为矢量微分算子( 即泊松算子) ，h 为磁场强度，b 为磁感应强度，p 为 电荷密度，d 为电位移，e 为电场强度，以为电流密度矢量。 由静电场的相关特性可得各上述物理量之间的关系，另外电场和磁场间各物 理量还存在着如下关系【2 7 】 d = s e 蔓篓， 【e = 一v 其中，g 、孝、y 、矽分别表示介电常数、磁导率、电导率、场域内电势分布函数。 对于实际的e c t 问题而言，传感器在施加激励电压或电流后构成一静电场， 因此它可用如下方程来描述【3 5 】 v d = v 6 e = 刃e + e v s = p ( 2 4 ) 对于均匀介质，v s = 0 ，则 刃e = 一刃v 矽= p ，( 2 5 ) 得到如下p o s s i o n 方程 v 2 矽= - p 6 ( 2 6 ) 若场中无自由体电荷( 即p = 0 ) ，则有如下l a p l a c e 方程成立 v 2 矽= 0 ，( 2 7 ) 其中，v 2 为l a p l a c e 算子。 在实际的应用中，e c t 正问题属于第一类边值问题( 即d i r i c h l e t 边界条件) p 4 1 ，其a n * 边界条件后可用以下的方程表达 v s ( 五少) v ( 五y ) 】= 0 ( 2 1 1 ) 基于贝叶斯理论的电容层析成像算法研究 基于以上推导，当电极( i = l ，2 ，1 1 ) 依次作为激励电极时，可用静电场的数学 模型表示如下： v s ( x ，y ) v c i ( x ，j ，) 】= 0 ， 矽7 = 圪 x x ，, y y ) ) cf 10 f 。( j | ：l ，2 ，1 2 ，k f ) 且( 五y ) s ( r ，+ f 昭) ，( 2 1 2 ) 7 【 ( x ，y ) i ( j | = l ，2 ， ， f ) 且( 五y ) s ( r ，+理) 。7 其中，圪为施加的边界正弦激励电压；8 ( x ，y ) 为管内介质的相对介电常数分布； r lr ：，r 。：为1 2 个表面电极的位置；f ，为屏蔽层的位置；f 昭为保护电极的位置。 根据( 2 1 1 ) ，( 2 1 2 ) 计算出电位分布矽( x ，y ) 后，利用下式即可求解i j 极板 对间的电容值 q 广一古叮m ，j ，) v 矽( x , y ) d r ， ( 2 1 3 ) 其中，y 为两电极板之间的电位差。 2 3e c t 正问题 2 3 1 正问题求解方法 e c t 重建包括两个计算过程：即正演过程和反演过程。其中正演过程由已知 场域内介电常数分布求解电极对间的电容值。e c t 正演过程的求解方法主要包括 数值计算法和解析法【2 8 3 2 1 ；其中，解析法主要包括复位函数法、分离变量法等； 数值计算法 1 7 , 3 3 - 3 6 】包括边界元法( b e m ) 、有限差分法( f d m ) 、有限元法( f e m ) 等。 。， 解析法过程较为复杂，仅适用于二维均匀场的规则分布情况下的计算。由于 被测场域中的介质分布比较复杂，且多为三维不均匀区域，所以解析法不太适合 e c t 正问题的研究。 有限差分法较为简单，只需针对每一节点写微分方程，并且用差分方程代替 导数，求解该线性方程组便可得到问题的近似数值解。有限差分法对于简单工程 问题的求解效果比较好，但是应用于具有复杂几何分布和边界条件的实际工程问 题时，计算过程相对复杂。 第2 章e c t 系统及其工作原理 边界元法【蚓只在定义域的边界上划分单元，用满足控制方程的插值函数去逼 近边界条件。但是边界元法很难适用于非均匀介质等问题，故其适用范围远不如 有限元法广泛。 有限元法【3 7 】是目前求解各类实际工程问题应用最为广泛的一种数值模拟计算 方法。它基于变分原理，把连续的几何结构离散成有限个单元，并且在每一个单 元中设定有限元节点，从而将连续体看作仅在节点处相连接的一组单元的集合体， 同时选定场函数的节点值作为基本未知量，并在每一单元中假设一个近似插值函 数以表示单元中场函数的分布规律，再建立用于求解节点未知量的有限元方程组， 从而将一个连续域中的无限自由度问题转化为离散域中的有限自由度问题，求解 得到节点值后就可以通过设定的插值函数确定单元上以至整个集合体上的场函 数。 有限元法的最大优点是对不规则区域的适应性好，而e c t 问题场域几何特征 不规则，目前基于有限元法的商业软件如c o m s o l 、a n s y s 等已经被广泛应用【3 8 】。 本文在求解e c t 正问题的过程中，通过自己编写有限元剖分程序，对两相流典型 流型的敏感场场域进行建模仿真，得到的电容值作为图像重构的样本数据。但有 限元法计算的工作量一般比较大，因此改进有限元求解方法，为e c t 的反演结果 的准确性提供保障，也是以后e c t 问题的研究方向之一。 2 3 2 有限元法 本文e c t 正问题所采用的有限元法的步骤为： 1 ) 基于变分原理，将求解域离散化为有限单元，即将连续体问题分解为节点 和单元等个体问题： 2 ) 假设代表单元物理行为的行函数，即假设代表单元解得近似连续函数，建 立单元方程； 3 ) 归一化各单位元方程，组合成整体方程，施加边界条件、初始条件和载荷； 4 ) 求解线性或非线性的微分方程组，得到节点求解结果。 基于贝叶斯理论的电容层析成像算法研究 ( 一) 变分原理 变分问题【5 4 1 是指研究泛函极值问题。如前所述，在对e c t 正问题的研究中， 将被测场域简化为静电场，利用静电场的汤姆逊定理，式( 2 1 2 ) 可转化为 删) = 三小( x , y ) e 2 d s = 哗( ( 掣) 2 + ( 掣) 2 ) d x d y ：m i n ， 。dd 。“ 。j 1 = ： x x , ，y y ) ，c r f 。7 。后：。，2 ，。2 ，七。且。x ，y ，。r ，+ r 馏，c 2 4 ， 式中，d 为被测管内所包含的区域。 ( 二) 三角形剖分 二维有限元有多种剖分和插值的方法，剖分的单元结构包括三角形、四边形 等，由于本文研究的主要是二维简单区域的e c t 问题，因此选择三角形剖分和相 应的三个顶点线性插值方法。 在进行剖分时，应遵循以下规贝1 j 3 7 , 1 3 】 ( 1 ) 应避免出现太大或太小的三角形； ( 2 ) 在电位变化比较大的区域，剖分更为细密； ( 3 ) 任意剖分三角形的顶点，必须是相邻其它单位元的顶点，而不能是其它 元的内点； ( 4 ) 为方便数据统计，对剖分单元按从内到外、逆时针方向依次连续编号， 对单位元的节点编号也是逆时针的顺序标记。 ( 三) 线性插值函数 取一个三角元，对其内部进行线性插值，即用插值函数 矽( x ，y ) = 口+ f i x + y y ， ( 2 1 5 ) 近似代替该三角元内的待求电位分布函数( z ，y ) ，其中三个顶点的电位值分别为 办，唬，九，节点坐标分别记为( 五，m ) ，( 乇，此) ，( 毛，乃) ，节点编号为( 1 ，2 ，3 ) 。 具体如下图所示 第2 章e c t 系统及其工作原理 由插值函数可得， 有上式解得口，7 为， 图2 3 一个剖分单元 f i g 2 3as p l i tu n i t 么= ( 五，乃) = 口+ + 7 y i ， 欢= 妒( 恐，y 2 ) = 口+ 恐+ 7 儿， 死= 痧( 玛，y 3 ) = 口+ p x 3 + y y 3 ， ( 2 1 6 ) 口：竺! 垒丝丝堡垒 ：丝掣， ( 2 1 7 ) 2 、 y ：刍亟垒丝鱼垒 7 2 式中，= 吾( 6 l c 2 6 2 q ) 为取得的三角元的面积 a l = 恐乃一毛咒，6 l = y 2 一y 3 ，c l = 恐一x 2 ， a 2 = 为m x t y 3 ，6 2 = 乃一y l ，c 2 = 五一毛， ( 2 18 ) a 32 五虼一x 2 y l ，岛2 y l y 2 ，c 3 = x 2 一玉 因此，该单位元上的电位分布函数为 = 矛( x ，y ) = i t ( 口+ 岛z + q j ，) 办+ ( 口：+ 包工+ c 2 j ，) 唬+ ( 口3 + 岛x + c j y ) 红】( 2 1 9 ) 在此引入一个基函数畔( z ，少) = 2 - - 去( a , + 1 5 i x + q y ) ，f = 1 ，2 ，3 则( 2 1 8 ) 式改写成如 下形式 基于贝叶斯理论的电容层析成像算法研究 = 矛( z ，y ) = m ，m ，m 门【磊唬九】r = m l l ， ( 2 2 0 ) 其中，【办欢织】7 = 【】。 由式( 2 1 7 ) 可得， v 孑= c 罢考，r = 去匿乏乏卜， c 2 ， 优 卵z ic lc c 1l 即】= 越乏斗则( 2 加) 可写成如下矩阵形式 v = b l 【矽】 ( 2 2 2 ) 由汤姆逊定理可得，单位元的能量积分为 w ) 以( 醅三肛( 础妒出= 庐笋( ( 掣) 2 + ( 掣) 2 ) 蚴，( 2 2 3 ) 。d 。id 。 。 。j 以( 矽) 厶( 多) = 毒舻( x ，y ) 孑 r v 孑】出咖 ( 2 2 4 ) 记 】。为单位元上的电场能系数矩阵，且设 【】= 肛【曰】7 b ld x d y ， | 6 f + c f 6 l 包+ c 。c 2 6 1 6 3 + c i c 3 i = 去协+ c 2 q 包吃+ c 2 c 26 2 6 3 + c 2 c 3l， ( 2 2 5 ) l6 3 6 i + 乞c lb a + c 3 c 2 岛6 3 + c 3 c 3l 因此，( 2 2 3 ) 转化为如下矩阵形式， 厶( ) 厶( 矛) = 】r 旭 ( 2 2 6 ) ( 四) 总体方程组集成 为了整个被测场域内所有节点电位的离散方程组，将( 2 2 6 ) 式中单位元的插 值函数进行扩展，可写为 ，( 矽) 厶( 孑) = 皋卯【】【矽】- m i n ， ( 2 2 7 ) 第2 章e c t 系统及其工作原理 由函数极值理论可得，警：0 ，即 d 口 】 】_ 0 ， ( 2 2 8 ) 下面对其进行强加边界条件处理，则( 2 2 8 ) 变为 】【矽】= b 】，( 2 2 9 ) 其中， 】为，z 刀阶的对称正定矩阵，我们定义为刚度矩阵；依次对每个单元算 出该单元对应的系数矩阵 l ，最后再累加到刚度矩w n 】中。 曰】为边界条件离散后的 1 ) 的向量；眵 为各节点的待求点位值；，z 为三角 形剖分的节点个数。 因为刚度矩阵为满秩、对称的良态矩阵，故e c t 正问题的，z 个结点电位为 矽】= 】叫 b 】( 2 3 0 ) ( 五) 离散电容值计算 在获得传感器内电位分布数值后，现有的方法主要是通过高斯通量积分法或 电荷法计算待测极板上的电荷量q ，进而用静电场的相关理论求得激励电极和待 测电极间的电容值c 。例如，电极1 与电极i 之间的电容q ，为 c i - - 害， ( 2 - 3 1 ) 其中，i = 2 ，3 ，m ，m 为系统电极数。 上节推论可知，z 个结点电位分布求得后，可根据高斯通量定理求出系统f j 极板对间的电容值为式( 2 1 3 ) 。 由于高斯通量法需要计算待测极板的曲面积分，因此造成计算量很大，且随 着系统中电极数目的增加，计算所需时间会呈现指数倍增长，甚至大于有限元计 算电位分布所用时间。因此从时效性考虑，本文采用相对简单的电荷法来计算待 测极板上的电容值c 。 先计算待测极板上的电荷量q g = ( 删) 。， ( 2 3 2 ) 其中，( 删) 。是刚度矩阵 】的第n 行，即对应于第，个待测电极上的节点n 与场 域内剖分节点的离散电位值 矽 的标量积，等于节点n 上的电荷量，再由式( 2 3 1 ) 求得激励电极和待测电极之间的电容值。 ( 六) e c t 正问题仿真 如下图2 4 所示的有限元剖分结构，分别为8 电极、1 2 电极、1 6 电极。 - i 舳舶“00 2 0 4 0 e 0 t l 图2 4 ( a ) 电极数= 8 ，管内剖分层数= 8 ；剖分单元数= 2 5 6 f i g 2 4 ( a ) e l e c t r o d e sn u m b e r = 8 ，t h et u b es p l i tl a y e r s = 8 ；s p l i tu n i tn u m b e r = 2 5 6 4 棚0 1 “m00 20 40 6n l1 图2 4 ( b ) 电极数= 1 2 ，管内剖分层数= 1 2 ，剖分单元数= 5 7 6 f i g 2 4 ( b ) e l e c t r o d e sn u m b e r = 12 ，t h et u b es p l i tl a y e r s = 12 ；s p l i tu n i tn u m b e r = 5 7 6 第2 章e c t 系统及其工作原理 - i 0 j - 0 6 也4 0 200 20 40 , 6 图2 4 ( c ) 电极数= 1 6 ，管内剖分层数= 1 6 ，剖分单元数= 1 0 2 4 f i g 2 4 ( c ) e l e c t r o d e sn u m b e r = 16 ，t h em b es p l i tl a y e r s = 16 ；s p l i tu n i tn u m b e r = 10 2 4 图2 - 4 所示的是采用三角形剖分的有限元分割模型，对于8 电极的仿真模型 1 6 l ，管内剖分单元总数为2 5 6 个，剖分节点总数为3 6 1 ，管内区域介质相对介电常 数设为1 ，管壁相对介电常数设为3 ；则矩阵 为( 3 6 1 x 3 6 1 ) 的方阵， b 是3 6 1 1 的列向量。由于管道内无自由电荷，即边界条件 曰】是一零向量，只需求出系数矩 阵 】的值。按照剖分单元的编号，依次对单元f ( 1 ，2 ，2 5 6 ) 算出该单元对应的系 数n ，( x = l ，2 ，3 ；y = 1 ，2 ，3 ) ( 其中每个单元的三个顶点编号逆时针为1 ，2 ，3 ) ，并 将这些系数累加到刚度矩阵 中；最后所有单元都计算完毕时，刚度矩阵 】就 确定了。下表给出了8 电极实际测量值和有限元计算值对比， 基丁= 贝叶斯理论的电容层析成像算法研究 表2 5 有限元计算的理论数据和实际测量数据对照表 t a b 2 5t h et h e o r yo ff i n i t ed e m e n td a t aa n dt h ea c t u a lm e a s u r e m e n td a t a t a b l e 全油全水 序号电极对 测量值理论值测量值理论值 11 20 4 6 2 20 4 6 7 30 9 8 7 20 9 8 6 2 21 30 0 5 6 4 0 0 5 6 2o 1 7 3 lo 1 7 5 8 31 - 40 0 3 2 20 0 3 5 30 0 9 6 4 0 0 8 9 8 4 1 50 0 2 6 30 0 2 7 7o 0 8 1 50 0 8 5 3 51 6 0 0 3 2 50 0 3 3 70 0 9 6 10 0 9 4 7 61 70 0 5 8 20 0 5 8 40