CN110706298B - 正则化加权最小二乘的透射反射双模式超声成像重建方法 - Google Patents

Abstract

本发明涉及一种正则化加权最小二乘的透射反射双模式超声成像重建方法，包括下列步骤：步骤一：获取边界测量值：在被测场域外表面均匀布置一定数量的超声换能器，并分别使用正弦连续波电压与两周期脉冲波电压对发射探头进行激励。步骤二：基于收发超声探头相对几何位置进行直线投影。在透射模式下使用同步代数重建方法SART进行透射预成像。步骤三：使用透射预成像结果对全部收发探头的边界反射点进行聚类分析和曲线拟合，得到光滑内含物边界拟合轮廓，根据轮廓构建边界点约束方程。步骤四：对目标函数使用预优近端梯度下降法进行求解。

Description

正则化加权最小二乘的透射反射双模式超声成像重建方法

技术领域

本发明属于超声层析成像技术领域，涉及实现一种利用正则化加权最小二乘的透射反射双模式超声层析成像重建方法，用于实现被测物场两相介质分布的重建与可视化表征。

背景技术

超声层析成像技术(Ultrasonic Tomography,UT)是一种结构性成像技术，其通过对布置于被测场域外的超声传感器阵列施加一定的激励，并记录传感器接收到的时变电压信号，结合图像重建算法对被测场域内折射系数、衰减系数或声阻抗分布等参数进行可视化表征。相比其他成像方法如电阻抗层析成像(Electrical Impedance Tomography,EIT)、电磁层析成像(Magnetic Impedance Tomography,MIT)、X射线层析成像(X-ray ComputedTomography,X-CT)及光学层析成像方法(Optical Coherence Tomography,OCT)，UT具有成本低、使用安全、结构简单、成像速度快等优势，是一种较为理想的过程可视化检测手段，在多相流可视化检测、化工石油输送、航空发动机探查以及生物医学诊断中均有广泛的应用。

超声层析成像主要包含三个部分：超声换能器阵列、信号激励与采集***、超声成像图像重建算法。其中超声成像算法通过对从采集***得到的换能器接收信号进行处理，通过解调提取测量幅值或渡越时间，得到某个确定激励下的全部换能器的有效测量数据，进一步通过图像重建方法得到被测物场两相介质分布的合理估计。目前，超声成像重建算法作为一种主要利用硬场特性的成像方法，超声成像方法严重依赖于场域边界换能器的数量，其逆问题求解具有严重的病态性(对测量值得微小扰动会导致重建结果的大幅度变化)和欠定性(所需求解的方程数远小于未知量的数目，方程有无穷多解)。为克服这个问题，学者们提出了许多图像重建算法，超声层析成像图像重建方法根据其敏感效应的不同主要分为透射模式和反射模式：其中透射模式记录激励源对侧超声探头接收到的时变电压信号实现渡越时间、衰减幅值的提取，进而对介质分布进行重建；反射模式记录激励源及其同侧超声探头接收到的时变电压信号，提取渡越时间以重建两相介质中离散相边界位置，提取信号幅值以确定边界声阻抗差异。在透射模式超声成像算法研究中，典型的图像重建方法有徐立军等人1998年在《仪器仪表学报》第17卷，1-7页，发表的题为《气液两相泡状流体监测用超声层析成像***的研究》(Investigation of Ultrasound Tomography System usedfor Monitoring Bubbly Gas/Liquid Two-phase fluid)的文章中提到的二值反投影方法、Rahim等人2007年在《传感器与执行器》(Sensors and Actuators)第135卷，337-345页发表的题为《超声透射层析成像技术对液体、气体的非侵入性成像》(Non-invasiveimaging of liquid/gas flow using ultrasonic transmission mode tomography)的文章中提到的采用滤波函数的线性反投影方法、Rahim等人2017年在《IEEE工业电子汇刊》(IEEETransactions on Industrial Electronics)第64卷，第12期，9636-9645页发表的题为《导电管壁内气泡柱上升过程的相含率检测：仿真和实验研究》(Simulative andExperimental Studies:Void Fraction Detection in Bubble Column of a ConductingPipe)的文章中提出的超声阈值滤波方法、苏邦良等人2000年在《化学工程》(ChemicalEngineering Journal)第77卷，37-41页发表的题为《同步迭代重建技术在电容层析成像中的应用》(The use of simultaneous iterative reconstruction technique forelectrical capacitance tomography)的文章中提出的同步迭代重建方法、等；在反射模式超声成像算法研究中，典型的图像重建算法有徐春广等人在《仪表技术与传感器》2013年第4期中发表的题为《基于VC++的超声探伤仪发射端信号测量与处理***》的文章中提出的反射式超声椭圆成像算法、Langener等人2017年在《流动测量与仪表》(Flow Measurementand Instrumentation)第53卷，107-115页发表的题为《一种基于超声反射重建技术的实时超声过程层析成像***》(Areal-time ultrasound process tomography system usingareflection mode reconstruction technique)的文章中提出的基于路径追踪的超声反射图像重建算法。

上述超声成像算法研究中，往往根据被测场域内内含物分布特性单独使用透射模式或反射模式其中一种，透射模式和反射模式。实际上，透射模式与反射模式有较为互补的成像效果，如反射模式对靠近边界的内含物重建效果较好，透射模式对靠近中央的内含物重建效果较好，单一模式下不同位置内含物的成像效果差别较大，无法做到对任意位置离散相介质进行快速、准确的可视化表征。

此外，探头数量对超声成像算法重建精度及分辨率有至关重要的影响：投影路径越多，成像精度越高。但在UPT实际应用过程中，受限于场域尺寸及测试实时性要求，超声换能器数目不能无限增加。单独使用透射法或反射法进行图像重建均仅能利用超声信号中包含的部分信息，无法在保证较快的数据成像速度前提下实现较高精度的图像重建。针对该问题，2017年谭超等人在《IEEE工业电子汇刊》(IEEE Transactions on IndustrialElectronics)第66卷，9539-9548页发表的题为《针对工业多相流的超声透反射层析成像***》(An Ultrasonic Transmission/Reflection Tomography System for IndustrialMultiphase Flow Imaging)中提出在超声层析成像中将透射反射模式进行图像融合，即分别在透射模式与反射模式下进行图像重建，将成像结果对应位置像素相加进行图像融合。该方法在探头数量较低时受反射模式渡越时间提取误差影响较大，其图像重建精度相比于透射模式及反射模式提升较为有限。

发明内容

本发明的目的在于提供一种需要一种在低投影数量下透射反射多模式的基于信息融合的图像重建算法，以期在保证实时性的同时实现高精度的图像重建，用于实现被测物场两相介质分布的重建和可视化表征。本发明在透射及反射模式下分别使用连续波与脉冲波两种激励方式，对透射信号和反射信号分别进行幅值和渡越时间提取，获得透射衰减系数及反射点位置。使用透射预成像结果对全部接收信号提取出的反射点进行分类和轮廓拟合，进而得到用于约束联合图像重建的边界约束方程。使用正则化加权最小二乘法对透射模式下的同步代数重建方法进行改造，进而构建适用于拉格朗日乘子框架下等式约束目标函数。使用预优近端梯度下降法对等式约束目标函数进行迭代求解，得到内含物分布重建结果。技术方案如下：

一种正则化加权最小二乘的透射反射双模式超声成像重建方法，包括下列步骤：

步骤一：获取边界测量值：在被测场域外表面均匀布置一定数量的超声换能器，并分别使用正弦连续波电压与两周期脉冲波电压对发射探头进行激励。

使用正弦连续波电压激励时，对除发射外的其他超声探头，记录其接收到正弦信号连续多个单峰值的平均值，记为透射边界电压测量值V，据此获取重建所需的投影衰减测量值τ：

其中，f_c是激励信号的中心频率，V_s为场域内充满单一背景介质时的边界电压测量值，V_r为场域内存在离散的内含物介质情况下的边界电压测量值，ln表示自然对数符号。

使用两周期脉冲波电压激励时，对发射探头及其相邻两探头，记录其接收到时变电压信号并使用小波阈值收缩法进行去噪，使用赤池信息准则对时变信号的各个回波进行渡越时间计算，并根据收发探头几何位置转换为边界反射点。

步骤二：基于收发超声探头相对几何位置进行直线投影，构建系数矩阵R，其元素R_i,j表示场域内第i条投影路径穿过场域内第j个像素的相对长度，同时对应系数矩阵中第i行、第j列的元素。在透射模式下使用同步代数重建方法SART进行透射预成像，具体计算方式为：

a^k+1＝a^k+λD_p(D_rR)^T(τ-Ra^k)

其中，λ表示迭代步长，k表示迭代次数，a表示被测场域内像素单元衰减系数分布，a^k为对a进行求解时的第k步迭代值，D_p＝diag(1/R_+,1,1/R_+,2,····,1/R_+,N)，D_r＝diag(1/R_1,+,1/R_2,+,····,1/R_M,+)，R_+,1表示对系数矩阵内第一列所有元素进行求和，R_1,+表示对系数矩阵内第一行所有元素进行求和，N表示场域内像素个数，M表示场域内投影路径数量，diag()表示对角矩阵，()^T表示矩阵转置。

步骤三：使用透射预成像结果对全部收发探头的边界反射点进行聚类分析和曲线拟合，得到光滑内含物边界拟合轮廓，根据轮廓构建边界点约束方程，要求边界的梯度下降值取为较大的数值，其他区域的梯度下降值取为较小的数值，以此构建反射模式约束方程G(a^k)表示为：

G(Δa^k)＝H(a^k)-J_H(a^k)·Δa^k＝0

其中J_H(a^k)表示在第k次迭代过程中，方程H(a^k)的一阶偏微分矩阵，Δa^k＝a^k+1-a^k，H(a^k)为：

H(a)＝[a(p_i)-α·a_b,...,a(p_b)-a_b,...,a(p_o)-α^-1·a_b,...,a(p)-a]＝0

其中，p_i表示拟合轮廓上像素内法向向量指向的像素单元，p_b表示在拟合轮廓上的像素单元，p_o表示拟合轮廓上像素外法向向量指向的像素单元，a_b表示在拟合轮廓上的像素单元的衰减系数，α为梯度下降取值，选取为10⁵。

步骤四：使用正则化加权最小二乘法对透射模式下的同步代数重建方法进行改造，构建适用于拉格朗日乘子框架下等式约束目标函数，并结合反射模式约束方程对目标函数使用预优近端梯度下降法进行求解。

[1]使用正则化加权最小二乘法将同步代数重建方法改造为待优化目标函数，具体计算方式为：

其中，a^*为待求解的最优衰减系数分布，

表示满足表达式取得最小值时a的取值，

表示D_r阶范数的平方，

表示二范数的平方，β表示预先选定的正则化参数，取为0.001，Q为正则化矩阵。

[2]构建拉格朗日乘子函数并在其极小值点求解Δa^k，具体计算方式为：

其中，

表示偏微分求解。根据上述目标函数和约束方程，上式可以展开并整理为：

使用

表示上式，并使用预优近端梯度下降法进行求解，得到Δa^k。

[3]根据求得的Δa^k更新各个像素单元衰减系数分布，具体计算方式为：

a^k+1＝a^k+Δa^k

步骤五：重复步骤四直至残差满足要求

其中，Rea^k＝||S·a^k-τ||表示残差值,ε为人为设定的残差阈值。

本发明的有益效果是将超声透射模式与反射模式进行基于信息融合的双模式图像重建，提出了一种正则化加权最小二乘的透射反射双模式超声层析成像重建方法，对超声层析成像图像重建算法进行了补充和发展，解决了传统超声透射成像或超声反射成像针对局部位置成像失败的问题，实现了测量中对透射及反射信息的同步提取及联合利用。与此同时，所提方法将逆问题的典型迭代逆推求解问题转化成为目标函数的等式约束优化问题，拓展了拉格朗日乘子框架在多模态图像重建中的应用，大幅提高超声层析成像的重建精度和图像重建质量。

附图说明

图1为正则化加权最小二乘的透射反射双模式超声层析成像重建方法基本流程图；

图2为本发明中采用透射模式和反射模式激励测量原理示意图；

图3为本发明中针对离散相介质分布(Distribution)的实验图像重建结果，并分别给出了相应的透射法成像结果(T-mode)、反射法成像结果(R-mode)、本算法所提出双模式成像结果(Dual-mode)；

图4为针对本发明的十组仿真模型不同成像结果的相对误差(Relative Error,RE)和相关系数(Correlation Coefficient,CC)对比。

具体实施方式

结合附图和实施例对本发明的基于非均匀形状约束的像素值域滤波超声成像算法加以说明。

本发明的正则化加权最小二乘的透射反射双模式超声层析成像重建方法，实施例中针对工业输油管道中气水两相流的成像这一UT技术的常见应用形式，使用基于投影路径长度的系数矩阵构建方法表征油水两相流超声可视化检测的正问题模型，同时将图像重建逆问题的求解过程分解为透射模式预成像、反射边界点提取及拟合、透反射目标函数及约束方程构建、拉格朗日乘子函数构建及预优近端梯度下降法求解四个部分，达到透射反射双模态信息同步提取及联合利用的目的，以提高低投影数量下的对离散相介质位置、轮廓的准确描述。

图1为本发明的基于非均匀形状约束的像素值域滤波超声成像算法完整流程图。算法主要分为透射模式预成像、反射边界点提取及拟合、透反射目标函数及约束方程构建、拉格朗日乘子函数构建及预优近端梯度下降法求解四个部分。透射模式所需系数矩阵在已知场域分布、探头尺寸和探头布置位置的基础上，由数值计算方法获得。

图2为透射反射双模式超声层析成像激励测量基本原理示意图，在对进行测量时，共计16个超声换能器均匀的沿管壁安装负责激励、接收超声波。透射模式采用循环激励、一发全收的测量模式，探头按顺时针方向均匀分布。16个超声探头按顺序接入峰峰值50V、频率1MHz的方波电压激励，探头通道切换时间间隔2.5ms。四与此同时，16个通道同步接收稳态时刻的电压正弦信号持续1ms，并通过正交解调得到接收电压有效值。每次透射模式测量总计获得16×15＝240组边界电压测量数据。反射模式采用两周期正弦脉冲波对每个超声探头进行循环激励，每次激励后，激励探头和其相邻两探头记录时变电压信号。每次反射模式测量总计获得16×3＝48组时变电压测试数据。

图3中分别给出了模型1-模型10对应的成像结果。可以看出，传统透射模式或反射模式均无法有效实现不同位置、不同数目、不同尺寸离散相介质的有效重建。本发明所提出的算法能够较准确地重建出离散相介质的准确位置和完整形状。重建结果中，离散相介质边界清晰，图像无过多伪影及噪声，图像分辨率、成像精度均有明显提高。

图4中分别给出了传统超声层析成像图像重建算法的成像结果与本算法的成像结果重建指标对比，重建指标包括相对误差(RE)和图像相关系数(CC)两种，其计算方法表示为：

其中，σ表示重建的像素单元衰减分布，σ^*表示真实情况下衰减分布，σ_j和σ_j ^*表示第j个像素单元重建的和真实的衰减分布，

和

表示重建的和真实的衰减分布的平均值。

本算法实施例包括如下具体步骤：

(1)：获取边界测量值，在被测场域外表面均匀布置一定数量的超声换能器，并分别使用正弦连续波电压与两周期脉冲波电压对发射探头进行激励。

使用正弦连续波电压激励时，对除发射外的其他超声探头，记录其接收到正弦信号连续50个单峰值的平均值，记为透射边界电压测量值V。据此获取重建所需的投影衰减测量值τ，具体计算方式为：

其中，f_c是激励信号的中心频率，V_s为场域内充满单一背景介质时的边界电压测量值，V_r为场域内存在离散的内含物介质情况下的边界电压测量值，ln表示自然对数符号。

使用两周期脉冲波电压激励时，对发射探头及其相邻两探头，记录其接收到时变电压信号并使用小波阈值收缩法进行去噪，使用赤池信息准则对时变信号的各个回波进行渡越时间计算，并根据收发探头几何位置转换为边界反射点。

(2)：基于收发超声探头相对几何位置进行直线投影，构建系数矩阵R，其元素R_i,j表示场域内第i条投影路径穿过场域内第j个像素的相对长度，同时对应系数矩阵中第i行、第j列的元素。在透射模式下使用同步代数重建方法(SART)进行透射预成像，具体计算方式为：

a^k+1＝a^k+λD_p(D_rR)^T(τ-Ra^k)

其中，λ表示迭代步长，k表示迭代次数，a表示被测场域内像素单元衰减系数分布，a^k为对a进行求解时的第k步迭代值，D_p＝diag(1/R_+,1,1/R_+,2,····,1/R_+,N)，D_r＝diag(1/R_1,+,1/R_2,+,····,1/R_M,+)，R_+,1表示对系数矩阵内第一列所有元素进行求和，R_1,+表示对系数矩阵内第一行所有元素进行求和，N表示场域内像素个数，M表示场域内投影路径数量，diag()表示对角矩阵，()^T表示矩阵转置。

(3)：使用透射预成像结果对全部收发探头的边界反射点进行聚类分析和曲线拟合，得到光滑内含物边界拟合轮廓。根据轮廓构建边界点约束方程，要求边界的梯度下降值取为较大的数值，其他区域的梯度下降值取为较小的数值。以此构建反射模式约束方程G(a^k)表示为：

G(Δa^k)＝H(a^k)-J_H(a^k)·Δa^k＝0

其中J_H(a^k)表示在第k次迭代过程中，方程H(a^k)的一阶偏微分矩阵，Δa^k＝a^k+1-a^k，H(a^k)为：

H(a)＝[g(a_i)-α·a_b,...,g(a_b)-a_b,...,g(a_o)-α^-1·a_b,...,a(p)-a]＝0

其中，p_i表示拟合轮廓上像素内法向向量指向的像素单元，p_b表示在拟合轮廓上的像素单元，p_o表示拟合轮廓上像素外法向向量指向的像素单元，a_b表示在拟合轮廓上的像素单元的衰减系数，α为梯度下降取值，选取为10⁵。

(4)：使用正则化加权最小二乘法对透射模式下的同步代数重建方法进行改造，构建适用于拉格朗日乘子框架下等式约束目标函数，并结合反射模式约束方程对目标函数使用预优近端梯度下降法进行求解。

a.使用正则化加权最小二乘法将同步代数重建方法改造为待优化目标函数，具体计算方式为：

其中，a^*为待求解的最优衰减系数分布，

表示满足表达式取得最小值时a的取值，

表示D_r阶范数的平方，

表示二范数的平方，β表示预先选定的正则化参数，一般取为0.001，Q为正则化矩阵。

b.构建拉格朗日乘子函数并在其极小值点求解Δa^k，具体计算方式为：

其中，

表示偏微分求解。根据上述目标函数和约束方程，上式可以展开并整理为：

使用

表示上式，并使用预优近端梯度下降法进行求解，得到Δa^k。

c.根据求得的Δa^k更新各个像素单元衰减系数分布，具体计算方式为：

a^k+1＝a^k+Δa^k

(5)：重复步骤四直至残差满足要求

其中，Rea^k＝||S·a^k-τ||表示残差值,ε为人为设定的残差阈值。

Claims (1)

1.一种正则化加权最小二乘的透射反射双模式超声成像重建方法，包括下列步骤：

步骤一：获取边界测量值：在被测场域外表面均匀布置一定数量的超声换能器，并分别使用正弦连续波电压与两周期脉冲波电压对发射探头进行激励；

使用正弦连续波电压激励时，对除发射外的其他超声探头，记录其接收到正弦信号连续多个单峰值的平均值，记为透射边界电压测量值V，据此获取重建所需的投影衰减测量值τ：

其中，f_c是激励信号的中心频率，V_s为场域内充满单一背景介质时的边界电压测量值，V_r为场域内存在离散的内含物介质情况下的边界电压测量值，ln表示自然对数符号；

使用两周期脉冲波电压激励时，对发射探头及其相邻两探头，记录其接收到时变电压信号并使用小波阈值收缩法进行去噪，使用赤池信息准则对时变信号的各个回波进行渡越时间计算，并根据收发探头几何位置转换为边界反射点；

步骤二：基于收发超声探头相对几何位置进行直线投影，构建系数矩阵R，其元素R_i,j表示场域内第i条投影路径穿过场域内第j个像素的相对长度，同时对应系数矩阵中第i行、第j列的元素；在透射模式下使用同步代数重建方法SART进行透射预成像，具体计算方式为：

a^k+1＝a^k+λD_p(D_rR)^T(τ-Ra^k)

其中，λ表示迭代步长，k表示迭代次数，a表示被测场域内像素单元衰减系数分布，a^k为对a进行求解时的第k步迭代值，D_p＝diag(1/R_+，1，1/R_+，2，...，1/R_N，+),，D_r＝diag(1/R_1,+,1/R_2,+,...,1/R_M,+)，R_+,1表示对系数矩阵内第一列所有元素进行求和，R_1,+表示对系数矩阵内第一行所有元素进行求和，N表示场域内像素个数，M表示场域内投影路径数量，diag()表示对角矩阵，()^T表示矩阵转置；

步骤三：使用透射预成像结果对全部收发探头的边界反射点进行聚类分析和曲线拟合，得到光滑内含物边界拟合轮廓，根据轮廓构建边界点约束方程，要求边界的梯度下降值大于其他区域的梯度下降值，以此构建反射模式约束方程G(a^k)表示为：

G(Δa^k)＝H(a^k)-J_H(a^k)·Δa^k＝0

其中J_H(a^k)表示在第k次迭代过程中，方程H(a^k)的一阶偏微分矩阵，Δa^k＝a^k+1-a^k，H(a^k)为：

H(a)＝[a(p_i)-α·a_b,...,a(p_b)-a_b,...,a(p_o)-α^-1·a_b,...,a(p)-a]＝0

其中，p_i表示拟合轮廓上像素内法向向量指向的像素单元，p_b表示在拟合轮廓上的像素单元，p_o表示拟合轮廓上像素外法向向量指向的像素单元，a_b表示在拟合轮廓上的像素单元的衰减系数，α为梯度下降取值，选取为10⁵；

步骤四：使用正则化加权最小二乘法对透射模式下的同步代数重建方法进行改造，构建适用于拉格朗日乘子框架下等式约束目标函数，并结合反射模式约束方程对目标函数使用预优近端梯度下降法进行求解；

[1]使用正则化加权最小二乘法将同步代数重建方法改造为待优化目标函数，具体计算方式为：

其中，a^*为待求解的最优衰减系数分布，表示满足表达式取得最小值时a的取值，表示D_r阶范数的平方，表示二范数的平方，β表示预先选定的正则化参数，取为0.001，Q为正则化矩阵；

[2]构建拉格朗日乘子函数并在其极小值点求解Δa^k，具体计算方式为：

其中，表示偏微分求解；根据上述目标函数和约束方程，上式可以展开并整理为：

使用表示上式，并使用预优近端梯度下降法进行求解，得到Δa^k；

[3]根据求得的Δa^k更新各个像素单元衰减系数分布，具体计算方式为：

a^k+1＝a^k+Δa^k

步骤五：重复步骤四直至残差满足要求

其中，Rea^k＝||S·a^k-τ||表示残差值，ε为人为设定的残差阈值。

Images