研究在體角膜材料的力學特性，是進一步借助有限元方法研究角膜生理與病理現象的重要基礎。本文借助可視化角膜生物力學分析儀所提供的數據，使用標準線性固體模型計算正常角膜和圓錐角膜在脈沖氣壓作用下的彈性系數（E）和粘性系數（η）。研究發現，正常角膜和圓錐角膜的 E 與 η 之間的差異具有統計學意義（P < 0.05）。受試者工作特征曲線分析發現，E、η 及其聯合指標的曲線下面積（AUC）分別為 0.776、0.895、0.948，說明該指標可以對圓錐角膜進行較好的預測。本文研究結果或可為圓錐角膜的早期診斷提供參考，以避免術后圓錐角膜的發生，因此具有一定的臨床實用價值。
引用本文： 趙科超, 王曉君, 陳維毅, 賀瑞, 李曉娜, 高妍. 正常角膜與圓錐角膜粘彈性的對比研究. 生物醫學工程學雜志, 2019, 36(4): 613-618. doi: 10.7507/1001-5515.201812052 復制
引言
角膜是人體眼部組織的重要組成部分，它在人眼屈光系統中起著十分重要的作用。據統計，在我國眼疾患者致盲原因中，角膜病位居第二位，是致盲者中最常見的病因之一[1]。圓錐角膜是一種原發性角膜變性疾病，主要表現為患者的角膜中央區變薄呈現前凸圓錐狀，并伴隨著不規則散光和高度近視[2]。角膜基質中膠原纖維和纖維間物質的丟失或滑移可導致角膜生物力學不穩定，這也是誘發圓錐角膜的一個因素[3-4]。針對在體角膜材料的力學特性開展的細致研究，可為圓錐角膜的早期診斷提供幫助，同時也可在此基礎上借助有限元方法設計和評估角膜手術，具有重要的臨床實用價值和意義。
角膜具有復雜的生物力學特性，是典型的粘彈性組織[5]。針對角膜的研究主要有離體實驗和在體研究兩種情況。離體實驗主要包括角膜軸向拉伸實驗、角膜膨脹實驗和離體全眼球測量[6]。由于角膜具有其非線性、各向異性、粘彈性等生物力學性能，其彈性模量在不同的生理條件和實驗條件下的測量值差異很大[7]。不同的加載速率下，角膜也表現出不同的彈性模量[8-11]。在體研究則主要是通過臨床檢測儀器直接對在體角膜進行測量，利用儀器所提供的參數和數據進一步分析在體角膜的生物力學性能。在體測量更好地保留了角膜的生理環境和結構特點，測量結果能夠更加真實地反映角膜力學性能。眼反應分析儀（ocular response analyze，ORA）和可視化角膜生物力學分析儀（corneal visualization Scheimpflug technology，Corvis ST）是目前臨床主要采用的角膜測量儀器，兩者都是對角膜施加脈沖氣流，角膜凹陷之后回復至初始狀態，在此過程中記錄角膜變形等參數，進而分析角膜的生物力學特性。
Glass 等[8]使用三參量標準線性固體模型分析了粘彈性系數的改變對角膜 ORA 檢測指標的影響。Tian 等[12]通過使用 Corvis ST 對比分析了圓錐角膜和正常角膜的參數差異，發現變形幅度能夠較好地預測圓錐角膜。Wang 等[13]借助 Corvis ST 的檢測數據，提出了瞬時剛度系數和能量吸收面積兩個參數，并發現二者能較好地預測圓錐角膜。但上述研究工作均針對角膜的整體力學性能，沒有排除角膜幾何參數的影響，也未直接給出角膜材料的粘彈性參數。
鑒于此，本文借助 Corvis ST 和三維眼前節分析儀所提供的脈沖壓力、角膜弧長變化（delta arclength，Darc）、中央角膜厚度（central corneal thickness，CCT）、角膜曲率半徑（corneal radius，CR）等數據，求解角膜在脈沖氣壓作用下的應力、應變隨時間的變化關系，使用標準線性固體模型計算角膜的彈性系數（E）和粘性系數（η）。此外，本文還分析了正常角膜與圓錐角膜的彈性系數和粘性系數的差異，以期通過對比粘彈性系數之間的差異實現對圓錐角膜的臨床早期診斷，避免術后圓錐角膜的發生。
1 方法
從山西省眼科醫院采集正常角膜受試者和圓錐角膜患者的 Corvis ST 與三維眼前節分析儀臨床檢測數據，包括 CCT 值、CR 值、Darc 值。采集數據包括正常角膜 31 例、圓錐角膜 24 例，不包含患者的任何其他個人信息。信息采集工作由山西省眼科醫院準分子激光科相關人員協助完成，該過程由山西省眼科醫院倫理委員會審核批準通過。
首先，將作用在角膜表面的脈沖壓力轉化為脈沖壓強，進而求得角膜的周向應力。借助 Corvis ST 測得的 Darc 值求得角膜變形過程中的應變。然后，選用標準線性固體模型作為角膜本構模型，利用多參數擬合求解角膜的 E 與 η。
1.1 儀器設備
本文在研究過程中使用了如下的測量儀器和計算分析軟件：
（1）Corvis ST 分析儀（G/72100，Oculus，德國）；
（2）三維眼前節分析儀（Pentacam，E/70700，Oculus，德國）；
（3）綜合優化軟件包 1stOpt（5.0，7D-Soft High Technology Inc.，中國）；
（4）科學計算軟件 Matlab（2014a，MathWorks，美國）；
（5）統計產品與服務解決方案軟件 SPSS（17.0，SPSS Inc.，美國）
1.2 脈沖氣流壓強
Corvis ST 在檢測過程可以記錄脈沖氣流壓強數據，但由于射流氣體在空氣中存在一定的能量耗散，該數據無法表示作用在角膜表面的實際壓強，需要對脈沖壓強進行重新計算。
采用 Wang 等[13]實驗測量后給出的擬合公式計算作用在角膜上的脈沖壓力，如式（1）所示：
 |
其中，r 表示脈沖氣流作用半徑，t 表示脈沖氣流作用時間，a1 = 36.04、a2 = 9.799、b1 = 15.54、b2 = 9.642、c1 = 5.119、c2 = 3.777。
脈沖氣壓的時間空間分布之間是相互獨立的[14]，即 P（r, t）= P（r）·P（t）。Joda 等[15]通過實驗測量給出 Corvis ST 脈沖氣壓在角膜表面的空間分布，該分布符合高斯分布。本文據此進行高斯擬合得到脈沖氣壓的空間分布函數，如式（2）所示：
 |
將 P（r, t）= P（r）·P（t）代入式（1）和式（2），可求得作用在角膜表面中央處的脈沖壓強隨時間變化的計算公式，如式（3）所示：
 |
Corvis ST 測量的 Darc 值反映了直徑 7 mm 范圍內的角膜弧長變化，如圖1 所示。圖中紅線表示角膜初始狀態時的弧長，藍線表示角膜在最大凹陷時的弧長，后者與前者的差值即為該時刻角膜弧長變化。本文對應地選擇直徑 7 mm 的范圍對角膜的受力情況進行分析。假設脈沖氣壓作用到角膜表面時是均勻分布的，則平均脈沖氣壓
的計算方法如式（4）所示：
圖1
角膜最大凹陷時的弧長變化示意圖
Figure1.
Diagram of dilated arclength of cornea in deepest concave
 |
1.3 角膜應力應變
本文假定角膜為厚度均勻、各向同性的球形膜，所受脈沖氣壓
和眼內壓分別均勻分布在角膜的前后表面，且角膜變形關于脈沖氣流中軸線對稱。由于角膜厚度很小，因此，當考慮角膜整體弧長變化時可以忽略彎曲應力對變形的影響，而主要研究作用在板層平面的拉壓應力[8]。角膜縱向切面的受力示意圖如圖2 所示，其中板層平面的周向應力為 σ，脈沖氣壓作用范圍為 2r1，眼內壓作用范圍為 2r2。
圖2
角膜在脈沖氣流作用下所受壓力示意圖
Figure2.
Diagram of pressure on cornea under air pulse
將角膜僅在眼內壓作用下的狀態作為初始狀態，此時 CR 以符號 R 表示，CCT 以符號 T 表示。根據玻爾茲曼疊加原理，可以對脈沖氣壓的作用過程進行單獨分析。當給角膜施加如圖2 所示脈沖氣壓
時，可計算角膜在脈沖氣流作用下的應力，如式（5）所示：
 |
假設角膜為線性粘彈性材料，角膜在周向應力作用下的應變為 ε，計算方法如式（6）所示：
 |
其中，J 為蠕變模量，ν 表示角膜材料泊松比，取 ν = 0.49。
在圖2 中取微元 dθ，則微元 dθ 所對應圓弧的弧長變化量（以符號 Δ 表示）的計算方法如式（7）所示：
 |
根據式（7）可以求解 2r2 范圍內角膜整體應變，如式（8）、（9）所示。
 |
 |
脈沖氣壓作用過程中，直徑 7 mm 范圍內的角膜整體應變的計算方法如式（10）所示。其中，脈沖氣壓為零時直徑 7 mm 范圍內所對應的角膜初始弧長用 arc 表示，Corvis ST 測量的 Darc 值用 darc 表示。
 |
1.4 角膜本構模型
本文選擇標準線性固體模型來描述角膜的粘彈性力學行為，粘彈性系數用 E1、E2、η 表示，如圖3 所示。其本構方程用式（11）表示，σ 表示平均應力、ε 表示整體應變。
圖3
標準線性固體模型
Figure3.
Standard linear solid model
 |
其中，
。
正常狀態下角膜后表面受眼內壓作用，角膜板層產生拉伸應力，角膜伸長；Corvis ST 檢測時，脈沖氣壓作用在角膜前表面，使得角膜板層拉伸應力減小，角膜縮短。Corvis ST 檢測過程中，作用在角膜前表面的脈沖氣壓先增大后減小，中央處脈沖氣壓最大值（38 mm Hg）高于眼內壓，角膜在脈沖氣流的作用下局部被壓平（壓平長度為 2 mm 左右），然后凹陷，之后反彈經第二次壓平后恢復至初始狀態。Corvis ST 檢測的是直徑 7 mm 范圍內的 Darc 值，該段角膜前表面脈沖氣壓的平均值最大為 13 mm Hg，小于眼內壓的平均值，因此該段角膜在脈沖氣壓的作用下會經歷兩個狀態：整體縮短至最大凹陷狀態，之后伸長恢復至初始狀態。
1.5 參數求解
由式（9）計算角膜應力隨時間變化的函數關系，進而求得應力率；對 Corvis ST 所測角膜應變進行數據擬合，獲得角膜應變隨時間變化的函數關系，從而求得角膜應變率。將應力、應變、應力率、應變率代入式（11），通過多參數擬合求解粘彈性系數 E1、E2、η。
由于本構模型的局限性，擬合結果無法與角膜實際情況完全匹配，在多參數擬合過程中會出現多個局部最優解的情況。為了使多參數擬合得到更加穩定的解，本文參照 Glass 等[8]將三參數擬合改為雙參數擬合，即假設 E1 = E2 = E。利用綜合優化軟件包 1stOpt 中的多參數擬合對角膜的 E 和 η 進行求解。
1.6 數據處理及統計學方法
Corvis ST 檢測時施加在眼角膜上的脈沖氣流具有高度的一致性[16]，脈沖壓力選用式（1）進行計算。在 Corvis ST 輸出參數的重復性研究中，Hon 等[17]發現 CCT 值的重復性最好，CCT 值及 Darc 值采用 Corvis ST 所提供的數據。三維眼前節分析儀 Pentacam 所測得的角膜曲率重復性較好[18]，本文選取 Pentacam 測得的水平和垂直曲率半徑的平均值作為 CR 值。
使用綜合優化軟件包 1stOpt 進行多參數擬合，算法采用麥考特法，模式為標準+通用全局優化法。角膜在脈沖氣流作用下的應力、應變及其所對應的導數隨時間的變化關系通過 Matlab 軟件進行編程計算。
為了比較正常角膜和圓錐角膜粘彈性之間的差異，利用 SPSS 軟件對正常角膜和圓錐角膜的相關數據進行統計學分析。使用柯爾莫可洛夫-斯米洛夫檢驗 (Kolmogorov-Smirnov，K-S test) 法對所測數據的正態分布情況進行檢驗，數據用“均值 ± 標準差”表示。采用獨立樣本 t 檢驗對正常角膜和圓錐角膜組進行分析，P < 0.05 則認為差異有統計學意義。采用 SPSS 軟件中的受試者工作特征曲線（receiver operating characteristic curve，ROC）對 E、η 及其聯合指標的預測精度進行綜合評價。
2 結果
對 31 例正常角膜、24 例圓錐角膜數據分別進行多參數擬合，擬合結果的決定系數均大于 0.984，擬合效果較好。
正常角膜和圓錐角膜的 CCT 值、CR 值、E 和 η 的分布均屬于正態分布（P > 0.05）。正常角膜和圓錐角膜的 CCT 值、CR 值、E 和 η 的“均值 ± 標準差”如表1 所示。正常角膜和圓錐角膜的 CCT 值、CR 值、E 和 η 之間差異均具有統計學意義（P < 0.05）。
表1
正常角膜和圓錐角膜 CCT 值、CR 值、E 和 η 的對比
Table1.
Comparison of CCT, CR, E and η in healthy group and keratogonus group
根據 ROC 曲線分析結果，繪制了 E、η 及其聯合指標的 ROC 曲線，如圖4 所示。三者的曲線下面積（area under curve，AUC）分別為 0.776、0.895、0.948，表明該方法對圓錐角膜的預測達到了一個較好的水平。E 的截斷值為 1.51 MPa，敏感性為 0.75，特異性為 0.874；η 的截斷值為 1.29 MPa·s，敏感性為 0.917，特異性為 0.874；聯合指標的截斷值為 E = 1.44 MPa，η = 1.425 MPa·s，敏感性為 0.833，特異性為 0.968。
圖4
E、η 及其聯合指標的 ROC 曲線
Figure4.
ROC curves of E, η and their combined indicators
3 討論
本文提出了一種求解在體角膜 E 和 η 的計算方法，選用標準線性固體模型作為角膜的本構模型，將 Corvis ST 提供的脈沖壓力數據和弧長變化數據轉化為角膜固定弧段的應力曲線與應變曲線，通過多參數擬合求解在體角膜 E 和 η。現有研究一般是通過對比頂點位移、變形幅度以及角膜整體剛度等的差異來區分圓錐角膜與正常角膜，而本文求解了在體角膜材料的 E 和 η，去除了角膜厚度等幾何因素的影響，對圓錐角膜與正常角膜的材料力學性能進行了比較。
本文研究發現，采用標準線性固體模型進行多參數擬合，擬合的決定系數均大于 0.984，說明該模型能夠較好地反映角膜的粘彈性行為。擬合結果表明圓錐角膜 E 和 η 都大于正常角膜，ROC 曲線分析發現 E、η 及其聯合指標對于判斷圓錐角膜具有較高的敏感性和特異性，能夠為圓錐角膜的臨床診斷提供一定的參考。
本文發現，圓錐角膜 E 大于正常角膜，這一變化可能與圓錐角膜微觀結構的改變有關。White 等[19]研究發現，正常角膜在基質層和后彈力層之間存在網狀的微纖維束，而圓錐角膜的微纖維束則主要出現在上皮層下方。Muller 等[20]發現由于角膜曲率的改變，圓錐角膜基質中基質神經密度有所增加。這些研究均表明，圓錐角膜和正常角膜在微觀結構上存在差異，這一差異可能是引起圓錐角膜 E 增加的原因。本文通過簡化計算發現，圓錐角膜 E 雖有所增加，但由于角膜厚度的降低，其整體剛度卻有所降低，這與 Wang 等[13]發現的圓錐角膜的整體剛度減小相一致。因此本文將其簡化計算，即當考慮寬度相等的矩形截面時，圓錐角膜與正常角膜剛度之比，等于彈性模量之比與厚度三次方之比的乘積。代入圓錐角膜和正常角膜彈性模量與厚度，計算可得圓錐角膜與正常角膜剛度之比為 0.84。這一結果表明圓錐角膜彈性模量雖大于正常角膜，但其剛度卻小于正常角膜。
本文在研究過程中也存在一些局限性。雖然角膜厚度在中央區域變化較小，但仍舊有所差別，尤其對于圓錐角膜，等厚度假設可能使得計算結果準確度有所降低。即使 Corvis ST 提供的脈沖壓力高度重復，但仍不能確保作用在每個角膜上的脈沖壓力完全一致。本文對脈沖氣壓在角膜表面均勻分布的假設與實際情況存在一定差異。同時本文選擇的本構方程是線性的，但角膜具有非線性性質，這可能也使得計算結果存在一定誤差。其中三參數變為雙參數的處理方式也可能引起一定的誤差。此外，由于數據量有限，本文沒有考慮年齡和性別這兩種因素對計算結果的影響。
綜上所述，本文提出了一種求解在體角膜 E 和 η 的計算方法，并發現圓錐角膜與正常角膜的 E 和 η 之間的差異具有統計學意義。這一發現有利于實現對圓錐角膜的早期診斷，具有重要的臨床實用價值。同時，本方法計算所得的 E 和 η 或可用于角膜有限元模型的建立，從而為術后角膜變形預測、角膜個性化治療等臨床研究提供支持。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
角膜是人體眼部組織的重要組成部分，它在人眼屈光系統中起著十分重要的作用。據統計，在我國眼疾患者致盲原因中，角膜病位居第二位，是致盲者中最常見的病因之一[1]。圓錐角膜是一種原發性角膜變性疾病，主要表現為患者的角膜中央區變薄呈現前凸圓錐狀，并伴隨著不規則散光和高度近視[2]。角膜基質中膠原纖維和纖維間物質的丟失或滑移可導致角膜生物力學不穩定，這也是誘發圓錐角膜的一個因素[3-4]。針對在體角膜材料的力學特性開展的細致研究，可為圓錐角膜的早期診斷提供幫助，同時也可在此基礎上借助有限元方法設計和評估角膜手術，具有重要的臨床實用價值和意義。
角膜具有復雜的生物力學特性，是典型的粘彈性組織[5]。針對角膜的研究主要有離體實驗和在體研究兩種情況。離體實驗主要包括角膜軸向拉伸實驗、角膜膨脹實驗和離體全眼球測量[6]。由于角膜具有其非線性、各向異性、粘彈性等生物力學性能，其彈性模量在不同的生理條件和實驗條件下的測量值差異很大[7]。不同的加載速率下，角膜也表現出不同的彈性模量[8-11]。在體研究則主要是通過臨床檢測儀器直接對在體角膜進行測量，利用儀器所提供的參數和數據進一步分析在體角膜的生物力學性能。在體測量更好地保留了角膜的生理環境和結構特點，測量結果能夠更加真實地反映角膜力學性能。眼反應分析儀（ocular response analyze，ORA）和可視化角膜生物力學分析儀（corneal visualization Scheimpflug technology，Corvis ST）是目前臨床主要采用的角膜測量儀器，兩者都是對角膜施加脈沖氣流，角膜凹陷之后回復至初始狀態，在此過程中記錄角膜變形等參數，進而分析角膜的生物力學特性。
Glass 等[8]使用三參量標準線性固體模型分析了粘彈性系數的改變對角膜 ORA 檢測指標的影響。Tian 等[12]通過使用 Corvis ST 對比分析了圓錐角膜和正常角膜的參數差異，發現變形幅度能夠較好地預測圓錐角膜。Wang 等[13]借助 Corvis ST 的檢測數據，提出了瞬時剛度系數和能量吸收面積兩個參數，并發現二者能較好地預測圓錐角膜。但上述研究工作均針對角膜的整體力學性能，沒有排除角膜幾何參數的影響，也未直接給出角膜材料的粘彈性參數。
鑒于此，本文借助 Corvis ST 和三維眼前節分析儀所提供的脈沖壓力、角膜弧長變化（delta arclength，Darc）、中央角膜厚度（central corneal thickness，CCT）、角膜曲率半徑（corneal radius，CR）等數據，求解角膜在脈沖氣壓作用下的應力、應變隨時間的變化關系，使用標準線性固體模型計算角膜的彈性系數（E）和粘性系數（η）。此外，本文還分析了正常角膜與圓錐角膜的彈性系數和粘性系數的差異，以期通過對比粘彈性系數之間的差異實現對圓錐角膜的臨床早期診斷，避免術后圓錐角膜的發生。
1 方法
從山西省眼科醫院采集正常角膜受試者和圓錐角膜患者的 Corvis ST 與三維眼前節分析儀臨床檢測數據，包括 CCT 值、CR 值、Darc 值。采集數據包括正常角膜 31 例、圓錐角膜 24 例，不包含患者的任何其他個人信息。信息采集工作由山西省眼科醫院準分子激光科相關人員協助完成，該過程由山西省眼科醫院倫理委員會審核批準通過。
首先，將作用在角膜表面的脈沖壓力轉化為脈沖壓強，進而求得角膜的周向應力。借助 Corvis ST 測得的 Darc 值求得角膜變形過程中的應變。然后，選用標準線性固體模型作為角膜本構模型，利用多參數擬合求解角膜的 E 與 η。
1.1 儀器設備
本文在研究過程中使用了如下的測量儀器和計算分析軟件：
（1）Corvis ST 分析儀（G/72100，Oculus，德國）；
（2）三維眼前節分析儀（Pentacam，E/70700，Oculus，德國）；
（3）綜合優化軟件包 1stOpt（5.0，7D-Soft High Technology Inc.，中國）；
（4）科學計算軟件 Matlab（2014a，MathWorks，美國）；
（5）統計產品與服務解決方案軟件 SPSS（17.0，SPSS Inc.，美國）
1.2 脈沖氣流壓強
Corvis ST 在檢測過程可以記錄脈沖氣流壓強數據，但由于射流氣體在空氣中存在一定的能量耗散，該數據無法表示作用在角膜表面的實際壓強，需要對脈沖壓強進行重新計算。
采用 Wang 等[13]實驗測量后給出的擬合公式計算作用在角膜上的脈沖壓力，如式（1）所示：
|
其中，r 表示脈沖氣流作用半徑，t 表示脈沖氣流作用時間，a1 = 36.04、a2 = 9.799、b1 = 15.54、b2 = 9.642、c1 = 5.119、c2 = 3.777。
脈沖氣壓的時間空間分布之間是相互獨立的[14]，即 P（r, t）= P（r）·P（t）。Joda 等[15]通過實驗測量給出 Corvis ST 脈沖氣壓在角膜表面的空間分布，該分布符合高斯分布。本文據此進行高斯擬合得到脈沖氣壓的空間分布函數，如式（2）所示：
|
將 P（r, t）= P（r）·P（t）代入式（1）和式（2），可求得作用在角膜表面中央處的脈沖壓強隨時間變化的計算公式，如式（3）所示：
|
Corvis ST 測量的 Darc 值反映了直徑 7 mm 范圍內的角膜弧長變化，如圖1 所示。圖中紅線表示角膜初始狀態時的弧長，藍線表示角膜在最大凹陷時的弧長，后者與前者的差值即為該時刻角膜弧長變化。本文對應地選擇直徑 7 mm 的范圍對角膜的受力情況進行分析。假設脈沖氣壓作用到角膜表面時是均勻分布的，則平均脈沖氣壓的計算方法如式（4）所示：
圖1
角膜最大凹陷時的弧長變化示意圖
Figure1.
Diagram of dilated arclength of cornea in deepest concave
|
1.3 角膜應力應變
本文假定角膜為厚度均勻、各向同性的球形膜，所受脈沖氣壓和眼內壓分別均勻分布在角膜的前后表面，且角膜變形關于脈沖氣流中軸線對稱。由于角膜厚度很小，因此，當考慮角膜整體弧長變化時可以忽略彎曲應力對變形的影響，而主要研究作用在板層平面的拉壓應力[8]。角膜縱向切面的受力示意圖如圖2 所示，其中板層平面的周向應力為 σ，脈沖氣壓作用范圍為 2r1，眼內壓作用范圍為 2r2。
圖2
角膜在脈沖氣流作用下所受壓力示意圖
Figure2.
Diagram of pressure on cornea under air pulse
將角膜僅在眼內壓作用下的狀態作為初始狀態，此時 CR 以符號 R 表示，CCT 以符號 T 表示。根據玻爾茲曼疊加原理，可以對脈沖氣壓的作用過程進行單獨分析。當給角膜施加如圖2 所示脈沖氣壓時，可計算角膜在脈沖氣流作用下的應力，如式（5）所示：
|
假設角膜為線性粘彈性材料，角膜在周向應力作用下的應變為 ε，計算方法如式（6）所示：
|
其中，J 為蠕變模量，ν 表示角膜材料泊松比，取 ν = 0.49。
在圖2 中取微元 dθ，則微元 dθ 所對應圓弧的弧長變化量（以符號 Δ 表示）的計算方法如式（7）所示：
|
根據式（7）可以求解 2r2 范圍內角膜整體應變，如式（8）、（9）所示。
|
|
脈沖氣壓作用過程中，直徑 7 mm 范圍內的角膜整體應變的計算方法如式（10）所示。其中，脈沖氣壓為零時直徑 7 mm 范圍內所對應的角膜初始弧長用 arc 表示，Corvis ST 測量的 Darc 值用 darc 表示。
|
1.4 角膜本構模型
本文選擇標準線性固體模型來描述角膜的粘彈性力學行為，粘彈性系數用 E1、E2、η 表示，如圖3 所示。其本構方程用式（11）表示，σ 表示平均應力、ε 表示整體應變。
圖3
標準線性固體模型
Figure3.
Standard linear solid model
|
其中，。
正常狀態下角膜后表面受眼內壓作用，角膜板層產生拉伸應力，角膜伸長；Corvis ST 檢測時，脈沖氣壓作用在角膜前表面，使得角膜板層拉伸應力減小，角膜縮短。Corvis ST 檢測過程中，作用在角膜前表面的脈沖氣壓先增大后減小，中央處脈沖氣壓最大值（38 mm Hg）高于眼內壓，角膜在脈沖氣流的作用下局部被壓平（壓平長度為 2 mm 左右），然后凹陷，之后反彈經第二次壓平后恢復至初始狀態。Corvis ST 檢測的是直徑 7 mm 范圍內的 Darc 值，該段角膜前表面脈沖氣壓的平均值最大為 13 mm Hg，小于眼內壓的平均值，因此該段角膜在脈沖氣壓的作用下會經歷兩個狀態：整體縮短至最大凹陷狀態，之后伸長恢復至初始狀態。
1.5 參數求解
由式（9）計算角膜應力隨時間變化的函數關系，進而求得應力率；對 Corvis ST 所測角膜應變進行數據擬合，獲得角膜應變隨時間變化的函數關系，從而求得角膜應變率。將應力、應變、應力率、應變率代入式（11），通過多參數擬合求解粘彈性系數 E1、E2、η。
由于本構模型的局限性，擬合結果無法與角膜實際情況完全匹配，在多參數擬合過程中會出現多個局部最優解的情況。為了使多參數擬合得到更加穩定的解，本文參照 Glass 等[8]將三參數擬合改為雙參數擬合，即假設 E1 = E2 = E。利用綜合優化軟件包 1stOpt 中的多參數擬合對角膜的 E 和 η 進行求解。
1.6 數據處理及統計學方法
Corvis ST 檢測時施加在眼角膜上的脈沖氣流具有高度的一致性[16]，脈沖壓力選用式（1）進行計算。在 Corvis ST 輸出參數的重復性研究中，Hon 等[17]發現 CCT 值的重復性最好，CCT 值及 Darc 值采用 Corvis ST 所提供的數據。三維眼前節分析儀 Pentacam 所測得的角膜曲率重復性較好[18]，本文選取 Pentacam 測得的水平和垂直曲率半徑的平均值作為 CR 值。
使用綜合優化軟件包 1stOpt 進行多參數擬合，算法采用麥考特法，模式為標準+通用全局優化法。角膜在脈沖氣流作用下的應力、應變及其所對應的導數隨時間的變化關系通過 Matlab 軟件進行編程計算。
為了比較正常角膜和圓錐角膜粘彈性之間的差異，利用 SPSS 軟件對正常角膜和圓錐角膜的相關數據進行統計學分析。使用柯爾莫可洛夫-斯米洛夫檢驗 (Kolmogorov-Smirnov，K-S test) 法對所測數據的正態分布情況進行檢驗，數據用“均值 ± 標準差”表示。采用獨立樣本 t 檢驗對正常角膜和圓錐角膜組進行分析，P < 0.05 則認為差異有統計學意義。采用 SPSS 軟件中的受試者工作特征曲線（receiver operating characteristic curve，ROC）對 E、η 及其聯合指標的預測精度進行綜合評價。
2 結果
對 31 例正常角膜、24 例圓錐角膜數據分別進行多參數擬合，擬合結果的決定系數均大于 0.984，擬合效果較好。
正常角膜和圓錐角膜的 CCT 值、CR 值、E 和 η 的分布均屬于正態分布（P > 0.05）。正常角膜和圓錐角膜的 CCT 值、CR 值、E 和 η 的“均值 ± 標準差”如表1 所示。正常角膜和圓錐角膜的 CCT 值、CR 值、E 和 η 之間差異均具有統計學意義（P < 0.05）。
表1
正常角膜和圓錐角膜 CCT 值、CR 值、E 和 η 的對比
Table1.
Comparison of CCT, CR, E and η in healthy group and keratogonus group
根據 ROC 曲線分析結果，繪制了 E、η 及其聯合指標的 ROC 曲線，如圖4 所示。三者的曲線下面積（area under curve，AUC）分別為 0.776、0.895、0.948，表明該方法對圓錐角膜的預測達到了一個較好的水平。E 的截斷值為 1.51 MPa，敏感性為 0.75，特異性為 0.874；η 的截斷值為 1.29 MPa·s，敏感性為 0.917，特異性為 0.874；聯合指標的截斷值為 E = 1.44 MPa，η = 1.425 MPa·s，敏感性為 0.833，特異性為 0.968。
圖4
E、η 及其聯合指標的 ROC 曲線
Figure4.
ROC curves of E, η and their combined indicators
3 討論
本文提出了一種求解在體角膜 E 和 η 的計算方法，選用標準線性固體模型作為角膜的本構模型，將 Corvis ST 提供的脈沖壓力數據和弧長變化數據轉化為角膜固定弧段的應力曲線與應變曲線，通過多參數擬合求解在體角膜 E 和 η。現有研究一般是通過對比頂點位移、變形幅度以及角膜整體剛度等的差異來區分圓錐角膜與正常角膜，而本文求解了在體角膜材料的 E 和 η，去除了角膜厚度等幾何因素的影響，對圓錐角膜與正常角膜的材料力學性能進行了比較。
本文研究發現，采用標準線性固體模型進行多參數擬合，擬合的決定系數均大于 0.984，說明該模型能夠較好地反映角膜的粘彈性行為。擬合結果表明圓錐角膜 E 和 η 都大于正常角膜，ROC 曲線分析發現 E、η 及其聯合指標對于判斷圓錐角膜具有較高的敏感性和特異性，能夠為圓錐角膜的臨床診斷提供一定的參考。
本文發現，圓錐角膜 E 大于正常角膜，這一變化可能與圓錐角膜微觀結構的改變有關。White 等[19]研究發現，正常角膜在基質層和后彈力層之間存在網狀的微纖維束，而圓錐角膜的微纖維束則主要出現在上皮層下方。Muller 等[20]發現由于角膜曲率的改變，圓錐角膜基質中基質神經密度有所增加。這些研究均表明，圓錐角膜和正常角膜在微觀結構上存在差異，這一差異可能是引起圓錐角膜 E 增加的原因。本文通過簡化計算發現，圓錐角膜 E 雖有所增加，但由于角膜厚度的降低，其整體剛度卻有所降低，這與 Wang 等[13]發現的圓錐角膜的整體剛度減小相一致。因此本文將其簡化計算，即當考慮寬度相等的矩形截面時，圓錐角膜與正常角膜剛度之比，等于彈性模量之比與厚度三次方之比的乘積。代入圓錐角膜和正常角膜彈性模量與厚度，計算可得圓錐角膜與正常角膜剛度之比為 0.84。這一結果表明圓錐角膜彈性模量雖大于正常角膜，但其剛度卻小于正常角膜。
本文在研究過程中也存在一些局限性。雖然角膜厚度在中央區域變化較小，但仍舊有所差別，尤其對于圓錐角膜，等厚度假設可能使得計算結果準確度有所降低。即使 Corvis ST 提供的脈沖壓力高度重復，但仍不能確保作用在每個角膜上的脈沖壓力完全一致。本文對脈沖氣壓在角膜表面均勻分布的假設與實際情況存在一定差異。同時本文選擇的本構方程是線性的，但角膜具有非線性性質，這可能也使得計算結果存在一定誤差。其中三參數變為雙參數的處理方式也可能引起一定的誤差。此外，由于數據量有限，本文沒有考慮年齡和性別這兩種因素對計算結果的影響。
綜上所述，本文提出了一種求解在體角膜 E 和 η 的計算方法，并發現圓錐角膜與正常角膜的 E 和 η 之間的差異具有統計學意義。這一發現有利于實現對圓錐角膜的早期診斷，具有重要的臨床實用價值。同時，本方法計算所得的 E 和 η 或可用于角膜有限元模型的建立，從而為術后角膜變形預測、角膜個性化治療等臨床研究提供支持。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
