CN113034635B - 一种抑制工业成像阶梯伪影的图像重建方法 - Google Patents

Abstract

本发明公开一种抑制工业成像阶梯伪影的图像重建方法，结合一种新颖的处理函数对电导率分布系数进行处理，最后利用重建电导率分布的算子对最优电导率分布系数进行重建得到最优电导率分布。稀疏正则项可以用更少的数据来还原最真实的解，将优化模型中无约束的问题转化为相应的约束问题，并对其子问题进行迭代求解，直到满足指定迭代次数即可结束迭代，求得此时的最优电导率分布系数。最后利用重建电导率分布的算子得到最优电导率分布并进行成像。本发明有效地减少了TV正则化方法重建图像背景中的阶梯伪影，提高了重建图像的质量，并且具有相对良好的抗噪性，为工业成像领域提供了新的可能。

Description

一种抑制工业成像阶梯伪影的图像重建方法

技术领域

本发明属于电学层析成像技术领域，具体涉及一种抑制工业成像阶梯伪影的图像重建方法，以实现工业两相流/多相流的电导率分布图像重建。

背景技术

电学层析成像(Electrical Tomography,ET)的工作原理是：利用基于电极敏感阵列获取被测区域物质空间分布信息，以电学信号作为载体进行处理与传输，采用信息重建算法，重建被测区域物质空间分布的全部信息。ET技术包括三种工作模态：电阻层析成像(Electrical Resistance Tomography,ERT)、电容层析成像(Electrical CapacitanceTomography,ECT)和电磁层析成像(Electromagnetic Tomography,EMT)，其研究对象涉及了电导率、电容率、磁导率等主要电学参数。ET作为一种可视化测量手段，以其非侵入、易携带、低成本、响应快等技术优势，受到了科研人员的广泛关注。其中ERT可视化技术能够使工业测量信息二维化，进一步促进了过程参数检测技术的发展。

ERT图像重建是一个病态的逆问题，它导致了重建图像质量低的问题，这一问题阻碍了ERT的应用和发展。正则化方法是解决逆问题病态性的有效方法，其中基于L2范数的Tikhonov正则化(Tikhonov Regularization,TR)方法和基于L1范数的总变分(Totalvariation regularization method,TV)正则化方法是两种常用的正则化方法。Tikhonov正则化方法针对电导率连续分布的被测物体具有较好的性能。但是当被测物体边界的电导率发生急剧变化时，会出现重建图像质量较差的问题。相比之下，TV正则化方法具有良好的保边性，可以有效地提高重建图像的质量。然而，TV方法的重建图像会在背景区域产生“阶梯伪影”。阶梯伪影的出现，严重影响了重建图像的质量。因此，本发明提出了一种抑制工业成像阶梯伪影的图像重建方法，以减少重建图像中的阶梯伪影，进而提高重建图像的质量。

发明内容

本发明的目的是针对TV正则化方法的ERT图像重建问题，提出了一种抑制工业成像阶梯伪影的图像重建方法，用于优化ERT重建工业两相流/多相流图像的电导率分布。为了去除TV正则化方法产生的阶梯伪影，提出处理函数对电导率分布系数进行处理。另外，为了促进处理函数的稀疏性，给出非凸函数作为处理函数的非凸正则项；同时，为了保证处理函数有唯一的最优解，对非凸函数中限制参数的范围进行了约束，以保证处理函数的凸性；更进一步，为了求解处理函数的最优解，将无约束的优化问题转化为约束问题并对其子问题进行迭代求解；最后，对求解得出的最优电导率分布系数进行逆变换，以得到最优电导率分布。与Tikhonov正则化方法和TV正则化方法相比，本发明提出的一种抑制工业成像阶梯伪影的图像重建方法有效地改善了重建图像背景中的阶梯伪影，明显地提高了重建图像的质量。

本发明为解决上述技术问题采用如下技术方案：一种抑制工业成像阶梯伪影的图像重建方法，其特征在于具体过程为：首先利用TV正则化方法求得初始电导率分布

TV正则化方法将工业两相流/多相流电导率分布图像重建过程看作一个线性不适定的逆问题Sg＝b，式中，S为灵敏度矩阵，b为边界相对电压测量值，g为电导率分布，求解初始电导率分布

的目标函数如下：

式中，λ是正则化参数，用来平衡保真项

和惩罚项

β是一个正常数，L_l为对应的电导率分布g的稀疏矩阵，采用牛顿法求解J(g)的极小值，以得到初始电导率分布

然后，利用处理函数中的分解电导率分布的算子将

分解为电导率分布系数，并进行处理，提出的处理函数F(x)为：

式中，

是数据保真项，x是电导率分布系数；

是初始电导率分布；W是分解电导率分布的算子，W^T是重建电导率分布的算子，是W的逆运算，

非凸正则项，其中λ_e是非凸正则化参数，

是非凸函数，x_e,w是当分解尺度为e，时间标度为w时的电导率分布系数，a_e是约束非凸函数

的限制参数；λ_o||DW^Tx||₁是凸正则项，λ_o是凸正则化参数，D是一个(n-1)×n的差分矩阵；

非凸函数

需要具有较强的稀疏性，并通过对限制参数a_e的约束，以保证处理函数F(x)的凸性，从而确保处理函数有唯一且最小的解。本发明提出的非凸函数

的具体形式为：

而后，为了获得最优电导率分布

需要将F(x)最小化来求得最优电导率分布系数

并利用重建电导率分布的算子W^T重建最优电导率分布

其优化模型为：

为了求解最优电导率分布系数

需要将优化模型中无约束的问题转化为相应的约束问题：

s.t.r＝x

式中，定义r为约束变量；U₁(x)和U₂(r)分别为：

U₂(r)＝λ_o||DW^Tr||₁

约束问题的增广拉格朗日函数可以表示为：

式中，μ是拉格朗日参数，μ＝0.7；d为迭代变量；

根据增广拉格朗日函数，将约束问题分解成下列子问题：

式中，上标K代表迭代次数，K的初始值设为0，设定K的最大值为50；r^K表示r第K次迭代的值，r的初始值

d^K表示d第K次迭代的值，d的初始值d⁰＝0。

重建方法的整体步骤为：

(1)首先在没有物体的电阻层析成像测量场域，采用相邻模式对电极对进行电流激励和电压测量，以获得空场边界电压测量值b₁。再将物体置于测量场域内，获得此时的满场边界电压测量值b₂。根据测量得到的b₁和b₂，可以求得计算所需的边界相对电压测量值b,即b＝b₂-b₁。根据灵敏度理论得到灵敏度矩阵S。

(2)求解初始电导率分布

其中正则化参数λ＝0.68；正常数β＝0.0001；k为电导率分布g求解的循环次数，初始值为0，设定k＝10。

(3)提出处理函数F(x)，给出非凸函数作为非凸正则项，同时约束限制参数的范围(0≤a_e＜1/λ_e)以保证F(x)的凸性。

(4)对F(x)中的参数进行赋值。非凸正则化参数λ_e＝0.95，凸正则化参数λ_o＝0.16，限制参数a_e＝0.99/λ_e。

(5)给出求解最优电导率分布系数

的优化模型。

(6)将无约束的优化问题转化为有约束的问题，并对其子问题进行迭代求解，直到满足迭代次数，停止迭代，求得最优解。

(7)用求解得到的最优电导率分布

进行成像。

本发明的显著效果是：一种抑制工业成像阶梯伪影的图像重建方法，该方法在求解初始电导率分布

的基础上，通过优化初始电导率分布

减少阶梯伪影以提高重建图像的质量。相比于Tikhonov正则化方法和TV正则化方法，本发明提出的一种抑制工业成像阶梯伪影的图像重建方法在减少阶梯伪影和提升图像质量等方面有着良好的效果。

附图说明

图1为一种改善工业两相流/多相流重建图像背景的正则化方法的流程图。

图2为本发明的电阻层析成像系统圆形单截面被测场域、激励电流、测量电压的模式以及电极分布。

图3、图4为本发明的实例，选取了2种典型模型(模型A-B)，分别用三种方法：Tikhonov正则化方法、TV正则化方法和本发明的一种抑制工业成像阶梯伪影的图像重建方法，设置不同的噪声环境(无噪环境和5％噪声水平的环境)，进行图像重建，重建的图像结果如图所示。

图5、图6分别为无噪环境和5％噪声水平的环境下相对重建模型空间分辨率的定量评价——相对误差和相关系数的值。

具体实施方式

结合附图和实例对本发明的一种抑制工业成像阶梯伪影的图像重建方法加以说明。

本发明的一种抑制工业成像阶梯伪影的图像重建方法，用于优化重建工业两相流/多相流的电导率分布。在求解初始电导率分布

的基础上，通过优化初始电导率分布

在减少重建图像背景的阶梯伪影、提高重建图像质量的同时，保留重建图像的锐利边缘。

如图1所示，为本发明的一种抑制工业成像阶梯伪影的图像重建方法流程图。

如图2所示，为电学层析成像中的电阻层析成像系统圆形单截面被测场域、激励电流和测量电压的模式以及电极分布，16电极均匀分布在场域外壁。

为了直观地体现本发明一种抑制工业成像阶梯伪影的图像重建方法的效果，分别用Tikhonov正则化方法、TV正则化和本方法，对两种典型模型进行了图像重建。

工业两相流/多相流电导率分布图像重建过程是一个严重不适定性的逆问题。逆问题可以用最小二乘优化形式的目标函数表示：

式中，f(g)为目标函数。正则化方法是解决逆问题不适定性的有效方法。正则化方法的一般形式可以表述为：

其中，λ是一个正则化参数，用来平衡保真项

和正则化项R(g)。

标准的Tikhonov正则化方法可将正则化项R(g)替换为正则化项

可以描述为：

其中N为单位矩阵。然而，Tikhonov方法在被测介质不连续分布时，边界上会出现过度平滑的现象，导致重建图像不够准确，空间分辨率降低。为了解决这一问题，提出了TV正则化方法。TV正则化方法可以表述为：

这一方法有效地改善了重建图像中锐利物体边缘过于平滑的现象。与此同时，TV正则化方法中产生了分段常数解，导致重建图像的光滑区域产生“阶梯伪影”，导致重建图像的质量受到影响。为改善TV正则化方法的阶梯伪影，提高重建图像质量，本发明在求解初始电导率分布

的基础上，进一步优化初始电导率分布

提出了一种抑制工业成像阶梯伪影的图像重建方法。

具体实施步骤如下：

步骤一：首先在没有物体的电阻层析成像测量场域，采用相邻模式对电极对进行电流激励和电压测量，以获得空场边界电压测量值b₁。再将物体置于测量场域内，获得此时的满场边界电压测量值b₂。根据测量得到的b₁和b₂，可以求得计算所需的边界相对电压测量值b,即b＝b₂-b₁。

灵敏度矩阵S是根据灵敏度理论得到的，灵敏度矩阵中的元素称为灵敏度系数，灵敏度系数的计算公式为：

式中，s_ij是第j个电极对对第i个电极对的灵敏度系数，φ_i、φ_j分别为第i个电极对及第j个电极对在激励电流为I_i、I_j时的场域电势分布。

步骤二：初始电导率分布

计算：

式中，λ是正则化参数(取λ＝0.68)，可以平衡保真项

和惩罚项

β是一个正常数(取β＝0.00001)，L_l为对应的电导率分布g的稀疏矩阵。

采用牛顿法求解J(g)的极小值，J(g)的梯度函数可以表示为：

J′(g)＝S^T(Sg-b)+λL_β(g)g

其中

式中，S^T是S的转置，L^T是L的转置，diag表示对角矩阵。

J(g)的Hessian矩阵可以表示为：

H(g)＝S^TS+λL_β(g)

则电导率分布g的解为：

g_k+1＝g_k-H(g)^-1J′(g)

式中，k为电导率分布g求解循环次数，初始值为0，设定k最大值为10；H(g)^-1是H(g)的逆矩阵。通过循环求解，可得初始电导率分布

为：

步骤三：提出处理函数F(x)：

式中，

是数据保真项，x是电导率分布系数；

是初始电导率分布；W是分解电导率分布的算子，W^T是重建电导率分布的算子，是W的逆运算。

非凸正则项，其中

是非凸函数，x_e,w是当分解尺度为e，时间标度为w时的电导率分布系数，a_e是约束非凸函数

的限制参数；λ_o||DW^Tx||₁是凸正则项，D是一个(n-1)×n的差分矩阵；λ_e和λ_o分别是非凸正则化参数和凸正则化参数。

非凸函数

需要具有较强的稀疏性，并通过对限制参数a_e的约束(约束范围为0≤a_e＜1/λ_e)，以保证处理函数F(x)的凸性，从而确保处理函数有唯一且最小的解。本发明给出的非凸函数

为：

步骤四：对F(x)中的参数进行赋值。非凸正则化参数λ_e＝0.95，凸正则化参数λ_o＝0.16，限制参数a_e＝0.99/λ_e。

步骤五：为了获得最优电导率分布

需要将F(x)最小化来求得最优电导率分布系数

并利用重建电导率分布的算子W^T重建最优电导率分布

其优化模型为：

步骤六：对优化模型进行求解。为了求解最优电导率分布系数

需要将优化模型中无约束的问题转化为相应的约束问题：

s.t.r＝x

式中，定义r为约束变量；U₁(x)和U₂(r)分别为：

U₂(r)＝λ_o||DW^Tr||₁

约束问题的增广拉格朗日函数可以表示为：

式中，μ是拉格朗日参数，μ＝0.7；d为迭代变量。

根据增广拉格朗日函数，将约束问题分解成下列子问题：

式中，K是迭代次数，K的初始值设为0，K的最大值为50；r^K表示r第K次迭代的值，r的初始值

d^K表示d第K次迭代的值，d的初始值d⁰＝0。

由增广拉格朗日函数，再结合L2范数的定义，a^K+1的解可以写为如下形式：

式中，p_e,w是第一辅助变量，

为了求解r^K+1，利用邻近算子及其半正交线性变换，定义第二辅助变量v，用v^K表示算法第K次迭代时v值，并设v^K＝x^K+1-d^K。根据子问题r^K+1，v^K的邻近算子可定义为：

式中，h(r)＝φ(W^Tr)，

z是第三辅助变量。利用邻近算子及其半正交线性变换，W^Tv^K的邻近算子定义为：

则r^K+1的解最终表示为：

式中，TVD表示TV去噪方法，可通过绷弦算法计算。

当上述迭代次数K满足K＝50时，迭代过程终止。根据最终求出的解，对最优电导率分布系数

进行逆变换，即

从而得到最优电导率分布

步骤七：用求解得到的最优电导率分布

进行成像。

图3、图4为本发明的实例，选取2种典型的模型，分别利用Tikhonov正则化方法、TV正则化方法和本方法，对图像进行重建。从结果中可以看出，Tikhonov正则化方法重建的图像锐利边缘过于平滑，呈现出的图形形状和大小不够准确。相比之下，TV方法的重建图像保留了锐利边缘，重建的形状和大小也更加准确，重建效果更好。然而，重建图像的背景中却出现了冗余的阶梯伪影。对比Tikhonov正则化方法和TV正则化方法，本方法在重建过程中能够有效地减少阶梯伪影，提高重建图像质量。该方法的重建图像能够更加准确地检测工业两相流/多相流。

在工业成像测量领域中，用ERT技术重建工业两相流/多相流具有非常重要的意义。实际环境中，各种噪声对测量的影响是不可避免的。从图4所示的5％噪声水平情况下的重建图像可以看出，相比于Tikhonov正则化方法和TV正则化方法，本发明所提出的方法抗噪性更强，在噪声环境下依然能够呈现出相对质量较高的重建图像。这也表明本方法在实际应用中的可行性和实用性更好。

为了更好地评价该方法的性能效果，采用相对误差(Relative Error，RE)和相关系数(Correlation Coefficient，CC)对重建图像的质量进行定量评估：

图像的相对误差越小、相关系数越大，表明重建图像的空间分辨率越高，质量越好。式中，g′是重建区域的计算电导率分布，g^*是实际电导率分布，g′_q和

分别是第q个网格的计算电导率分布和实际电导率分布，

和

分别是第q个网格的平均计算电导率分布和平均实际电导率分布。图5和图6给出了Tikhonov、TV和本发明所提出的三种不同的正则化方法在无噪环境下和5％噪声水平环境下，2种典型模型重建图像的相对误差和相关系数。数据结果表明，与Tikhonov正则化方法和TV正则化方法相比，本发明提出的一种抑制工业成像阶梯伪影的图像重建方法，具有最低的相对误差和最高的相关系数，在无噪环境和5％噪声水平情况下，能够较为准确地重建出被测区域内的电导率分布，有效地减少了TV正则化方法重建图像中的阶梯伪影，改善了重建图像的质量。

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

Claims (1)

1.一种抑制工业成像阶梯伪影的图像重建方法，其特征在于具体过程如下：

步骤一：首先在没有物体的电阻层析成像测量场域，采用相邻模式对电极对进行电流激励和电压测量，以获得空场边界电压测量值b₁，再将物体置于测量场域内，获得此时的满场边界电压测量值b₂，根据测量得到的b₁和b₂，求得计算所需的边界相对电压测量值b，即b＝b₂-b₁；

灵敏度矩阵S是根据灵敏度理论得到的，灵敏度矩阵中的元素称为灵敏度系数，灵敏度系数的计算公式为：

式中，s_ij是第j个电极对对第i个电极对的灵敏度系数，φ_i、φ_j分别为第i个电极对及第j个电极对在激励电流为I_i、I_j时的场域电势分布；

步骤二：计算初始电导率分布

式中，λ是正则化参数，取λ＝0.68，能够平衡保真项

和惩罚项

β是一个正常数，取β＝0.00001，L_l为对应的电导率分布g的稀疏矩阵；

采用牛顿法求解J(g)的极小值，J(g)的梯度函数表示为：

J'(g)＝S^T(Sg-b)+λL_β(g)g

其中

式中，S^T是S的转置，L^T是L的转置，diag表示对角矩阵；

J(g)的Hessian矩阵表示为：

H(g)＝S^TS+λL_β(g)

则电导率分布g的解为：

g_k+1＝g_k-H(g)^-1J'(g)

式中，k为电导率分布g求解的循环次数，k初始值为0，设定k的最大值为10；H(g)^-1是H(g)的逆矩阵，通过循环求解，得初始电导率分布

为：

步骤三：设计处理函数F(x)：

式中，

是数据保真项，W是分解电导率分布的算子，x是电导率分布系数；W^T是重建电导率分布的算子，是W的逆运算，

非凸正则项，其中λ_e是非凸正则化参数，

是非凸函数，x_e,w是当分解尺度为e，时间标度为w时的电导率分布系数，a_e是约束非凸函数

的限制参数；

是凸正则项，λ_o是凸正则化参数，D是一个(n-1)×n的差分矩阵；

非凸函数

需要具有较强的稀疏性，并通过对限制参数a_e的约束，约束范围为0≤a_e<1/λ_e，以保证处理函数F(x)的凸性，从而确保处理函数有唯一且最小的解，给出非凸函数

为：

步骤四：对F(x)中的参数进行赋值，非凸正则化参数λ_e＝0.95，凸正则化参数λ_o＝0.16，限制参数a_e＝0.99/λ_e；

步骤五：为了获得最优电导率分布

需要将F(x)最小化来求得最优电导率分布系数

并利用重建电导率分布的算子W^T重建最优电导率分布

其优化模型为：

步骤六：对优化模型进行求解，为了求解最优电导率分布系数

需要将优化模型中无约束的问题转化为相应的约束问题：

s.t.r＝x

式中，定义r为约束变量；U₁(x)和U₂(r)分别为：

U₂(r)＝λ_o||DW^Tr||₁

约束问题的增广拉格朗日函数表示为：

式中，μ是拉格朗日参数，μ＝0.7；d为迭代变量；

根据增广拉格朗日函数，将约束问题分解成下列子问题：

式中，上标K表示迭代次数，K的初始值设为0，设定K最大值为50；r^K表示r第K次迭代的值，r的初始值

d^K表示d第K次迭代的值，d的初始值d⁰＝0；

由增广拉格朗日函数，再结合L2范数的定义，x^K+1的解写为如下形式：

式中，p_e,w是第一辅助变量，

为了求解r^K+1，利用邻近算子及其半正交线性变换，定义第二辅助变量v，用v^K表示算法第K次迭代时v值，并设v^K＝x^K+1-d^K，根据子问题r^K+1，v^K的邻近算子定义为：

式中，h(r)＝φ(W^Tr)，

z是第三辅助变量，利用邻近算子及其半正交线性变换，W^Tv^K的邻近算子定义为：

则r^K+1的解最终表示为：

式中，TVD表示TV去噪方法，通过绷弦算法计算；

当上述迭代次数K满足K＝50时，迭代过程终止，根据最终求出的解，对最优电导率分布系数

进行逆变换，即

从而得到最优电导率分布

步骤七：用求解得到的最优电导率分布

进行成像。

Images