CN114596377A - 一种多指数特征提取的ect图像重建方法 - Google Patents

Abstract

一种多指数特征提取的ECT图像重建方法，首先将图像灰度矩阵和灵敏度矩阵经过零向量扩展来重组，构建成一个稀疏化的稀疏观测方程，其次将Tikhonov正则化算法得到灰度向量建模成一个有限新息率信号，用指数再生采样核对其滤波后均匀采样，提取特征信息，从采样样本中计算获取测量值后，把稀疏观测矩阵与FRI观测矩阵相结合，稀疏观测向量和FRI观测向量相结合,随机重组综合观测矩阵和综合观测向量的行向量，获得综合观测方程；将电容层析图像重构的问题转化为L0范数最小化问题，求解L0范数最小化问题得到稀疏向量的估计，计算出灰度向量的估计，最后重构原始图像。本发明具有良好的图像质量和较高的重建精度。

Description

一种多指数特征提取的ECT图像重建方法

技术领域

本发明涉及图像处理技术领域，具体涉及一种多指数特征提取的 ECT图像重建方法。

背景技术

电容层析成像技术(Electrical Capacitance Tomography，ECT)，是一种以电容传感器作为测量电极的过程成像技术(Process Tomography)。随着工业的发展，工业上对管道内部截面两相流的流动参数可视化检测的需求越来越大，近年来ECT技术已经在工业上的过程成像中广泛应用，是过程层析成像技术的研究热点。并且ECT技术是基于电容传感器的成像技术，具有设备成本低、非侵入性、结构简单、安全等许多的优点，能有效的解决两相流的测量问题。

电容层析成像的基本思想是：在实际情况中，管道中的流体不是单一的，通常由两种或两种以上的物质混合而成，称为两相流和多相流。每一种流体有不同的介电常数。两相流的介电常数会随着流体分布的改变而发生变化。在被测管道周围设置一系列电容传感器，用于获取测量值。数据采集系统把电容传感器获取的电容测量值上传到成像计算机，通过图像重建算法将管道内部截面的介电常数分布转换为图像。典型的ECT系统结构图如1所示，100为计算机，200是检测线圈，300是屏蔽层，400是励磁线圈。

电容传感器i-j之间的灵敏度矩阵可以用公式(1)描述为：

其中，C_i,j(k)表示的是两相流管道截面上第k个像素上的介电常数为ε_h，剩下的像素上介电常数为ε_l时，电容传感器i-j之间的电容测量值， k＝1,2,...,m。μ(k)是修正因子。

分别代表的是两相流管道截面上所有像素都填入介电常数ε_h和ε_l时，电容传感器i-j之间的电容测量值。

被测管道横截面两相流体介电常数分布ε(x,y)与从电容传感器获取的不同测量电极之间的电容向量C_i之间的关系可以表示为：

C_i＝∫∫_DS_i(x,y)ε(x,y)dxdy (2)

其中，D表示被检测的管道区域，(x,y)∈D表示像素点的取值都在被检测管道截面上；ε(x,y)为介电常数分布函数，C_i(i＝1,2...m)为电容传感器测量所得的m组电容测量值，m为电容传感器之间两两组合数， m为整数。S_i(x,y)是灵敏度分布函数，对应于C_i。

由公式(2)可以看出，该公式中的C_i和ε(x,y)不是线性关系，只有当介电常数区别不大的时候，才可以近似为线性关系。首先对式(2)进行离散处理，再对其线性处理，最后对其归一化，可以得到ECT的矩阵形式数学模型：

C＝Sg (3)

其中，C是归一化的电容测量值，是一个大小为m×1的电容向量。g 是对应于管道截面介电常数分布的大小为n×1归一化灰度值向量。m 为测量极板间两两组合个数，m为整数，n为所剖分的像素个数，n为整数。S是大小为m×n归一化的灵敏度矩阵。ECT图像重建的问题是已知测量电容向量C求解灰度向量g的过程，从而得出管道内的介电常数分布。结果通常以图像呈现，因此这种过程称为图像重构。

由公式(3)可以看出，只有当m＝n时，S的逆矩阵才会存在。在实际的ECT系统中，电容传感器的数量一般很少，导致获得的测量值数量m远远小于剖分的像素个数n，导致S不可逆，要在测量电容值不足的情况下求解方程。这是一个欠定的问题，欠定问题不具有唯一解。电容层析成像的关键是图像重构算法，目前已经提出了很多比较成熟的图像重构算法，如线性反投影法(LBP)、Landweber算法、Tikhonov正则化算法等。LBP算法对于复杂的流型，重建图像的质量相对较低。但是一种简单快速的算法；与LBP算法相比，Landweber 算法具有更好的图像质量，近几年来得到了广泛的应用，但是这种算法核心是最速下降法，在步长参数不适当的情况下，会导致结果的分歧。步长过大。可能会找不到最优点或者陷入局部极值；步长过小，需要迭代的次数过多，耗费的时间太长。Tikhonov正则化算法可以得到一幅近似稳定的图像，但图像的边缘并不理想。迄今为止，提高电容层析图像重建的质量和精度仍然是个关键技术问题。

发明内容

为了克服已有技术的不足，针对电容层析成像问题，本发明提出一种多指数特征提取的ECT图像重建方法，将图像灰度矩阵和灵敏度矩阵经过零向量扩展来重组，构建一个稀疏化的稀疏观测方程。然后将ECT的图像的灰度向量建模成一个有限新息率(FiniteRate of Innovation，FRI)信号，用指数再生采样核对其滤波后均匀采样，提取输入信号中的特征信息，从采样样本中计算获取测量值后，把稀疏观测矩阵与FRI观测矩阵相结合，稀疏观测向量和FRI观测向量相结合，并随机重组综合观测矩阵和综合观测向量的行向量，获得综合观测方程。将电容层析图像重构的问题转化为L0范数的最小化问题，为了得到稀疏向量的估计，最后通过求解L0范数最小化问题，估计出原始图像的灰度。

本发明解决其技术问题所采用的技术方案是：

一种多指数特征提取的ECT图像重建方法，包括以下步骤：

步骤一，初始化数据，获取并归一化数据，ECT图像重构的问题是通过测量电容向量C求解灰度向量g的过程，估计管道内的介电常数分布；

步骤二，补零重组，构建稀疏化的稀疏观测方程。在实际的ECT 系统中，电容传感器的数量一般不会很多，所剖分的像素个数n远大于测量电容值个数m，从而导致采样率太低，降低重构的精度；通过零向量扩展的方式增加虚拟电极，以获取更多的测量电容值，构建一个稀疏化的稀疏观测方程，ECT系统的数学模型写为：

λ₀＝S₀g (4)

其中，λ₀是通过补零重组的方法扩展后的电容向量；S₀是通过补零重组的方法扩展后的灵敏度矩阵；

步骤三，有限新息率FRI信号建模，由于灰度向量取值范围为[0,1]，是一个典型的离散FRI信号，因此把Tikhonov正则化算法所得到的灰度向量g_v建模为狄拉克脉冲序列信号g(x)：

g_v＝(S^T·S+α·I)^-1·S^TC (5)

其中，g_v为大小为n×1灰度值向量，α是正则化参数，I是一个大小为 n×n的单位阵，由于S^T·S存在不可逆的可能，所以加入参数α用于解决S^T·S不适定的问题，正则化算法可以有效解决灵敏度矩阵不可逆的病态问题；

其中，x＝1,2...,N是像素的位置，L是g(x)非零灰度值的个数，a_l∈[0,1] 是像素x_l的灰度值，x_l∈{1,2,....,N}是非零灰度值的像素位置；

步骤四，有限新息率FRI采样，首先使用指数再生采样核

对狄拉克脉冲序列信号g(x)滤波，再对其均匀采样，可以得到K个样本值y_k(k＝1,2,...,K)；

均匀采样后得到K个样本y_k：

其中，y(x)＝g(x)*h(x)是FRI信号g(x)滤波后的结果，＜·,·＞表示内积， y_k是滤波后均匀采样的样本值，T是采样间隔，K是采样值的数量；

步骤五，从FRI采样样本中获取测量值，构建FRI观测方程。用样条系数C_m,k对FRI的采样值y_k进行加权后求和，可以计算出M个复指数

K个测量值可以由K个样本y_k计算得到：

公式(8)写成矩阵形式：

因为像素位置x_l属于集合{1,2,...,N}，a_l∈[0,1]是像素x_l的灰度值；所以公式(8)写成：

U＝Ag (11)

其中，U＝[τ₀,τ₁,...,τ_M-1]^T∈R^M×1是FRI测量向量，g＝[g₁,g₂,...,g_N]^T是 ECT系统中大小为N×1的灰度向量，A是由集合

构成M×N的FRI测量矩阵；

步骤六，构建综合观测方程，步骤如下：

步骤6.1：稀疏观测矩阵结合FRI观测矩阵，由于公式(4)和(10) 都是对g的一个多维观测，所以两个观测方程可以结合，得到一个新的综合观测方程，结合公式(4)和(11)得到新的综合观测方程：

步骤6.2：以相同规则随机重组综合观测向量和综合观测矩阵的行向量。因此ECT系统的数学模型重写为：

λ_new＝S_newg (13)

其中，λ_new是综合观测向量；S_new是综合观测矩阵；

步骤七，求综合观测方程的稀疏解，步骤如下：

步骤7.1：稀疏变换，只有信号是稀疏的，才能满足压缩感知的要求，对于ECT系统，信号稀疏程度大部分情况下不能满足压缩感知的要求，为了确保输入信号能够适应压缩感知，需要对输入信号进行正交变换，将输入信号转化为稀疏信号，正交变换由以下公式表示：

g＝ψs (14)

其中，大小为N×N的矩阵ψ为稀疏基，原始信号g投影到稀疏基上得到大小为N×1的稀疏向量s；

步骤7.2：求稀疏解，电容层析图像重构问题转化为L0范数优化问题，由于信号s是稀疏的，获得稀疏解

最直接的方法就是求解L0 范数优化问题，用OMP算法：

其中，L0范数||s||₀表示向量s中非零系数的数目，通过正交匹配追踪算法OMP来求解，大小为N×N的矩阵ψ为稀疏基；

步骤八，图像重建，原始图像信号被估计为：

进一步，所述步骤四中，用指数再生采样核

对狄拉克脉冲序列信号g(x)滤波后均匀采样，得到样本y_k(k＝1,2,...,K)；

M-1阶指数再生采样核的公式如下：

其中，M-1＝3是再生采样核的阶数，M-1阶再生采样核可以再生出M 个指数

C_m,k由下式得出：

其中

是

的对偶函数，属于准正交函数。

再进一步，所述步骤五中，从FRI采样样本中计算测量值，构建 FRI观测方程，指数再生采样核的阶数和FRI信号的关系为M-1＝2L-1， L为狄拉克脉冲个数，M-1阶指数再生采样核可以再生出M个复指数

其中，α_m＝a₀+jmλ，a₀为自由设定的复指数，λ为自由设定的实数。

更进一步，所述步骤7.2中，为了保证ECT系统的原始信号能够重构，需要对原始信号进行正交变换，使输入信号变得稀疏化，稀疏基为离散余弦变换基。

本发明提出了一种多指数特征提取的ECT图像重建方法。首先，将图像灰度矩阵和灵敏度矩阵经过零向量扩展来重组，构建一个稀疏化的稀疏观测方程，通过这种方法来增加测量电容值的数量，可以提高信号的采样率，提高重构信号的精度。其次，将Tikhonov正则化算法得到的灰度向量建模成一个有限新息率(Finite Rate of Innovation， FRI)信号，用指数再生采样核对其滤波后均匀采样，提取输入信号中的特征信息，从采样样本中计算获取测量值后，把稀疏观测矩阵与 FRI观测矩阵相结合，稀疏观测向量和FRI观测向量相结合，并随机重组综合观测矩阵和综合观测向量的行向量，获得综合观测方程。将图像重构的问题转化为L0范数最小化问题，求解L0范数最小化问题得到稀疏向量的估计。最后，重构出原始图像的灰度向量。

本发明的有益效果主要表现在：提高信号的采样率，提高重构信号的精度。

附图说明

图1是ECT系统结构图。

图2是多指数特征提取的ECT图像重建方法的原理框图。

图3是经典的FRI采样结构图。

图4是二种流型的原始图像，分别为核心流和半管流。

图5是LBP算法重建的图像效果图。

图6是Landweber算法重建的图像效果图。

图7是Tikhonov算法重建的图像效果图。

图8是本发明方法重建的图像效果。

具体实施方式

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

参照图2～图8，一种多指数特征提取的ECT图像重建方法，包括以下步骤：

步骤一，初始化数据，获取并归一化数据，ECT图像重构的问题是通过测量电容向量C求解灰度向量g的过程，估计管道内的介电常数分布；

步骤二，补零重组，构建稀疏化的稀疏观测方程。在实际的ECT 系统中，电容传感器的数量一般不会很多，所剖分的像素个数n远大于测量电容值个数m，从而导致采样率太低，降低重构的精度；通过零向量扩展的方式增加虚拟电极，以获取更多的测量电容值，构建一个稀疏化的稀疏观测方程，ECT系统的数学模型写为：

λ₀＝S₀g (4)

其中，λ₀是通过补零重组的方法扩展后的电容向量；S₀是通过补零重组的方法扩展后的灵敏度矩阵；

步骤三，有限新息率FRI信号建模，由于灰度向量取值范围为[0,1]，是一个典型的离散FRI信号，因此把Tikhonov正则化算法所得到的灰度向量g_v建模为狄拉克脉冲序列信号g(x)：

g_v＝(S^T·S+α·I)^-1·S^TC (5)

其中，g_v为大小为n×1灰度值向量，α是正则化参数，I是一个大小为 n×n的单位阵，由于S^T·S存在不可逆的可能，所以加入参数α用于解决S^T·S不适定的问题，正则化算法可以有效解决灵敏度矩阵不可逆的病态问题；

其中，x＝1,2...,N是像素的位置，L是g(x)非零灰度值的个数，a_l∈[0,1] 是像素x_l的灰度值，x_l∈{1,2,....,N}是非零灰度值的像素位置；

步骤四，有限新息率FRI采样，首先使用指数再生采样核

对狄拉克脉冲序列信号g(x)滤波，再对其均匀采样，可以得到K个样本值y_k(k＝1,2,...,K)；

均匀采样后得到K个样本y_k：

其中，y(x)＝g(x)*h(x)是FRI信号g(x)滤波后的结果，＜·,·＞表示内积， y_k是滤波后均匀采样的样本值，T是采样间隔，K是采样值的数量；

步骤五，从FRI采样样本中获取测量值，构建FRI观测方程。用样条系数C_m,k对FRI的采样值y_k进行加权后求和，可以计算出M个复指数

K个测量值可以由K个样本y_k计算得到：

公式(8)写成矩阵形式：

因为像素位置x_l属于集合{1,2,...,N}，a_l∈[0,1]是像素x_l的灰度值；所以公式(8)写成：

U＝Ag (11)

其中，U＝[τ₀,τ₁,...,τ_M-1]^T∈R^M×1是FRI测量向量，g＝[g₁,g₂,...,g_N]^T是ECT系统中大小为N×1的灰度向量，A是由集合

构成M×N的FRI测量矩阵；

步骤六，构建综合观测方程，步骤如下：

步骤6.1：稀疏观测矩阵结合FRI观测矩阵，由于公式(4)和(10) 都是对g的一个多维观测，所以两个观测方程可以结合，得到一个新的综合观测方程，结合公式(4)和(11)得到新的综合观测方程：

步骤6.2：以相同规则随机重组综合观测向量和综合观测矩阵的行向量。因此ECT系统的数学模型重写为：

λ_new＝S_newg (13)

其中，λ_new是综合观测向量；S_new是综合观测矩阵；

步骤七，求综合观测方程的稀疏解，步骤如下：

步骤7.1：稀疏变换，只有信号是稀疏的，才能满足压缩感知的要求，对于ECT系统，信号稀疏程度大部分情况下不能满足压缩感知的要求，为了确保输入信号能够适应压缩感知，需要对输入信号进行正交变换，将输入信号转化为稀疏信号，正交变换由以下公式表示：

g＝ψs (14)

其中，大小为N×N的矩阵ψ为稀疏基，原始信号g投影到稀疏基上得到大小为N×1的稀疏向量s；

步骤7.2：求稀疏解，电容层析图像重构问题转化为L0范数优化问题，由于信号s是稀疏的，获得稀疏解

最直接的方法就是求解L0 范数优化问题，用OMP算法：

其中，L0范数||s||₀表示向量s中非零系数的数目，通过正交匹配追踪算法OMP来求解，大小为N×N的矩阵ψ为稀疏基；

步骤八，图像重建，原始图像信号被估计为：

进一步，所述步骤四中，用指数再生采样核

对狄拉克脉冲序列信号g(x)滤波后均匀采样，得到样本y_k(k＝1,2,...,K)；

M-1阶指数再生采样核的公式如下：

其中，M-1＝3是再生采样核的阶数，M-1阶再生采样核可以再生出M 个指数

C_m,k由下式得出：

其中

是

的对偶函数，属于准正交函数。

再进一步，所述步骤五中，从FRI采样样本中计算测量值，构建 FRI观测方程，指数再生采样核的阶数和FRI信号的关系为M-1＝2L-1， L为狄拉克脉冲个数，M-1阶指数再生采样核可以再生出M个复指数

其中，α_m＝a₀+jmλ，a₀为自由设定的复指数，λ为自由设定的实数。

更进一步，所述步骤7.2中，为了保证ECT系统的原始信号能够重构，需要对原始信号进行正交变换，使输入信号变得稀疏化，稀疏基为离散余弦变换基。

为了验证本发明方法的性能，进行了仿真实验验证。利用 COMSOL Multiphysics软件建立三维仿真模型并生成仿真数据m×1的电容向量C和m×n的灵敏度矩阵S。利用Matlab R2019a仿真工具对数据进行处理并重构ECT图像。三维仿真模型建立如下:管道为方形， 8个电容传感器，两相流的管道内截面剖分为20×20个像素，两相流的介电常数分别设置为3和1。指数再生采样核阶数设为3。即M-1＝3。

二种流型的原始图像如图4所示，分别为核心流和半管流，LBP 算法重建的图像效果如图5所示、Landweber算法重建的图像效果如图6所示、Tikhonov算法重建的图像效果图7所示、本发明方法重建的图像效果图8所示。从重构结果可以看出，本发明方法重建的图像与原始图像接近，图像边缘清晰平滑。本发明方法重建的图像质量明显优于传统算法。

为了对图像重建方法的性能进行数据评价，将误差和相关系数作为性能指标。误差定义为重建后的图像灰度向量与原始图像灰度向量之间的差异度，公式为:

其中,g是原始灰度向量，

是重建的灰度向量。

相关系数定义为重建后的图像灰度向量与原始图像灰度向量之间的线性相关度，可以定义为:

其中,

是向量g的平均值，而

是向量

的平均值。

表1给出了三个传统算法LBP算法、Landweber算法、Tikhonov 正则化算法和本发明方法在两种流型(核心流和半管流)下的重建图像误差和相关系数。

表1

由表1可以得出，对于两种不同的流型(核心流和半管流)，本发明方法与传统算法相比，本发明方法重建的图像误差远小于传统算法，相关系数都要高于传统算法。这意味着用本发明方法重建的图像与原始图像最相似。本发明方法获得了较好的图像，效果优于传统算法，具有良好的图像质量和较高的重建精度。

Claims (4)

1.一种多指数特征提取的ECT图像重建方法，其特征在于，所述方法包括以下步骤：

步骤一，初始化数据，获取并归一化数据，ECT图像重构的问题是通过测量电容向量C求解灰度向量g的过程，估计管道内的介电常数分布；

步骤二，补零重组，构建稀疏化的稀疏观测方程，在实际的ECT系统中，电容传感器的数量一般不会很多，所剖分的像素个数n远大于测量电容值个数m，从而导致采样率太低，降低重构的精度；通过零向量扩展的方式增加虚拟电极，以获取更多的测量电容值，构建一个稀疏化的稀疏观测方程，ECT系统的数学模型写为：

λ₀＝S₀g (4)

其中，λ₀是通过补零重组的方法扩展后的电容向量；S₀是通过补零重组的方法扩展后的灵敏度矩阵；

步骤三，有限新息率FRI信号建模，由于灰度向量取值范围为[0,1]，是一个典型的离散FRI信号，因此把Tikhonov正则化算法所得到的灰度向量g_v建模为狄拉克脉冲序列信号g(x)：

g_v＝(S^T·S+α·I)^-1·S^TC (5)

其中，g_v为大小为n×1灰度值向量，α是正则化参数，I是一个大小为n×n的单位阵，由于S^T·S存在不可逆的可能，所以加入参数α用于解决S^T·S不适定的问题，正则化算法可以有效解决灵敏度矩阵不可逆的病态问题；

其中，x＝1,2...,N是像素的位置，L是g(x)非零灰度值的个数，a_l∈[0,1]是像素x_l的灰度值，x_l∈{1,2,....,N}是非零灰度值的像素位置；

步骤四，有限新息率FRI采样，首先使用指数再生采样核

对狄拉克脉冲序列信号g(x)滤波，再对其均匀采样，可以得到K个样本值y_k(k＝1,2,...,K)；

均匀采样后得到K个样本y_k：

其中，y(x)＝g(x)*h(x)是FRI信号g(x)滤波后的结果，<·,·>表示内积，y_k是滤波后均匀采样的样本值，T是采样间隔，K是采样值的数量；

步骤五，从FRI采样样本中获取测量值，构建FRI观测方程，用样条系数C_m,k对FRI的采样值y_k进行加权后求和，可以计算出M个复指数

个测量值可以由K个样本y_k计算得到：

公式(8)写成矩阵形式：

因为像素位置x_l属于集合{1,2,...,N}，a_l∈[0,1]是像素x_l的灰度值；

所以公式(8)写成：

U＝Ag (11)其中，U＝[τ₀,τ₁,...,τ_M-1]^T∈R^M×1是FRI测量向量，g＝[g₁,g₂,...,g_N]^T是ECT系统中大小为N×1的灰度向量，A是由集合

构成M×N的FRI测量矩阵；

步骤六，构建综合观测方程，步骤如下：

步骤6.1：稀疏观测矩阵结合FRI观测矩阵，由于公式(4)和(10)都是对g的一个多维观测，所以两个观测方程可以结合，得到一个新的综合观测方程，结合公式(4)和(11)得到新的综合观测方程：

步骤6.2：以相同规则随机重组综合观测向量和综合观测矩阵的行向量，因此ECT系统的数学模型重写为：

λ_new＝S_newg (13)

其中，λ_new是综合观测向量；S_new是综合观测矩阵；

步骤七，求综合观测方程的稀疏解，步骤如下：

步骤7.1：稀疏变换，只有信号是稀疏的，才能满足压缩感知的要求，对于ECT系统，信号稀疏程度大部分情况下不能满足压缩感知的要求，为了确保输入信号能够适应压缩感知，需要对输入信号进行正交变换，将输入信号转化为稀疏信号，正交变换由以下公式表示：

g＝ψs (14)

其中，大小为N×N的矩阵ψ为稀疏基，原始信号g投影到稀疏基上得到大小为N×1的稀疏向量s；

步骤7.2：求稀疏解，电容层析图像重构问题转化为L0范数优化问题，由于信号s是稀疏的，获得稀疏解

最直接的方法就是求解L0范数优化问题，用OMP算法：

其中，L0范数||s||₀表示向量s中非零系数的数目，通过正交匹配追踪算法OMP来求解，大小为N×N的矩阵ψ为稀疏基；

步骤八，图像重建，原始图像信号被估计为：

2.如权利要求1所述的一种多指数特征提取的ECT图像重建方法，其特征在于，所述步骤四中，用指数再生采样核

对狄拉克脉冲序列信号g(x)滤波后均匀采样，得到样本y_k(k＝1,2,...,K)；

M-1阶指数再生采样核的公式如下：

其中，M-1＝3是再生采样核的阶数，M-1阶再生采样核可以再生出M个指数

C_m,k由下式得出：

其中

是

的对偶函数，属于准正交函数。

3.如权利要求1或2所述的一种多指数特征提取的ECT图像重建方法，其特征在于，所述步骤五中，从FRI采样样本中计算测量值，构建FRI观测方程，指数再生采样核的阶数和FRI信号的关系为M-1＝2L-1，L为狄拉克脉冲个数，M-1阶指数再生采样核可以再生出M个复指数

其中，α_m＝a₀+jmλ，a₀为自由设定的复指数，λ为自由设定的实数。

4.如权利要求1或2所述的一种多指数特征提取的ECT图像重建方法，其特征在于，所述步骤7.2中，为了保证ECT系统的原始信号能够重构，需要对原始信号进行正交变换，使输入信号变得稀疏化，稀疏基为离散余弦变换基。

Images