修正路径追踪描述的连续波超声层析成像重建方法与流程

本发明属于超声层析成像技术领域，涉及实现一种采用修正路径追踪方法描述的连续波超声成像重建方法，用于实现场域内低对比度两相介质的重建。

背景技术：

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

完整的ut系统主要包含三个部分：超声换能器设计、制造与安装；信号激励与采集系统；超声成像图像重建算法。其中超声成像算法通过对从采集系统得到的换能器接收信号进行处理，通过解调提取测量幅值或渡越时间，得到某个确定激励下的全部换能器的有效测量数据，进一步通过图像重建方法得到场域内含物介质分布的合理估计。目前，超声成像算法作为一种主要利用硬场特性的成像方法，超声成像方法严重依赖于场域边界换能器的数量，其逆问题求解具有严重的病态性(对测量值得微小扰动会导致重建结果的大幅度变化)和欠定性(所需求解的方程数远小于未知量的数目，方程有无穷多解)。为克服这个问题，学者们提出了许多图像重建算法，其中，基于路径的投影重建算法是一种克服病态性的有效手段。该方法通过计算激励、接收换能器间的路径，将收发探头间的时间延迟或幅值衰减均匀的分配给所计算路径上的每一个像素，通过对不同收发探头间的路径进行计算并对同一像素在不同路径上的估计值进行叠加，得到场域内每个像素值的有效估计，以达到可视化测量和图像重建的目的。典型的投影重建方法有徐立军等人1998年在《仪器仪表学报(chinesejournalofscientificinstrument)》第17卷，1-7页，发表的“气液两相泡状流体监测用超声层析成像系统的研究(investigationofultrasoundtomographysystemusedformonitoringbubblygas/liquidtwo-phasefluid)”中提到的二值反投影方法、rahim等人在《sensorsandactuators(传感器与执行器)》第135卷，337-345页发表的“non-invasiveimagingofliquid/gasflowusingultrasonictransmissionmodetomography(超声对液体、气体的非侵入性成像)”中提到的采用r-l函数的线性滤波反投影方法、gordon等人在《journaloftheoreticalbiology(理论生物学期刊)》第29卷，第3期，471-481页发表的“algebraicreconstructiontechniques(art)forthree-dimensionalelectronmicroscopyandx-rayphotography(用于三维电子显微镜和x射线ct的代数重建技术)”中提出的代数重建方法、苏邦良等人在《chemicalengineeringjournal(化学工程期刊)》第77卷，37-41页发表的“theuseofsimultaneousiterativereconstructiontechniqueforelectricalcapacitancetomography(同步迭代重建技术在电容层析成像中的应用)”中提出的同步迭代重建方法、anderson等人在《ultrasoundimaging(超声成像)》第6卷，81-94页发表的“simultaneousalgebraicreconstructiontechnique(sart):asuperiorimplementationoftheartalgorithm(同步代数重建技术(sart)：art算法更优越的实现)”中提出的同步代数重建方法等。其中，sart算法以其收敛快速、残差较小的优点被广泛使用在超声成像的研究中。

不同的应用对象对超声成像图像重建方法提出了不同的要求，面向工业两相流动过程尤其是以油水两相流为代表的弱声阻抗比介质分布，超声成像重建算法需要满足实时性的要求(即成像速度每秒大于24幅)。此外，针对弱声阻抗比介质，除基于已有线性成像模型设计高精度成像算法外，还应对其正问题中的非线性现象进行有效刻画。在满足实时性要求方面：鲍勇等人在《measurementscienceandtechnology(测量科学和技术)》第28卷，074002页发表的“real-timetemperaturefieldmeasurementbasedonacoustictomography(基于超声层析成像的实时温度场测量)”中提出了多探头同步激励方法；在正问题的非线性现象描述方面，anders等人在《ultrasoundimaging(超声成像)》第12卷，268-291页发表的“aray-tracingapproachtorestorationandresolutionenhancementinexperimentalultrasoundtomography(实验超声断层扫描中恢复和分辨率增强的射线追踪方法)”中推导了基于程函方程的超声传播路径追踪模型，瞿晓磊等人在《japanesejournalofappliedphysics(日本应用物理学报)》第56卷，07jf14页发表的“computationalcostreductionbyavoidingray-linkingiterationinbent-raymethodforsoundspeedimagereconstructioninultrasoundcomputedtomography(使用超声成像进行声速图像重建时通过避免在弯曲射线方法中进行射线链接迭代来降低计算成本)”中提出了虚拟接收换能器的概念以提升路径追踪的计算时间，贾佳斌等人在《2017ieeeinternationalconferenceonimagingsystemandtechniques(2017ieee国际成像系统与技术大会)》发表的“nonlineartemperaturefieldreconstructionusingacoustictomography(超声层析成像非线性温度场重建)”中应用非线性路径追踪方法进行温度场分布的有效重建。

在超声成像图像重建算法中，基于路径追踪的正问题非线性描述方法与超声成像实时性的要求间存在矛盾。首先，基于路径追踪的正问题非线性描述方法的目的在于优化确定发射探头与接收探头的最优传播路径，其本质为迭代优化过程，若所有发射与接收探头间传播路径均采用路径追踪方法进行计算，则重建时间会大幅提高，无法满足实时性成像的要求。其次，场域内用于成像的像素数量的多少会直接影响重建算法的计算量及计算时间，随着场域内像素数量提高，重建算法对计算资源的占用较大，在迭代重建算法中，过多的像素无法满足超声成像的实时性要求。

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

技术实现要素：

本发明针对超声层析成像逆问题图像重建，提出一种修正路径追踪描述的连续波超声成像重建方法，可以获得高精度、高信噪比的边界测量数据，提高重建结果的求解精度和图像分辨率。技术方案如下：

一种修正路径追踪描述的连续波超声层析成像重建方法，包括如下步骤：

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

其中，fc是激励信号的中心频率，as为背景介质(水)下的边界电压测量值，ar为存在内含物(水中的离散油泡)情况下的边界电压测量值，ln表示对数符号。

步骤二：基于收发探头相对几何位置进行直线线性投影，构建系数矩阵r，其元素表示为ri,j，表示场域内第i条投影路径穿过场域内第j个像素的相对长度，同时对应系数矩阵中第i行、第j列的元素。

步骤三：使用同步代数重建方法(sart)进行成像迭代计算求得场域内衰减系数的分布，具体计算方式为：

a^(k+1)＝a^(k)+αdp(drr)^t(τ-ra^(k))

其中，α表示迭代步长，k表示迭代次数，a^(k)表示第k次迭代时的像素单元衰减系数分布，dp＝diag(1/r+,1,1/r+,2,····,1/r+,n)，dr＝diag(1/r1,+,1/r2,+,····,1/rm,+)，r+,1表示对系数矩阵内第一列所有元素进行求和，r1,+表示对系数矩阵内第一行所有元素进行求和，n表示场域内像素个数，m表示场域内投影路径数量，diag()表示对角矩阵，()^t表示矩阵转置。

步骤四：使用高斯滤波函数对步骤三的结果进行图像滤波。

步骤五：使用直方图均衡化映射对步骤四的结果进行图像增强。

步骤六：使用差分曲率方法对步骤五的结果进行内含物轮廓提取，该方法通过计算能量熵函数判断图像各像素点是否对应内含物轮廓，在内含物轮廓对应像素点，能量熵函数较大；在非内含物轮廓对应像素点，能量熵函数约等于零。能量熵函数计算方式表示为：

其中，ax表示成像结果水平方向的一阶导数，ay表示成像结果垂直方向的一阶导数，axx表示对ax求水平方向的偏导数，ayy表示对ay求垂直方向的偏导数，axy表示对ax求垂直方向的偏导数，|·|表示绝对值。

步骤七：根据步骤六中提取的内含物轮廓求取声速分布，轮廓内部声速为内含物(油)的声速，轮廓外部声速为背景介质(水)的声速。

步骤八：根据步骤七中得到的声速分布，对给定发射超声探头和接收超声探头进行路径追踪计算。路径追踪计算过程的输入变量为起始点坐标和起始角，起始点坐标为发射超声探头坐标，起始角为发射与接收探头间连线的角度。路径追踪的输出为该条路径上的坐标点向量和路径追踪与场域边界的交点坐标(终止点)。路径追踪的具体计算方式为：

其中，表示下一个像素路径点的方向矢量，表示当前像素路径点的路径矢量，表示路径方向矢量对当前像素求一阶偏导数，表示路径方向矢量对当前像素求二阶偏导数，o((δs)³)表示以像素尺寸为自变量的高阶无穷小，s表示路径追踪计算当前所在的像素点，δs表示像素点的尺寸，n表示场域内折射系数分布，表示场域内折射系数分布的梯度分布，≈为约等于号。

步骤九：根据步骤八中的路径追踪终止点坐标和不同路径的坐标点向量构建新的系数矩阵投影衰减测量值τ和系数矩阵r。投影衰减测量值τ的计算方法为根据路径追踪终止点坐标的位置对边界测量值进行线性插值，系数矩阵r的计算方法为计算不同路径的坐标点向量构成的近似曲线同场域内不同像素相交的长度。

步骤十：重复步骤三至步骤九直至残差满足要求其中，rea^(k)＝||r·a^(k)-τ||，ε为人为设定的残差阈值，m为步骤三至步骤九迭代过程中的迭代次数。

本发明提供的修正路径追踪描述的连续波超声成像重建方法。所提出修正路径追踪描述的连续波超声成像重建方法，实质性特点是：采用路径追踪方法精确表征超声层析成像正问题中的非线性特性，采用低精度迭代、高精度可视化的策略进行逆问题迭代重建。在正问题描述方面,使用连续波超声激励以获得弱声阻抗比两相介质准确的衰减信息，通过衰减分布的成像结果映射声速分布以进行路径追踪计算；提出使用虚拟接收换能器进行路径追踪计算以大幅减小路径追踪计算时间；在逆问题重建方面，采用低分辨率成像迭代计算以满足算法实时性要求，通过使用基于周期性协方差矩阵的高斯过程回归方法在保证成像精度的同时实现高分辨率图像重建。所提出的超声成像重建方法可以获得高精度、高信噪比的边界测量数据，同时给出了超声传播过程中非线性特性的准确描述，在满足实时性要求的情况下，成像结果可以有效解决传统反投影类算法在进行油水两相分布逆问题求解时成像精度差、成像分辨率低、成像伪影严重的问题。所提算法有效扩展了超声透射成像算法的应用，在满足工业流动过程实时性成像要求的基础上，显著提高重建结果的求解精度和图像分辨率。

附图说明

图1为修正路径追踪描述的连续波超声层析成像重建方法基本流程图；

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

图3为本发明中提出的利用虚拟接收探头进行路径追踪计算的示意图；

图4为本发明的四个典型仿真模型，并分别给出了相应的总变差正则化(tv)；

图5给出传统算法与本算法的重建指标对比。

具体实施方式

结合附图和实施例对本发明的修正路径追踪描述的连续波超声层析成像重建方法加以说明。

本发明的修正路径追踪描述的连续波超声层析成像重建方法，实施例中针对工业输油管道中油水两相流的成像这一upt技术的常见应用形式使用本发明所提算法进行计算。在正问题描述方面，使用连续波超声激励以获得弱声阻抗比两相介质准确的衰减信息，通过衰减分布的成像结果映射声速分布以进行路径追踪计算。在逆问题重建计算方面，提出使用虚拟接收换能器进行路径追踪计算以大幅减小路径追踪计算时间，与此同时，采用低分辨率成像迭代计算以满足算法实时性要求，通过使用基于周期性协方差矩阵的高斯过程回归方法在保证成像精度的同时实现高分辨率图像重建。所提出的超声成像重建方法可以获得高精度、高信噪比的边界测量数据，同时给出了超声传播过程中非线性特性的准确描述，在满足实时性要求的情况下，成像结果可以有效解决传统反投影类算法在进行油水两相分布逆问题求解时成像精度差、成像分辨率低、成像伪影严重的问题。所提算法有效扩展了超声透射成像算法的应用，在满足工业流动过程实时性成像要求的基础上，显著提高重建结果的求解精度和图像分辨率。

图1为本发明所提算法的基本流程图，主要分为sart迭代成像、声衰减分布映射声速分布、路径追踪三部分

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

图3为使用虚拟接收探头的路径追踪方法示意，路径追踪起始点为发射超声探头位置坐标，路径追踪起始角为发射探头与接收探头连线角度，将路径追踪终止点作为虚拟接收探头，根据其位置坐标进行边界测量值线性插值，得到每条路径追踪路线对应的测试值。

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

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

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

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

其中，fc是激励信号的中心频率，as为单一背景介质下(水)的边界电压测量值，ar为存在内含物(为离散油泡)情况下的边界电压测量值，ln表示对数符号。

(2)基于收发探头相对几何位置进行直线线性投影，构建系数矩阵r，其元素表示为ri,j，表示场域内第i条投影路径穿过场域内第j个像素的相对长度，同时对应系数矩阵中第i行、第j列的元素。

(3)：使用同步代数重建方法(sart)进行成像迭代计算求得场域内衰减系数的分布，具体计算方式为：

a^(k+1)＝a^(k)+αdp(drr)^t(τ-ra^(k))

其中，α表示迭代步长，采用单因素变量法进行优化，最优取值为1.75，k表示迭代次数，a^(k)表示第k次迭代时的像素单元衰减系数分布，dp＝diag(1/r+,1,1/r+,2,····,1/r+,n)，dr＝diag(1/r1,+,1/r2,+,····,1/rm,+)，r+,1表示对系数矩阵内第一列所有元素进行求和，r1,+表示对系数矩阵内第一行所有元素进行求和，n表示场域内像素个数，m表示场域内投影路径数量，diag()表示对角矩阵，()^t表示矩阵转置。

(4)：使用高斯滤波函数对(3)的结果进行图像滤波。

(5)：使用直方图均衡化映射对(4)的结果进行图像增强。

(6)：使用差分曲率方法对(5)的结果进行内含物轮廓提取，该方法通过计算能量熵函数判断图像各像素点是否对应内含物轮廓，在内含物轮廓对应像素点，能量熵函数较大；在非内含物轮廓对应像素点，能量熵函数约等于零。能量熵函数计算方式表示为：

其中，ax表示成像结果水平方向的一阶导数，ay表示成像结果垂直方向的一阶导数，axx表示对ax求水平方向的偏导数，ayy表示对ay求垂直方向的偏导数，axy表示对ax求垂直方向的偏导数，|·|表示绝对值。

(7)：根据(6)中提取的内含物轮廓求取声速分布，轮廓内部声速为内含物(油)的声速，轮廓外部声速为背景介质(水)的声速。

(8)：根据(7)中得到的声速分布，对给定发射超声探头和接收超声探头进行路径追踪计算。路径追踪计算过程的输入变量为起始点坐标和起始角，起始点坐标为发射超声探头坐标，起始角为发射与接收探头间连线的角度。路径追踪的输出为该条路径上的坐标点向量和路径追踪与场域边界的交点坐标(终止点)。路径追踪的具体计算方式为：

其中，表示下一个像素路径点的方向矢量，表示当前像素路径点的路径矢量，表示路径方向矢量对当前像素求一阶偏导数，表示路径方向矢量对当前像素求二阶偏导数，o((δs)³)表示以像素尺寸为自变量的高阶无穷小，s表示路径追踪计算当前所在的像素点，δs表示像素点的尺寸，n表示场域内折射系数分布，表示场域内折射系数分布的梯度分布，≈为约等于号。

(9)：根据(8)中的路径追踪终止点坐标和不同路径的坐标点向量构建新的系数矩阵投影衰减测量值τ和系数矩阵r。投影衰减测量值τ的计算方法为根据路径追踪终止点坐标的位置对边界测量值进行线性插值，系数矩阵r的计算方法为计算不同路径的坐标点向量构成的近似曲线同场域内不同像素相交的长度。

(10)：重复(3)～(9)直至残差满足要求其中，rea^(k)＝||r·a^(k)-τ||，ε为人为设定的残差阈值，这里选定为0.0001，m为(3)～(9)迭代过程中的迭代次数。

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