一种超声心动图粒子图像测速速度场修正方法

摘要

一种超声心动图粒子图像测速速度场修正方法，它涉及速度场修正方法。本发明解决了解决现有超声医学图像进行噪声抑制的方法中空域滤波法会引起图像模糊和细节丢失、小波变换法对超声图像的乘性噪声效果不佳及各向异性扩散降噪法算法难度大的技术问题。方法：1、读入t时刻的超声心动图的初始速度场；2、计算各速度矢量对应的协方差矩阵；3、对各速度矢量进行各向异性扩散；4、各速度矢量进行中值滤波；5、各速度矢量移流计算；6、判断是否是最后的循环；如果是，则执行步骤9的结束；如果不是，则执行步骤7的计算下一循环的初始速度场；步骤8、令循环次数加1后，返回至步骤2；本方法能保持速度场细节，可对各种粒子图像测速场进行处理。

说明书

技术领域

本发明涉及速度场修正方法。

背景技术

超声心动图粒子图像测速技术是一种新兴的医学诊断方法，应用于对心腔内血流流场 的成像观察。由于超声成像的物理原因，相干波叠加形成的散斑噪声严重影响了成像质量， 尤其是使速度场的提取变得十分困难，降低了血流速度场的定量分析精度，限制了粒子图 像测速在超声心动图临床诊断中的进一步应用。对超声心动图粒子图像测速速度场进行修 正处理，可以为定量分析等后续处理提供更有利的条件，提高诊断的准确性。因此，研究 对超声心动图粒子图像测速速度场的修正，有助于提高临床诊断的准确性，具有重要的意 义。

超声心动图粒子图像测速的基本原理是基于灰度互相关的匹配算法。其应用的主要弊 端是在超声图像特有的乘性噪声下导致匹配错误率很高，从而导致匹配后得到的速度矢量 场杂乱无章，难以探寻其中的规律。从目前在该领域内对超声医学图像进行噪声抑制的方 法来看，主要有空域滤波、小波变换和各向异性扩散算法。空域滤波主要指邻域平均和中 值滤波等，本质是根据像素点邻域窗口信息计算出一个新值来取代原值，在去噪同时会引 起图像模糊和细节丢失。小波变换法是将超声图像变换到小波域，通过阈值法去除某些小 波系数，再逆变换以去除噪声，但该方法对加性噪声更有效，对于超声图像的乘性噪声效 果不佳。各向异性扩散降噪使用选择性扩散方式，在图像的噪声处有较大的扩散系数，有 利于图像的平滑，而在边缘处有较小的扩散系数，保持了图像的细节，但超声图像中，由 于乘性噪声的影响，难以从灰度变化上准确的区分图像的边缘和噪声，加大了算法的难度。

发明内容

本发明是为了解决现有超声医学图像进行噪声抑制的方法中空域滤波法会引起图像模 糊和细节丢失、小波变换法对超声图像的乘性噪声效果不佳及各向异性扩散降噪法算法难 度大的技术问题，而提供一种超声心动图粒子图像测速速度场修正方法。

本发明的一种超声心动图粒子图像测速速度场修正方法按以下步骤进行：

步骤一：用超声仪器检测超声心动图，然后读入用粒子图像测速法(Particle Image  Velocimetry)测量的t时刻的超声心动图的初始速度场该初始速度场u₁由m×n个网 格组成，每个网格中有一个速度矢量其中1≤i≤m，1≤j≤n；用k表示循环次数，K 为循环次数的最大值，K为3～5；初始化k＝1；

步骤二：计算初始速度场中2≤i≤m-1且2≤j≤n-1的各个速度矢量对应的协 方差矩阵S_c，其中，$S_c=\left(\begin{array}{cc}{\psi }_{\mathrm{xx}}& {\psi }_{\mathrm{xy}}\\ {\psi }_{\mathrm{yx}}& {\psi }_{\mathrm{yy}}\end{array}\right),$其${\psi }_{\mathrm{xx}}=\frac{{\partial }^2{u}_x}{{\partial x}^2},$${\psi }_{\mathrm{yy}}=\frac{{\partial }^2{u}_y}{{\partial x}^2},$${\psi }_{\mathrm{xy}}=\frac{{\partial }^2{u}_x}{\partial x\partial y},$${\psi }_{\mathrm{yx}}=\frac{{\partial }^2{u}_y}{\partial y\partial x};$u_x为速度矢量在x轴方向上的分量，u_y为速度矢量在y轴方向上的分量；

步骤三：以步骤二中得到的协方差矩阵作为扩散核，对初始速度场中2≤i≤m-1且 2≤j≤n-1的各个速度矢量分别按

$u_2^k(i,j)=\frac{\underset{-1\le \beta \le 1}{\underset{-1\le \alpha \le 1}{\Sigma }}{u}_1^k(i+\alpha ,j+\beta )\exp [{-(\alpha ,\beta )}^T·S_c^{-1}·(\alpha ,\beta )]}{\underset{-1\le \beta \le 1}{\underset{-1\le \alpha \le 1}{\Sigma }}\exp [{-(\alpha ,\beta )}^T·S_c^{-1}·(\alpha ,\beta )]}$进行各向异性扩散， 得到速度场α、β为整数；S_c指步骤二中计算得到的对应的扩散核，即协方 差矩阵；为S_c的逆阵；(α，β)^T为向量(α，β)的转置；

步骤四、对步骤三得到的速度场中各速度矢量实施中值滤波，中值滤波的 范围是当前计算点的3×3邻域，得到速度场

步骤五：将步骤四的速度场中各速度矢量进行移流计算，反推出Δt时刻前 该速度矢量所在的位置以及在该位置的速度值再用这个速度矢量值替代得到末态速度矢量从而得到末态速度场其中Δt为时间间隔，Δt取1/300秒 ～1/30000秒；

步骤六：判断k＝K是否成立，如果否，执行步骤七；如果是，执行步骤九；

步骤七：按计算第k+1次循环的初始速度场ε为比例系 数，取0.05～0.15；

步骤八；令k加1，执行步骤二；

步骤九：结束。

本发明吸收了光流场处理中Navier-Stokes方程的优点，同时继承了各向异性扩散能 够保持图像细节的特点，不会引起图像模糊和细节丢失，而且排除了超声图像固有的噪声 影响，可以从定性和定量上更加准确地分析流场的整体趋势和具体参数指标，算法简单。 本发明可广泛适用于各种需要对粒子图像测速速度场进行处理的场合。

附图说明

图1是具体实施方式一的流程图；图2是具体实施方式三中称移计算的图解图；图3 是具体实施方式三中双线性插值法图解图；图4是具体实施方式四中步骤二中速度场的网 格示意图；图5是具体实施方式四中t时刻的超声心动图的血流速度初始速度场图；图6 是具体实施方式四中经修正后的血流速度速度场图；图7是具体实施方式五中经修正后的 血流速度速度场图。

具体实施方式

具体实施方式一：(请参考附图1)本实施方式的一种超声心动图粒子图像测速速度场 修正方法按以下步骤进行：

步骤一：用超声仪器检测超声心动图，然后读入用粒子图像测速法(Particle Image  Velocimetry)测量的t时刻的超声心动图的初始速度场该初始速度场u₁由m×n个网 格组成，每个网格中有一个速度矢量其中1≤i≤m，1≤j≤n；用k表示循环次数，K 为循环次数的最大值，K为3～5；初始化k＝1；

步骤二：计算初始速度场中2≤i≤m-1且2≤j≤n-1的各个速度矢量对应的协 方差矩阵S_c，其中，$S_c=\left(\begin{array}{cc}{\psi }_{\mathrm{xx}}& {\psi }_{\mathrm{xy}}\\ {\psi }_{\mathrm{yx}}& {\psi }_{\mathrm{yy}}\end{array}\right),$其${\psi }_{\mathrm{xx}}=\frac{{\partial }^2{u}_x}{{\partial x}^2},$${\psi }_{\mathrm{yy}}=\frac{{\partial }^2{u}_y}{{\partial x}^2},$${\psi }_{\mathrm{xy}}=\frac{{\partial }^2{u}_x}{\partial x\partial y},$${\psi }_{\mathrm{yx}}=\frac{{\partial }^2{u}_y}{\partial y\partial x};$u_x为速度矢量在x轴方向上的分量，u_y为速度矢量在y轴方向上的分量；

步骤三：以步骤二中得到的协方差矩阵作为扩散核，对初始速度场中2≤i≤m-1且 2≤j≤n-1的各个速度矢量分别按

$u_2^k(i,j)=\frac{\underset{-1\le \beta \le 1}{\underset{-1\le \alpha \le 1}{\Sigma }}{u}_1^k(i+\alpha ,j+\beta )\exp [{-(\alpha ,\beta )}^T·S_c^{-1}·(\alpha ,\beta )]}{\underset{-1\le \beta \le 1}{\underset{-1\le \alpha \le 1}{\Sigma }}\exp [{-(\alpha ,\beta )}^T·S_c^{-1}·(\alpha ,\beta )]}$进行各向异性扩散， 得到速度场α、β为整数；S_c指步骤二中计算得到的对应的扩散核，即协方 差矩阵；为S_c的逆阵；(α，β)^T为向量(α，β)的转置；

步骤四、对步骤三得到的速度场中各速度矢量实施中值滤波，中值滤波的 范围是当前计算点的3×3邻域，得到速度场

步骤五：将步骤四的速度场中各速度矢量进行移流计算，反推出Δt时刻前 该速度矢量所在的位置以及在该位置的速度值，再用这个速度矢量值替代得到 末态速度矢量从而得到末态速度场其中Δt为时间间隔，Δt取1/300秒 ～1/30000秒；

步骤六：判断k＝K是否成立，如果否，执行步骤七；如果是，执行步骤九；

步骤七：按计算第k+1次循环的初始速度场ε为比例系 数，取0.05～0.15；

步骤八；令k加1，执行步骤二；

步骤九：结束。

本实施方式吸收了光流场处理中Navier-Stokes方程的优点，同时继承了各向异性扩 散能够保持图像细节的特点，不会引起图像模糊和细节丢失，而且排除了超声图像固有的 噪声影响，可以从定性和定量上更加准确地分析流场的整体趋势和具体参数指标，算法简 单。本实施方式可广泛适用于各种需要对超声粒子图像测速进行处理的场合。

具体实施方式二：本实施方式与具体实施方式一不同的是；步骤四中中值滤波时速度 矢量的排序规则为边缘中值法或矢量中值法；其它与具体实施方式一相同。

本实施方式中边缘中值法是将速度矢量的各分量分别比较大小排序后，取各分量的中 值得到的新矢量；而矢量中值法则是计算矢量与其3×3邻域内各矢量的欧式距离，根据 欧式距离的大小进行排序。进行中值滤波是为了减小步骤三中对速度矢量实施各向异性扩 散带来的局部误差。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤五中移流计算 的方法如下：

移流计算的图解图如图2所示，当前时刻点a网格的速度矢量为计算出在Δt 时刻前该速度矢量所在的位置为b点，用双线性插值法计算b点处的速度矢量的值 再将赋给a网格中的速度矢量，即

双线性插值法计算b点处的速度矢量的值的方法如下：

双线性插值法图解如图3所示，将b点周围的四个速度矢量u_A、u_B、u_C和u_D所在 的点联成一个正方形，四个顶点分别用A、B、C和D表示，其中f_x＝AP/AB，f_y＝FD/BD， 则$u_3^k{(i,j)}^{\prime }=(1-f_x)f_y·u_A+f_x{f}_y·u_B+(1-f_x)(1-f_y)·u_C+f_x(1-f_y)·u_D;$其中为b点处的速度矢量的值。其它与具体实施方式一或二相同。

本实施方式的移流是为了使速度矢量值更精确。

具体实施方式四：本实施方式的一种超声心动图粒子图像测速速度场修正方法按以下 步骤进行：

步骤一：用超声仪器检测超声心动图，然后读入用粒子图像测速法(Particle Image  Velocimetry)测量的t时刻的超声心动图的初始速度场该初始速度场u₁由m×n个网 格组成，每个网格中有一个速度矢量其中，m＝40，n＝50，1≤i≤m，1≤j≤n；用k表 示循环次数，K为循环次数的最大值，K为3；初始化k＝1；

步骤二：计算初始速度场中2≤i≤39且2≤j≤49的各个速度矢量对应的协方 差矩阵S_c，其中$S_c=\left(\begin{array}{cc}{\psi }_{\mathrm{xx}}& {\psi }_{\mathrm{xy}}\\ {\psi }_{\mathrm{yx}}& {\psi }_{\mathrm{yy}}\end{array}\right),$其${\psi }_{\mathrm{xx}}=\frac{{\partial }^2{u}_x}{{\partial x}^2},$${\psi }_{\mathrm{yy}}=\frac{{\partial }^2{u}_y}{{\partial x}^2},$${\psi }_{\mathrm{xy}}=\frac{{\partial }^2{u}_x}{\partial x\partial y},$${\psi }_{\mathrm{yx}}=\frac{{\partial }^2{u}_y}{\partial y\partial x};$u_x为速度矢量在x轴方向上的分量，u_y为速度矢量在y轴方向上的分量；

步骤三：以步骤二中得到的协方差矩阵作为扩散核，对初始速度场中2≤i≤39且 2≤j≤49的各个速度矢量分别按

$u_2^k(i,j)=\frac{\underset{-1\le \beta \le 1}{\underset{-1\le \alpha \le 1}{\Sigma }}{u}_1^k(i+\alpha ,j+\beta )\exp [{-(\alpha ,\beta )}^T·S_c^{-1}·(\alpha ,\beta )]}{\underset{-1\le \beta \le 1}{\underset{-1\le \alpha \le 1}{\Sigma }}\exp [{-(\alpha ,\beta )}^T·S_c^{-1}·(\alpha ,\beta )]}$进行各向异性扩散， 得到速度场α、β为整数；S_c指步骤二中计算得到的对应的扩散核，即协方 差矩阵；为S_c的逆阵；(α，β)^T为向量(α，β)的转置；

步骤四、对步骤三得到的速度场中各速度矢量实施中值滤波，中值滤波的 范围是当前计算点的3×3邻域，得到速度场

步骤五：将步骤四的速度场中各速度矢量进行移流计算，反推出Δt时刻前 该速度矢量所在的位置以及在该位置的速度值再用这个速度矢量值替代得到末态速度矢量从而得到末态速度场其中Δt为时间间隔，Δt取1/300秒；

步骤六：判断k＝K是否成立，如果否，执行步骤七；如果是，执行步骤九；

步骤七：按计算第k+1次循环的初始速度场ε为比例系 数，ε取0.10；

步骤八；令k加1，执行步骤二；

步骤九：结束。

本实施方式中步骤二中速度场的网格如图4所示，${\psi }_{\mathrm{xx}}=\frac{{\partial }^2{u}_x}{{\partial x}^2},$${\psi }_{\mathrm{yy}}=\frac{{\partial }^2{u}_y}{{\partial x}^2},$${\psi }_{\mathrm{xy}}=\frac{{\partial }^2{u}_x}{\partial x\partial y},$${\psi }_{\mathrm{yx}}=\frac{{\partial }^2{u}_y}{\partial y\partial x}$的计算方法如下：

${\psi }_{\mathrm{xy}}=\frac{{\partial }^2{u}_x}{{\partial x}^2}=u_x(i+1,j)-2u_x(i,j)+u_x(i-1,j)$

${\psi }_{\mathrm{xx}}=\frac{{\partial }^2{u}_x}{{\partial x}^2}=u_x(i+1,j)-{2u}_x(i,j)+u_x(i-1,j)$

${\psi }_{\mathrm{yy}}=\frac{{\partial }^2{u}_y}{{\partial y}^2}=u_y(i,j+1)-2u_y(i,j)+u_y(i,j-1)$

${\psi }_{\mathrm{yx}}=\frac{{\partial }^2{u}_y}{{\partial y}^2}=u_y(i,j+1)-{2u}_y(i,j)+u_y(i,j-1)$

本实施方式中步骤一用粒子图像测速法(Particle Image Velicimetry)测量的t时刻的 超声心动图的血流速度初始速度场如图5所示，

经本实施方式的方法修正后，得到的超声心动图的血流速度速度场如图6所示，比较 图5和图6可知，经过本实施方式的方法修正后，可以清楚看到血流的方向、涡流等局部 特征。

具体实施方式五：本实施方式与具体实施方式四不同的是步骤五中的时间间隔Δt取 1/30000。其它与具体实施方式四相同。

经本实施方式的方法修正后，得到的超声心动图的血流速度速度场如图6所示，比较 图5和图7可知，经过本实施方式的方法修正后，可以清楚看到血流的方向、涡流等局部 特征，图像清晰，细节保留，而且排除了超声图像固有的噪声影响，可以从定性和定量上 更加准确地分析流场的整体趋势和具体参数指标，算法简单。