面向磁共振弹性成像的数值建模与仿真研究

第６期２０１７年６月电　　子　　学　　报ＡＣＴＡＥＬＥＣＴＲＯＮＩＣＡＳＩＮＩＣＡＶｏｌ．４５　Ｎｏ．６Ｊｕｎ．　２０１７收稿日期：２０１５１００３；修回日２０１６１０１０；责任编辑：梅志强基金项目：国家自然科学基金（Ｎｏ．６１２７１１２３，Ｎｏ．６１５７１１７６，Ｎｏ．６１５１１１４００９９）；安徽省杰出青年科学基金（Ｎｏ．１６０８０８５Ｊ０４）；安徽省国际科技合作计划（Ｎｏ．１５０３０６２０１５）面向磁共振弹性成像的数值建模与仿真研究朱霞丽１，李炳南１，向　馗２（１．合肥工业大学生物医学工程系，安徽合肥２３０００９；２．武汉理工大学自动化学院，湖北武汉４３００７０）　　摘　要：　磁共振弹性成像技术能够定量可视化人体组织的生物力学属性，正逐步成长为一种新型的医学成像手段．然而，其研发不仅要深入探索弹性波的产生、传播及耗散特性，而且需要建立特殊的磁共振成像装备．建立数值模型并开发仿真实验系统，可以有效促进磁共振弹性成像研究工作的开展．本文首先对人体组织进行弹性力学分析，得到外力作用下人体组织的运动微分方程，并利用有限元法求解该方程．同时建立了不同结构和弹性分布的人体组织数值模型．采用剪切弹性波作为成像探针，并根据人体组织的运动微分方程，仿真计算得到弹性波在组织中的传播信息．通过比较不同模型的估算值和理论值，验证了仿真实验系统的有效性和可靠性，说明建立的数值模型可以有效支撑磁共振弹性成像研究．关键词：　弹性波；磁共振弹性成像；软组织；组织弹性；数值仿真中图分类号：　Ｒ８１４　　　文献标识码：　Ａ　　　文章编号：　０３７２２１１２（２０１７）０６１４８３０７电子学报ＵＲＬ：ｈｔｔｐ：／／ｗｗｗ．ｅｊｏｕｒｎａｌ．ｏｒｇ．ｃｎ　ＤＯＩ：１０．３９６９／ｊ．ｉｓｓｎ．０３７２２１１２．２０１７．０６．０２９ＮｕｍｅｒｉｃａｌＭｏｄｅｌｉｎｇｏｆＭａｇｎｅｔｉｃＲｅｓｏｎａｎｃｅＥｌａｓｔｏｇｒａｐｈｙＺＨＵＸｉａｌｉ１，ＬＩＢｉｎｇｎａｎ１，ＸＩＡＮＧＫｕｉ２（１．ＤｅｐａｒｔｍｅｎｔｏｆＢｉｏｍｅｄｉｃａｌＥｎｇｉｎｅｅｒｉｎｇ，ＨｅｆｅｉＵｎｉｖｅｒｓｉｔｙｏｆＴｅｃｈｎｏｌｏｇｙ，Ｈｅｆｅｉ，Ａｎｈｕｉ２３０００９，Ｃｈｉｎａ；２．ＳｃｈｏｏｌｏｆＡｕｔｏｍａｔｉｏｎ，ＷｕｈａｎＵｎｉｖｅｒｓｉｔｙｏｆＴｅｃｈｎｏｌｏｇｙ，Ｗｕｈａｎ，Ｈｕｂｅｉ４３００７０，Ｃｈｉｎａ）Ａｂｓｔｒａｃｔ：　Ｍａｇｎｅｔｉｃｒｅｓｏｎａｎｃｅｅｌａｓｔｏｇｒａｐｈｙ（ＭＲＥ）ｉｓｅｍｅｒｇｉｎｇｔｏｖｉｒｔｕａｌｌｙｐａｌｐａｔｅｈｕｍａｎｂｏｄｙａｎｄｖｉｓｕａｌｉｚｅｔｉｓｓｕｅｅｌａｓｔｉｃｉｔｙ．ＡｌｔｈｏｕｇｈＭＲＥｒｅｃｅｉｖｅｓｍｏｒｅａｎｄｍｏｒｅａｔｔｅｎｔｉｏｎｉｎｃｌｉｎｉｃ，ｔｈｅｓｔｕｄｙｉｓｌｉｍｉｔｅｄｄｕｅｔｏｔｈｅｉｎａｃｃｅｓｓｉｂｉｌｉｔｙａｎｄｅｘｐｅｎｓｉｖｅｎｅｓｓｏｆｍａｇｎｅｔｉｃｒｅｓｏｎａｎｃｅｓｃａｎｎｉｎｇ．Ａｎｕｍｅｒｉｃａｌｍｏｄｅｌｉｓｔｈｕｓｐｒｏｐｏｓｅｄｂｙａｎａｌｙｚｉｎｇｔｉｓｓｕｅｅｌａｓｔｉｃｉｔｙａｎｄｔｈｅｆｏｒｃｅｄｍｏｖｅｍｅｎｔｄｉｆｆｅｒｅｎｔｉａｌｅｑｕａｔｉｏｎ．Ｔｈｅｆｉｎｉｔｅｅｌｅｍｅｎｔｍｅｔｈｏｄ（ＦＥＭ）ｉｓｉｎｔｒｏｄｕｃｅｄｔｏｒｅｓｏｌｖｅｔｈｉｓｎｕｍｅｒｉｃａｌｍｏｄｅｌｏｆＭＲＥ．Ｉｔｓｐｅｒｆｏｒｍａｎｃｅｉｓｖａｌｉｄａｔｅｄｗｉｔｈｄｉｆｆｅｒｅｎｔｍｏｄｅｌｓｏｆｔｉｓｓｕｅｓｔｒｕｃｔｕｒｅａｎｄｅｌａｓｔｉｃｉｔｙｃｏｍｐｏｓｉｔｉｏｎ．ＴｈｅｒｅｓｕｌｔｓｏｆｑｕａｎｔｉｔａｔｉｖｅｅｘｐｅｒｉｍｅｎｔｓｃｏｎｆｉｒｍｔｈａｔｔｈｅｎｅｗｍｏｄｅｌｉｓｅｆｆｅｃｔｉｖｅｔｏｐｒｏｍｏｔｅｎｕｍｅｒｉｃａｌＭＲＥｓｔｕｄｙ．Ｋｅｙｗｏｒｄｓ：　ｅｌａｓｔｉｃｗａｖｅ；ｍａｇｎｅｔｉｃｒｅｓｏｎａｎｃｅｅｌａｓｔｏｇｒａｐｈｙ；ｓｏｆｔｔｉｓｓｕｅ；ｔｉｓｓｕｅｅｌａｓｔｉｃｉｔｙ；ｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎ１　引言　　弹性是人体组织的一个重要生物力学属性：首先，不同组织的弹性不同，人体不同组织之间弹性模量的差异程度明显大于其他物理特性；其次，同一组织在不同生理或者病理状态下，其弹性也会不同，软组织弹性的变化与其病理过程有着密切的关系［１］．临床上已经广泛利用触诊来感知组织弹性的变化，进行健康评估和疾病诊断．但是，触诊面临着不少限制和不便：一是只适用于表面器官，难以获取身体内部组织二维或者三维的弹性分布；二是诊断效果取决于医生的主观感受和个人经验．因此，越来越多的医生希望发展一种能够客观的、无创的人体组织弹性成像技术［２，３］．弹性成像是通过提取组织弹性特征进行成像的一种新技术，主要通过对组织施加外界激励，并应用超声、磁共振等影像技术对组织响应情况进行分析，获得反映组织弹性的图像［４］．磁共振成像（ＭａｇｎｅｔｉｃＲｅｓｏｎａｎｃｅＩｍａｇｉｎｇ，ＭＲＩ）有着良好的软组织分辨率及空间分辨率，是一种重要的影像学手段［５，６］．磁共振弹性成电　　子　　学　　报２０１７年像（ＭａｇｎｅｔｉｃＲｅｓｏｎａｎｃｅＥｌａｓｔｏｇｒａｐｈｙ，ＭＲＥ）能够定量可视化人体组织的弹性属性，是一种动态弹性成像技术［７］．由于磁共振弹性成像为无创伤性检查，具有良好的临床应用前景，已经被应用于多种器官研究，例如脑功能［８］和肝脏纤维化［９］等．人体是一个复杂的不均匀体，使得不能准确估算组织的弹性属性［１０］．但是，要进行磁共振弹性成像扫描，不仅需要对磁共振系统作软硬件改造，技术门槛较高，而且磁共振成像系统一般机时紧张，费用比较昂贵，这都给开展磁共振弹性成像研究造成比较大的不便．因此，有必要在实验前首先开展有针对性的理论研究和仿真分析．实际上，已经有一些采用数值仿真及分析处理平台来推动磁共振弹性成像研究的案例［１１］．在本研究中，将致力于建立一套磁共振弹性成像数值模型，仿真计算弹性波探针在人体组织中的传播模型，分析不同组织对弹性波探针的响应情况，从而为开展磁共振弹性成像系统开发和实验研究打下良好的基础．２　磁共振弹性成像原理　　磁共振弹性成像系统主要包括弹性波产生装置、磁共振扫描系统以及弹性重构算法三个部分，其成像的基本过程是：产生剪切弹性波作用于人体，使体内的组织或器官产生质点位移；通过特殊编码的磁共振脉冲序列获取组织位移相位场；对获得的磁共振相位图进行处理和弹性重构，得到组织器官的弹性属性分布图，供临床医学参考［３］．在进行分析时，通常先假设被测组织是一个均匀的无限介质［３］．当剪切应力以一个固定的频率振动，振动籍由组织内部质点的来回振动得以传播，伴随的波动会传入组织内部．这时，在组织中产生的应变波动有２种：纵波和横波．在同等条件下，纵波的传播速度较快，波长较长，对弹性检测精度较低．因此，磁共振弹性成像主要关注的是剪切横波，称为弹性波［１２］．考虑弹性波的传播速度和弹性常数的关系：υｔ（ｘ，ｙ）＝Ｇ（ｘ，ｙ）槡ρ（１）其中ρ是材料密度，υｔ是弹性波的波速．并且，考虑波速υ能由振动频率ｆ和波长λ表达：υ（ｘ，ｙ）＝ｆ·λ（ｘ，ｙ）（２）将式（２）代入式（１），得到：Ｇ（ｘ，ｙ）＝ρ（ｆｔλｔ（ｘ，ｙ））２（３）式中ｆｔ，λｔ表示弹性波的频率和波长，Ｇ（ｘ，ｙ）表示点（ｘ，ｙ）处的剪切模量．剪切模量和杨氏模量Ｅ存在以下关系，其中ζ表示泊松比：Ｅ（ｘ，ｙ）＝２Ｇ（ｘ，ｙ）（１＋ζ）（４）将式（３）代入式（４），得到求解弹性波波长的公式：λ（ｘ，ｙ）＝１ｆＥ（ｘ，ｙ）２ρ（１＋ζ槡）（５）式（５）说明应变波在较硬的物质中传播时波长较长；在较软的物质中传播时波长较短．在磁共振弹性成像中，驱动装置产生的弹性波从组织器官的表面进入后，在组织内部进行传播，会在垂直于波的传播路径上产生周期性位移，而位移大小与质点的弹性相关．但是，弹性波在组织局部产生的位移幅度比较小（一般在微米级别），故很难用传统的磁共振成像方法来进行检测．所以，必须在成像平面增加双极性对称梯度场，利用运动敏感脉冲序列来检测弹性波引起的微小位移，并编码生成与弹性波振幅正相关的相位信息φ：φ（ｒ，θ）＝２γＮＴ（Ｇ·ξ）πｓｉｎ（ｋ·ｒ＋θ）（６）式中，γ为旋磁比，Ｎ为波动编码梯度的周期总数，Ｔ为外源激励振动的周期，Ｇ为磁场梯度矢量，ξ为振幅正方向矢量，ｋ为弹性波波数，ｒ为质点所处位置矢量，θ为周期性梯度场和振动源之间的相位差．从式（６）可以看出，磁共振信号的相位差跟振幅空间分布有关，呈现周期性的变化．而且，可以增加总的编码梯度周期数Ｎ，从而提高对微小振动检测的灵敏度（可达０１μｍ）［３］．将式（６）变换到频域，得到相位、波动、梯度三者的关系如下：φ（ｆ）＝２γＳ（ｆ）［Γ（ｆ）Ｗ（ｆ）］（７）其中，ｆ表示频率、Ｓ（ｆ）表示傅里叶变换后的时空位移、Г（ｆ）表示傅里叶变换后的波动编码梯度、Ｗ（ｆ）表示傅里叶变换后的弹性波激励函数，而“”表示卷积运算．在磁共振弹性成像中，通过获得相位分布φ（ｆ）、已知的编码函数Г（ｆ）和弹性波激励函数Ｗ（ｆ），根据式（７）可以计算出弹性波在组织中产生的位移函数Ｓ（ｆ）．经过处理和增强后［１３］，可以将Ｓ（ｆ）作为输入，通过波动方程求逆算法，如局部频率估计、偏微分方程或者迭代反演等，求解出弹性波的波长λ（ｘ，ｙ）．然后根据式（３）和（４），计算得出人体组织器官的剪切模量或者杨氏模量．３　弹性力学有限元法　　在研究磁共振弹性成像的过程中，研究弹性波在人体组织中的传播特性是一个非常重要的问题．这里涉及两个方面：一类是正问题，指已知弹性波参数以及组织器官的密度、力学属性等，数值求解弹性波动场的时空变化［１４］；另一类是逆问题，指已知弹性波激励源和弹性波传播的时空分布，估算重构组织器官的弹性属性［１５］．就磁共振弹性成像而言，对获得的波动编码相位信息进行弹性重构，属于逆问题范畴．但是，为了更好地４８４１第　６　期朱霞丽：面向磁共振弹性成像的数值建模与仿真研究求解逆问题，必须优化设计正问题，两者往往相辅相成．３．１　人体组织的弹性力学建模在磁共振弹性成像研究中，一般将人体组织器官假定为异构弹性体，因此适合应用弹性力学的一些基本假设，即人体组织是一个连续的、各项同性的、局部均匀的线弹性体。在此基础上求解磁共振弹性成像的力学问题时，先确定已知条件，例如组织器官的形态分布、力学属性、边界条件等，再求解出应力、形变和位移等．整个过程分两部分建立方程进行求解．３．１．１　运动方程根据微分线段上应变位移的几何条件，建立几何方程，再根据应力应变间的物理条件建立物理方程，即平衡微分方程（应力部分）、几何关系（应变部分）、本构方程（应力与应变的关系），得到平衡状态下组织器官的弹性力学基本方程：（λ＋μ）ｇｒａｄθ＋μ２Ｕ＋ｆ＝０（８）其中，θ为体应变，θ＝ｕｘ＋ｖｙ＋ｗｚ，Ｕ是位移分量ｕ、ｖ、ｗ的函数，ｆ表示外力，其外力分量为ｆｘ、ｆｙ、ｆｚ，由于假设组织器官是一个局部均匀体，可以将λ、μ表达成与位置无关的拉梅常数．在运动状态下（微小扰动），位移Ｕ不仅是点（ｘ，ｙ，ｚ）的函数，还是时间ｔ的函数，而质点的位置也是时间ｔ的函数．在有外源驱动的情况下，需要添加一个微量惯性力∫∫∫ｖ－ρｄ２Ｕｄｔ２，又因为有θｘ＝２ｕｘ２＋２ｖｘｙ＋２ｗｘｚ＝（·Ｕ），可以得到运动状态下用位移表示的各向弹性力学基本方程：（λ＋μ）（·ｕ）＋μ２ｕ＋ｆｘ－ρｄ２ｕｄｔ２＝０（λ＋μ）（·ｖ）＋μ２ｖ＋ｆｙ－ρｄ２ｖｄｔ２＝０（λ＋μ）（·ｗ）＋μ２ｗ＋ｆｚ－ρｄ２ｗｄｔ２＝０（９）其中，ｄｕｄｔ＝ｕｘｄｘｄｔ＋ｕｙｄｙｄｔ＋ｕｚｄｚｄｔ＋ｕｔ，而（ｄｘｄｔ，ｄｙｄｔ，ｄｚｄｔ）为质点的运动速度．在微小形变的情况下，ｕｘｄｘｄｔ、ｕｙｄｙｄｔ、ｕｚｄｚｄｔ都是微量乘积项，可忽略不计．所以，ｄｕｄｔ＝ｕｔ，从而得到ｄ２ｕｄｔ２＝２ｕｔ２，ｄ２ｖｄｔ２＝２ｖ