CN101515266A - 一种利用计算机对心脏运动序列图像的处理方法 - Google Patents

Abstract

本发明公开了一种根据心脏运动序列图像进行心脏运动和材料特性联合估计的处理方法。通过获得的一个心动周期的心脏运动序列图像，得到采样节点的位移和速度。采用图像分割方法分割出序列图像上的心脏边界，匹配连续两帧之间的边界得到边界位移量。根据心脏边界将心脏图像生成有限元模型，输入生物力学模型参数，通过扩展卡尔曼滤波方法对心脏的运动和材料特性进行迭代估计。本发明基于状态空间的扩展卡尔曼滤波理论，结合有限元模型，实现在已知部分节点位移和边界条件的情况下，方便准确地得到心脏运动参数和心肌组织材料特性。

Description

一种利用计算机对心脏运动序列图像的处理方法

技术领域

本发明涉及一种利用计算机对心脏运动序列图像的处理方法，具体地说是一种根据心脏运动序列图像进行心脏运动和心肌材料特性联合估计的处理方法。

背景技术

人体的血液循环系统负责向全身各个器官提供氧气和营养物质，并带走代谢产生的废物。心脏是循环系统的动力器官，由于心脏的“泵”的作用，血液循环才得以维持。心脏的泵血能力使得循环系统保持着足够的压力。如果心脏不能实现泵血功能，动脉血压即迅速下降，使全身各器官的供血不足，从而发生功能障碍以至危及生命。在心脏疾病的诊断和治疗中，如何从定量的角度对心脏功能和结构状况进行分析，是国内外医学图像分析领域和生物力学领域一直关注的重要课题。理论研究和大量临床实验表明，心室的形状变化与某些心脏疾病有着直接联系，通过分析心肌的运动状况，可以得到大量关于心脏结构和功能的重要信息，如应力应变状态。心肌的应力应变状态具有重要的生理、病理启示：它影响冠脉血流的分布和耗氧状况；也可作为心肌变形、肥大和纤维化、心肌梗塞的诊断和修复的信号等等。

现代医学成像技术的发展，使得对心脏形态变化、运动情况等的非入侵观察和分析成为可能。研究人员和医生们可以利用MRI、CT、超声波以及PET等多种成像技术来取得心脏的结构或功能图。

在心脏运动分析领域，通常的研究方法是利用MRI或者其他成像工具采集一个心脏循环的多幅图象；然后进行图像分割，接着得到幅-幅之间的直接匹配位移数据，即特征点的位移场；最后利用这些稀疏的部分数据来重建整体的光滑位移场。

其实这类问题都可归结为已知部分数据来估计整体的逆问题，用微分方程可以表示为S(q，u(q))＝F，q为材料模型参数，u为系统的状态比如位移量，S为系统微分操作算子，F为负载，目的是从部分u值中估计出整体的u值。为了使逆问题得到唯一解，需要加上约束模型即正则化处理。由于心脏信号的往复运动特性，人们也意识到进行多幅(multi-frame)分析的重要性。Meyer等人假设心脏组织元素的运动轨迹为椭圆，使用Kalman滤波器来估计左心室的形变；McEachen等人使用自适应滤波器；Prince等人使用Kriging方法来跟踪心脏的运动。但是，无论是frame-to-frame的分析还是multi-frame的分析，无论使用那种模型作为约束条件，这些模型在求解过程中都是固定不变的(即q值是固定不变的)。而研究人员早已经意识到这些模型的参数选择将对分析的最终结果产生重要影响。

连续生物力学模型因其可以提供丰富的力学参数，在心脏分割和心脏运动分析领域都得到了广泛的应用。连续生物力学是应用力学原理和方法对生物体中的力学问题进行定量研究的生物物理学分支，其重点是研究与生理学、医学有关的力学问题。连续生物力学的基础是能量守恒、动量定律，质量守恒三定律。利用生物力学模型可以将复杂的生物组织属性转换为数学上可以求解的形式。

生物力学中研究心脏的力学和材料属性是假设心脏的运动信息已知(比如依靠植入种子获得或者通过成像方法测量运动信息，即u值为完全已知)，然后从这些已知的信息中来估计本构关系中的材料参数如杨氏模量、泊松比等等。标识MRI技术最近被广泛用于心脏的研究，心脏材料参数数据通常是通过极小化标识MRI的观测值与有限元模型的预测值来获得。但是，问题是如何从标识MRI或者其他成像手段中获得可靠的位移运动信息依然是个待解的问题。

总结起来，心脏运动是非刚性的运动过程，从心脏力学(或者数学)模型出发，给定材料(模型)参数和边界约束条件，求解心肌的位移、应变应力场，是心脏运动分析的研究内容。而把该问题进行逆向分析，求解心脏的材料特性，是心脏材料分析的研究内容。尽管现有的分析方法都取得了一定的成效，但是它们都是割裂地看待这两个密切相关的问题，一方面，对于某一具体的研究对象(比如特定的病人)，我们只知道心脏材料模型的先验数据但并不知道符合该特定对象的约束模型及其参数，因而所求得的运动分析结果往往并不能完全反应该特定对象的心脏运动状态，如何选定合适的模型和参数成为问题的瓶颈。另一方面，从图像中观测到的运动数据通常是稀疏和存在噪声的，因而得到的心脏材料特性往往不够理想。

发明内容

本发明提供一种利用计算机，根据心脏运动序列图像，对心脏运动和心肌材料特性联合估计的处理方法。

本发明在充分分析心脏运动图像的成像特点和退化问题的基础上，针对不同模态下图像的特征，提出相应的解决方案和初步建立适应于医学图像的分割方法与理论，在此基础上完成对心脏图像的分割，然后从序列图像信息中同时估计心脏的动力学参数和材料属性参数。即我们的最终目标就是利用含有噪声的不完全数据，确定最佳的材料模型参数，同时在物理空间内求解最佳的状态参数。

本发明的处理方法包括以下步骤：

(1)输入包含至少一个心动周期的心脏运动序列图像；

(2)在心脏运动序列图像上设置采样节点，测量采样节点的位移和速度；

(3)采用图像分割方法分割出序列图像上的心脏边界，匹配连续两帧之间的边界得到边界位移量；

(4)根据心脏内外边界第一帧的图像，在内外两条边界之间使用区域填充算法均匀布点，得到心脏的有限元模型；

(5)输入生物力学模型参数；

(6)通过扩展卡尔曼滤波方法对心脏的运动和材料特性进行迭代估计。

所述的扩展卡尔曼滤波器，是通过将期望和方差线性化来对心脏的运动过程进行估计，同时在其线性化的过程中也估计出了作为卡尔曼方程参数一部分的心肌材料参数，从而实现了两者的联合估计。

上述步骤(1)中的心脏运动图像可以是MR tagging图像、MR相对比图像、普通MR图像、CT序列图像或者PET序列图像。

上述的有限元模型通过Delaunay剖分或者其他有限元网格生成。

本发明还可通过计算机可读介质，存储上述技术方案的计算机可执行程序代码。

本发明的优点是：基于状态空间的扩展卡尔曼滤波理论，结合有限元模型，实现在已知部分节点位移和边界条件的情况下，方便准确地得到心脏运动参数和心肌组织材料特性。

附图说明

图1是本发明的流程框图；

图2a、2b心脏图像示意图；图2c为现代医学的染色技术标记出来的受伤位置；图2d为边界点的心脏周期的位移轨迹；

图3a、3b、3c为心脏左心室的有限元网格生成结果示意图；

图4a、4b、4c、4d、4e、4f为采用状态空间方法与生物力学模型相结合，估计出的各帧之间的位移场的示意图；

图5a、5b为利用本发明方法估计得到的心脏左心室杨氏模量和泊松比分布图。

具体实施方式

本发明的处理方法包括以下步骤：

(1)输入包含至少一个心动周期的心脏运动序列图像；

(2)在心脏运动序列图像上设置采样节点，测量采样节点的位移和速度；

(3)采用图像分割方法分割出序列图像上的心脏边界，匹配连续两帧之间的边界得到边界位移量；

(4)根据心脏内外边界第一帧的图像，在内外两条边界之间使用区域填充算法均匀布点，得到心脏的有限元模型；

(5)输入生物力学模型参数；

(6)通过扩展卡尔曼滤波方法对心脏的运动和材料特性进行迭代估计。

根据心脏力学模型，心肌的动态方程为：

$M\stackrel{\·\·}{U}+C\stackrel{\·}{U}+\mathrm{KU}=R---\left(1\right)$

其中M，C和K分别是质量，阻尼和刚度矩阵，R为载荷向量，U是位移向量。M是材料密度的函数，不随时间和空间变化。K是材料结构规律(应变-应力)的函数，与材料特定的杨氏模量和泊松比有关，而这些材料参数会随时间和空间变化。在本发明中，这两个量被当作仅知道部分先验统计值的随机变量，需要根据运动函数来估计。C与频率有关，这里假设为瑞利阻尼，满足公式C＝αM+βK。

将公式(2)变形，可以得到连续时间线性随机系统的状态空间表达式：

$\stackrel{\·}{x}\left(t\right)=A_c\left(\mathrm{\θ}\right)x\left(t\right)+B_{c}w\left(t\right)---\left(2\right)$

y(t)＝Dx(t)+e(t)                                (3)

其中：材料向量θ＝[E  v]^T，状态向量 $x\left(t\right)={\left[\begin{array}{cc}U\left(t\right)& \stackrel{\·}{U}\left(t\right)\end{array}\right]}^T,$ 控制量w(t)＝[0   R]^T。

系统矩阵是材料向量，即杨氏模量E和泊松比v的函数：

$A_c=\left[\begin{array}{cc}0& I\\ -M^{-1}K& -M^{-1}C\end{array}\right],$ $B_c=\left[\begin{array}{cc}0& 0\\ 0& M^{-1}\end{array}\right]$

测量矩阵D为单位矩阵，e(t)代表观测噪声，是加性，零均值的白噪声。

方程(2)和(3)描述了一个连续时间系统的离散时间测量，或者说采样系统。方程的输入由系统方程计算而来，是在采样间隔T内分段常数。由此得到系统离散时间方程：

x(t+1)＝A(θ)x(t)+B(θ)w(t)+v(t)                            (4)

其中： $A=e^{A_{c}T},$ $B=A_c^{-1}(e^{A_{c}T}-I)B_c.$

为了利用扩展卡尔曼滤波方法对材料参数进行估计，这里将状态向量x和材料参数向量θ组成一个新的状态向量z＝[x θ]^T，由此根据方程(3)和(4)推出一组新的状态与测量方程：

$z(t+1)=\left[\begin{array}{c}A\left(\mathrm{\θ}\right)x\left(t\right)+B\left(\mathrm{\θ}\right)w\left(t\right)\\ \mathrm{\θ}\end{array}\right]+\left[\begin{array}{c}v\left(t\right)\\ 0\end{array}\right]=f(z\left(t\right),w\left(t\right))+\left[\begin{array}{c}v\left(t\right)\\ 0\end{array}\right]---\left(5\right)$

$y\left(t\right)=\left[\begin{array}{cc}D& 0\end{array}\right]\left[\begin{array}{c}x\left(t\right)\\ \mathrm{\θ}\left(t\right)\end{array}\right]+\left[\begin{array}{c}e\left(t\right)\\ 0\end{array}\right]=h\left(z\left(t\right)\right)+\left[\begin{array}{c}e\left(t\right)\\ 0\end{array}\right]---\left(6\right)$

运动和材料参数联合估计问题可以被理解为一个非线性系统的状态估计问题。方程(5)和(6)可以引出一个基于连续系统离散测量的扩展卡尔曼滤波(EKF)框架的滤波问题的解。扩展卡尔曼滤波(EKF)框架是将每一时刻的状态方程线性化而得到。由此得到对运动状态和材料参数联合估计的一组递归公式如下：

状态向量更新方程：

$\hat{x}(t+1)=A\left(\widehat{\mathrm{\θ}}\left(t\right)\right)\hat{x}\left(t\right)+B\left(\widehat{\mathrm{\θ}}\left(t\right)\right)w\left(t\right)+L\left(t\right)[y\left(t\right)-D\hat{x}\left(t\right)]---\left(7\right)$

$\widehat{\mathrm{\θ}}(t+1)=\widehat{\mathrm{\θ}}\left(t\right)+G\left(t\right)[y\left(t\right)-D\hat{x}\left(t\right)]---\left(8\right)$

状态向量初值： $\hat{x}\left(0\right)={\hat{x}}_0,\widehat{\mathrm{\θ}}\left(0\right)={\widehat{\mathrm{\θ}}}_0.$

方程(7)和(8)中卡尔曼增益系数L(t)，G(t)的更新方程：

$L\left(t\right)=[A\left(\widehat{\mathrm{\θ}}\left(t\right)\right)P_1\left(t\right)D^T+M_t{P}_2^T\left(t\right)D^T]S^{-1}\left(t\right)---\left(9\right)$

$G\left(t\right)=P_2^T\left(t\right)D^T{S}^{-1}\left(t\right)---\left(10\right)$

状态向量误差协方差矩阵P时间更新方程：

$P\left(t\right)=\left[\begin{array}{cc}{P}_1\left(t\right)& P_2\left(t\right)\\ P_2^T\left(t\right)& P_3\left(t\right)\end{array}\right]$

$P_1(t+1)=A\left(\widehat{\mathrm{\θ}}\left(t\right)\right)[P_1\left(t\right)A^T\left(\widehat{\mathrm{\θ}}\left(t\right)\right)+P_2\left(t\right)M_t^T]$

$+M_t[P_2^T\left(t\right)A^T\left(\widehat{\mathrm{\θ}}\left(t\right)\right)+P_3\left(t\right)M_t^T]+Q_v-L\left(t\right)S\left(t\right)G^T\left(t\right)---\left(11\right)$

$P_2(t+1)=A\left(\widehat{\mathrm{\θ}}\left(t\right)\right)P_2\left(t\right)+M_t{P}_3\left(t\right)-L\left(t\right)S\left(t\right)G^T\left(t\right)---\left(12\right)$

P₃(t+1)＝P₃(t)-G(t)S(t)G^T(t)                            (13)

其中：

S(t)＝DP₁(t)D^T+R_e                                            (14)

$M_t=\frac{\∂}{\∂\mathrm{\θ}}(A\left(\mathrm{\θ}\right)\hat{x}+\mathrm{Bw}\left(t\right)){|}_{\mathrm{\θ}=\widehat{\mathrm{\θ}}}---\left(15\right)$

R_e是观测误差协方差矩阵。

将图像数据和生物力学模型设定的初始条件和边界约束作为初始值代入方程，迭代运算直到得到收敛的最优结果。

设迭代过程循环次数为j，帧数为i。j＝1，i＝1时，初始位移为0，否则初始位移是直到第j次循环，第i-1帧之前所有帧计算得到的预测值。

如果有其他测量信息，比如MR相对比速度图像，那就用图像中的信息来提供心肌壁间点x和y方向的瞬时速度。对于其它点，我们使用直到第j次循环，第i-1帧之前所有帧计算得到的预测值；

初始所有点的加速度使用直到第j次循环，第i-1帧之前所有帧计算得到的预测值；

j＝1，i＝1时，初始的杨氏模量和泊松比根据先验值来设定(对于人体来讲，可以分别采用E＝75000Pacscal，γ＝0.47)。否则使用直到第j次循环，第i-1帧之前所有帧计算得到的预测值；

初始的等效总载荷通过系统方程根据其它初始条件计算得到。

将初始条件代入方程(7)-(15)，进行迭代运算，直到得到收敛的最优结果，完成对心脏运动和材料参数的联合估计。

图2a为相位对比MR心脏强度图像；图2b为相位对比MR心脏速度矢量图；图2c为现代医学的染色技术标记出来的受伤位置(绿色线包围的地区)；图2d为边界点的心脏周期的位移轨迹。

图3a、3b、3c为心脏左心室的有限元网格生成结果示意图，从左到右分别是初始网格(舒张末期)，变形后网格(收缩末期)和从舒张末期到收缩末期之间的位移向量图。

图4a、4b、4c、4d、4e、4f为采用状态空间方法与生物力学模型相结合，估计出的各帧之间的位移场的示意图。请注意受伤区域的明显的运动滞后性；

图5a、5b为利用本发明方法估计得到的心脏左心室杨氏模量和泊松比分布图。

Claims (4)

1、一种利用计算机对心脏运动序列图像的处理方法，其特征在于：包括以下步骤：

(1)输入包含至少一个心动周期的心脏运动序列图像；

(2)在心脏运动序列图像上设置采样节点，测量采样节点的位移和速度；

(3)采用图像分割方法分割出序列图像上的心脏边界，匹配连续两帧之间的边界得到边界位移量；

(4)根据心脏边界将心脏图像生成有限元模型；

(5)输入生物力学模型参数；

(6)通过扩展卡尔曼滤波方法对心脏的运动和材料特性进行迭代估计。

2、根据权利要求1所述的方法，其特征在于：步骤(1)中心脏运动图像是MR tagging图像、MR相对比图像、普通MR图像、CT序列图像或者PET序列图像。

3、根据权利要求1所述的方法，其特征在于：所述的有限元模型通过Delaunay剖分生成。

4、一种计算机可读介质，包括存储有执行权利要求1所述步骤的计算机可执行程序代码。

Images