《电阻抗断层成像的求解方法分析》4700字

电阻抗断层成像的求解方法分析综述目录TOC\o"1-2"\h\u10363电阻抗断层成像的求解方法分析综述 1312731.1EIT技术的正问题 15666Σ∂φ∂n=J (2-3) 1306861.1.1EIT正问题的基本实现原理 2142391.2.2有限元法求解EIT 4306352.2EIT技术的逆问题 5284282.2-1EIT逆问题 55251.2.2EIT的图像重建算法 5160141.2.3牛顿-拉夫逊算法 6178361.2.4正则化修正技术 81.1EIT技术的正问题EIT研究具有特殊边界条件的电场，其数学模型可以等效于Laplace方程： ∆∙σ∆ϕ=0 (2-1)加强边界条件：ϕ=ϕ0 自然边界条件： Σ∂φ∂n=J 其中，ϕ是物体内部的电压分布，σ为电导率，ϕ0和J分别为边界区域的电压和电流密度REF_Ref1122394767\n\h[13].根据驱动信号和测量模型的阻抗分布计算内部电流和电压分布，这在EIT的研究中统称为正问题(ForwardProblem)，为了解决正问题，可以通过求解拉普拉斯方程得到该区域的节点电压REF_Ref1119739261\n\h[12]，然后利用已知的条件计算出内部电流密度，从而达到综合分析电场的目的.当前科学研究团队中，用于解决正问题的数值方法是有限元法，有限体积元法，边界元法等.有限元法(FiniteElementMethod，FEM)，通过将离散和连续求解区域简化为有限数量的元素节点，该方法适用于分割复杂字段，该方法更全面，更灵活.边界元法(BoundaryElementMethod，BEM)适用于开域问题.用高斯定理降阶复杂高阶问题，使得待求解方程被转化成边界积分方程，从而将边界细分成有限个单元，这种方法由于Gauss定理，保证了精确度，减少了计算工作量REF_Ref745565020\n\h[9].有限体元法(FiniteVolumeMethod，FVM)，它将场域划分为大量的主单元元素，并围绕单元主元素的每个节点构建一个封闭的控制体，称为次单元元素.电流通量由主单元中每个节点的电势来计算和表示.这样就可以形成一个矩阵方程来求解每个节点的电势.该方法精度高，计算量少，可以容易地解决电导率突然变化的精确模拟问题.图1.2.1-1EIT正问题EIT问题无法通过解析方法解决.这与许多其他非线性逆问题相同.必须借助迭代计算.首先，假设电导率分布为已知的感应区域，解决正问题以获取场边界处的电势分布.根通过基于已知边界条件迭代地校正假定的电导率分布，并进行无限次，直到正问题的结果的误差小于某个值，也就是说，此时的正问题解已经被认定为收敛到满足已知边界条件，从而获得导电率的分布情况，然后在利用现有的图像重建算法对求解出的结果进行成像处理.虽然已经能将这个过程通过有限元法变形，将式(2-1)转化成求解线性方程组，但是，正问题求解过程迭代一直是最大的问题，它会花费很长的时间.有限元法求解过程：开始，场域剖分，取得各种信息(节点总数、强加边界条件节点编号、单元节点编号、三角元总数以及单元电导率)，形成系数矩阵，构成有限元方程，边界条件、激励模式，求解方程，结束REF_Ref1128932690\n\h[14].1.1.1EIT正问题的基本实现原理EIT技术的主要研究问题是EIT正问题计算，EIT逆问题计算和硬件系统设计REF_Ref1151807017\n\h[15].EIT正向问题的数学模型主要从麦克斯韦的电磁场方程开始，推导它，然后使用有限元方法求解.现有研究已经证实，生物组织中的载体材料具有非均匀分布的各向异性，并且生物组织的电学性质不仅仅是电阻，而是电阻和电容，这证实了生物组织中存在着虚部信息.电磁场中的电流和电导率分布之间基础约束关系可通过麦克斯韦方程组进行求解表达： ∇∙L=ρ (2-4) ∇×E=∂B∂t ∇∙B=0 (2-6) ∇×M=J+∂D∂t其中： B=μM ( D=εE (2-9) I=σE (2-10)电位移——L.电场强度——E.磁感应强度——B，磁场强度——M.电流密度——I.电荷密度——ρ。介质的磁导率——μ.介质的介电常数——ϵ.介质电导率——σ.由于电阻抗成像问题中敏感场特性包含实部信息和虚部信息，当向场域内加入频率为𝜔的电流激励时，其电场强度为： E=E可得： ∇×M=σE+ε∂E∂t又因为：E=−φ是场内电势的分布，故上式可以变形为： ∇×M=−γ∇Φ (2-13)将上式两端取散度得到的Laplace方程如下： ∇∙(γ∇Φ=0) (2-14)其中：γ=σ(x,y)+jωε(x,y)，γ是敏感长内复电导率.一旦确定了边界条件，它也可以由上式唯一地确定：对于EIT问题，由于激励源是电流源，因此满足Neumann边界条件，即，给定边界. γ∙∂Φ∂其中，n为边界外法线向量，S为场域边界，I为电流密度.1.2.2有限元法求解EIT当给定了σ和在实际的测量中，因为代替电位来测量物体上各点的电压，所以以O点为基准点，其电势为0.如果将其他边界点和点O之间的电位差设定为电位值.将点O的坐标设为(x0，y0)，则测得边界电势分布函数φ(x,y)(首先将问题化为变分形式.记Ω=Ω∪∂Ω，假设C’C`(Ω)={u|u及uU=u(x,y)是方程(N’)的解，则对任意的v=v(x,y)∈ ν∙(∇∙ρ−1即∇∙(v在Ω上积分后再利用高斯公式可得：Ωρ−1∇ν∙∇udxdy=Ω∇∙(νρ−1∇u)dxdy记D(u,ν)=Ωρ−1∇ν∙∇udxdy，G(ν)=(G)求u∈在原始问题中，需要u的二次偏导数的存在，但是在积分形式中，积分和格林公式的应用大大减少了解的平滑度，ux面对一般的u∈C2.2EIT技术的逆问题1.2.1EIT逆问题的基本实现原理2.2-1EIT逆问题EIT的逆问题是一个病态的非线性问题，由图2.2-1以及2.1节公式2-1，2-2，2-3可知，逆问题就是给出边界的电压测量数值和激励电流反解模型内的电阻抗分布，即由ϕ0，J求解ρ1.2.2EIT的图像重建算法EIT最关键的一环，就是最终的成像，成像的好坏往往决定着它的使用价值.EIT图像重建算法主要是根据边界测量值来确认测量体内部电导率分布情况，随着国际上对EIT技术的关注，EIT图像重建算法逐渐成为EIT技术中最频繁的研究方向.当前，智能算法和非智能算法是EIT图像重建算法的两个主要类别.非智能算法通常通过使用线性模型模拟逆EIT问题来解决，该问题可大致分为迭代算法和非迭代算法.非迭代类算法不需要经历迭代来优化目标函数就可以获得EIT逆问题的解.1985年，Muria等人提出了一种基于Geselowitz的灵敏度算法，以通过灵敏度矩阵对线性方程进行求解，从而简化了问题.1990年，Cheny等人创新出基于最小二乘法的单步牛顿残差重构法.1991年，Somersalo等人提出层剥法.1998年，Vauhkonen提出了对应阻抗分布先验假设的Tikhonov正则化矩阵的构造法，等等.迭代算法是通过更新参数的连续迭代来使目标函数达到最小值的算法的一种REF_Ref1161487849\n\h[16].1987年，Yorkey提出了改进Newton-Raphson法，并提出用标准方法和补偿定理计算雅可比矩阵.1991年，Huda率领众人基于改进的Newton-Raphson法，提出一种先验信息融入图像重建的正则化法，等等.而在具体的实现里，众多科学家都在致力于研究出一个有效的求解器.目前，EIT的大多数研究小组使用的是在MATLAB中电阻抗层析成像和漫反射光学层析成像重建软件EIDORS.EIDORS提供了一组有用的性质，例如二维和三维正演模拟，以及一组广泛的重建算法、可视化功能等等.Horesh等人使用不同的预处理条件和更有效的例程对EIDORS进行了调整，产生了一个称为SuperSolver的版本，目前仍在UCL的小组中使用.然而，对于大网格，MATLAB缺乏高效的并行编程能力，这使得计算正解成为一项冗长的任务.而后，Borsic等人将Jacobian计算(但不是系统矩阵的集合)转移到稀疏并行直接求解器库PARDISO中，以超越这些限制.与Horesh等人相比，他们能够将正向模拟的速度提高约5.3倍.他们把它用在含有大约50万个元素的网格上.在较大的网格上，直接求解需要大量内存，这些内存通常会限制可计算的网格大小.此外，在本文中，直接求解器的装配速度要比一个好的预处理器慢得多，这使得迭代方法的执行速度更快，而这取决于唯一电流注入模式的数目.特别地，代数多重网格预处理程序已经被证明可以显著地提高求解时间.基于图形处理单元(GPU)的计算已经成功地应用于对内存的快速访问至关重要的Jacobian矩阵的计算中.1.2.3牛顿-拉夫逊算法当使用牛顿迭代法作为EIT图像重建算法时，它基于最优化理论的思想，以重建模型边界的电压的计算值与测量得到的边缘电压的值两者间差的二范数最小，即得到一个合适的静态电导率分布REF_Ref1164244197\n\h[17]σ∗，使得(2-18)最小化 f(σ)=1其中σ作为电导率分布是一个非线性函数.F(σ):Cn→假设T(σ)=F(σ)−V，由于电阻抗是非线性的，有牛顿迭代法，在点( T(σ+h)=T(σ)+T`采用线性最小二乘法，找到适合的h，让目标函数达到最小，即T(σ+h)≈0， T(σ+h)=T(σ)+T`(σ)h≈0求解： h=−(T`(σ))−1T同时： h=σk+1−σk=[F即： σk+1=σk+[F`基于之前的理论，更异于使趋近于0，也就是可以使得目标函数达到最小值.而通常就是Jacobian矩阵，定义为[F`由上述推导可知，最终目标是最小化的F(σ)−V，但是，在展开之后，不计算高阶导数项，从而引起较大的误差，导致最终结果的质量.所以公式(2-18)是最小化的研究对象，令∇f(σ)= ∇f(σ)=对f(σ) f``(σ)=F`(σ)(2-25)中的F`(σ)TF将公式(2-18)在(σ+h)处展开成泰勒级数，可得： f(σ+h)=f(σ)+f`求上式的梯度并将其设为0，就可以使其最小化，可得： h∙∇f=求解得： h=(F`(σ)T综合上述推导，可得： h=(F在上述的推导中，本文还是做了简化计算，例如不计算高阶变量，在实际应用中，由于Hessian阵的条件数一般都是比较大的，所以这样的处理会使得误差较大，对最终结果影响也很大.因此，还要引入正则化法来修正.1.2.4正则化修正技术Hessian矩阵的广义逆条件的数目较大，导致该解的病态需要一些额外的先验和解的附加约束条件，为了将不适当的问题转变为适当的问题，可以通过归一化方法获得稳定的近似解并进行校正.常见的归一化算法包括Tikhonov归一化和Marquatdt方法.Marquatdt方法是通过引入和归一化单位矩阵来实现的.此参数减少了Hessian矩阵中的广义逆条件的数量，以减少不良条件.本文主要利用Tikhonov正则化规则来修改Newton-Raphson算法，这也是目前使用的最成熟，最常用的修改方法，Tikhonov正则化的基本原理是最小二乘问题.在实际中，求解方法是在目标函数中加入一定的约束元素，在一定程度上提高了解的稳定性.把最小化函数看作： Z(a,L)=argmin||h(Z)−d|(2-30)中，Z0为参数Z的初始值，α是修正参数，F(Z)=||L(Z−Z0 Z(a,L)(k+1)=Z(a,L)通过上述公式，原来求Jaccobian矩阵Jk的逆的问题转化成了求JkTJk所以将这种修正方法使用到逆问题求解中，对f(σ) f(σ)=||F(σ)−V||其中λ为正则化参数，R是一个调制矩阵，它包含一些关于电