CN109598769B - 总变差正则化约束的超声成像同步代数迭代重建方法 - Google Patents

Abstract

本发明涉及一种总变差正则化约束的超声成像同步代数迭代重建方法，用于超声层析成像，包含：获取重建所需的投影衰减测量值；构建系数矩阵，在考虑投影路径与场域像素重叠面积比例的同时考虑重叠区域相对于场域像素的几何位置；使用同步代数重建方法进行成像预处理；搭建正则化加权最小二乘法框架并将总变差正则化作为正则化项，将同步代数迭代重建方法作为代数项，进行逆问题迭代重建计算；迭代直至残差满足要求。

Description

总变差正则化约束的超声成像同步代数迭代重建方法

技术领域

本发明属于超声层析成像技术领域，用于实现场域内弱声阻抗比介质分布衰减系数的重建。

背景技术

超声层析成像技术(Ultrasonic Tomography,UT)是一种结构性成像技术，其通过在被测场域外布置超声传感器阵列并施加一定的激励以得到边界电压测量数据，以此来重建被测场域内部的折射系数、衰减系数或声阻抗分布情况。相比软场成像技术例如电阻抗层析成像(Electrical Impedance Tomography,EIT)和电磁层析成像(MagneticImpedance Tomography,MIT)，UT具有非侵入、分辨率高的优点，相比精度较高的硬场成像技术如X射线层析成像(X-ray Computed Tomography,X-CT)及光学层析成像方法(OpticalCoherence Tomography,OCT)，UT使用安全、结构简单、可以实现实时成像。此外UT还有着非接触、方向性好、成本低等优势，是一种较为理想的过程可视化检测监视手段。UT作为一种层析成像技术手段，在多相流可视化检测、化工石油输送、航空发动机探查以及生物医学诊断中均有广泛的应用。

完整的UT系统主要包含三个部分：传感器阵列设计与换能器安装；信号激励、采集系统；超声成像重建算法。其中超声成像算法通过对从采集系统得到的换能器接收信号进行处理，通过解调提取测量幅值或渡越时间，得到某个确定激励下的全部换能器的有效测量数据，进一步通过图像重建方法得到场域内含物介质分布的合理估计。目前，超声成像重建算法主要存在成像分辨率低、成像精度差、图像伪影严重等三个方面的问题。此外，作为一种主要利用硬场特性的成像方法，超声成像方法严重依赖于场域边界换能器的数量，其逆问题求解具有严重的病态性(对测量值得微小扰动会导致重建结果的大幅度变化)和欠定性(所需求解的方程数远小于未知量的数目，方程有无穷多解)。为克服这个问题，学者们提出了许多图像重建算法，其中，基于路径的投影重建算法是一种克服病态性的有效手段。该方法通过计算激励、接收换能器间的路径，将收发探头间的时间延迟或幅值衰减均匀的分配给所计算路径上的每一个像素，通过对不同收发探头间的路径进行计算和对同一像素在不同路径上的估计值进行叠加，得到场域内每个像素值的有效估计，以达到可视化测量和图像重建的目的。典型的投影重建方法有徐立军等人1998年在《仪器仪表学报(ChineseJournal of Scientific Instrument)》第17卷，1-7页，发表的“气液两相泡状流体监测用超声层析成像系统的研究(Investigation of Ultrasound Tomography System used forMonitoring Bubbly Gas/Liquid Two-phase fluid)”中提到的二值反投影方法、Rahim等人在《Sensors and Actuators(传感器与执行器)》第135卷，337-345页发表的“Non-invasive imaging of liquid/gas flow using ultrasonic transmission modetomography(超声对液体、气体的非侵入性成像)”中提到的采用R-L函数的线性滤波反投影方法、Gordon等人在《Journal of theoretical biology(理论生物学杂志)》第29卷，第3期，471-481页发表的“Algebraic reconstruction techniques(ART)for three-dimensional electron microscopy and X-ray photography(用于三维电子显微镜和X射线CT的代数重建技术)”中提出的代数重建方法、苏邦良等人在《Chemical EngineeringJournal(化学工程期刊)》第77卷，37-41页发表的“The use of simultaneous iterativereconstruction technique for electrical capacitance tomography(同步迭代重建技术在电容层析成像中的应用)”中提出的同步迭代重建方法、Anderson等人在《UltrasoundImaging(超声成像)》第6卷，81-94页发表的“Simultaneous algebraic reconstructiontechnique(SART):a superior implementation of the ART algorithm(同步代数重建技术(SART)：ART算法更优越的实现)”中提出的同步代数重建方法等。其中，SART算法以其收敛快速、残差较小的优点被广泛使用在生物组织超声成像的研究中。

正则化方法是一种克服病态性的有效手段，通过构造正则化惩罚项的方式，将一定的先验信息融入到逆问题的求解过程中，约束解的搜索空间，指引解的优化方向，以达到改善病态性的目的。正则化方法常用于电阻抗层析成像逆问题重建算法中，如M.Vauhkonen等人1998年在《IEEE Transactions on Medical Imaging(IEEE医学成像期刊)》第17卷，285-293页，发表的“Tikhonov regularization and prior information in electricalimpedance tomography(在电阻抗成像中的吉洪诺夫正则化和先验信息)”中提到的吉洪诺夫正则化方法、J.Zhao等人在《IEEE International Conference on Imaging Systemsand Techniques(IEEE国际成像系统与技术会议)》第25-30页发表的“Sparseregularization for small objects imaging with electrical resistancetomography(电阻层析成像中小物体成像的稀疏正则化方法)”中提到的L1正则化方法、A.Borsic等人2007年在《Inverse Problems(逆问题)》第99卷，A12-A12页发表的“Totalvariation regularization in electrical impedance tomography(电阻抗成像中的总变差正则化方法)”中提出的总变差正则化(TV)方法等。不同的正则化方法往往引入不同类型的先验信息，比如Tikhonov先验引入了均匀分布信息，Laplace先验引入的光滑性信息以及TV正则化引入的边界梯度信息。正则化方法与超声成像重建算法的结合将有效提高对场域内含物边界进行重建的能力，并对成像结果中的伪影和噪声进行有效抑制。

在众多超声重建算法的应用对象中，以油水两相流为代表的液液两相流在石油化工、航空航天等领域广泛存在，是层析成像的常见对象之一。该类对象在声学特性上表现为密度相近，声速差异小，由于其黏性差异声衰减系数有较大区别。这对超声层析成像的边界电压测量和图像重建算法的准确性都提出了更高的要求，如较深的传播距离、较窄的激励信号频带，场域内部尽可能少的多重散射等等。与此同时，在UT应用于实际生产过程中时，由于采集速度和对象尺寸的关系，超声换能器数目一般较少。较高精度的超声成像需求和较弱声阻抗比的内含物分布形成了较大的矛盾。因此需要一种在低投影数量下针对弱声阻抗比介质的高精度超声图像重建算法。

发明内容

本发明的目的是提供一种能对弱声阻抗比介质分布进行较高精度有效重建的总变差正则化约束的超声成像同步代数迭代重建方法。技术方案如下：

一种总变差正则化约束的超声成像同步代数迭代重建方法，用于超声层析成像，包含以下步骤：

步骤一：在被测场域外表面均匀的布置一定数量的超声换能器，并使用连续正弦波电压对发射探头进行激励，对除发射外的其他超声探头，记录其接收到正弦信号连续多个单峰值的平均值，记为边界电压测量值，从而获取重建所需的投影衰减测量值τ；

步骤二：构建系数矩阵R，在考虑投影路径与场域像素重叠面积比例的同时考虑重叠区域相对于场域像素的几何位置，计算公式为：

其中，R_i,j是场域内第i条投影路径穿过场域内第j个像素的相对长度，同时对应系数矩阵中第i行、第j列的元素，d_diag表示像素单元中心相对于投影面积几何中心在对角线方向上的距离，l_pixel为像素对角线长度，若第i条投影路径与第j个像素无投影重叠，则R_ij＝0，area_overlap表示重叠区域面积，area_pixel表示像素单元面积；

步骤三：使用同步代数重建方法进行成像预处理：

[1]使用总变差正则化方法计算迭代重建正则化约束项Q(a)，其计算方式可以表示为：

其中，a表示场域内像素单元的衰减系数，

表示场域内第j个像素在横向方向上的灰度变化特征，

表示场域内第j个像素在纵向方向上的灰度变化特征，β表示为防止像素值为零时的不可微现象而加入的一个较小数，Δp_j表示场域内第j个像素单元的面积，N代表所有像素单元的个数；

[2]计算同步代数迭代重建过程中的系数矩阵R，具体计算过程中可以表示为：

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

其中，α表示迭代步长，k表示迭代次数，a^(k)表示第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列的全部元素之和，R_1,+表示系数矩阵R中第1行的全部元素之和，N表示场域内像素个数，M表示场域内投影路径数量，diag()表示对角矩阵，()^T表示矩阵转置；

[3]列出预优牛顿最速下降法的基本框架，具体计算方法表示为：

其中，P为非负约束矩阵，D为预优牛顿最速下降法中的权值矩阵，

表示第k次迭代中求解目标函数的梯度，H(Ψ(a^(k)))表示第k次迭代中目标函数的二阶偏微分矩阵，即海瑟矩阵，Ψ(a^(k))表示第k次迭代的目标矩阵；

步骤四：搭建正则化加权最小二乘法框架并将总变差正则化作为正则化项，将同步代数迭代重建方法作为代数项，进行逆问题迭代重建计算；

[1]给出基于正则化加权最小二乘法框架的求解方程和待优化目标函数Ψ(a)：

其中，a^*表示带求解的像素单元的衰减系数分布，

表示满足目标函数取最小值时a的取值；

表示W-范数的平方，W表示预先确定的加权范数矩阵，μ表示正则化参数，

表示二范数的平方，Q(a)表示正则化项；

[2]进行非负约束矩阵P的计算，其计算方式表示为；

其中，a_j表示第j个像素单元的重建衰减系数值；

[3]根据步骤二中的预处理结果，进行基于总变差正则化加权最小二乘框架下同步代数重建算法迭代计算，其具体计算方式可以表示为：

其中，L_β(a^(k))＝D_x ^TLD_x+D_y ^TLD_y，L表示一个N×N维，对角线元素为

的对角矩阵，D_x和D_y分别表示将

和

两个列向量转换为N×N维对角矩阵的形式；

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

其中，Rea^(k)＝||R·a^(k)-τ||，ε为人为设定的残差阈值。

本发明在给出透射模态下超声衰减系数重建逆问题模型的基础上，提出一种基于几何位置修正的重叠面积系数矩阵构建方法表征弱声阻抗比介质分布超声可视化检测模型，提高正问题模型的精确性并显著减小非有效投影路径上像素参与逆问题计算的规模；同时在图像重建逆问题的迭代求解过程中搭建正则化加权最小二乘法框架，将同步代数迭代方法纳入其中的代数项，将总变差正则化方法纳入其中的正则化项，以达到在较高精度重建内含物位置的同时，有效进行灵敏度均一化约束，实现低投影数量下弱声阻抗比介质轮廓的准确重建。所提出总变差正则化约束的超声成像同步代数迭代重建方法，其核心思想是“搭建正则化加权最小二乘法框架并将总变差正则化和同步代数重建方法融合”其中：通过对加权最小二乘的变形推导，提出了融合正则化方法和代数重建方法的基本框架；通过将同步代数迭代方法纳入其中的代数项，将总变差正则化方法纳入其中的正则化项，以达到在以较高精度重建出内含物位置的同时，有效进行灵敏度均一化约束，以提高低投影数量下对内含物位置、轮廓的准确重建。该算法可以在重建结果中保留较为清晰和准确的内含物形状结构，在保证实时成像的基础上显著提升UT的成像精度。

附图说明

图1为总变差正则化约束的超声成像同步代数迭代重建方法基本流程图；

图2为本发明中采用连续波成像的超声层析成像测量方法示意图；

图3为本发明的系数矩阵构建方式示意图，主要分为投影路径与像素单元重叠面积的计算与重叠区域同像素单元几何位置关系判定两部分；

图4为本发明的四个典型仿真模型，并分别给出了相应的总变差正则化(TV)成像结果、滤波反投影(FBP)成像结果、同步代数重建(SART)成像结果及本发明算法(RW-SART)的最终成像结果。

图5为针对本发明四组仿真模型不同成像结果的相对误差和图像相关系数对比。

具体实施方式

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

本发明的总变差正则化约束的超声成像同步代数迭代重建方法，实施例中针对工业输油管道中油水两相流的成像这一UT技术的常见应用形式，使用基于几何位置修正的重叠面积系数矩阵构建方法表征弱声阻抗比两相流超声可视化检测的正问题模型，同时在图像重建逆问题的迭代求解过程中搭建正则化加权最小二乘法框架，将同步代数迭代方法纳入其中的代数项，将总变差正则化方法纳入其中的正则化项，以达到在以较高精度重建出内含物位置的同时，有效进行灵敏度均一化约束，以提高低投影数量下对内含物位置、轮廓的准确重建。

图1为总变差正则化约束的超声成像同步代数迭代重建方法所提算法的基本流程图，主要分为基于几何位置修正的重叠面积系数矩阵构建、正则化与同步迭代预计算、预优牛顿法下降计算及循环迭代等部分。

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

图3为基于几何位置修正的重叠面积系数矩阵构建方法，如图所示，在计算投影面积和像素单元重叠面积的同时，考虑重叠区域中心相对于像素单元中心的几何位置并以其对角线投影距离作为几何位置修正系数。

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

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

和

表示重建的和真实的电导率分布的平均值。

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

(1)在被测场域外表面均匀的布置一定数量的超声换能器，并使用正弦波电压对发射探头进行激励。对除发射外的其他超声探头，记录其接收到正弦信号连续50个单峰值的平均值，记为边界电压测量值。据此，获取重建所需的投影衰减测量值τ，具体计算方式为

式中f_c是激励信号的中心频率，A_s为空场下的边界电压测量值，A_r为存在内含物情况下的边界电压测量值，ln表示对数符号。

(2)构建系数矩阵，本发明中提出全新基于几何位置修正的系数矩阵构建方法。即在考虑投影路径与场域像素重叠面积比例的同时考虑重叠区域相对于场域像素的几何位置，具体计算公式为：

式中R_i,j是场域内第i条投影路径穿过场域内第j个像素的相对长度，同时对应系数矩阵中第i行、第j列的元素，d_diag表示像素单元中心相对于投影面积几何中心在对角线方向上的距离，l_pixel为像素对角线长度，若第i条投影路径与第j个像素无投影重叠，则R_ij＝0；

(3)使用同步代数重建方法进行成像预处理：

a.使用总变差正则化方法计算迭代重建正则化项；(其中k表示当前迭代次数)，具体计算方式表示为：

其中，a表示场域内像素单元的衰减系数，

表示场域内第j个像素在横向方向上的灰度变化特征，

表示场域内第j个像素在纵向方向上的灰度变化特征，β表示为防止像素值为零时的不可微现象而加入的一个较小数，Δp_j表示场域内第j个像素单元的面积。

b.计算同步代数迭代重建过程中的相关矩阵，具体计算过程中可以表示为：

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

其中，α表示迭代步长，k表示迭代次数，a^(k)表示第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,+)，N表示场域内像素个数，M表示场域内投影路径数量，diag()表示对角矩阵，()^T表示矩阵转置。

c.计算预优牛顿最速下降法的基本框架，具体计算方法表示为：

其中，P为非负约束矩阵，D为预优牛顿最速下降法中的权值矩阵，

表示第k次迭代中求解目标函数的梯度，H(Ψ(a^(k)))表示第k次迭代中目标函数的二阶偏微分矩阵，即海瑟矩阵，Ψ(a^(k))表示第k次迭代的目标矩阵。

(4)根据(3)的计算结果，搭建正则化加权最小二乘法框架并将总变差正则化作为正则化项，将同步代数迭代重建方法作为代数项，进行逆问题迭代重建计算。

a.给出基于正则化加权最小二乘法框架的求解方程和待优化目标函数，其计算方式表示为：

其中，a^*表示带求解的像素单元的衰减系数分布，

表示满足目标函数取最小值时a的取值。

表示W-范数的平方，W表示预先确定的加权范数矩阵，μ表示正则化参数，

表示二范数的平方，Q(a)表示正则化项。

b.进行非负约束矩阵P的计算，其计算方式表示为；

c.根据步骤二中的预处理结果，进行基于总变差正则化加权最小二乘框架下同步代数重建算法迭代的计算，其具体计算方式可以表示为：

其中，L_β(a^(k))＝D_x ^TLD_x+D_y ^TLD_y，L表示一个N×N维，对角线元素为

的对角矩阵，D_x和D_y分别表示将

和

两个列向量转换为N×N维对角矩阵的形式。

(5)重复(4)直至残差满足要求

其中，Rea^(k)＝||R·a^(k)-τ||，ε为人为设定的残差阈值。。

以上所述实施例为本发明的几个示例模型，本发明不局限于该实施例和附图所公开的内容。凡是不脱离本发明所公开的精神下完成的等效或修改，都在本发明保护的范围。

Claims (2)

1.一种总变差正则化约束的超声成像同步代数迭代重建方法，用于超声层析成像，包含以下步骤：

步骤一：在被测场域外表面均匀的布置一定数量的超声换能器，并使用连续正弦波电压对发射探头进行激励，对除发射外的其他超声探头，记录其接收到正弦信号连续多个单峰值的平均值，记为边界电压测量值，从而获取重建所需的投影衰减测量值τ；

步骤二：构建系数矩阵R，在考虑投影路径与场域像素重叠面积比例的同时考虑重叠区域相对于场域像素的几何位置，计算公式为：

其中，R_i,j是场域内第i条投影路径穿过场域内第j个像素的相对长度，同时对应系数矩阵中第i行、第j列的元素，d_diag表示像素单元中心相对于投影面积几何中心在对角线方向上的距离，l_pixel为像素对角线长度，若第i条投影路径与第j个像素无投影重叠，则R_ij＝0，area_overlap表示重叠区域面积，area_pixel表示像素单元面积；

步骤三：使用同步代数重建方法进行成像预处理：

[1]使用总变差正则化方法计算迭代重建正则化约束项Q(a)，其计算方式可以表示为：

其中，a表示场域内像素单元的衰减系数，

表示场域内第j个像素在横向方向上的灰度变化特征，

表示场域内第j个像素在纵向方向上的灰度变化特征，β表示为防止像素值为零时的不可微现象而加入的一个较小数，Δp_j表示场域内第j个像素单元的面积，N代表所有像素单元的个数；

[2]计算同步代数迭代重建过程中的系数矩阵R，具体计算过程中可以表示为：

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

其中，α表示迭代步长，k表示迭代次数，a^(k)表示第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列的全部元素之和，R_1,+表示系数矩阵R中第1行的全部元素之和，N表示场域内像素个数，M表示场域内投影路径数量，diag( )表示对角矩阵，( )^T表示矩阵转置；

[3]列出预优牛顿最速下降法的基本框架，具体计算方法表示为：

其中，P为非负约束矩阵，D为预优牛顿最速下降法中的权值矩阵，

表示第k次迭代中求解目标函数的梯度，H(Ψ(a^(k)))表示第k次迭代中目标函数的二阶偏微分矩阵，即海瑟矩阵，Ψ(a^(k))表示第k次迭代的目标矩阵；

步骤四：搭建正则化加权最小二乘法框架并将总变差正则化作为正则化项，将同步代数迭代重建方法作为代数项，进行逆问题迭代重建计算；

[1]给出基于正则化加权最小二乘法框架的求解方程和待优化目标函数Ψ(a)：

其中，a^*表示带求解的像素单元的衰减系数分布，

表示满足目标函数取最小值时a的取值；

表示W-范数的平方，W表示预先确定的加权范数矩阵，μ表示正则化参数，

表示二范数的平方，Q(a)表示正则化项；

[2]进行非负约束矩阵P的计算，其计算方式表示为；

其中，a_j表示第j个像素单元的重建衰减系数值；

[3]根据步骤二中的预处理结果，进行基于总变差正则化加权最小二乘框架下同步代数重建算法迭代计算，其具体计算方式可以表示为：

其中，L_β(a^(k))＝D_x ^TLD_x+D_y ^TLD_y，L表示一个N×N维，对角线元素为

的对角矩阵，D_x和D_y分别表示将

和

两个列向量转换为N×N维对角矩阵的形式；

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

其中，Rea^(k)＝||R·a^(k)-τ||，ε为人为设定的残差阈值。

2.根据权利要求1所述的超声成像同步代数迭代重建方法，其特征在于：所述步骤一中边界测量值的获取，是指将被测对象置于超声层析成像被测场域中，场域外均匀分布多个超声换能器，对发射探头进行连续波正弦信号电压激励，采用电压激励电压测量的基本测量模式，循环激励、一发全收的测量策略，投影衰减测量值为不含内含物的空场边界测量电压值和含有内含物的有物场的边界测量电压值的比值。

Images