（计算数学专业论文）电阻抗断层成像与核磁共振电阻抗断层成像中的若干数学问题.pdf

中文摘要 摘要 本篇论文主要讨论电阻抗断层成像( e i t ) 与核磁共振电阻抗断层成像 ( m r e r d 中的若干数学问题。在第一章中，我们介绍电阻抗断层成像的数学模 型和应用前景，并简要回顾它在理论和数值算法上的发展现状。最后，对于分 片电导率的情形，我们提出了基于水平集方法的数值算法，通过数值例子说明 算法的准确性和稳定性。第二章主要对核磁共振电阻抗断层成像方法的相关数 学问题进行讨论，我们证明了介质电导率为正交各向异性时的某种唯一性，此 外，对于正交各向异性的情形提出相应的重构算法。通过数值例子说明，我们 所采用的间断识别方法和数值微分方法是有效的。 关键词：电阻抗断层成像，核磁共振电阻抗断层成像，反问题，唯一性，水平 集方法，重构算法 中图法分类号：0 2 4 1 8 2 一一 英文摘要 a b s t r a c t i nt h i s t h e s i s ，w em a i n l yd i s c u s ss o m em a t h e m a t i c a lp r o b l e m si ne l e c t r i e a l i m p e d a n c et o m o g r a p h y ( e r r ) a n dm a g n e t i cr 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 ( m r e i t ) i nc h a p t e r1 ，w ef i r s t i yf o r m u l a t et h em a t h e m a t i c a lm o d e l so fe i t , a n dr e v i e wb r i e f l ys o m et h e o r c c i c a la sw e l la sn u m e r i c a lw o r k so nt h ep r o b l e m t h e n w ep r o p o s ean u m e r i c a lm e t h o db a s e do nt h el e v c ls e tm e t h o dt or e c o n s t r u c td i e c e w i s e c o n s t a n tc o n d u c t i v i t i e s n u m e r i c a le x a m d l e sa r cm a d et oi l l u s t r a t et h ep e r f o r m a n c eo f o u ra l g o r i t h m i nc h a p t e r2 ，w ed e a lw i t ht h em a t h e m a t i c a lp r o b l e mo fm p e i t w e p r o v eau n i q u e n e s sr e s u l to fm r e i tp r o b l e mw i t ho r t h o t r o p i cc o n d u c t i v i t yu n d e rt h e h y p o t h e s i st h a tt h er a t i o so f c o n d u c t i v i t i e sa r ck n o w n m o r e o v e r , an u m e r i c a la l g o r i t h m i sp r o p o s e dt or e c o l l s t l m c tt h eo r t h o t r o p i cc o n d u c t i v i t y w ed e m o n s t r a t eb yn u m e r i c a l e x a m p l e st h a tt h en u m e r i c a ld i f f e r e n t i a t i o na sw e l la st h ed i s c o n t i n u i t yd e t e c t i o nm e t h o d sa l ee f f e c t i v ef o rt h er e c o f l s t r n c t i o n k e yw o r d s ： 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 , 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 ，i n v e r s ep r o b l e m ，u n i q u e n e s s ，l e v e ls e tm e t h o d ，r e c o n s t r u c t i o na l g o r i t h m c h i n e s el i b r a r yc l a s s i f i c a t i o nn u m b e r ：0 2 4 1 8 2 一m 一 别舌 电阻抗断层成像( 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 ) 是通过物体表面 的电压和电流来重构内部的电导率。这种技术具有很多优势和广阔的应用前 景，目前它在医学、地质勘探和无损探伤等方面有重要的应用。自从1 9 8 0 年 c a l d e r 6 n 对这个问题进行数学上描述以来 9 1 ，e i t 已被广泛地研究，并取得了 长足的进展。如今这仍旧是个非常活跃的研究领域，因为有许多理论问题和数 值重构算法都具有很大的挑战性。对于实际应用来说，最为关心的是能否有一 种有效、快速的算法来重构电导率。由于问题本身的不适定性与非线性性，目 前重构算法的分辨率还不理想。 在本文中，我们提出基于水平集方法的重构算法。水平集方法是o s h e r 和 s e t h i a n 于1 9 8 8 年提出【3 3 】，现在已经是模拟界面运动的一种有力工具，但关 于它在反问题方面的应用研究才刚刚开始。s a n t o s a 最早将此方法用于反问 题 3 7 1 ，此后很多学者相继展开研究 1 i ，1 2 ，1 6 】。本文主要将这种方法应用e i t 的一个简化模型，即电导率为分片常数的情形，并且通过数值例子说明它的准 确性和稳定性，与此同时，我们也指出该方法的不足之处。 为了克服e i t 问题的不适定性和非线性性，近几年来韩国的一个研究小组 开始致力于寻找一种新的成像方法，这就是核磁共振电阻抗断层成像( m r e i t ) 提出的背景。很多学者都已经对m r e i t 进行了研究，包括直接重构算法和迭代 算法 2 l ，2 3 2 4 ，2 7 】，从数值例子和实际实验可以看出，这些算法很大程度上改 善了e i t 的重构图像。 然而，至今为止关于m r e i t 的理论和算法大部分都是基于“电导率是各 向同性”这一基本假设，而绝大多数的生物体组织的电导率却为各向异性。因 此，为了使得这种成像方法能用于实际的临床实验，需要对各向异性m r e i t 的 唯一性做出回答，并有相应的数值算法来准确并稳定的重构各向异性电导率。 另一方面，我们发现在实际中很多各向异性物体都可以近似地看成是正交各向 异性的导体，即电导率里除对角线外全为0 。本文主要探讨这种假设下的唯一 性和重构算法。我们将证明如果电导率的比率已知，那么在至多两组输入电流 的条件下，电导率可以被唯一性重构。此外，对于正交各向异性的情形，我们 提出相应的重构算法，并通过数值例子说明所采用的间断识别方法和数值微分 方法是有效的。 一一 论文独创性声明 本论文是我个人在导师指导下进行的研究工作及取得的研究成果论文中 除了特另4 加以标注和致谢的地方外不包含其他入或其它机构已经发表或撰写 过的研究成果其他同志对本研究的启发和所设的贡献均已在论文中作了明确 的声明并表示了谢意 作者釜名：跑7 日期：三12 ：竺口l 论文使用授权声明 本人完全了解复旦大学有关保留、使用学位论文的规定。即：学校有权保 留送交论文的复印件允许论文被查阅和借阋：学校可以公布论文的全部或部 分内容可以采用影印、缩印或其它复制手段保存论文保密的论文在解密后 遵守此规定 作者签名 日期：坦! z 鱼少躲獬 删 第一章电阻抗断层成像( e i t ) 1 1 引言 第一章电阻抗断层成像( e l t ) 电导率是材料的一个重要物理特性易导电物体的电导率较高，且极易通直 流电和交流电，相反绝缘体电导率通常低于l o “西门子，而且只有在通交流电时 才有电流通过物体。不同的物体都会呈现出不同的导电性质，如表1 11 4 ，2 2 ，所 以可以根据不同的电导率来推测物体内部的结构。 表1 1 人体组织、岩石以及流体的导电性 人体组织 1 口( a ) 岩石和流体 1 口( 1 2 ) 肺部 9 5 0 海砂 卜l o 肌肉7 6 0 花岗岩 5 0 0 一2 0 0 0 心脏 6 0 0 油田中的氯化物 o 1 6 脂肪 1 0 0 0 油田中的硫酸盐i 2 电阻抗断层成像( 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 ) 正是基于这一物理 现象发展起来的。通过给物体注入小的安全电流，我们可以测量表面的电位势 来重构物体内部的电导率或者电阻抗的分布图像，从而来确定物体内部的结 构，如图1 11 1 0 】。e i t 技术具有很多优势和广阔的应用前景。研究表明某些人体 组织病理改变( 如癌变) 能引起组织电导率的变化，这些信息将会在e 1 t 图像中 体现出来，目前它在医学方面应用的例子包括乳腺癌探测 1 0 l ，心脏功能的监 测【2 0 等。此外，e i t 在材料的无损探伤和地质勘探方面有重要应用，比如管 道腐蚀、地下油罐渗漏的探铡以及矿床和地下水的寻找 3 4 _ 3 6 】。 图1 1 ( a ) 电阻抗断层成像的实验模型( b ) 重构图像 第一章电阻抗断层成像f e l t ) 在过去的二十几年中。e i t 被广泛地研究，并且在理论和实际应用的问题 都取得了长足的进展。如今这仍旧是个非常活跃的研究领域，因为有许多理论 问题和数值重构算法都具有很大的挑战性。在这一章的后续小节，我们将分别 介绍e i t 的数学模型，简单回顾它的唯一性问题，最后，我们对一种分片常数 的情形提出基于水平集方法的重构算法。 1 2 数学模型 1 2 1 正问题的连续模型与电极模型 假设有一导电体qc p ( n = 2 ，3 ) 的电导率为口，如果在边界注入电流g 那么电位势“满足边值问题 v 和v 札) = o 在n 内，嘉= 夕a o 上 ( 1 1 ) 其中礼外法向方向。当盯l ”( q ) ，并且g h 一 ( a q ) ，那么上述n e u m a n n 问 题在相差一个常数的意义下有唯一解“1 ( q ) 【18 】其中一 ( a q ) 与1 ( q ) 都 是标准的s o b o l e v 空间【l 】。通常情况下我们取解u 满足u d s = 0 。 ( 1 1 ) 是e i t 的连续模型，它可以由m a x w e l l 方程推导得到【8 】。它的不足之 处在于要求边界上每一点的电流g 都己知，但在实际实验中，我们只知道边界 电极上的电流，所以需要一种更为合理的数学语言来描述。下面叙述的模型是 s o m e r s a l o 等所建立的更加精确的完备电极模型 3 8 】 在完备电极模型中，设m 个电极所占区域为最，z = 1 ，一，m 。那么在没有 安放电极的地方，即x a n u 圣l s f 时，无法向电流，即嚣= o ；在电极处， 法向电流满足 卜尝d s ：h ， ( 1 2 ) 厶盯瓦幽2 ( 1 2 ) 而电压u 满足 “+ 盈盯褰= v t ( 1 3 ) 咖 其中m 是电极处的电压，匀称为表面阻抗，它考虑了电极和物体之间的化学反 一2 一 第一章电阻抗断层成像( e i t ) 应，该反应的结果是在物体之间形成一层很薄的高电阻层，该层阻抗正是由句 来刻画。最后，由电荷守恒定律和固定参考电压，电极上的电流和电压满足 ( 1 4 ) 虽然这个完备电极模型比较合理，但是相对连续型模型1 1 来说它要复杂的多。 基于这个模型的反问题在理论上和数值重构算法方面还没有本质上的突破，因 此在这篇论文中我们只讨论连续模型。 1 2 2 反问题模型 在e i t 中电导率盯未知，我们希望通过测量边界上的电压u 和输入电流g 来重构它的信息。电压u 和电流夕之间的关系包含了电导率的信息，通常可以 用映射来建立它们之间的关系：d i r i c h l c t t o n e u m a n n ( n t d ) 映射。 定义1 1 ：d i r i c h l e t - t o - n e u m a n n ( d t n ) 映射a ，：日 ( a q ) 一h 一( a q ) ，且 埘( x ) 水) 掣当x m ( 1 5 ) 其中f 日( a q ) ，u 是边值问题 v ( a v u ) = 0 在q 内， = ，在a q 上 的解。 值得一提的是映射如具有一些基本性质【8 】： ( 1 ) 非平凡零空间：n u l l ( a ，) = ff ( z ) 在a q 上为常数) ( 2 ) 自共轭：( a ，g ) = ( ，a ，夕) ，害e 中( a ，f ，g ) = 厶n a ，f ( x ) g ( x ) d s ； ( 3 ) 变分原理：对任意，日 ( a q ) ，( “，) = r a i n 。j 。：，矗盯i v “1 2 d x 对于反问题模型的定义始于c a l d e r 6 n ，他在【9 】中的叙述如下： 定义l ，2 ：给定d i r i c h i l e t - t o - n e u m a n n 映射& ，求电导率函数口( x ) 。( q ) 满足 矿( x ) 0 ，v x q 一3 一 o l i 所 m | i “ m 第一章电阻抗断层成像( e r r ) 实际中通常我们不直接用d t n 映射，而是用它的逆映射( a ，) ( 。) ，这是因 为( a 。) 卜1 ) 光滑性比也好，但两者在理论上是等价的。另一方亟，一般来说实 际中我们无法得知( a ，) ( 以) 或者“的所有信息，通常只能做组实验，每组实 验给定特定的激励电流g t ，并测得相应的边界电压，l 。这是理论与实际应用的 差异之处，所以更为合理的反问题模型定义为： 定义1 3 ：给定n e u m a n n t - t o d i r i c h l e t 映射( 虬) ( 一1 ) 的部分信息，重构电导率函数 盯( x ) l o o ( q ) 本篇论文有关理论上的结果都是假定知道屯的所有信息。即定义1 2 ，而对于数 值重构算法都是采用定义1 3 。 1 3e i t 问题唯一性简单回顾 给定d t n 映射k ，是否能唯一重构物体q 的电导率矿? 这就是e i t 的唯 一性问题。这一节简单回顾了关于唯一性的一些结论，包括盯在边界上的唯一 性重构与物体内部的唯一性重构。 唯一性理论的第一个重要结果是由l ( 0 l l i l 与v o g e l i u s 提出【2 5 1 ，他们证明 了如果区域边界是g * 的，盯l e e ( f 2 ) 而且在边界附近口g * ，那么d t n 映 射k 唯一决定了边界上的电导率口。为了证明在x o a q 的唯一性，k o t m 与 v o g e l i u s 构造了高度振荡函数v 日m ( q ) ，并且矿在x o 的一个小邻域外为0 。 这样的边界条件就保证了边值问题 v ( c v u ) = 0i n q ，t l = fo na q 的解在殉附近快速趋近于0 ，因此可以在】c 0 处可以对电导率盯进行讨论。 s y l v e s t e r 和u h l m a r m 4 1 】通过微局部分析的技巧给出了盯在a q 上的显示表 达式。证明是基于假设“区域q 为g ”而且矿c 。( q ) ”，因为这个假设保证 了为一阶拟微分算子。设x 为x o a q 附近的一个点，令x x o = ( z 7 ，z n ) ， 其中z 7 p 那么a ，的象征就是 a ，( z 7 ，7 ) = 盯( z 7 ，0 ) j 亭7 i + a ( z 7 ，专) + r ( z ，毒7 ) ， 其中a ( z 7 ，关于7 是零阶算子，由口在( z ，o ) 的一阶切向和法向导数所决 定。r ( z 7 ，) 是1 阶算子，由盯在( z ，0 ) 的高阶切向和法向导数所决定。因此 4 - 第一章电阻抗断层成像( e r r ) 盯在边界的重构公式就是 巾，0 ) _ i 器箐舻_ ( 1 6 ) n a c h m a n 【3 0 】也给出了口在边界的重构公式，因为他在证明中并没有要求 人。是一个拟微分算子，所以对区域q 和电导率仃的光滑性要求降低了许多，实 际上在a q 为c 1 ，1 并且盯c 1 , 1 ( q ) 的光滑性要求下，就有重构公式 i 概f f 2 e - i k s o a c e i x k - - 1 7 ( z ) 恢啦2 o ，( 1 7 ) 其中s o 为单层位势算子岛，( x ) = g o ( x ，y ) f ( y ) d s ，g o & l a p l a c e 方程的格林函数。 上述结果讨论的都是关于边界上g r 能否被唯一重构，对于物体内部的唯一 性问题，我们有 定理1 4 ：如果qcr 3 为有界的单连通区域，并且a q 具有e * 光滑性，假设 0 1 和o 2 是正的e 。( q ) 函数。那么a n = a ，：寺在q 上0 1 = o 2 。 这个定理的证明是由s y l v e s t e r 和u h l m a n n 4 0 】给出的，主要的证明思想是把 e i t 问题转化为s c h r o d i n g e r 方程 妒一q 妒= 0 ， 其中庐= t t l 2 u ， q = 盯1 2 盯1 2 通过证明逆散射问题的唯一性来说明口可以被儿唯一确定，因此由椭圆型边值 问题 盯m 一妒m = 0 ， 盯在a q 上给定 ( 1 8 ) 的唯一性可以推出仃被a 。唯一确定。这是三维问题的结论，对于二维问题有同 样的结论，n a c h m a n n 【3 1 】在1 9 9 6 年给出它的证明。 1 4 重构电导率间断界面的一种水平集方法 1 4 1 引言 在这一节，我们考虑e i t 的一个简化模型。假设在( 1 1 ) 中qcr 2 ，电导 一5 一 第一章电阻抗断层成像( e i t ) 图1 2 ( a ) 电导率为分片常数的简化模型( b ) 界面运动轨迹 率是盯是分片常数，如图1 2 ( a ) 所示，即 f 矿x d i 盯一。a 一d 一 进一步假定区域d 的电导率矿d r 一。这种情况在提炼油的时候会出现， d 表示水所占的区域，而q 西中为油。那么根据实验发现，在d 中近似有 t 0 4 5 】。另一方面，我们假设边界的激励电流为g ，相应的在边界所测量的 电压为，。那么简化后的模型为 甚纂1 a q 上 ( 1 9 ) 对于实际应用来说，最为关心的是能否有一种有效、快速的算法来重构电 导率的间断界面o d 。至今为止，对于e f t 问题的数值重构算法主要有两类： 直接法和迭代法。直接法包括由h a n k 和b r i i h l 发展的分解算法【5 - 7 】和基于 n a c h m a n n 对e l f 唯一性证明的重构公式 2 8 ，3 9 】，这些算法的优点是速度快， 但最大的缺陷是要求k 的所有信息都给定，这给实际应用造成了很大限制。迭 代法包括最小二乘法【1 5 ，1 6 ，3 9 】和勋l l i l 提出的变分法【2 6 ，它们是基于非线性 优化问题的求解。由于问题本身的不适定性与非线性性，目前的这些算法的分 辨率还不理想。 这一节我们探讨基于水平集方法的重构算法。水平集方法是o s h e r 和 一6 一 第一章电阻抗断埕成像( e i t ) s e t h i a n 于1 9 8 8 年提出【3 3 ，最初的目的在于追踪不同流体界面的运动，如今 它已经是模拟界面运动的一种有力工具，特别是界面发生拓扑变化的情况，在 图像处理和计算物理方面有重要应用【3 2 ，但关于它在反问题方面的应用研究 才刚刚开始。s a n t o s a 最早将此方法用于一个积分方程反问题 3 7 】，此后很多 学者相继对它在偏微分方程反问题的应用展开研究，如d o r a 等对逆散射问题 提出水平集方法 1 6 】，而t o n yc h a n 对类似( 1 1 ) 的问题也提出相应的水平集算 法【1 1 ，1 2 】，但与我们所讨论的问题不同之处在于，他要求给定u 或者v u 在q 内的值，并且将问题直接转化为关于水平集函数西的优化问题。这种方法的不 足之处在于对同样一个区域d ，可能对应于不同的水平集函数西，这会给优化 问题的求解带来困难。在这里，我们则是通过给d 。在t 时刻的边界o d 。选取适 当的速度v ( x ，t ) ，满足 磊d ( 夙) 一川：( 删 0 是正则化参数，也就是说把边界o d 的长 一7 一 第一章电阻抗断层成像( e i t ) 度作为惩罚项。这种正则化方法在图像处理中经常使用【2 9 。从某种意义上来 讲，上述问题是一个形状优化闯题( s h a p eo p t i m i z a t i o np r o b l e m ) 。先考虑它的子 问题，即o = 0 的情形， r l 【d ) ：= 丐1min m i n j j 廿( d ) 一，各( 触) ( 1 1 1 ) r l 【d ) ：= 写j j 廿( d ) 一，刍触) ( 1 1 1 ) 要求解上述优化问题，我们假设初始时刻t = 0 猜测为风，并且在t 时刻 所求的解为d 。，如图1 2 ( b ) 所示，那么按照优化理论，通常希望找一个 “下降方向”，或者等价地，给觑的边界o d c 选取一个速度v t ( x ) ，使得 d 磊r l ( d t ) 0 。当然v t ( x ) 的选取方法有很多，我们所选取的速度如下述定理所 述： 定理1 5 ( 速度的选取) ：设d tcq ，n ( t ，x ) 为功的外法向向量。如果令v t ( x ) = u ( ，x ) 礼( t ，x ) ，其中u ( ，x ) = v u l ( t ，x ) v u 2 ( t ，x ) ，且u l ( t ，x ) 、u 2 ( t ，x ) 分别满足 与 当x q 夙 当x o d ，( i 1 2 ) 当x a q 。t 上2 ( ，x ) = 0当x n d t 抛( t ，x ) = 0当x o d t ，( 1 1 3 ) 豪啦x ) = 州淞) 一舷) 当x m d 其中n 为d t 的外法向导数，那么磊r 1 ( d t ) 0 在证明这个定理之前，我们需要一个引理。 引理1 6 ：【4 3 】设t i ( t ，x ) 满足边值问题( 1 1 2 ) ， i , 1 ( t + 乳，x ) 满足 a 。l ( + f 5 t ，x ) = 0当x q 西件疵 牡l ( t + 疵，x ) = 0当x c o d t + 缸， 未“m + 阮x ) = 出) 当x m 一8 一 x 0 “ = 0 = 、，、， x = x 船神m u 屯 u & “旦饥 ，、l 第一章电阻抗断层成像( e i t ) 我们定义e u l e r 导数似，x ) = l i m 以埘竺堡坐弓掣，那么船，x ) 满 足边值问题 i 。“7 ( ，z ) = 0 当x a z 3 t c ( t ，x ) + v u l ( t ，x ) v ( t ，x ) = 0 当x o d ， ( 1 1 4 ) 【未让协，x ) = o 当x 讹 定理1 5 的证明： 我们注意到 d f 五r t ( d t ) 2 上n ( “1 ( 。，x ) 一，( x ) ) 牡船，x ) ， ( 1 1 5 ) 因此引入中间变量啦( t ，x ) ，它满足边值问题( i 1 3 ) 。根据g r e e n 公式，有 f 。d t a x u l 锄一“让7 奴= ( 厶筹砌一丽o u 2 u 幽) 一( z 仇筹锄一鬻d s ) 由( 1 1 3 ) ，( 1 1 4 ) ，上式转化为 ( u ，一，) ，如： j m| a d t 丽o u z 憾 一厶鬻c v u l ，v ，如 注意到1 1 9 d r = “, 2 l o d t = 0 ，v = 吼= ( v n l v u 口) n ，那么 ( 钍l g ) c d s = 一( v u l v 地) ( v n ) d s 0 ( 1 1 6 ) j a n j o d t 结合( i 1 5 ) 定理得证。 对于问题( 1 1 0 ) 的速度( 下降方向) 选取，我们有下面的推论。 推论1 7 ：设d t ，n ( t ，x ) 如定理1 5 所定义。如果令v t ( x ) = ”( t ，x ) n ( ，x ) ，其中 u ( ，x ) = v u t ( t ，x ) v 地( t ，x ) 一a k ( ，x ) ，托为d 的曲率，且u 1 ( t ，x ) 、坳( 以x ) 分 d 别满足( 1 1 2 ) ，( 1 - 1 3 ) ，那么瓦n ( o t ) 0 9 第一章电阻抗断层成像( e i t ) 这个推论的证明非常简单，只要注意到 4 3 】 那么结合定理1 5 有 丢z 仇- d s = z 皿一c v 礼，幽 d r ( 厶x ) = 一z 。( v u x ，v 坳一。k ) ( v ，佗) 如 ( 1 1 7 ) 1 4 3 水平集方法 对于优化问题( 1 1 0 ) ，假定t 时刻速度v ( t ，x ) 如推论( 1 7 ) 选定，我们希望 对于0 d t 上的每个点都以此速度运动。最简单的方法是在o d 。上求解常微分方 程 譬= y ( f ，x ) ( 1 1 8 ) 面2 。i ，x j l 1 这是对界面的一种l a g r a n g c 描述方式。要描述它的运动轨迹，只能对( 1 1 8 ) 进 行数值求解，这就意味着要对o d t 进行划分。如果o d ：的拓扑没有发生变化或 者变形不大，这种方法完全没有问题。但是，如果速度场v ( t ，x ) 会导致边界发 生拓扑变化( 如图3 ( 上) 所示) ，这种描述方式的就完全不可靠，因为这种情况 下对同一个边界的离散点集合可以有不同的连接方式，也就是说，可以有完全 不同的拓扑结构。 图1 3 界面的拓扑变化( 上) 与相应的水平集描述( 下) 一1 0 一 第一章电f r 抗断层成像( e i t ) 为了克服拓扑变化和变形带来的困难，我们用水平集函数( t ，x ) 来表示界 面o d 。的演化，即定义 妒：( 0 ，o 。) r 2 _ r 1 ， 满足o d t = ( ，x ) i 咖( t ，x ) = o ( 1 1 9 ) 通常要求妒在o d t 所包含的区域内大于0 ，在区域外小于0 ( 如图3 ) 。注意到当 界面从左演化到右时，尽管它的拓扑发生了变化，我们只需要把水平集函数西 往下移即可，而不用直接去追踪界面本身的运动轨迹。这种通过水平集函数西 来描述界面运动的方法称为水平集方法，它是一种e u l e r 描述方式，由o s h e r 和 s e t h i a n 于1 9 8 8 年提出【3 3 】。 对于水平集函数西，容易看出它满足偏微分方程 掣+ v ( t ，x ) - v ，x ) = o 这个方程通常称为水平集方程，它是一个h a m i l t o n - j a e o b i 方程。对这类方程的 数值求解，时间方向的离散格式可以用e u l e r 法，空间方向的离散可以用迎风格 式，它的精度为一阶，或者用高精度格式e n o ，w e n o 格式【3 2 1 。 1 4 4 重构算法与数值例子 对于问题( 1 9 ) ，给定边界上的测量值，和g ，结合前两小节的讨论，我们 给出基于水平集方法的重构算法： 1 设定初始时刻的水平集函数护，令c o d o = ( t ，x ) fo o ( t ，x ) = o ) 。 2 对于后= 0 ，1 ，2 ，求解边值问题 憎亍。耋裳筹一眦 z - ， 滢a x u l k 。薹裳筹刊一舭z z ， i 舡。在眠上，筹刊一，在别二 q 2 动 第一章屯阻抗断层成像( e i t ) 如果l i u 1 1 l z ( 船) e ，其中f 是给定的误差精度，那么停止迭代，否则 转至下一步。 3 在区域f l d 上，令速度v 女= k n ，其中仇= v 艟- v 谚。在区域d k 内， 对于每一行网格采用线性插值得到速度v k 。 a ，始 4 求解水平集方程：善+ v k v 扩= 0 。时间方向采用e u l e r 格式离散，即 z 七+ 1詹 口 兰+ v v 扩= 0 ，对于空间离散采用e n o 或者w e n o 格式，其 中时间步长a t 与空间步长x x 的选取应该满足差分格式的稳定性条件。 5 如果需要，对西。十1 进行周期性重新初始化( r e i n i t i a l i z a t i o n ) 。 6 求扩+ 1 的零水平集o d k + l ，即o d k + l = ( t ，x ) i + 1 ( t ，x ) = o ，。并转到 第2 步。 在每步迭代过程中，第2 、3 步通过求解正问题( 1 - 2 1 ) ( 1 2 2 ) 来确定边界 c g d k 的速度v ，这是因为这种速度选取方法可以保证i i u 一，儿弘旧氆随k 单 调下降。另一方面，我们引入水平集函数扩，通过求解h a m i l t o n - j a c o b i 方程 ：善+ v k v 扩= 0 来模拟边界的运动，这正是水平集方法。此外，迭代停止 的准则是判断让f 与测量值，之间的差的l 2 范数是否小于给定精度e 。值得一 提的是，在数值求解h a m i l t o a - j a c o b i 方程后( 第3 步) ，如果需要，我们对咖o + 1 进行周期性重新初始化，它的目的在于让西1 尽量接近一个有符号的距离函 数( s i g n e dd i s t a n c ef u n c t i o n ) 。这是因为，我们最为关心的是扩+ 1 = 0 的点的集 合o d k + 1 ，而扩“本身的值并不那么重要，而转化成一个有符号的距离函数 ( s i g n e dd i s t a n c ef u n c t i o n ) 可以保证o d k + l 能够更为准确地求出。关于重新初始 化的具体实现可以参考 3 2 】的第七章。 对于上述重构算法，我们希望通过数值例子来说明它的优点和不足。下面 的例子将分别考虑区域d 为单连通、多连通的情形。以及测量数据带误差和不 带误差的情况。在所有例子中，我们都假设区域n ；( 0 ，1 ) ( 0 ，1 ) ，在边界上的 激励电流g = 1 ，边界上的电压值，通过正闯题求解得到，正问题的求解基于 二阶有限元方法。 例1 8 ( 单连通区域) ：真实区域d 的边界是中心在( 0 4 ，o 5 ) 、长短半轴分别为 0 2 与o 1 的椭圆，如图中实线所示，虚线所示为数值解( 下同) 。我们假设初始 一1 2 一 第一章电f r 抗断层成像( e t ) 猜测d o 为以( 0 5 ，o 5 ) 为圆心、半径是o 4 的圆。可以看到当k = 5 0 时，数值解 玖和真实解已经d 十分接近，丽当k = 2 0 0 时，两者的边界基本吻合。因此， 基于水平集方法的重构算法在处理d 为单连通区域时精度高，速度快。 ( c )( d ) 图1 4 ( a ) 初始猜测d o ( b ) 岛= 5 0 的重构边界形状( c ) k = i 0 0 ( d ) k = 2 0 0 例1 , 9 ( 不连通区域) ：真实区域d 由两部分d 1 和d 2 组成，其中d 1 为圆，其 半径为0 , i ，d 2 是长、短轴分别为0 i 与0 0 5 的椭圆。初始猜测 o 与例i 8 相同，我们发现当女= 2 0 0 时，d 的拓扑与d o 相比已经发生了非常大的变 化。k = 5 0 0 时，图( b ) 中集中在中间的小区域已经消失，而且数值解逼近真实 解。之后，速度v 的速度非常小，所以收敛较慢，七= i 0 0 0 时，数值解基本 上能勾勒出j d 的形状。 第一章电阻抗断层成像饵i t ) ( c )( d ) 图1 5 ( a ) 初始猜测d o ( b ) k = i 0 0 的重构边界形状( c ) k = 2 0 0 ( d ) k = 5 0 0 例1 1 0 ( 不连通区域) ：真实区域d 由三个部分d 1 、d 2 和d 3 组成，与例1 9 不 同的是它们的距离较近。我们注意到当k = 1 0 0 时，a d k 在d 的外日形成一个 包络，而且k = 2 0 0 ，5 0 0 时，这个包络并没有发生太大的变化，只是在形状上 与d 的三个小区域的分布有些类似。我们在数值实验中发现，即使当k = 2 0 0 0 时，所求的数值解0 d k 仍为一个包络，而且，即使在边界上有多组激励电流、 或者改变电流值，的取值，仍无法准确重构出真实解。我们推测，出现这种情 况的原因是所求的解只是优化问题( 1 1 0 ) 的一个局部极小值点，而不是整体最 小值点，而这是求解优化问题一个本质性的困难，这就说明重构算法在处理距 离较近的不连通区域时分辨率还不够。 一1 4 第一章电阻抗断层成像( e i t ) ( c )( d ) 图1 6 ( a ) 初始猜测d o ( b ) k = 1 0 0 的重构边界形状( c ) k = 2 0 0 ( d ) k = 5 0 0 例1 1 1 ( 初始猜测的选取) ：为了解决上述例子中对于距离较近的不连通区域分辨 率不高的困难，我们希望初始猜测d o 也是个不连通区域，如图1 7 ( a ) 所示。 这时重构图像已不再是真实解d 的包络，它已经基本上能勾勒出d 1 、d 2 和d 3 的轮廓，如图1 7 ( b ) 。因此，这可能是解决不连通区域的有效方法。但是，是 否可以对于其他问题也用类似的初始猜测? 我们讨论例1 8 的情况。不幸的是， 这种选法并不能保证它收敛于精确解，我们所得到的数值结果与真实解完全不 符合，如图1 8 ( b ) 。因此，在不知道真实解d 的情况下，如何选取适当的初始 猜测是一个困难的问题。 一1 5 第一章电阻抗断层成像( e r r ) 图1 7 ( a ) 初始猜测d o( b ) k = 1 0 0 0 的重构边界形状 图1 8 ( a ) 初始猜测d o( b ) k = 2 0 0 的重构边界形状 例1 1 2 ( 稳定性) ：这个例子考虑数据带误差的情况。对于例1 8 、1 9 ，我们把有 限元方法求解出来的边界电压g 分别加上1 0 的均匀误差。在数值计算中，取 正则化参数口= 1 0 一，重构的结果图1 9 ( b ) ，1 1 0 ( b ) ，可以看到这种算法是稳定 的。另外，我们发现正则化参数无需选取得太大，这也确保在求解水平集方程 时步长a t 不会太小 图1 9 ( a ) 初始猜测d o( b ) = 2 0 0 的重构边界形状 一1 6 一 第一章电阻抗断层成像( e r r ) 图1 1 0 ( a ) 初始猜测 o ( b ) k = 1 0 0 0 的重构边界形状 通过以上的数值例子，我们可以看出，基于水平集方法的重构算法在d 为 单连通区域、或者由距离较远的几个部分组成不连通区域时非常有效，而且数 值稳定性高。但是，如果对于例1 1 1 的情况，即d 1 、d 2 和d 3 距离较近，我 们还无法准确重构真实值，这是由求解优化问题( 1 1 0 ) 本身的困难所导致。选 择合理的初始猜测可能是种有效的办法。但在不知道真实值d 的情况下，如 何选取初始猜测仍是个困难的问题。此外，还需要提到的是这种方法属于迭代 法，而且每次迭代需要求解两个正问题，所以相对直接重构算法来说，它的速 度是一个劣势。用边界元来求解正问题可以提高算法的速度，但这并不能解决 问题的本质，我们希望能找到一种更为合理的速度来v 减少迭代次数，也就是 说，在求解优化问题( 1 1 0 ) 时需要寻求更为合理的“下降方向”。 一1 7 第二章正交各向异性核磁共振电阻抗断层成像( o m r e i t ) 第二章正交各向异性核磁共振电阻抗断层成像( o m r e i t ) 2 。1 引言 在第一章中我们已经提到，由于传统e i t 问题的不适定性和非线性性，对 于它的重构算法分辨率都比较差。近几年来韩国的一个研究小组开始致力于寻 找一种新的成像方法，从而来克服上述困难，这就是核磁共振电阻抗断层成像 ( m r e r r ) 提出的背景。很多学者都已经对m r e i t 进行了研究，包括直接重构算 法和迭代算法【2 1 ，2 3 ，2 4 ，2 7 1 ，另外，也有些重构算法则直接利用导体内磁感应 强度b 在z 方向的信息b ： 4 2 1 ，这是由于在实际实验中，b 的另外两个方向 往往很难测量到。从数值例子和实际实验可以看出，这些算法很大程度上改善 了e i t 的重构图像。 然而，至今为止关于m r e i t 的理论和算法大部分都是基于“电导率是各 向同性”这一基本假设，而绝大多数的生物体组织的电导率却为各向异性。因 此，为了使得这种成像方法能用于实际的临床实验，需要对各向异性m r e i t 的 唯一性做出回答，此外，仍没有相应的数值算法来准确并稳定的重构各向异性 电导率。另一方面，我们发现在实际中很多各向异性物体都可以近似地看成是 正交各向异性的导体，即电导率堡除对角线外金为0 ，这种情况下问题的处理 相对不是那么复杂，这一章主要探讨这种假设下的唯一性和重构算法。我们发 现这时会得到与各向同性m r e l t 类似的一些结论和方法。在2 2 中，我们将证 明如果电导率的比率已知，那么在至多两组输入电流的条件下，电导率可以被 唯一性重构，2 3 是正交各向异性电导率的重构算法，我们通过数值例子来说 明算法的准确性和稳定性。在进行这些讨论之前，先简单介绍一下o m r e i t 的 模型。 假设导电物体q cr 2 的电导率里为正交各向异性，即 ，0 、 驴io 观。| ( 2 1 ) 。e 1 ( 囝& = l ，2 ) ，且恒正。边界o f t 上的电流g 通过电极p 和q 输入输 出( 如图2 1 ) ，我们假定每个电极是一段长度为2 e 的弧，所以电流9 可以通过下 一1 8 一 第二章正交各向异性核磁共振屯阴抗断层成像( o m r e f r ) 图2 1 边界a q 上的电极p 、9 以及电流9 输入方式 面的式子来表示，即 gcx，=gi，尸；qjcx，=；羔 同第一章，导体的电位势t 满足边值问题 x j x p i e n o g f 2 ( 2 2 ) x l x o i 0 是g 1 ( 而) 中的一个已知函数，它表示0 - 1 1 和0 2 2 之间的比率。另 外，0 - e 1 ( q ) 恒正。因此，在这种情况下，正交异性电导率0 - 中只有0 - 未知。 我们将证明电导率具( 2 5 ) 形式下问题( 2 4 ) 的唯一性。值得一提的是，这个结 论其实是对 2 4 】中唯一性定理在正交各向异性问题上的推广。【2 4 】中作者证明 了各向同性m r e i t 问题的唯一性，也就是( 2 5 ) 中k = 1 的情形，而这里我们 假定片为c 1 ( q ) 中任意的正函数。对于o m r e i t 问题的唯一性，我们主要有下 面的结论： 定理2 1 ( 唯一性定理) ：假定电流9 l 、9 2 如( 2 2 ) 所给定，相应的电流密度 j 1 、j 2 已知。如果( u 1 ，u , 2 ，盯) 与( 面1 ，豆2 ，子) 分别是电导率具( 2 5 ) 形式下边值 问题( 2 4 ) 的解。并且0 - ，矛c 1 ( q ) ，那么t l = 豇l ，坳= 奶，而且0 - = 矛。 一2 0 一 尝等“= 如娑嚣 第二章正交各向异性核磁共振电阻抗断层成像( o m r e i t ) 在证明这个定理之前，我们需要证明下面的引理： 引理2 2 ：假设盯e 1 ( q ) 恒正，吼= g kp j q 】( i 1 ，2 ) 如( 2 2 ) 所定义如果 ( u 1 ，i 2 ) 是电导率具( 2 5 ) 形式下边值问题( 2 4 ) 的