基于Pentacam和Corvis ST測量近視患者角膜的彈性模量

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

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

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

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

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

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

1 方法

1.1 數據采集

經山西省眼科醫(yī)院倫理委員會批準通過，選取山西省眼科醫(yī)院準分子激光科確診的近視患者20例。選取的20例受試者中，男12例(12眼)，女8例(8眼)，受試者年齡均值為25.02±4.60歲，均為大專以上受教育程度；健康狀況良好，無慢性病、免疫性疾病、精神疾病及其他疾??；無眼部手術史或外傷史，無圓錐角膜、青光眼等其他眼部疾?。磺舛确€(wěn)定兩年及以上，屈光度年增長不超過0.5 D；雙眼柱鏡度數不超過-4.0 D，等效球鏡度數不高于-10.0 D；角膜厚度大于470 μm，眼壓小于3.192 kPa；停戴各類接觸鏡和塑形鏡一月以上。所有患者經三維眼前節(jié)分析儀和可視化角膜生物力學分析儀進行檢測，檢查均由同一專業(yè)人員操作，每眼檢測3次，取圖像最佳、重復性最好的結果分析。從眼前節(jié)分析儀中獲取患者角膜前后表面高度圖、厚度圖等地形數據；從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檢測的是角膜在眼內壓作用下的坐標位置，并非角膜在無眼壓狀態(tài)下的位置，因此需要通過迭代的方法找到角膜在無眼壓狀態(tài)下的坐標點，進而建立無應力狀態(tài)下的幾何模型，參考ARIZA-GRACIA et al[15]的方法進行計算。首先對角膜上各點所有Z坐標值減去一個微小量δ，δ初始值設置為0.005，在其基礎上施加眼內壓并進行運算，通過Abaqus運算后的*inp文件導出角膜表面各點的Z坐標值并求和，求和結果與儀器檢測值相減，當差值小于10-3時計算結束，通過不斷嘗試找到角膜零應力狀態(tài)下的坐標，本文采集的患者數據較少，故所有模型均采用嘗試法計算角膜的初始狀態(tài)坐標值。

圖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所示。模型網格采用四面體單元，考慮到計算規(guī)模及收斂性，對網格密度進行了驗證。

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

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

1.4 統(tǒng)計學分析

使用SPSS 20.0對結果進行統(tǒng)計學分析，描述性統(tǒng)計結果用均值±標準差進行描述。

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

3 討論

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

不同方法、取材和檢測手段得出的角膜彈性模量值均有差別。黃來鑫等[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，也略高于本文計算值；這可能是由于其實驗對象為角膜基質部分，而非全角膜材料，且在離體環(huán)境下實驗也會產生一定的誤差。

角膜組織本身為超彈性材料，以往研究中有用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檢測為動態(tài)實驗，角膜材料具有粘彈性特性，粘彈性對實驗結果存在一定的影響，但由于目前對角膜粘彈性的本構模型研究較少，本文未考慮角膜粘彈性對試驗結果的影響，是本研究的一個局限。

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

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