CN108711178A - 一种基于闭环控制原理的电容层析成像图像重建方法 - Google Patents

Abstract

本发明涉及一种基于闭环控制原理的电容层析成像图像重建方法，主要包括以下步骤：首先，利用电容层析成像系统通过测量分别得到空场条件下的N(N‑1)/2个独立测量值和有物体条件下的N(N‑1)/2个独立测量值，利用两者的差建立N×N的电容变化量矩阵；然后，将Calderon算法作为被控对象，利用该算法反演出的介电常数分布求解新的电容变化量矩阵作为负反馈，通过调整PID控制器的参数来修正电容变化量的偏差量，使整个闭环收敛；最后，若迭代次数达到预设迭代次数，则结束迭代并输出重建结果，否则继续进行迭代。该方法首次将闭环控制原理与Calderon算法相结合，能有效减小图像重建误差，显著提高图像重建质量，在电学成像领域具有重要的实用价值和很好的应用前景。

Description

一种基于闭环控制原理的电容层析成像图像重建方法

技术领域

本发明涉及电学成像领域，尤其涉及一种基于闭环控制原理的电容层析成像图像重建方法。

背景技术

电容层析成像(Electrical Capacitance Tomography，简称ECT)技术是20世纪80年代发展起来的一种基于电磁场敏感机理的层析成像技术，具有无辐射、非侵入性、便携性、响应速度快、价格低廉等优势，在工业、化学和医学等领域具有重要的应用前景。电容层析成像系统采用特殊设计的敏感空间阵列电极，以非侵入的方式获取被测敏感场信息，根据边界电容信息，利用特定的图像重建算法重建被测对象内部介电常数的分布，从而得到物体内部的分布情况(王化祥.电学层析成像技术[J].自动化仪表,2017,38(5):1-6.)。ECT系统主要包括三部分：空间敏感电极阵列、多通道电容采集系统以及图像重建。

图像重建是涉及到非线性偏微分方程的逆问题求解过程，是ECT的关键核心技术。对于ECT逆问题的求解，难点主要表现在：(1)欠定性。通过测量得到场域的边界信息，得到的信息量远远小于所要求解的未知量，所以ECT逆问题的解不唯一。(2)“软场”特性。场域内的电位是介质分布的函数，通过电位的测量值求解场域内部介质分布是一个非线性问题。(3)病态性。求解ECT逆问题涉及到严重的病态性问题，这导致在求解ECT逆问题时，测量边界电位值的微小变化会导致求解出的场域内部介质分布有很大的变化，实际测量中存在的噪声会对图像重建产生很大的影响。(王化祥等.电学层析成像[M].科学出版社,2013.)。

在文献中，大多数的图像重建算法采用了灵敏度矩阵原理，该原理是Geselowitz在1971年于《IEEE生物医学工程汇刊》(IEEE Transactions on Biomedical Engineering)第18卷第1期38-41页发表的题为《心电图的导联理论在阻抗测量成像中的应用》(Anapplication of electrocardiographic lead theory to impedance plethysmography)和Lehr在1972年于《IEEE生物医学工程汇刊》(IEEE Transactions on BiomedicalEngineering)第19卷第2期156-157页发表的题为《一种用于阻抗测量成像的向量推导方法》(A vector derivation useful in impedance plethysmographic fieldcalculations)的文章中提出的，是一种采用摄动原理的线性化重建方法。常见的图像重建算法有线性反投影算法(linear back projection,简称LBP)，兰德韦伯算法(Landweber)，牛顿-拉弗森法(Newton-Raphson)等。LBP算法是由Xie等人在1992年于《IEE会议论文集：电路设备与系统》(IEE Proceedings G:Circuits Devices and Systems)第139卷第1期89-98页发表的题为《用于流体成像的电容层析成像：用于图像重建算法和传感器设计的系统模型》(Electrical capacitance tomography for flow imaging:system model fordevelopment of image reconstruction algorithms and design of primary sensors)的文章中提出的，这是一种相对简单的成像方法，它将通过某点的所有投影线进行累加，再反向估算出该点的密度值，成像精度不高，多用于定性分析。Landweber算法由Yang等人在1999年于《测量科学与技术》(Measurement Science and Technology)第10卷第11期第1065页发表的题为《一种基于Landweber迭代方法的电容层析成像图像重建算法》(Animage-reconstruction algorithm based on Landweber's iteration method forelectrical-capacitance tomography)文章中提出的，用迭代的方式逼近灵敏度矩阵，由于Landweber算法沿负梯度方向进行搜索，计算收敛速度较慢。牛顿-拉弗森(Newton-Raphson)法的基本思想是将目标函数作泰勒展开，利用二次函数近似代替目标函数，求出该二次函数的极小点，一般情况下收敛性好，但运算量相对较大，且函数一阶导数有时不易求出，Edic等人于1998年在《IEEE生物医学工程汇刊》(IEEE Transactions on BiomedicalEngineering)第45卷第7期899页发表的题为《采用迭代牛顿-拉弗森方法求解导纳逆问题》(An iterative Newton-Raphson method to solve the inverse admittivity problem)的文章中对该方法进行了详细的讨论。这些方法的基本思想是假设物质分布变化较小时，场域内敏感场的分布近似不变，即忽略了“软场”特性。实际上，灵敏度矩阵需要根据真实的分布进行更新，然而，真实的分布是未知的，因此，Mirkowski等人于2008年在《IEEE仪器和测量汇刊》(IEEE Transactions on Instrumentation and Measurement)第57卷第5期973-980页发表的题为《一种基于电容-网格模型的电容层析成像正问题求解方法》(A newforward-problem solver based on a capacitor-mesh model for electricalcapacitance tomography)的文章中，提出利用重建出的图像更新灵敏度矩阵；Soleimani等人于2009年在《IEEE仪器和测量汇刊》(IEEE Transactions on Instrumentation andMeasurement)第59卷第1期78-83页发表的题为《利用电容层析成像实验数据通过亥姆霍兹正则化方法对介电常数分布进行重建》(Helmholtz-type regularization method forpermittivity reconstruction using experimental phantom data of electricalcapacitance tomography)的文章中，通过求解有限元方程更新灵敏度矩阵。虽然基于灵敏度矩阵的图像重建方法已经成功应用，但它们也有自己的局限性，这些方法将整个区域所有点的灰度值计算出来，若整个区域中只有一小部分需要重建，则会增加许多不必要的计算量。

近年来，学者们提出了一些直接图像重建算法，如Calderon算法，分解方法，D-bar方法等，直接图像重建算法可以重建出场域内部任一点的灰度值。分解方法只适用于重建具有连续背景的扰动区域，若待重建区域是环形的，则分解方法失效，分解方法多用于定性分析，Harrach等人于2010年在《IEEE医学成像汇刊》(IEEE Transactions on MedicalImaging)第29卷第11期1918-1926页发表的题为《频率差分电阻抗层析成像中的分解方法及其物理验证》(Factorization method and its physical justification infrequency-difference electrical impedance tomography)的文章中，介绍了利用分解方法进行电阻抗层析成像的基本原理及应用。D-bar方法是将非线性傅里叶变换应用到电阻抗成像中，Murphy等人于2009年在《IEEE医学成像汇刊》(IEEE Transactions onMedical Imaging)第28卷第10期1576-1584页发表的题为《场域形状建模和测量误差对EIT二维D-bar方法的影响》(Effect of domain shape modeling and measurement errorson the 2-D D-bar method for EIT)的文章中对该方法进行了讨论，利用D-bar方法进行重建计算量较大。Calderon算法是由Calderon于1980年在《计算应用和数学》(Computational and Applied Mathematics)发表的题为《关于逆边界问题》(On aninverse boundary value problem)的文章中提出的，提出了一种用于求解低电导率对比度二维逆电导率问题的线性化准则。Bikowski等人于2008年在《Inverse Problems andImaging》(逆问题与成像)第2卷第1期43-61页发表的题为《利用Calderon方法进行二维EIT重建》(2DEIT reconstructions using Calderon's method)的文章中将Calderon方法应用到二维电阻抗层析成像来进行电导率分布的重建。Cao等人于2011年在《IEEE仪器与测量汇刊》(IEEE Transactions on Instrumentation and Measurement)第60卷第3期900-907页发表的题为《利用Calderon算法针对方形截面传感器的电容层析成像》(Electricalcapacitance tomography for sensors of square cross sections using Calderon'sMethod)的文章中，提出了狄利克雷-诺依曼边界映射的构造方法，将Calderon方法应用到电容层析成像，实现场域内部介电常数的重建。

如果在电极数足够多的情况下精确得到边界电容信息，并且介电常数的扰动比较微小，Calderon算法有非常良好的重建效果，然而，在实际测量中会存在电极数量、电容测量电路分辨率和传感器精度等多方面的限制。为了满足实际应用中重建图像更高空间分辨率的要求，本发明提出一种基于闭环控制原理的电容层析成像图像重建方法，将Calderon算法与闭环控制理论相结合，对比利用Calderon算法的图像重建结果，应用该方法得到的重建结果产生的伪影少，视觉效果好，重建精度更高，有效减小了图像重建误差，也可用于其他成像技术的逆问题求解过程。

发明内容

本发明的目的在于提出一种基于闭环控制原理的电容层析成像图像重建方法，所述方法重建精度高，能有效改善电容层析成像图像重建质量。

本发明的技术方案是：

步骤一、对同一圆形截面上有等面积均匀分布的N个电极的电容层析成像传感器，将N个电极逆时针标号为j(1≤j≤N)，采用传统的相邻激励测量模式，即，一次测量过程中，第1步，电极1上施加幅值为V的交流电压，其他N-1个电极均接地或与地保持同电位，测量得到电极1与电极2到电极N共N-1个电容值；第2步，电极2上施加幅值为V的交流电压，其他N-1个电极均接地或与地保持同电位，测量得到电极2与电极3到电极N共N-2个电容值；第3步，电极3上施加幅值为V的交流电压，其他N-1个电极均接地或与地保持同电位，测量得到电极3与电极4到电极N共N-3个电容值；以此类推，第N-1步，电极N-1上施加幅值为V的交流电压，其他N-1个电极均接地或与地保持同电位，测量得到电极N-1与电极N共1个电容值；通过测量分别得到空场条件下的N(N-1)/2个独立测量值和有物场条件下的N(N-1)/2个独立测量值，所述空场条件表示场域内部只有介电常数均匀分布的背景介质，所述有物场条件表示在均匀分布的背景介质上存在与背景介质介电常数不同的物体，空场条件下与有物场条件下的测量值可分别建立N×N的电容矩阵C₀和C_m，C₀表示空场电容矩阵，C_m表示有物场电容矩阵，C_m与C₀的差表示为：

其中，△C_s,t为第s个与第t个电极之间的电容的变化量(s≠t)，△C_s,t＝△C_t,s，△C_s,s为第s个电极上的自电容变化量，所有电极的电容变化量的关系为：

利用Calderon算法计算介电常数分布。即利用公式(3)计算散射变换t(re^iθ)：

其中，(r,θ)表示极坐标参数，r为极径，θ为极角，π表示圆周率，e表示自然常数，θ_s表示第s个电极的中心位置，θ_s＝2πs/N，A表示电极面积。

利用公式(4)可反演出场域内部物质介电常数分布情况：

δε(x,y)表示坐标为(x,y)的点处介电常数的变化值，R为数值积分区域的半径。

步骤二、根据所求得的分布得到新的电容矩阵C'_m，将该电容矩阵与步骤一中测量得到的空场电容矩阵C₀相减得到△C'_N×N(k)作为负反馈信息：

其中，△C'_N×N(k)表示第k次迭代计算的负反馈电容变化量矩阵。

步骤三、计算第k次偏差p_N×N(k)＝△C_N×N-△C'_N×N(k)，利用公式(6)得到PID控制器第k次输出信息q_N×N(k)，接着利用公式(3)和(4)所描述的Calderon算法更新重建的分布。选取的数字式PID控制器为：

q(k)＝q(k-1)+(K_p+K_i+K_d)p(k)-(K_p+2K_d)p(k-1)+K_dp(k-2) (6)

其中，K_p,K_i,K_d分别为比例、积分、微分系数；q(k),q(k-1)分别为PID控制器第k次迭代计算和第(k-1)次迭代计算的输出，p(k),p(k-1),p(k-2)分别为PID控制器第k次迭代计算、第(k-1)次迭代计算和第(k-2)次迭代计算的输入，k＝1,2,3,…,m，m为预设迭代次数；设定初值为q_N×N(0)＝O_N×N,p_N×N(0)＝O_N×N,p_N×N(-1)＝O_N×N，△C'_N×N(1)＝O_N×N,O_N×N表示N×N的全零矩阵。

步骤四、判断迭代计算次数k是否大于预设迭代次数m，如果是，则执行步骤五，如果否，则令k＝k+1并返回步骤二。

步骤五、结束迭代，输出成像结果。

本发明的有益效果是：以Calderon算法为基础，将闭环控制原理引入电容层析成像图像重建算法中，以反演误差为控制对象，以图像重建算法为被控对象，显著提高了图像重建的质量。

附图说明

图1是实施流程图。

图2是16电极电容层析成像传感器。

图3是具体实施原始模型。

图4是利用Calderon算法重建效果图。

图5是将闭环控制原理与Calderon算法结合重建效果图。

具体实施方式

参见附图1，一种基于闭环控制原理的电容层析成像图像重建方法实施流程图。结合如附图2所示的16电极电容层析成像传感器及附图3所示的原始模型对本算法进行说明。

步骤一、对如附图2所示的同一圆形截面上等面积均匀分布有16个电极的电容层析成像传感器，对如附图3所示的原始分布进行测量，其中白色部分相对介电常数为3，黑色部分相对介电常数为1，将16个电极逆时针标号为j(1≤j≤16)，采用传统的相邻激励测量模式，即，一次测量过程中，第1步，电极1上施加幅值为V的交流电压，其他15个电极均接地或与地保持同电位，测量得到电极1与电极2到电极16共15个电容值；第2步，电极2上施加幅值为V的交流电压，其他15个电极均接地或与地保持同电位，测量得到电极2与电极3到电极16共14个电容值；第3步，电极3上施加幅值为V的交流电压，其他15个电极均接地或与地保持同电位，测量得到电极3与电极4到电极16共13个电容值；以此类推，第15步，电极15上施加幅值为V的交流电压，其他15个电极均接地或与地保持同电位，测量得到电极15与电极16共1个电容值；通过测量分别得到空场条件下的120个独立测量值和有物场条件下的120个独立测量值，空场表示场域内部的介质为空气，相对介电常数为1，有物场表示在均匀分布的背景介质上存在物体，相对介电常数为3，空场条件下与有物场条件下的测量值可分别建立16×16的电容矩阵C₀和C_m，C₀表示空场电容矩阵，C_m表示有物体时的电容矩阵，C_m与C₀的差表示为：

其中，△C_s,t为第s个与第t个电极之间的电容的变化量(s≠t)，△C_s,t＝△C_t,s，△C_s,s为第s个电极上的自电容变化量，所有电极的电容变化量的关系为：

利用Calderon算法计算介电常数分布。即利用公式(9)计算散射变换t(re^iθ)：

其中，(r,θ)表示极坐标参数，r为极径，θ为极角，π表示圆周率，e表示自然常数，θ_s表示第s个电极的中心位置，θ_s＝2πs/16，A表示电极面积。

利用公式(10)可反演出场域内部物质介电常数分布情况：

δε(x,y)表示坐标为(x,y)的点处介电常数的变化值，R为数值积分区域的半径。

步骤二、根据所求得的分布得到新的电容矩阵C'_m，将该电容矩阵与步骤一中测量得到的空场电容矩阵C₀相减作为负反馈信息：

其中，△C'_16×16(k)表示第k次迭代计算的负反馈电容变化量矩阵。

步骤三、计算第k次偏差p_16×16(k)＝△C_16×16-△C'_16×16(k)，该偏差作为PID控制器的输入从而得到PID控制器第k次输出信息q_16×16(k)，接着利用公式(9)和(10)所描述的Calderon算法更新重建的分布。选取的数字式PID控制器为：

q(k)＝q(k-1)+(K_p+K_i+K_d)p(k)-(K_p+2K_d)p(k-1)+K_dp(k-2) (12)

其中，K_p,K_i,K_d分别为比例、积分、微分系数；q(k),q(k-1)分别为PID控制器第k次迭代计算和第(k-1)次迭代计算的输出，p(k),p(k-1),p(k-2)分别为PID控制器第k次迭代计算、第(k-1)次迭代计算和第(k-2)次迭代计算的输入，k＝1,2,3,...,m，m为预设迭代次数；设定初值为q_16×16(0)＝O_16×16,p_16×16(0)＝O_16×16，p_16×16(-1)＝O_16×16，△C'_16×16(1)＝O_16×16,O_16×16表示16×16的全零矩阵。PID控制器参数通过遍历寻优设置为K_p＝1,K_i＝100,K_d＝0。

步骤四、预设迭代次数设置为m＝100，判断迭代计算次数k是否大于预设迭代次数m，如果是，则执行步骤五，如果否，则令k＝k+1并返回步骤二。

步骤五、结束迭代，输出成像结果。利用Calderon算法反演出场域内部物质介电常数分布情况如附图4所示，将闭环控制原理与Calderon算法相结合得到的反演结果如附图5所示。

通过比较附图3、附图4和附图5，将闭环控制原理和Calderon算法相结合可有效提高图像重建质量，重建出的图像更加接近原始分布，边缘更加清晰。

以上对本发明及其实施方式的描述，并不局限于此，附图中所示仅是本发明的实施方式之一。在不脱离本发明创造宗旨的情况下，不经创造地设计出与该技术方案类似的结构或实施例，均属于本发明保护范围。

Claims (1)

1.一种基于闭环控制原理的电容层析成像图像重建方法，其特征在于，该方法具体操作步骤如下：

步骤一、对同一圆形截面上有等面积均匀分布N个电极的电容层析成像传感器，采用传统的相邻激励测量模式，通过测量分别得到空场条件下的N(N-1)/2个独立测量值和有物场条件下的N(N-1)/2个独立测量值，所述空场条件表示场域内部只有介电常数均匀分布的背景介质，所述有物场条件表示在均匀分布的背景介质上存在与背景介质介电常数不同的物体，空场条件下与有物场条件下的测量值可分别建立N×N的电容矩阵C₀和C_m，C₀表示空场电容矩阵，C_m表示有物场电容矩阵，C_m与C₀的差表示为：

其中，ΔC_s,t为第s个与第t个电极之间的电容的变化量(s≠t)，ΔC_s,t＝ΔC_t,s，ΔC_s,s为第s个电极上的自电容变化量，所有电极的电容变化量的关系为：

利用Calderon算法计算介电常数分布，即利用公式(3)计算散射变换t(re^iθ)：

其中，(r,θ)表示极坐标参数，r为极径，θ为极角，π表示圆周率，e表示自然常数，θ_s表示第s个电极的中心位置，θ_s＝2πs/N，A表示电极面积；

利用公式(4)可反演出场域内部物质介电常数分布情况：

δε(x,y)表示坐标为(x,y)的点处介电常数的变化值，R为数值积分区域的半径；

步骤二、根据所求得的分布得到新的电容矩阵C'_m，将该电容矩阵与步骤一中测量得到的空场电容矩阵C₀相减得到ΔC'_N×N(k)作为负反馈信息：

其中，ΔC'_N×N(k)表示第k次迭代计算的负反馈电容变化量矩阵；

步骤三、计算第k次偏差p_N×N(k)＝ΔC_N×N-ΔC'_N×N(k)，利用公式(6)得到PID控制器第k次输出信息q_N×N(k)，接着利用公式(3)和(4)所描述的Calderon算法更新重建分布，选取的数字式PID控制器为：

q(k)＝q(k-1)+(K_p+K_i+K_d)p(k)-(K_p+2K_d)p(k-1)+K_dp(k-2) (6)

其中，K_p,K_i,K_d分别为比例、积分、微分系数；q(k),q(k-1)分别为PID控制器第k次迭代计算和第(k-1)次迭代计算的输出，p(k),p(k-1),p(k-2)分别为PID控制器第k次迭代计算、第(k-1)次迭代计算和第(k-2)次迭代计算的输入，k＝1,2,3,...,m，m为预设迭代次数,设定初值为q_N×N(0)＝O_N×N,p_N×N(0)＝O_N×N,p_N×N(-1)＝O_N×N，ΔC'_N×N(1)＝O_N×N,O_N×N表示N×N的全零矩阵；

步骤四、判断迭代计算次数k是否大于预设迭代次数m，如果是，则执行步骤五，如果否，则令k＝k+1并返回步骤二；

步骤五、结束迭代，输出成像结果。