CN109102552B - 一种非均匀形状约束的像素值域滤波超声成像重建方法 - Google Patents

Abstract

本发明涉及一种基于非均匀形状约束的像素值域滤波超声成像重建方法，包括：根据被测场域，获取重建所需的投影衰减测量值；构建基于路径长度的系数矩阵，即在不含内含物的空场内根据激励探头和接收探头的相对位置计算其连线即投影路径穿过场域内像素的情况；进行成像迭代计算：针对场域内每个像素计算基于非均匀滤波核的像素值滤波模板：将场域内所有像素对应滤波器模板进行卷积，得到滤波后场域内每个位置的像素值；计算滤波后边界测量值和重建估计值之间的残差，迭代直至残差满足要求。

Description

一种非均匀形状约束的像素值域滤波超声成像重建方法

技术领域

本发明属于超声层析成像技术领域，涉及实现一种采用非均匀形状约束的像素值域滤波超声成像方法并用于实现场域内超声层析成像的图像重建。

背景技术

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

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

上述超声成像重建算法及其改进方法中，探头数量的多少对重建图像精度及重建图像分辨率的有至关重要的影响，即超声成像图像重建的优劣与探头间的有效投影路径数目紧密相关：投影路径越多，成像精度越高，伪影越少。但在UPT实际应用过程中，受限于场域尺寸及信号激励幅值限制，场域边界的探头数目不能无限增加；另一方面，超声作为机械波在场域内的传播需要一定的渡越时间，换能器数目过多无法满足可视化监测的实时性要求。在UPT应用于实际生产过程中时，超声换能器的数目一般不超过32个。较高精度的超声成像要求和较快的数据成像速度需要形成了较大的矛盾。因此需要一种在低投影数量下以较高精度和较少伪影进行超声成像逆问题计算的图像重建算法。

发明内容

本发明针对超声过程层析成像逆问题图像重建中，传统成像方法不能同时满足较高测量速度(低投影数量)和较高重建精度(高投影数量)的问题，提出一种基于非均匀滤波核的像素值域滤波超声逆问题图像重建算法。该算法可以在重建结果中保留较为清晰和准确的内含物形状结构，在保证实时成像的基础上显著提升UPT的成像精度。技术方案如下：

一种基于非均匀形状约束的像素值域滤波超声成像重建方法，包括如下步骤：

步骤一：根据被测场域，获取重建所需的第i条投影路径上衰减测量值τ_i，计算方式为

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

步骤二：构建基于路径长度的系数矩阵，即在不含内含物的空场内根据激励探头和接收探头的相对位置计算其连线即投影路径穿过场域内像素的情况，计算公式为：

式中R_i,j是场域内第i条投影路径穿过场域内第j个像素的相对长度，同时对应系数矩阵中第i行、第j列的元素，l_ij为第i条投影路径穿过第j个像素的长度，l_pixel为像素对角线长度，若第i条投影路径未经过第j个像素，则R_ij＝0；

步骤三：使用同步代数重建方法进行成像迭代计算：

[1]给出上一次迭代结果得到的像素值分布a_j ^(k-1)，其中k表示当前迭代次数；

[2]根据系数矩阵R计算投影衰减估计值

n表示场域内所有像素的个数；

[3]计算投影衰减测量值与投影衰减估计值的偏差

[4]计算需要更新的像素值分布

其中

m表示场域内总的投影路径数；

步骤四：根据步骤二的计算结果，针对场域内每个像素计算基于非均匀滤波核的像素值滤波模板：

[1]对全场域像素进行高斯滤波权值计算

[2]对全场域像素进行基于像素值的滤波权值计算

[3]对全场域像素进行非均匀滤波核形状约束权值计算

[4]根据上述计算的权值进行整体滤波器模板设计

式中χ表示滤波中目标像素的位置，ξ表示滤波中场域任一像素的位置，f(χ)表示目标像素的像素值，f(ξ)表示场域内任一像素的像素值，σ_g表示高斯滤波中的位置约束系数，σ_b表示像素值滤波中的像素值约束系数，s表示基于灵敏度先验拟合得到的幂指数因子，σ_s为非均匀滤波核形状约束中的松弛因子，N表示长约内所有像素的集合，sgn表示符号函数，||·||₂表示元素的二范数；

步骤五：将场域内所有像素对应滤波器模板进行卷积，得到滤波后场域内每个位置的像素值，在对目标像素进行滤波计算时，所采用的滤波器窗口包含场域内所有像素，为全尺寸模板；

步骤六：计算滤波后边界测量值和重建估计值之间的残差；

步骤七：重复步骤二～步骤四直至残差满足要求。

本发明在给出透射模态下超声衰减系数重建逆问题模型的基础上，提出一种点到点的系数矩阵构建方式，减小非有效投影路径上像素参与逆问题计算的规模；以基于像素空间域和值域的双边滤波为基础，提取场域尺寸和探头数目先验下表征成像灵敏度信息的特征向量，设计双边滤波核的非均匀约束项，并利用基于同步代数重建方法的逆问题计算框架，实现低投影数量下较高精度的超声衰减图像重建。所提出基于非均匀双边滤波的超声衰减成像算法，其核心思想是“在成像迭代中去除成像伪影噪声并保留内含物边界高频信息和成像中各位置灵敏度均一化表征”其中：通过双边滤波及其全尺寸滤波模板，实现了对重建图像的保边去噪；通过基于灵敏度分布表征的滤波器模板非均匀形状约束，实现了成像中各点位置灵敏度的均一化表征，在低投影数量下的成像结果中保留较完整准确的内含物边界，在给出内含物位置准确位置的同时有效减少了成像伪影，显著提高了UPT图像的重建质量。

附图说明

图1为本发明的基于非均匀形状约束的像素值域滤波算法完整流程图，主要分为同步代数重建计算与非均匀形状约束滤波器模板设计两部分；

图2为本发明中基于投影路径长度的系数矩阵构建方法示意图；

图3为本发明所使用的超声过程层析成像(UPT)系统结构图；

图4为本发明的四个典型仿真模型，并分别给出了相应的传统Tikhonov先验成像结果、滤波反投影(LBP)成像结果、同步代数重建(SART)成像结果及本发明算法(BF-SART)的最终成像结果。

具体实施方式

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

本发明的基于非均匀形状约束的像素值域滤波超声成像算法，实施例中针对工业输油管道中油水两相流的成像这一UPT技术的常见应用形式，使用基于投影路径长度的系数矩阵构建方法表征油水两相流超声可视化检测的正问题模型，同时将图像重建逆问题的迭代求解过程分解为同步迭代重建和非均匀形状约束的像素值域滤波，同步迭代重建部分给出内含物的准确位置和大致轮廓，非均匀形状约束像素值域滤波部分负责去除计算结果的噪声伪影、给出内含物轮廓的准确描述，进行灵敏度均一化约束，以提高低投影数量下对内含物位置、轮廓的准确重建。

如图1所示，为本发明的基于非均匀形状约束的像素值域滤波超声成像算法完整流程图。算法主要分为基于投影路径长度的系数矩阵构建、同步迭代重建与非均匀形状约束的像素值域滤波三部分，系数矩阵在已知场域分布、探头尺寸和探头布置位置的基础上，由数值计算方法获得，图2是其构建方式的基本示意。

图3为超声过程层析成像系统的基本原理图示意，在对油水两相流进行测量时，共计16个超声换能器均匀的沿管壁安装负责激励、接收超声波。采用循环激励、一发全收的测量模式，探头按顺时针方向均匀分布。16个超声探头按顺序接入峰峰值50V、频率1MHz的方波电压激励，探头通道切换时间间隔2.5ms。四与此同时，16个通道同步接收稳态时刻的电压正弦信号持续1ms，并通过正交解调得到接收电压有效值。每次测量总计获得16×15＝240个边界电压测量数据。图4中分别给出了传统UPT成像算法的成像结果与本算法的成像结果。本算法实施例包括如下具体步骤：

(1)系数灵敏度矩阵构建：在不含内含物的空场内根据激励探头和接收探头的相对位置计算其连线(投影路径)穿过场域内像素的情况，其具体计算公式为：

式中R_i,j是场域内第i条投影路径穿过场域内第j个像素的相对长度，若第i条投影路径未经过第j个像素，则R_ij＝0；

(2)获得空场下的边界测量值：在管道中充满背景介质(纯水)的情况下，采用16个超声探头均匀安装在管壁四周，在循环激励、一发全收的测量模式下，获得空场下的240个边界电压测量数据，记为A_s。

(3)针对图4中模型1-模型4的内含物分布情况，分别获取各自重建所需的边界电压测量数据记为A_r，使用空场下边界电压测量数据对内含物仿真模型的边界电压测量数据进行处理，公式表示如下

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

(4)根据噪声阈值选取原则，对处理后的边界测量值向量τ中小于50mV的元素进行剔除，并在删除系数矩阵中对应行的所有数据。

(5)采用线性反投影方法，获得衰减系数分布的初始估计：

a_j ⁽⁰⁾＝R^T·τ

(6)采用同步代数重建计算方法，计算单次迭代结果：

a^(k+1)＝a^(k)+αS_p(S_rR)^T(τ-Ra(^k))

其中，S_p＝diag(1/R_+,1,1/R_+,2,…,1/R_+,N)，S_r＝diag(1/R_1,+,1/R_2,+,…,1/R_M,+)，其中α表示松弛因子，它的计算方式由经验公式给出：

(7)遍历同步代数重建方法所得结果中的每个像素值，计算其相应位置对应的滤波模板：

式中χ表示滤波中目标像素的位置，ξ表示滤波中场域任一像素的位置，f(χ)表示目标像素的像素值，f(ξ)表示场域内任一像素的像素值，σ_g表示高斯滤波中的位置约束系数，σ_b表示像素值滤波中的像素值约束系数，s表示基于灵敏度先验拟合得到的幂指数因子，σ_s为非均匀滤波核形状约束中的松弛因子，N表示长约内所有像素的集合，sgn表示符号函数，||·||₂表示元素的二范数。

(8)对测量区域内每个像素，将其对应的滤波模板与同步代数重建得到的衰减值向量进行卷积运算，得到相应位置衰减系数的滤波值：

式中

表示经第k次滤波后像素χ的取值，

表示经k次同步代数迭代计算后像素χ的取值，p(χ,ξ)为场域内任一像素ξ对目标像素χ的滤波权值。

(9)根据滤波后的衰减系数向量，计算理论边界测量值向量，计算其与实际边界测量值的残差：

Rea_g,k＝||R·a_g,k-τ||

(10)判断前后两次迭代的残差相对变化值是否小于误差允许范围或者达到迭代次数，如果是，算法结束输出成像结果，否则回到步骤(6)。

图4中分别给出了模型1-模型4对应的成像结果。可以看出，传统UPT成像算法结果只能大概地反映出内含物的相对位置，无法较好的提供边界和轮廓信息。本发明所提出的算法能够较准确地重建出内含物的准确位置和完整形状。重建结果中，内含物边界清晰，图像无过多伪影及噪声，图像分辨率、成像精度均有明显提高。

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

Claims (1)

1.一种基于非均匀形状约束的像素值域滤波超声成像重建方法，包括如下步骤：

步骤一：根据被测场域，获取重建所需的第i条投影路径上衰减测量值τ_i，计算方式为

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

步骤二：构建基于路径长度的系数矩阵，即在不含内含物的空场内根据激励探头和接收探头的相对位置计算其连线即投影路径穿过场域内像素的情况，计算公式为：

式中R_i,j是场域内第i条投影路径穿过场域内第j个像素的相对长度，同时对应系数矩阵中第i行、第j列的元素，l_ij为第i条投影路径穿过第j个像素的长度，l_pixel为像素对角线长度，若第i条投影路径未经过第j个像素，则R_ij＝0；

步骤三：使用同步代数重建方法进行成像迭代计算：

[1]给出上一次迭代结果得到的像素值分布a_j ^(k-1)，其中k表示当前迭代次数；

[2]根据系数矩阵R计算投影衰减估计值

n表示场域内所有像素的个数；

[3]计算投影衰减测量值与投影衰减估计值的偏差

[4]计算需要更新的像素值分布

其中

m表示场域内总的投影路径数；

步骤四：根据步骤二的计算结果，针对场域内每个像素计算基于非均匀滤波核的像素值滤波模板：

[1]对全场域像素进行高斯滤波权值计算

[2]对全场域像素进行基于像素值的滤波权值计算

[3]对全场域像素进行非均匀滤波核形状约束权值计算

[4]根据上述计算的权值进行整体滤波器模板设计

式中χ表示滤波中目标像素的位置，ξ表示滤波中场域任一像素的位置，f(χ)表示目标像素的像素值，f(ξ)表示场域内任一像素的像素值，σ_g表示高斯滤波中的位置约束系数，σ_b表示像素值滤波中的像素值约束系数，s表示基于灵敏度先验拟合得到的幂指数因子，σ_s为非均匀滤波核形状约束中的松弛因子，N表示长约内所有像素的集合，sgn表示符号函数，||·||₂表示元素的二范数；

步骤五：将场域内所有像素对应滤波器模板进行卷积，得到滤波后场域内每个位置的像素值，在对目标像素进行滤波计算时，所采用的滤波器窗口包含场域内所有像素，为全尺寸模板；

步骤六：计算滤波后边界测量值和重建估计值之间的残差；

步骤七：重复步骤二～步骤四直至残差满足要求。

Images