CN112798654A - 用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法 - Google Patents

Abstract

本发明公开了一种用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法，分析目标物体的形状数据，结合目标物体场域内的电阻抗信息，在计算机上完成标准三维有限元模型的构建。基于构建的模型，在相对电流激励模式下，获得16组场域电势值和一组参考时刻边界电压测量值，在相邻电流激励模式下，获得16组场域电势值。本发明在EIT的逆问题求解过程中，使用改进的快速梯度法求解电导率分布，降低了算法求解的复杂度，提高了迭代计算速度。该快速梯度法和自适应构建雅可比矩阵方法可以应用于其他电阻抗层析成像算法中，适用性广泛。

Description

用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构 方法

技术领域

本发明属于电阻抗层析成像领域，具体涉及一种用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法。

背景技术

电阻抗层析成像(EIT)是一种重建物体内部电导率特性的成像方法，根据不同介质具有不同电阻抗这一基本物理性质，在一定的激励和测量模式下，通过电极对测量目标施加激励信号，然后测量边界信号得出被测区域内介质的分布信息，进而对目标物体被测平面的电阻抗特性分布信息进行成像。与CT和MRI相比，电阻抗层析成像具有响应时间快、低成本、无辐射、便携和非侵入等优点，在工业过程、大地勘探和医学成像等领域中得到了广泛的应用。然而，电阻抗层析成像的逆问题天然是病态的，也就是说边界电压测量值的数目远小于需要重建的像素数目，所以逆问题的解不唯一，而且测量电压对噪声很敏感，很难得到稳定的解。电阻抗层析成像的病态性导致成像的空间分辨率通常较低，不能得到满意的效果。为提高电阻抗层析成像技术的成像精度，开展了一系列的研究，各种各样的算法被用来解决逆问题的病态性，在这些算法中正则化方法是最常用的有效方法。正则化方法的思想是选取一个先验信息(惩罚项)对最小化误差函数(保真项)上加约束来逼近真实解，先验信息的选取和最小化误差函数的各种形式产生了不同的正则化方法，比如：Tikhonov正则化算法、TV正则化算法。为了获得最优解并且提高解的稳定性，这些算法通常需要进行迭代求解，但是在迭代求解的过程中包含大量的矩阵求逆和转置的运算，极大地增加了运算量，使EIT的实时性降低。此外，电阻抗层析成像还受“软场效应”的影响，即物体的场域分布受物体场域内介质分布的影响，这会导致计算得到的雅可比矩阵与图像重建时物体的场域分布不匹配，边界测量信号存在一定的误差，从而导致重建的图像易产生伪影且重建目标不准确。

针对现有电阻抗层析成像方法运算量大造成的成像耗时较长以及雅可比矩阵与测量场域分布不匹配的问题，有必要发明一种适用于电阻抗层析成像的快速计算方法和自适应构建雅可比矩阵方法。本发明方法不但可以有效减少算法迭代过程运算量，还可通过实时重构雅可比矩阵，减少边界测量信号的建模计算误差。

发明内容

本发明解决的技术问题是提供了一种适用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法。该方法通过准确构建目标物体三维有限元模型来计算测量场域的雅可比矩阵。将计算得到的雅可比矩阵进行初始标准化处理，同时在迭代计算过程中自适应重构雅可比矩阵来减少建模误差，降低软场效应的影响。在EIT的逆问题求解过程中，使用改进的快速梯度法求解电导率分布，降低了算法求解的复杂度，提高了迭代计算速度。该快速梯度法和自适应构建雅可比矩阵方法可以应用于其他电阻抗层析成像算法中，适用性广泛。

本发明解决其技术问题所采用的技术方案：分析目标物体的形状数据，结合目标物体场域内的电阻抗信息，在计算机上完成标准三维有限元模型的构建。基于构建的模型，在相对电流激励模式下，获得16组场域电势值和一组参考时刻边界电压测量值，在相邻电流激励模式下，获得16组场域电势值。在t时刻对于真实测量物体，利用相对电流激励模式，获得一组边界电压测量值。利用获得的场域电势值计算雅可比矩阵A，并将雅可比矩阵进行标准化，标准化方法为:

式中，A_m,n是位于雅可比矩阵m行n列的元素，

是标准化之后的位于雅可比矩阵m行n列的元素。

在逆问题中建立的最小化目标函数为：

式中，||||表示欧几里得范数；||表示绝对值；g为电导率变化量，g_n∈g；N为图像重建时测量场域划分的像素数；b为边界电压测量变化值，即参考时刻边界电压测量值与t时刻边界电压测量值的差值。在目标函数中引入辅助变量p＝Ag和q＝g，其增广拉格朗日函数表示为：

式中，T表示转置，λ为正则化参数，β₁为惩罚项参数Ⅰ、β₂为惩罚项参数Ⅱ，γ为增广拉格朗日乘子Ⅰ、δ为增广拉格朗日乘子Ⅱ，γ和δ的更新方法如下：

对辅助变量p，q和电导率变化量g进行求解：

迭代求解辅助变量p，q和电导率变化量g。迭代求解过程如下：

(1)初始化各个参数及变量，设置迭代次数k，最大迭代次数为k_max；

(2)设置C^k＝L_A(g,p,q；γ,δ)，C表示函数平均值，然后进入内部迭代求解；

(3)通过计算步长与梯度直接求解得到电导率变化量g，避免大量的矩阵运算，并且加入收敛判断条件Ⅰ，确保步长不会太大，使迭代计算收敛，该计算方法还可以应用到其他电阻抗层析成像算法的求解过程中。电导率变化量g的求解过程如下：

更新步长μ：

上标k为第k次迭代；

式中，D^k＝g^k-g^k-1，y^k＝d^k(g^k)-d^k(g^k-1)，d是目标函数的梯度方向；

判断μ是否满足收敛条件Ⅰ：G_k(g^k-μ^kd^k)≤C^k-ωμ^k(d^k)^Td_k

式中，ω为权重参数。若μ不满足该条件，则令μ^k＝ρμ^k，ρ为收缩参数。

更新电导率变化量g：g^k+1＝g^k-μ^kd^k；

(4)更新辅助变量p：

(5)更新辅助变量q：

式中，max表示最大值函数；

(6)自适应重构雅可比矩阵来匹配测量场域介质的分布，从而减少建模误差、降低软场效应的影响，该重构方法也可应用于其他线性化迭代的电阻抗成像算法中。自适应雅可比矩阵重构方法如下：

式中，E为对角矩阵，e_n,n为重构因子，A^*是重构之后的雅可比矩阵。

自适应重构因子e_n,n表示如下，式中，e为重构参数，exp表示以自然常数e为底的指数函数。

(7)若

则结束内部迭代过程，否则转到(3)；

(8)更新增广拉格朗日乘子Ⅰ、Ⅱ，迭代次数k＝k+1；

(9)若k小于最大迭代次数k_max，则转到(2)，否则停止迭代。

(10)更新电导率变化量

本发明的有益效果是：本发明提出了用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法。首先将计算得到的雅可比矩阵进行标准化处理，同时在迭代计算电导率变化量的过程中自适应重构雅可比矩阵来减少误差，降低软场效应的影响。在EIT的逆问题求解过程中，使用快速梯度法，降低了算法求解的复杂度，提高了迭代计算速度。该快速梯度法和自适应构建雅可比矩阵方法可以应用于多种电阻抗层析成像求解过程中，适用性广泛，在减少重建图像伪影、提高重建图像的准确度、加快成像速度等方面均具有很好的效果。

附图说明

图1为本发明的流程框图；

图2为快速梯度法流程框图；

图3为雅可比矩阵自适应重构流程框图；

图4为在四种模型下，Tikhonov正则化方法和本发明所提方法的图像重建结果图；

图5为图4中四种模型下重建结果的模糊半径(BR)和结构相似度(SSIM)。

具体实施方式

下面结合附图对本发明作进一步说明。

结合图1的本发明流程框图对本发明进行说明：

步骤一：分析目标物体的形状数据，结合目标物体场域内的电阻抗信息，在计算机上完成标准三维有限元模型的构建。

步骤二：在三维有限元模型中，将16个电极等距围绕贴合在检测平面与模型表面相交的闭合曲线上，闭合曲线所围平面区域即为测量场域。在相对电流激励模式下，获得16组场域电势值

和一组参考时刻边界电压测量值b₀

在相邻电流激励模式下，获得16组场域电势值

步骤三：由步骤二中得到的场域电势值计算雅可比矩阵A，雅可比矩阵包含M行、N列，M表示对所有电极对依次激励时获得的测量值数量之和，N表示测量场域划分的像素数。第m行的雅可比向量计算公式为：

式中，Ω表示测量场域，上标1、2分别表示相对激励模式和相邻激励模式，

分别是在第i次相对激励和第j次相邻激励条件下N个场域单元的电势值，

和

是在相对激励和相邻激励条件下的激励电流，1≤i≤16、1≤j≤16。

步骤四：将步骤三中计算得到的雅可比矩阵进行标准化，以提高成像质量：

式中，A_m,n是位于雅可比矩阵m行n列的元素，

是标准化之后的位于雅可比矩阵m行n列的元素。

步骤五：在真实测量物体表面等距围绕贴合16个电极，贴合的位置与步骤二中电极贴合位置相同。在相对电流激励模式下，获得一组边界电压测量值b_t

步骤六：将电阻抗层析成像的图像重建过程当作一个病态的逆问题，在逆问题中建立最小化目标函数为：

式中，|| ||表示欧几里得范数；| |表示绝对值；g为电导率变化量，g_n∈g；N为图像重建时测量场域划分的单元数；b为边界电压测量变化值。

步骤七：在目标函数中引入辅助变量p＝Ag和q＝g，其增广拉格朗日函数表示为：

式中，T表示转置，λ为正则化参数，β₁为惩罚项参数Ⅰ、β₂为惩罚项参数Ⅱ，γ为增广拉格朗日乘子Ⅰ、δ为增广拉格朗日乘子Ⅱ，γ和δ的更新方法如下：

步骤八：对辅助变量p，q和电导率变化量g进行求解：

其迭代求解过程如下：

(1)初始化各个参数及变量，设置迭代次数k，最大迭代次数为k_max；

(2)设置C^k＝L_A(g,p,q；γ,δ)，C表示函数平均值，然后进入内部迭代求解；

(3)更新步长μ：

上标k为第k次迭代；

式中，D^k＝g^k-g^k-1，y^k＝d^k(g^k)-d^k(g^k-1)，d是目标函数的梯度方向；

判断μ是否满足收敛条件Ⅰ：

G_k(g^k-μ^kd^k)≤C^k-ωμ^k(d^k)^Td_k

式中，ω为权重参数。若μ不满足该条件，则令μ^k＝ρμ^k，ρ为收缩参数，再转到(4)；若μ满足该条件，则直接转到(4)；

(4)更新电导率变化量g：g^k+1＝g^k-μ^kd^k；

(5)更新辅助变量p：

(6)更新辅助变量q：

式中，max表示最大值函数；

(7)对雅可比矩阵进行自适应重构：

式中，E为对角矩阵，e_n,n为重构因子，A^*是重构之后的雅可比矩阵。

自适应重构因子e_n,n表示如下，式中，e为重构参数，exp表示以自然常数e为底的指数函数。

(8)若

结束内部迭代过程，否则转到(3)；

(9)更新增广拉格朗日乘子Ⅰ、Ⅱ，迭代次数k＝k+1；

(10)若k小于最大迭代次数k_max则转到(2)，否则停止迭代。

步骤八：更新电导率变化量

将更新得到的电导率变化量按坐标信息进行图像重建。

图4用Tikhonov算法和所提方法进行对比，比较了不同数量、不同位置的重建图像。可以看出，本方法重建的图像背景清晰、基本没有伪影，重建结果准确，而传统Tikhonov算法所重建的图像伪影较多、重建目标不准确。

同时，为了定量分析上述重建图像，采用模糊半径(Blur Radius，BR)和结构相似性(Structural Similarity Index，SSIM)对图像进行比较。图像的模糊半径越接近于0越好，结构相似性越接近于1越好。图5给出了上述重建图像的模糊半径和结构相似性。可以看出，本发明所提方法重建图像的模糊半径远远小于Tikhonov算法重建图像的模糊半径，而本发明所提方法重建图像的结构相似性均大于Tikhonov算法重建图像的结构相似性，验证了本发明方法的优越性。

以上所述仅为本发明的较佳实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法，其特征在于：在电阻抗层析成像逆问题的计算中引入自适应雅可比矩阵重构和快速梯度法，所述方法按照以下步骤进行：

步骤一：分析目标物体的形状数据，结合目标物体场域内的电阻抗信息，在计算机上完成标准三维有限元模型的构建；

步骤二：在三维有限元模型中，将16个电极等距围绕贴合在检测平面与模型表面相交的闭合曲线上，闭合曲线所围平面区域即为测量场域，在相对电流激励模式下，获得16组场域电势值φ_i ¹和一组参考时刻边界电压测量值b₀，

在相邻电流激励模式下，获得16组场域电势值

步骤三：由步骤二中得到的场域电势值计算雅可比矩阵A，雅可比矩阵包含M行、N列，M表示对所有电极对依次激励时获得的测量值数量之和，N表示测量场域划分的像素数，第m行的雅可比向量计算公式为：

式中，Ω表示测量场域，上标1、2分别表示相对激励模式和相邻激励模式，

分别是在第i次相对激励和第j次相邻激励条件下N个场域单元的电势值，

和

是在相对激励和相邻激励条件下的激励电流，1≤i≤16、1≤j≤16；

步骤四：将步骤三中计算得到的雅可比矩阵进行标准化，以提高成像质量：

式中，A_m,n是位于雅可比矩阵m行n列的元素，

是标准化之后的位于雅可比矩阵m行n列的元素；

步骤五：在真实测量物体表面等距围绕贴合16个电极，贴合的位置与步骤二中电极贴合位置相同，t时刻在相对电流激励模式下，获得一组边界电压测量值b_t，

步骤六：将电阻抗层析成像的图像重建过程当作一个病态的逆问题，在逆问题中建立的最小化目标函数为：

式中，|| ||表示欧几里得范数；| |表示绝对值；g为电导率变化量，g_n∈g；b为边界电压测量变化值，b＝b_t-b₀；

步骤七：在目标函数中引入辅助变量p＝Ag和q＝g，其增广拉格朗日函数表示为：

式中，T表示转置，λ为正则化参数，β₁为惩罚项参数Ⅰ、β₂为惩罚项参数Ⅱ，γ为增广拉格朗日乘子Ⅰ、δ为增广拉格朗日乘子Ⅱ，γ和δ的更新方法如下：

步骤八：对辅助变量p，q和电导率变化量g进行求解：

其迭代求解过程如下：

(1)初始化各个参数及变量，设置迭代次数k，最大迭代次数为k_max；

(2)设置C^k＝L_A(g,p,q；γ,δ)，C表示函数平均值，然后进入内部迭代求解；

(3)更新步长μ：

上标k为第k次迭代；

式中，D^k＝g^k-g^k-1，y^k＝d^k(g^k)-d^k(g^k-1)，d是目标函数的梯度方向；

判断μ是否满足收敛条件Ⅰ：

式中，ω为权重参数，若μ不满足该条件，则令μ^k＝ρμ^k，ρ为收缩参数，再转到(4)；若μ满足该条件，则直接转到(4)；

(4)更新电导率变化量g：g^k+1＝g^k-μ^kd^k；

(5)更新辅助变量p：

(6)更新辅助变量q：

式中，max表示最大值函数；

(7)对雅可比矩阵进行自适应重构：

式中，E为对角矩阵，e_n,n为重构因子，A^*是重构之后的雅可比矩阵；

自适应重构因子e_n,n表示如下，式中，e为重构参数，exp表示以自然常数e为底的指数函数；

(8)若

结束内部迭代过程，否则转到(3)；

(9)更新增广拉格朗日乘子Ⅰ、Ⅱ，迭代次数k＝k+1；

(10)若k小于最大迭代次数k_max则转到(2)，否则停止迭代；

步骤九：更新电导率变化量

将更新得到的电导率变化量按坐标信息进行图像重建。

Images