孟喬宇，王曉君，李曉娜，賀 瑞，孟 瑩

(1.北京理工大學(xué) 先進(jìn)結(jié)構(gòu)技術(shù)研究院，北京 100081；2.太原理工大學(xué) a.機(jī)械與運(yùn)載工程學(xué)院，b.生物醫(yī)學(xué)工程學(xué)院，太原 030024；3.山西省眼科醫(yī)院 準(zhǔn)分子激光科，太原 030002；4.天津城建大學(xué) 控制與機(jī)械工程學(xué)院，天津 300384)

隨著電子產(chǎn)品的日益增多，我國近視人群數(shù)量龐大。治療近視非常流行的一種方法是角膜屈光手術(shù)，如準(zhǔn)分子激光原位角膜磨鑲術(shù)(LASIK手術(shù))、飛秒激光輔助準(zhǔn)分子激光原位角膜磨鑲術(shù)(Flex手術(shù))、小切口飛秒激光基質(zhì)透鏡取出術(shù)(SMILE手術(shù))等[1]。角膜的生物力學(xué)性能對手術(shù)方案的選取有著重要的影響，同時對術(shù)后角膜變形的預(yù)測以及安全性的評估有著重要的意義。

國內(nèi)外學(xué)者使用多種方法來測量角膜生物力學(xué)性能，諸如軸向拉伸實驗、膨脹實驗和壓痕實驗等[2-3]。由于人眼角膜的獲取較為困難，大多數(shù)人眼角膜實驗中的材料均取自手術(shù)剩余邊角料，且都是在角膜脫離人體之后進(jìn)行的測量，實驗過程全部伴隨著破壞性手段[4-5]，這樣的實驗條件不能代表角膜真實的生理環(huán)境和力學(xué)環(huán)境，測得的力學(xué)參數(shù)也很難應(yīng)用于臨床當(dāng)中。目前在臨床中應(yīng)用廣泛的在體檢測設(shè)備是可視化角膜生物力學(xué)分析儀(Corvis ST)，該設(shè)備系統(tǒng)利用高速成像儀記錄角膜在氣流作用下的變形過程，可重復(fù)且有效測得角膜厚度、眼內(nèi)壓、角膜頂點位移隨時間的變化曲線、角膜變形幅度、壓平時間等多個角膜在體生物力學(xué)參數(shù)[6]，是目前十分推崇的在體檢測設(shè)備。

有限元的發(fā)展對角膜生物力學(xué)研究起到了重要的推動作用，在屈光手術(shù)、眼外傷、青光眼以及角膜擴(kuò)張等眼部疾病的研究中往往需要有限元模型進(jìn)行研究。角膜的屈光力和角膜形態(tài)有很大關(guān)系，國內(nèi)外研究中常借助有限元模擬的方法研究角膜的生物力學(xué)性能。PANDOLFI et al[7]將角膜視作橢球面，在柱坐標(biāo)系中用橢球面近似近視眼角膜和散光角膜表面，用球面方程計算角膜屈光度；NGUYEN et al[8]建立了二維軸對稱角膜模型，通過comsol模擬了Corvis ST檢測過程研究角膜變形；GEFEN et al[9]建立了理想的完全軸對稱正常角膜和圓錐角膜，分析了不同眼壓作用下角膜屈光和表面應(yīng)力的變化。以往的角膜幾何模型常根據(jù)角膜的解剖學(xué)參數(shù)建立[10-11]，即角膜的前后表面曲率、高度、橫徑等參數(shù)采用解剖學(xué)平均值，角膜厚度方向上處理為均勻厚度；但根據(jù)解剖學(xué)參數(shù)建立的角膜模型與真實角膜形態(tài)有一定差距，無法準(zhǔn)確還原角膜表面不規(guī)則的曲率變化。

三維眼前節(jié)分析儀(Pentacam)是廣泛應(yīng)用于眼部疾病檢查的非接觸式檢測儀，檢測系統(tǒng)通過旋轉(zhuǎn)掃描角膜，測得角膜前后表面形態(tài)、曲率數(shù)據(jù)以及厚度等參數(shù)，并生成彩色圖像，檢測過程快捷、無損傷，臨床應(yīng)用十分廣泛。該系統(tǒng)測量覆蓋范圍廣，有較高的準(zhǔn)確度和良好的重復(fù)性，對角膜形態(tài)學(xué)特征描述比較全面，并且可定點描述角膜某處的形態(tài)學(xué)特征，借助三維眼前節(jié)分析儀的檢測數(shù)據(jù)可以較準(zhǔn)確地還原角膜的幾何形態(tài)[12]。本文結(jié)合Corvis ST檢測數(shù)據(jù)和Pentacam地形數(shù)據(jù)對20例近視眼角膜進(jìn)行有限元分析計算，求解出對應(yīng)角膜彈性模量。

1 方法

1.1 數(shù)據(jù)采集

經(jīng)山西省眼科醫(yī)院倫理委員會批準(zhǔn)通過，選取山西省眼科醫(yī)院準(zhǔn)分子激光科確診的近視患者20例。選取的20例受試者中，男12例(12眼)，女8例(8眼)，受試者年齡均值為25.02±4.60歲，均為大專以上受教育程度；健康狀況良好，無慢性病、免疫性疾病、精神疾病及其他疾病；無眼部手術(shù)史或外傷史，無圓錐角膜、青光眼等其他眼部疾??；屈光度穩(wěn)定兩年及以上，屈光度年增長不超過0.5 D；雙眼柱鏡度數(shù)不超過-4.0 D，等效球鏡度數(shù)不高于-10.0 D；角膜厚度大于470 μm，眼壓小于3.192 kPa；停戴各類接觸鏡和塑形鏡一月以上。所有患者經(jīng)三維眼前節(jié)分析儀和可視化角膜生物力學(xué)分析儀進(jìn)行檢測，檢查均由同一專業(yè)人員操作，每眼檢測3次，取圖像最佳、重復(fù)性最好的結(jié)果分析。從眼前節(jié)分析儀中獲取患者角膜前后表面高度圖、厚度圖等地形數(shù)據(jù)；從Corvis ST中獲取檢測過程中角膜真實頂點位移、中央角膜厚度和眼內(nèi)壓等生物力學(xué)參數(shù)。Corvis ST檢測結(jié)果中提供了IOPnct和bIOP(biomechanical corrected IOP)兩種眼壓值，bIOP為修正眼內(nèi)壓，排除了其他因素對眼壓測量的影響，更接近真實眼壓值。本文采用修正眼內(nèi)壓bIOP值，眼壓單位用kPa表示；Corvis ST還同時提供了角膜頂點位移曲線(deflection amplitude)和包含了眼球位移的角膜頂點位移曲線(deformation amplitude)，本文所建模型僅包含角膜部分，不包含除角膜外的其他部分，因此角膜頂點位移選用deflection amplitude值。

1.2 建立角膜幾何模型

采集Pentacam地形儀8 mm×8 mm范圍內(nèi)角膜表面上各點與最適球面間的徑向距離差值(ddiff)、角膜前表面上各點處對應(yīng)的角膜厚度值和最適球面半徑等數(shù)據(jù)。建立如圖1(a)所示三維空間直角坐標(biāo)系[13-14]，將角膜頂點置于空間坐標(biāo)系原點處。Pentacam提供了角膜前后表面8 mm×8 mm(XOY平面)范圍內(nèi)各點橫、縱坐標(biāo)(x,y)以及對應(yīng)ddiff值，但未提供各點Z方向坐標(biāo)。需通過幾何關(guān)系計算角膜前、后表面各點Z坐標(biāo)值。圖1(b)、(c)中，最適球面半徑(AD段)為r，A點坐標(biāo)為(0,0,-r)，角膜表面上某點C與最適球面球心A的連線或連線的延長線與最適球面相交于D點，CD間長度即ddiff值。ddiff值小于0表示角膜表面低于最適球面，對應(yīng)圖1(b)；ddiff值大于0表示角膜表面高于最適球面，對應(yīng)圖1(c).d為角膜上某點到Z軸的距離(BC段)，角膜前表面上某點的Z坐標(biāo)值表示為：

Z=-(r-((r+ddiff)2-d2)1/2) .

式中d2=x2+y2.將求得的前表面Z坐標(biāo)值減去對應(yīng)點的厚度值即為后表面Z坐標(biāo)值。將計算得到的前后表面三維坐標(biāo)點數(shù)據(jù)導(dǎo)入逆向建模軟件Geomagic 2013中，形成角膜前后表面點云輪廓圖，對點云封裝生成*iges格式曲面片，對前后表面曲面片縫合并填充生成實體，將實體導(dǎo)入有限元計算軟件Abaqus 6.14中，建模過程如圖2所示。Pentacam檢測的是角膜在眼內(nèi)壓作用下的坐標(biāo)位置，并非角膜在無眼壓狀態(tài)下的位置，因此需要通過迭代的方法找到角膜在無眼壓狀態(tài)下的坐標(biāo)點，進(jìn)而建立無應(yīng)力狀態(tài)下的幾何模型，參考ARIZA-GRACIA et al[15]的方法進(jìn)行計算。首先對角膜上各點所有Z坐標(biāo)值減去一個微小量δ，δ初始值設(shè)置為0.005，在其基礎(chǔ)上施加眼內(nèi)壓并進(jìn)行運(yùn)算，通過Abaqus運(yùn)算后的*inp文件導(dǎo)出角膜表面各點的Z坐標(biāo)值并求和，求和結(jié)果與儀器檢測值相減，當(dāng)差值小于10-3時計算結(jié)束，通過不斷嘗試找到角膜零應(yīng)力狀態(tài)下的坐標(biāo)，本文采集的患者數(shù)據(jù)較少，故所有模型均采用嘗試法計算角膜的初始狀態(tài)坐標(biāo)值。

圖1 三維空間直角坐標(biāo)系Fig.1 Three dimensional rectangular coordinates

圖2 角膜模型建模過程Fig.2 Processing of cornea model modeling

1.3 角膜模擬檢測過程

角膜材料采用線彈性本構(gòu)模型描述，在Abaqus中需設(shè)置彈性模量和泊松比。角膜通常視為不可壓縮性材料[16]，在有限元計算中泊松比設(shè)置為0.49，彈性模量逆向求解。模擬角膜在Corvis ST檢測時的變形過程[17]，角膜內(nèi)部施加眼內(nèi)壓均布載荷，外部施加氣流均布載荷，氣流載荷采用WANG et al[18]測得的數(shù)值，固定角膜邊緣的轉(zhuǎn)角和位移(如圖3所示)。角膜先后歷經(jīng)壓平、凹陷、再次壓平和回彈的過程，模擬變形過程如圖4所示。模型網(wǎng)格采用四面體單元，考慮到計算規(guī)模及收斂性，對網(wǎng)格密度進(jìn)行了驗證。

圖3 角膜約束條件和載荷Fig.3 Corneal constraints and loads

圖4 角膜變形模擬過程Fig.4 Simulation processing of corneal deformation

1.4 統(tǒng)計學(xué)分析

使用SPSS 20.0對結(jié)果進(jìn)行統(tǒng)計學(xué)分析，描述性統(tǒng)計結(jié)果用均值±標(biāo)準(zhǔn)差進(jìn)行描述。

2 結(jié)果

模擬20例角膜的Corvis ST檢測過程，通過對比模擬與檢測的角膜頂點位移曲線，反推出每個角膜的彈性模量值，圖5給出了部分角膜頂點位移與時間關(guān)系的模擬與檢測曲線對比圖，圖中紅色曲線為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，眼內(nèi)壓均值為1.70 kPa±0.29 kPa.通過模擬反推出角膜彈性模量均值為0.822 MPa±0.146 MPa，20例角膜彈性模量結(jié)果如表1所示。

表1 不同患者眼壓和角膜彈性模量Table 1 IOP and corneal elastic modulus of different patients

3 討論

角膜材料的力學(xué)性能對于研究角膜屈光術(shù)等角膜病變至關(guān)重要。對于結(jié)構(gòu)十分精密的角膜組織而言，角膜表面曲率的微小變化都會嚴(yán)重影響視力，若角膜前后表面的幾何參數(shù)設(shè)置為固定值，將與角膜實際情況產(chǎn)生很大誤差，可能對角膜材料性能分析產(chǎn)生較大影響。本文通過Pentacam地形數(shù)據(jù)建立角膜幾何模型，較為準(zhǔn)確地還原了近視眼角膜真實幾何形態(tài)，在一定程度上降低了因角膜幾何形狀產(chǎn)生的模擬誤差。

不同方法、取材和檢測手段得出的角膜彈性模量值均有差別。黃來鑫等[19]通過建立理想的近視患者眼球有限元模型，結(jié)合Corvis ST結(jié)果計算了8例近視眼角膜彈性模量，均值為0.68 MPa，與本文均值0.822 MPa的結(jié)果較接近。薛超等[5]對SMILE手術(shù)取出的基質(zhì)透鏡進(jìn)行單軸拉伸實驗，當(dāng)基質(zhì)應(yīng)力為0.02 MPa時，對應(yīng)的彈性模量為1.26MPa± 0.71 MPa，略高于本文計算值；XIANG et al[20]對10例SMILE手術(shù)的角膜基質(zhì)透鏡在水平和豎直兩個方向上進(jìn)行單軸拉伸測試，測得彈性模量分別為1.300 MPa±0.508 MPa、1.142 MPa± 0.238 MPa，也略高于本文計算值；這可能是由于其實驗對象為角膜基質(zhì)部分，而非全角膜材料，且在離體環(huán)境下實驗也會產(chǎn)生一定的誤差。

角膜組織本身為超彈性材料，以往研究中有用HGO模型描述角膜的材料屬性，其中需要待定的參數(shù)有4個，分別為：k1(代表膠原纖維剛度)、k2(代表膠原纖維非線性的無量綱常數(shù))、κ(表示膠原纖維的分散程度，0≤κ<1/3)和C10(表示基質(zhì)相關(guān)材料參數(shù))。材料參數(shù)的確定是研究者根據(jù)實驗擬合出的結(jié)果，但是對于某個實驗結(jié)果，可以擬合出許多種與實驗曲線相符合的參數(shù)組合結(jié)果。天津大學(xué)向堯齊[21]用SMILE手術(shù)取出的角膜透鏡進(jìn)行單軸拉伸實驗，對34組實驗結(jié)果通過最小二乘法擬合出HGO模型中的參數(shù)，C10、k1、k2的均值分別為0.220、0.615、121.633，與ARIZA-GRACIA[15]擬合的結(jié)果(C10=0.05、k1=130.9、k2=2490)相差較大；向堯齊[21]還提出κ值的改變對角膜應(yīng)力-應(yīng)變的影響非常大，也會影響k1和k2的取值。而每個角膜的κ值均不同，這就造成了擬合出的其他參數(shù)結(jié)果可能并不準(zhǔn)確。不同的角膜取材、實驗方法和擬合方法都有可能對擬合參數(shù)的結(jié)果產(chǎn)生非常大的影響，以至于不同研究者的研究結(jié)果相差甚大。目前通過拉伸實驗測定角膜線彈性材料屬性的研究較成熟，在仿真模擬中待定參數(shù)少。故本文采用線彈性本構(gòu)描述角膜材料屬性。此外由于Corvis ST檢測為動態(tài)實驗，角膜材料具有粘彈性特性，粘彈性對實驗結(jié)果存在一定的影響，但由于目前對角膜粘彈性的本構(gòu)模型研究較少，本文未考慮角膜粘彈性對試驗結(jié)果的影響，是本研究的一個局限。

有限元模型是理想化的幾何和力學(xué)模型，與人體真實情況有一定的差別。首先，在角膜幾何模型方面，因很難獲取全部鞏膜和眼球內(nèi)容物的準(zhǔn)確幾何形狀及位置信息，建模時僅考慮了角膜組織，忽略了眼內(nèi)容物等其他組織成分對邊界條件以及模擬結(jié)果產(chǎn)生的影響；角膜中包含的5層組織為各向異性材料，且力學(xué)性能各不相同，根據(jù)DIAS et al[4]的研究，角膜前、后基質(zhì)的彈性模量并不相同，本模型未對角膜各部分分別賦予材料屬性；不同幾何形態(tài)的角膜表面在氣流載荷作用時壓力并不完全相同，本文假定了所有角膜受到的氣流壓力值均相同。其次，同一患者在進(jìn)行多次檢測時，眼壓、角膜頂點位移等檢測結(jié)果也有一定程度差別，這也可能造成模擬結(jié)果與真實情況之間的誤差。最后，本文選取的20例近視患者樣本數(shù)量有限，并不能代表全部近視患者人群，且考慮到其他個體差異原因，本研究所求角膜彈性模量值并不具有普適性，僅可作為參考值。

根據(jù)本文所求近視眼角膜彈性模量，可為建立近視角膜有限元模型提供一定的參考，以期為后續(xù)更準(zhǔn)確地模擬角膜屈光手術(shù)及預(yù)測術(shù)后角膜力學(xué)性能提供理論基礎(chǔ)。