CN112798654B - 用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法 - Google Patents
用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法 Download PDFInfo
- Publication number
- CN112798654B CN112798654B CN202110108632.XA CN202110108632A CN112798654B CN 112798654 B CN112798654 B CN 112798654B CN 202110108632 A CN202110108632 A CN 202110108632A CN 112798654 B CN112798654 B CN 112798654B
- Authority
- CN
- China
- Prior art keywords
- jacobian matrix
- electrical impedance
- reconstruction
- adaptive
- field
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N27/00—Investigating or analysing materials by the use of electric, electrochemical, or magnetic means
- G01N27/02—Investigating or analysing materials by the use of electric, electrochemical, or magnetic means by investigating impedance
- G01N27/04—Investigating or analysing materials by the use of electric, electrochemical, or magnetic means by investigating impedance by investigating resistance
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/15—Correlation function computation including computation of convolution operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Chemical & Material Sciences (AREA)
- General Engineering & Computer Science (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Algebra (AREA)
- Computing Systems (AREA)
- Biochemistry (AREA)
- Life Sciences & Earth Sciences (AREA)
- Electrochemistry (AREA)
- General Health & Medical Sciences (AREA)
- Chemical Kinetics & Catalysis (AREA)
- Health & Medical Sciences (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Analytical Chemistry (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
本发明公开了一种用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法,分析目标物体的形状数据,结合目标物体场域内的电阻抗信息,在计算机上完成标准三维有限元模型的构建。基于构建的模型,在相对电流激励模式下,获得16组场域电势值和一组参考时刻边界电压测量值,在相邻电流激励模式下,获得16组场域电势值。本发明在EIT的逆问题求解过程中,使用改进的快速梯度法求解电导率分布,降低了算法求解的复杂度,提高了迭代计算速度。该快速梯度法和自适应构建雅可比矩阵方法可以应用于其他电阻抗层析成像算法中,适用性广泛。
Description
技术领域
本发明属于电阻抗层析成像领域,具体涉及一种用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法。
背景技术
电阻抗层析成像(EIT)是一种重建物体内部电导率特性的成像方法,根据不同介质具有不同电阻抗这一基本物理性质,在一定的激励和测量模式下,通过电极对测量目标施加激励信号,然后测量边界信号得出被测区域内介质的分布信息,进而对目标物体被测平面的电阻抗特性分布信息进行成像。与CT和MRI相比,电阻抗层析成像具有响应时间快、低成本、无辐射、便携和非侵入等优点,在工业过程、大地勘探和医学成像等领域中得到了广泛的应用。然而,电阻抗层析成像的逆问题天然是病态的,也就是说边界电压测量值的数目远小于需要重建的像素数目,所以逆问题的解不唯一,而且测量电压对噪声很敏感,很难得到稳定的解。电阻抗层析成像的病态性导致成像的空间分辨率通常较低,不能得到满意的效果。为提高电阻抗层析成像技术的成像精度,开展了一系列的研究,各种各样的算法被用来解决逆问题的病态性,在这些算法中正则化方法是最常用的有效方法。正则化方法的思想是选取一个先验信息(惩罚项)对最小化误差函数(保真项)上加约束来逼近真实解,先验信息的选取和最小化误差函数的各种形式产生了不同的正则化方法,比如:Tikhonov正则化算法、TV正则化算法。为了获得最优解并且提高解的稳定性,这些算法通常需要进行迭代求解,但是在迭代求解的过程中包含大量的矩阵求逆和转置的运算,极大地增加了运算量,使EIT的实时性降低。此外,电阻抗层析成像还受“软场效应”的影响,即物体的场域分布受物体场域内介质分布的影响,这会导致计算得到的雅可比矩阵与图像重建时物体的场域分布不匹配,边界测量信号存在一定的误差,从而导致重建的图像易产生伪影且重建目标不准确。
针对现有电阻抗层析成像方法运算量大造成的成像耗时较长以及雅可比矩阵与测量场域分布不匹配的问题,有必要发明一种适用于电阻抗层析成像的快速计算方法和自适应构建雅可比矩阵方法。本发明方法不但可以有效减少算法迭代过程运算量,还可通过实时重构雅可比矩阵,减少边界测量信号的建模计算误差。
发明内容
本发明解决的技术问题是提供了一种适用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法。该方法通过准确构建目标物体三维有限元模型来计算测量场域的雅可比矩阵。将计算得到的雅可比矩阵进行初始标准化处理,同时在迭代计算过程中自适应重构雅可比矩阵来减少建模误差,降低软场效应的影响。在EIT的逆问题求解过程中,使用改进的快速梯度法求解电导率分布,降低了算法求解的复杂度,提高了迭代计算速度。该快速梯度法和自适应构建雅可比矩阵方法可以应用于其他电阻抗层析成像算法中,适用性广泛。
本发明解决其技术问题所采用的技术方案:分析目标物体的形状数据,结合目标物体场域内的电阻抗信息,在计算机上完成标准三维有限元模型的构建。基于构建的模型,在相对电流激励模式下,获得16组场域电势值和一组参考时刻边界电压测量值,在相邻电流激励模式下,获得16组场域电势值。在t时刻对于真实测量物体,利用相对电流激励模式,获得一组边界电压测量值。利用获得的场域电势值计算雅可比矩阵A,并将雅可比矩阵进行标准化,标准化方法为:
在逆问题中建立的最小化目标函数为:
式中,||||表示欧几里得范数;||表示绝对值;g为电导率变化量,gn∈g;N为图像重建时测量场域划分的像素数;b为边界电压测量变化值,即参考时刻边界电压测量值与t时刻边界电压测量值的差值。在目标函数中引入辅助变量p=Ag和q=g,其增广拉格朗日函数表示为:
式中,T表示转置,λ为正则化参数,β1为惩罚项参数Ⅰ、β2为惩罚项参数Ⅱ,γ为增广拉格朗日乘子Ⅰ、δ为增广拉格朗日乘子Ⅱ,γ和δ的更新方法如下:
对辅助变量p,q和电导率变化量g进行求解:
迭代求解辅助变量p,q和电导率变化量g。迭代求解过程如下:
(1)初始化各个参数及变量,设置迭代次数k,最大迭代次数为kmax;
(2)设置Ck=LA(g,p,q;γ,δ),C表示函数平均值,然后进入内部迭代求解;
(3)通过计算步长与梯度直接求解得到电导率变化量g,避免大量的矩阵运算,并且加入收敛判断条件Ⅰ,确保步长不会太大,使迭代计算收敛,该计算方法还可以应用到其他电阻抗层析成像算法的求解过程中。电导率变化量g的求解过程如下:
式中,Dk=gk-gk-1,yk=dk(gk)-dk(gk-1),d是目标函数的梯度方向;
判断μ是否满足收敛条件Ⅰ:Gk(gk-μkdk)≤Ck-ωμk(dk)Tdk
式中,ω为权重参数。若μ不满足该条件,则令μk=ρμk,ρ为收缩参数。
更新电导率变化量g:gk+1=gk-μkdk;
式中,max表示最大值函数;
(6)自适应重构雅可比矩阵来匹配测量场域介质的分布,从而减少建模误差、降低软场效应的影响,该重构方法也可应用于其他线性化迭代的电阻抗成像算法中。自适应雅可比矩阵重构方法如下:
式中,E为对角矩阵,en,n为重构因子,A*是重构之后的雅可比矩阵。
自适应重构因子en,n表示如下,式中,e为重构参数,exp表示以自然常数e为底的指数函数。
(8)更新增广拉格朗日乘子Ⅰ、Ⅱ,迭代次数k=k+1;
(9)若k小于最大迭代次数kmax,则转到(2),否则停止迭代。
本发明的有益效果是:本发明提出了用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法。首先将计算得到的雅可比矩阵进行标准化处理,同时在迭代计算电导率变化量的过程中自适应重构雅可比矩阵来减少误差,降低软场效应的影响。在EIT的逆问题求解过程中,使用快速梯度法,降低了算法求解的复杂度,提高了迭代计算速度。该快速梯度法和自适应构建雅可比矩阵方法可以应用于多种电阻抗层析成像求解过程中,适用性广泛,在减少重建图像伪影、提高重建图像的准确度、加快成像速度等方面均具有很好的效果。
附图说明
图1为本发明的流程框图;
图2为快速梯度法流程框图;
图3为雅可比矩阵自适应重构流程框图;
图4为在四种模型下,Tikhonov正则化方法和本发明所提方法的图像重建结果图;
图5为图4中四种模型下重建结果的模糊半径(BR)和结构相似度(SSIM)。
具体实施方式
下面结合附图对本发明作进一步说明。
结合图1的本发明流程框图对本发明进行说明:
步骤一:分析目标物体的形状数据,结合目标物体场域内的电阻抗信息,在计算机上完成标准三维有限元模型的构建。
步骤二:在三维有限元模型中,将16个电极等距围绕贴合在检测平面与模型表面相交的闭合曲线上,闭合曲线所围平面区域即为测量场域。在相对电流激励模式下,获得16组场域电势值和一组参考时刻边界电压测量值b0/>在相邻电流激励模式下,获得16组场域电势值/>
步骤三:由步骤二中得到的场域电势值计算雅可比矩阵A,雅可比矩阵包含M行、N列,M表示对所有电极对依次激励时获得的测量值数量之和,N表示测量场域划分的像素数。第m行的雅可比向量计算公式为:
式中,Ω表示测量场域,上标1、2分别表示相对激励模式和相邻激励模式,分别是在第i次相对激励和第j次相邻激励条件下N个场域单元的电势值,/>和/>是在相对激励和相邻激励条件下的激励电流,1≤i≤16、1≤j≤16。
步骤四:将步骤三中计算得到的雅可比矩阵进行标准化,以提高成像质量:
步骤六:将电阻抗层析成像的图像重建过程当作一个病态的逆问题,在逆问题中建立最小化目标函数为:
式中,|| ||表示欧几里得范数;| |表示绝对值;g为电导率变化量,gn∈g;N为图像重建时测量场域划分的单元数;b为边界电压测量变化值。
步骤七:在目标函数中引入辅助变量p=Ag和q=g,其增广拉格朗日函数表示为:
式中,T表示转置,λ为正则化参数,β1为惩罚项参数Ⅰ、β2为惩罚项参数Ⅱ,γ为增广拉格朗日乘子Ⅰ、δ为增广拉格朗日乘子Ⅱ,γ和δ的更新方法如下:
步骤八:对辅助变量p,q和电导率变化量g进行求解:
其迭代求解过程如下:
(1)初始化各个参数及变量,设置迭代次数k,最大迭代次数为kmax;
(2)设置Ck=LA(g,p,q;γ,δ),C表示函数平均值,然后进入内部迭代求解;
式中,Dk=gk-gk-1,yk=dk(gk)-dk(gk-1),d是目标函数的梯度方向;
判断μ是否满足收敛条件Ⅰ:
Gk(gk-μkdk)≤Ck-ωμk(dk)Tdk
式中,ω为权重参数。若μ不满足该条件,则令μk=ρμk,ρ为收缩参数,再转到(4);若μ满足该条件,则直接转到(4);
(4)更新电导率变化量g:gk+1=gk-μkdk;
式中,max表示最大值函数;
(7)对雅可比矩阵进行自适应重构:
式中,E为对角矩阵,en,n为重构因子,A*是重构之后的雅可比矩阵。
自适应重构因子en,n表示如下,式中,e为重构参数,exp表示以自然常数e为底的指数函数。
(9)更新增广拉格朗日乘子Ⅰ、Ⅱ,迭代次数k=k+1;
(10)若k小于最大迭代次数kmax则转到(2),否则停止迭代。
图4用Tikhonov算法和所提方法进行对比,比较了不同数量、不同位置的重建图像。可以看出,本方法重建的图像背景清晰、基本没有伪影,重建结果准确,而传统Tikhonov算法所重建的图像伪影较多、重建目标不准确。
同时,为了定量分析上述重建图像,采用模糊半径(Blur Radius,BR)和结构相似性(Structural Similarity Index,SSIM)对图像进行比较。图像的模糊半径越接近于0越好,结构相似性越接近于1越好。图5给出了上述重建图像的模糊半径和结构相似性。可以看出,本发明所提方法重建图像的模糊半径远远小于Tikhonov算法重建图像的模糊半径,而本发明所提方法重建图像的结构相似性均大于Tikhonov算法重建图像的结构相似性,验证了本发明方法的优越性。
以上所述仅为本发明的较佳实施例而已,并不用于限制本发明,凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (1)
1.用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法,其特征在于:在电阻抗层析成像逆问题的计算中引入自适应雅可比矩阵重构和快速梯度法,所述方法按照以下步骤进行:
步骤一:分析目标物体的形状数据,结合目标物体场域内的电阻抗信息,在计算机上完成标准三维有限元模型的构建;
步骤二:在三维有限元模型中,将16个电极等距围绕贴合在检测平面与模型表面相交的闭合曲线上,闭合曲线所围平面区域即为测量场域,在相对电流激励模式下,获得16组场域电势值φi 1和一组参考时刻边界电压测量值b0,在相邻电流激励模式下,获得16组场域电势值/>
步骤三:由步骤二中得到的场域电势值计算雅可比矩阵A,雅可比矩阵包含M行、N列,M表示对所有电极对依次激励时获得的测量值数量之和,N表示测量场域划分的像素数,第m行的雅可比向量计算公式为:
式中,Ω表示测量场域,上标1、2分别表示相对激励模式和相邻激励模式,分别是在第i次相对激励和第j次相邻激励条件下N个场域单元的电势值,/>和/>是在相对激励和相邻激励条件下的激励电流,1≤i≤16、1≤j≤16;
步骤四:将步骤三中计算得到的雅可比矩阵进行标准化,以提高成像质量:
步骤六:将电阻抗层析成像的图像重建过程当作一个病态的逆问题,在逆问题中建立的最小化目标函数为:
式中,|| ||表示欧几里得范数;| |表示绝对值;g为电导率变化量,gn∈g;b为边界电压测量变化值,b=bt-b0;
步骤七:在目标函数中引入辅助变量p=Ag和q=g,其增广拉格朗日函数表示为:
式中,T表示转置,λ为正则化参数,β1为惩罚项参数Ⅰ、β2为惩罚项参数Ⅱ,γ为增广拉格朗日乘子Ⅰ、δ为增广拉格朗日乘子Ⅱ,γ和δ的更新方法如下:
步骤八:对辅助变量p,q和电导率变化量g进行求解:
其迭代求解过程如下:
(1)初始化各个参数及变量,设置迭代次数k,最大迭代次数为kmax;
(2)设置Ck=LA(g,p,q;γ,δ),C表示函数平均值,然后进入内部迭代求解;
式中,Dk=gk-gk-1,yk=dk(gk)-dk(gk-1),d是目标函数的梯度方向;
判断μ是否满足收敛条件Ⅰ:
式中,ω为权重参数,若μ不满足该条件,则令μk=ρμk,ρ为收缩参数,再转到(4);若μ满足该条件,则直接转到(4);
(4)更新电导率变化量g:gk+1=gk-μkdk;
式中,max表示最大值函数;
(7)对雅可比矩阵进行自适应重构:
式中,E为对角矩阵,en,n为重构因子,A*是重构之后的雅可比矩阵;
自适应重构因子en,n表示如下,式中,e为重构参数,exp表示以自然常数e为底的指数函数;
(9)更新增广拉格朗日乘子Ⅰ、Ⅱ,迭代次数k=k+1;
(10)若k小于最大迭代次数kmax则转到(2),否则停止迭代;
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110108632.XA CN112798654B (zh) | 2021-01-26 | 2021-01-26 | 用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110108632.XA CN112798654B (zh) | 2021-01-26 | 2021-01-26 | 用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112798654A CN112798654A (zh) | 2021-05-14 |
CN112798654B true CN112798654B (zh) | 2023-06-20 |
Family
ID=75811987
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110108632.XA Active CN112798654B (zh) | 2021-01-26 | 2021-01-26 | 用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112798654B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114224313B (zh) * | 2021-12-27 | 2023-01-17 | 深圳融昕医疗科技有限公司 | 电阻抗成像方法及计算机可读存储介质 |
CN116524123B (zh) * | 2023-04-20 | 2024-02-13 | 深圳市元甪科技有限公司 | 一种三维电阻抗断层扫描图像重建方法及相关设备 |
CN116824048B (zh) * | 2023-06-05 | 2024-01-30 | 南京航空航天大学 | 一种传感器、雅可比矩阵求解方法、三维成像系统及方法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2012061475A2 (en) * | 2010-11-02 | 2012-05-10 | University Of Florida Research Foundation, Inc. | Systems and methods for fast magnetic resonance image reconstruction |
CN111616708A (zh) * | 2020-05-25 | 2020-09-04 | 中国人民解放军第四军医大学 | 一种精准识别脑卒中颅内病变区域的图像重建方法 |
-
2021
- 2021-01-26 CN CN202110108632.XA patent/CN112798654B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2012061475A2 (en) * | 2010-11-02 | 2012-05-10 | University Of Florida Research Foundation, Inc. | Systems and methods for fast magnetic resonance image reconstruction |
CN111616708A (zh) * | 2020-05-25 | 2020-09-04 | 中国人民解放军第四军医大学 | 一种精准识别脑卒中颅内病变区域的图像重建方法 |
Non-Patent Citations (1)
Title |
---|
动态电阻抗成像时空相关性重建方法研究;王琦;彭圆圆;汪剑鸣;连志杰;李秀艳;陈庆良;王化祥;窦汝振;;电子测量与仪器学报(第02期);全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN112798654A (zh) | 2021-05-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN112798654B (zh) | 用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法 | |
CN109919844B (zh) | 一种高分辨率的电学层析成像电导率分布重建方法 | |
CN109584323B (zh) | 超声反射信息约束的腹部病变电阻抗图像重建方法 | |
CN107146228B (zh) | 一种基于先验知识的大脑磁共振图像超体素生成方法 | |
CN109859285B (zh) | 基于空洞卷积网络的电阻抗图像重建方法 | |
CN109035352B (zh) | L1-l2空间自适应电学层析成像正则化重建方法 | |
CN109934885B (zh) | 一种锐利边缘保持的电阻层析成像图像重建方法 | |
Mair et al. | Estimation of images and nonrigid deformations in gated emission CT | |
CN108830875B (zh) | 一种基于残差最小的电阻抗层析成像图像分割方法 | |
CN110934586B (zh) | 一种灰度值矩阵快速分解重建的正则化方法 | |
CN111047663A (zh) | 电学层析成像伪影抑制图像重建方法 | |
WO2015112804A1 (en) | System and method for generating magnetic resonance imaging (mri) images using structures of the images | |
CN111161182B (zh) | Mr结构信息约束的非局部均值引导的pet图像部分容积校正方法 | |
CN110910466B (zh) | 一种新型多频差分电阻抗层析成像重建算法 | |
CN114549682A (zh) | 一种用于电阻抗肺部成像图像的优化方法 | |
CN112581385B (zh) | 基于多先验约束扩散峰度成像张量估计方法、介质和设备 | |
CN116843679B (zh) | 基于深度图像先验框架的pet图像部分容积校正方法 | |
CN114052701A (zh) | 一种电容耦合电阻层析成像图像重建方法 | |
CN117197159A (zh) | 一种高效且高质量的高温物体点云分割方法 | |
Yan et al. | A multi-modal medical image fusion method in spatial domain | |
CN110992385B (zh) | 抑制伪影与保护边缘的颅内图像重建方法 | |
CN113034632A (zh) | 一种用于检测工业两相流的图像重建方法 | |
CN113870377A (zh) | 基于V-ResNet肺部成像方法 | |
CN113034633B (zh) | 一种目标函数准确限定支配的图像重建方法 | |
CN111476888B (zh) | 一种基于三维空间体拟合的医学图像层间插值方法,设备及可读存储介质 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |