CN104398271A - 血管与斑块的三维力学及组织特性成像检测方法 - Google Patents

Abstract

血管与斑块的三维力学特性及组织特性成像检测方法，基于多波束聚焦波与超快速平面波交替发射序列又可扩展至逐线扫描成像方式，分别从血管的径向、圆周及轴向三维力学特性、血管壁剪切率、斑块形态以及组织定征等方面对血管斑块形态与功能进行评价与成像，为颈动脉易损斑块的检测提供新的方法来作为现有方法的提升。

Description

血管与斑块的三维力学及组织特性成像检测方法

技术领域

本发明涉及超声成像技术领域，特别涉及一种基于多波束聚焦波与超快速平面波交替发射序列的颈动脉血管与斑块的三维力学特性及组织特性成像检测方法。

背景技术

动脉粥样硬化是一种全身性的系统疾病，定期对颈部动脉进行超声检查成为临床上动脉粥样硬化斑块普查的基本方法。现有临床血管斑块的诊断技术主要依靠B超图像、脉冲多普勒、彩色多普勒等成像方法来确定斑块的位置、大小及血管的狭窄程度。虽然狭窄程度经常作为动脉粥样硬化诊疗过程中一个重要的常规测量参数，但是众多研究表明大部分急性脑血管事件的发生与狭窄程度并不存在直接的联系，而是与斑块的实际组成成分以及所处的血流动力学环境有密切的关系。因此，有必要发展一种成像序列分别从血管径向、圆周以及横向的三维力学特性、血管壁剪切率、斑块的形态以及组织定征等方面等对血管斑块形态与功能进行评价与成像。这将大大提升现有的斑块检测方法并提高对易损斑块的检出率，降低急性心血管事件的发生。

血管在血压作用下沿径向与圆周方向的形变以及脉搏波沿血管轴向的传播速度可以分别用于表征血管在三个方向上的力学特性。同时，血管壁剪切率的变化也是影响血管斑块发展与破裂的重要血流动力学因素。然而，传统B超成像模式是由逐条扫描线依次聚焦发射接收实现的，这种成像模式的成像帧频只有几十赫兹，而人体脉搏波波速可达到数米/秒，因此传统超声成像的帧率难以捕捉脉搏波传播。除此之外，传统B超中脉冲多普勒发射波束与血流之间存在一定的夹角，只能获取血管腔内沿声束剖面不同深度处的流速分布图，降低了血管壁剪切率的估计精度。同时，常规脉冲多普勒技术测量流速时，存在频谱展宽效应，受声束几何特性、采样单元内流速梯度等因素的影响，造成频谱模糊，所测得的血流流速不准确。如果需要在一次发射中实现血管与斑块沿径向、圆周以及横向的三维力学特性评价、同时获取时空匹配的血管壁剪切率、以及斑块的形态与组织特征，就有必要发展一种复合成像序列，来克服传统超声成像模式的局限性，通过多角度的参数分析，实现对血管与斑块的形态与功能的同时评价。

发明内容

针对高帧率与高脉冲重复频率的要求，为了实现对颈动脉血管与斑块的形态和功能的成像检测及参数提取，本发明提供一种血管与斑块的三维力学特性及组织特性成像检测方法，基于多波束聚焦波与超快速平面波交替发射序列又可扩展至逐线扫描成像方式，分别从血管的径向、圆周及轴向三维力学特性、血管壁剪切率、斑块形态以及组织定征等方面对血管斑块形态与功能进行评价与成像，实现颈动脉血管与斑块的三维力学特性检测以及多点血管壁剪切率的同时测量，可显著提高血管壁周围位移估计的精度。

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

血管与斑块的三维力学特性及组织特性成像检测方法，包括以下步骤：

步骤一、使用多波束聚焦波与超快速平面波交替发射序列，从颈动脉血管与斑块的径向与圆周方向进行弹性成像，求取脉搏波沿血管轴向的局部传播速度，获得血管径向、圆周及轴向的三维力学特性；计算斑块部位多点血管壁剪切率随心动周期的变化曲线；根据膨胀波形提取血管壁的扩张系数、顺应系数以及硬度指数表征血管硬化程度的参数；

步骤二、使用Nakagami分布模型估计血管及斑块回波包络的概率密度分布曲线并使用由粗到精的M参量成像方法估计Nakagami模型中形状参数M的大小并成像，用于斑块中内部成分的组织定征；

步骤三、在B超图像的血管壁周围选取感兴趣区域，利用梯度及模糊C均值聚类对血管内-中膜以及中-外膜进行初始提取，之后经分段哈夫变换与GVF-Snake活动轮廓模型最终实现血管及斑块内中膜厚度的半自动测量；

步骤一所述的使用多波束聚焦波与超快速平面波交替发射序列，从颈动脉血管与斑块的径向与圆周方向进行弹性成像，测量脉搏波沿血管轴向的局部传播速度，获得血管径向、圆周及轴向的三维力学特性；计算斑块部位多点血管壁剪切率随心动周期的变化曲线；根据膨胀波形提取血管壁的扩张系数、顺应系数以及硬度指数表征血管硬化程度的参数，具体步骤为：

(1)设计多波束聚焦波与超快速平面波交替发射序列：根据线阵换能器的阵元的个数及尺寸，多波束聚焦波由四条间隔均匀的扫描线依次发射构成，每根扫描线由32个阵元聚焦发射，扫描线声束方向与血管内血流方向垂直；超快速平面波由超声探头128个阵元无聚焦无延时发射构成；多波束聚焦波与超快速平面波交替发射序列由四条聚焦波依次发射完成后发射一次平面波实现；利用延时叠加法交替获得的平面波与聚焦波的射频回波数据，分别用于血管和斑块三维力学特性的检测以及血管壁剪切率的测量；

(2)选用交替接收的平面波数据，用于计算血管和斑块沿径向、圆周及轴向的三维力学特性：首先，利用自相关方法计算平面波射频数据的自相关系数用于区分血液和周围组织；之后，对血液区域乘以校正系数，使两区域超声射频数据的幅度相近，实现血管腔内回波幅度变换；使用幅度变换后的射频数据，利用归一化互相关的方法计算二维位移场；采用阈值判别方法判断位移场中较大值奇异点，对于奇异点处，增大纵向与横向窗长再次计算，实现奇异值的去除；为了减小血管壁与血液交界区域不连续运动场的干扰，使用之前计算得到的自相关系数去除位移场中的血液信息，之后利用二维最小二乘应变估计器估计沿横向、轴向以及切向的应变，最后，通过矩阵变换来计算血管沿径向与圆周的弹性图像；

手动选择参考帧第一条扫描线前后血管壁的信号，之后每帧数据每条扫描线均自动选择该位置信号进行处理，利用互相关法计算每一根扫描线处相邻两帧之间血管前后壁的运动波形，积分并相减得到血管壁的膨胀波形即血管内径变化曲线，确定该膨胀波形中收缩期脚点所对应的时间，通过对每一根扫描线所在血管的位置与其膨胀波形收缩期脚点的时间进行线性回归分析求取脉搏波沿血管轴向的局部传播速度，最终获取血管及斑块沿血管轴向的力学特性；

(3)选用步骤(2)中计算的一个心动周期的相邻帧之间的二维位移场与径向弹性图像，手动选取感兴趣区域，计算该区域位移与应变的均值，并按时间绘制成曲线，即渐进位移与应变曲线；测量渐进位移与应变曲线的峰峰值用于表征血管与斑块在心脏快速收缩期的瞬时最大位移与应变；之后对渐进位移与应变曲线进行积分得到累加位移与应变曲线，测量其峰值，用于表征血管与斑块在一个心动周期的最大位移与形变；

(4)选用步骤(2)中计算的血管壁的膨胀波形，记录的臂部收缩压与舒张压按照式(7)计算颈部血管压力波形p(t)，之后，按照式(9)分别计算血管的顺应系数(CC)、扩张系数(DC)以及硬度指数(β)；

$\mathrm{\α}=\frac{A_d\ln \left(\frac{p_s}{p_d}\right)}{A_s-A_d},p\left(t\right)=p_d{e}^{\mathrm{\α}(\frac{A\left(t\right)}{A_d}-1)}---\left(7\right)$

其中，p_s是收缩期臂部血压，p_d是舒张期臂部血压，A(t)是根据膨胀波形计算的颈动脉血管横截面积变化波形，A_s是收缩期血管横截面积，A_d是舒张期血管横截面积，(式(7)中没有这2个变量，请确认)；将计算得到的p(t)按照下式(8)迭代得到优化的α值：

$\mathrm{\α}=\frac{{\stackrel{\‾}{p}}_{\mathrm{BrA}}}{{\stackrel{\‾}{p}}_{\mathrm{CCA}}}\mathrm{\α}---\left(8\right)$

其中， ${\stackrel{\‾}{p}}_{\mathrm{BrA}}=p_d+(p_s-p_d)/3,{\stackrel{\‾}{p}}_{\mathrm{CCA}}=\mathrm{average}\left(p\left(t\right)\right),$ 将优化后的α值带入(7)式重复上述过程，直到与之间的差异小于0.01；

分别使用颈动脉的膨胀波形与血压波形计算收缩期与舒张期的血管横截面变化ΔA、内径变化Δd与血压变化ΔP，按照式(9)分别计算血管的顺应系数CC、扩张系数DC与硬度指数β，其中d_d是舒张期血管内径；

顺应系数 $\mathrm{CC}=\frac{\mathrm{\ΔA}}{\mathrm{\ΔP}},$ 扩张系数 $\mathrm{DC}=\frac{\mathrm{\ΔA}}{A_d\mathrm{\ΔP}},$ 硬度指数 $\mathrm{\β}=\frac{\ln (p_s/p_d)}{\left[(\mathrm{\Δd}/d_d)\right]}---\left(9\right)$

(5)利用横向脉冲多普勒技术，选用交替接收到的多个聚焦波束信号，经正交解调、低通滤波、壁滤波、加窗及FFT变换后得到横向多普勒谱；由式(10)计算多普勒谱频宽所对应的流速v_max；根据不同深度所测的流速绘制血管剖面的分布曲线，血管壁剪切率由流速剖面对血管壁位置求导得到；之后，绘制血管壁剪切率随时间的变化曲线；

$\mathrm{Bd}=\frac{v_{\max }\·f_0}{c}(\frac{W_1}{\sqrt{{(W_1/2)}^2+F^2}}+\frac{W_2}{\sqrt{{(W_2/2)}^2+F^2}})---\left(10\right)$

其中，W₁与W₂分别为发射与接收孔径、F为焦点深度，Bd为横向多普勒谱宽，f₀为换能器实际工作频率，c为声速；

步骤二所述的使用Nakagami分布模型估计血管及斑块回波包络的概率密度分布曲线并使用由粗到精M参量成像方法估计Nakagami模型中形状参数M的大小并成像，用于斑块中内部成分的组织定征，具体步骤为：

(1)利用希尔伯特变换对射频数据取包络，将包络检波之后的数据通过降2采样进行金字塔分层，其中第零层为最底层，第三层为最高层；

(2)从金字塔最高层开始，使用式(11)中二阶导逼近的方法，计算每个窗长内回波包络数据的M值作为其中心点的值，通过滑动窗长计算该层的M参量图像，其中窗长尺寸为三个脉冲长度：

${\tilde{m}}_{\mathrm{TP}}=\frac{1+\sqrt{1+(4y/3)}}{4y}---\left(11\right)$

式中，y＝ln(μ₂/G)，μ₂为窗长内信号的二阶矩，N为窗长内包含的数据点数，x为每个数据点的回波包络值；之后，使用作为式(12)递归迭代的初始值，进一步计算该层的M参量图像：

${\tilde{m}}_i=\frac{{\tilde{m}}_{i-1}\{\ln \left({\tilde{m}}_{i-1}\right)-\mathrm{\ψ}\left({\tilde{m}}_{i-1}\right)\}}{y}---\left(12\right)$

(3)将第三层计算得到的M参量图像按照第二层图像的尺寸进行插值作为第二层递归迭代的初始值，使用公式(12)计算该层的M参量图像；之后，重复步骤三逐层向下直到最底层得到最终的M参量图像；

(4)在回波包络图像中选择感兴趣血管斑块区域，绘制该区域回波包络的概率分布曲线，使用Nakagami模型拟合该概率分布曲线并记录感兴趣区域的M参量；

步骤三所述在B超图像的血管壁周围选取感兴趣区域，通过梯度及模糊C均值聚类对血管内-中膜以及中-外膜进行初始提取，之后经分段Hough变换与GVF-Snake活动轮廓生长最终实现血管及斑块内中膜厚度的半自动测量；具体步骤为：

(1)图像预处理：在颈总动脉CCA的B超图像上手动选取包含血管腔、血管壁及部分周围组织的感兴趣区域ROI；利用非线性各项异性扩散滤波器去除感兴趣区域ROI的散斑噪声；

(2)计算感兴趣区域内图像梯度，根据内膜梯度阈值寻找内膜边界所在位置，作为初始内-中膜边界线；将初步搜索到的内-中膜边界划分为若干个连通域线段，若某一连通域的起始点距前一连通域终点的纵向距离较大，则去除该连通域线段，实现轮廓线的修正；将修正过的轮廓线分成若干段，对每段内轮廓线离散点应用Hough变换检测直线，将各段检测的直线插值连接成完整的一条轮廓线，作为内-中膜边界线；

(3)利用模糊C均值聚类方法将血管壁区域划分为三个灰度级的区域块，灰色区域与白色区域相交的边界即为初始中-外膜轮廓线；将初步搜索到的中-外膜轮廓线分为若干个连通域，若某一连通域的起始点距前一连通域终点的纵向距离较大，则去除该连通域线段；将修正过的轮廓线分成若干段，每段内轮廓线离散点应用Hough变换检测直线，将各段检测的直线插值连接成完整的一条轮廓线，作为中-外膜边界线；

(4)利用参数活动轮廓GVF-snake模型在保持曲线平滑性和连续性的基础上，根据图像梯度力进一步演变上述步骤(2)和步骤(3)分别得到的内-中膜以及中-外膜轮廓曲线，使之收敛于真实的边界位置；每次迭代时轮廓线上单个离散点下一时刻位置是由该点以及其前后四个点的位置决定；每次迭代后，轮廓线的前后端点位置由轮廓曲线的中间点位置线性外插进行修正，再将修正后的轮廓曲线进入下次迭代，迭代一定次数后停止，从而获得最终的内-中膜以及中-外膜轮廓曲线；

(5)测量内-中膜以及中-外膜两条轮廓线上各离散点之间的距离，求内中膜的平均厚度，最小厚度以及最大厚度参数。

所述校正系数的求取方法是：分别在原始超声射频数据中的血液和周围组织部分取大小相同的矩形区域，求取相应的均值比，即为校正系数。

本发明的优点：

1.多波束聚焦波与超快速平面波交替发射序列结合了多点聚焦波的高脉冲重复频率与平面波高帧率的优点，实现颈动脉血管与斑块的三维力学特性检测以及多点血管壁剪切率的同时测量。

2.血管壁径向及圆周弹性成像通过血管腔内回波幅度变换来降低血液信号突变引起的互相关估计误差，利用自相关矩阵去除血液信息对血管壁位移的干扰，结合阈值判断的反向由粗到精分层算法去除大奇异点，可显著提高血管壁周围位移估计的精度。

3.脉搏波沿血管轴向传播波速测量中应用交替接收的平面波数据进行脉搏波波速测量，平面波发射时，128个阵元无聚焦无延时发射，扫描线之间不存在时间延迟，显著提高了成像的帧率，可实现局部超短距脉搏波传播速度的测量。

4.多点血管壁剪切率估计应用超声横向脉冲多普勒技术，声束与流速夹角为90°，生成的频谱宽度不受采样单元内流速梯度影响，只由采样单元内的最大流速决定，流速估计精度提高，可获得较为准确的沿血管横断面不同深度的流速分布曲线。此外，多点聚焦波的发射模式有利于同时检测同一段血管多个位置处的血管壁剪切率大小，获取同一段血管不同位置处的血液动力学特性。

5.多分辨由粗到精的Nakagami-M参量成像方法，通过将上一级大窗长的计算结果传递到下一级作为初始值来保证计算的稳定性，同时通过逐级减小计算窗长来保证图像较好的分辨性能。

6.血管壁内中膜厚度测量方法根据图像梯度及模糊C均值聚类方法提取感兴趣区域内血管壁内-中膜和中-外膜初始轮廓线，经连通域划分，根据前后连通域之间的距离去除距前一连通域距离较大的连通域线段，得到修正后的轮廓线，再利用Hough变换分段检测直线，将检测到的各段直线进行插值后通过GVF-snake模型进一步优化，实现血管内中膜厚度及血管斑块厚度的半自动测量。

附图说明

图1是多波束聚焦波与超快速平面波交替发射序列示意图。

图2是颈动脉血管与斑块的三维力学特性与组织特性的成像检测示意图。

图3是血管沿径向与圆周的弹性成像流程图。

图4是脉搏波沿血管轴向的局部传播速度测量流程图。

图5是基于横向脉冲多普勒技术的血管壁多点剪切率估计流程图。

图6是多分辨由粗到精的Nakagami-M参量成像流程图。

图7是血管及斑块内中膜厚度的半自动测量流程图。

具体实施方式

下面结合附图对本发明做详细叙述。

本发明血管与斑块的三维力学特性及组织特性成像检测方法，其具体实施方式，包括以下几个步骤：

步骤一、使用多波束聚焦波与超快速平面波交替发射序列，从颈动脉血管与斑块的径向与圆周方向进行弹性成像，求取脉搏波沿血管轴向的局部传播速度，获得血管径向、圆周及轴向的三维力学特性。计算斑块部位多点血管壁剪切率随心动周期的变化曲线，根据膨胀波形提取血管壁的扩张系数、顺应系数以及硬度指数等表征血管硬化程度的参数，具体步骤为：

(1)设计多波束聚焦波与超快速平面波交替发射序列如图1所示。根据线阵换能器的阵元的个数及尺寸，多波束聚焦波由四条间隔均匀的扫描线依次发射构成，每根扫描线由32个阵元聚焦发射，扫描线声束方向与血管内血流方向垂直；超快速平面波由超声探头128个阵元无聚焦无延时发射构成。多波束聚焦波与超快速平面波交替发射序列由四条聚焦波依次发射完成后发射一次平面波实现。复合发射序列中每条聚焦波的PRF以及平面波的成像帧频FR都为：

$\mathrm{PRF}=\mathrm{FR}=\frac{c}{2d}*\frac{1}{5}---\left(1\right)$

式中，c为声速，d为成像深度，其余成像参数见表1：

表1 复合发射序列的成像参数

(2)选用交替接收的平面波数据，用于颈动脉血管与斑块的三维力学特性与组织特性的成像检测，示意图如图2所示。接收交替发射序列的回波数据时，对平面波和聚焦波的回波数据进行交替接收。其中平面波发射时128个阵元的单通道数据经延时叠加法得到平面波的射频回波数据，可用于血管三维力学特性的分析以及M参量成像。同时，多波束聚焦波射频回波数据可用于多点血管壁剪切率的估计。

血管沿径向与圆周的弹性成像流程如图3所示。首先，在位移估计前对原始射频数据进行预处理，包括射频数据滤波和幅度变换。然后，进行位移场估计，包括基于归一化互相关的二维位移获取、奇异值去除、血液信息去除和二维应变场计算。最后，使用矩阵变换来计算血管沿径向与圆周的弹性图像。

数据预处理：读取相邻两帧的射频回波数据对其进行滤波处理，滤波器设置为以超声发射频率为中心频率，带宽为70％的巴特沃斯带通滤波器。对滤波后的射频数据进行数据幅度变换，调整血液与周围组织的超声RF数据幅度后再进行位移估计，降低数据幅度突变所产生的位移场误差。幅度变换的规则是：使用射频数据的自相关矩阵区分血液和周围组织；对血液区域乘以一定校正系数，使两区域超声RF数据的幅度相近。其校正系数的求取方法是：分别在原始超声射频数据中的血液和周围组织部分取大小相同的矩形区域，求取相应的均值比，即为校正系数。

基于归一化互相关的二维位移估计：对于幅度变换后的相邻两帧射频数据设置一定的数据窗长和重叠率，对对应数据窗内的信号应用二维互相关算法计算它们之间发生的位移，其归一化二维互相关系数计算公式如式(2)。二维互相关矩阵获得后，利用二维亚采样时延技术提高位移估计精度，通过对二维互相关矩阵两个方向上分别进行抛物线插值来实现纵向、横向位移场的亚采样间距位移估计。

$R(u,v)=\frac{{\mathrm{\Σ}}_{x,y}[\mathrm{Pre}(x,y)-{\stackrel{\‾}{\mathrm{Pre}}}_{u,v}][\mathrm{Post}(x-u,y-v)-{\stackrel{\‾}{\mathrm{Post}}}_{u,v}]}{\left\{{\mathrm{\Σ}}_{x,y}{[\mathrm{Pre}(x,y)-{\stackrel{\‾}{\mathrm{Pre}}}_{u,v}]}^2{\mathrm{\Σ}}_{x,y}{[\mathrm{Post}(x-u,y-v)-{\stackrel{\‾}{\mathrm{Post}}}_{u,v}]}^2\right\}}---\left(2\right)$

其中Pre表示形变前RF数据窗，Post表示形变后RF数据窗，为Pre的均值，为Post的均值。

由粗到精分层法去除奇异值：位移场中大部分奇异点的幅值接近于窗长点数的幅值，与实际位移相差较大，故主要是针对大值奇异点进行去除，可以先采用阈值判别方法判断较大值奇异点，使用反向由粗到精分层算法结合校正算法来实现大奇异值的去除，具体方法为对于奇异点处，增大纵向与横向窗长，并使用周围非奇异值点的位移去校正二维互相关位移估计，使用估算到的“粗”位移来替代奇异值点。最后，使用局部窗的差分中值滤波实现位移场剩余奇异值和噪声的去除。

基于自相关矩阵的血液信息去除：由于血液的流动性，前后两帧超声射频数据的相关性较低，血液处位移场存在大量的奇异值与噪声点，这会干扰血管壁位移场显示与应变估计，故需要使用射频数据的自相关矩阵区分血液和周围组织，并对位移场中血管内血液信息进行去除。

应变成像：以计算得到的二维位移场(u,v)作为输入，利用二维的最小二乘应变估计器(LSQSE)可以估计沿横向(E_xx)、轴向(E_yy)以及切向(E_xy)的应变。对于中心点为(x₀,y₀)，大小为(2M+1)×(2M+1)矩形窗内的二维位移场可以用以下拟合公式来表达：

$\left\{\begin{array}{c}u(x,y)=a_0+a_{1}x+a_{2}y\\ v(x,y)=b_0+b_{1}x+b_{2}y\end{array}\right.,\mathrm{for}\begin{array}{c}x=x_0-M..x_0+M\\ y=y_0-M..y_0+M\end{array}---\left(3\right)$

式(3)中的a_i,b_i(i＝0,1,2)分别为待求的多项式系数。该式最小二乘解矩阵表达为：

$\left[\begin{array}{c}{\hat{a}}_0\\ {\hat{a}}_1\\ {\hat{a}}_2\end{array}\right]=X^T\·{\left[{\mathrm{XX}}^T\right]}^{-1}\·u,\left[\begin{array}{c}{\hat{b}}_0\\ {\hat{b}}_1\\ {\hat{b}}_2\end{array}\right]=X^T\·{\left[{\mathrm{XX}}^T\right]}^{-1}\·v---\left(4\right)$

式(4)中是a_i,b_i(i＝0,1,2)的最小二乘估计结果。X是一个大小为(2M+1)²×3的矩阵，其第一列的值为1，第二列的值为横坐标，第三列值为纵坐标。通过矩阵运算可以求得多项式拟合的系数。小变形条件下的二维应变分量可以用这些系数表示为^[68]：

${\mathrm{\ϵ}}_x=\frac{\∂u}{\∂x}=a_1,{\mathrm{\ϵ}}_y=\frac{\∂v}{\∂y}=b_2,2{\mathrm{\ϵ}}_{\mathrm{xy}}=\frac{\∂u}{\∂y}+\frac{\text{\∂v}}{\∂y}=a_2+b_1---\left(5\right)$

直角坐标系下的血管径向及圆周应变可通过二维应变分量乘以以下矩阵变换得到^[69]：

$\left[\begin{array}{cc}{E}_{\mathrm{radial}}& E_{\mathrm{shear}}\\ E_{\mathrm{shear}}& E_{\mathrm{circum}}\end{array}\right]={\left[\begin{array}{cc}\cos \mathrm{\θ}& -\sin \mathrm{\θ}\\ \sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right]}^T\·\left[\begin{array}{cc}{E}_{\mathrm{xx}}& E_{\mathrm{xy}}\\ E_{\mathrm{xy}}& E_{\mathrm{yy}}\end{array}\right]\·\left[\begin{array}{cc}\cos \mathrm{\θ}& -\sin \mathrm{\θ}\\ \sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right]---\left(6\right)$

式中，θ就是直角坐标系与极坐标系之间的转换角度，通过手动选择血管腔内的中心点位置可以计算直角坐标系中任意点做对应的θ。式(6)中矩阵对角线上的E_radial与E_circum分别代表了血管壁沿径向及圆周方向的形变大小，也即径向与圆周弹性图像。

脉搏波沿血管轴向的局部传播速度的测量流程图如图4所示。手动选择参考帧第一条扫描线前后血管壁的信号(图中矩形框)，之后每帧数据每条扫描线均自动选择该位置信号进行处理。利用一维互相关技术对相邻两帧之间血管壁位置处的信号进行运动估计，可获得血管前后壁的运动波形。使用截止频率为80Hz的2阶巴特沃斯低通滤波器对血管壁运动波形进行滤波，目的是为了去除系统高频噪声所造成的运动波形中的“毛刺”现象。之后，通过血管壁运动曲线减去其平均值的方法去除呼吸所带来的基线漂移。对血管前后壁的运动波形积分并相减，可得到血管壁的膨胀波波形。最后，由各扫描线位置处膨胀波形上的收缩期脚点可求出各扫描线处膨胀波形之间的时延T。另外，由于局部血管的长度很难测量，因此我们通过扫描线之间的间距(阵元间距)来确定膨胀波形之间的距离L，最后对膨胀波形间的距离与时延进行线性回归分析求局部超短距脉搏波传播速度PWV＝dL/dT。

(3)选用步骤(2)中计算的一个心动周期的相邻帧之间的二维位移场与径向应变图，手动选取感兴趣区域，计算该区域位移与应变的均值，并按时间绘制成曲线，即渐进位移与应变曲线。测量渐进位移与应变曲线的峰峰值用于表征血管与斑块在心脏快速收缩期的瞬时最大位移与应变。之后对渐进位移与应变曲线进行积分得到累加位移与应变曲线，测量其峰值，用于表征血管与斑块在一个心动周期的最大位移与形变。

(4)使用血管壁的膨胀波形与臂部测量压力按照式(7)计算颈部血管压力波形p(t)，之后，按照式(9)分别计算血管的顺应系数(CC)、扩张系数(DC)与硬度指数(β)。

$\mathrm{\α}=\frac{A_d\ln \left(\frac{p_s}{p_d}\right)}{A_s-A_d},p\left(t\right)=p_d{e}^{\mathrm{\α}(\frac{A\left(t\right)}{A_d}-1)}---\left(7\right)$

其中，p_s是收缩期臂部血压，p_d是舒张期臂部血压，A(t)是根据膨胀波形计算的颈动脉血管横截面积变化波形，A_s是收缩期血管横截面积，A_d是舒张期血管横截面积，d_s是收缩期血管内径，d_d是舒张期血管内径。将计算得到的p(t)按照下式迭代得到优化的α值：

$\mathrm{\α}=\frac{{\stackrel{\‾}{p}}_{\mathrm{BrA}}}{{\stackrel{\‾}{p}}_{\mathrm{CCA}}}\mathrm{\α}---\left(8\right)$

其中， ${\stackrel{\‾}{p}}_{\mathrm{BrA}}=p_d+(p_s-p_d)/3,{\stackrel{\‾}{p}}_{\mathrm{CCA}}=\mathrm{average}\left(p\left(t\right)\right),$ 将优化后的α值带入(7)式重复上述过程，直到与之间的差异小于0.01。

分别使用颈动脉的膨胀波形与血压波形计算收缩期与舒张期的血管横截面变化ΔA、内径变化Δd与血压变化ΔP，按照式(9)分别计算血管的顺应系数(CC)、扩张系数(DC)与硬度指数(β)，其中A_d为舒张期血管横截面积

(5)利用横向脉冲多普勒技术，选用交替接收到的多个聚焦波束信号，根据血管腔内感兴趣区域选择距离选通的位置及窗长，距离选通后的信号经正交解调得到其I、Q分量，并经过低通滤波滤除正交解调过程中产生的二倍频分量，得到包含频移信息的I、Q分量。为了减小血管壁及组织运动对频移的影响，需要使用壁滤波器去除这部分低频成分。将I、Q分量进行复信号合成，通过采样保持获得采样单元处位移曲线，并做FFT变换获得该距离选通的血流信号的横向多普勒谱，信号处理过程见图5。根据式(10)计算多普勒谱频宽所对应的流速。之后，根据沿血管腔横截面不同深度所测的流速分布计算血管剪切率的分布，选取血管壁处的剪切率沿时间轴展开即可获取血管壁剪切率随时间的变化曲线。

$\mathrm{Bd}=\frac{v_{\max }\·f_0}{c}(\frac{W_1}{\sqrt{{(W_1/2)}^2+F^2}}+\frac{W_2}{\sqrt{{(W_2/2)}^2+F^2}})---\left(10\right)$

其中，W₁与W₂分别为发射与接收孔径，F为焦点深度，Bd为横向多普勒谱宽，f₀为换能器实际工作频率，c为声速。

步骤二、使用Nakagami分布模型估计血管及斑块回波包络的概率密度分布曲线并使用由粗到精的M参量成像方法估计Nakagami模型中形状参数M的大小并成像，用于斑块中内部成分的组织定征，多分辨由粗到精的Nakagami-M参量成像流程如图6所示，具体实现步骤为：

(1)利用希尔伯特变换对射频数据取包络，将包络检波之后的数据通过降2采样进行金字塔分层，其中第零层为最底层，第三层为最高层。

(2)从金字塔最高层开始，使用式(11)中二阶导逼近的方法，计算每个窗长内回波包络数据的M值作为其中心点的值，通过滑动窗长计算该层的M参量图像，其中窗长尺寸为三个脉冲长度：

${\tilde{m}}_{\mathrm{TP}}=\frac{1+\sqrt{1+(4y/3)}}{4y}---\left(11\right)$

式中，y＝ln(μ₂/G)，μ₂为窗长内信号的二阶矩，N为窗长内包含的数据点数，x为每个数据点的回波包络值。之后，使用作为式(12)递归迭代的初始值，进一步计算该层的M参量图像，具体计算公式为：

${\tilde{m}}_i=\frac{{\tilde{m}}_{i-1}\{\ln \left({\tilde{m}}_{i-1}\right)-\mathrm{\ψ}\left({\tilde{m}}_{i-1}\right)\}}{y}---\left(12\right)$

(3)将第三层计算得到的M参量图像按照第二层图像的尺寸进行插值作为第二层递归迭代的初始值，使用公式(12)计算该层的M参量图像。之后，重复步骤三逐层向下直到最底层得到最终的M参量图像。

(4)在回波包络图像中选择感兴趣血管斑块区域，绘制该区域回波包络的概率分布曲线，使用Nakagami模型拟合该概率分布曲线并记录感兴趣区域的M参量。

步骤三、在B超图像的血管壁周围选取感兴趣区域，通过梯度及模糊C均值聚类对血管内-中膜以及中-外膜进行初始提取，之后经分段Hough变换与GVF-Snake活动轮廓生长最终实现血管及斑块内中膜厚度的半自动测量。信号处理流程如图7所示，具体步骤为：

(1)图像预处理：在颈总动脉(CCA)B超图像上手动选取包含血管腔、血管壁及部分周围组织的感兴趣区域(ROI)。由于超声图像存在散斑噪声，影响后续的图像处理，故对ROI区域首先进行散斑去噪处理，此处采用非线性各项异性扩散SRAD滤波器，在滤除噪声的同时能有效保持图像边缘信息。

(2)内-中膜初始轮廓提取

使用图像预处理后的ROI图像计算梯度矩阵，对梯度矩阵从上至下逐列搜索，当搜索到某一位置处的梯度值大于该列的梯度阈值(本研究中各列的梯度阈值为每列梯度最大值的20％)时，则初步认为该位置为该列的内-中膜边界位置，依次搜索直至所有列都搜索完毕，即可得到初步的内-中膜轮廓点。

将初步搜索到的内-中膜轮廓点划分为若干个连通域线段，若某一连通域的起始点相距前一连通域终点之间的距离较大，则去除该连通域线段，直至判断完所有连通域线段以去除与真实内-中膜边界偏离较大的奇异点。将上步修正过的轮廓点分成若干段，每段内轮廓点利用Hough变换检测直线，将各段检测的直线插值连接成一条完整的轮廓线，即内-中膜初始轮廓。

(3)中-外膜初始轮廓提取：利用模糊C均值聚类分割方法将ROI区域划分为三个灰度级的区域块，灰色区域与白色区域相交的边界即为中-外膜轮廓线。确定外膜的初始位置后，通过去除奇异点与分段Hough变换修正中-外膜轮廓线，此部分同步骤(2)。

(4)基于GVF-snake模型的轮廓优化

本研究采用的参数活动轮廓GVF-Snake模型是能量方程最小化模型，能量方程如式(13)。基于GVF-snake模型的轮廓曲线优化过程首先计算得到GVF梯度矢量流场，建立GVF-snake模型，最后将内-中膜与中-外膜两条初始轮廓线代入GVF-snake模型反复迭代至收敛。

$E\left(r\left(s\right)\right)=\underset{0}{\overset{1}{\∫}}[E_{\mathrm{int}}\left(r\left(s\right)\right)+E_{\mathrm{image}}\left(r\left(s\right)\right)]\mathrm{ds}=\underset{0}{\overset{1}{\∫}}[\mathrm{\α}\·{\left|r_s\left(s\right)\right|}^2/2+\mathrm{\β}\·{{|r}_{\mathrm{ss}}\left(s\right)|}^2/2+k\·E_{\mathrm{image}}\left(r\left(s\right)\right)]\mathrm{ds}---\left(13\right)$

GVF梯度矢量流场的计算需要先对图像预处理后的ROI区域根据梯度算子得到其边缘图。在计算边缘图时，只计算图像矩阵的垂直梯度，剔除掉较小的梯度值和负梯度值，保留较大的正梯度值，减小其它非内中膜边缘对梯度矢量流场GVF产生的影响。利用式(14)可将梯度矢量延伸到离边界较远的区域和同质区域。GVF矢量场V(x,y)＝[u(x,y)v(x,y)]定义为由最小化以下能量泛函得到：

$\begin{array}{c}{E}_{\mathrm{GVF}}=\∫\∫\mathrm{\μ}({|\▿u|}^2+{|\▿v|}^2)+{|\▿f|}^2{|V-\▿f|}^2\mathrm{dxdy}\\ =\∫\∫\mathrm{\μ}({u_x}^2+{u_y}^2+{v_x}^2+{v_y}^2)+{|\▿f|}^2[{(u-f_x)}^2+{(v-f_y)}^2]\mathrm{dxdy}\end{array}---\left(14\right)$

将GVF梯度矢量流场代入GVF-snake模型后，就可对轮廓曲线进行反复迭代，在每次迭代后，需要处理轮廓曲线两端的边界点，根据轮廓曲线的中间点线性拟合插值形成新的边界点，再进入下一次迭代，反复迭代至停止。

(5)IMT计算：测量内-中膜以及中-外膜两条轮廓线上各离散点之间的距离，求内中膜的平均厚度，最小厚度以及最大厚度等参数。

Claims (2)

1.血管与斑块的三维力学特性及组织特性成像检测方法，其特征在于，包括以下步骤：

步骤一、使用多波束聚焦波与超快速平面波交替发射序列，从颈动脉血管与斑块的径向与圆周方向进行弹性成像，求取脉搏波沿血管轴向的局部传播速度，获得血管径向、圆周及轴向的三维力学特性；计算斑块部位多点血管壁剪切率随心动周期的变化曲线；根据膨胀波形提取血管壁的扩张系数、顺应系数以及硬度指数表征血管硬化程度的参数；

步骤二、使用Nakagami分布模型估计血管及斑块回波包络的概率密度分布曲线并使用由粗到精的M参量成像方法估计Nakagami模型中形状参数M的大小并成像，用于斑块中内部成分的组织定征；

步骤三、在B超图像的血管壁周围选取感兴趣区域，利用梯度及模糊C均值聚类对血管内-中膜以及中-外膜进行初始提取，之后经分段哈夫变换与GVF-Snake活动轮廓模型最终实现血管及斑块内中膜厚度的半自动测量；

步骤一所述的使用多波束聚焦波与超快速平面波交替发射序列，从颈动脉血管与斑块的径向与圆周方向进行弹性成像，测量脉搏波沿血管轴向的局部传播速度，获得血管径向、圆周及轴向的三维力学特性；计算斑块部位多点血管壁剪切率随心动周期的变化曲线；根据膨胀波形提取血管壁的扩张系数、顺应系数以及硬度指数表征血管硬化程度的参数，具体步骤为：

(1)设计多波束聚焦波与超快速平面波交替发射序列：根据线阵换能器的阵元的个数及尺寸，多波束聚焦波由四条间隔均匀的扫描线依次发射构成，每根扫描线由32个阵元聚焦发射，扫描线声束方向与血管内血流方向垂直；超快速平面波由超声探头128个阵元无聚焦无延时发射构成；多波束聚焦波与超快速平面波交替发射序列由四条聚焦波依次发射完成后发射一次平面波实现；利用延时叠加法交替获得的平面波与聚焦波的射频回波数据，分别用于血管和斑块三维力学特性的检测以及血管壁剪切率的测量；

(2)选用交替接收的平面波数据，用于计算血管和斑块沿径向、圆周及轴向的三维力学特性：首先，利用自相关方法计算平面波射频数据的自相关系数用于区分血液和周围组织；之后，对血液区域乘以校正系数，使两区域超声射频数据的幅度相近，实现血管腔内回波幅度变换；使用幅度变换后的射频数据，利用归一化互相关的方法计算二维位移场；采用阈值判别方法判断位移场中较大值奇异点，对于奇异点处，增大纵向与横向窗长再次计算，实现奇异值的去除；为了减小血管壁与血液交界区域不连续运动场的干扰，使用之前计算得到的自相关系数去除位移场中的血液信息，之后利用二维最小二乘应变估计器估计沿横向、轴向以及切向的应变，最后，通过矩阵变换来计算血管沿径向与圆周的弹性图像；

手动选择参考帧第一条扫描线前后血管壁的信号，之后每帧数据每条扫描线均自动选择该位置信号进行处理，利用互相关法计算每一根扫描线处相邻两帧之间血管前后壁的运动波形，积分并相减得到血管壁的膨胀波形即血管内径变化曲线，确定该膨胀波形中收缩期脚点所对应的时间，通过对每一根扫描线所在血管的位置与其膨胀波形收缩期脚点的时间进行线性回归分析求取脉搏波沿血管轴向的局部传播速度，最终获取血管及斑块沿血管轴向的力学特性；

(3)选用步骤(2)中计算的一个心动周期的相邻帧之间的二维位移场与径向弹性图像，手动选取感兴趣区域，计算该区域位移与应变的均值，并按时间绘制成曲线，即渐进位移与应变曲线；测量渐进位移与应变曲线的峰峰值用于表征血管与斑块在心脏快速收缩期的瞬时最大位移与应变；之后对渐进位移与应变曲线进行积分得到累加位移与应变曲线，测量其峰值，用于表征血管与斑块在一个心动周期的最大位移与形变；

(4)选用步骤(2)中计算的血管壁的膨胀波形，记录的臂部收缩压与舒张压按照式(7)计算颈部血管压力波形p(t)，之后，按照式(9)分别计算血管的顺应系数(CC)、扩张系数(DC)以及硬度指数(β)；

$\mathrm{\α}=\frac{A_d\ln \left(\frac{p_s}{p_d}\right)}{A_s-A_d},p\left(t\right)=p_d{e}^{\mathrm{\α}(\frac{A\left(t\right)}{A_d}-1)}---\left(7\right)$

其中，p_s是收缩期臂部血压，p_d是舒张期臂部血压，A(t)是根据膨胀波形计算的颈动脉血管横截面积变化波形，A_s是收缩期血管横截面积，A_d是舒张期血管横截面积；将计算得到的p(t)按照下式(8)迭代得到优化的α值：

$\mathrm{\α}=\frac{{\stackrel{\‾}{p}}_{\mathrm{BrA}}}{{\stackrel{\‾}{p}}_{\mathrm{CCA}}}\mathrm{\α}---\left(8\right)$

其中， ${\stackrel{\‾}{p}}_{\mathrm{BrA}}=p_d+(p_s-p_d)/3,{\stackrel{\‾}{p}}_{\mathrm{CCA}}=\mathrm{average}\left(p\left(t\right)\right),$ 将优化后的α值带入(7)式重复上述过程，直到与之间的差异小于0.01；

分别使用颈动脉的膨胀波形与血压波形计算收缩期与舒张期的血管横截面变化ΔA、内径变化Δd与血压变化ΔP，按照式(9)分别计算血管的顺应系数CC、扩张系数DC与硬度指数β，其中d_d是舒张期血管内径；

顺应系数 $\mathrm{CC}=\frac{\mathrm{\ΔA}}{\mathrm{\ΔP}},$ 扩张系数 $\mathrm{DC}=\frac{\mathrm{\ΔA}}{A_d\mathrm{\ΔP}},$ 硬度指数 $\mathrm{\β}=\frac{\ln (p_s/p_d)}{\left[(\mathrm{\Δd}/d_d)\right]}---\left(9\right)$

(5)利用横向脉冲多普勒技术，选用交替接收到的多个聚焦波束信号，经正交解调、低通滤波、壁滤波、加窗及FFT变换后得到横向多普勒谱；由式(10)计算多普勒谱频宽所对应的流速v_max；根据不同深度所测的流速绘制血管剖面的分布曲线，血管壁剪切率由流速剖面对血管壁位置求导得到；之后，绘制血管壁剪切率随时间的变化曲线；

$\mathrm{Bd}=\frac{v_{\max }\·f_0}{c}(\frac{W_1}{\sqrt{{(W_1/2)}^2+F^2}}+\frac{W_2}{\sqrt{{(W_2/2)}^2+F^2}})---\left(10\right)$

其中，W₁与W₂分别为发射与接收孔径、F为焦点深度，Bd为横向多普勒谱宽，f₀为换能器实际工作频率，c为声速；

步骤二所述的使用Nakagami分布模型估计血管及斑块回波包络的概率密度分布曲线并使用由粗到精M参量成像方法估计Nakagami模型中形状参数M的大小并成像，用于斑块中内部成分的组织定征，具体步骤为：

(1)利用希尔伯特变换对射频数据取包络，将包络检波之后的数据通过降2采样进行金字塔分层，其中第零层为最底层，第三层为最高层；

(2)从金字塔最高层开始，使用式(11)中二阶导逼近的方法，计算每个窗长内回波包络数据的M值作为其中心点的值，通过滑动窗长计算该层的M参量图像，其中窗长尺寸为三个脉冲长度：

${\tilde{m}}_{\mathrm{TP}}=\frac{1+\sqrt{1+(4y/3)}}{4y}---\left(11\right)$

式中，y＝ln(μ₂/G)，μ₂为窗长内信号的二阶矩，N为窗长内包含的数据点数，x为每个数据点的回波包络值；之后，使用作为式(12)递归迭代的初始值，进一步计算该层的M参量图像：

${\tilde{m}}_i=\frac{{\tilde{m}}_{i-1}\left\{\ln \left({\tilde{m}}_{i-1}\right)\mathrm{\ψ}\left({\tilde{m}}_{i-1}\right)\right\}}{y}---\left(12\right)$

(3)将第三层计算得到的M参量图像按照第二层图像的尺寸进行插值作为第二层递归迭代的初始值，使用公式(12)计算该层的M参量图像；之后，重复步骤三逐层向下直到最底层得到最终的M参量图像；

(4)在回波包络图像中选择感兴趣血管斑块区域，绘制该区域回波包络的概率分布曲线，使用Nakagami模型拟合该概率分布曲线并记录感兴趣区域的M参量；

步骤三所述在B超图像的血管壁周围选取感兴趣区域，通过梯度及模糊C均值聚类对血管内-中膜以及中-外膜进行初始提取，之后经分段Hough变换与GVF-Snake活动轮廓生长最终实现血管及斑块内中膜厚度的半自动测量；具体步骤为：

(1)图像预处理：在颈总动脉CCA的B超图像上手动选取包含血管腔、血管壁及部分周围组织的感兴趣区域ROI；利用非线性各项异性扩散滤波器去除感兴趣区域ROI的散斑噪声；

(2)计算感兴趣区域内图像梯度，根据内膜梯度阈值寻找内膜边界所在位置，作为初始内-中膜边界线；将初步搜索到的内-中膜边界划分为若干个连通域线段，若某一连通域的起始点距前一连通域终点的纵向距离较大，则去除该连通域线段，实现轮廓线的修正；将修正过的轮廓线分成若干段，对每段内轮廓线离散点应用Hough变换检测直线，将各段检测的直线插值连接成完整的一条轮廓线，作为内-中膜边界线；

(3)利用模糊C均值聚类方法将血管壁区域划分为三个灰度级的区域块，灰色区域与白色区域相交的边界即为初始中-外膜轮廓线；将初步搜索到的中-外膜轮廓线分为若干个连通域，若某一连通域的起始点距前一连通域终点的纵向距离较大，则去除该连通域线段；将修正过的轮廓线分成若干段，每段内轮廓线离散点应用Hough变换检测直线，将各段检测的直线插值连接成完整的一条轮廓线，作为中-外膜边界线；

(4)利用参数活动轮廓GVF-snake模型在保持曲线平滑性和连续性的基础上，根据图像梯度力进一步演变上述步骤(2)和步骤(3)分别得到的内-中膜以及中-外膜轮廓曲线，使之收敛于真实的边界位置；每次迭代时轮廓线上单个离散点下一时刻位置是由该点以及其前后四个点的位置决定；每次迭代后，轮廓线的前后端点位置由轮廓曲线的中间点位置线性外插进行修正，再将修正后的轮廓曲线进入下次迭代，迭代一定次数后停止，从而获得最终的内-中膜以及中-外膜轮廓曲线；

(5)测量内-中膜以及中-外膜两条轮廓线上各离散点之间的距离，求内中膜的平均厚度，最小厚度以及最大厚度参数。

2.根据权利要求1所述的方法，其特征在于，所述校正系数的求取方法是：分别在原始超声射频数据中的血液和周围组织部分取大小相同的矩形区域，求取相应的均值比，即为校正系数。