有限元法心脏运动建模

第四节　有限元法心脏运动建模

基于有限元方法的模型更侧重于分析心室各个内部单元的应力应变分布情况，在局部运动分析方面有限元模型比可形变模型更细致精确。之所以采用四面体单元划分心脏表面实体模型，是因为用四面体单元划分三维实体很灵活，可以逼近较复杂的几何形状，并按需要在局部区域随意加密或变稀网格，而心脏的几何形状比较复杂、不规则，适合用四面体单元划分。

一、心脏有限元模型的应变分析

心肌是近似不可压缩的，心肌被视为是一种非刚性的会发生形变的弹性体材料，心脏在收缩、舒张时必然会产生体力或表面力；而应变分析是描述连续体内部形变的方法，是目前用于量化分析心脏形变的常用工具。许多心脏运动分析算法都使用应变张量来表示心肌的形变，应变张量是弹性力学中的概念，能包括形变的全部信息。该方法能对心肌组织的形变进行精确估计，因此被广泛应用于描述弹性塑性体的形变。在有限元模型中，重建心脏在整个心动周期中的3D形状以及任意心肌点的运动轨迹，并能对心肌的运动情况进行分析，是目前心脏图像分析领域中重要的研究方向。

图102为初始时刻的心脏表面有限元模型，全部由四面体单元组成，该模型由619个四面体单元组成，共771个节点。运用四面体单元剖分三维实体时，可以更灵活地进行剖分，逼近复杂的几何形状，也可以按需要在局部区域随意加密或变稀网格。通常在应力集中区域或应力梯度高的区域布置较密的网格，在应力变化平缓的区域布置较稀疏的网格。这样做可以同时满足精度和效率两方面的要求。

根据所建立的基于3D表面重建的动态有限元方程的模型以及基于形变模型的约束条件，就可得到整个结构的所有节点的位置矢量，由形态函数可得到表面上任一点的位置。考虑到心肌运动应该是连续的弹性运动，本书运用薄板样条对模型受力进行了优化，针对本书建立的四面体元心脏表面三维有限元模型，关键是导入边界条件求得作用在各单元节点上的节点力，计算模型中各单元节点的位移。已知位移，计算单元节点的应变，其中某一点的应变状态可用六个独立的应变分量表示：

其中：ε_x，ε_y和ε_z是正应变，γ_xy，γ_yz和γ_zx是切应变。用u，v，w表示四面体各节点沿x，y，z方向的位移，这些位移一般应为各节点坐标的函数，即

为了反映某节点总的应变状态，把公式（1016）中6个应变分量用米泽斯屈服条件转换为米泽斯应变，即：

所以，如何求左心室表面各单元节点的位移是求解心脏表面应变分布的关键。

心脏表面的几何形状比较复杂、不规则，适于用四面体单元划分，在所建立的有限元模型中，目前非刚体运动对应点估计方法已成功地运用于心脏心室的表面运动跟踪以及大脑皮质表面点匹配。

UNGCK算法是一种结合了单位法线（unit normal，UN）算法和高斯曲率（Gaussian curvature，GC）算法各自优点的非刚体运动对应点估计方法，而且已应用在测量犬心脏的左心室体积膨胀系数上。有限元方法和UNGCK算法结合分析心脏表面应变，其适用于心脏三维实体的有限元单元网格模型。

在曲面造型中，曲面在任一点附近的形状与在该点曲面的主曲率的乘积即高斯曲率有关。该点与附近点的高斯曲率比较可以反映出该点附近的形状变化，故可以用高斯曲率表达该点的形状信息，对该点附近的形状质量进行评判。高斯曲率计算原理图图10-3所示。三维空间的3点确定1个平面，在这一平面内的这3点确定1个球，即确定该球的半径和沿着球的切矢。设表示曲面S在一点P上的单位法矢，切S且经过的平面与曲面相交成一条曲线。同样，不经过但经过P点的平面与曲面同样也可以相交成一条曲线。让每一个法平面与一个方向及单位切矢对应，即在曲线P点，一个法曲面曲率k_n，对应一个位置。这个法曲面曲率随着切平面绕的旋转而变化。k_n存在最大和最小值，即为P点的主曲率。

图10-3　高斯曲率计算原理图

UNGCK算法具有较高的准确率和较好的数值稳定性，但要直接用在心脏有限元模型上还应做相应的改进。我国学者刘复昌、夏德深在相关的论文中提出了一种基于左心室有限元模型的α剪枝平均UNGCK算法。论文的α剪枝平均UNGCK算法考虑到左心室有限元模型表面由许多三角形拼成，拓扑关系不规则，因此取在表面上与节点邻近的一些点代替节点邻域拟合曲面。在网格节点所在的拟合曲面上，以该网格节点为中心的邻域内重新找运动误差最小的点作为对应点，这样找到的对应点不再局限于网格节点。采用经典的刚体配准算法——最近点迭代（iterative closest point，ICP）算法来估计相邻时刻对应点的候选点，在一对多的对应关系下，它选取距离最近的且该点法线方向与待匹配点法线方向夹角最小的点作为对应点。α剪枝平均的好处就是可以去除数据集中的一些不正常点（即过大或过小点），而对正常的点改变很小，这样可避免把不正常点以相同的权值和正常点作平均，比直接对数据点作平均要合理得多。他们提出基于有限元和α剪枝平均UNGCK算法的左心室应变分析，利用α剪枝平均对UNGCK算法存在的局部运动函数与物体整体运动不一致的情况作修正，能直观、有效地反映左心室在整个心脏收缩期的形变情况及应变分布。

二、有限元模型的心脏运动建模

如何求心脏表面各单元节点的位移是求解心脏表面应变分布的关键。显然四面体元有12个自由度，每个节点有3个自由度，因此，对一个有n个节点的结构而言，其整体刚度矩阵K是3n×3n的，可求解δ的结构节点位移矢量，F的结构节点载荷矢量。一旦求得未知的位移与节点力矢量，就可以用公式（1019）求得每个单元的应力矢量。

公式（1019）中，σ是6×1的单元应力矢量，δ是12×1的单元节点位移矢量，每个单元的σ形式为：

其中前三项为正应力，后三项为切应力，在实际的程序设计实验中，这些参数的选择将根据要解决的问题进行调整和取值。

在建立基于有限元方法的心脏实体模型中，初始时刻的心脏表面有限元模型全部由四面体单元组成，该模型由619个四面体单元组成，共771个节点。在实际进行的程序检测中，我们发现771个节点中的12个点是在心脏表面有限元模型的上边缘，形成锯齿，它们即所谓的孤立点，有的形成的四面体体积V为零的，有的没有成为在我们所生成的四面体中的单元节点。因此必须把这些点先处理好，使其不得参与后面的有限单元的计算，以免影响求位移与应力的正确性。通过应用程序我们检测到如果有四面体单元的体积为零的四面体，必须排除出去，否则会出现分母为零时得到无穷大的结果，显然这是不合理的。经调试程序时发现生成的四面体单元中共有4个四面体单元的体积为零。我们必须要把这四个四面体体积为零的四面体作上标记，使其不参加后面的有限元方法的计算，以免影响整个程序计算的精度与速度。综合上述考虑，那么对于心肌这种被认为是近似不可压缩的物质，我们参考人体血流动力学参数取其泊松比v都取0.4，弹性模量即杨氏模量E取50kPa，质量密度ρ取1.07g/cm³。我们在有限元模型的心脏表面位移计算算法研究中，实际程序中采用的点是771个点减去这四个四面体体积为零的12个点后的759个节点，619个常应变四面体单元来计算的。

求节点位移矢量的算法思路设计

在进行有限元法求表面节点的位移矢量时，其中我们就要利用虚功原理与变分原理以及至关重要的边界条件来进行各初始条件的导入，进而求得在载荷的边界条件下心脏表面节点的位移矢量。在相关文献中的论述中计算单位法线及高斯曲率时，采用等间隔分布的网格数据点（如与坐标轴对应的点）作为节点邻域拟合曲面，也就是先将每层短轴方向的轮廓点拟合，得到短轴平面上形式统一的轮廓点数据；然后将长轴轮廓点数据转换成短轴平面数据形式，得到形式完全一致的心脏内轮廓表面线模型，考虑了肌纤维方向及其横截面方向的差异，运用了运动心脏模拟器验证本书模型的结果，简言之就是等参元、等间隔分层长、短轴法。由于先前采用的四面体算法不是按分层等间隔来进行考虑的，因此上述方法并不太适用于我们所建的心脏表面有限元模型，那么我们采用的是近似椭球体的各参数来进行模拟。其思想是：α表示沿模型的基部到顶部方向的取值为-π/2≤α≤π/2，β表示沿模型圆周方向的取值，-π≤β≤π，其边界载荷力F是基于系数函数cosβ，cosα，sinα的形式，向外、向内的F分别表示舒张、收缩的载荷力均由椭球体的参数函数来确定。考虑到进行的是心脏表面的四面体模型，并由于心脏本身的功能特点及人体血流动力学参数，心脏的收缩与舒张主要发生在下半部，即左、右心室，故应着重使得分析中的中心坐标点下移到心脏下半部的左、右心室的中心，使其满足并接近真实的心脏收缩与舒张形态。