基于MRI图像的心肌三维弹性成像

金杰锋刘华锋*胡红杰

1(浙江大学现代光学仪器国家重点实验室，杭州 310027)2(浙江大学医学院附属邵逸夫医院，杭州 310016)

基于MRI图像的心肌三维弹性成像

金杰锋1刘华锋1*胡红杰2

1(浙江大学现代光学仪器国家重点实验室，杭州 310027)2(浙江大学医学院附属邵逸夫医院，杭州 310016)

材料参数比起运动参数能够更显著地反映物体的内在属性变化，定量获得病变组织的材料参数对心脏疾病诊断具有非常重要的意义。本研究提出一种在利用MRI图像获得部分形变位移等运动信息和约束条件的情况下，依据连续力学模型，使用有限元分析结合H∞滤波估计方法实现三维弹性成像的算法。模拟数据仿真给出了与最小二乘法结果的对比，验证了方法的有效性。最后给出了对病人心脏MRI图像三维弹性成像结果。

左心室;三维弹性成像;有限元方法;H∞滤波

Abstract:Quantification of organic tissue’s material parameters has significant implication in clinic diagnosis of cardiac diseases，since material parameters are far more sensitive measures of the object intrinsic property changes.With known nodal kinematic data and boundary constrains，a new 3D elastic modulus estimation algorithm was based on finite-element method and H∞ filter.Simulation results proved the efficiency of the method compared with the least square method.A myocardiac 3D elastography using cardiac MRI images was presented as well.

Key words:left ventricle;3D elastography;finite elements method;H∞filter

引言

心脏是人体最重要的器官之一，它运作的方式和机理复杂而精巧，对心脏的研究一直是医学上一个非常重要的领域。早期认识心脏主要是根据解剖学的知识，由于技术手段的限制人们无法了解其在人体内部工作的细节。随着现代医学影像技术的发展，人们已经可以无创地观察人体内部组织器官，利用CT、MRI等手段，不仅可以获得二维的人体断层图像，而且可以利用先进的计算机技术生成心脏的三维模型，更加有效地对心脏的运作规律进行研究，从而将对心脏疾病的研究和诊断提升到一个新的高度。

很多研究表明，生物组织的病变会引起弹性系数的变化，尤其是对于像心脏这样由肌肉组织构成的器官。因此在临床上，对心肌组织的材料系数进行量化的测定有着非常重要的实际意义。目前，超声应变成像与弹性成像［1］以及 MR弹性成像(MRE)［2］都已经在心脏运动研究领域得到了应用。但超声应变(应变率)成像只能提供一维信息，而且受到操作手法、成像角度等因素影响。MRE成像需要外部激发装置，目前只能成平面图像，而且现有的图像处理方法也存在分辨率不高或者易受噪声影响等缺点。

连续生物力学模型将通过实验得到的心脏材料属性先验知识引入到心脏图像分析中，取得了比较好的效果，很多研究者在研究中采用了这种方法［3-5］。利用连续生物力学模型研究心脏器官的材料属性的方法是假设心脏的运动信息已知(比如依靠植入种子获得或者通过图像处理的方法测量运动信息，如位移值或瞬时速度值)，然后从这些运动信息中来估计材料参数如杨氏模量、泊松比等，由此可以实现对材料的成像。

以往基于医学图像信息求解组织的弹性模量分布信息的工作，通常只对某一断层图像进行单独处理［6］，忽略了垂直于图像平面方向上的运动信息，且所采用的最小二乘法重建无法应用于更复杂的三维数据，容易受到噪声的影响。随着三维成像和基于图像的三维轮廓运动追踪方法的不断发展［7-8］，已经可以从多层断层图像中提取三维运动信息，利用三维运动信息估计整个心肌的弹性模量，从而更准确地反映实际情况。因此，本研究提出了一种基于有限元与H∞滤波估计的方法，利用MRI心脏图像信息对心肌实现三维弹性成像，并通过实验证明了其有效性。

1 方法

1.1 心肌的有限元表达

使用有限元方法来表达心肌求解域。有限元方法是一种在工程问题中使用非常广泛的近似数值解法，其基本思想将求解域分解为较小单元来处理，在三维问题中，常见的单元类型是四面体单元和六面体单元。由于四面体单元对于复杂形状的网格化有更好的适应性，网格划分算法也比较成熟，而且计算量小，因此本研究采用这种单元形状。对左心室划分网格的方法是，先对MRI等断层成像图像使用图像分割方法获得左心室心肌表层的内外轮廓边界，并在内外轮廓线界定的区域均匀布点，然后使用三维Delaunay四面体剖分算法生成所需的四面体网格。

当四面体单元形变时，四个顶点都有沿 x，y，z三个方向的位移，结合起来用一个矩阵表示:

四面体单元内的位移u可以仅用其四个顶点的位移ue来表示

其中 N为四面体单元形状函数矩阵，具体形式:

其中 Ni，Nj，Nk，Nm分别是四面体单元四个顶点对应的形状函数，可以由下式计算得到:

记由四面体单元四个顶点坐标构成的矩阵为

则ap，n为矩阵 Λ 的第 p行，第 n列元素的代数余子式，单元体积Ve=|Λ|/6。

1.2 心肌三维连续生物力学模型

实际心肌的力学属性非常复杂，从计算的角度考虑，假设心脏是各向同性的线性弹性体，其应力向量σ与应变向量ε成线性关系:

式中，材料矩阵

式中，E是弹性模量，表示物体产生形变的难易度;ν是泊松比，表示横向应变与纵向应变之比。这两项都属于物体的材料性质，实际中泊松比的变化较小，可以认为是常量。

空间中的应变与位移的关系由三维几何方程表示式中，ρ是心肌组织平均密度，I是12×12的单位矩阵。

1.3 从图像获取运动信息

建立序列图像不同时间帧之间像素点的对应关系一直是计算机视觉领域的重要研究课题。在心脏图像研究领域中，对于普通 MRI图像，所能依据的仅仅是心脏内外边界的匹配。一种方法是基于几何形态标记匹配和定位的心肌壁轮廓运动追踪方法［7］。该方法基于左心室轮廓形变在连续的两个时间帧之间尽可能小的假设，利用曲率能量作为匹配原则，来得到节点在前后两个时间帧中的匹配关系。如果利用MRI标记成像和MRI相对比成像技术，那就可以匹配更多像素点，得到更丰富的运动信息，使弹性成像的结果更好。

在获取左心室内外壁三维空间位移场时，本研究利用文献［12］中提出的心脏左心室形状和运动联合估计框架，其主要思想是通过使内能与外能达到平衡获得帧与帧之间的边界点的对应关系，其中内能由采用点的力学特性所决定，外能为边界力、先验力、时间相关力和保形力的组合。

1.4 重建算法

基于有限元方法原理得到心脏运动的系统动

由式(2)和式(8)可以得到:

由此可见四面体单元内各点的应变是一样的，是常应变单元，应变大小由四个节点的位移确定，B是单元应变矩阵。四面体单元刚度矩阵的计算如下式:

四面体单元的质量矩阵采用集中质量矩阵的形式，即认为质量集中在单元的节点上:态方程的矩阵形式如下［9］

式中，R是系统所有节点的载荷向量，反映物体所受作用力的大小和分布;U是系统所有节点的位移向量，其一阶微分和二阶微分分别表示速度向量和加速度向量。M和K分别是系统的质量矩阵和刚度矩阵，均利用相应的单元矩阵(10)和(11)通过有限元方法构造而成。C为阻尼矩阵，假设为瑞利阻尼，满足公式

α，β均取 0.01。

本研究的目标是重建弹性模量，即通过重组系统动态方程(式(12))，利用已有的位移信息和边界力信息来求得弹性模量E的分布。因此我们需要把式(12)的左边变换成 E向量的函数，即 KU=G1E，K=G2E，由方程(7)可知这样的变换是可行的0。Ke可以表示成含Ee项的形式:

由此，根据刚度矩阵构造的方法，可以得到构造G1和 G2的方法:

1)初始化两个N×Ne的空矩阵 G1和 G2，N为系统节点变量数，在三维空间中，每个节点有三个方向的自由度，即系统总节点数×3，Ne为系统总单元数。

3)对于一个编号为j的单元，利用方程(14)构造不包含该单元弹性模量Ej的刚度矩阵。

4)按照系统生成网格时单元和节点的编号规则，把K′j中的元素插入系统临时刚度矩阵中相应的位置中，并将其单元内下标改成相应的系统下标。

6)返回到第2)步直到所有节点遍历。

由此构造出 N×Ne的矩阵 G1，G2。式(12)经过变形:

为了求解E，将式(15)整理成如下

令 G=βG2+G1，R′=R-MÜ- αM，得到简化形式:

这样利用一定的边界条件，系统弹性模量向量E就能直接利用最小二乘法得到近似的解:

在三维空间中，单元的数量常常大于节点自由度。

1.5 H∞滤波估计

由成像手段获得的运动信息必然含有一定噪声，最小二乘法不能很好的消除噪声带来的影响，本研究采用一种基于H∞的噪声扰动全状态信息方法(NPFSI)来实现估计过程。

为了应用H∞滤波估计方法，首先必须建立系统状态方程。设状态向量为x(t)=［U(t)(t)］T，控制向量为w(t)=［0 R］T，将式(11)变形，建立连续时间线性随机系统状态空间表达式:

式中，D为测量矩阵，D的行数为三维运动信息可以测量得到的有限元节点自由度，D的列数为系统总节点自由度。D的每行中，索引位置对应的有限元节点的三维运动信息如果可以测量得到，该索引位置的元素设为1，否则设为0。v(t)代表过程噪声，e(t)代表观测噪声，本研究均假设为0均值的白噪声(E(v(t))=0，E(v(t)v(s)′)=Qv，E(e(t))=0，E(e(t)e(s)’)=Re)。

由于弹性模量E的是包含在参数矩阵中的，为了实现估计的目的，设材料向量为θ=［E］，再将原运动状态量x和材料参数θ组成一个新的状态向量，根据式(15)将原状态方程变形，得到如下一组新的状态方程，这里假设材料参数θ不随时间变化

根据文献［11］由状态方程得到如下一组连续状态估计方程:

为了使估计算法达到预期的效果，在开始估计之前，必须先设定合适的噪声干扰衰减水平，权重矩阵 Q，Q0，Q1和 γ 值。通常很难决定最优的 γ*，它可能是无穷的，意味着无法达到期望的噪声控制水平。然而，实际上如果Q被选定了，γ*能够通过分析获得［12］。

在计算过程中，弹性模量E的初值可以根据对心肌的经验知识，设定与所成像对象的真实弹性模量接近的某一个先验的值，然后利用式(25)～(28)对预设的弹性模量分布E进行递归估计修正，直到获得最优解。

2 实验与讨论

为了验证本方法，首先采用了类似心脏左心室结构的仿真模型进行了实验。仿真模型高100 mm，内径60 mm，外径100 mm，深灰色部分表示正常组织，弹性模量取正常心肌组织的先验值75 kPa，浅灰色代表较硬组织，弹性模量为105 kPa，黑色代表较软组织，弹性模量为45 kPa，其形变状态与位移场分布如图2所示。

将模型的底面固定，顶面施加1 000 N垂直向下的压力，使用 ABAQUS软件生成节点位移数据，作为用于重建的理想数据。

从图2中可以看到不同材料之间的位移差别非常小。这一方面说明从位移数据中重建物体的材料参数非常困难，尤其当数据包含噪声时。另一方面，如果在这样的条件下可靠的重建出了物体的弹性模量分布，那么这说明物体的材料参数能够比运动参数更显著地反映物体属性细微变化。

图3是在不同噪声程度影响下，使用H∞滤波估计方法得到三维可视化结果。

图1 三维仿真模型Fig.1 3D synthetic model

图2 施加力后物体的形变状态和位移场。(a)形变状态;(b)单元节点位移向量Fig.2 Deformed state and displacement field after applying forces.(a)Deformedstate;(b)Nodal displacement field.

图3 不同噪声情况H∞滤波估计三维可视化结果(左列为外侧正面视图，右列为内侧反向视图)。(a)和(b)不含噪声;(c)和(d)SNR=40 dB;(e)和(f)SNR=30 dBFig.3 Visual results of reconstruction from the synthetic data corrupted by noise using H∞filter algorithm.(The left column is front view，the right column is semi-sectional view from opposite direction)(a)and(b)noise free;(c)and(d)SNR=40 dB;(e)and(f)SNR=30 dB

表1是加入不同程度和种类的噪声后重建的结果，可以发现信噪比在30 dB的水平时，最小二乘法已经不能重建出真实分布，但使用H∞滤波估计仍然能得到可以接受的接近真实弹性模量分布结果，而且对不同类型的噪声都有很好的适应性。

真实心脏三维弹性成像实验采用来自医院病人心脏数据，整个心脏分10个断层，每层原始图像大小512像素(512像素，层内分辨率0.704 4 mm/像素，层间隔为8 mm，整个心脏运动周期内包含20帧，从心脏舒张末期开始。

为了获得心肌三维运动信息量，本研究利用文献［12］中提出的心脏左心室形状和运动联合估计框架，其主要思想是利用连续生物力学心脏有限元模型(如图4所示)，每个采样点的受力表示为边界力、先验力、时间相关力和保形力的组合，从而同时求得左心室心肌的运动轨迹和序列图像分割的结果。使用该框架得到的分割结果如图5所示。图6是得到左心室舒张末期到收缩末期的位移值。

表1 H∞滤波重建和LS重建在不同噪声情况下的结果对比 (均值±标准差，kPa)Tab.1 Results of reconstruction using H∞and LS from the synthetic data corrupted by different kinds and degrees of noise(mean±SD，kPa)

图4 心脏左心室有限元模型Fig.4 Left ventricular finite elements model

图5 左心室MRI图像序列的分割结果。(a)舒张末期 (b)收缩末期Fig.5 Segmentation of LV from cardiac MR image sequence(a)end-diastole(ED);(b)end-systole(ES).

图6 左心室从舒张末期到收缩末期的三维位移信息Fig.6 The displacement field of LV model between ED and ES

该框架中加入了几何形态匹配的约束条件，因此对于左心室图像的分割和心肌壁上节点的运动跟踪能够得到比较好的效果，但由于其框架假设心肌材料是均匀一致的，而且所采用的作用力仅是来自图像信息的虚拟力，所以对于壁间节点的运动重建与实际有一定的差异。因此仅采用其框架得到的内外心肌壁上的节点舒张末期至收缩末期的图像帧的三维空间位移值，作为观测信息即式(19)的观测值y。由于没有相应的相对比MRI图像，而心脏舒张末期和收缩末期速度值比较小，y中的速度观测量可以设为0。采用病人血压的收缩压和舒张压之差作为心肌内壁上的压力载荷R。利用本算法进行估计，最终得到弹性成像的结果图7所示。从图6可以看出，上部分心肌在同样压力作用下产生的位移较小，其弹性模量应该比较大，图7的结果正好反映了这一点，说明本算法得到的结果是符合实际意义的。

图7 左心室三维弹性成像结果Fig.7 Visual results of LV elastic modulus distribution

3 结论

本研究提出了利用MRI图像信息，使用有限元结合H∞的滤波估计方法来进行心肌组织的三维弹性成像。仿真实验表明算法在有噪声干扰的情况下，仍可以估计出接近真实的结果，验证了本方法的有效性。进一步的实验利用真实病人心脏MRI图像数据得到了弹性成像结果。未来在对本方法加以一定的完善和实用化之后，可以利用MRI图像数据给医生提供辅助诊断所需的组织材料信息，对于心脏疾病的诊断具有重要参考意义。

［1］Konofagou E，D'hooge J，Ophir J.Cardiac elastography-A feasibility study［J］.2000 IEEE Ultrasonics Symposium Proceedings，2000，1-2:1273-1276.

［2］Elgeti T.Cardiac magnetic resonance elastography initial results［J］.Investigative Radiology，2008;43(11):762-772.

［3］Frangi AF，Niessen WJ，Viergever MA.Three-dimensional modeling for functional analysis of cardiac images:A review［J］.IEEE Transactions on Medical Imaging，2001，20(1):2-25.

［4］Papademetris X，Sinusas AJ，Dione DP，Duncan JS.Estimation of 3D left ventricular deformation from echocardiography［J］.Medical Image Analysis，2001，5(1):17-28.

［5］Shi Pengcheng，Sinusas AJ，Constable RT，et al.Volumetric deformation analysis using mechanics-based data fusion:Applications in cardiac motion recovery［J］.International Journal of Computer Vision，1999，35(1):87-107.

［6］Zhu Yanning，HallTJ，Jiang Jingfeng.A finite-element approach forYoung's modulus reconstruction［J］.IEEE Transactions on Medical Imaging.2003，22(7):890-901.

［7］Shi Pengcheng，Sinusas AJ，Constable RT，et al.Point-tracked quantitative analysis of left ventricular surface motion from 3-D image sequences［J］.IEEE Transactions on Medical Imaging，2000，19(1):36-50.

［8］庄玲，刘华锋，鲍虎军.基于图像的心脏左心室形状和运动参数联合估计框架的研究［J］.自然科学进展，2007，17(3):396-400.

［9］ShiPengcheng，Liu Huafeng.Stochastic finite element framework forsimultaneousestimation ofcardiackinematic functions and material parameters［J］.Medical Image Analysis，2003，7(4):445-464.

［10］Shi Pengcheng，Liu Huafeng.Robust identification of object elasticity［J］.Computer Vision and Mathematical Methods in Medical and Biomedical Image Analysis，2004，3117:423-435.

［11］Didinsky G，Pan ZG，Basar T.Parameter-identification for uncertain plants using H-infinity methods［J］.Automatica，1995，31(9):1227-1250.

［12］Didinsky G，Basar T，Bernhard P.Structural-properties of minimax controllers for a class of differential-games arising in nonlinear H-infinity control［J］.Systems ＆ Control Letters，1993，21(6):433-441.

Myocardiac 3D Elastography Using MRI Images

JIN Jie-Feng1LIU Hua-Feng1＊HU Hong-Jie2
1(State Key Laboratory of Modern Optical Instrumentation，Zhejiang University，Hangzhou 310027，China)2(The Sir Run Run Shaw Hospital，College of Medical Sciences Zhejiang University，Hangzhou 310016，China)

TP391.41

A

0258-8021(2010)04-0564-07

10.3969/j.issn.0258-8021.2010.04.014

2009-11-20，

2010-02-20

国家自然科学基金资助项目(60872068);浙江省自然基金项目(R207119)

*通讯作者。 E-mail:liuhf@zju.edu.cn