[发明专利]血管与斑块的三维力学及组织特性成像检测方法

摘要

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

主权项

血管与斑块的三维力学特性及组织特性成像检测方法，其特征在于，包括以下步骤：步骤一、使用多波束聚焦波与超快速平面波交替发射序列，从颈动脉血管与斑块的径向与圆周方向进行弹性成像，求取脉搏波沿血管轴向的局部传播速度，获得血管径向、圆周及轴向的三维力学特性；计算斑块部位多点血管壁剪切率随心动周期的变化曲线；根据膨胀波形提取血管壁的扩张系数、顺应系数以及硬度指数表征血管硬化程度的参数；步骤二、使用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)测量内‑中膜以及中‑外膜两条轮廓线上各离散点之间的距离，求内中膜的平均厚度，最小厚度以及最大厚度参数。