（计算数学专业论文）介质导电率成像的反演问题及算法实现.pdf

摘要 电阻抗成像( e i t ，e l e c t r i c a li m p e d a n c et o m o g r a p h y ) 是当今生物医学工程领域的 重大研究课题之一，数学上对应于一类经典的椭圆型方程的系数反问题传统的 e i t 技术由于仅利用了介质外部的直接测量数据，得到的图像的分辨率较低，很 难满足应用问题的需要近年来发展起来的核磁共振成像技术( m r e i t ，m a g n e t i c r e s o n a n c ee l e c t r i c a li m p e d a n c et o m o g r a p h y ) 是改进成像精度的一个新的研究方向，它 利用间接测量到的介质内部的磁场信息来反演介质内部的导电率，由此有望得到 稳定的高精度的成像结果数学上而言，它是一类由椭圆型方程边值问题解的内部 信息来反演方程的待定系数的新的反问题由于利用了解在内部的有关数据，而 不仅是边界数据，因此问题的不适定性有所减弱 本文讨论该反问题在一个具有明确医学背景的模型上的算法一调和b ；算法的 实现问题该模型可以用下面的带积分边值条件的椭圆型方程定解问题来描述： 上一瞎d s r n r u 一 r 0 n s + u 一 在分析了算法基本思想和引起算法不适定性的因素的基础上，本文重点考虑了利 用正则化方法处理该算法中扰动输入数据的求导不适定性的问题，在此基础上提 出了b 反演方法对扰动数据的数值实现通过计算机模拟实验验证，取得了令人 满意的结果 本文的工作包含了以下四个方面的内容首先，与目前已有的模型问题相比， 我们讨论了一个更为现实的模型，利用的是边界电极上的总的输入电流的强度，数 学上这是一个带有积分边界条件的定解问题其次，我们把此问题转化为一个标准 的椭圆型方程混合边值问题，自编了有限元程序得到了有效的正问题的数值解， 实现了反演数据的模拟，同时也为尻反演方法中方程的迭代求解提供了一个必要 的工具再次，我们考虑了利用扰动输入数据b 。时，算法的稳定实现问题调和 最反演方法中的关键一步是要对输入数据进行l a p l a c e 运算众所周知要对扰动数 据求数值微分是个不适定的问题，对此我们应用基于t i k h o n o v e 正则化方法的五次 样条函数的数值微分方法，提出了一个稳定的数值求导方法，大大减弱了迭代算 法中的不适定性，并在适当的假设下得到了误差估计最后，我们给出了所提出反 演方法的数值实现，尤其是对扰动数据的实现效果，验证了本文的理论结果 关键词；反问题，椭圆型方程，有限元法，调和b ：算法，样条函数，正则化 数值微分，迭代算法，数值解 i 鬻掌 呻饕畦 a b s t r a c t e l e c t r i c a li m p e d a n c et o m o g r a p h y ( e i t ) i so n eo ft h em o s ti m p o r t a n tr e s e a r c ha r e a si nt h e b i o m e d i c a le n g i n e e r i n g w h e r e a se i t t e c h n i q u ch a s t h ed i s a d v a n t a g e so fs e v e r ei l l - p o s e d n c s sa d l o ws p a t i a lr e s o l u t i o n ，m u c ha t t e n t i o na l ep a i dt ot h em r e i t ( m a g n e t i cr e s o n a n c ee l e c t r i c a l i m p e d a n c et o m o g r a p h y ) t e c h n i q u e ，w h i c hi sar e c e n t l y - d e v e l o p e dd y n a m i cc o n d u c t i v i t yi m a g i n g m o d a l i t y t h eb a s i ci d e ao fm r e i ti st or e c o n s t r u c tt h em e d i u mc o n d u c t i v i t yu s i n gt h ei n t e r n a l m a g n e t i cf l u xi n f o r m a t i o n d u et ot h ea p p l i c a t i o no fi n t e r n a le l e c t r i c a lc u r r e n td i s t r i b u t i o n ，h e t h t h ca c c u r a c ya n dt h es p a t i a lr e s o l u t i o no ft h ei m a g e sa r ee x p e c t e d t h i sp a p e rc o n s i d e r st h er e c o n s t r u c t i o ns c h e m eb a s e do nt h ep a r t i a l l ym e a s u r e dm a g n e t i c f l u xd a t a 0 u rm o d e lp r o b l e mc a nb ed e s c r i b e db yt h ef o l l o w i n gs y s t e m ： iv - ( a v u ) = 0 r n jj = j ：+ a 舞d s = 一止一a 舞d s lv uxn = 0r + u 5 一 i 一口器= 0 r o n i 干i i i 3 w ea d o p tt i l er e c e n t l y - d e v e l o p e dh a r m o n i c b za l g o r i t h mt or e c o n s t r u c tt h ec o n d u c t i v i t y 口 c o n s i d e r i n gt h cu n a v o i d a b l en o i s ei nt h ei n p u tb ；d a t a ，o u rm a i nt a s ki st op r o p o s ea ne f f i c i e n t c o m p u t a t i o n a ls c h e m ef o rv 2 l ，s u c ht h a tt h eh a r m o n i c - b za l g o r i t h mc a ng e n e r a t eac o n v e r g e n t a n ds t a b l es o l u t i o n o u rr e s e a r c hw o r ki nt h i sp a p e rc o n f a i n sf o u rp a r t s f i r s to fa l l ，w ec o n s i d e ram o r e p r a c t i c a lm o d e lu s i n gt h et o t a la m o u n to fi n p u tc u r r e n tf r o mt h ec l c c t r o d e s b yi n t r o d u c i n ga f u n c t i o nt r a n s f o r m ，w ec o n v e r tt h i se l l i p t i cp r o b l e mw i t hn o n l o c a lb o u n d a r yc o n d i t i o ni n t oo n e w i t hs t a n d a r dm i x e db o u n d a r yc o n d i t i o n s e c o n d l y , w es o l v et h i sn o n s t a n d a r dp r o b l e mb yf e m m e t h o d t h ed i r e c ts o l v e rd e v e l o p e di nt h i sp a p e ri sr e q u i r e df o rt h es i m u l a t i o no fi n v e r s i o n i n p u td a t a i ta l s op r o v i d e sa ne f f i c i e n tt o o lf o rt h ei t e r a t i v ep r o c e d u r ci nt h eh a r m o n i c - b = r e c o n s t r u c t i o na l g o r i t h m s t h et h i r dw o r ki nt h i st h e s i si st od i s c u s st h en u m e r i c a ld i f f e r e n t i a t i o n p r o b l e mf o rn o i s yd a t a w ec o n s t r u c ta na p p r o x i m a t ef u n c t i o nf r o mt h et i k h o n o vr e g u l a r i z i n g f u n c t i o n a lw i t he x p l i c i te x p r e s s i o nb a s e do nt h es p l i n ef u n c t i o n s t h i st r e a t m e n tw e a k e n st h ei l l - p o s e d n e s so ft h eh a r m o n i c - b za l g o r i t h ma n dp r o v i d e sac o n v e r g e n ti t e r a t i o ns e q u e n c e f i n a u y , n u m e r i c a le x a m p l e sa r eg i v e nt os h o wt h ev a l i d i t yo ft h ep r o p o s e dm e t h o d k e y w o r d s ：i n v e r s ep r o b l e m ，e l l i p t i ce q u a t i o n ，f e m ，h a r m o n i e - b za l g o r i t h m ，s p l i n ef u n c t i o n s ，r e g u l a r i z a t i o n ，n u m e r i c a ld i f f e r e n t i a t i o n ，i t e r a t i o n ，n u m e r i c s 东南大学学位论文 独创性声明及使用授权的说明 一学位论文独创性声明 本人声明所呈交的学位论文是我个人在导师指导下进行的研究工作及取得的 研究成果尽我所知，除了文中特别加以标明和致谢的地方外，论文中不包含其他 人已经发表或撰写过的研究成果，也不包含为获得东南大学或其它教育机构的学 位或证书而使用过的材料与我一同工作的同志对本研究所做的任何贡献均已在 论文中作了明确的说明并表示了谢意 ：，关于学位论文使用授权的说明 东南大学、中国科学技术信息研究所、国家图书馆有权保留本人所送交学位 论文的复印件和电子文档，可以采用影印、缩印或其他复制手段保存论文本人电 子文档的内容和纸质论文的内容相一致除在保密期内的保密论文外，允许论文 被查阅和借阅，可以公布( 包括刊登) 论文的全部或部分内容论文的公布( 包括刊 登) 授权东南大学研究生院办理 期， 第一章引言 1 1 介质成像的有关背景 数学物理反问题研究的内容非常广泛，涉及到医学成像，热流逆传导，地球物 理勘探。遥感技术，信号( 图象) 处理，工业控制乃至经济决策等众多的科学技术领 域当我们把具体的物理现象用微分方程定解问题来描述时，反问题一般是重建 微分方程中未知的系数区域形状或一些边界及初始条件，求解反问题是较困难 的，因为很多反问题都是不适定的并且通常是非线性的这意味着当我们利用的原 始数据是有扰动的测量数据时，在经典意义下问题的解可熊不存在解即使存在 也可能不唯一，且可能不连续依赖于输入数据，即原始数据小的观测误差会导致 近似解与真解的严重偏离处理不适定问题的基本方法是正则化理论对一般的 不适定问题，典型的求解方法有t i k h o n o v 正则化，l a n d w e b e r 迭代方法，及基于随 机分布理论的随机正则化方法等 因为很多不适定问题的困难本质上是由给定的输入数据的不足或不精确而产 生的，于是人们试图通过提高测量仪器的精度或采用新的技术来获得更多的信息 从而减少问题的不适定性这是处理一般不适定舟题的从输入信息角度可以采用 的一个改进途径本论文考虑具有医学成像背景的由核磁共振技术( m r e i t ) 重建 生物体内部的介质导电率的问题。数学上该问题由一类椭圆型方程的系数反问题 来描述通过间接测量介质内部磁场分布的相关信息，m r e i t 技术基于m a x w e l l 方 程组试图利用不完全的磁场数据重建介质导电率该技术采用的数学模型与经典 的e i t 三维介质成像的数学描述有所不同通常的基于微分方程系数反演的e i t 介 质成像方法，所利用的输入数据都是解在介质外部的观测信息，解对输入数据的 敏感性太强，因而是一类强不适定问题上世纪九十年代新发展的核磁共振技术 ( m r e i t ) ，利用技术上可以实际测量得到介质内部的一个方向的磁场分布( 而不仅 是边界测量数据) 作为成像的输入数据，大大减弱了问题的不适定性，更可能得到 稳定的数值结果 本章首先对传统阻抗成像技术及其研究意义进行了简单的介绍，并对其技术 的难点进行了分析，在此基础上引出了基于磁共振成像系统的阻抗成像技术，分 析了该技术的优点和算法上还存在的困难最后介绍了本硕士学位论文为克服此 困难而进行的数值计算方面的主要研究工作，它主要包含了一类带有非局部边界 条件的椭圆型方程的数值求解和噪音数据的数值微分的问题在适当的假设条件 下，我们还建立了利用噪音输入数据时调和乜算法的收敛速度的估计 1 1 1 电阻抗成像( e i t ) 技术 电阻抗成像技术( e l e c t r i c a li m p e d a n c et o m o g r a p h y , e i t ) 是当今生物医学工程重 l 要研究课题之一e i t 是在医学检测c t ( c o m p u t e rt o m o g r a p h y ) 成像技术的基础上 进一步发展起来的( 【l 】) 它根据人体内部不同组织具有不同的电阻抗这一物理原 理，借助于附加电极对生物体外加一定的电流，通过测得生物体外界的电量信息 来重构生物体的阻抗分布 电阻抗成像技术是通过物体表面电极的测量值来对物体内部导电率的分布进 行重建，数学上对应于一类由d i r i c h l e t - t o - n e u m a n n 映射作为反演数据的椭圆型方程 的系数反问题这种成像方式本身决定了e i t 技术目前的发展遇到了一些根本的困 难这主要是由于以下的两个原因( 1 2 1 ) 首先是问题的不适定性通过物体表面电极 的测量值对物体内部导电率的分布进行重建，由于边界测量值对物体内部导电率 的变化十分不敏感，边界测量值的微小变化会导致求解出的内部导电率有很大变 化，使得求解过程很不稳定其次是得到的导电率分布的图像分辨率不高通过物 体表面电极的测量值对物体内部导电率的分布进行重建，获得的信息量是有限的， 因为测量上不可能给出无穷组的边界电流电位数据，因此由d i r i c h l e t - t o - n e u m a n n 映 射的有限的近似数据重建得到的物体导电率图像内部易出现模糊 1 1 2 基于磁共振原理的电阻抗成像( m r e i t ) 如上所述，目前e i t 技术面临求解过程的不适定性和得到的物体导电率分布 图像分辨率不高这两个难点这两个问题是e i t 本身所固有的，因此，有必要对传 统生物电阻抗成像技术的数据获取方式与图像重建技术进行改进，以得到高分辨 率，高精确度、适用于临床医学诊断的生物电阻抗分布图像基于磁共振成像原理 的生物阻抗成像技术( m r m i t ) 就是一个主要的研究方向，数学上它是一类椭圆型 方程利用解的内部信息重建系数的反问题 m r e i t 是传统的电阻抗成像技术( e i t ) 与电流密度成像( c u r r e n t d e n s i t yi m a g i n g ， c d i ) 的结合，最早由加拿大t o r o n t o 大学的j o ym ，p e i k a np 和s c o t tg 等人于1 9 8 9 年 提出其基本原理( 【3 】，【4 】) 该技术通过物体表面的电极向物体注入电流，利用磁共振 成像系统( m r i ) 测量介质内部电流在物体内部产生的磁感应强度分布由著名的 m a x w e n 方程组描述的电磁场的关系，这实际上是给出了介质内部电流密度分布的 一种间接测量方法而m r e i t 技术上的概念则是在1 9 9 2 年的一篇学位论文中提出 的( f 5 】) 由于c d i 技术利用m r i 系统得到物体内部电流密度分布，因此物体内部精 确的导电率可以通过m r 图像间接得到同时由于利用c d i 得到的物体内部电流 密度分布其分辨率在整个物体范围内是一致的，因此解决了传统e i t 边界测量值 对物体内部导电率的变化不敏感的问题，因而更可能得到物体高分辨率、高精确 度的导电率分布图目前该技术受到国际生物医学界越来越多的关注 由于该技术出现得比较晚，目前国际上只有英国，韩国，德国，土耳其、加 拿大等国的几个研究小组在进行该技术的研究，且均处于理论研究阶段。尚未将 m r e i t 技术应用于人体具体部位的阻抗成像。与真正的医学临床应用尚有一段距 离该问题数学上的困难在于。利用m r i 测量得到的的数据在某些部分真有较强 的噪音，基于椭圆型方程系数反演的成像算法得到的数值解难以达到可以应用的 分辨率国内大多是对传统e i t 的研究，m r e i t 的研究还很少 1 2m r e i t 巳有的相关工作 m r e i t 技术最早的提出者之一是nz h a n g ( s ) ，其基本思想是利用物体内电流 密度分布与边界电压测量值来确定物体内导电率的分布由于边界任何两点之间 的电压差是连接这两点的任何一条线的电场强度的积分，而电场强度是电阻率与 电流密度分布的乘积，即e = p j ，其中p 是电阻率假定物体内的电流密度分布 可以间接测量得到，则对连接两点的不同的线和不同的边界点对，可以建立一系 列线性方程，通过求解该方程组则可得到物体内的电阻率分布这个方法的缺点 是需要许多边界电压测量值来提高图像的精确度与空间分辨率m r e i t 经过了十 多年的发展，已经形成了不少标准的结果和成像算法，下面我们简单总结一下这 些方法 1 基于等势线的m r e i t 算法 o z d e m i r 和e y i i b o g l u ( 2 0 0 2 ) ( 6 ) ，及k w o ne ta l ( 2 0 0 2 ) ( 7 d 】提出了基于等势线的 m r e i t 算法物体内电流密度分布为己知，而且电流密度分布与等势线是垂直的， 寻找等势线并应用边界电压测量值。则物体内电场分布可以求得，利用求得的电 场分布、测量得到的电流密度分布及欧姆定律则可以求得物体内电阻率分布这 是一种非迭代方法。但需要许多边界电压测量值与单个输入电流 2 电流迭代算法 ejw o oe ta 1 ( 1 9 9 4 ) ( 8 ) 与ok w o ne ta l ( 2 0 0 2 ) ( 【9 】) 等提出了电流迭代( j - s u b s t i t u t i o n ) 算法这种方法利用有限元方法，使得对每个有限单元来讲，计算得到的电流密度 分布的幅值与测量得到的电流密度分布的幅值最小，是一种迭代算法但为了得 到电流密度分布，必须利用三个方向的磁感应强度，因此成像物体需要在m 砒系 统中做三维转动，这在技术上难度较大 3 调和磁场( h a r m o n i c - b z ) 算法 jks e o 等人于2 0 0 3 年提出h a r m o n i c - 色算法( 1 0 】) 该算法利用两个注入电流测 量得到的单个方向的磁感应强度对物体内的导电率分布进行估计该算法的原理 是将内部电位叫口l 和内部磁通量向量b ( r ) 的关系 v 2 口= 瑚vx ( a v u ) = 一i i o v u v 叮 与描述电位分布的椭圆型方程相结合，对物体导电率在每一个垂直于。轴的介质 截面上逐层成像，从而本质上是一个二维椭圆型方程的系数反问题该反演方法的 基本思想是由l a p l a c e 方程的基本癣把介质内部的待求的导电率表示为毋的二阶 导数和相应定解问题的解的一个积分，于是可以得到该算法的迭代公式 jjl i u 等人在这方面做了稳定性和收敛性的分析( 【1 l 】， 1 2 1 ) ，为之提供基本的数学理论 这种算法需要对磁感应强度求两次微分，因此对噪音比较敏感对忍数据有噪音 的情况，如何去噪是关键，该问题也是本文的主要研究内容 4 梯度玩分解算法 cp a r k 等人于2 0 0 4 年提出了梯度玩分解算法( g r a d i e n t 最d e c o m p o s i t i o n ) ( 1 3 ) 该算法首先将所求区域分成许多薄层，对每个层分别求导电率分布对每个层m 而言。方向导电率是不变的，并且设定一个比薄层范围更宽的层n 6 ，并假设薄层 n ，内的电势分布只受薄层n 。与n d 内的导电率分布的影响，n j 之外的导电率对其 电势分布没有影响在上述假设条件下。对q 内的单方向磁感应强度岛进行分 解，使得尾成为n ，内电势的函数，然后结合正问题中的椭圆型方程。对导电率分 布进行估计 综合上面对物体内部导电率已有的重构算法的介绍，根据获取信息类型的不 同，对m r e i t 技术已有的算法大体上分为两类 第一类是利用物体内部电流密度分布对物体内导电率进行重建nz h a n g 在 其硕士论文中提出的m r e i t 算法，基于等势线的m r e i t 算法、电流迭代算法，均 属于此类算法这种类型的m r e i t 算法需要的输入数据信息量相对较多，理论上 该类反问题具有很高的反演精度和稳定性但是这种算法的一个最根本的弱点就 是为了得到介质内部完整的电流向量的分布，根据标准的m a x w e l l 方程组的关系， 必须测最介质内部三个方向的磁感应强度，因此成像物体需要在m r i 系统中做三 维转动这在实际临床应用中是难以实现的 第二类是根据物体内部磁感应强度的的部分分布数据对物体内导电率的分布 进行重建j ks e o 等人于2 0 0 3 年提出调和见算法，c p a r k 等人于2 0 0 4 年提出的 梯度玩分解算法均属于此类算法上述几种m r e i t 算法由于仅利用了单个方向 的磁感应强度测量值，因而都有效地解决了成像介质在m r i 系统中旋转的问题， 较第一类算法有很大的进步与更大的应用空间，但此类算法也有缺点，就是要对 实验数据毋求数值微分，其中调和展算法需要两次微分，梯度岛分解算法需要 一次微分众所周知数值微分是典型的不适定问题，当忍带有噪音时，需要发展 有效的数值反演算法 目前m r e i t 技术的发展大多处于理论研究阶段。距真正的临床应用尚有距 离m r e i t 技术，尤其是有效的反演算法，还需要进一步深入的研究 1 3 本文的工作 目前已有的关于利用一个磁场方向的反演算法的工作大都是从工程的观点以 圣曼查篓塑圭耋堡墼圣 董三塞 丝堡 5 及实验上出发的，考虑的是具有如下n e u m a n n 边值条件的椭圆型方程模型( 【m 1 ) v一(evu)=0ttvug ：曼， s 【一 n =r 踟 7 其中u ( r ) 为物体内电势分布$ j o 为真空磁导率矿( r ) 为物体内导电率分布。n 为 所要成像的三维物体而物体内部磁场。内部电流和电位满足 肿) ；等) x 高，e l - v u ，= 以 该种边界条件相当于电极上的输入电流的密度分布要逐点已知。这是难以做 到的，实际上只能给出通过每个电极的总电流 本文采用h a r m o n i c - b z 算法，从实际出发考虑了一种更为接近现实的模型。该 模型是由下面的带积分边值条件的椭圆型方程来描述 iv p 钆) = 0 j ，= 正+ a 舞幽；一正一瞎r l d s ， lv “n = 0 【一砖= o r n r e - i - u r ( 1 删 r a q e + u 暑一 在给定一个方向的磁通量数据鼠的条件下来重建导电率a ( r ) 的问题，其中的一个 主要工作就是要处理调和兄算法中对噪音输入数据的二阶微分问题 本文的主要工作如下 在第二章我们考虑了m r e i t 技术中电位分布的正问题的数值求解，我们通过 引理2 1 先将( 1 3 2 ) 转化为一个标准的椭圆型方程混合边值问题，再用有限元法求 解正问题由给定的已知的精确导电率矿代入( 1 3 2 ) 求解得到介质内部的电势分 布，模拟产生边界电势差以及介质内部磁通量分布b 该工作为我们模拟产生反 演输入数据和在调和b 算法中迭代求解导电率时提供了必要的工具 第三章中讨论了m r e i t 技术中的调和e 反演算法介绍了调和e 算法及其 实现的具体流程，其中迭代算法是求解关于导电率的非线性方程的的有效工具， 并给出了算法具体的流程图，该项工作清楚地阐明了调和且反演算法的优点及影 响算法精度的主要原因 第四章讨论调和玩反演算法中的不适定性的处理方法，即对扰动输入数据毋 的二次微分因为每一迭代过程都要利用v 2 b 作为积分核，所以我们认为高精度 的数值微分算法是保证调和兄反演算法迭代序列收敛的关键我们用五次样条函 数和正则化方法相结合，得到了稳定的的数值微分。并在适当的假设下得到了误 差估计该处理过程大大减弱了调和风反演算法中的不适定性数值饲子表明了 五次样条函数方法对离散数据。特别是扰动的离散数据的数值微分的不适定性有 很大改善 在第五章讨论调和且算法的数值实现用f o r t r a n 和m a t l a b 语言混合自编了数 值反演的有限元程序，进行了细致的仿真实验验证数值试验表明，这是一个迭代 速度较快的算法，同时，该成像算法具有较高的成像精确度和一定的抗噪音能力 第二章m r e i t 成像技术中的正问题 给定介质内的导电率的分布，对介质内的电势分布。边界电位和介质内部磁 场分布的求解一般称为m r e i t 技术中的正问题众所周知，正问题对相关的反问 题的求解起着基本的作用，反问题的求解离不开正问题对我们所讨论的由椭圆 型方程描述的介质成像的反问题，它同样如此一方面，在反问题的迭代求解中， 在每一迭代步，我们都要解一个正问题以产生对应的电位分布从而由此更新得到 导电率下一个迭代值另一方面，为了模拟产生反问题需要的磁场输入数据，也需 要解一个正问题如果模拟磁场输入数据已包含了较大的误差，则对反演方法的 评价就会产生误导，或者难以评价因此针对具体模型的一个高效的正问题的求 解方法是反问题算法实现并进行评价的关键对本文讨论的带有积分边界条件的 电位分布的正问题的数值求解，在本章我们先把原问题转化为一个等价的标准的 椭圆型方程的混合边值问题，采用有限元的数值方法来求解得到电势分布。本章 建立的正问题求解方案，为反问题的迭代求解和模拟得到反演输入数据( 内部磁场 分布) 提供了保证 2 1m r e i t 正问题的数学描述 众所周知，电流在场域内的流动受场域内导电率分布的影响。导电率分布和 电流之间的关系可用m a x w e l l 方程组的形式表达针对m r e i t 的特殊情况，可做 如下假设，场域内的电流场满足似稳条件；场域内没有频率和外加电流频率相同 的电流源；物体的导电率a 是各向同性的；导电率与电流密度无关给定以上条件 后。导电率分布就可以表示为一个标量函数一( r ) 在以上假设条件下。m r e i t 技术中电位分布的正问题可以用以下简化了的椭 圆方程的边值问题来描述( 1 4 d v ( 7 乳) _ 0睡! ， ( 2 1 1 ) 1 一# v u n 29 r 砚 ”7 其中u ( r ) 为物体内电势分布，( r ) 为物体内导电率分布。n 为所要成像的三维物 体边界条件相当于要求给定边界的电流密度分布该问题的解在相差一个常数 的意义下是唯一的，这对应于电位本身只有相对意义一旦得到物体内部的电势 分布u ( r ) ，则物体内电场强度e 电流密度，可通过以下两式求得； e = 一v 口j = e 给定物体内导电率的分布与边界条件，求解物体内电位分布与磁感应强度的 分布称之为m r e i t 技术中的正问题上面给出的是m r e i t 技术中的一个简化的模 7 型，但在现实中边界电流密度并不能被逐点给出，只能知道通过每个电极的总电 流即解在边界的n e u m a n n 数据并不是真正给出的因此我们从实际出发，考虑一 个更为实际的m r e i t 成像模型我们设介质为nc 帮，边界加是光滑的，导电率 分布口是各向同性的通过勰上一对电极矿，e 一输入电流总量为j ，则边界输入 电流j 在介质内部产生电流密度分布j = ( 五，山，以) = 一f v u 对这样的模型，基于 m a x w e l l 方程组，电势分布t 与物体内导电率的关系可用以下带积分边值条件的椭 圆方程来描述 f v ( 盯v u ) = o r n v 7 i 曼霎0 拈以鼬+ u e - ， ( 2 - 2 ) i “n =r s + 。 i 一口器= 0 r 踟而 根据m a x w e l l 定律。物体内磁通量b 和电流密度分布，有如下关系 伽，一v b ，( 2 1 3 ) 或者由b i o t - s a v a t 定律 b ( r ) = 等j ( r ，) x 品 ( 2 1 4 ) 来描述，其中p o 为真空磁导率 对一般的区域q ，上述边值问题不存在解析解，只能进行数值求解。这不是一 个标准的带局部边界条件的边值问题，但可以利用下面的函数变换将其转化为一 个标准的椭圆型方程的混合边值问题来求解 引理2 1 设豇为混合边值问豚 f v ( 矿v 豇) = o r n 矗l 一= l ，乩一= 0 ( 2 1 5 ) 【一f 豢= 0 r 鲫丽 的解，则t = 赤矗+ g 为( 2 1 2 ) 在n 内的解，其中e 为常数 证明，设= t 一c 豇，则t t ，满足v ( a v 曲= 0 于是由g r e e n 公式和边界电流守恒 性质 厶，蒙幽；一二，蒙幽 并利用( 2 1 2 ) ，( 2 1 5 ) 中的边界条件得 上a 1 1 2 打；厶一筹”幽= 厶，筹幽c 吨+ 一。+ c 丘一蒙幽m 。一 = 一k c ) ( 卜c z + ，蒙凼) 塞童奎堂堑圭堂堡垒塞董三塞些塑坚童墼薹堡星9 所以，对常数c 2 西刁i 甄，可得到在n 上l v 训- = 0 即 在n 上是常数，得证 由上面的引理，我们只要求解定解问题( 2 l 5 ) 得到豇，再由“2 南豇+ g 即 可得到正问题的解，其中c 可由指定介质内( 或边界上) 一点的电位来确定 2 2 正问题的数值解 数值计算偏微分方程边值问题的方法一般包括边界元法( b o u n d a r ye l e m e n t m e t h o d ，b e m ) ，有限元法( f i n i t ee l e m e n tm e t h o d ，f e m ) 和有限差分法( f i n i t ed i f f e r e n c e m e t h o d ，f d m ) 等几种边界元法是形成较晚的一种数值计算方法，是把边界积分方 程法与有限元法的离散方程式组合起来的产物有限元法基本思想是将偏微分方 程边值问题化为一个等价的变分问题，通过有限元剞分将其离散化到各个单元， 并最终将其转化为关于节点处函数值的一个线性代数方程组来求解与有限差分 法耜比，有限元法能更好地模拟电极尺寸和电极与介质之间的接触阻抗边界元法 虽然兼有边界积分法和有限元法的优点，但计算量很大在电磁场计算( 【1 5 】) 中， 有限元法使用最为广芝，在这里我们也采用有限元法来求解m r e i t 中的正问题 2 2 1 有限元方法 应用有限元法求解m r e i t 成像技术中正同题的步骤遵循一般的有限元方法的 基本步骤首先给出与待求边值问题相应的变分问题；进而剖分场域并选取相应 的插值基函数；在每个单元上计算单元刚度矩阵和单元荷载向量；合成单元刚度 矩阵和单元荷载向量得到总刚度矩阵和总荷载向量；处理边值条件( 主要是第一 边值条件的处理) ；导出一组联立的代数方程( 有限元方程，其系数矩阵是大型稀 疏矩阵) ；选择适当的代数解法求解有限元方程，最后得到待求边值问题的数值解 由于对一般的三维介质q ，m r e i t 成像技术也是化为在每一个垂直于：方向的= 维 截面上来求饵的，因此我们这里仅器要考虑如下二维区域带混合边界条件的求解 1 青况。设n = n n q ，弘。) ：= 为) lv p v = 0r n j = 五+ 口舞幽= 一正一一舞幽，( 2 2 1 ) lv “n = 0 r + u e 一 、7 l 一仃甓= 0 r a n 句+ u 根据引理2 1 ，求解此问题本质上只要求解如下椭圆型方程混合边值问题 窖= 。 r n 如 r 留可f 坠翼垄兰墅坠奏! l 冀：i 坠一： 董三塞些g i 堡塞墼垂塑矍 1 0 对此混合边值问题，取试验函数空间为矿= 扣日1 ( q 。) ：口i d ；o 。并记 g d 1 ( n ) 满足g d l t + = l ，g d l = 0 是声上u 的d i r i c h l e t 数据在n 上。上的一个已知 箜要拓函数则上述向题的解等价于满足下面的变分问题；求函效。满足。一妨y 僖得 d ( u ，一j ( 2 0 ，坳m ( z z s l 其中 、 d 们。j ( k 神 珊蛐，p 厶埘矿撕一) 础 ( 2 2 2 ) 和( 2 2 3 ) 的等价性可以证明如下 如果“满足( 2 2 2 ) ，则显然一卵y 并且对一切的t ，v ，有 忙厶玑c 两咖2 五k 。象俨。，a v u v v = - 上褂妊+ 恤一，鲫一z 。，讥珊， 即( 2 2 3 ) 成立 反之，假定( 2 2 3 ) 对一切的 y 成立则我们先取 y 为任意满足8 卿t ，c c n 目 的函数，则( 2 2 3 ) 和。的紧支集性质表明 上，。玑p 乳”一z k a v u , v v 厶埘一) 鲫2 o 由口的任意性即得v ( a v u ) t0 在0 。上成立，由此对任意的p 瘩y 由( 2 2 3 ) 得 o - 厶仉c 内咖；k ，囊”一上。概渤 2 五褂。+ 。一，囊”十五埘扛+ 。一，= 上址+ 。，p 蒙+ 印” 下面要来说明对任意的口n 训出k 、( 。+ ) 也可以是任意的该结果不是显然的。 注意这里鼬m p + u r ) 只是a n 。的一个闭集其证明需要在并的一定条件下的 坐枥 变换，具体细节见( f 1 6 j ，瓶1 9 ) $ j p o i s s o n 方程混合边值条件的证明( p 1 3 1 ) ) 从 而由口j 房“”一) 的任意性得 f 象( r ) + 口( r ) 一o ，r a f k 忙+ u e 一) 从而t 满足n e u m a u n 边界条件另一方面，由于“一窖de n 故t 满足s 上的d i 鼬l e 边界条件 由此等价性，就可以对变分形式用有限元方法求解来得到问题( 2 2 2 ) 的解式 ( 2 2 。3 ) 的有限元近似为，对给定的有限维子空间h c y ，求满足一l h 卵h 使 得 d ，v h ) 一f ( 鲰) 一0 ，怕h e ， 其中霸是把9 d ( z ) 投影到玩上的全局插值算子对f 的有限单元离散。及对应 的的基函数的构造，可由标准的方法如下进行令单元e ( 图2 1 ) 内电势分布函 数近似为； “0 ，f ) u h ( z ，) = + 6 l ，+ c 图2 1 ：三角形单元 ( 2 2 4 ) 设t “，坳，分别为顶点最，弓，p m 的电位，并记t 。= h ，叼，l r 利用标准的线 性有限元的方法，可以求出在单元e 上有 撕铋功；，= 蕃 = 玩， 其中一m 扣，妁，m 仁，) n m ( ：r ，口) 】， m = 去 三曼；j 吩= 击j 三墨汁m = 击l 三兰；| c z 卫s ， b = 蓦喜蕃】= 去 ：z 乏 ，色；| 主曼 i ，q ，b 轨 新h 茁t 甜勺= j ： l ，：l 矾 ii 珊 i , b m = f liq 跏ii 雏l 向。i 斗 ：l ， ( 2 。6 ) 以弘l 句珊l 。 我们考虑求解区域o 。为矩形的情况先把矩形等分成c = m 个小的矩 形，再把每个小矩形分成两个直角三角形元，记剖分后得到的区域为n 乞我们记第n 个元为e ，i ，它的边界记为钯n 如果是n 乞边界上的元( 的边界是a q 乞的一部分) 且a n 鼬。不在产上，则把这部分边界记成q n ，否则7 n = 0 记巩= v h g ( 彘) ： 在每个e 上 p h 是线性函数，蜥i 一= o ) ，这样( 2 2 3 ) 式可以写成，求鲰一z h g d u h ， 使得 登厂，守”v t ，h d z d y - - - 釜g v h d 8 ， 玩7 ) 卜守” ，玩7 ) n - = l 。厶一 下面戎们来形成有限元方程 设单元e 。为b 毋f m ，在其三个顶点上t h ( $ ，”) ， h ( $ ，) 的函数值分别为啦，畔，“。 和地，哆，并记向量t = h ，吩，t 训7 ，= h ，吩，于是在单元岛上有 a v v u h d t , d y = 一( v ) r ( v “一) d z d u = ，( b y e ) t ( 口。) 蛐= t 嚣耳。u “ 其中k 为单元刚度矩阵 ， 一匙雠忙：r 不妨设为p - z 。设j p - ：e 弓i = z ，参数t 为再巧上的线段弧长，于是可以得到 t i 再巧= 1 一 ，i 丽= ，i 研= 0 可以计算 g v h d 3 ；( n v e 。) t g d 8 = 咳k ， 盖盖 舞( 1 一 ) 拙) 1 其中如= ，勺幽；1 名 删 l 为单元荷载向量 “ 【0j 甜帅彩舳彩 = = = m 以 q 分别把每个单元刚度矩阵和单元荷载向量合成到总刚度矩阵和总荷载向量 中，由式( 2 2 7 ) 可得 u t k u = 斧f 其中口= ( ”l ，啦，) t ，“= ( u l ，t | 2 ，螂；) t ，n p 为节点数，k 和f 为总刚度矩阵 和总荷载向量 于是可得有限元方程为 耳 一f = 0 由于我们考虑的定解问题( 2 2 2 ) 在一的函数值是已知的，即要求一z d v a 在程序编制上是通过对刚度矩阵的装配来实现此要求的对这一都分的第一边值 情况“= 面，对某个“= 砺，我们的处理方法是对总刚度矩阵耳令 j ；1 ，2 ，8 1 ，8 + 1 ，p l = 1 ，2 ，o 一1 ，j 十1 ，n p ，2 1 2 8 ) 对总荷载向量f 令 墨2 0 一 2 1 , 2 ，- s - l ，5 十1 ，n p ，( 2 2 1 9 ) le = 舀 最后就得到了同题( 2 2 2 ) 的有限元方程，求解就得到了( 2 2 2 ) 的数值解 2 2 2电位q 和磁通量密度岛的计算 现在就可以用上一小节的结果求解我们的带有积分边界条件的的m r e i t 的正 问题首先由上一小节叙述的有限元法在每个截面n 。求解带有混合边界条件的同 题( 2 2 2 ) 得到解蟊，再由= f ；残f i + c 得到正问题( 2 2 1 ) 的解，即为电势分布仉 在计算“的时候，需要计算豇在边界的法向导数我们的计算方法是用三次b 一样条数值微分方法求解，其基本思想是和后面要用到的五次样条求二阶导数的 方法是一样的 对我们