技术领域
[0001] 本
发明涉及电磁
层析成像(EMT)技术领域，具体涉及一种基于有限新息率的EMT图像重构方法。
背景技术
[0002] 电
磁层析成像(EMT)技术是一种基于
电磁感应原理的新型
过程层析成像技术，用于检测和重现可用电导率或磁导率表示的
多相流物场空间，如工业管道内的多相流检测。EMT系统通过激励线圈内的交变
电流产生激励
磁场，激励磁场与被测区域内
导电性或导
磁性物质相互作用形成物场。当被测物体的电导率或磁导率发生改变时，物场内的磁感应
密度也会随之变化，进而影响检测线圈检测到的
电压值。最后，利用检测线圈检测到的电压值重建出物场空间内电导率或磁导率的分布情况。EMT技术具有非
接触、非侵入、安全性高、成本低和实时成像等优点，目前已成为国内外的研究热点。
[0003] 典型的EMT成像系统由三部分组成：
传感器阵列(即激励和检测线圈)、
接口与控制单元和计算机。传感器阵列用于产生激励磁场，并检测被测物场内电导率或磁导率的变化，然后将检测到的电压
信号传送给接口与控制
电路；接口与控制单元一方面控制激励线圈产生磁场，另一方面实时地采集检测线圈电压信号，并完成相应的信号解调与滤波等处理，同时负责与上位机的通讯；计算机即上位机主要运行Matlab等
软件，通过一定的EMT图像重建
算法重构出物场空间的电导率和磁导率的分布图像。
[0004] EMT技术基于电磁感应原理，其正问题实际上是求解一个
角频率为ω的时谐开域
涡流场问题，可以用麦克斯韦(Maxwell)方程组来描述，正弦时变
电磁场的Maxwell方程为：
[0005]
[0006] 其中，B是磁感应强度，E是
电场强度，D是位移电流，σ是介质的电导率，μ是介质的磁导率，Js是激励线圈的电流密度。
[0007] 在忽略位移电流D的情况下，引入矢量磁位A，即 采用二维近似求解，并引入库仑规范 上述时谐电磁场方程最终变为：
[0008]
[0009] 利用有限元方法求解公式二，即可获得被测区域内的矢量磁位A。
[0010] 在进行图像重建之前，需要获取检测线圈的测量电压V和被测物场的敏感场S。其中测量电压可以通过矢量磁位A计算：
[0011] V＝jω∫A·dl
[0012] 敏感场S的计算是实现图像重建的一个重要环节，它表征着被测物场内介质电导率的微小变化与检测线圈上感生电动势变化的一种映射关系。根据互易原理，定义第k个
位置处导体的电导率发生变化时的敏感矩阵为：
[0013]
[0014] 其中， 为单元k的体积，Ai和Aj分别为当激励线圈和检测线圈通以电流I0时的正问题计算结果。
[0015] EMT逆问题的实质就是在正问题分析所得检测线圈的测量电压V和敏感场S的
基础之上，求出被测对象的电导率或磁导率分布信息，即EMT图像重建。定义敏感场的分布函数为Sed[(x,y),μ(x,y),σ(x,y)]，其物理含义是第e个投影方向下，第d个检测线圈输出的电压变化率。类似于敏感场的分布，第e个投影方向下，第d个检测线圈的测量电压值为Ved。当介质的电导率为一定值时，敏感场的分布函数Sed[(x,y),μ(x,y),σ(x,y)]和检测电压Ved具有以下关系：
[0016]
[0017] 其中Φ为被测区域。
[0018] 在EMT系统中，测量电压值和被测区域内电导率或磁导率的分布之间的关系为非线性的。当网格剖分的单元足够多时，非线性关系可以近似为线性关系，表示为：
[0019] V＝Sg
[0020] 其中，向量V∈RD×1表示检测线圈获取的测量电压值，D是投影的数目；矩阵S∈RD×NN×1表示敏感场，N为被测区域被剖分的网格单元数；向量g∈R 表示被测区域电导率或磁导率分布。EMT图像重建问题即是在已知测量电压V和敏感场S的情况下，求解被测区域电导率或磁导率分布g的问题，由于D＜＜N，该问题是病态的。
[0021] EMT图像重建问题是EMT技术的核心，其最终目标是重构出多相流二维或三维的分布信息，重建图像的
质量将直接影响到管道内电导率和磁导率分布的视觉效果。目前国内外的学者提出了许多EMT图像重建算法：线性反投影(LBP)算法是最早提出且被广泛应用的EMT图像重建算法，该算法编程特别简单，成像速度快，但成像
精度不高，适合于准确度要求不高的工业在线检测；Landweber算法的成像质量比LBP算法要好一些，然而该算法基于最速下降法的思想，沿着负梯度方向搜索最优解的，该方向并不是理想的搜索方向，容易收敛到局部极值；Tikhonov正则化算法是处理病态问题的有效方法，虽然能得到稳定的近似图像重建结果，然而该算法的解过于光滑，导致图像边缘的重建效果不是很好；Newton-Raphson算法通过求解非线性方程的方式获得最小二乘解，然而该方法的收敛速度很慢，且容易局部收敛；神经网络算法利用大量的样本进行训练，从而对网络的
阈值和权值进行调整以确立其映射关系模型，只适用于处理一些简单的流型；代数重建算法(ART)和同步
迭代重建算法(SIRT)当敏感场矩阵较大时尤其适用，但是图像重建质量并不能够满足工业检测的需求；其他的卡尔曼滤波方法和粒子滤波方法通过迭代的思想，使用一些先验信息来估计图像，然而卡尔曼滤波在迭代过程中要求矩阵的逆，粒子滤波方案则需要大量的样本，这将导致很长的运算时间。迄今为止，高精确度且高效率的EMT图像重建，仍是一个关键问题。
发明内容
[0022] 针对现有EMT图像重建问题，本发明提出一种基于有限新息率的EMT图像重建算法。首先，将EMT技术重建的图像建模为有限新息率(FRI)信号，并在FRI
采样框架下用多项式再生采样核采样，通过这种采样的方式提取出原始图像信号的一些特征信息。然而，EMT图像信号并不能够直接被采样，所以需要使用一个LBP线性反投影信号来近似。在经过FRI采样之后，根据多项式再生核的性质用采样样本重新构建测量方程，弥补EMT系统测量信息不足的
缺陷。最后，将EMT图像重建问题转换为一个最小L0范数下的优化问题，采用
正交匹配追踪(OMP)算法来求解。
[0023] 本发明解决其技术问题所采用的技术方案是：
[0024] 一种基于有限新息率的EMT图像重构方法，所述方法包括以下步骤：
[0025] 步骤一，初始化：考虑电磁层析成像EMT系统检测工业管道内多相流分布的场景，从检测线圈上获取的测量电压为V＝[V1,V2,…,VD]T，其中D为管道周围的线圈数量；管道被D剖分的网格数为N＝P×P，其中P为横轴或纵轴的网格数；管道内物场的敏感场矩阵为S∈R×N；
[0026] 步骤二，有限新息率FRI信号建模：用一个大小为N×1的向量g来表示EMT图像信号，由于该图像灰度值的取值范围为[0,1]，图像信号g建模为有限长的Diracs序列，这是一个典型的离散FRI信号：
[0027]
[0028] 其中，x＝1,2,…,N是
像素位置，L是g(x)中非零灰度值的数量，xl∈{1,2,...,N}是g(x)中非零灰度值的像素位置，l＝0,1,…,L-1，al∈[0,1]是像素位置xl处的灰度值，由于图像信号不能被直接采样，用一个LBP线性反投影信号g0＝STV来近似g；
[0029] 步骤三，有限新息率FRI采样：对于离散FRI信号g(x)，采用M-1阶的多项式再生采样核 对其进行滤波后均匀采样，得到K个采样值：
[0030]
[0031] 其中，k＝1,2,…,K是采样值的数量，T是采样间隔，K＝N/T是总的采样值数量；
[0032] 步骤四，信号观测：根据多项式再生核的性质，M-1阶再生核能够再生出M个多项式，用系数cm,k(m＝0,1,…,M-1；k＝1,2,…,K)对FRI采样值yk(k＝1,2,…,K)进行加权并求和，即得到M个测量值τm(m＝0,1,…,M-1)：
[0033]
[0034] 令um＝τm·Tm，公式(3)写成矩阵的形式：
[0035]
[0036] 步骤五，测量值的完备线性表示：由于像素的位置xl∈{1,2,...,N},测量值um表示为集合{1,2,...,N}中所有元素的完备线性组合：
[0037]
[0038] 其中，[g1,g2,…,gN]T是一个N×1的向量，由L个灰度值 和N-L个零元素组成。令U＝[u0,u1,…,uM-1]T∈RM×1，g＝[g1,g2,…,gN]T，那么公式(5)简化为：
[0039] U＝Ag (6)
[0040] 其中，A是一个M×N大小的矩阵，由集合 中的元素组成；
[0041] 步骤六，重新生成测量方程：结合EMT系统中获得的检测线圈测量电压V，得到新的测量向量为λ＝[V；U]；同理，结合EMT系统中管道内物场的敏感场矩阵S，可得到新的敏感场矩阵Φ＝[S；A]，最后，新的测量方程变为：
[0042] λ＝Φg (7)
[0043] 步骤七，求稀疏解：一般情况下，向量g是非稀疏的，采用正交基Ψ对其作
稀疏变换，即g＝Ψs，这样，对公式(7)的测量方程的求解就转换为求解一个最小L0范数下的优化问题：
[0044]
[0045] 其中，s是N×1的稀疏向量，是原始信号g在稀疏基ψ上的投影，公式(8)是一个病态问题，采用正交匹配追踪OMP算法进行求解；
[0046] 步骤八，EMT图像重建：在求得稀疏解 之后，原始图像信号能够用公式(9)来估计；
[0047]
[0048] 其中，Ψ′是正交基矩阵Ψ的逆。
[0049] 进一步，所述步骤三中，所述的多项式再生核 的阶数为M-1，能够得到M个测量值，M的取值越大，能够提取到的特征信息越多，图像重构的精度也就越高，但是算法运行时间也会随着增加。
[0050] 所述步骤四中，所述的系数cm,k的计算方式为：
[0051]
[0052] 其中，函数 为多项式再生核 的对偶函数，属于准正交函数。
[0053] 所述步骤七中，当被测管道内电导率非零的介质较少时，可以不用稀疏变换，即将正交基Ψ用单位矩阵I来代替。
[0054] 本发明的有益效果主要表现在：将EMT图像建模为FRI信号，并用多项式再生核进行采样。FRI采样的目的是为了从原始的图像信号中提取出一些特征参数作为测量值，弥补EMT图像重构问题中测量信息不足的缺陷，从而增强图像重构的效果。用正交基对图像信号做稀疏变换，从而将新的测量方程转换为求解一个最小L0范数下的优化问题，可以用OMP算法求解。仿真实验结果显示，本发明方法的图像重建精度较高，是一种高效的EMT图像重建算法。
附图说明
[0055] 图1是FRI采样结构图。
[0056] 图2是四种流型的原始图像。
[0057] 图3是LBP算法的图像重构结果。
[0058] 图4是Landweber迭代算法的图像重构结果。
[0059] 图5是TV正则化算法的图像重构结果。
[0060] 图6是本发明方法的图像重构结果。
具体实施方式
[0061] 下面结合附图对本发明作进一步描述。
[0062] 参照图1，一种基于有限新息率的EMT图像重构方法，其特征在于，所述方法包括以下步骤：
[0063] 步骤一，初始化：考虑电磁层析成像EMT系统检测工业管道内多相流分布的场景，从检测线圈上获取的测量电压为V＝[V1,V2,…,VD]T，其中D为管道周围的线圈数量；管道被剖分的网格数为N＝P×P，其中P为横轴或纵轴的网格数；管道内物场的敏感场矩阵为S∈RD×N；
[0064] 步骤二，有限新息率FRI信号建模：用一个大小为N×1的向量g来表示EMT图像信号，由于该图像灰度值的取值范围为[0,1]，图像信号g建模为有限长的Diracs序列，这是一个典型的离散FRI信号：
[0065]
[0066] 其中，x＝1,2,…,N是像素位置，L是g(x)中非零灰度值的数量，xl∈{1,2,...,N}是g(x)中非零灰度值的像素位置，l＝0,1,…,L-1，al∈[0,1]是像素位置xl处的灰度值，由于图像信号不能被直接采样，用一个LBP线性反投影信号g0＝STV来近似g；
[0067] 步骤三，有限新息率FRI采样：参照图1，对于离散FRI信号g(x)，采用M-1阶的多项式再生采样核 对其进行滤波后均匀采样，得到K个采样值：
[0068]
[0069] 其中，k＝1,2,…,K是采样值的数量，T是采样间隔，K＝N/T是总的采样值数量；
[0070] 步骤四，信号观测：根据多项式再生核的性质，M-1阶再生核能够再生出M个多项式，用系数cm,k(m＝0,1,…,M-1；k＝1,2,…,K)对FRI采样值yk(k＝1,2,…,K)进行加权并求和，即得到M个测量值τm(m＝0,1,…,M-1)：
[0071]
[0072] 令um＝τm·Tm，公式(3)写成矩阵的形式：
[0073]
[0074] 步骤五，测量值的完备线性表示：由于像素的位置xl∈{1,2,...,N},测量值um表示为集合{1,2,...,N}中所有元素的完备线性组合：
[0075]
[0076] 其中，[g1,g2,…,gN]T是一个N×1的向量，由L个灰度值 和N-L个零元素组成。令U＝[u0,u1,…,uM-1]T∈RM×1，g＝[g1,g2,…,gN]T，那么公式(5)简化为：
[0077] U＝Ag (6)
[0078] 其中，A是一个M×N大小的矩阵，由集合 中的元素组成；
[0079] 步骤六，重新生成测量方程：结合EMT系统中获得的检测线圈测量电压V，得到新的测量向量为λ＝[V；U]；同理，结合EMT系统中管道内物场的敏感场矩阵S，可得到新的敏感场矩阵Φ＝[S；A]，最后，新的测量方程变为：
[0080] λ＝Φg (7)
[0081] 步骤七，求稀疏解：一般情况下，向量g是非稀疏的，采用正交基Ψ对其作稀疏变换，即g＝Ψs，这样，对公式(7)的测量方程的求解就转换为求解一个最小L0范数下的优化问题：
[0082]
[0083] 其中，s是N×1的稀疏向量，是原始信号g在稀疏基ψ上的投影，公式(8)是一个病态问题，采用正交匹配追踪OMP算法进行求解；
[0084] 步骤八，EMT图像重建：在求得稀疏解 之后，原始图像信号能够用公式(9)来估计；
[0085]
[0086] 其中，Ψ′是正交基矩阵Ψ的逆。
[0087] 进一步，所述步骤三中，所述的多项式再生核 的阶数为M-1，能够得到M个测量值。M的取值越大，能够提取到的特征信息越多，图像重构的精度也就越高，但是算法运行时间也会随着增加。
[0088] 所述步骤四中，所述的系数cm,k的计算方式为：
[0089]
[0090] 其中，函数 为多项式再生核 的对偶函数，属于准正交函数。
[0091] 所述步骤七中，当被测管道内电导率非零的介质较少时，可以不用稀疏变换，即将正交基Ψ用单位矩阵I来代替。
[0092] 实验对比：为了验证本发明方法的性能，进行了仿真实验。使用COMSOL Multiphysics软件建立一个3维的仿真模型并生成仿真数据；使用MATLAB R2013a仿真工具处理数据和重构EMT图像。3维仿真模型的设置如下：被测管道具有8个检测线圈的方形管，可以检测到的测量电压数量为 管道横截面被分成40×40个像素元，即N＝1600；管道半径为0.05m，屏蔽层半径为0.07m，检测线圈半径为0.01m；管道内的被测介质是电导率为σ＝3.774×107S/m的
铝和电导率为σ＝0的空气。FRI的多项式再生采样核的阶数是7，即M＝8。
[0093] 四种流型的原始图像如图2所示，不同EMT重构算法的仿真结果如图3-6所示。将这些重建图像进行比较，可以看出本发明方法重构的图像非常接近原始图像，并且图像的边缘清晰且光滑，散粒噪声较少。本发明方法重构图像的质量明显优于用LBP算法、Landweber迭代算法和TV正则化算法。
[0094] 为了定量地评价本发明方法的性能，采用图像误差和相关系数作为评价指标。图像误差是指重构图像和原始图像之间的差异程度，被定义为：
[0095]
[0096] 其中g是原始图像的灰度向量， 是重构图像的灰度向量。相关系数指的是重构图像和原始图像之间的线性相关程度，定义为：
[0097]
[0098] 其中 是向量g的平均值； 是向量 的平均值。
[0099] 表1和表2分别展示了使用不同算法重构的图像误差和相关系数。从表1的图像误差数据对比可以看出，对于四种不同的流型，本发明方法重构的图像误差要小于LBP算法、Landweber迭代算法和TV正则化算法。类似的从表2的相关系数数据对比可以看出，本发明方法重构的图像相关系数要大于其他算法。综上所述，本发明方法是一种有效的EMT图像重构方法，而且重建图像的精度较高。表1是不同算法重建的图像误差
[0100]
[0101] 表1
[0102] 表2是不同算法重建的图像相关系数
[0103]
[0104] 表2。
