CN112712573B - 一种基于频谱特性的电磁层析成像图像重建方法 - Google Patents

Abstract

一种基于频谱特性的电磁层析成像图像重建方法，包括以下步骤：步骤一，电磁层析成像系统原始信号的有限新息率信号建模；步骤二，Sinc采样核设计和均匀采样；步骤三，电磁层析成像系统一维图像信号的有限新息率信号采样；步骤四，电磁层析成像系统图像信号的重构。本发明利用有限新息率信号模型将管道内目标物的电导率分布建模为有限长Diracs流，然后用Sinc采样核进行均匀采样，获得样本，提取样本的频谱信息建立新的测量方程从而利用稀疏解决方案进行图像重建。理论分析及实验结果表明，该方法在图像误差和相关系数等指标上优于原有方法，为电磁层析成像研究提供了一种新的方法和手段。

Description

一种基于频谱特性的电磁层析成像图像重建方法

技术领域

本发明涉及图像处理技术领域，具体涉及一种基于频谱特性的电磁层析成像图像重建方法。

背景技术

图像重建是通过测量物体外部数据再经数字处理获得三维被测物体的形状信息的技术，最早应用于医学领域，近年来已逐渐扩展到工业应用领域。电学层析成像是通过对被测目标施加电激励检测其电学参数的变化情况以根据所得信息逆推被测目标分布情况的工业过程层析成像技术。电磁层析成像即是基于图像重建原理和电磁学定律进行参数检测分析目标情况的一种电学层析成像技术。由于电磁层析成像兼具无需接触，检测速度快，成像效果较好等优点，因此被广泛应用于工业生产检测环节及医学检测过程中。

电磁层析成像系统中，励磁线圈通过发射励磁信号穿透被测物体产生时变磁场，时变磁场产生的涡流场再产生二次磁场，通过检测线圈将检测信号传输至电子控制单元，经电子控制单元进行信号处理后发送到计算机，由计算机进行数据处理及图像重建工作。

V＝Sg

输入信号g经过合适的采样和建模后可基于其模型结构通过检测单元检测到的电压信号V推演出被测物场的灵敏度分布情况S，从而进一步还原被测目标的电导率分布情况。作为电磁层析成像中最为复杂及关键的部分，如何快速准确重建被研究目标的复杂电导率分布是研究的主要问题。

近年来，已有多种基于电磁层析成像的重建算法被提出，诸如线性反投影(Linearback-projection,LBP)算法、Tikhonov正则化算法、Landweber迭代算法等。其中LBP图像重建算法提出时间最早，由于其简单快速的优点应用也最广泛，但其局限性在于对复杂流型的重建图像质量较低。而Tikhonov正则化算法的重建图像质量优于LBP算法，能有效应用于病态问题，但在边缘效果方面不如人意。与前两者相比，Landweber迭代算法适用性更加广泛，但该算法基于梯度下降法，所以在参数不合适的情况下结果会有所不同。因此，由于电磁层析成像系统的病态性及现有算法的局限性，如何高质量地重建图像仍是目前研究的重难点问题。

发明内容

为了克服已有电磁层析成像图像重建方法的不足，本发明提供了一种基于频谱特性的电磁层析成像图像重建方法，该方法能够高质量、快速地重建图像，首先将有限新息率信号模型应用在电磁层析成像系统中，其中两相流参数的图像信号可以看作有限长Diracs流，Diracs的时延和振幅分别表征了密相介质的像素位置和电导率；然后对有限新息率采样核进行设计，利用低通滤波方法对输入信号进行均匀采样，得到信号的傅里叶级数系数，这些系数经由零化滤波器处理获得未知的时间延迟和振幅，从而完成频谱特征信息提取工作；最后提出信号联合重建算法，利用信号的频谱特性重建原始图像。

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

一种基于频谱特性的电磁层析成像图像重建方法，包括以下步骤：

步骤一，电磁层析成像系统原始信号的有限新息率信号建模：对于两相流参数检测问题，电磁层析成像系统的原始图像为二值图像，灰度值为0或1。假设管道内的截面被分割成d×d像素，其中d为整数，其取值由网格剖分数决定；二值图像表示为二维d×d矩阵或一维N×1向量g，其中N＝d×d。由于EMT系统结构的限制，原始信号不能直接采样或测量，先将原始信号建模为可采样的具有有限新息率的信号，有限新息率即指单位时间内具有有限个参数个数，具有有限新息率的信号可在单位时间内由有限个参数确定，将信号建模为有限新息率信号使其能被均匀采样从而提取出输入信号的特征参数信息；过程为：

步骤1.1：对离散图像信号向量进行有限新息率信号建模为有限长的Diracs流：

其中T是一个周期长度，默认为1秒，便于分析处理；L₀是g₀(t)的非零灰度值个数，t_l∈{0,1,...,N-1}·(T/N)是密相介质的像素位置，n＝0,1,…,N-1是像素位置；

步骤1.2：原始信号由于其局限性，通过如LBP算法、Landweber迭代算法和TV正则化算法重建的图像来代替。这些重建的信号是灰度图像，灰度值范围从0到1，所以一个更普遍适用的电磁层析成像信号模型为：

其中L为g(n)的非零灰度值个数，t_l∈{0,1,...,N-1}·(T/N)为密相介质的像素位置，n＝0,1,…,N-1为像素位置，a_l∈[0,1]是像素t_l的灰度值，n(t)是噪声，通过公式(2)得到的g(t)为具有有限新息率的信号，其新息率为ρ＝2L/T；

步骤二，Sinc采样核设计和均匀采样：对有限长度的信号g(t)进行周期延拓使其适用于Sinc采样核，假设周期为T的信号是由Diracs产生的，这种信号通过傅立叶级数自然表示为：

Sinc采样核设计为：

其中，选择B大于或等于ρ＝2L/T的新息率进行均匀采样。通过Sinc采样核可以获取被采样信号的基带频谱信息从而为后续信号重建工作做准备；

步骤三，电磁层析成像系统一维图像信号的有限新息率信号采样过程：在对每个输入信号g(t)用Sinc采样核h(t)进行滤波之后，以t_s为采样间隔进行均匀采样，得到样本y_m，然后按照如下公式从样本y_m中获得傅立叶级数系数G[k]：

其中k＝-K,…,K，K＝BT/2；m＝0,1,…,M-1，M为样本总数，t_s＝T/M是均匀采样的时间间隔，原始信号经过采样得到样本y_m，y_m通过离散傅立叶变换提取傅立叶级数系数G[k]得到信息；

步骤四，电磁层析成像系统图像信号的重构过程：通过对原始信号建模及采样获取频谱信息后对所得信息进行分析得到新的测量方程并通过求解稀疏解得到重建后的图像信号。

步骤4.1，根据采样所得频谱信息建立新测量方程：采样阶段之后，从样本y_m中得到输入信号的傅立叶级数系数G[k]，按单周期[0,T)计算，G[k]表示为：

由于实际中噪声的存在，通过对模拟时间轴以Δ＝T/N的时间间隔进行量化，将模拟时间t近似为t＝nΔ，其中n＝0,1,…,N-1，未知参数的时间延迟t_l近似为t_l＝n_lΔ，则公式(6)近似为：

其中，n_l∈{0，1,…,N-1}为时间延迟t_l的离散数值，则测量值看作：

步骤4.2，根据信号的稀疏性建立新的测量方程：在一个周期[0,T)内，忽略量化噪声的情况下，得到一个完整的模拟时间集U＝{0,Δ,2Δ,…,(N-1)Δ}，则时间延迟参数集V＝{n₀Δ,n₁Δ,…,n_L-1Δ}构成了一个更小的子集U，即其中L＜＜N。最终公式(8)改写为：

其中[g₀,g₁,…,g_N-1]^T是一个N×1矢量，由L振幅参数和N-L零值组成。我们的目标是从测量值中找到向量[g₀,g₁,…,g_N-1]^T的非零部分，为了简单起见，公式(9)改写为：

U＝Cg (10)

其中U＝[G₁[-K],…,G₁[K]]^T∈R^Q×1是测量向量，Q＝2K+1；C是由集合组成的Q×N矩阵，g＝[g₁,g₂,…,g_N]^T是电磁层析成像的一个N×1灰度向量。

步骤4.3，给定公式(11)：

V＝Sg (11)

将两个测量方程进行重组：结合公式(10)和(11)，电磁层析成像系统的测量方程重新描述为：

设新的测量向量为λ＝[V；U]，新的灵敏度矩阵为Φ＝[S；C]，公式(5)可以简化为：

λ＝Φg (13)

其中，λ是(D+Q)×1的向量，Φ是(D+Q)×N的矩阵；

步骤4.4，针对测量方程解的稀疏性选择OMP算法进行稀疏解求解过程得到重构后的图像信号：正交变换可以将信号g转化为稀疏信号，正交变换过程描述为：

g＝Ψs (14)

结合公式(13)，公式(14)改写为：

λ＝ΦΨs (15)

由于信号是稀疏的，求解稀疏解的最直接的方法就是求解L0范数最佳化问题：

其中L0范数||s||₀表示向量s中非零系数的个数，再利用正交匹配追踪算法(OMP)迭代寻找λ和ΦΨ的列之间的最大相关性来找到s的非零部分以求解L0范数最优化问题；

步骤4.4，利用稀疏解重构原始图像：在求解稀疏信号s之后，原始图像信号可以估计为：

其中Ψ′是矩阵Ψ′的逆。

本发明的有益效果主要表现在：利用有限新息率信号模型将管道内目标物的电导率分布建模为有限长Diracs流，然后用Sinc采样核进行均匀采样，获得样本，提取样本的频谱信息建立新的测量方程从而利用稀疏解决方案进行图像重建。理论分析及实验结果表明，该方法在图像误差和相关系数等指标上优于原有方法，为电磁层析成像研究提供了一种新的方法和手段。

附图说明

图1是基于频谱特性的电磁层析成像图像重建方法的系统结构图。

图2是EMT有限新息率信号采样结构。

具体实施方式

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

参照图1和图2，一种基于频谱特性的电磁层析成像图像重建方法，包括以下步骤：

步骤一，电磁层析成像系统原始信号的有限新息率信号建模：对于两相流参数检测问题，电磁层析成像系统的原始图像为二值图像，灰度值为0或1。假设管道内的截面被分割成d×d像素，其中d为整数，其取值由网格剖分数决定；二值图像表示为二维d×d矩阵或一维N×1向量g，其中N＝d×d。由于EMT系统结构的限制，原始信号不能直接采样或测量，先将原始信号建模为可采样的具有有限新息率的信号，有限新息率即指单位时间内具有有限个参数个数，具有有限新息率的信号可在单位时间内由有限个参数确定，将信号建模为有限新息率信号使其能被均匀采样从而提取出输入信号的特征参数信息；过程为：

步骤1.1：对离散图像信号向量进行有限新息率信号建模为有限长的Diracs流：

其中T是一个周期长度，默认为1秒，便于分析处理；L₀是g₀(t)的非零灰度值个数，t_l∈{0,1,...,N-1}·(T/N)是密相介质的像素位置，n＝0,1,…,N-1是像素位置；

步骤1.2：原始信号由于其局限性，通过如LBP算法、Landweber迭代算法和TV正则化算法重建的图像来代替。这些重建的信号是灰度图像，灰度值范围从0到1，所以一个更普遍适用的电磁层析成像信号模型为：

其中L为g(n)的非零灰度值个数，t_l∈{0,1,...,N-1}·(T/N)为密相介质的像素位置，n＝0,1,…,N-1为像素位置，a_l∈[0,1]是像素t_l的灰度值，n(t)是噪声，通过公式(2)得到的g(t)为具有有限新息率的信号，其新息率为ρ＝2L/T。

步骤二，Sinc采样核设计和均匀采样：对有限长度的信号g(t)进行周期延拓使其适用于Sinc采样核。假设周期为T的信号是由Diracs产生的，这种信号通过傅立叶级数自然表示为：

Sinc采样核设计为：

其中，选择B大于或等于ρ＝2L/T的新息率进行均匀采样。通过Sinc采样核可以获取被采样信号的基带频谱信息从而为后续信号重建工作做准备；

步骤三，电磁层析成像系统一维图像信号的有限新息率信号采样过程：输入信号g(t)的有限新息率信号采样结构如图2所示，在对每个输入信号g(t)用Sinc采样核h(t)进行滤波之后，以t_s为采样间隔进行均匀采样，得到样本y_m，然后按照如下公式从样本y_m中获得傅立叶级数系数G[k]：

其中k＝-K,…,K，K＝BT/2；m＝0,1,…,M-1，M为样本总数，t_s＝T/M是均匀采样的时间间隔，原始信号经过采样得到样本y_m，y_m通过离散傅立叶变换提取傅立叶级数系数G[k]得到信息；

步骤四，电磁层析成像系统图像信号的重构过程：通过对原始信号建模及采样获取频谱信息后对所得信息进行分析得到新的测量方程并通过求解稀疏解得到重建后的图像信号。

步骤4.1，根据采样所得频谱信息建立新测量方程：采样阶段之后，从样本y_m中得到输入信号的傅立叶级数系数G[k]，按单周期[0,T)计算，G[k]表示为：

由于实际中噪声的存在，通过对模拟时间轴以Δ＝T/N的时间间隔进行量化，将模拟时间t近似为t＝nΔ，其中n＝0,1,…,N-1，未知参数的时间延迟t_l近似为t_l＝n_lΔ，则公式(6)近似为：

其中，n_l∈{0，1,…,N-1}为时间延迟t_l的离散数值，则测量值看作：

步骤4.2，根据信号的稀疏性建立新的测量方程：在一个周期[0,T)内，忽略量化噪声的情况下，得到一个完整的模拟时间集U＝{0,Δ,2Δ,…,(N-1)Δ}，则时间延迟参数集V＝{n₀Δ,n₁Δ,…,n_L-1Δ}构成了一个更小的子集U，即其中L＜＜N。最终公式(8)改写为：

其中[g₀,g₁,…,g_N-1]^T是一个N×1矢量，由L振幅参数和N-L零值组成。我们的目标是从测量值中找到向量[g₀,g₁,…,g_N-1]^T的非零部分，为了简单起见，公式(9)改写为：

U＝Cg (10)

其中U＝[G₁[-K],…,G₁[K]]^T∈R^Q×1是测量向量，Q＝2K+1；C是由集合组成的Q×N矩阵，g＝[g₁,g₂,…,g_N]^T是电磁层析成像的一个N×1灰度向量。

步骤4.3，给定公式(11)：

V＝Sg (11)

将两个测量方程进行重组：结合公式(10)和(11)，电磁层析成像系统的测量方程重新描述为：

设新的测量向量为λ＝[V；U]，新的灵敏度矩阵为Φ＝[S；C]，公式(5)可以简化为：

λ＝Φg (13)

其中，λ是(D+Q)×1的向量，Φ是(D+Q)×N的矩阵；

步骤4.4，针对测量方程解的稀疏性选择OMP算法进行稀疏解求解过程得到重构后的图像信号：正交变换可以将信号g转化为稀疏信号，正交变换过程描述为：

g＝Ψs (14)

结合公式(13)，公式(14)改写为：

λ＝ΦΨs (15)

由于信号是稀疏的，求解稀疏解的最直接的方法就是求解L0范数最佳化问题：

其中L0范数||s||₀表示向量s中非零系数的个数，再利用正交匹配追踪算法(OMP)迭代寻找λ和ΦΨ的列之间的最大相关性来找到s的非零部分以求解L0范数最优化问题；

步骤4.4，利用稀疏解重构原始图像：在求解稀疏信号s之后，原始图像信号可以估计为：

其中Ψ′是矩阵Ψ′的逆。

本说明书的实施例所述的内容仅仅是对发明构思的实现形式的列举，仅作说明用途。本发明的保护范围不应当被视为仅限于本实施例所陈述的具体形式，本发明的保护范围也及于本领域的普通技术人员根据本发明构思所能想到的等同技术手段。

Claims (2)

1.一种基于频谱特性的电磁层析成像图像重建方法，其特征在于，所述方法包括以下步骤：

步骤一，电磁层析成像系统原始信号的有限新息率信号建模：对于两相流参数检测问题，电磁层析成像系统的原始图像为二值图像，灰度值为0或1；设管道内的截面被分割成d×d像素，其中d为整数，其取值由网格剖分数决定；二值图像表示为二维d×d矩阵或一维N×1向量g，其中N＝d×d；由于EMT系统结构的限制，原始信号不能直接采样或测量，先将原始信号建模为可采样的具有有限新息率的信号，有限新息率即指单位时间内具有有限个参数个数，具有有限新息率的信号可在单位时间内由有限个参数确定，将信号建模为有限新息率信号使其能被均匀采样从而提取出输入信号的特征参数信息；过程为：

步骤1.1：对离散图像信号向量进行有限新息率信号建模为有限长的Diracs流：

其中T是一个周期长度，默认为1秒，便于分析处理；L₀是g₀(t)的非零灰度值个数，t_l∈{0,1,...,N-1}·(T/N)是密相介质的像素位置，n＝0,1,…,N-1是像素位置；

步骤1.2：原始信号由于其局限性，通过如LBP算法、Landweber迭代算法和TV正则化算法重建的图像来代替，这些重建的信号是灰度图像，灰度值范围从0到1，所以一个更普遍适用的电磁层析成像信号模型为：

其中L为g(n)的非零灰度值个数，t_l∈{0,1,...,N-1}·(T/N)为密相介质的像素位置，n＝0,1,…,N-1为像素位置，a_l∈[0,1]是像素t_l的灰度值，n(t)是噪声，通过公式(2)得到的g(t)为具有有限新息率的信号，其新息率为ρ＝2L/T；

步骤二，Sinc采样核设计和均匀采样：对有限长度的信号g(t)进行周期延拓使其适用于Sinc采样核，设周期为T的信号是由Diracs产生的，这种信号通过傅立叶级数自然表示为：

Sinc采样核设计为：

其中，选择B大于或等于ρ＝2L/T的新息率进行均匀采样，通过Sinc采样核可以获取被采样信号的基带频谱信息从而为后续信号重建工作做准备；

步骤三，电磁层析成像系统一维图像信号的有限新息率信号采样过程：在对每个输入信号g(t)用Sinc采样核h(t)进行滤波之后，以t_s为采样间隔进行均匀采样，得到样本y_m，然后按照如下公式从样本y_m中获得傅立叶级数系数G[k]：

其中k＝-K,…,K，K＝BT/2；m＝0,1,…,M-1，M为样本总数，t_s＝T/M是均匀采样的时间间隔，原始信号经过采样得到样本y_m，y_m通过离散傅立叶变换提取傅立叶级数系数G[k]得到信息；

步骤四，电磁层析成像系统图像信号的重构过程：通过对原始信号建模及采样获取频谱信息后对所得信息进行分析得到新的测量方程并通过求解稀疏解得到重建后的图像信号。

2.如权利要求1所述的一种基于频谱特性的电磁层析成像图像重建方法，其特征在于，所述步骤四的过程如下：

步骤4.1，根据采样所得频谱信息建立新测量方程：采样阶段之后，从样本y_m中得到输入信号的傅立叶级数系数G[k]，按单周期[0,T)计算，G[k]表示为：

由于实际中噪声的存在，通过对模拟时间轴以Δ＝T/N的时间间隔进行量化，将模拟时间t近似为t＝nΔ，其中n＝0,1,…,N-1，未知参数的时间延迟t_l近似为t_l＝n_lΔ，则公式(6)近似为：

其中，n_l∈{0，1,…,N-1}为时间延迟t_l的离散数值，则测量值看作：

步骤4.2，根据信号的稀疏性建立新的测量方程：在一个周期[0,T)内，忽略量化噪声的情况下，得到一个完整的模拟时间集U＝{0,Δ,2Δ,…,(N-1)Δ}，则时间延迟参数集V＝{n₀Δ,n₁Δ,…,n_L-1Δ}构成了一个更小的子集U，即其中L＜＜N，最终公式(8)改写为：

其中[g₀,g₁,…,g_N-1]^T是一个N×1矢量，由L振幅参数和N-L零值组成，目标是从测量值中找到向量[g₀,g₁,·,g_N-1]^T的非零部分，为了简单起见，公式(9)改写为：

U＝Cg (10)

其中U＝[G₁[-K],·,G₁[K]]^T∈R^Q×1是测量向量，Q＝2K+1；C是由集合组成的Q×N矩阵，g＝[g₁,g₂,…,g_N]^T是电磁层析成像的一个N×1灰度向量；

步骤4.3，给定公式(11)：

V＝Sg (11)

将两个测量方程进行重组：结合公式(10)和(11)，电磁层析成像系统的测量方程重新描述为：

设新的测量向量为λ＝[V；U]，新的灵敏度矩阵为Φ＝[S；C]，公式(5)简化为：

λ＝Φg (13)

其中，λ是(D+Q)×1的向量，Φ是(D+Q)×N的矩阵；

步骤4.4，针对测量方程解的稀疏性选择OMP算法进行稀疏解求解过程得到重构后的图像信号：正交变换可以将信号g转化为稀疏信号，正交变换过程描述为：

g＝Ψs (14)

结合公式(13)，公式(14)改写为：

λ＝ΦΨs (15)

由于信号是稀疏的，求解稀疏解的最直接的方法就是求解L0范数最佳化问题：

其中L0范数||s||₀表示向量s中非零系数的个数，再利用正交匹配追踪算法OMP迭代寻找λ和ΦΨ的列之间的最大相关性来找到s的非零部分以求解L0范数最优化问题；

步骤4.4，利用稀疏解重构原始图像：在求解稀疏信号s之后，原始图像信号估计为：

其中Ψ′是矩阵Ψ′的逆。