CN103860201B - 一种基于宽波束造影成像及提取灌注时间强度曲线的方法 - Google Patents

Abstract

本发明公开了一种基于宽波束造影成像及提取灌注时间强度曲线的方法，通过基于PIWSSD的宽波束造影成像方法以提高CTR，并提供可快速准确提取宽波束造影成像TIC趋势的自适应分析方法，克服宽波束造影成像CTR降低及灌注评价曲线TIC的SCR降低的限制。本发明对于有效降低超声造影成像声功率与造影微泡灌注浓度、降低对人体可能损伤的可能、得到高CTR的造影图像并进行准确的血流灌注评价与诊断具有重要意义。

Description

一种基于宽波束造影成像及提取灌注时间强度曲线的方法

技术领域

本发明涉及超声成像技术领域，特别涉及一种基于宽波束造影成像及提取灌注时间强度曲线的方法。

背景技术

传统超声B模式成像通过逐条扫描线声聚焦实现成像，其帧率可达到60Hz左右。同时，临床为了得到较高分辨率的超声图像和成像景深，通常选用较高的声功率。这不仅可能会对组织造成损伤，也增加了对微泡的破坏，一定程度降低了造影组织比（CTR），而且传统的逐条扫描模式对微泡存在多次破坏，这使得临床为了保持较高的CTR值进一步增强造影剂灌注浓度。而且临床已发现偏高的声功率与造影剂浓度可能对人体存在潜在危害，因此也限制了超声造影成像在临床医学上的进一步应用。同时，在活体组织低帧率B超造影成像中，呼吸运动、脏器不自主蠕动等都会产生运动伪迹，导致前后帧图像具有明显差异，对成像质量和后续时间强度曲线（TIC）的准确提取及灌注的可靠评价造成影响。此外，低帧率对快速运动脏器的捕捉能力差，导致对其评价缺失，对瞬态运动信息的感知不敏感。为此快速高帧率造影成像应需出现。

目前，已广泛应用的超声造影成像技术基本都仅利用微泡振动的线性或者非线性特性，如谐波成像、谐波功率多普勒成像、脉冲逆转成像、脉冲逆转谐波成像等。而高帧率平面波成像技术使微泡变化的瞬时观测及灌注区域的快速运动监测成为可能，并有助于发现新的造影机制。

同时，其中脉冲逆转是一种常用的利用微泡非线性特性与组织线性特征的明显回波差异，获得较高SNR的方法；而微泡小波变换则是在尽可能多的利用微泡回波信息的基础上，将组织的非线性声学特性影响进行抑制，进一步增加CTR值，增强造影微泡的检测灵敏度。解相关技术是基于视频信号下，相邻回波信号之间造影微泡与组织背景的解相关性差异，设定解相关阈值，从而提高CTR的方法。它不仅对组织的运动有较高的灵敏度，而且该技术基于视频信号的处理，弥补了前述脉冲逆转和微泡小波变换基于射频信号的局限性。

通过造影成像得到的散射回波信号包含一系列生理信息，在具体操作中，通常通过提取TIC曲线得到如血流速度信息、血容量时间分布信息、血流动力学等时间分布信息。同时由于微泡引起的回波信号增强与浓度成比例，因而也与组织血容量成比例，这将TIC所表征的时间分布信息转化为空间分布，推动了灌注参量估计的发展。为此TIC趋势的准确快速提取对于参量成像图像的质量以及血流灌注评价和诊断的真实性具有十分重要的意义。

传统TIC曲线提取后有多种处理方法，如基线归零、插值处理、滤波处理、公式或灌注模型拟合等。由于宽波束造影成像具有与传统窄波束条件下不同的特征，这使得传统TIC处理方法，特别是TIC灌注模型拟合的应用受到极大限制。为此本发明采用的DFA（Detrended Fluctuation Analysis）拟合方法，不仅将TIC数据进行分段、兼顾拟合片段前后的数据对重叠片段的影响权重，而且进行自适应非线性多项式拟合。而传统TIC处理技术多基于整体数据拟合和逐点处理，虽能完整提取TIC整体趋势但对于局部动态趋势难以准确拟合；同时，因低帧率下原始TIC数据少，单个点对最终拟合结果的权重较大，后续处理难以全部剔除局部奇异点对整体TIC趋势拟合的影响，而且为保证TIC拟合的准确度，逐点处理的时间成本难以下降。

发明内容

本发明解决的问题在于提供一种基于宽波束造影成像及提取灌注时间强度曲线的方法，通过基于PIWSSD的宽波束造影成像方法以提高CTR，并提供可快速准确提取宽波束造影成像TIC趋势的自适应分析方法，克服宽波束造影成像CTR降低及灌注评价曲线TIC的SCR降低的限制。

为达到上述目的，本发明采用的技术方案是：

一种基于宽波束造影成像的方法，包括以下操作：

1）基于Morgan模型修正的Herring-Trilling微泡振动方程，构建适用于宽波束造影成像下的0/180相微泡母小波；

2）超声阵列探头自发射宽波束射频信号，接收对象反射回来的射频信号；对超声阵列探头得到的0/180相的宽波束射频信号进行基于步骤1）得到的微泡母小波的小波相关性分析，在最佳小波变换尺度下提取最大小波相关系数，构建脉冲逆转微泡小波变换图像矩阵PIW；

3）构造脉冲逆转插值平方和解相关阈值矩阵PISSD，对步骤2）得到的PIW矩阵进行解相关阈值判定，得到宽波束的造影成像图像PIWSSD。

所述的适用于宽波束造影成像下的0/180相微泡母小波的构建为：

分别以宽波束0/180相扫描声压曲线驱动（1）式所示的Herring-Trilling微泡振动方程，得到微泡振动半径随时间的变化曲线及其径向振动速度和加速度，然后带入（2）式中得到宽波束下微泡振动辐声压力波形，再经过带通滤波、归一化处理后作为适用于宽波束造影成像下的0/180相微泡母小波；

$\begin{array}{c}\mathrm{\ρ}\stackrel{\·}{R}\stackrel{\·\·}{R}+\frac{3}{2}\mathrm{\ρ}{\stackrel{\·}{R}}^2=[P_0+\frac{2\mathrm{\σ}}{R_0}+\frac{2\mathrm{\χ}}{R_0}]{\left[\frac{R_0}{R}\right]}^{3\mathrm{\γ}}[1-\frac{3\mathrm{\γ}}{c}\stackrel{\·}{R}]-\frac{4\mathrm{\μ}\stackrel{\·}{R}}{R}-\frac{2\mathrm{\σ}}{R}[1-\frac{1}{c}\stackrel{\·}{R}]\\ -\frac{2\mathrm{\χ}}{R}{\left(\frac{R_0}{R}\right)}^2(1-\frac{3}{c}\stackrel{\·}{R})-12{\mathrm{\μ}}_{\mathrm{sh}}\mathrm{\ϵ}\frac{\stackrel{\·}{R}}{R(R\·\mathrm{\ϵ})}-[P_0+P_{\mathrm{driv}}\left(t\right)]\end{array}---\left(1\right)$

$P={\mathrm{\ρr}}^{-1}(R^2\stackrel{\·\·}{R}+R{\stackrel{\·}{R}}^2)---\left(2\right).$

所述的脉冲逆转微泡小波变换图像矩阵PIW的构建为：

（1）在宽波束射频信号采集完成后，选择合适的尺度，基于步骤1）所构建的0/180相微泡母小波，分别对原始0相和180相射频回波信号进行小波相关性分析，得到一系列小波相关性系数；

（2）对0相和180相回波信号，分别选取最大的小波相关系数所在尺度的所有小波系数作为返回值以表征原始0相和180相射频回波信号；

（3）将0相和180相回波信号中表征原始射频回波信号的小波相关系数相加，进行希尔伯特包络检波后得到脉冲逆转微泡小波变换图像矩阵PIW。

所述的造影成像PIWSSD图像的构建为：

（1）将由步骤2）得到的由0相和180相小波相关系数表征的回波信号划分为长度为L_w的采样窗，采样窗彼此重叠长度L_o；

（2）PISSD阈值矩阵由采样点总数N、采样窗序数n、采样窗长度L_w、采样窗重叠长度L_o、采样窗函数w[]及0相和180相小波相关系数决定，PISSD算法如式3所示：

$\mathrm{PISSD}\left(n\right)=\frac{1}{L_w}\underset{k=1+(n-1)\mathrm{\Δ}}{\overset{L_0+\mathrm{n\Δ}}{\mathrm{\Σ}}}{\left(w\right[s_{f1}\left(k\right)+w\left[s_{f2}\left(k\right)\right]\left]\right)}^2---\left(3\right)$

n=1,2,...,(N-L_w)Δ+1,Δ=L_w-L_o,L_o<L_w

其中s_f1(k)和s_f2(k)分别为步骤3）得到的由0相和180相小波相关系数表征的回波信号；

（3）利用步骤（2）的PISSD算法构建描述组织背景解相关程度的二维PISSD阈值矩阵；

（4）基于PISSD阈值矩阵设定四种PISSD阈值：统一阈值、分段阈值、低幅联合阈值与高幅联合阈值；

所述的统一阈值为PISSD矩阵的均值与标准偏差的和；

所述的分段阈值为沿扫描波束方向将PISSD矩阵分为M个等尺寸的片段，计算每个片段内PISSD数值的均值与标准偏差的和，将该数值作为对应片段的分段阈值；

所述的低幅联合阈值LCT如式4所示：

$\mathrm{LCT}(i,j)=\left\{\begin{array}{cc}\mathrm{DT}(i,j),& \mathrm{DT}(i,j)\&le;\mathrm{UT}\\ \mathrm{UT},& \mathrm{DT}(i,j)>\mathrm{UT}\end{array}\right.---\left(4\right)$

i=1,2,...,N_s,j=1,2,...,N_l

所述的高幅联合阈值HCT如式5所示：

$\mathrm{HCT}(i,j)=\left\{\begin{array}{cc}\mathrm{DT}(i,j)& \mathrm{DT}(i,j)\&GreaterEqual;\mathrm{UT}\\ \mathrm{UT},& \mathrm{DT}(i,j)<\mathrm{UT}\end{array}\right.---\left(5\right)$

i=1,2,...,N_s,j=1,2,...,N_l

（5）基于步骤（4）的四种阈值矩阵，将PISSD阈值矩阵，PISSD图像和PIW图像矩阵输入PISSD阈值判定模块，得到PIHSSD图像。

进一步，本发明提供一种基于所述的宽波束造影成像的提取灌注时间强度曲线的方法，还包括以下操作：

4）提取步骤3）的PIWSSD图像的感兴趣区域的灌注评价曲线TIC；

5）构造适合DFA分析的窗长阶数矩阵X和多项式拟合系数矩阵Fit；

6）对步骤4）获取的TIC曲线分别进行DFA分析提取TIC趋势。

所述的窗长阶数矩阵X和多项式拟合系数矩阵Fit的构建为：

（1）适合DFA分析的窗长阶数矩阵X由原始TIC数据长度Num、DFA分析半窗长N及多项式拟合阶数K决定；

所述适合DFA分析的窗长阶数矩阵X如式（6）所示；

（2）适合DFA分析的多项式拟合系数矩阵Fit由原始TIC数据长度Num和DFA分析半窗长N决定；

所述适合DFA分析的多项式拟合系数矩阵Fit由最小二乘多元回归分析方法获取TIC数据片段的最优多项式拟合系数矩阵P如式7所示；

适合DFA分析的多项式拟合系数矩阵Fit如式8所示；

所述的DFA分析提取TIC趋势包括以下操作：

（1）首先对原始TIC数据进行判断，若行数小于列数，则进行矩阵转秩，使行数大于列数；

（2）将TIC分为2N+1个相同长度的片段，并将片段长度记为窗长，相邻片段之间有N+1个点重叠；

（3）对于每一个片段寻找最优的阶数为K的多项式拟合，当K为0或1的时候，分别对应于常数和线性拟合；

（4）构造适合DFA分析的窗长阶数矩阵X和多项式拟合系数矩阵Fit；

（5）由最小二乘多元回归分析方法获取TIC数据片段的最优多项式拟合系数矩阵P；

（6）判断整段TIC数据能否被窗长整除，若能够整除则球多项式线性拟合的加权系数，若不能则首先对最后一个长度小于窗长的片段进行加权拟合；

（7）将第i个片段和第i+1个片段的经拟合的多项式分别标记为yⁱ(l₁)和yⁱ⁺¹(l₂)，l₁,l₂=1,...,2N+1；

重叠片段的拟合多项式y^o(l)为：y^o(l)=w₁yⁱ(l+N)+w₂yⁱ⁺¹(l),l=1,...,N+1，其中w₁=1-(l-1)/N,w₂=(l-1)/N，为多项式线性拟合的加权系数；

（8）依次对重叠部分数据进行加权拟合；

（9）将最后一个长度小于窗长的片段及其之前的部分的加权拟合结果合并，输出拟合TIC曲线。

还对拟合TIC曲线还进行评价拟合质量：

（1）选择拟合质量QOF作为时间强度曲线整体趋势的评价指标；

（2）选择细节评价指数DEI作为时间强度曲线整体趋势的评价指标；

（3）利用matlab软件编程实现QOF、DEI计算；

（4）分别对TIC数据进行DFA拟合，计算QOF、DEI和运算时间。

所述的拟合质量QOF达到最大值100%时表明TIC的拟合质量最佳，当QOF低于80%时认为TIC的拟合不准确，不具备可信性；

拟合质量QOF如式（9）、（10）、（11）所示：

$\mathrm{QOF}=(1-\frac{{\mathrm{SSR}}_f}{{\mathrm{VAR}}_d})\&times;100\%---\left(9\right)$

${\mathrm{SSR}}_f=\underset{t=t_{\mathrm{start}}}{\overset{t_{\mathrm{end}}}{\mathrm{\&Sigma;}}}{({\mathrm{TIC}}_f\left(t\right)-{\mathrm{TIC}}_d\left(t\right))}^2---\left(10\right)$

${\mathrm{VAR}}_d=\underset{t=t_{\mathrm{start}}}{\overset{t_{\mathrm{end}}}{\mathrm{\&Sigma;}}}{({\mathrm{TIC}}_d\left(t\right)-\frac{\underset{t=t_{\mathrm{start}}}{\overset{t_{\mathrm{end}}}{\mathrm{\&Sigma;}}}{\mathrm{TIC}}_d\left(t\right)}{N})}^2---\left(11\right)$

其中SSR_f为拟合区间内拟合TIC与滤波TIC的残差平方和，VAR_d为滤波TIC在拟合区间内的方差，TIC_f为拟合TIC，TIC_d为滤波TIC，t_start,t_end为拟合区间的起始点和终止点，N为拟合区间内的样本点数。

所述的细节评价指数DEI越大表明拟合后的TIC对原始数据的细节保留能力越强，当DEI小于0.85时认为拟合TIC对灌注细节失真，不具备可信性；

DEI的定义为：

$\mathrm{DEI}=\frac{\mathrm{\&Sigma;}\left(\right|{\mathrm{TIC}}_f\left(t\right)-{\mathrm{TIC}}_f(t-1)|+|{\mathrm{TIC}}_f\left(t\right)-{\mathrm{TIC}}_f(t+1)\left|\right)}{\mathrm{\&Sigma;}({\mathrm{TIC}}_d\left(t\right)-{\mathrm{TIC}}_d(t-1)|+|{\mathrm{TIC}}_d\left(t\right)-{\mathrm{TIC}}_d(t+1)|)}---\left(12\right)$

式中：TIC_f(t)和TIC_d(t)分别为拟合前后的TIC值，DEI最大值为1，最小值为0。

与现有技术相比，本发明具有以下有益的技术效果：

本发明提供的基于宽波束造影成像方法，是将脉冲逆转、微泡小波变换与解相关技术相结合，而提出的针对宽波束造影成像的造影成像算法。该方法结合了脉冲逆转、小波变换技术与解相关算法的优点，可大幅提高宽波束超声造影成像的CTR，可在保证超声造影图像CTR需求的同时，大幅降低超声探头的声功率，在一定程度可降低超声高声功率对人体的潜在损害。

而且与传统模式逐条扫描的方式不同，宽波束是多条扫描线同时发射和接收。同等帧率时，宽波束对微泡的破坏次数减少，作用于微泡的声功率也降低，这使得探头对微泡的破坏程度降低，客服了传统B超高声功率和高造影浓度的缺点。本发明提供的脉冲逆转小波变换解相关宽波束造影成像方法一改传统超声探头逐条扫描对造影微泡多次破坏的缺点，为临床在保证图像质量的同时降低灌注造影微泡浓度提供成像方法。

本发明提供的提取灌注时间强度曲线的方法，是一种更为准确的自适应提取TIC整体趋势及局部灌注趋势的方法，基于数据分段、局部交叉和高阶次拟合的原理，综合QOF（Quality of Fitting）整体趋势和DEI（DetailedEvaluation Index）局部细节的评价指标。

宽波束条件下的TIC曲线虽然点数多，局部奇异点也多，但因DFA具有局部交叉、高阶次拟合的特征，使得在有效剔除奇异点、降低运算时间的同时，可以最大限度提取整体趋势及局部动态趋势，这对于二次循环或灌注和灌注时间延长的TIC更加有利。为此，在基于脉冲逆转小波变换解相关造影成像的宽波束平面波条件下实现快速提取TIC曲线趋势、获取准确的灌注信息对血流灌注评价和诊断具有重要意义。

本发明提供的提取灌注时间强度曲线的方法，是一种更为快速的提取宽波束造影成像-TIC趋势的方法，相比于传统逐点移动窗长进行拟合处理的方法，DFA以较大的半窗长为步长进行移动，使得其运算时间成本降低。

本发明提供的基于宽波束造影成像及提取灌注时间强度曲线的方法，对于有效降低超声造影成像声功率与造影微泡灌注浓度、降低对人体可能损伤的可能、得到高CTR的造影图像并进行准确的血流灌注评价与诊断具有重要意义。

附图说明

图1为脉冲逆转小波变换解相关宽波束造影成像及灌注评价方法框图；

图2为脉冲逆转小波变换解相关宽波束造影成像原理流程图；

图3为准确提取宽波束造影成像-时间强度曲线趋势(DFA)的快速计算流程图；

图4为宽/窄波束激励下造影微泡静态衰减比较图

图5为活体兔肾脉冲逆转小波变换解相关宽波束造影成像结果

图6为管TIC的DFA与中值滤波结果比较及拟合系数与运算时间评价。

具体实施方式

下面结合具体的实施例对本发明做进一步的详细说明，所述是对本发明的解释而不是限定。

针对宽波束造影成像CTR降低及灌注评价曲线TIC的SCR降低等限制，本发明提供的基于宽波束造影成像及提取灌注时间强度曲线的方法，首先通过基于PIWSSD宽波束造影成像方法以提高CTR，然后基于该成像方法进行可快速准确提取宽波束造影成像TIC趋势的分析，这对于有效降低超声造影成像声功率与造影微泡灌注浓度、降低对人体可能损伤的可能、得到高CTR的造影图像并进行准确的血流灌注评价与诊断具有重要意义。

参见图1所示，脉冲逆转小波变换解相关宽波束造影成像，先利用超声成像系统得到射频数据，进行小波分析得到脉冲逆转微泡小波变换（PIW）图像矩阵，在此基础上结合差值平方和解相关（PISSD）技术，从而得到宽波束的PIWSSD图像；

然后，提取PIWSSD图中感兴趣区域（ROI）的灌注评价曲线TIC，准确快速提取该TIC趋势，进行宽波束造影成像灌注评价。

本发明提供的基于宽波束的造影成像方法，包括以下步骤：

1）基于Morgan模型修正的Herring–Trilling微泡振动方程（如公式1、2所示），构建适用于宽波束造影成像下分别适合的压缩/扩张相（0/180相）母小波。

$\begin{array}{c}\mathrm{\&rho;}\stackrel{\&CenterDot;}{R}\stackrel{\&CenterDot;\&CenterDot;}{R}+\frac{3}{2}\mathrm{\&rho;}{\stackrel{\&CenterDot;}{R}}^2=[P_0+\frac{2\mathrm{\&sigma;}}{R_0}+\frac{2\mathrm{\&chi;}}{R_0}]{\left[\frac{R_0}{R}\right]}^{3\mathrm{\&gamma;}}[1-\frac{3\mathrm{\&gamma;}}{c}\stackrel{\&CenterDot;}{R}]-\frac{4\mathrm{\&mu;}\stackrel{\&CenterDot;}{R}}{R}-\frac{2\mathrm{\&sigma;}}{R}[1-\frac{1}{c}\stackrel{\&CenterDot;}{R}]\\ -\frac{2\mathrm{\&chi;}}{R}{\left(\frac{R_0}{R}\right)}^2(1-\frac{3}{c}\stackrel{\&CenterDot;}{R})-12{\mathrm{\&mu;}}_{\mathrm{sh}}\mathrm{\&epsiv;}\frac{\stackrel{\&CenterDot;}{R}}{R(R\&CenterDot;\mathrm{\&epsiv;})}-[P_0+P_{\mathrm{driv}}\left(t\right)]\end{array}---\left(1\right)$

$P={\mathrm{\&rho;r}}^{-1}(R^2\stackrel{\&CenterDot;\&CenterDot;}{R}+R{\stackrel{\&CenterDot;}{R}}^2)---\left(2\right)$

具体操作如下：

分别以宽波束0/180相扫描声压曲线驱动Herring–Trilling微泡振动方程（1式所示），得到微泡振动半径随时间的变化曲线及其径向振动速度和加速度，带入（2）式中得到宽波束下微泡振动辐声压力波形，然后经过带通滤波、归一化处理后作为适用于宽波束造影成像下的0/180相微泡母小波。

上述方程中各参量定义与取值参见表1

表1：Morgan模型中符号的含义与仿真计算时使用的参数

2）超声阵列探头自发射宽波束射频信号，接收对象反射回来的宽波束射频信号；

对超声阵列探头得到的0/180相的宽波束射频信号进行基于步骤1）得到的母小波的小波相关性分析，在最佳小波变换尺度下提取最大小波相关系数，构建脉冲逆转微泡小波变换（PIW）图像矩阵；

参见图2所示，具体步骤为：

（1）在宽波束射频信号采集完成后，根据经验选择合适的尺度，基于步骤1）所构建的合适的0/180相微泡母小波，分别对原始0相和180相射频回波信号进行小波相关性分析，得到一系列小波相关性系数；

（2）对0相和180相回波信号，分别选取最大的小波相关系数所在尺度的所有小波系数作为返回值以表征原始0相和180相射频回波信号；

（3）将0相和180相回波信号中表征原始射频回波信号的小波相关系数相加，进行希尔伯特包络检波后得到PIW图像矩阵。

3）构造脉冲逆转插值平方和解相关（PISSD）阈值矩阵，对步骤2）得到的PIW图像矩阵进行解相关阈值判定，得到宽波束下高CTR的PIWSSD图像；

参见图2，具体操作如下：

（1）将由步骤2）得到的由0相和180相小波相关系数表征的回波信号划分为长度为L_w的采样窗，采样窗彼此重叠长度L_o；

（2）PISSD阈值矩阵由采样点总数N、采样窗序数n、采样窗长度L_w、采样窗重叠长度L_o、采样窗函数w[]及0相和180相小波相关系数决定，PISSD算法如式3所示：

$\mathrm{PISSD}\left(n\right)=\frac{1}{L_w}\underset{k=1+(n-1)\mathrm{\&Delta;}}{\overset{L_0+\mathrm{n\&Delta;}}{\mathrm{\&Sigma;}}}{\left(w\right[s_{f1}\left(k\right)+w\left[s_{f2}\left(k\right)\right]\left]\right)}^2---\left(3\right)$

n=1,2,...,(N-L_w)Δ+1,Δ=L_w-L_o,L_o<L_w

其中s_f1(k)和s_f2(k)分别为步骤3）得到的由0相和180相小波相关系数表征的回波信号；

（3）利用上述PISSD算法构建描述组织背景解相关程度的二维PISSD阈值矩阵。

（4）基于PISSD阈值矩阵设定四种PISSD阈值：统一阈值、分段阈值、低幅联合阈值与高幅联合阈值；

所述的统一阈值（Uniform threshold,UT）为PISSD矩阵的均值与标准偏差的和；

所述的分段阈值（Divided threshold,DT）为沿扫描波束方向将PISSD矩阵分为M个等尺寸的片段，计算每个片段内PISSD数值的均值与标准偏差的和，将该数值作为对应片段的分段阈值。

所述的低幅联合阈值（Lower combined threshold,LCT）如式4所示：

$\mathrm{LCT}(i,j)=\left\{\begin{array}{cc}\mathrm{DT}(i,j),& \mathrm{DT}(i,j)\&le;\mathrm{UT}\\ \mathrm{UT},& \mathrm{DT}(i,j)>\mathrm{UT}\end{array}\right.---\left(4\right)$

i=1,2,...,N_s,j=1,2,...,N_l

所述的高幅联合阈值（Higher combined threshold,HCT）如式5所示：

$\mathrm{HCT}(i,j)=\left\{\begin{array}{cc}\mathrm{DT}(i,j)& \mathrm{DT}(i,j)\&GreaterEqual;\mathrm{UT}\\ \mathrm{UT},& \mathrm{DT}(i,j)<\mathrm{UT}\end{array}\right.---\left(5\right)$

i=1,2,...,N_s,j=1,2,...,N_l

（5）基于上述四种阈值矩阵，将PISSD阈值矩阵，PISSD图像和PIW图像矩阵输入PISSD阈值判定模块，得到PIHSSD图像。

进一步，在基于所得的PIHSSD图像，本发明提出的提取灌注时间强度曲线的方法包括以下操作：

4）提取步骤3）PIWSSD图像的感兴趣区域（ROI）的灌注评价曲线TIC；

5）构造适合DFA（Detrended Fluctuation Analysis）分析的窗长阶数矩阵X和多项式拟合系数矩阵Fit；

其步骤如下：

（1）适合DFA分析的窗长阶数矩阵X由原始TIC数据长度Num、DFA分析半窗长N及多项式拟合阶数K决定；

所述适合DFA分析的窗长阶数矩阵X如式（6）所示；

（2）适合DFA分析的多项式拟合系数矩阵Fit由原始TIC数据长度Num和DFA分析半窗长N决定；

所述适合DFA分析的多项式拟合系数矩阵Fit由最小二乘多元回归分析方法获取TIC数据片段的最优多项式拟合系数矩阵P如式7所示；

适合DFA分析的多项式拟合系数矩阵Fit如式8所示；

6）对步骤4）获取的TIC曲线分别进行DFA分析提取TIC趋势，同时为了便于对所提取的效果进行评价，同时进行中值滤波、多项式拟合以作为对照；

参见图5，所述的DFA（Detrended Fluctuation Analysis）分析，包括以下步骤：

（1）首先对原始TIC数据进行判断，若行数小于列数，则进行矩阵转秩，使行数大于列数；

（2）将TIC分为若干个相同长度（2N+1）（记为窗长）的片段，相邻片段之间有N+1个点重叠；

（3）对于每一个片段寻找最优的阶数为K的多项式拟合，当K为0或1的时候，分别对应于常数和线性拟合；

（4）构造适合DFA分析的窗长阶数矩阵X和多项式拟合系数矩阵Fit（如式4所示）；

（5）由最小二乘多元回归分析方法获取TIC数据片段的最优多项式拟合系数矩阵P（如式5所示）；

（6）判断整段TIC数据能否被窗长整除，若能够整除则球多项式线性拟合的加权系数，若不能则首先对最后一个长度小于窗长的片段进行加权拟合。

（7）将第i个片段和第i+1个片段的经拟合的多项式分别标记为yⁱ(l₁)和yⁱ⁺¹(l₂)，l₁,l₂=1,...,2N+1；

重叠片段的拟合多项式y^o(l)定义为：y^o(l)=w₁yⁱ(l+N)+w₂yⁱ⁺¹(l),l=1,...,N+1，其中w₁=1-(l-1)/N,w₂=(l-1)/N，为多项式线性拟合的加权系数。

（8）依次对重叠部分数据进行加权拟合，该加权方式在具备了对称性的同时，能够较好的去除相邻片段重叠部分（片段边界）的跳跃或者不连续。

（9）将上述最后一个长度小于窗长的片段及其之前的部分的加权拟合结果合并，输出拟合TIC曲线。

对步骤6）得到的TIC定义QOF及DEI，评价拟合质量；在matlab中采用tic与toc方法计算运行时间评价运算速度。

具体步骤为：

（1）选择拟合质量（QOF）作为时间强度曲线整体趋势的评价指标。QOF的定义如式9、10、11所示：

$\mathrm{QOF}=(1-\frac{{\mathrm{SSR}}_f}{{\mathrm{VAR}}_d})\&times;100\%---\left(9\right)$

${\mathrm{SSR}}_f=\underset{t=t_{\mathrm{start}}}{\overset{t_{\mathrm{end}}}{\mathrm{\&Sigma;}}}{({\mathrm{TIC}}_f\left(t\right)-{\mathrm{TIC}}_d\left(t\right))}^2---\left(10\right)$

${\mathrm{VAR}}_d=\underset{t=t_{\mathrm{start}}}{\overset{t_{\mathrm{end}}}{\mathrm{\&Sigma;}}}{({\mathrm{TIC}}_d\left(t\right)-\frac{\underset{t=t_{\mathrm{start}}}{\overset{t_{\mathrm{end}}}{\mathrm{\&Sigma;}}}{\mathrm{TIC}}_d\left(t\right)}{N})}^2---\left(11\right)$

其中SSR_f为拟合区间内拟合TIC与滤波TIC的残差平方和，VAR_d为滤波TIC在拟合区间内的方差，TIC_f为拟合TIC，TIC_d为滤波TIC，t_start,t_end为拟合区间的起始点和终止点，N为拟合区间内的样本点数。QOF是个百分数，它与SSR_f呈反比，SSR_f越小，QOF越大。当SSR_f为零时，QOF达到最大值100%，表明TIC的拟合质量最佳。当QOF低于80%时，认为TIC的拟合不准确，不具备可信性；

（2）选择细节评价指数（DEI）作为时间强度曲线整体趋势的评价指标。DEI的定义如方程12所示：

$\mathrm{DEI}=\frac{\mathrm{\&Sigma;}\left(\right|{\mathrm{TIC}}_f\left(t\right)-{\mathrm{TIC}}_f(t-1)|+|{\mathrm{TIC}}_f\left(t\right)-{\mathrm{TIC}}_f(t+1)\left|\right)}{\mathrm{\&Sigma;}({\mathrm{TIC}}_d\left(t\right)-{\mathrm{TIC}}_d(t-1)|+|{\mathrm{TIC}}_d\left(t\right)-{\mathrm{TIC}}_d(t+1)|)}---\left(12\right)$

式中：TIC_f(t)和TIC_d(t)分别为拟合前后的TIC值。DEI最大值为1，最小值为0。DEI越大，表明拟合后的TIC对原始数据的细节保留能力越强，当DEI小于0.85时，认为拟合TIC对灌注细节失真，不具备可信性。

（3）利用matlab软件编程实现QOF、DEI计算；

（4）分别对TIC数据进行DFA拟合和中值滤波处理，计算QOF、DEI和运算时间。

参见图4所示，在磁力搅拌子的作用下发射中心频率5MHz的超声波作用于600ml稀释1000倍的浓度为1%SonoVue，分别得到逐条扫描的窄波束（图a）和宽波束（图b）条件下的造影图像，并提取两种模式造影图中同样大小ROI的TIC曲线（图c所示）。可以看到两种曲线随时间都呈下降趋势，这是由声束对微泡的破坏及微泡的自由破裂引起的。但相比于宽波束，窄波束曲线的起始点高，下降速度快。这说明窄波束的发射功率更大，而且存在多次破坏，故对微泡的破坏作用更强。

参见图5所示，为体重3.3kg兔子肾脏脉冲逆转小波变换解相关宽波束造影成像结果。其中图（a）、（b）为构建的适用于宽波束造影成像下分别适合的压缩/扩张相（0/180相）母小波，图（c）、（d）、（e）、（f）分别为原始脉冲逆转成像、脉冲逆转微泡小波变换成像、脉冲逆转解相关成像及脉冲逆转小波变换解相关成像。可以看到，四种方法对三个红色圈内的区域都有一定程度的保留，这应该是具有血流信息的肾脏内血管区和肾脏外血管区，其他区域从（c）到（f）逐渐被去掉，这应该是运动信息小的组织区。可以看出从（c）到（f），造影组织比逐渐提高，说明在宽波束条件下脉冲逆转小波变换解相关的成像质量优于前三种方法。

参见图6所示，图（a）为原始TIC数据、十阶多项式拟合曲线、窗长为60的DFA拟合后及中值滤波后曲线对比显示，整体多项式拟合提取只能提取出TIC十分粗糙的趋势，几乎没有展示细节信息；而DFA拟合与中值滤波都能提取出TIC的整体趋势，但DFA的细节展示能力更强，而且中值滤波在重要的波峰和波谷处都产生不真实的平滑数据。从图（b）可以看到，在10到70的窗长范围内，同样大小窗长下DFA的QOF值平均高出中值滤波的6%，且在窗长为10时中值滤波的QOF低于80%，不具备可信性。图（c），在10到70的窗长范围内，DFA的运算时间也比中值滤波平均缩短6.3ms。从这两方面说明DFA在拟合效果和运行速度方面都要优于中值滤波。

Claims (10)

1.一种基于宽波束造影成像的方法，其特征在于，包括以下操作：

1)基于Morgan模型修正的Herring-Trilling微泡振动方程，构建适用于宽波束造影成像下的0/180相微泡母小波；

2)超声阵列探头自发射宽波束射频信号，接收对象反射回来的射频信号；对超声阵列探头得到的0/180相的宽波束射频信号进行基于步骤1)得到的微泡母小波的小波相关性分析，在最佳小波变换尺度下提取最大小波相关系数，构建脉冲逆转微泡小波变换图像矩阵PIW；

3)构造脉冲逆转插值平方和解相关阈值矩阵PISSD，对步骤2)得到的PIW矩阵进行解相关阈值判定，得到宽波束的造影成像图像PIWSSD。

2.如权利要求1所述的基于宽波束造影成像的方法，其特征在于，所述的适用于宽波束造影成像下的0/180相微泡母小波的构建为：

分别以宽波束0/180相扫描声压曲线驱动(1)式所示的Herring-Trilling微泡振动方程，得到微泡振动半径随时间的变化曲线及其径向振动速度和加速度，然后带入(2)式中得到宽波束下微泡振动辐声压力波形，再经过带通滤波、归一化处理后作为适用于宽波束造影成像下的0/180相微泡母小波；

$\begin{array}{c}\mathrm{\&rho;}\stackrel{.}{R}\stackrel{..}{R}\frac{3}{2}\mathrm{\&rho;}{\stackrel{.}{R}}^2=[P_0+\frac{2\mathrm{\&sigma;}}{R_0}+\frac{2\mathrm{\&chi;}}{R_0}]{\left[\frac{R_0}{R}\right]}^{3\mathrm{\&gamma;}}[1-\frac{3\mathrm{\&gamma;}}{c}\stackrel{.}{R}]-\frac{4\mathrm{\&mu;}\stackrel{.}{R}}{R}-\frac{2\mathrm{\&sigma;}}{R}[1-\frac{1}{c}\stackrel{.}{R}]\\ -\frac{2\mathrm{\&chi;}}{R}{\left(\frac{R_0}{R}\right)}^2(1-\frac{3}{c}\stackrel{.}{R})-12{\mathrm{\&mu;}}_{\mathrm{sh}}\mathrm{\&epsiv;}\frac{\stackrel{.}{R}}{R(R\&CenterDot;\mathrm{\&epsiv;})}-[P_0+P_{\mathrm{driv}}\left(t\right)]\end{array}---\left(1\right)$

$P={\mathrm{\&rho;r}}^{-1}(R^2\stackrel{..}{R}+R{\stackrel{.}{R}}^2)---\left(2\right);$

其中，ρ为液体密度，P₀为液体静态压，σ为液体表面张力，χ为薄包膜的弹性模量，R₀为微泡的静态半径，R为瞬态微泡半径，γ为气体的多归系数，c为液体中声速，μ为液体的粘滞系数，ε为薄包膜的厚度，μ_shε为膜厚与材料粘滞系数的积，P_driv(t)为驱动声压。

3.如权利要求1所述的基于宽波束造影成像的方法，其特征在于，所述的脉冲逆转微泡小波变换图像矩阵PIW的构建为：

(1)在宽波束射频信号采集完成后，选择合适的尺度，基于步骤1)所构建的0/180相微泡母小波，分别对原始0相和180相射频回波信号进行小波相关性分析，得到一系列小波相关性系数；

(2)对0相和180相回波信号，分别选取最大的小波相关系数所在尺度的所有小波系数作为返回值以表征原始0相和180相射频回波信号；

(3)将0相和180相回波信号中表征原始射频回波信号的小波相关系数相加，进行希尔伯特包络检波后得到脉冲逆转微泡小波变换图像矩阵PIW。

4.如权利要求1所述的基于宽波束造影成像的方法，其特征在于，所述的造影成像PIWSSD图像的构建为：

(1)将由步骤2)得到的由0相和180相小波相关系数表征的回波信号划分为长度为L_w的采样窗，采样窗彼此重叠长度L_o；

(2)PISSD阈值矩阵由采样点总数N、采样窗序数n、采样窗长度L_w、采样窗重叠长度L_o、采样窗函数w[·]及0相和180相小波相关系数决定，PISSD算法如式3所示：

$\begin{array}{c}\mathrm{PISSD}\left(n\right)=\frac{1}{L_w}\underset{k=1+(n-1)\mathrm{\&Delta;}}{\overset{L_o+\mathrm{n\&Delta;}}{\mathrm{\&Sigma;}}}{\left(w\right[s_{f1}\left(k\right)+w\left[s_{f2}\left(k\right)\right]\left]\right)}^2\\ n=\mathrm{1,2},...,(N-L_w)/\mathrm{\&Delta;}+1,\mathrm{\&Delta;}=L_w-L_o,L_o<L_w\end{array}---\left(3\right)$

其中s_f1(k)和s_f2(k)分别为步骤3)得到的由0相和180相小波相关系数表征的回波信号；

(3)利用步骤(2)的PISSD算法构建描述组织背景解相关程度的二维PISSD阈值矩阵；

(4)基于PISSD阈值矩阵设定四种PISSD阈值：统一阈值、分段阈值、低幅联合阈值与高幅联合阈值；

所述的统一阈值为PISSD矩阵的均值与标准偏差的和；

所述的分段阈值为沿扫描波束方向将PISSD矩阵分为M个等尺寸的片段，计算每个片段内PISSD数值的均值与标准偏差的和，将该数值作为对应片段的分段阈值；

所述的低幅联合阈值LCT如式4所示：

$\begin{array}{c}\mathrm{LCT}(i,j)=\left\{\begin{array}{cc}\mathrm{DT}(i,j),& \mathrm{DT}(i,j)\&le;\mathrm{UT}\\ \mathrm{UT},& \mathrm{DT}(i,j)>\mathrm{UT}\end{array}\right.\\ i=\mathrm{1,2}...,N_s,j=\mathrm{1,2},...,N_l\end{array}---\left(4\right)$

其中，UT表示统一阈值，统一阈值为PISSD矩阵的均值与标准偏差的和；DT表示分段阈值，分段阈值为沿扫描波束方向将PISSD矩阵分为M个等尺寸的片段，计算每个片段内PISSD数值的均值与标准偏差的和，将该数值作为对应片段的分段阈值；

所述的高幅联合阈值HCT如式5所示：

$\begin{array}{c}\mathrm{HCT}(i,j)=\left\{\begin{array}{cc}\mathrm{DT}(i,j),& \mathrm{DT}(i,j)\&GreaterEqual;\mathrm{UT}\\ \mathrm{UT},& \mathrm{DT}(i,j)<\mathrm{UT}\end{array}\right.\\ i=\mathrm{1,2}...,N_s,j=\mathrm{1,2},...,N_l\end{array}---\left(5\right)$

(5)基于步骤(4)的四种阈值矩阵，将PISSD阈值矩阵，PISSD图像和PIW图像矩阵输入PISSD阈值判定模块，得到PIHSSD图像。

5.一种基于权利要求1所述的宽波束造影成像的提取灌注时间强度曲线的方法，其特征在于，还包括以下操作：

4)提取步骤3)的PIWSSD图像的感兴趣区域的灌注评价曲线TIC；

5)构造适合DFA分析的窗长阶数矩阵X和多项式拟合系数矩阵Fit；

6)对步骤4)获取的TIC曲线分别进行DFA分析提取TIC趋势。

6.如权利要求5所述的提取灌注时间强度曲线的方法，其特征在于，所述的窗长阶数矩阵X和多项式拟合系数矩阵Fit的构建为：

(1)适合DFA分析的窗长阶数矩阵X由原始TIC数据长度Num、DFA分析半窗长N及多项式拟合阶数K决定；

所述适合DFA分析的窗长阶数矩阵X如式(6)所示；

(2)适合DFA分析的多项式拟合系数矩阵Fit由原始TIC数据长度Num和DFA分析半窗长N决定；

所述适合DFA分析的多项式拟合系数矩阵Fit由最小二乘多元回归分析方法获取TIC数据片段的最优多项式拟合系数矩阵P如式7所示；

适合DFA分析的多项式拟合系数矩阵Fit如式8所示；

7.如权利要求5或6所述的提取灌注时间强度曲线的方法，其特征在于，所述的DFA分析提取TIC趋势包括以下操作：

(1)首先对原始TIC数据进行判断，若行数小于列数，则进行矩阵转秩，使行数大于列数；

(2)将TIC分为2N+1个相同长度的片段，并将片段长度记为窗长，相邻片段之间有N+1个点重叠；

(3)对于每一个片段寻找最优的阶数为K的多项式拟合，当K为0或1的时候，分别对应于常数和线性拟合；

(4)构造适合DFA分析的窗长阶数矩阵X和多项式拟合系数矩阵Fit；

(5)由最小二乘多元回归分析方法获取TIC数据片段的最优多项式拟合系数矩阵P；

(6)判断整段TIC数据能否被窗长整除，若能够整除则求多项式线性拟合的加权系数，若不能则首先对最后一个长度小于窗长的片段进行加权拟合；

(7)将第i个片段和第i+1个片段的经拟合的多项式分别标记为yⁱ(l₁)和yⁱ⁺¹(l₂)，l₁,l₂＝1,...,2N+1；

重叠片段的拟合多项式y^o(l)为：y^o(l)＝w₁yⁱ(l+N)+w₂yⁱ⁺¹(l),l＝1,...,N+1，其中w₁＝1-(l-1)/N,w₂＝(l-1)/N，为多项式线性拟合的加权系数；

(8)依次对重叠部分数据进行加权拟合；

(9)将最后一个长度小于窗长的片段及其之前的部分的加权拟合结果合并，输出拟合TIC曲线。

8.如权利要求7所述的提取灌注时间强度曲线的方法，其特征在于，对拟合TIC曲线还进行评价拟合质量：

(1)选择拟合质量QOF作为时间强度曲线整体趋势的评价指标；

(2)选择细节评价指数DEI作为时间强度曲线整体趋势的评价指标；

(3)利用matlab软件编程实现QOF、DEI计算；

(4)分别对TIC数据进行DFA拟合，计算QOF、DEI和运算时间。

9.如权利要求8所述的提取灌注时间强度曲线的方法，其特征在于，所述的拟合质量QOF达到最大值100％时表明TIC的拟合质量最佳，当QOF低于80％时认为TIC的拟合不准确，不具备可信性；

拟合质量QOF如式(9)、(10)、(11)所示：

$\mathrm{QOF}=(1-\frac{{\mathrm{SSR}}_f}{{\mathrm{VAR}}_d})\&times;100\%---\left(9\right)$

${\mathrm{SSR}}_f=\underset{t=t_{\mathrm{start}}}{\overset{t_{\mathrm{end}}}{\mathrm{\&Sigma;}}}{({\mathrm{TIC}}_f\left(t\right)-{\mathrm{TIC}}_d\left(t\right))}^2---\left(10\right)$

${\mathrm{VAR}}_d=\underset{t=t_{\mathrm{start}}}{\overset{t_{\mathrm{end}}}{\mathrm{\&Sigma;}}}{({\mathrm{TIC}}_f\left(t\right)-\frac{\underset{t=t_{\mathrm{start}}}{\overset{t_{\mathrm{end}}}{\mathrm{\&Sigma;}}}{\mathrm{TIC}}_d\left(t\right)}{N})}^2---\left(11\right)$

其中SSR_f为拟合区间内拟合TIC与滤波TIC的残差平方和，VAR_d为滤波TIC在拟合区间内的方差，TIC_f为拟合TIC，TIC_d为滤波TIC，t_start,t_end为拟合区间的起始点和终止点，N为拟合区间内的样本点数。

10.如权利要求8所述的提取灌注时间强度曲线的方法，其特征在于，所述的细节评价指数DEI越大表明拟合后的TIC对原始数据的细节保留能力越强，当DEI小于0.85时认为拟合TIC对灌注细节失真，不具备可信性；

DEI的定义为：

$\mathrm{DEI}=\frac{\mathrm{\&Sigma;}\left(\right|{\mathrm{TIC}}_f\left(t\right)-{\mathrm{TIC}}_f(t-1)|+|{\mathrm{TIC}}_f\left(t\right)-{\mathrm{TIC}}_f(t+1)\left|\right)}{\mathrm{\&Sigma;}\left(\right|{\mathrm{TIC}}_d\left(t\right)-{\mathrm{TIC}}_d(t-1)|+|{\mathrm{TIC}}_d\left(t\right)-{\mathrm{TIC}}_d(t+1)\left|\right)}---\left(12\right)$

式中：TIC_f(t)和TIC_d(t)分别为拟合前后的TIC值，DEI最大值为1，最小值为0。