CN103340620B - 一种管壁应力相位角的测量方法和系统 - Google Patents

Abstract

本申请公开了一种管壁应力相位角的测量方法，包括：在流体中加入超声造影微泡；使用高频超声成像装置采集连续采集多帧感兴趣区域的超声图像；对相邻两幅超声图像对应区域进行互相关分析，得到管壁的周向应力和流场剪切应力；绘制管壁周向应力和流场剪切应力相对于时间变化的波形图，测量应力相位角。本申请还公开了一种管壁应力相位角的测量系统。在本申请的具体实施方式中，通过高频超声成像装置和超声微泡成像，对相邻帧图像进行互相关分析，能同步、实时、精确地获得周向应力和流场剪切应力相对于时间的波形图，从而测量应力相位角。

Description

一种管壁应力相位角的测量方法和系统

技术领域

本申请涉及管壁应力的测量技术，尤其涉及一种管壁应力相位角的测量方法和系统。

背景技术

心血管疾病（CVD）已成为人类死亡病因最高的“头号杀手”，也是人们健康的“无声凶煞”，亟需发展新型早期诊断方法。血管内皮细胞衬于整个心血管系统的内表面，是血管壁与血液之间的分界细胞，是形成心血管封闭管道系统的形态基础。血管内皮细胞的功能复杂多样，在体内作为重要的“调节组织”维持着心血管系统的稳态。研究已证实，血管内皮功能障碍与多种心血管疾病的发生和发展密切相关，既是心血管疾病的早期表现形式，又是早于动脉硬化形态学发展的第一步。因此正确评价血管内皮功能，敏锐识别血管早期病变，可为早期诊断和防治高血压、冠心病、糖尿病等常见心血管疾病及其相应的靶器官病变提供有效手段，为临床治疗提供新方向。

在体血管内皮细胞处于一个复杂的力学环境中，作用于血管内皮细胞的机械应力有二种：血液流动产生的剪切应力(wallshearstress，WSS)和血压导致的周向应力(circumferentialstress，CS)。这二种力同时作用于血管内皮细胞，从而调节其表型。最近的研究指出，最容易出现动脉粥样硬化和内膜增生的循环区域，如弯曲血管内壁，分支血管外壁，也是血流剪切应力与血管壁周向应力最不同步的部位。为了说明这个不同相(asynchronous)的现象，美国生物医学工程学会前主席JohnM.Tarbell教授引入了应力相位角(stressphaseangle，SPA)的概念，应力相位角是血管壁周向应力与血流剪切应力之间的相位角。应力相位角反映压力和流体波形的不同步性，它是结合动脉壁力学因素和流体力学因素对血管内皮细胞影响的唯一参数。因此，如果能够明确应力相位角在检测血管病变部位中的作用，可以为血管早期病变提供新的诊断方法和依据。

目前的研究主要是通过计算机仿真，从细胞水平上验证了应力相位角理论的正确性。然而，尚无合适的方法对应力相位角进行定量分析。因此，为了研究弹性管腔内应力相位角的变化规律，必须有一种方法可以同步测量血管壁周向应力和血流剪切应力。

关于同步探测血管壁运动和血流动力学参数方法的研究国际上已有少量报道。RichardG.Wise等人利用相位对比磁共振成像技术同时测量大鼠心脏的流速分布图和心室壁的运动。S.Sampath等人利用磁化空间调制的磁共振成像技术测量健康志愿者心肌位移，利用梯度编码磁共振成像技术测量心肌血流速度。但是，磁共振设备价格昂贵、体积庞大、实时性差，使其应用受到了限制。而超声方法因其方便、实时等优点受到普遍应用。P.Tortoli等开发一种系统采用两个超声束来同时监测动脉的直径变化和流场速度分布，该系统的复杂性限制了它的广泛应用。H.Hasegawa等采用高帧频获得的射频回波信号同时测量动脉壁的径向应变和流场速度。但该方法计算流场速度依赖超声束与血流方向的夹角。J.Luo等应用斑点追踪方法同时测量动脉壁轴向速度和血流速度。然而，该方法要求较高的帧频(8kHz)。综上所述，目前文献中提到的磁共振成像无法满足测量实时性的要求，而报道的超声方法对成像系统有特殊的要求，都存在一定的局限性。另外，流道或流体的非透明性，使得目前应用较广的光学方法很难适用。

发明内容

本申请要解决的技术问题是针对现有技术的不足，提供一种实时、准确管壁应力相位角的测量方法。

本申请要解决的另一技术问题是提供一种基于上述方法的测量系统。

本申请要解决的技术问题通过以下技术方案加以解决：

一种管壁应力相位角的测量方法，包括：

在流体中加入超声造影微泡；

使用高频超声成像装置采集连续采集多帧感兴趣区域的超声图像；

对相邻两幅超声图像对应区域进行互相关分析，得到管壁的周向应力和流场剪切应力；

绘制管壁周向应力和流场剪切应力相对于时间变化的波形图，测量应力相位角。

其中所述管壁包括血管壁；所述流体包括血液；所述使用高频超声成像装置采集连续采集多帧感兴趣区域的超声图像包括：使用高频超声成像装置采集连续至少三个心动周期内的感兴趣区域的超声图像；所述绘制管壁周向应力和流场剪切应力相对于时间变化的波形图，测量应力相位角包括：绘制血管壁周向应力和血流剪切应力随心动周期变化的波形图，测量应力相位角。

所述使用高频超声成像装置采集连续采集多帧感兴趣区域的超声图像包括：使用高频超声成像装置采集连续采集获得N幅感兴趣区域的超声图像g⁽¹⁾(r,θ)、g⁽²⁾(r,θ)…g^(N)(r,θ)，其中(r,θ)为图像平面内的像素坐标；

所述对相邻两幅超声图像对应区域进行互相关分析，得到管壁的周向应力和流场剪切应力包括：

对于感兴趣区域的超声图像g⁽ⁿ⁾(r,θ)，n∈｛1，N｝；

n＝1时开始，将g⁽ⁿ⁾(r,θ)和g⁽ⁿ⁺¹⁾(r,θ)划分为小的网格区域，计算每个网格区域的平动位移；

用位移梯度估计管壁的旋转和形变以及高速度梯度流；

减小网格区域提高结果的空间分辨率；

对于每个网格区域进行计算，得到管壁位移和流体速度矢量分布图；

一直到n＝N－1重复以上过程；

通过弹性模量与位移梯度的乘积计算管壁的周向应力；

通过流体粘度与流场速度梯度的乘积计算流场剪切应力。

所述将g⁽ⁿ⁾(r,θ)和g⁽ⁿ⁺¹⁾(r,θ)划分为小的网格区域，计算每个网格区域的平动位移包括：将g⁽ⁿ⁾(r,θ)和g⁽ⁿ⁺¹⁾(r,θ)划分为小的网格区域，利用二维互相关算法结合亚像素算法和滤波插值算法计算每个网格区域的平动位移。

所述用位移梯度估计管壁的旋转和形变以及高速度梯度流包括：使用多次迭代算法用位移梯度估计管壁的旋转和形变以及高速度梯度流。

所述对于每个网格区域进行计算，得到管壁位移和流体速度矢量分布图之前，还包括：利用基于频域滤波和连续性方程的错误矢量剔除算法来提高结果的精确性。

一种管壁应力相位角的测量系统，包括超声成像模块和测量模块，

所述超声成像模块包括超声微泡和高频超声成像装置，所述超声微泡用于流体中的示踪；所述高频超声成像装置用于采集连续采集多帧感兴趣区域的超声图像；

所述测量模块包括互相关分析单元和应力相位角测量单元，所述互相关分析单元用于对相邻两幅超声图像对应区域进行互相关分析，得到管壁的周向应力和流场剪切应力；所述应力相位角测量单元用于绘制管壁周向应力和流场剪切应力相对于时间变化的波形图，测量应力相位角。

其中所述管壁包括血管壁；所述流体包括血液；所述超声成像装置还用于采集连续至少三个心动周期内的感兴趣区域的超声图像；所述应力相位角测量单元还用于绘制血管壁周向应力和血流剪切应力随心动周期变化的波形图，测量应力相位角。

其中所述超声成像装置还用于连续采集获得N幅感兴趣区域的超声图像g⁽¹⁾(r,θ)、g⁽²⁾(r,θ)…g^(N)(r,θ)，其中(r,θ)为图像平面内的像素坐标；

所述互相关分析单元还用于：

对于感兴趣区域的超声图像g⁽ⁿ⁾(r,θ)，n∈｛1，N｝；

n＝1时开始，将g⁽ⁿ⁾(r,θ)和g⁽ⁿ⁺¹⁾(r,θ)划分为小的网格区域，计算每个网格区域的平动位移；

用位移梯度估计管壁的旋转和形变以及高速度梯度流；

减小网格区域提高结果的空间分辨率；

对于每个网格区域进行计算，得到管壁位移和流体速度矢量分布图；

一直到n＝N－1重复以上过程；

通过弹性模量与位移梯度的乘积计算管壁的周向应力；

通过流体粘度与流场速度梯度的乘积计算流场剪切应力。

所述互相关分析单元还用于：将g⁽ⁿ⁾(r,θ)和g⁽ⁿ⁺¹⁾(r,θ)划分为小的网格区域，利用二维互相关算法结合亚像素算法和滤波插值算法计算每个网格区域的平动位移。

所述互相关分析单元还用于：使用多次迭代算法用位移梯度估计管壁的旋转和形变以及高速度梯度流。

所述互相关分析单元还用于：利用基于频域滤波和连续性方程的错误矢量剔除算法来提高结果的精确性。

由于采用了以上技术方案，使本申请具备的有益效果在于：

⑴在本申请的具体实施方式中，通过高频超声成像装置和超声微泡成像，对相邻帧图像进行互相关分析，能同步、实时、精确地获得周向应力和流场剪切应力相对于时间的波形图，从而测量应力相位角。

⑵在本申请的具体实施方式中，通过减少网格区域，提高了结果的空间分辨率。

⑶在本申请的具体实施方式中，利用基于频域滤波和连续性方程的错误矢量剔除算法来提高结果的精确性。

⑷在本申请的具体实施方式中，应用于血管壁的应力相位角测量，临床医生易掌握实施，且不易受其他因素影响。

⑸在本申请的具体实施方式中，该测量模块可作为一个图像后处理软件模块集成到现有的成像装置中，以提升成像系统的功能无需对现有临床成像系统进行硬件升级，升级成本低，容易被医院和医生接受，方便临床推广。

附图说明

图1为根据本申请管壁应力相位角的测量方法一个实施例的流程图；

图2为根据本申请管壁应力相位角的测量系统一个实施例的结构示意图。

具体实施方式

下面通过具体实施方式结合附图对本发明作进一步详细说明。

图1示出根据本申请管壁应力相位角的测量方法一个实施例的流程图，包括：

步骤102：在流体中加入超声造影微泡；

步骤104：使用高频超声成像装置采集连续采集多帧感兴趣区域的超声图像；

步骤106：对相邻两幅超声图像对应区域进行互相关分析，得到管壁的周向应力和流场剪切应力；

步骤108：绘制管壁周向应力和流场剪切应力相对于时间变化的波形图，测量应力相位角。

一种实施方式，其中管壁包括血管壁；流体包括血液；步骤102包括：使用高频超声成像装置采集连续三到五个心动周期内的感兴趣区域的超声图像；步骤108包括：绘制血管壁周向应力和血流剪切应力随心动周期变化的波形图，测量应力相位角。

一种实施方式，步骤104包括：使用高频超声成像装置采集连续采集获得N幅感兴趣区域的超声图像g⁽¹⁾(r,θ)、g⁽²⁾(r,θ)…g^(N)(r,θ)，其中(r,θ)为图像平面内的像素坐标；

一种实施方式，步骤106包括：

步骤S02：对于感兴趣区域的超声图像g⁽ⁿ⁾(r,θ)，n∈｛1，N｝；n＝1时开始，将g⁽ⁿ⁾(r,θ)和g⁽ⁿ⁺¹⁾(r,θ)划分为小的网格区域，计算每个网格区域的平动位移；

步骤S04：用位移梯度估计管壁的旋转和形变以及高速度梯度流；

步骤S06：减小网格区域提高结果的空间分辨率；

步骤S08：对于每个网格区域进行计算，得到管壁位移和流体速度矢量分布图；

步骤S10：一直到n＝N－1重复步骤S02到步骤S08的过程；

步骤S12：通过弹性模量与位移梯度的乘积计算管壁的周向应力；

步骤S14：通过流体粘度与流场速度梯度的乘积计算流场剪切应力。

一种实施方式，S02包括：将g⁽ⁿ⁾(r,θ)和g⁽ⁿ⁺¹⁾(r,θ)划分为小的网格区域，利用二维互相关算法结合亚像素算法和滤波插值算法计算每个网格区域的平动位移。

一种实施方式，步骤S04包括：使用多次迭代算法用位移梯度估计管壁的旋转和形变以及高速度梯度流。

一种实施方式，S06和S08之间还包括：

S07：利用基于频域滤波和连续性方程的错误矢量剔除算法来提高结果的精确性。

本实施例应用于血管壁应力相位角测量的更详细实施方式如下：

步骤202：在流体中加入示踪粒子—超声造影微泡，利用高频超声成像系统高帧速采集连续多帧感兴趣区域(血管壁和造影血流区域)的超声造影图像。

帧频（FR）为：FR=1/(BLD*T_t*N_ele)其中BLD为超声束线密度，T_t为超声声束的扫描时间，N_ele为激励状态的阵元数。另外，V_max=FR*W/4。W为选择的诊断窗口，若W一定，FR必须保证可以测量ROI的最大速度。“连续多帧”是3－5个心动周期内的图像帧。

步骤204：设f⁽¹⁾(r,θ)…f^(N)(r,θ)为N幅连续的超声造影图像，(r,θ)为图像平面内的像素坐标。在f⁽ⁿ⁾(r,θ)上选择感兴趣区域g⁽ⁿ⁾(r,θ)，当n=1时计算开始。随后在f⁽ⁿ⁺¹⁾(r,θ)上选择感兴趣区域g⁽ⁿ⁺¹⁾(r,θ)。每个感兴趣区域被划分为小的网格区域-诊断窗口。

首先将极坐标系变换为直角坐标系，利用二维互相关算法结合亚像素算法和滤波插值算法来计算每个诊断窗口的平动位移，用多次迭代算法、错误矢量剔除算法，提高测量精度，减小测量误差。

（1）二维标准互相关算法是利用图像的灰度分布相似性，来计算两个分析窗口（即k1在第1帧和第2帧中）的位移。其表达式为：

$R(p,q)=\frac{\underset{x=1}{\overset{M_1}{\mathrm{\Σ}}}\underset{y=1}{\overset{N_1}{\mathrm{\Σ}}}\left(f(x,y){\mathrm{\μ}}_f\right)\left(g(x+p,y+q){\mathrm{\μ}}_g\right)}{\sqrt{\underset{x=1}{\overset{M_1}{\mathrm{\Σ}}}\underset{y=1}{\overset{N_1}{\mathrm{\Σ}}}{\left(f(x,y){\mathrm{\μ}}_f\right)}^2\underset{i=1}{\overset{M_1}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_1}{\mathrm{\Σ}}}{\left(g(x+p,y+q){\mathrm{\μ}}_g\right)}^2}}$

其中f和g为两个分析窗口内的灰度值，μ_f和μ_g为对应的灰度平均值。

为了提高计算精度采用亚像素方法，可以精确到小数级像素，在该算法中采用了高斯峰值拟合公式：

$\begin{array}{c}\mathrm{\Δx}=i+\frac{\ln \left(R_{\mathrm{pq}}(i-1,j)\right)-\ln \left(R_{\mathrm{pq}}(i+1,j)\right)}{2\ln \left(R_{\mathrm{pq}}(i-1,j)\right)-4\ln \left(R_{\mathrm{pq}}(i,j)\right)+2\ln \left(R_{\mathrm{pq}}(i+1,j)\right)}\\ \mathrm{\Δy}=j+\frac{\ln \left(R_{\mathrm{pq}}(i,j-1)\right)-\ln \left(R_{\mathrm{pq}}(i,j+1)\right)}{2\ln \left(R_{\mathrm{pq}}(i,j-1)\right)-4\ln \left(R_{\mathrm{pq}}(i,j)\right)+2\ln \left(R_{\mathrm{pq}}(i,j+1)\right)}\end{array}$

R_pq(i,j)表示互相关函数的最大值，i和j表示互相关函数取得最大值时相应的坐标，R_pq(i-1,j),R_pq(i+1,j),R_pq(i,j-1)和R_pq(i,j+1)分别表示互相关函数数组中与R_pq(i,j)最邻近的周围四个网格点上的互相关函数的数值。

对于滤波与插值算法，本实施方式采用的是中值滤波和双线性插值，其具体公式如下：

$\begin{array}{c}U(x,y)>\mathrm{median}\left(U(x-1:x+1,y-1:y+1)\right)+\mathrm{threshold}*\mathrm{std}\left(U(x-1:x+1,y-1:y+1)\right)\\ U(x,y)<\mathrm{median}\left(U(x-1:x+1,y-1:y+1)\right)-\mathrm{threshold}*\mathrm{std}\left(U(x-1:x+1,y-1:y+1)\right)\end{array}$

U(x,y)＝(1-x)(1-y)U(x-1,y-1)+x(1-y)U(x+1,y-1)+(1-x)yU(x-1,y+1)+xyU(x+1,y+1)U=（u，v）为二维平动位移矢量。

（2）采用迭代算法通过二维平动位移的位移梯度进一步计算感兴趣区域的旋转和变形，得到几何变换的二维位移。

为了考虑感兴趣区域的旋转和变形，引入泰勒级数：

$U\left(r\right)=U\left(r_0\right)+{\left(\frac{\&PartialD;U}{\&PartialD;r}\right)}_{r=r_0}(r-r_0)+\frac{1}{2!}{\left(\frac{{\&PartialD;}^{2}U}{{\&PartialD;r}^2}\right)}_{r=r_0}{(r-r_0)}^2+o{(r-r_0)}^3$

r＝(x,y)。

迭代算法的具体过程为：首先，设定迭代次数K(一般为2－3次)。然后，对于第k次迭代，利用第k-1次的位移梯度U′_k-1来重置分析窗口的灰度值。

$f\left(r\right)=f(r-\frac{U_{k-1}}{2}-\frac{U_{k-1}^{\&prime;}r}{2})$

$g\left(r\right)=g(r+\frac{U_{k-1}}{2}+\frac{U_{k-1}^{\&prime;}r}{2})$

然后再用标准互相关公式计算互相关平面内位移最大值的位置R_max，因此最后计算的几何变换的位移为：

U_k＝U_k-1+R_max

（3）对几何变换的二维位移采用错误矢量剔除算法以提高精度，得到第1帧和第2帧之间的相对位移矢量。

错误矢量剔除算法是基于连续性方程的修正方法，具体为：首先根据最初的位移矢量，设置一个初始值vall，并通过下面公式计算每一点的值val：

$\mathrm{val}=\frac{\mathrm{\&Sigma;}|u_{x-1:x+1,y-1:y+1}-u_{x,y}|+\mathrm{\&Sigma;}|v_{x-1:x+1,y-1:y+1}-v_{x,y}|}{\mathrm{\&Sigma;}\left|u_{x-1:x+1,y-1:y+1}\right|+\mathrm{\&Sigma;}\left|v_{x-1:x+1,y-1:y+1}\right|}$

其中u_x,y和v_x,y应用高斯模板计算：

$\left\{\begin{array}{c}{u}_{x,y}=\frac{u_{x-1,y+1}+2u_{x,y+1}+u_{x+1,y+1}+{2u}_{x-1,y}+{2u}_{x+1,y}+u_{x-1,y-1}+{2u}_{x,y-1}+u_{x+1,y-1}}{12}\\ v_{x,y}=\frac{v_{x-1,y+1}+{2v}_{x,y+1}+v_{x+1,y+1}+{2v}_{x-1,y}+{2v}_{x+1,y}+v_{x-1,y-1}+{2v}_{x,y-1}+v_{x+1,y-1}}{12}\end{array}\right.$

如果val大于初始值vall，则这个值应用下面公式纠正：

$u_{x,y}=\frac{\mathrm{\&Sigma;}{u}_{x-1:x+1,y-1:y+1}}{8}$ $v_{x,y}=\frac{\mathrm{\&Sigma;}{v}_{x-1:x+1,y-1:y+1}}{8}$

由此，得到了高精度的相对位移矢量u和v。从而可以计算应变、应力以及血流动力学参数(血流量、血流速度、剪切应力)。

本实施方式通过多次迭代算法用位移梯度来估计血管壁的旋转和形变以及高速度梯度流

本实施方式减小诊断窗口利用互相关算法计算位移，提高结果的空间分辨率

对图像的每个诊断窗口进行计算，最后得到血管壁位移和血流速度矢量分布图。对n=2的f⁽ⁿ⁾(r,θ)和f⁽ⁿ⁺¹⁾(r,θ)重复以上过程，直到n=N-1。

步骤206：血管壁的周向应力通过弹性模量与位移梯度的乘积计算，流场剪切应力通过血液粘度与流场速度梯度的乘积计算。因此，可以绘制血管壁周向应力和血流剪切应力随心动周期的变化曲线图，从而测量应力相位角。

应力相位角是血管壁周向应力和血流剪切应力之间的相位角。研究指出，血管壁周向应力与直径的改变是同步的，因此，应力相位角可表示为：

φ(D-τ)＝φ(D-Q)-φ(τ-Q)

≈φ(P-Q)-φ(τ-Q)(1)

D是直径，Q是血流量，P是压力，τ是血流剪切应力。对于一个弹性血管，直径D和压力P几乎是同相位的，因此，φ(D-Q)近似等于阻力相位角，φ(P-Q)。阻力相位角由外端阻力、弹性和波反射决定。φ(τ-Q)是血流剪切应力和血流量之间的相位角(WQPA)，WQPA与血管的几何形状有密切关系。由此可见，应力相位角与血管几何形状、血管壁弹性密切相关。

本实施方式采用的计算公式，SPA=φ(θ-τ)。θ是周向应力，τ是血流剪切应力。

图2示出根据本申请管壁应力相位角的测量系统一个实施例的结构示意图，包括：包括超声成像模块和测量模块，超声成像模块包括超声微泡和高频超声成像装置，超声微泡用于流体中的示踪；高频超声成像装置用于采集连续采集多帧感兴趣区域的超声图像；测量模块进一步包括互相关分析单元和应力相位角测量单元。该互相关分析单元用于对相邻两幅超声图像对应区域进行互相关分析，得到管壁的周向应力和流场剪切应力；该应力相位角测量单元用于绘制管壁周向应力和流场剪切应力相对于时间变化的波形图，测量应力相位角。

一种实施方式，管壁包括血管壁；流体包括血液；超声成像装置还用于采集连续三－五个心动周期内的感兴趣区域的超声图像；应力相位角测量单元还用于绘制血管壁周向应力和血流剪切应力随心动周期变化的波形图，测量应力相位角。

一种实施方式，超声成像装置还用于连续采集获得N幅感兴趣区域的超声图像g⁽¹⁾(r,θ)、g⁽²⁾(r,θ)…g^(N)(r,θ)，其中(r,θ)为图像平面内的像素坐标；

该互相关分析单元还用于：

对于感兴趣区域的超声图像g⁽ⁿ⁾(r,θ)，n∈｛1，N｝；

n＝1时开始，将g⁽ⁿ⁾(r,θ)和g⁽ⁿ⁺¹⁾(r,θ)划分为小的网格区域，计算每个网格区域的平动位移；

用位移梯度估计管壁的旋转和形变以及高速度梯度流；

减小网格区域提高结果的空间分辨率；

对于每个网格区域进行计算，得到管壁位移和流体速度矢量分布图；

一直到n＝N－1重复以上过程；

通过弹性模量与位移梯度的乘积计算管壁的周向应力；

通过流体粘度与流场速度梯度的乘积计算流场剪切应力。

一种实施方式，互相关分析单元还用于：使用多次迭代算法用位移梯度估计管壁的旋转和形变以及高速度梯度流。

一种实施方式，互相关分析单元还用于：使用多次迭代算法用位移梯度估计管壁的旋转和形变以及高速度梯度流。

一种实施方式，互相关分析单元还用于：利用基于频域滤波和连续性方程的错误矢量剔除算法来提高结果的精确性。

采用本实施例的测量系统进行血管壁应力相位角测量的更详细实施方式如下：

超声成像模块用于在流体中加入示踪粒子—超声造影微泡，利用高频超声成像系统高帧速采集连续多帧感兴趣区域(血管壁和造影血流区域)的超声造影图像。

帧频（FR）为：FR=1/(BLD*T_t*N_ele)其中BLD为超声束线密度，T_t为超声声束的扫描时间，N_ele为激励状态的阵元数。另外，V_max=FR*W/4。W为选择的诊断窗口，若W一定，FR必须保证可以测量ROI的最大速度。“连续多帧”是3－5个心动周期内的图像帧。

在测量模块中，设f⁽¹⁾(r,θ)…f^(N)(r,θ)为N幅连续的超声造影图像，(r,θ)为图像平面内的像素坐标。在f⁽ⁿ⁾(r,θ)上选择感兴趣区域g⁽ⁿ⁾(r,θ)，当n=1时计算开始。随后在f⁽ⁿ⁺¹⁾(r,θ)上选择感兴趣区域g⁽ⁿ⁺¹⁾(r,θ)。每个感兴趣区域被划分为小的网格区域-诊断窗口。

本实施例中，互相关分析单元首先将极坐标系变换为直角坐标系，利用二维互相关算法结合亚像素算法和滤波插值算法来计算每个诊断窗口的平动位移，用多次迭代算法、错误矢量剔除算法，提高测量精度，减小测量误差。

（1）互相关分析单元的二维标准互相关算法是利用图像的灰度分布相似性，来计算两个分析窗口（即k1在第1帧和第2帧中）的位移。其表达式为：

$R(p,q)=\frac{\underset{x=1}{\overset{M_1}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{N_1}{\mathrm{\&Sigma;}}}\left(f(x,y){\mathrm{\&mu;}}_f\right)\left(g(x+p,y+q){\mathrm{\&mu;}}_g\right)}{\sqrt{\underset{x=1}{\overset{M_1}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{\left(f(x,y){\mathrm{\&mu;}}_f\right)}^2\underset{i=1}{\overset{M_1}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{\left(g(x+p,y+q){\mathrm{\&mu;}}_g\right)}^2}}$

其中f和g为两个分析窗口内的灰度值，μ_f和μ_g为对应的灰度平均值。

为了提高计算精度采用亚像素方法，可以精确到小数级像素，在该算法中采用了高斯峰值拟合公式：

$\begin{array}{c}\mathrm{\&Delta;x}=i+\frac{\ln \left(R_{\mathrm{pq}}(i-1,j)\right)-\ln \left(R_{\mathrm{pq}}(i+1,j)\right)}{2\ln \left(R_{\mathrm{pq}}(i-1,j)\right)-4\ln \left(R_{\mathrm{pq}}(i,j)\right)+2\ln \left(R_{\mathrm{pq}}(i+1,j)\right)}\\ \mathrm{\&Delta;y}=j+\frac{\ln \left(R_{\mathrm{pq}}(i,j-1)\right)-\ln \left(R_{\mathrm{pq}}(i,j+1)\right)}{2\ln \left(R_{\mathrm{pq}}(i,j-1)\right)-4\ln \left(R_{\mathrm{pq}}(i,j)\right)+2\ln \left(R_{\mathrm{pq}}(i,j+1)\right)}\end{array}$

R_pq(i,j)表示互相关函数的最大值，i和j表示互相关函数取得最大值时相应的坐标，R_pq(i-1,j),R_pq(i+1,j),R_pq(i,j-1)和R_pq(i,j+1)分别表示互相关函数数组中与R_pq(i,j)最邻近的周围四个网格点上的互相关函数的数值。

对于滤波与插值算法，本实施例采用的是中值滤波和双线性插值，其具体公式如下：

$\begin{array}{c}U(x,y)>\mathrm{median}\left(U(x-1:x+1,y-1:y+1)\right)+\mathrm{threshold}*\mathrm{std}\left(U(x-1:x+1,y-1:y+1)\right)\\ U(x,y)<\mathrm{median}\left(U(x-1:x+1,y-1:y+1)\right)-\mathrm{threshold}*\mathrm{std}\left(U(x-1:x+1,y-1:y+1)\right)\end{array}$

U(x,y)＝(1-x)(1-y)U(x-1,y-1)+x(1-y)U(x+1,y-1)+(1-x)yU(x-1,y+1)+xyU(x+1,y+1)U=（u，v）为二维平动位移矢量。

（2）采用迭代算法通过二维平动位移的位移梯度进一步计算感兴趣区域的旋转和变形，得到几何变换的二维位移。

为了考虑感兴趣区域的旋转和变形，引入泰勒级数：

$U\left(r\right)=U\left(r_0\right)+{\left(\frac{\&PartialD;U}{\&PartialD;r}\right)}_{r=r_0}(r-r_0)+\frac{1}{2!}{\left(\frac{{\&PartialD;}^{2}U}{{\&PartialD;r}^2}\right)}_{r=r_0}{(r-r_0)}^2+o{(r-r_0)}^3$

r＝(x,y)。

迭代算法的具体过程为：首先，设定迭代次数K(一般为2－3次)。然后，对于第k次迭代，利用第k-1次的位移梯度U′k-1来重置分析窗口的灰度值。

$f\left(r\right)=f(r-\frac{U_{k-1}}{2}-\frac{U_{k-1}^{\&prime;}r}{2})$

$g\left(r\right)=g(r+\frac{U_{k-1}}{2}+\frac{U_{k-1}^{\&prime;}r}{2})$

然后再用标准互相关公式计算互相关平面内位移最大值的位置R_max，因此最后计算的几何变换的位移为：

U_k＝U_k-1+R_max

（3）本实施例中，对几何变换的二维位移采用错误矢量剔除算法以提高精度，得到第1帧和第2帧之间的相对位移矢量。

错误矢量剔除算法是基于连续性方程的修正方法，具体为：首先根据最初的位移矢量，设置一个初始值vall，并通过下面公式计算每一点的值val：

$\mathrm{val}=\frac{\mathrm{\&Sigma;}|u_{x-1:x+1,y-1:y+1}-u_{x,y}|+\mathrm{\&Sigma;}|v_{x-1:x+1,y-1:y+1}-v_{x,y}|}{\mathrm{\&Sigma;}\left|u_{x-1:x+1,y-1:y+1}\right|+\mathrm{\&Sigma;}\left|v_{x-1:x+1,y-1:y+1}\right|}$

其中u_x,y和v_x,y应用高斯模板计算：

$\left\{\begin{array}{c}{u}_{x,y}=\frac{u_{x-1,y+1}+2u_{x,y+1}+u_{x+1,y+1}+{2u}_{x-1,y}+{2u}_{x+1,y}+u_{x-1,y-1}+{2u}_{x,y-1}+u_{x+1,y-1}}{12}\\ v_{x,y}=\frac{v_{x-1,y+1}+{2v}_{x,y+1}+v_{x+1,y+1}+{2v}_{x-1,y}+{2v}_{x+1,y}+v_{x-1,y-1}+{2v}_{x,y-1}+v_{x+1,y-1}}{12}\end{array}\right.$

如果val大于初始值vall，则这个值应用下面公式纠正：

$u_{x,y}=\frac{\mathrm{\&Sigma;}{u}_{x-1:x+1,y-1:y+1}}{8}$ $v_{x,y}=\frac{\mathrm{\&Sigma;}{v}_{x-1:x+1,y-1:y+1}}{8}$

由此，得到了高精度的相对位移矢量u和v。从而可以计算应变、应力以及血流动力学参数(血流量、血流速度、剪切应力)。

本实施例通过多次迭代算法用位移梯度来估计血管壁的旋转和形变以及高速度梯度流

本实施例减小诊断窗口利用互相关算法计算位移，提高结果的空间分辨率。

对图像的每个诊断窗口进行计算，最后得到血管壁位移和血流速度矢量分布图。对n=2的f⁽ⁿ⁾(r,θ)和f⁽ⁿ⁺¹⁾(r,θ)重复以上过程，直到n=N-1。

在本实施例的应力相位角测量单元中，血管壁的周向应力通过弹性模量与位移梯度的乘积计算，流场剪切应力通过血液粘度与流场速度梯度的乘积计算。因此，可以绘制血管壁周向应力和血流剪切应力随心动周期的变化曲线图，从而测量应力相位角。

应力相位角是血管壁周向应力和血流剪切应力之间的相位角。研究指出，血管壁周向应力与直径的改变是同步的，因此，应力相位角可表示为：

φ(D-τ)＝φ(D-Q)-φ(τ-Q)

≈φ(P-Q)-φ(τ-Q)(1)

D是直径，Q是血流量，P是压力，τ是血流剪切应力。对于一个弹性血管，直径D和压力P几乎是同相位的，因此，φ(D-Q)近似等于阻力相位角，φ(P-Q)。阻力相位角由外端阻力、弹性和波反射决定。φ(τ-Q)是血流剪切应力和血流量之间的相位角(WQPA)，WQPA与血管的几何形状有密切关系。由此可见，应力相位角与血管几何形状、血管壁弹性密切相关。

本实施例采用的计算公式，SPA=φ(θ-τ)。θ是周向应力，τ是血流剪切应力。

利用上述实施例的方法和系统可以进行应力相位角的测量研究：

（1）血液仿体制备：直径为5μm、10μm、或20μm的微小尼龙粒子作为散射子。通过磁力搅拌器混合适量的尼龙粒子和适当重量的液体(纯水85.41%，纯甘油10.25%，右旋糖酐3.42%，表面活性剂0.92%)。用滤网过滤血液仿体中的残余小块。通过减压除去血液仿体中的气泡。水和甘油的比例决定血液仿体的声速和密度。

（2）弹性血管仿体制备：将一定质量分数的聚乙烯醇(polyvinylalcohol,PVA)粉末溶于水，添加一定质量分数的散射子，加热混合均匀。将PVA、散射子的水溶液灌注到模具的空隙，然后放入冷冻装置，使其经历一系列的“冷冻-解冻”周期，最终形成固态、弹性凝胶PVA-c，从模具中取出，即为弹性血管仿体。通过调控模具的几何形状、PVA质量分数、散射子质量分数、“冷冻-解冻”次数，可以制备不同几何形状(直管、弯曲、分叉，狭窄)、不同弹性(弹性模量50kPa～800kPa)、声速和声衰减系数的血管仿体。利用机械压缩方法和超声辐射力弹性成像方法测量不同“冷冻-解冻”周期血管仿体的弹性模量和剪切模量；脉冲声波插入法测量仿体的声速和声衰减系数。仿真血液循环系统搭建：脉动循环系统由Harvard的血流脉动泵(55-3305,美国)、压力传感器、弹性血管仿体、水槽、高频超声成像系统组成。脉动泵产生脉动流，经过传输管道到达血管仿体，使血管仿体产生周期性的收缩舒张，调节脉动泵的泵血量、频率、收缩期/舒张期时间比，使其模仿人体血管脉动和血流运动。基于该体外仿真血液循环系统开展实验，研究血管几何形状和血管壁弹性对应力相位角的影响。

（3）研究血管几何形状和血管壁弹性改变时应力相位角的演变规律：

体外脉动管流实验：在流体循环系统中采用不同几何形状(直管、弯曲、分叉，狭窄)，相同弹性的血管仿体，保证血管仿体入口处流量或压力恒定，保持脉动泵频率、每搏输出量、“收缩期/舒张期”时间比等参数不变，利用高分辨率超声成像方法获得血管壁周向应力与血流剪切应力随心动周期的变化曲线，测量应力相位角，研究血管几何形状改变时应力相位角的演变规律。然后，对不同弹性(弹性模量50kPa～800kPa)，相同几何形状的血管仿体进行实验，研究血管壁弹性改变时应力相位角的演变规律。

活体动物实验：初步的动物实验中，将对常规饲料喂养的正常小鼠和高脂饲料喂养的ApoE基因敲除小鼠的颈总动脉，颈内动脉和颈外动脉进行分析，比较组间相同位置及组内不同位置处周向应力、剪切应力、应力相位角的变化。ApoE基因敲除小鼠的结果也将与“金标准”病理切片对比，验证应力相位角在血管病变的检测和定位中的作用。

现有技术主要通过计算机仿真，从细胞水平上验证了应力相位角理论的正确性。然而，尚无合适的方法对应力相位角进行定量分析。因此，为了研究弹性管腔内应力相位角的变化规律，必须有一种方法可以同步测量血管壁周向应力和血流剪切应力。本申请的测量方法能同步精确获得血管壁周向应力和血流剪切应力随心动周期变化的波形图，测量应力相位角，可以此深入研究血管几何形状和血管壁弹性改变时应力相位角的演变规律，为心血管疾病早期诊断提供新的手段。临床医生以及具有一定成像设备使用经验的操作者都可以采集到所需的图像，不易受其他因素影响。而且该方法可作为一个图像后处理软件模块集成到现有的成像系统中，以提升成像系统的功能。由于无需对现有临床成像系统进行硬件升级，升级成本低，容易被医院和医生接受，方便临床推广。

以上内容是结合具体的实施方式对本申请所作的进一步详细说明，不能认定本申请的具体实施只局限于这些说明。对于本申请所属技术领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干简单推演或替换。

Claims (10)

1.一种管壁应力相位角的测量方法，其特征在于，包括：

在流体中加入超声造影微泡；

使用高频超声成像装置采集连续采集多帧感兴趣区域的超声图像；

对相邻两幅超声图像对应区域进行互相关分析，得到管壁的周向应力和流场剪切应力；

绘制管壁周向应力和流场剪切应力相对于时间变化的波形图，测量应力相位角；

所述使用高频超声成像装置采集连续采集多帧感兴趣区域的超声图像包括：

使用高频超声成像装置采集连续采集获得N幅感兴趣区域的超声图像g⁽¹⁾(r,θ)、g⁽²⁾(r,θ)…g^(N)(r,θ)，其中(r,θ)为图像平面内的像素坐标；

所述对相邻两幅超声图像对应区域进行互相关分析，得到管壁的周向应力和流场剪切应力包括：

对于感兴趣区域的超声图像g⁽ⁿ⁾(r,θ)，n∈{1，N}；

n＝1时开始，将g⁽ⁿ⁾(r,θ)和g⁽ⁿ⁺¹⁾(r,θ)划分为小的网格区域，计算每个网格区域的平动位移；

用位移梯度估计管壁的旋转和形变以及高速度梯度流；

减小网格区域提高结果的空间分辨率；

对于每个网格区域进行计算，得到管壁位移和流体速度矢量分布图；

重复以上过程，直到n＝N－1；

通过弹性模量与位移梯度的乘积计算管壁的周向应力；

通过流体粘度与流场速度梯度的乘积计算流场剪切应力。

2.如权利要求1所述的方法，其特征在于，所述使用高频超声成像装置采集连续采集多帧感兴趣区域的超声图像包括：使用高频超声成像装置采集连续至少三个心动周期内的感兴趣区域的超声图像；所述绘制管壁周向应力和流场剪切应力相对于时间变化的波形图，测量应力相位角包括：绘制管壁周向应力和流场剪切应力随心动周期变化的波形图，测量应力相位角。

3.如权利要求1所述的方法，其特征在于，所述将g⁽ⁿ⁾(r,θ)和g⁽ⁿ⁺¹⁾(r,θ)划分为小的网格区域，计算每个网格区域的平动位移包括：

将g⁽ⁿ⁾(r,θ)和g⁽ⁿ⁺¹⁾(r,θ)划分为小的网格区域，利用二维互相关算法结合亚像素算法和滤波插值算法计算每个网格区域的平动位移。

4.如权利要求1所述的方法，其特征在于，所述用位移梯度估计管壁的旋转和形变以及高速度梯度流包括：

使用多次迭代算法用位移梯度估计管壁的旋转和形变以及高速度梯度流。

5.如权利要求1至4中任一项所述的方法，其特征在于，所述对于每个网格区域进行计算，得到管壁位移和流体速度矢量分布图之前，还包括：

利用基于频域滤波和连续性方程的错误矢量剔除算法来提高结果的精确性。

6.一种管壁应力相位角的测量系统，其特征在于，包括超声成像模块和测量模块，

所述超声成像模块包括超声微泡和高频超声成像装置，所述超声微泡用于流体中的示踪；所述高频超声成像装置用于采集连续采集多帧感兴趣区域的超声图像；

所述测量模块包括互相关分析单元和应力相位角测量单元，所述互相关分析单元用于对相邻两幅超声图像对应区域进行互相关分析，得到管壁的周向应力和流场剪切应力；所述应力相位角测量单元用于绘制管壁周向应力和流场剪切应力相对于时间变化的波形图，测量应力相位角；

其中所述超声成像装置还用于连续采集获得N幅感兴趣区域的超声图像g⁽¹⁾(r,θ)、g⁽²⁾(r,θ)…g^(N)(r,θ)，其中(r,θ)为图像平面内的像素坐标；

所述互相关分析单元还用于：

对于感兴趣区域的超声图像g⁽ⁿ⁾(r,θ)，n∈{1，N}；

n＝1时开始，将g⁽ⁿ⁾(r,θ)和g⁽ⁿ⁺¹⁾(r,θ)划分为小的网格区域，计算每个网格区域的平动位移；

用位移梯度估计管壁的旋转和形变以及高速度梯度流；

减小网格区域提高结果的空间分辨率；

对于每个网格区域进行计算，得到管壁位移和流体速度矢量分布图；

重复以上过程，直到n＝N－1；

通过弹性模量与位移梯度的乘积计算管壁的周向应力；

通过流体粘度与流场速度梯度的乘积计算流场剪切应力。

7.如权利要求6所述的系统，其特征在于，所述超声成像装置还用于采集连续至少三个心动周期内的感兴趣区域的超声图像；所述应力相位角测量单元还用于绘制管壁周向应力和流场剪切应力随心动周期变化的波形图，测量应力相位角。

8.如权利要求6所述的系统，其特征在于，所述互相关分析单元还用于：将g⁽ⁿ⁾(r,θ)和g⁽ⁿ⁺¹⁾(r,θ)划分为小的网格区域，利用二维互相关算法结合亚像素算法和滤波插值算法计算每个网格区域的平动位移。

9.如权利要求6所述的系统，其特征在于，所述互相关分析单元还用于：使用多次迭代算法用位移梯度估计管壁的旋转和形变以及高速度梯度流。

10.如权利要求6至9中任一项所述的系统，其特征在于，所述互相关分析单元还用于：利用基于频域滤波和连续性方程的错误矢量剔除算法来提高结果的精确性。