摘要：弹性成像在医学成像领域具有广阔的应用前景.在本文中，有限元法被用于模拟角膜组织的超声弹性成像，通过改变角膜病变组织与正常组织杨氏模量的数值大小及分布，模拟角膜病变的不同情形，分析计算各种不同情形时的应变、应力和位移分布，并对仿真结果进行分析.当角膜正常组织与病变组织的杨氏模量之比为1:4时，病变区域中心的应变为-0.00854，而对应正常角膜相同位置的应变为-0.02277，为病变区域中心的应变值的2.67倍.病变区域中心的应力为0.04337，而对应正常角膜相同位置的应力为0.02729，相当于病变区域中心的应力值的0.64倍.当角膜正常组织与病变组织的杨氏模量之比为4:1时，情况刚好相反.基于角膜组织的粘弹性，优化了角膜组织模型.结果表明，利用应变压缩方法可以实现离体组织的弹性成像，超声弹性成像中病变组织与正常组织的生物力学响应差异明显，研究结果可为弹性成像在眼科临床应用上提供理论指导.
关键词：生物力学;超声弹性成像;有限元全眼模型;粘弹性;角膜
在传统意义上，超声成像使用声阻抗作为其成像参数，能够反映组织的声学特性的改变，但是，在某些组织产生病理变化以后，声阻抗改变很小，但弹性变化很大，这需要超声成像可以通过成像显示出力学性质的变化[1].自从Ophir于1991年提出了超声弹性成像技术的概念[2]以来，它就一直受到医学界关注并迅速发展，从那时起，弹性成像成为了超声成像医学应用领域的热门研究课题.
弹性成像基本的原理是在生物组织上提供内部(包括其自身)或者外部、动态或者静态(准静态)激励，在物理定律作用下，组织将产生位移、应变和速度等响应，并且它们的分布会随着激励而改变[3-4].这些变化可以通过一系列数字信号处理方法转化为图像，从而获得组织内部力学性质的弹性分布.超声弹性成像主要包括传统的弹性成像，剪切波成像和声辐射力成像[5].传统的弹性成像是一种静态压缩方法，可以获得组织的整体弹性分布.已有众多学者对超声弹性成像进行理论及临床研究.先前的研究主要集中在原理验证、信号处理、噪声抑制、块匹配和相移的混合算法、应变估计算法以及二维弹性成像算法追踪盆底肌位移场等方面[6-10]，并且超声弹性成像逆问题的研究也得到了发展[11].而梁铭等[12]利用超声弹性成像技术对乳腺病灶切除的190例患者共212个病灶进行检测，结果表明超声弹性成像技术对早期乳腺癌有较好的诊断价值.Cassinotto等[13]对349例需要进行穿刺的慢性肝病患者利用超声弹性成像技术进行肝硬度测量，结果显示超声弹性成像技术更有助于精确评估慢性肝病肝纤维化.张明博等[14]通过对34例前列腺增生患者进行超声弹性成像后，结合相关指标分析得出，前列腺组织弹性模量与前列腺增生导致的膀胱出口梗阻(bladderoutletobstruction,BOO)程度具有良好的相关性，结果表明超声弹性成像可以为前列腺增生患者提供无创评价BOO程度的新方法.局部弹性信息可以通过声辐射力成像和剪切波成像获得.为了能够获得弹性特性的整体分布，对不同的区域进行聚焦是必要的.目前，还没有研究出低噪声、实时、快速、高分辨率的成像方法.
有限元法是研究人体组织机理的重要方法[15]，它可以模拟生物软组织的应变分布，并研究能对弹性成像产生影响的因素，如软组织的不同力学和几何参数等[16].离体组织的弹性成像是研究弹性成像的基础.尽管前人在超声弹性成像技术上做了很多模拟和实验工作[16]，但是大多数模拟是基于二维模型的，并且没有眼组织的相关模拟研究.本文主要基于人眼三维有限元模型，考虑了粘弹性，对角膜疾病进行超声弹性成像的仿真研究.
1研究方法
本研究通过建立人眼有限元全眼模型，并利用COMSOLMultiphysics软件对角膜组织进行有限元分析.通过提取角膜后表面Y方向的应变、应力和位移数据，对比分析了正常和病变角膜应变、应力和位移的差异性，进而确定角膜病变与超声弹性成像结果的关联性，实现超声弹性成像在角膜疾病诊断中的仿真研究.
1.1几何模型
人眼模型的精确构建是剖析问题的关键.在本文中，使用了经典Gullstrand精确模型眼的数据[17]，角膜位于眼球的前部，约占眼睛外纤维膜的1/6.角膜前表面的曲率半径为7.7mm，角膜后表面的曲率半径为6.8mm，并且角膜中央厚度为0.5mm.巩膜构成眼睛外纤维膜的后5/6，主要由胶原纤维组成.根据解剖和生理的正常值，巩膜的水平直径为23.5mm，垂直直径为23mm，前后直径为24.4mm.巩膜厚度有区域差异，通常赤道部分较薄，后极部最厚[18].考虑到计算成本和眼部组织的复杂性，在该模型中忽略了晶状体，虹膜和其他组织的建模.根据Gullstrand精确眼模型的数据，在合理简化和假设的基础上，采用参数化建模方法，并使用三维建模软件SiemensNX来构建眼球模型.将构建的三维实体模型导入COMSOL软件中的三维-固体力学-稳态分析功能模块进行相关分析.
1.2模型选择
角膜具有复杂的生物力学特性，是一种典型的粘弹性组织[19].因而加载后，角膜表现出一些瞬时变形(纯弹性行为)，然后是渐进变形(粘弹性行为)[20-21].这些行为可以用相当简单的弹簧和阻尼系统进行建模.减震器表示对施加力的时滞粘性阻力，而弹簧仅表示纯弹性行为.常用的粘弹性模型是Kelvin-Voigt，Maxwell和StandardLinearSolid模型[22].这种带有附加弹簧的Kelvin-Voigt模型，不仅是因为其简单性，而且还因为它可以清晰，直观地分离角膜的纯弹性和粘弹性应变响应.这类似于Kobayashi等[21]的配置，他们使用两个Kelvin-Voigt模型与弹簧串联来比较快速和慢速粘弹性响应.因为它不具有纯弹性分量，不允许任何瞬时变形，所以我们未选择纯Kelvin-Voigt模型.而未选择Maxwell模型是因为它在施加载荷时会无限期地蠕变.在角膜中，蠕变接近渐近线[23].Maxwell模型的另一个局限性在于它不能从纯粘性成分引起的变形中完全恢复.
1.3材料属性
1.3.1未增加粘弹性
为显示方便，角膜中Y轴方向上的圆形部分是区域Ⅰ，为病变组织;角膜剩余部分是区域Ⅱ，为正常组织;巩膜部分是区域Ⅲ，为正常组织.整个眼球组织网格密度经过放大处理，如图2所示.在本文中，第2部分为模拟结果图，自变量都是Y轴坐标，最小是-5mm，最大是5mm.通过改变角膜正常和病变两个区域的杨氏模量比值，可观察到角膜病变组织弹性差异导致的应力与应变的变化.
材料属性的设置如表1所示.角膜病变组织是表1中的区域I，区域I-1是角膜病变组织杨氏模量大于角膜正常组织的情况、区域I-2是角膜病变组织杨氏模量小于角膜正常组织的情况和区域I-3是角膜组织杨氏模量按高斯分布的情况;区域2是角膜正常组织，杨氏模量的值保持不变.特别要注意的是，生物软组织的杨氏模量在实际中不会有很大的突变，因此高斯分布更符合实际.角膜正常组织和病变组织的密度均设置为1076kg·m-3，泊松比均设为0.49.巩膜组织的密度设为1243kg·m-3，泊松比设为0.49.
1.3.2增加粘弹性
求解结果是非对称的，因此选择三维-固体力学-瞬态分析功能模块进行模拟.几何模型设置、材料属性设置和边界设定与第1.3.1小节设置情况基本相同，但采用线性增大的瞬态激励且材料上增加了粘弹性[24]，如表2所示.
1.4边界设定
未增加粘弹性时，自角膜顶部向下施加1%的应变，为了阻止眼球转动，在巩膜底部(直径约为4mm)进行零位移约束.本文所有模型均只对角膜后表面进行相关应变、应力和位移分析，且均在15mmHg眼压下模拟，贴合正常人眼的眼内压.在本文中，三维情况描述的是正常角膜组织中含有一个病变的组织，且据此设置对应的材料参数.根据几何模型的特点，采用四面体网格进行划分，网格密度满足计算的精度和准确性.增加粘弹性后求解器的设置，时间范围(0，0.1，12)，相对容差0.01，采取MUMPS求解器.为了防止计算时内存溢出，网格划分精度相比未增加粘弹性时略有降低，以满足正常计算需要.
1.5有限元求解
超声弹性成像正向问题是了解组织内部弹性分布和边界条件，并计算组织中的应力应变分布[11].在本文中，有限元法被用来模拟弹性成像.有限元法是将弹性体划分为有限个有限大小的构件，并把有限个节点相互连接组成的有限大小的构件称为单元，这些单元之间通过有限个节点相互连接，单元之间的力通过节点传递[25].操作过程包括网格生成、单元计算、全局计算、代码存储和其他步骤.其中，单元分析指的是建立节点位移与单元中任何点之间的关系：位移、单元应力、单元应变和单元节点力之间的关系，这是整个有限元计算过程中最重要的部分[1].
选择COMSOL软件进行有限元分析.COMSOLMultiphysics是一款先进的基于有限元方法的数值模拟软件，通过求解过程中的偏微分方程(单场)或偏微分方程组(多场)，可以实现物理现象的精确数值模拟.
2结果
2.1病变组织为硬块时的模拟结果
图3是角膜病变组织杨氏模量值为4.44MPa(硬块)时的模拟结果.图3(a)是角膜后表面二维应变图、图3(b)是角膜后表面Y轴上的应变图、图3(c)是角膜后表面二维应力图、图3(d)是角膜后表面Y轴上的应力图、图3(e)是角膜后表面二维位移图、图3(f)是角膜后表面Y轴上的位移图.二维图中位移单位为毫米(mm)。
从图3可以看出，当角膜正常组织与病变组织的杨氏模量之比为1:4时，即区域I的杨氏模量为4.44MPa，区域II的杨氏模量为1.11MPa.从图中也能明显的看出有几个应变、应力高峰值，这是由于外部激励的存在引起的.对应于角膜组织病变为硬块时，结果表明，当施加在角膜顶部的固定应变为1%时，角膜正常组织与病变组织的应变和应力差异较明显，位移则差异较小.从图3(a)和图3(b)中可以看出，角膜病变组织的应变绝对值小于正常组织的应变.病变区域中心的应变为-0.00854，而对应正常角膜相同位置的应变为-0.02277，为病变区域中心应变值的2.67倍.从图3(c)和图3(d)中可以看出，角膜病变组织的应力大于正常组织的应力.病变区域中心的应力为0.04337MPa，而对应正常角膜相同位置的应力为0.02729MPa，相当于病变区域中心的应力值的0.64倍.在图3(e)中能够很清楚的观察到组织内部位移分布不均匀.同样的，在图3(f)中，角膜后表面Y轴上的位移分布也是不均匀的，且病变区域的组织位移值比正常区域的要小一些.角膜病变组织的边缘的应力和应变都有很大的跳变.该模拟结果是弹性成像的一个正向问题，它验证了胡克定律，并阐明了将应力应变关系应用于弹性成像的可行性.
2.2病变组织为囊肿时的模拟结果
图4是角膜病变组织杨氏模量值为0.278MPa(囊肿)时的模拟结果.图4(a)是角膜后表面二维应变图、图4(b)是角膜后表面Y轴上的应变图、图4(c)是角膜后表面二维应力图、图4(d)是角膜后表面Y轴上的应力图、图4(e)是角膜后表面二维位移图、图4(f)是角膜后表面Y轴上的位移图.
从图4可以看出，当正常组织与病变组织的杨氏模量之比为4:1时，即区域I的杨氏模量为0.278MPa，区域II的杨氏模量为1.11MPa.对应于角膜组织病变为囊肿时，结果表明，当施加在角膜顶部的固定应变为1%时，角膜正常组织与病变组织的应变、应力差异较明显，而位移的差异仍较小.从图4(a)和图4(b)中可以看出，角膜病变组织的应变绝对值大于正常组织的应变.病变区域中心的应变为-0.03949，而对应正常角膜相同位置的应变为-0.02277，相当于病变区域中心的应变值的0.58倍.从图4(c)和图4(d)中可以看出，角膜病变组织的应力小于正常组织的应力.病变区域中心的应力为0.01152MPa，而对应正常角膜相同位置的应力为0.02729MPa，为病变区域中心的应力值的2.37倍.从图4(e)和图4(f)中可以看出，病变区域的组织位移值比正常区域的要大一些.在角膜病变组织的边缘的应力和应变仍都有很大的跳变.——论文作者：陶兴明1,2,3，方利华1,2,3，林鲁超1,2,3，何兴道1,2,3，伏燕军1,2,3，黄采敏4
相关期刊推荐：《力学季刊》(季刊)(原《上海力学》)创刊于1980年，反映我国高等院校、研究院所和工程技术单位的前沿力学研究的成果和动态，报道力学在解决我国重大工程设计和技术问题中的进展和作用。主要刊登：理论研究、实验研究、计算研究和应用研究等的学术论文及简报，以及国际力学研究新方向的综述和展望;同时刊登学术讨论、新书评介和教学研究。
