随着电子产品的日益增多，我国近视人群数量庞大。治疗近视非常流行的一种方法是角膜屈光手术，如准分子激光原位角膜磨镶术(LASIK手术)、飞秒激光辅助准分子激光原位角膜磨镶术(Flex手术)、小切口飞秒激光基质透镜取出术(SMILE手术)等[1]。角膜的生物力学性能对手术方案的选取有着重要的影响，同时对术后角膜变形的预测以及安全性的评估有着重要的意义。
国内外学者使用多种方法来测量角膜生物力学性能，诸如轴向拉伸实验、膨胀实验和压痕实验等[2-3]。由于人眼角膜的获取较为困难，大多数人眼角膜实验中的材料均取自手术剩余边角料，且都是在角膜脱离人体之后进行的测量，实验过程全部伴随着破坏性手段[4-5]，这样的实验条件不能代表角膜真实的生理环境和力学环境，测得的力学参数也很难应用于临床当中。目前在临床中应用广泛的在体检测设备是可视化角膜生物力学分析仪(Corvis ST)，该设备系统利用高速成像仪记录角膜在气流作用下的变形过程，可重复且有效测得角膜厚度、眼内压、角膜顶点位移随时间的变化曲线、角膜变形幅度、压平时间等多个角膜在体生物力学参数[6]，是目前十分推崇的在体检测设备。
有限元的发展对角膜生物力学研究起到了重要的推动作用，在屈光手术、眼外伤、青光眼以及角膜扩张等眼部疾病的研究中往往需要有限元模型进行研究。角膜的屈光力和角膜形态有很大关系，国内外研究中常借助有限元模拟的方法研究角膜的生物力学性能。PANDOLFI et al[7]将角膜视作椭球面，在柱坐标系中用椭球面近似近视眼角膜和散光角膜表面，用球面方程计算角膜屈光度；NGUYEN et al[8]建立了二维轴对称角膜模型，通过comsol模拟了Corvis ST检测过程研究角膜变形；GEFEN et al[9]建立了理想的完全轴对称正常角膜和圆锥角膜，分析了不同眼压作用下角膜屈光和表面应力的变化。以往的角膜几何模型常根据角膜的解剖学参数建立[10-11]，即角膜的前后表面曲率、高度、横径等参数采用解剖学平均值，角膜厚度方向上处理为均匀厚度；但根据解剖学参数建立的角膜模型与真实角膜形态有一定差距，无法准确还原角膜表面不规则的曲率变化。
三维眼前节分析仪(Pentacam)是广泛应用于眼部疾病检查的非接触式检测仪，检测系统通过旋转扫描角膜，测得角膜前后表面形态、曲率数据以及厚度等参数，并生成彩色图像，检测过程快捷、无损伤，临床应用十分广泛。该系统测量覆盖范围广，有较高的准确度和良好的重复性，对角膜形态学特征描述比较全面，并且可定点描述角膜某处的形态学特征，借助三维眼前节分析仪的检测数据可以较准确地还原角膜的几何形态[12]。本文结合Corvis ST检测数据和Pentacam地形数据对20例近视眼角膜进行有限元分析计算，求解出对应角膜弹性模量。
1 方法
1.1 数据采集
经山西省眼科医院伦理委员会批准通过，选取山西省眼科医院准分子激光科确诊的近视患者20例。选取的20例受试者中，男12例(12眼)，女8例(8眼)，受试者年龄均值为25.02±4.60岁，均为大专以上受教育程度；健康状况良好，无慢性病、免疫性疾病、精神疾病及其他疾病；无眼部手术史或外伤史，无圆锥角膜、青光眼等其他眼部疾病；屈光度稳定两年及以上，屈光度年增长不超过0.5 D；双眼柱镜度数不超过-4.0 D，等效球镜度数不高于-10.0 D；角膜厚度大于470 μm，眼压小于3.192 kPa；停戴各类接触镜和塑形镜一月以上。所有患者经三维眼前节分析仪和可视化角膜生物力学分析仪进行检测，检查均由同一专业人员操作，每眼检测3次，取图像最佳、重复性最好的结果分析。从眼前节分析仪中获取患者角膜前后表面高度图、厚度图等地形数据；从Corvis ST中获取检测过程中角膜真实顶点位移、中央角膜厚度和眼内压等生物力学参数。Corvis ST检测结果中提供了IOPnct和bIOP(biomechanical corrected IOP)两种眼压值，bIOP为修正眼内压，排除了其他因素对眼压测量的影响，更接近真实眼压值。本文采用修正眼内压bIOP值，眼压单位用kPa表示；Corvis ST还同时提供了角膜顶点位移曲线(deflection amplitude)和包含了眼球位移的角膜顶点位移曲线(deformation amplitude)，本文所建模型仅包含角膜部分，不包含除角膜外的其他部分，因此角膜顶点位移选用deflection amplitude值。
1.2 建立角膜几何模型
采集Pentacam地形仪8 mm×8 mm范围内角膜表面上各点与最适球面间的径向距离差值(ddiff)、角膜前表面上各点处对应的角膜厚度值和最适球面半径等数据。建立如图1(a)所示三维空间直角坐标系[13-14]，将角膜顶点置于空间坐标系原点处。Pentacam提供了角膜前后表面8 mm×8 mm(XOY平面)范围内各点横、纵坐标(x,y)以及对应ddiff值，但未提供各点Z方向坐标。需通过几何关系计算角膜前、后表面各点Z坐标值。图1(b)、(c)中，最适球面半径(AD段)为r，A点坐标为(0,0,-r)，角膜表面上某点C与最适球面球心A的连线或连线的延长线与最适球面相交于D点，CD间长度即ddiff值。ddiff值小于0表示角膜表面低于最适球面，对应图1(b)；ddiff值大于0表示角膜表面高于最适球面，对应图1(c).d为角膜上某点到Z轴的距离(BC段)，角膜前表面上某点的Z坐标值表示为：
Z=-(r-((r+ddiff)2-d2)1/2) .
式中d2=x2+y2.将求得的前表面Z坐标值减去对应点的厚度值即为后表面Z坐标值。将计算得到的前后表面三维坐标点数据导入逆向建模软件Geomagic 2013中，形成角膜前后表面点云轮廓图，对点云封装生成*iges格式曲面片，对前后表面曲面片缝合并填充生成实体，将实体导入有限元计算软件Abaqus 6.14中，建模过程如图2所示。Pentacam检测的是角膜在眼内压作用下的坐标位置，并非角膜在无眼压状态下的位置，因此需要通过迭代的方法找到角膜在无眼压状态下的坐标点，进而建立无应力状态下的几何模型，参考ARIZA-GRACIA et al[15]的方法进行计算。首先对角膜上各点所有Z坐标值减去一个微小量δ，δ初始值设置为0.005，在其基础上施加眼内压并进行运算，通过Abaqus运算后的*inp文件导出角膜表面各点的Z坐标值并求和，求和结果与仪器检测值相减，当差值小于10-3时计算结束，通过不断尝试找到角膜零应力状态下的坐标，本文采集的患者数据较少，故所有模型均采用尝试法计算角膜的初始状态坐标值。
图1 三维空间直角坐标系
Fig.1 Three dimensional rectangular coordinates
图2 角膜模型建模过程
Fig.2 Processing of cornea model modeling
1.3 角膜模拟检测过程
角膜材料采用线弹性本构模型描述，在Abaqus中需设置弹性模量和泊松比。角膜通常视为不可压缩性材料[16]，在有限元计算中泊松比设置为0.49，弹性模量逆向求解。模拟角膜在Corvis ST检测时的变形过程[17]，角膜内部施加眼内压均布载荷，外部施加气流均布载荷，气流载荷采用WANG et al[18]测得的数值，固定角膜边缘的转角和位移(如图3所示)。角膜先后历经压平、凹陷、再次压平和回弹的过程，模拟变形过程如图4所示。模型网格采用四面体单元，考虑到计算规模及收敛性，对网格密度进行了验证。
图3 角膜约束条件和载荷
Fig.3 Corneal constraints and loads
图4 角膜变形模拟过程
Fig.4 Simulation processing of corneal deformation
1.4 统计学分析
使用SPSS 20.0对结果进行统计学分析，描述性统计结果用均值±标准差进行描述。
2 结果
模拟20例角膜的Corvis ST检测过程，通过对比模拟与检测的角膜顶点位移曲线，反推出每个角膜的弹性模量值，图5给出了部分角膜顶点位移与时间关系的模拟与检测曲线对比图，图中红色曲线为Corvis ST检测曲线，黑色曲线为模拟曲线，所有角膜对比曲线均类似，不一一列举。图6给出了角膜在不同压平阶段的有限元模型轮廓与实际检测图像轮廓的对比图，图中可看出不同阶段角膜模型变形与实际变形基本一致，验证了模型的有效性。
图5 模拟与检测角膜顶点位移-时间曲线对比
Fig.5 Comparision of corneal vertex displacement-time curves between simulation and detection
图6 角膜模拟轮廓与检测轮廓对比图
Fig.6 Comparison of corneal simulated contour
and detected contour
20例角膜的中央角膜厚度均值为531.68 μm±32.38 μm，眼内压均值为1.70 kPa±0.29 kPa.通过模拟反推出角膜弹性模量均值为0.822 MPa±0.146 MPa，20例角膜弹性模量结果如表1所示。
表1 不同患者眼压和角膜弹性模量
Table 1 IOP and corneal elastic modulus
of different patients
NumberIOP/kPaYoung’s modulus/MPa11.50.75221.70.99031.30.64341.50.81151.90.76161.20.61371.50.99781.80.96091.60.694101.71.012NumberIOP/kPaYoung’s modulus/MPa111.70.866121.60.734131.51.104141.60.812151.50.940161.70.655171.90.857181.80.899191.70.635201.30.700
3 讨论
角膜材料的力学性能对于研究角膜屈光术等角膜病变至关重要。对于结构十分精密的角膜组织而言，角膜表面曲率的微小变化都会严重影响视力，若角膜前后表面的几何参数设置为固定值，将与角膜实际情况产生很大误差，可能对角膜材料性能分析产生较大影响。本文通过Pentacam地形数据建立角膜几何模型，较为准确地还原了近视眼角膜真实几何形态，在一定程度上降低了因角膜几何形状产生的模拟误差。
不同方法、取材和检测手段得出的角膜弹性模量值均有差别。黄来鑫等[19]通过建立理想的近视患者眼球有限元模型，结合Corvis ST结果计算了8例近视眼角膜弹性模量，均值为0.68 MPa，与本文均值0.822 MPa的结果较接近。薛超等[5]对SMILE手术取出的基质透镜进行单轴拉伸实验，当基质应力为0.02 MPa时，对应的弹性模量为1.26MPa± 0.71 MPa，略高于本文计算值；XIANG et al[20]对10例SMILE手术的角膜基质透镜在水平和竖直两个方向上进行单轴拉伸测试，测得弹性模量分别为1.300 MPa±0.508 MPa、1.142 MPa± 0.238 MPa，也略高于本文计算值；这可能是由于其实验对象为角膜基质部分，而非全角膜材料，且在离体环境下实验也会产生一定的误差。
角膜组织本身为超弹性材料，以往研究中有用HGO模型描述角膜的材料属性，其中需要待定的参数有4个，分别为：k1(代表胶原纤维刚度)、k2(代表胶原纤维非线性的无量纲常数)、κ(表示胶原纤维的分散程度，0≤κ<1/3)和C10(表示基质相关材料参数)。材料参数的确定是研究者根据实验拟合出的结果，但是对于某个实验结果，可以拟合出许多种与实验曲线相符合的参数组合结果。天津大学向尧齐[21]用SMILE手术取出的角膜透镜进行单轴拉伸实验，对34组实验结果通过最小二乘法拟合出HGO模型中的参数，C10、k1、k2的均值分别为0.220、0.615、121.633，与ARIZA-GRACIA[15]拟合的结果(C10=0.05、k1=130.9、k2=2490)相差较大；向尧齐[21]还提出κ值的改变对角膜应力-应变的影响非常大，也会影响k1和k2的取值。而每个角膜的κ值均不同，这就造成了拟合出的其他参数结果可能并不准确。不同的角膜取材、实验方法和拟合方法都有可能对拟合参数的结果产生非常大的影响，以至于不同研究者的研究结果相差甚大。目前通过拉伸实验测定角膜线弹性材料属性的研究较成熟，在仿真模拟中待定参数少。故本文采用线弹性本构描述角膜材料属性。此外由于Corvis ST检测为动态实验，角膜材料具有粘弹性特性，粘弹性对实验结果存在一定的影响，但由于目前对角膜粘弹性的本构模型研究较少，本文未考虑角膜粘弹性对试验结果的影响，是本研究的一个局限。
有限元模型是理想化的几何和力学模型，与人体真实情况有一定的差别。首先，在角膜几何模型方面，因很难获取全部巩膜和眼球内容物的准确几何形状及位置信息，建模时仅考虑了角膜组织，忽略了眼内容物等其他组织成分对边界条件以及模拟结果产生的影响；角膜中包含的5层组织为各向异性材料，且力学性能各不相同，根据DIAS et al[4]的研究，角膜前、后基质的弹性模量并不相同，本模型未对角膜各部分分别赋予材料属性；不同几何形态的角膜表面在气流载荷作用时压力并不完全相同，本文假定了所有角膜受到的气流压力值均相同。其次，同一患者在进行多次检测时，眼压、角膜顶点位移等检测结果也有一定程度差别，这也可能造成模拟结果与真实情况之间的误差。最后，本文选取的20例近视患者样本数量有限，并不能代表全部近视患者人群，且考虑到其他个体差异原因，本研究所求角膜弹性模量值并不具有普适性，仅可作为参考值。
根据本文所求近视眼角膜弹性模量，可为建立近视角膜有限元模型提供一定的参考，以期为后续更准确地模拟角膜屈光手术及预测术后角膜力学性能提供理论基础。
[1] GYLDENKERNE A,IVARSEN A,JESPER H.Comparison of corneal shape changes and aberrations induced by FS-LASIK and SMILE for myopia[J].Journal of Refractive Surgery,2015,31(4):1-7.
[2] ELSHEIKH A,WHITFORD C,HAMARASHID R,et al.Stress free configuration of the human eye[J].Medical Engineering & Physics,2013,35(2):211-216.
[3] WU K S,LI S Y,LO A C Y,et al.Micro-scale stiffness change of cornea tissues suffered from elevated intraocular pressure investigated by nanoindentation[J].Soft Materials, 2013,11(3):244-253.
[4] DIAS J M ,ZIEBARTH N M .Anterior and posterior corneal stroma elasticity assessed using nanoindentation[J].Experimental Eye Research,2013,115:41-46.
[5] 薛超,王雁,向尧齐,等.近视患者角膜在体与离体生物力学特性相关性初步研究[J].眼科新进展,2021,41(10):943-947.
XUE C,WANG Y,XIANG Y Q,et al.Correlation of corneal biomechanical properties of myopic patients in vivo and in vitro[J].Recent Advances in Ophthalmology,2021,41(10):943-947.
[6] MARIA L S,MARCO Z,CLAUDIA T,et al.Corneal deformation parameters provided by the corvis-ST pachy-tonometer in healthy subjects and glaucoma patients[J].Journal of Glaucoma,2015,24(8):568-574.
[7] PANDOLFI A,MANGANIELLO F.A model for the human cornea:constitutive formulation and numerical analysis[J].Biomechanics & modeling in mechanobiology,2006,5(4):237-246.
[8] NGUYEN B A,ROBERTS C J,REILLY M A.Biomechanical impact of the sclera on corneal deformation response to an air-puff:a finite-element study[J].Frontiers in Bioengineering & Biotechnology,2019,6:210.
[9] GEFEN A,SHALOM R,ELAD D,et al.Biomechanical analysis of the keratoconic cornea[J].Journal of the Mechanical Behavior of Biomedical Materials,2009,2(3):224-236.
[10] QIN X,YU M,ZHANG H,et al.The mechanical interpretation of ocular response analyzer parameters[J].BioMed Research International,2019(1):1-11.
[11] QIN X,TIAN L,ZHANG H,et al.Evaluation of corneal elastic modulus based on corneal visualization scheimpflug technology[J].Biomedical Engineering Online,2019,18(1):42.
[12] ZHENG T,XU J,LU Y.Comparison of the accuracy of four Pentacam corneal astigmatism values in non-toric pseudophakic eyes[J].Albrecht von Graaees Archiv für Ophthalmologie,2020,258(4)：795-803.
[13] BAO F,CHEN H,YU Y,et al.Evaluation of the shape symmetry of bilateral normal corneas in a Chinese population[J].PloS One,2013,8(8):e73412.
[14] 包芳军.LASIK数字预测模型的临床验证及眼生物力学测量平台的创建[D].温州：温州医科大学,2015.
[15] ARIZA-GRACIA M
,ZURITA J F,PI
ERO D P,et al.Coupled biomechanical response of the cornea assessed by non-contact tonometry.A simulation study[J].PloS One,2015,10(3):e0121486.
[16] LAND M F.The human eye:structure and function[J].Nature Medicine,1999,5(11):1229.
[17] 孟乔宇,王晓君,陈维毅,等.Corvis ST检测在体角膜生物力学性能的有限元模拟[J].生物医学工程学杂志,2020,37(4):608-613.
MENG Q Y,WANG X J,CHEN W Y,et al.Finite element analysis of determining corneal biomechanical properties in vivo based on Corvis ST[J].Journal of Biomedical Engineering,2020,37(4):608-613.
[18] WANG L K,TIAN L,ZHENG Y P.Determining in vivo elasticity and viscosity with dynamic scheimpflug imaging analysis in keratoconic and healthy eyes[J].Journal of Biophotonics,2016,9(5):454-463.
[19] 黄来鑫,沈珉,王雁.基于Corvis ST测量人角膜的弹性模量[J].应用力学学报,2021,38(3):861-867.
HUANG L X,SHEN M,WANG Y.Inverse solution of elastic modulus of human cornea based on Corvis ST[J].Chinese Journal of Applied Mechanics,2021,38(3):861-867.
[20] XIANG Y Q,SHEN M,XUE C,et al.Tensile biomechanical properties and constitutive parameters of human corneal stroma extracted by SMILE procedure[J].Journal of the mechanical behavior of biomedical materials,2018,85:102-108.
[21] 向尧齐.基于单轴拉伸实验的人角膜的生物力学性能的研究[D].天津：天津大学,2018.