基于压缩感知流型辨识的变化灵敏度场电容层析成像方法
1.本发明属于图像处理技术领域，具体涉及一种基于压缩感知流型辨识的变化灵敏度场电容层析成像方法。
背景技术：
2.在工业和生活领域中，两相流现象是普遍存在的，故过程层析成像技术在工业测量领域发展前景非常广阔。两相流系统如今广泛应用于工业中以适应现代化工业的需要。电容层析成像技术(electrical capacitance tomography，ect)是一种成熟的过程层析成像技术。ect技术是利用被测管道周围的传感器阵列获取被测物质的投影数据(电容值)，并用图像重建算法重建出被测区域的介电常数分布。ect技术相对于x射线层析成像、伽马射线层析成像、超声波层析成像等其他过程层析成像技术，具有硬件成本低、响应速度快、硬件结构简单、安全性能好(非辐射)等优点，其应用前景十分广阔。
3.ect技术的基本思想是：不同的流动物质具有不同的相对介电常数，当各流相物流动不均或浓度发生改变时，都将致使流体相对介电常数的变化，从而引起检测极板之间电容响应的改变。检测中采用阵列式的电容传感器，各电容传感器之间彼此组合可以产生多个独立测量值。最后通过成像计算机即可重建出两相流体的介电常数分布情况。典型的ect系统结构图如图1所示。
4.被测管道内部的流体介电常数分布ε(x,y)与从传感器阵列获取的一组电容向量ci之间的关系可以描述为：
5.ci＝∫∫dsi(x,y)ε(x,y)dxdy
ꢀꢀꢀ
(1)
6.其中，d是被测管道区域，(x,y)∈d表示像素点的位置；ε(x,y)为介电常数分布函数，ci(i＝1,2...m)为传感器阵列测量所得的m个电容测量值。si(x,y)是灵敏度的分布函数。为了快速计算灵敏度场，电势分布法是一个很好的方法。根据电势分布法，灵敏度场可以描述为以下公式：
[0007][0008]
其中v是激励电压，和分别是电极i和电极j为激励电极时的电势分布。
[0009]
ect的图像重构是利用测量获得的m个独立电容值来反演管道内介电常数分布。目前已经有多种图像重建方法，但是由于ect灵敏度的“软场”问题，并且观测数据极为有限，是一个欠定问题，因此很难重建出高质量的图像，所以图像重建算法一直是人们研究的热点。目前在ect系统中有两种典型的图像重建算法：线性反投影法(lbp)和landweber迭代法。lbp算法重建的图像质量相对较低，但lbp算法的优点是速度快；landweber迭代法算法是一种成像精度相对较高的迭代算法，它使用最速下降法来进行收敛迭代，具有更好的图像质量。这两种传统算法都是利用固定灵敏度场来求解正问题的，固定灵敏度场是假定被测电容是与整个管道内部变化没有关系的，这显然是不符合实际的假设。因此，为了减少
ect灵敏度场的“软场”影响，这里提出一种基于压缩感知流型辨识的变化灵敏度场电容层析成像方法。迄今为止，提高电容层析图像重建的质量和精度仍然是个关键技术问题。
技术实现要素：
[0010]
为了克服已有技术的不足，针对电容层析成像问题，提出一种基于压缩感知流型辨识的变化灵敏度场电容层析成像方法。首先本发明方法设定k类典型流型，将其归一化后的ect测量电容信号作为标准样本，用过调整介质分布和叠加随机噪声的方式来生成训练样本集。其次将随机高斯矩阵作为测量矩阵，分别求得测试样本和标准样本的观测值，并用压缩感知信号重构算法如正交匹配追踪算法(omp)分别求出测试样本和标准样本在过完备训练样本集e投影下的稀疏系数。最后计算测试样本与所有标准样本的稀疏解之间的线性相关系数，并选择最大相关系数所对应的流型为测试样本的归属流型，并用该归属流型的灵敏度场替代原始的固定灵敏度场进行图像重构；提高了电容层析成像的精度。
[0011]
本发明解决其技术问题所采用的技术方案是：
[0012]
一种基于压缩感知流型辨识的变化灵敏度场电容层析成像方法，包括以下步骤：
[0013]
步骤一，初始化，归一化数据，被测管道内部的流体介电常数分布ε(x,y)与从传感器阵列获取的一组电容向量ci之间的关系可以描述为：
[0014]ci
＝∫∫dsi(x,y)ε(x,y)dxdy
ꢀꢀꢀ
(1)
[0015]
其中，d是被测管道区域，(x,y)∈d表示像素点的位置；ε(x,y)为介电常数分布函数，ci(i＝1,2...m)为传感器阵列测量所得的m个电容测量值，si(x,y)是灵敏度的分布函数；
[0016]
根据电势分布法，灵敏度场描述为以下公式：
[0017][0018]
其中，v是激励电压，和分别是电极i和电极j为激励电极时的电势分布；
[0019]
由公式(1)看出，当介电常数相近时，电容向量ci与介电常数分布函数ε(x,y)之间的关系可近似为线性；为了方便计算，对其进行离散和归一化处理，得到ect矩阵形式的数学模型如下所示：
[0020]
c＝sg
ꢀꢀꢀ
(3)
[0021]
其中，c是归一化的m
×
1的电容测量值向量，g是大小为n
×
1的灰度向量，代表管道不同介电常数的分布；s是归一化的大小为m
×
n的灵敏度矩阵，m为测量电容值的个数。n的大小代表了成像区域像素个数，从公式(3)可以看出，ect图像重建的问题可以简化为通过测量电容向量c求解灰度向量g的过程；
[0022]
步骤二，生成训练样本集。给定k类流型的测量电容数据组成训练样本集，记第i类流型的训练集ei共有ni个样本，训练集ei表示为：
[0023][0024]
其中，vj∈rm×1(j＝1,2
…
ni)代表第i类流型里面的第j个训练样本，即ect系统从第i类流型获取的测量电容向量；m是测量电容值得个数，8极板ect系统有m＝28，通过调整介质分布和叠加随机噪声的方式来生成训练样本集，每类流型的训练样本数为100个，则所
有k类流型组成的训练样本集为：
[0025]
e＝[e1,e2…ek
]∈rm×n，
ꢀꢀꢀ
(5)
[0026]
其中n＝n1+n2+
…
nk；
[0027]
对于流型i，若训练集ei完备或者过完备，则此类流型的任意测量电容向量λ可由其训练集ei线性表示：
[0028][0029]
其中，系数αi(j＝1,2
…
ni)为实数，为系数向量，在实际应用中，由于无法预知任意测量电容向量λ的归属类，可以考虑将其重新表示为k类流型组成的训练样本集e的线性组合：
[0030]
λ＝ex0.
ꢀꢀꢀ
(7)
[0031]
式中，为稀疏向量，至此，由ect系统获得的测量电容信号λ就可由以训练样本集组成的过完备字典变换基e来稀疏表示，这样就使得ect测量电容信号满足了压缩感知理论的前提条件；
[0032]
步骤三，计算测试样本和标准样本的观测值，随机高斯矩阵b∈rd×m作为测量矩阵，则根据压缩感知思想，针对以测量电容向量λ作为原始信号的线性测量过程为：
[0033][0034]
其中，为观测值，是d
×
1的向量，标准样本的观测值按照同样方法计算；
[0035]
步骤四，计算测试样本和标准样本的稀疏系数，由于观测向量的维数d远小于稀疏信号x0的维数，直接求解公式(8)是一个欠定问题，将其转换为如公式(9)所示的数学模型：
[0036][0037]
通过压缩感知的信号重构算法求解此l0范数下的最优化问题，使用正交匹配追踪(omp)算法得到测量电容向量λ和标准电容向量λk对应于过完备字典变换基e的稀疏解x0和xk；
[0038]
步骤五，计算测试样本和标准样本稀疏解的线性相关系数，通过x0与xi的线性相关程度来确定测试样本的归属流型，达到流型辨识的目的；
[0039]
相关系数则是衡量变量之间线性相关程度的指标，其计算方法如式(10)所示：
[0040][0041]
其中，和分别为x0和xi的平均值。根据公式(10)分别计算测试样本与所有标准样本对应的稀疏解之间的线性相关系数；
[0042]
步骤六，确定流型，选择相关系数最大的归属流型灵敏度场，过程如下：将最大相关系数所对应的流型确定为测试样本的归属流型。选择该归属流型的灵敏度场确定为新的灵敏度矩阵；
[0043]
步骤七，图像重构，输出灰度向量g。
[0044]
进一步，所述步骤七中，用变化灵敏度矩阵和重构算法来进行图像重构。
[0045]
再进一步，所述步骤三中，使用随机高斯矩阵b∈rd×m作为测量矩阵，为了准确的重构出原始信号，测量矩阵必须满足rip准则，rip准则的等价条件是测量矩阵和稀疏基不相关，当以d
×
m维的随机高斯矩阵作为测量矩阵时，只要满足条件：d≥cklog(m/k)，其中c是一个很小的数，k为原始信号的稀疏度，此测量矩阵在极大概率下满足rip准则。测量矩阵可以替换，例如离散fourier矩阵和伯努利随机矩阵。
[0046]
更进一步，所述步骤六中，为了确定流型，选择相关系数最大的归属流型灵敏度场。可以用其他性能指标来选择流型，例如最小相对误差和最小相对电容残差等。
[0047]
本发明的有益效果主要表现在：提高了电容层析成像的精度。
附图说明
[0048]
图1是ect系统结构图。
[0049]
图2是基于压缩感知流型辨识的变化灵敏度场电容层析成像方法流程图。
[0050]
图3是七种流型的两种算法成像效果图。
具体实施方式
[0051]
下面结合附图对本发明作进一步描述。
[0052]
参照图2和图3，一种基于压缩感知流型辨识的变化灵敏度场电容层析成像方法，包括以下步骤：
[0053]
步骤一，初始化，归一化数据，被测管道内部的流体介电常数分布ε(x,y)与从传感器阵列获取的一组电容向量ci之间的关系可以描述为：
[0054]ci
＝∫∫dsi(x,y)ε(x,y)dxdy
ꢀꢀꢀ
(1)
[0055]
其中，d是被测管道区域，(x,y)∈d表示像素点的位置；ε(x,y)为介电常数分布函数，ci(i＝1,2...m)为传感器阵列测量所得的m个电容测量值，si(x,y)是灵敏度的分布函数；
[0056]
根据电势分布法，灵敏度场描述为以下公式：
[0057][0058]
其中，v是激励电压，和分别是电极i和电极j为激励电极时的电势分布；
[0059]
由公式(1)看出，当介电常数相近时，电容向量ci与介电常数分布函数ε(x,y)之间的关系可近似为线性；为了方便计算，对其进行离散和归一化处理，得到ect矩阵形式的数学模型如下所示：
[0060]
c＝sg
ꢀꢀꢀ
(3)
[0061]
其中，c是归一化的m
×
1的电容测量值向量，g是大小为n
×
1的灰度向量，代表管道不同介电常数的分布；s是归一化的大小为m
×
n的灵敏度矩阵，m为测量电容值的个数。n的大小代表了成像区域像素个数，从公式(3)可以看出，ect图像重建的问题可以简化为通过测量电容向量c求解灰度向量g的过程；
[0062]
步骤二，生成训练样本集。给定k类流型的测量电容数据组成训练样本集，记第i类
流型的训练集ei共有ni个样本，训练集ei表示为：
[0063][0064]
其中，vj∈rm×1(j＝1,2
…
ni)代表第i类流型里面的第j个训练样本，即ect系统从第i类流型获取的测量电容向量；m是测量电容值得个数，8极板ect系统有m＝28，通过调整介质分布和叠加随机噪声的方式来生成训练样本集，每类流型的训练样本数为100个，则所有k类流型组成的训练样本集为：
[0065]
e＝[e1,e2…ek
]∈rm×n，
ꢀꢀꢀ
(5)
[0066]
其中n＝n1+n2+
…
nk；
[0067]
对于流型i，若训练集ei完备或者过完备，则此类流型的任意测量电容向量λ可由其训练集ei线性表示：
[0068][0069]
其中，系数αi(j＝1,2ani)为实数，为系数向量，在实际应用中，由于无法预知任意测量电容向量λ的归属类，可以考虑将其重新表示为k类流型组成的训练样本集e的线性组合：
[0070]
λ＝ex0.
ꢀꢀꢀ
(7)
[0071]
式中，为稀疏向量，至此，由ect系统获得的测量电容信号λ就可由以训练样本集组成的过完备字典变换基e来稀疏表示，这样就使得ect测量电容信号满足了压缩感知理论的前提条件；
[0072]
步骤三，计算测试样本和标准样本的观测值，随机高斯矩阵b∈rd×m作为测量矩阵，则根据压缩感知思想，针对以测量电容向量λ作为原始信号的线性测量过程为：
[0073][0074]
其中，为观测值，是d
×
1的向量。标准样本的观测值按照同样方法计算；
[0075]
步骤四，计算测试样本和标准样本的稀疏系数，由于观测向量的维数d远小于稀疏信号x0的维数，直接求解公式(8)是一个欠定问题，将其转换为如公式(9)所示的数学模型：
[0076][0077]
通过压缩感知的信号重构算法求解此l0范数下的最优化问题，使用正交匹配追踪(omp)算法得到测量电容向量λ和标准电容向量λk对应于过完备字典变换基e的稀疏解x0和xk；
[0078]
步骤五，计算测试样本和标准样本稀疏解的线性相关系数，通过x0与xi的线性相关程度来确定测试样本的归属流型，达到流型辨识的目的；
[0079]
相关系数则是衡量变量之间线性相关程度的指标，其计算方法如式(10)所示：
[0080][0081]
其中，和分别为x0和xi的平均值。根据公式(10)分别计算测试样本与所有标准
样本对应的稀疏解之间的线性相关系数；
[0082]
步骤六，确定流型，选择相关系数最大的归属流型灵敏度场，过程如下：将最大相关系数所对应的流型确定为测试样本的归属流型。选择该归属流型的灵敏度场确定为新的灵敏度矩阵；
[0083]
步骤七，图像重构，输出灰度向量g。用变化灵敏度矩阵和重构算法来进行图像重构。
[0084]
为了验证本发明方法的性能，进行了仿真模拟实验验证。利用comsol multiphysics软件建立不同流型仿真模型并生成仿真数据m
×
1的电容向量c和m
×
n的灵敏度矩阵s并导出数据，8种不同流型包括1/3层流、1/2层流、2/3层流、满管、环流、核心流、泡状流、空管。通过调整介质分布和叠加随机噪声的方式来生成训练样本集，每类流型的训练样本数为100个。用压缩感知求稀疏解的方法，求出测试样本和标准样本的稀疏系数，最后通过计算测试样本和不同标准样本之间的线性相关系数，出最大的相关系数的归属流型，用该流型的灵敏度场替换原来的固定灵敏度场。用matlab r2019a仿真工具对数据进行处理，用图像重构算法重建ect的图像。三维仿真模型建立如下：管道为圆形，8个测量极板组成阵列，两相流管道的截面被分割为40
×
40个像素，较密介质的介电常数为3，较稀介质的介电常数为1。成像算法选择lbp算法和landweber迭代算法，信噪比为50db。分别使用变化灵敏度场和固定灵敏度场进行试验。
[0085]
测量的8种流型分别为1/2层流、1/3层流、2/3层流、核心流、环流、泡状流、满管流和空管流，其中空管流的变化灵敏度场就是固定灵敏度场。两种算法分别使用变化灵敏度场和固定灵敏度场的成像结果如图x所示。从图像上可以看出，相比于使用固定灵敏度场成像，使用本发明方法重建的图像会有一个清晰的轮廓边缘，例如，层流和环状流。对于泡状流来说，使用固定灵敏度场不能很好的分辨出泡状，但本发明方法中lbp算法可以呈现出泡状的轮廓，landweber迭代算法可以清晰的显示出泡状。故本发明方法重建的图像从视觉效果上优于传统的固定灵敏度场成像方法。
[0086]
为了对本发明方法图像重建的性能进行数值上的评价，分析重建的图像与原图像的图像相对误差和相关系数。图像相对误差可定义为重建的图像灰度向量与原始图像灰度向量的差异程度，图像相对误差的公式如下:
[0087][0088]
其中，g是原始灰度向量。表示g的估计，表示为重建后的灰度向量。
[0089]
图像的相关系数可定义为重建后的灰度向量与原始灰度向量之间的线性相关度，图像的相关系数的公式如下:
[0090][0091]
其中，是向量g的平均值，而是向量的平均值。
[0092]
表1给出了本发明方法和传统的固定灵敏度场成像方法在七种流型下的图像相对误差和相关系数。
[0093][0094][0095]
表1
[0096]
由图3可得，本发明方法重建的图像误差比传统固定灵敏度场成像方法小。这说明本发明方法重建的图像好于传统固定灵敏度场成像方法。本发明方法效果良好，获得了较好的图像。
[0097]
本说明书的实施例所述的内容仅仅是对发明构思的实现形式的列举，仅作说明用途。本发明的保护范围不应当被视为仅限于本实施例所陈述的具体形式，本发明的保护范围也及于本领域的普通技术人员根据本发明构思所能想到的等同技术手段。
技术特征：
1.一种基于压缩感知流型辨识的变化灵敏度场电容层析成像方法，其特征在于，所述方法包括以下步骤：步骤一，初始化，归一化数据，被测管道内部的流体介电常数分布ε(x,y)与从传感器阵列获取的一组电容向量c
i
之间的关系可以描述为：c
i
＝∫∫
d
s
i
(x,y)ε(x,y)dxdy
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)其中，d是被测管道区域，(x,y)∈d表示像素点的位置；ε(x,y)为介电常数分布函数，c
i
(i＝1,2...m)为传感器阵列测量所得的m个电容测量值，s
i
(x,y)是灵敏度的分布函数；根据电势分布法，灵敏度场描述为以下公式：其中，v是激励电压，和分别是电极i和电极j为激励电极时的电势分布；由公式(1)看出，当介电常数相近时，电容向量c
i
与介电常数分布函数ε(x,y)之间的关系可近似为线性；为了方便计算，对其进行离散和归一化处理，得到ect矩阵形式的数学模型如下所示：c＝sg
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)其中，c是归一化的m
×
1的电容测量值向量，g是大小为n
×
1的灰度向量，代表管道不同介电常数的分布；s是归一化的大小为m
×
n的灵敏度矩阵，m为测量电容值的个数，n的大小代表了成像区域像素个数，从公式(3)可以看出，ect图像重建的问题可以简化为通过测量电容向量c求解灰度向量g的过程；步骤二，生成训练样本集，给定k类流型的测量电容数据组成训练样本集，记第i类流型的训练集e
i
共有n
i
个样本，训练集e
i
表示为：其中，v
j
∈r
m
×1(j＝1,2
…
n
i
)代表第i类流型里面的第j个训练样本，即ect系统从第i类流型获取的测量电容向量；m是测量电容值得个数，8极板ect系统有m＝28，通过调整介质分布和叠加随机噪声的方式来生成训练样本集，每类流型的训练样本数为100个，则所有k类流型组成的训练样本集为：e＝[e1,e2…
e
k
]∈r
m
×
n
，
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)其中n＝n1+n2+
…
n
k
；对于流型i，若训练集e
i
完备或者过完备，则此类流型的任意测量电容向量λ可由其训练集e
i
线性表示：其中，系数α
i
(j＝1,2
…
n
i
)为实数，为系数向量，在实际应用中，由于无法预知任意测量电容向量λ的归属类，可以考虑将其重新表示为k类流型组成的训练样本集e的线性组合：λ＝ex0.
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)式中，为稀疏向量，至此，由ect系统获得的测量电容
信号λ就可由以训练样本集组成的过完备字典变换基e来稀疏表示，这样就使得ect测量电容信号满足了压缩感知理论的前提条件；步骤三，计算测试样本和标准样本的观测值，随机高斯矩阵b∈r
d
×
m
作为测量矩阵，则根据压缩感知思想，针对以测量电容向量λ作为原始信号的线性测量过程为：其中，为观测值，是d
×
1的向量，标准样本的观测值按照同样方法计算；步骤四，计算测试样本和标准样本的稀疏系数，由于观测向量的维数d远小于稀疏信号x0的维数，直接求解公式(8)是一个欠定问题，将其转换为如公式(9)所示的数学模型：通过压缩感知的信号重构算法求解此l0范数下的最优化问题，使用正交匹配追踪(omp)算法得到测量电容向量λ和标准电容向量λ
k
对应于过完备字典变换基e的稀疏解x0和x
k
；步骤五，计算测试样本和标准样本稀疏解的线性相关系数，通过x0与x
i
的线性相关程度来确定测试样本的归属流型，达到流型辨识的目的；相关系数则是衡量变量之间线性相关程度的指标，其计算方法如式(10)所示：其中，和分别为x0和x
i
的平均值，根据公式(10)分别计算测试样本与所有标准样本对应的稀疏解之间的线性相关系数；步骤六，确定流型，选择相关系数最大的归属流型灵敏度场，过程如下：将最大相关系数所对应的流型确定为测试样本的归属流型，选择该归属流型的灵敏度场确定为新的灵敏度矩阵；步骤七，图像重构，输出灰度向量g。2.如权利要求1所述的基于压缩感知流型辨识的变化灵敏度场电容层析成像方法，其特征在于，所述步骤七中，用变化灵敏度矩阵和重构算法来进行图像重构。3.如权利要求1或2所述的基于压缩感知流型辨识的变化灵敏度场电容层析成像方法，其特征在于，所述步骤三中，使用随机高斯矩阵b∈r
d
×
m
作为测量矩阵，为了准确的重构出原始信号，测量矩阵必须满足rip准则，rip准则的等价条件是测量矩阵和稀疏基不相关，当以d
×
m维的随机高斯矩阵作为测量矩阵时，只要满足条件：d≥ck log(m/k)，其中c是一个很小的数，k为原始信号的稀疏度，此测量矩阵在极大概率下满足rip准则。
技术总结
一种基于压缩感知流型辨识的变化灵敏度场电容层析成像方法，首先，本发明方法设定k类典型流型，将其归一化后的ECT测量电容信号作为标准样本并在此基础上以调整介质分布和叠加随机噪声来生成训练样本集，待测流型电容值的归一化数据作为测试样本。然后，通过计算测试样本与所有标准样本的稀疏解之间的线性相关系数，并将最大相关系数所对应的流型确定为测试样本的归属流型，选择该归属流型的灵敏度场替代原始的固定灵敏度场进行图像重构。通过仿真模拟实验表明，所得的性能指标和成像效果优于传统的固定灵敏度场成像方法，图像边缘清晰平滑。故本发明方法提高了电容层析成像的精度，为电容层析成像研究提供了一种新的方法和手段。手段。手段。
