基于统计模型的心脏力学分析方法

摘要

一种基于统计模型的心脏力学分析方法，通过活动形体模型在不同时刻的建立的病人三维立体心脏图，这些三维立体心脏是由点集组成，这些点分别形成了心脏的心内膜和心外膜，每个点在不同时刻的立体心脏中找到对应的点，按照运动学分析方法计算心脏运动时的受力情况。本发明提供一种计算精度高、运算速度快、可用于临床诊断所需的基于统计模型的心脏力学分析方法。

说明书

技术领域

本发明涉及一种心脏力学分析方法，该方法结合使用图像处理、计算机图形学、医学、数学、统计学、生物学、运动学等技术，主要用于分析心脏运动过程中的动态性能，协助医生对心脏病人进行诊断。

背景技术

心脏是集电生理学、动力学、血液流体力学以及神经、生化控制等于一身的极其复杂的综合系统。建模仿真是研究复杂生物学问题的有效手段。在过去的几年中，人们对心脏结构和功能的生理意义有了深入的了解，并且建立了许多数学模型，努力量化所观察到的心肌机械行为、电传导行为和生物化学行为。但是由于心脏生理病理系统的复杂性，总的来说这些模型是相互独立发展，目前尚没有人能够把心脏的各种机制集成在一起研究。

近些年兴起的虚拟心脏研究将虚拟现实的思想引入到心血管系统这样一个复杂的研究领域，它是利用计算机强大的计算能力和图形处理显示能力，建立虚拟的心脏模型为深入研究心脏活动机理提供了可能。模型不仅要从形态上仿真心脏，而且应能模拟真实心脏的运动过程，仿真心脏的心肌、瓣膜和心腔运动的力学特征、心脏的电传导特性、以及心腔内血液的流体力学特性，并且能够仿真心脏病理状态，为临床诊断疾病提供资料。

目前有一些学者提出了一些基于模型的方法，用于获得心脏的形体和运动的描述。Kyoungju Park，Dimitris Metaxas等学者提出了一种心脏功能分析新的理论。用MRI的图像建立了一个基本的心脏模型，提出了有限元分析的方法计算整体和局部的功能性参数。实验表明，基于这样的模型得出的结构可以表征心脏壁的运动和动态变化规律。而Taratorin和Sideman则把心肌层分割成大量的立方体微元片进行建模和分析，得到效果比较理想。

然而，这些基于模型的心脏运动分析方法由于方法本身的原因，有些在计算精度上还达不到临床诊断所需的要求，有些运算速度慢。

发明内容

为了克服已有的心脏力学分析方法的计算精度或速度达不到临床诊断要求、实用性差的不足，本发明提供一种计算精度高、运算速度快、符合临床诊断所需的要求的基于统计模型的心脏力学分析方法。

本发明解决其技术问题所采用的技术方案是：

一种基于统计模型的心脏力学分析方法，所述分析方法包括：

1)、建立统计模型：对心脏的心室设置所需的物理点，设χ＝{x_i；i＝1....n}是n个被标记过的形体，每一个形体用m个串联的三维标记点向量描述，p_j＝(p_1j，p_2j，p_3j)；j＝1.....m，χ是在3m维空间分布，统计形体采用以下算式表示(1)：

$x=\stackrel{-}{x}+\mathrm{\phi b}---\left(1\right)$

上式中，$\stackrel{-}{x}=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{x}_i$是平均标记点向量，b是形体模型参数向量，φ是协方差矩阵对应的特征向量组成的矩阵，协方差矩阵的算式(2)为：

$s=\frac{1}{n-1}\underset{i=1}{\overset{n}{\Sigma }}(x_i-\stackrel{-}{x}){(x_i-\stackrel{-}{x})}^T---\left(2\right);$

把多个形体按获得时间的顺序进行排列，并计算已排列形体的平均形体，再标准化平均形体，接着用标准化形体重新排列形体，判断是否正交化，如是正交化，形体排列完成，再次迭代，得到最终的心脏的平均形体，即点分布模型，将点分布模型上的统计点进行曲线曲面拟合，得到连续的心脏表面模型；

2)、心脏的运动分析：每个点在不同时刻的立体心脏中找到对应的点，按照以下步骤计算心脏运动的受力大小：

2.1)、通过点分布模型在两个相邻时刻的心脏切片图上建立病人在两个不同时刻的立体心脏图，然后计算模型上每个点在该时间段内的位移，其算式为(3)：

$H=\sqrt{{(x_0-x^{\prime })}^2+{(y_0-y^{\prime })}^2+{(z_0-z^{\prime })}^2}---\left(3\right)$

上式中，x₀、y₀、z₀是前一时刻的点分布模型中某一点P在三维空间中的坐标，x′、y′、z′是后一时刻的点分布模型中P点对应的点P′在三维空间的坐标，H是该点在该时间段运动的位移；

2.2)、根据算式(3)计算出心脏模型上每一点某一时间内的位移，计算在这段时间内的平均速度，其算式为(4)：

v＝H/t    (4)

上式中，t是该时间段的时间长度，v是该点在该时间段内的平均速度；

2.3)、根据相邻时间段计算得到的速度，根据算式(5)计算加速度：

a_i＝(v_i-v_i-1)/Δt；    (5)

上式中，a_i为加速度，v_i为该时间段的速度，v_i-1为相邻前一时间段的速度，Δt为两个时间段的间隔；

2.4)、根据牛顿第二运动定律进行计算，其算是为(6)：

F_t2＝r[S(t₂)-M(t₁)]＝a₁.m

上式中，F_t2为物理点的受力大小，γ是力的弹性系数，a₁是t₁时刻到t₂的平均加速度，m是物理点的质量。

上式中，γ是力的弹性系数，a₁是t₁时刻到t₂的加速度，m是物理点的质量。

作为优选的一种方案：在建立统计模型时，采用SPECT医学图像、核磁共振图像、CT图像、螺旋CT图像、超声图像或者PET图像。

根据具体病人对象的心脏情况，作为候选的方案：所述的心室为左心室或者右心室。

本发明的技术构思为：首先通过活动形体模型在不同时刻的精确的快速的建立病人三维立体心脏图，这些三维立体心脏是由点集组成，这些点分别形成了心脏的心内膜和心外膜，而且对于不同的立体心脏图，点集中点的数目是一样的，并且每个点都可以在不同时刻的立体心脏中找到对应的点。然后按照运动学计算心脏的运动的位移、速度、加速度和受力大小等参数。

本发明的有益效果主要表现在：1、计算精度高、运算速度快、符合临床诊断所需的要求；2、实用性好。

附图说明

图1是SPECT医学图像平面的示意图。

图2是SPECT切片图。

图3是位移H计算的示意图。

图4是位移H大小的数值变化曲线图。

图5是速度计算的示意图。

图6是速度瞬时变化的曲线图。

图7是受力的示意图。

图8是心脏某区域力的动态方向示意图。

图9是区域力的数值标量变化示意图。

具体实施方式

下面结合附图对本发明作进一步描述。

参照图1～图9，一种基于统计模型的心脏力学分析方法，所述分析方法包括以下步骤：

1)、建立统计模型：对心脏的心室设置所需的物理点，设χ＝{x_i；i＝1....n}是n个被标记过的形体，每一个形体用m个串联的三维标记点向量描述，p_j＝(p_1j，p_2j，p_3j)；j＝1.....m，χ是在3m维空间分布，统计形体采用以下算式表示(1)：

$x=\stackrel{-}{x}+\mathrm{\phi b}---\left(1\right)$

上式中，$\stackrel{-}{x}=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{x}_i$是平均标记点向量，b是形体模型参数向量，φ是协方差矩阵对应的特征向量组成的矩阵，协方差矩阵的算式(2)为：

$s=\frac{1}{n-1}\underset{i=1}{\overset{n}{\Sigma }}(x_i-\stackrel{-}{x}){(x_i-\stackrel{-}{x})}^T---\left(2\right);$

把多个形体按获得时间的顺序进行排列，并计算已排列形体的平均形体，再标准化平均形体，接着用标准化形体重新排列形体，判断是否正交化，如是正交化，形体排列完成，再次迭代，得到最终的心脏的平均形体，即点分布模型，将点分布模型上的统计点进行曲线曲面拟合，得到连续的心脏表面模型；

2)、心脏的运动分析：每个点在不同时刻的立体心脏中找到对应的点，按照以下步骤计算心脏运动的受力大小：

2.1)、通过点分布模型在两个相邻时刻的心脏切片图上建立病人在两个不同时刻的立体心脏图，然后计算模型上每个点在该时间段内的位移，其算式为(3)：

$H=\sqrt{{(x_0-x^{\prime })}^2+{(y_0-y^{\prime })}^2+{(z_0-z^{\prime })}^2}---\left(3\right)$

上式中，x₀、y₀、z₀是前一时刻的点分布模型中某一点P在三维空间中的坐标，x′、y′、z′是后一时刻的点分布模型中P点对应的点P′在三维空间的坐标，H是该点在该时间段运动的位移；

2.2)、根据算式(3)计算出心脏模型上每一点某一时间内的位移，计算在这段时间内的平均速度，其算式为(4)：

v＝H/t    (4)

上式中，t是该时间段的时间长度，v是该点在该时间段内的平均速度；

2.3)、根据相邻时间段计算得到的速度，根据算式(5)计算加速度：

a_i＝(v_i-v_i-1)/Δt；    (5)

上式中，a_i为加速度，v_i为该时间段的速度，v_i-1为相邻前一时间段的速度，Δt为两个时间段的间隔；

2.4)、根据牛顿第二运动定律进行计算，其算是为(6)：

F_t2＝r[S(t₂)-M(t₁)]＝a₁.m    (6)

上式中，F_t2为物理点的受力大小，γ是力的弹性系数，a₁是t₁时刻到t₂的平均加速度，m是物理点的质量。

上式中，γ是力的弹性系数，a₁是t₁时刻到t₂的加速度，m是物理点的质量。

在建立统计模型时，采用SPECT医学图像、核磁共振图像、CT图像、螺旋CT图像、超声图像或者PET图像。所述的心室为左心室或者右心室。

本实施例中，考虑左心室的情况，我们总共设置了2848个点，其中心外膜用1777个点表示，心内膜用1071个点表示。

设χ＝{x_i；i＝1....n}是n个被标记过的形体，每一个形体用m个串联的三维标记点向量描述，p_j＝(p_1j，p_2j，p_3j)；j＝1.....m。χ是在3m维空间分布。

目标是获得一个总体上的紧凑的统计形体，这个形体我们可以用下面这样的表达式：

$x=\stackrel{-}{x}+\mathrm{\phi b}---\left(1\right)$

这里的$\stackrel{-}{x}=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}x$是平均标记点向量，b是形体模型参数向量，φ是协方差矩阵：

$s=\frac{1}{n-1}\underset{i=1}{\overset{n}{\Sigma }}(x_i-\stackrel{-}{x}){(x_i-\stackrel{-}{x})}^T---\left(2\right)$

对应的特征向量组成的矩阵，它的几何意义是主成分的集合，但这些主成分必须是正交的，这里的成分是协方差矩阵S对应的特征向量，对应的特征值λ_i(λ_i≥λ_i+1)，如果φ是最大的非零特征值对应的t(t是min{m，n})个特征向量。这样我们就能用等式(1)表示任意一个形体x，这里的φ＝((φ₁|φ₂|…|φ_t)，b是t维向量，$b={\phi }^T(x-\stackrel{-}{x}).$假设我们认为标记点向量服从高维高斯分布，λ_i的变化使得参数b_i的变化。

正确描述形体的先决条件是，这些形体必须在相同的坐标框架下格式化。格式化的目的是消除由于形体在旋转、位移、缩放过程中引起的差异。形体被排列后，残存的差异仅仅是跟形体相关的，因而，可以消除形体投影和缩放等细微变化引起的不良结果。

向量b定义了可变形模型的一组参数．通过改变b我们就能够改变形体，确定x。对向量b中任一个元素b_i取值限制在的范围内，这样就能保证得到的形体和原来的形体相似。

在这种条件下，所有的点可以分布在主轴附近$x\approx \stackrel{-}{x}+\mathrm{\phi b},$这里的b是沿着主轴从到x最近的距离，一般的二维的数据使用单独的参数b即可，在原始的训练集中，那些控制许多模型点的相似模型可能只需要一些参数即可。

在这里，我们描述了一个统计的形体模型，它被用于以图像的方式来表示对象。该对象的形体将用一组N个点来表示，这些点可能有n个自由度。通常这些点有两个或三个自由度。

形体由一些高质量的配置点定义，这些点在形体变化过程中是不变的。在二维或者三维空间中，我们通常要考虑如下变换：平移、旋转和缩放。

我们的目标是通过该模型既可以综合训练集中的相似形体，又可以分析新的形体。首先这些集中的图像都需要进行手工标记标注点，用于表示图像中的对象。

假设现在有一组被排列成同等结构的点xi，这些点形成nd维空问分布。如果能模拟这个分布，我们能产生新的例子，与最初的训练组相似，而且我们检验新的形体是否满足要求。

我们找到一个参数化模型x＝(b)，这里的b是模型参数的一个向量。这个模型能被用来产生新的向量x，如果能模拟参数的分布，通过P＝(b)我们能限制这些向量，使新产生的向量x和训练的组相似。同样的，在用这个模型的时候，应该要去估计p(x)。

从点到单个形体，得到了心脏的平均形体，再次迭代，我们得到了最终的平均形体，这样的形体是具有统计意义的，它上面每一个点能统计得表示对应一般形体的点。由这些点组成了这样的模型的就称为点分布模型PDM，然后，我们将这上面的统计点进行曲面拟合，得到了连续的PDM模型。

在计算运动学的参数之前，我们需要一些标记在心脏内外心膜的点，而这些点是在匹配过程中建立的某个病人的心脏心外膜或心内膜上点。

我们采用的数据源是SPECT图像(如图1)，SPECT心脏图像物理点和标记点很清楚的可以看到，所以大量的物理点能在很短的时间间隔内被标记，而且可以在心脏周期内被全程跟踪，所以用本发明所述的模型技术分析SPECT是非常适合的。我们可以看到在图1中，以黑色的运动初始的交点M为物理点，S点为图像标记点，标记点S在相应的时间时刻对应着相应的物理点M，运动初始t＝1到运动终点t＝5就是一个完整的心脏运动周期(从舒张末期到收缩末期)。

通过分割算法，我们可以得到和建立匹配后的心脏统计模型。

首先我们要对这些相位的位移这个基本的参数进行分析计算。表1和表2列出了心外膜8个点第一相位和第二相位的三维坐标值，它们表示了8个点在空间中的位置和方向。那么，对应的心内膜也存在这样的8个相位的三维标注点序列。

在本实施例中，我们只考虑了左心室心外膜和心内膜总共1164个点，其中心外膜是700个点，心内膜是464个点。这1164个点的运动集合就构成了这整个心肌的运动。所以我们先选出相邻的8个点，这8个点构成了一个心肌面片，分析这面片上的8个点的各项参数。

    X注标值    Y标值    Z标值    171.705    176.967    45.697    171.784    182.399    44.555    168.831    182.377    56.863    170.168    187.293    54.716

    172.108    187.132    46.399    170.240    181.393    51.198    170.070    174.331    53.905    172.012    177.391    43.616

表1心外膜第一个相位8个点的初始位置三维注标；

    X注标值    Y标值    Z标值 172.496994 182.835999 43.898998 168.960007 180.688004 55.717400 168.585999 176.763000 54.618599 170.766006 185.865005 54.752899 172.966995 186.947998 46.797798 170.559006 180.582001 49.895599 169.630005 173.016006 50.903000 173.970993 177.057007 46.435200

表2心外膜第二个相位8个点的初始位置三维注标；

这样的相位一共有7个，相互间隔的时间是100ms，即心脏的严格跳动周期是700ms。显然，各相位间对应点的三维坐标的差就是位移的数值差，相应的向量方向就是位移的方向，即点的运动方向，参照图3、图4。

参照图5、图6，考虑一个心脏周期是由七个相位组成，每个相位持续时间是100ms。一开始我们假设模型数据点和图像数据点的位置是重合的，M(t₁)是物理点，M(t₂)是数据点，M(t₂)是相邻下一个数据点的位置。那么，我们就定义速度的计算方法如下：

v＝H/t    (4)

实际上心肌物理点的速度是变化的，不是匀速的，有加速和减速的过程。加速度是心肌壁运动重要的一个基本参数，加速的快慢是心肌收缩和舒张强度的重要指数。规定心脏模型参数的计算是严格符合牛顿第二运动定律的，其算式为(5)：

a_i＝(v_i-v_i-1)/Δt；    (5)

上式中，a_i为加速度，v_i为该时间段的平均速度，v_i-1为相邻前一时间段的平均速度，Δt为两个时间段的间隔；

根据牛顿第二运动定律，心脏上某点的受力大小为(6)：

F_ti＝r[S(t_i)-M(t_i-1)]＝a_i-1.m    (6)

上式中，γ是力的弹性系数，a₁是t₁时刻到t₂的平均加速度，m是物理点的质量。

参照图7，力F_ti的方向跟运动平面的方向是平行，实际上，力的计算过程可以近似成Spring弹簧模型，但是，我们目前难以确定弹性系数，因此我们采用了经典的牛顿运动定律。

参照图8、9，心脏外心膜的局部区域收缩情况，这八个点中，其中6个点的方向是基本指向内侧的。需要说明的是，小点是运动的初始时刻，大点是终点时刻，这中间的是运动的时间间隔，其值为100ms。