本文基于CT圖像數據結合圖像處理軟件建立人體腰椎L4~5節段有限元模型,進行前屈、后伸、側彎和軸向旋轉工況下的生物力學特性分析。首先,選擇一名正常志愿者腰椎L4~5節段的CT圖像數據建立三維模型,模型包括椎體、椎間盤、軟骨終板、連接韌帶和小關節突關節等部分,在前屈、后伸、側彎和軸向旋轉工況加載下進行生物力學性能的分析。結果表明,建立的腰椎模型幾何形態逼真且驗證有效,在500、1 000、1 500和2 000 N分布壓力下,腰椎L4~5有限元模型的軸向位移分別為0.23、0.47、0.76、1.02 mm,與前人在相同條件下離體實驗和有限元分析的研究相近,幾種情況下腰椎整體及椎間盤的應力分布符合腰椎的生物力學特性。因此,本文所建立的腰椎L4~5節段有限元模型能有效地模擬腰椎的生物力學特性,為后續腰椎植入物的生物力學性能研究奠定良好的基礎。
引用本文： 顏文濤, 趙改平, 方新果, 郭昊翔, 馬童, 凃意輝. 人體腰椎L4~5節段有限元建模及分析. 生物醫學工程學雜志, 2014, 31(3): 612-618. doi: 10.7507/1001-5515.20140115 復制
引言
人體腰椎活動節段是腰椎最小的功能單位，呈現出腰椎一般的生物力學特性，全面地認識腰椎活動節段在生理和不同病理條件下的受力狀況，一直是臨床工作者和理論研究者關注的熱點問題之一。自1974年Belytschko等[1]首次進行椎間盤及毗鄰椎骨的有限元分析以來，利用有限元建模的方法已廣泛應用于腰椎生物力學的研究。隨著計算機技術和圖像處理技術的不斷發展，Goel等[2]于1988年首次提出應用CT掃描幾何形狀輸入模型，建立幾何形狀逼真的復雜脊柱有限元模型。在此之后研究人員認為脊柱建模的眾多方法中，利用CT獲得直接的影像數據是最準確的途徑[3]。張德盛等[4]利用CATIA V5軟件提取斷層圖像邊界關鍵點的二維坐標，結合相應層面生成三維坐標，從而建立L3~5腰椎的三維實體模型。蘇晉等[5]在建立全腰椎有限元接觸模型中，利用自編軟件和Hypermesh軟件建立全腰椎有限元模型。傅棟等[6]在腰椎運動節段數字模型建立及三維可視化研究中，利用Mimics軟件對腰椎骨性結構及各種軟組織進行三維重建。在這些已有研究的基礎上，本文基于CT掃描圖像數據，對斷層圖像進行三維逆向重建，幾何模型進行有限單元化，賦予各部分材料屬性，建立腰椎L4~5有限元模型對其進行生物力學特性分析，建模及分析的流程如圖 1所示。
圖1
腰椎L4~5有限元建模及分析流程
Figure1.
Flow of construction and analysis of lumbar L4~5 finite element model
1 材料和方法
1.1 CT數據
選取一名男性志愿者，利用普通X線檢查排除脊柱病變、畸形和損傷等異常情況，采用64排螺旋CT在自然狀態下對其腰椎L4~5節段以層厚1.0 mm進行掃描，得到CT斷層掃描圖像173張，并以512×512像素DICOM3.0標準格式保存。
1.2 幾何模型
基于CT斷層掃描圖像，將圖片導入數字醫學圖像處理軟件Mimics10.01，界定閾值，進行區域增長，提取出椎骨輪廓線后以彩色填充的數據集表示，如圖 2所示。經軟件計算后將二維圖像三維重建，優化直接生成的腰椎L4、L5椎體三角片面網格。
圖2
腰椎L4~5椎體輪廓彩色數據集
Figure2.
Color data set of lumbar L4~5 vertebral body contour
基于腰椎解剖結構和曲率變化的特點，利用逆向工程軟件Geomagic9.0對模型進行填充孔、去除特征及NURBS曲面片擬合等處理，參考腰椎L4~5椎間盤生理參數，提取腰椎L4下表面和腰椎L5上表面，邊緣光順處理后，向各自內法線方向平移0.5 mm(留出終板的位置)，使用填充孔搭橋法創建椎間盤部分，最終擬合椎間盤曲面如圖 3所示。在三維制圖軟件中根據腰椎L4~5椎間盤的生理參數，纖維環占椎間盤的56%，髓核占44%，創建髓核的形狀尺寸，得到腰椎L4~5髓核和纖維環。
圖3
腰椎L4~5椎體及椎間盤實體模型
Figure3.
Entity model of lumbar L4~5 vertebral body and intervertebral disc
1.3 有限元模型的建立
1.3.1 幾何清理
將建立后的L4、L5椎體，L4~5纖維環和髓核模型導入有限元分析前處理軟件Hypermesh10.0進行編輯操作，L4、L5椎體由面形成體，將整個椎體分割為前方椎體與后方骨性單元[5]。
1.3.2 網格劃分
根據腰椎解剖形態，基于“網格包絡”思想進行腰椎L4~5網格劃分，即先劃分椎體表面形成封閉的面網格(L4椎體下表面和L5椎體上表面由軟骨終板面網格提供)，向內偏置1 mm，由此便可通過包絡的面網格分別形成松質骨、皮質骨以及后方骨性體網格。椎間盤則通過共用髓核側面面網格，各自形成纖維環和髓核的面網格包絡后劃分體網格。軟骨終板體網格則通過提取椎間盤上下面網格，向各自外法向偏置0.5 mm而得。同時提取軟骨終板面網格，將L4椎體下終板的上表面網格和L5椎體上終板的下表面網格分別作為L4椎體下表面網格和L5椎體上表面網格。這樣便實現了包括皮質骨、松質骨、纖維環、髓核及與軟骨終板等結構所有節點共節點，如圖 4所示。
圖4
腰椎L4~5椎體及椎間盤有限元模型
Figure4.
Finite element model of lumbar L4~5 vertebral body and intervertebral disc
根據文獻測量起始位置，建立包括前縱韌帶(ALL)、后縱韌帶(PLL)、黃韌帶( LF)、棘間韌帶(ISL)、棘上韌帶(SSL)、橫向韌帶(TL)和囊韌帶(CL)七種韌帶模型[7]。關節軟骨由小關節突接觸關系所定義的面-面接觸表示，接觸面初始間隙約為0.4 mm，摩擦系數為0 [8]。腰椎L4~5有限元模型網格劃分所用單元屬性及單元數目如表 1所示。
表1
腰椎L4~5有限元模型單元屬性
Table1.
Unit attribute of lumbar L4~5 finite element model
1.3.3 材料屬性
骨性結構和椎間盤均假設為各向同性的彈性材料，用彈性模量和泊松比兩個參數描述[9]，韌帶是一種只能承受拉伸載荷而在受壓時響應為零的材料[10]。腰椎L4~5有限元模型的材料特性如表 2所示。
1.3.4 邊界與負載條件
在L4椎體上表面施加500 N分布壓力用于模擬人體上半身重量；同時耦合L4椎體上表面所有節點于中性點，根據右手定則，在該中性點上施加7.5 N·m扭矩[17]，分別模擬前屈、后伸，左、右側彎，左、右軸向旋轉6種加載情況；同時固定L5椎體下表面所有節點，限制其所有的自由度。
2 結果
本文基于CT圖像并充分利用Hypermesh軟件建立了幾何外形真實、組織結構完整、材料屬性可靠的腰椎L4~5有限元模型，共58 590個單元，91 118個節點。將所建立的腰椎L4~5有限元模型導入Ansys12.0中進行計算后處理，分析腰椎L4~5有限元模型在相應邊界條件下椎體和椎間盤等結構的力學行為。
2.1 腰椎L4~5有限元模型的驗證
在500、1 000、 1 500和2 000 N的分布壓力下，腰椎L4~5有限元模型的軸向位移分別為0.23、0.47、0.76和1.02 mm，如圖 5所示，其分布壓力-軸向位移曲線與Brown等[18]、Markolf[19]和Kim等[20]在相同條件下離體實驗和有限元分析的研究結果相近。隨著分布壓力的增加，腰椎L4~5最大位移也增加，幾乎呈線性關系，說明腰椎L4~5在正常受力狀態下表現出彈性性質[5]。
圖5
不同分布壓力下椎體的軸向位移曲線
Figure5.
Axial displacement curve of vertebral body under different distribution pressure
2.2 腰椎L4~5椎體的應力分布
腰椎椎體的應力分析采用Von Mises等效應力作為主要指標，各工況下腰椎整體力學變化如圖 6所示。前屈時最大應力出現在椎體皮質骨，其中L5椎體皮質骨后方承受最大的應力為32.36 MPa，椎弓根和椎管周圍也有較大的應力分布；后伸時椎弓根處存在著明顯的應力集中，小關節承受最大的壓應力為95.08 MPa；左右側彎時最大應力均出現在受壓側椎體皮質骨、椎弓根以及椎管周圍，椎弓根處最大應力為40.07 MPa，同時兩側小關節均受到一定的壓應力；左右軸向旋轉時以椎弓根以及后部一側小關節應力最大，椎體峽部和椎管周圍也存在較大的應力分布，最大應力為48.01 MPa。
圖6
腰椎L4~5有限元模型在各工況下應力云圖
(a)前屈；(b)后伸；(c)左側彎；(d)右側彎；(e)左軸向旋轉；(f)右軸向旋轉
Figure6. Stress nephogram of lumbar L4~5 finite element model under different states(a) flexion； (b) extension； (c) left lateral bending； (d) right lateral bending； (e) left axial rotation； (f) right axial rotation
2.3 腰椎L4~5椎間盤的應力分布
椎間盤在力學性能和功能上相當特殊，其主要生物力學性能為對抗壓縮載荷，各工況下腰椎椎間盤力學變化如圖 7所示。前屈時椎間盤應力主要集中在椎間盤前部且為壓應力，其最大應力值為1.36 MPa；后伸時最大壓縮應力出現在椎間盤后部，其最大值為0.97 MPa；左右側彎時椎間盤受壓側存在明顯的應力集中，其最大值為2.19 MPa，且向椎間盤中心有逐漸減小的趨勢；左右軸向旋轉時纖維環受到扭力而發生傾斜至牽張，最大應力主要出現于軸向旋轉方向側后方，其最大應力值為1.42 MPa。
圖7
腰椎L4~5椎間盤在各工況下應力云圖
(a)前屈；(b)后伸；(c)左側彎；(d)右側彎；(e)左軸向旋轉；(f)右軸向旋轉
Figure7. Stress nephogram of lumbar L4~5 intervertebral disc under different states(a) flexion； (b) extension； (c) left lateral bending； (d) right lateral bending； (e) left axial rotation； (f) right axial rotation
3 討論
本文基于CT圖像數據，利用Mimics軟件無需對圖像進行任何形式的轉化，快速而準確地建立腰椎的三角面片模型，避免了信息的丟失。充分利用Hypermesh軟件結合由三角面片生成的四面體單元，最適用于具有復雜邊界曲面的不規則椎體結構的離散[21]，同時采用十節點保證了計算精度且很大程度上縮短了有限元前處理的時間，提高了建模的效率。基于“網格包絡”思想，巧妙運用共用面網格及網格偏置等網格劃分技巧，實現所有單元共節點，不需要在椎體與椎間盤處進行綁定約束，提高了模型的收斂性，大大提高了計算速度[22]。
腰椎整體的應力云圖反映了椎體在前屈、后伸、側彎和軸向旋轉時的力學分布情況，結果證明椎弓根附近是椎體多種應力集中的部位，后部結構所受應力只有通過椎弓根才能傳達到椎體。這就將有限元分析和臨床研究有效地結合在一起且結果能夠相互印證，再次說明模型的準確性和有效性。同時可以看出，椎體在各工況下受壓縮側應力均較大，在后伸、側彎和軸向旋轉時，可以明顯看到椎體后方小關節突承載著較大的接觸應力，這說明小關節突可以對抗后伸、側彎和軸向旋轉的作用很明顯，這也與實際研究相符[13]。另外，由研究結果可知，椎間盤的應力分布也反映了在各加載情況下椎間盤受壓側應力較為集中，并有向受拉側逐漸擴散消釋的趨勢。
本文有效地模擬了腰椎的生物力學特性，但也存在著一些不足，尚需要完善和改進。首先，CT圖像數據的精度和軟件處理的精度等制約了建模精度。其次，腰椎椎體、椎間盤及韌帶等結構材料的生物力學特性較為復雜，在材料賦值時往往需要簡便處理，因此使得材質特性與真實情況難以維系。最后，腰椎處于人體復雜多變的環境，不僅受到體重的影響，肌肉等軟組織在腰椎運動過程中的牽張作用也是不容忽視的。因此，完善腰椎的有限元模型進行深入的生物力學特性分析，是我們進一步探索的工作。
引言
人體腰椎活動節段是腰椎最小的功能單位，呈現出腰椎一般的生物力學特性，全面地認識腰椎活動節段在生理和不同病理條件下的受力狀況，一直是臨床工作者和理論研究者關注的熱點問題之一。自1974年Belytschko等[1]首次進行椎間盤及毗鄰椎骨的有限元分析以來，利用有限元建模的方法已廣泛應用于腰椎生物力學的研究。隨著計算機技術和圖像處理技術的不斷發展，Goel等[2]于1988年首次提出應用CT掃描幾何形狀輸入模型，建立幾何形狀逼真的復雜脊柱有限元模型。在此之后研究人員認為脊柱建模的眾多方法中，利用CT獲得直接的影像數據是最準確的途徑[3]。張德盛等[4]利用CATIA V5軟件提取斷層圖像邊界關鍵點的二維坐標，結合相應層面生成三維坐標，從而建立L3~5腰椎的三維實體模型。蘇晉等[5]在建立全腰椎有限元接觸模型中，利用自編軟件和Hypermesh軟件建立全腰椎有限元模型。傅棟等[6]在腰椎運動節段數字模型建立及三維可視化研究中，利用Mimics軟件對腰椎骨性結構及各種軟組織進行三維重建。在這些已有研究的基礎上，本文基于CT掃描圖像數據，對斷層圖像進行三維逆向重建，幾何模型進行有限單元化，賦予各部分材料屬性，建立腰椎L4~5有限元模型對其進行生物力學特性分析，建模及分析的流程如圖 1所示。
圖1
腰椎L4~5有限元建模及分析流程
Figure1.
Flow of construction and analysis of lumbar L4~5 finite element model
1 材料和方法
1.1 CT數據
選取一名男性志愿者，利用普通X線檢查排除脊柱病變、畸形和損傷等異常情況，采用64排螺旋CT在自然狀態下對其腰椎L4~5節段以層厚1.0 mm進行掃描，得到CT斷層掃描圖像173張，并以512×512像素DICOM3.0標準格式保存。
1.2 幾何模型
基于CT斷層掃描圖像，將圖片導入數字醫學圖像處理軟件Mimics10.01，界定閾值，進行區域增長，提取出椎骨輪廓線后以彩色填充的數據集表示，如圖 2所示。經軟件計算后將二維圖像三維重建，優化直接生成的腰椎L4、L5椎體三角片面網格。
圖2
腰椎L4~5椎體輪廓彩色數據集
Figure2.
Color data set of lumbar L4~5 vertebral body contour
基于腰椎解剖結構和曲率變化的特點，利用逆向工程軟件Geomagic9.0對模型進行填充孔、去除特征及NURBS曲面片擬合等處理，參考腰椎L4~5椎間盤生理參數，提取腰椎L4下表面和腰椎L5上表面，邊緣光順處理后，向各自內法線方向平移0.5 mm(留出終板的位置)，使用填充孔搭橋法創建椎間盤部分，最終擬合椎間盤曲面如圖 3所示。在三維制圖軟件中根據腰椎L4~5椎間盤的生理參數，纖維環占椎間盤的56%，髓核占44%，創建髓核的形狀尺寸，得到腰椎L4~5髓核和纖維環。
圖3
腰椎L4~5椎體及椎間盤實體模型
Figure3.
Entity model of lumbar L4~5 vertebral body and intervertebral disc
1.3 有限元模型的建立
1.3.1 幾何清理
將建立后的L4、L5椎體，L4~5纖維環和髓核模型導入有限元分析前處理軟件Hypermesh10.0進行編輯操作，L4、L5椎體由面形成體，將整個椎體分割為前方椎體與后方骨性單元[5]。
1.3.2 網格劃分
根據腰椎解剖形態，基于“網格包絡”思想進行腰椎L4~5網格劃分，即先劃分椎體表面形成封閉的面網格(L4椎體下表面和L5椎體上表面由軟骨終板面網格提供)，向內偏置1 mm，由此便可通過包絡的面網格分別形成松質骨、皮質骨以及后方骨性體網格。椎間盤則通過共用髓核側面面網格，各自形成纖維環和髓核的面網格包絡后劃分體網格。軟骨終板體網格則通過提取椎間盤上下面網格，向各自外法向偏置0.5 mm而得。同時提取軟骨終板面網格，將L4椎體下終板的上表面網格和L5椎體上終板的下表面網格分別作為L4椎體下表面網格和L5椎體上表面網格。這樣便實現了包括皮質骨、松質骨、纖維環、髓核及與軟骨終板等結構所有節點共節點，如圖 4所示。
圖4
腰椎L4~5椎體及椎間盤有限元模型
Figure4.
Finite element model of lumbar L4~5 vertebral body and intervertebral disc
根據文獻測量起始位置，建立包括前縱韌帶(ALL)、后縱韌帶(PLL)、黃韌帶( LF)、棘間韌帶(ISL)、棘上韌帶(SSL)、橫向韌帶(TL)和囊韌帶(CL)七種韌帶模型[7]。關節軟骨由小關節突接觸關系所定義的面-面接觸表示，接觸面初始間隙約為0.4 mm，摩擦系數為0 [8]。腰椎L4~5有限元模型網格劃分所用單元屬性及單元數目如表 1所示。
表1
腰椎L4~5有限元模型單元屬性
Table1.
Unit attribute of lumbar L4~5 finite element model
1.3.3 材料屬性
骨性結構和椎間盤均假設為各向同性的彈性材料，用彈性模量和泊松比兩個參數描述[9]，韌帶是一種只能承受拉伸載荷而在受壓時響應為零的材料[10]。腰椎L4~5有限元模型的材料特性如表 2所示。
1.3.4 邊界與負載條件
在L4椎體上表面施加500 N分布壓力用于模擬人體上半身重量；同時耦合L4椎體上表面所有節點于中性點，根據右手定則，在該中性點上施加7.5 N·m扭矩[17]，分別模擬前屈、后伸，左、右側彎，左、右軸向旋轉6種加載情況；同時固定L5椎體下表面所有節點，限制其所有的自由度。
2 結果
本文基于CT圖像并充分利用Hypermesh軟件建立了幾何外形真實、組織結構完整、材料屬性可靠的腰椎L4~5有限元模型，共58 590個單元，91 118個節點。將所建立的腰椎L4~5有限元模型導入Ansys12.0中進行計算后處理，分析腰椎L4~5有限元模型在相應邊界條件下椎體和椎間盤等結構的力學行為。
2.1 腰椎L4~5有限元模型的驗證
在500、1 000、 1 500和2 000 N的分布壓力下，腰椎L4~5有限元模型的軸向位移分別為0.23、0.47、0.76和1.02 mm，如圖 5所示，其分布壓力-軸向位移曲線與Brown等[18]、Markolf[19]和Kim等[20]在相同條件下離體實驗和有限元分析的研究結果相近。隨著分布壓力的增加，腰椎L4~5最大位移也增加，幾乎呈線性關系，說明腰椎L4~5在正常受力狀態下表現出彈性性質[5]。
圖5
不同分布壓力下椎體的軸向位移曲線
Figure5.
Axial displacement curve of vertebral body under different distribution pressure
2.2 腰椎L4~5椎體的應力分布
腰椎椎體的應力分析采用Von Mises等效應力作為主要指標，各工況下腰椎整體力學變化如圖 6所示。前屈時最大應力出現在椎體皮質骨，其中L5椎體皮質骨后方承受最大的應力為32.36 MPa，椎弓根和椎管周圍也有較大的應力分布；后伸時椎弓根處存在著明顯的應力集中，小關節承受最大的壓應力為95.08 MPa；左右側彎時最大應力均出現在受壓側椎體皮質骨、椎弓根以及椎管周圍，椎弓根處最大應力為40.07 MPa，同時兩側小關節均受到一定的壓應力；左右軸向旋轉時以椎弓根以及后部一側小關節應力最大，椎體峽部和椎管周圍也存在較大的應力分布，最大應力為48.01 MPa。
圖6
腰椎L4~5有限元模型在各工況下應力云圖
(a)前屈；(b)后伸；(c)左側彎；(d)右側彎；(e)左軸向旋轉；(f)右軸向旋轉
Figure6. Stress nephogram of lumbar L4~5 finite element model under different states(a) flexion； (b) extension； (c) left lateral bending； (d) right lateral bending； (e) left axial rotation； (f) right axial rotation
2.3 腰椎L4~5椎間盤的應力分布
椎間盤在力學性能和功能上相當特殊，其主要生物力學性能為對抗壓縮載荷，各工況下腰椎椎間盤力學變化如圖 7所示。前屈時椎間盤應力主要集中在椎間盤前部且為壓應力，其最大應力值為1.36 MPa；后伸時最大壓縮應力出現在椎間盤后部，其最大值為0.97 MPa；左右側彎時椎間盤受壓側存在明顯的應力集中，其最大值為2.19 MPa，且向椎間盤中心有逐漸減小的趨勢；左右軸向旋轉時纖維環受到扭力而發生傾斜至牽張，最大應力主要出現于軸向旋轉方向側后方，其最大應力值為1.42 MPa。
圖7
腰椎L4~5椎間盤在各工況下應力云圖
(a)前屈；(b)后伸；(c)左側彎；(d)右側彎；(e)左軸向旋轉；(f)右軸向旋轉
Figure7. Stress nephogram of lumbar L4~5 intervertebral disc under different states(a) flexion； (b) extension； (c) left lateral bending； (d) right lateral bending； (e) left axial rotation； (f) right axial rotation
3 討論
本文基于CT圖像數據，利用Mimics軟件無需對圖像進行任何形式的轉化，快速而準確地建立腰椎的三角面片模型，避免了信息的丟失。充分利用Hypermesh軟件結合由三角面片生成的四面體單元，最適用于具有復雜邊界曲面的不規則椎體結構的離散[21]，同時采用十節點保證了計算精度且很大程度上縮短了有限元前處理的時間，提高了建模的效率。基于“網格包絡”思想，巧妙運用共用面網格及網格偏置等網格劃分技巧，實現所有單元共節點，不需要在椎體與椎間盤處進行綁定約束，提高了模型的收斂性，大大提高了計算速度[22]。
腰椎整體的應力云圖反映了椎體在前屈、后伸、側彎和軸向旋轉時的力學分布情況，結果證明椎弓根附近是椎體多種應力集中的部位，后部結構所受應力只有通過椎弓根才能傳達到椎體。這就將有限元分析和臨床研究有效地結合在一起且結果能夠相互印證，再次說明模型的準確性和有效性。同時可以看出，椎體在各工況下受壓縮側應力均較大，在后伸、側彎和軸向旋轉時，可以明顯看到椎體后方小關節突承載著較大的接觸應力，這說明小關節突可以對抗后伸、側彎和軸向旋轉的作用很明顯，這也與實際研究相符[13]。另外，由研究結果可知，椎間盤的應力分布也反映了在各加載情況下椎間盤受壓側應力較為集中，并有向受拉側逐漸擴散消釋的趨勢。
本文有效地模擬了腰椎的生物力學特性，但也存在著一些不足，尚需要完善和改進。首先，CT圖像數據的精度和軟件處理的精度等制約了建模精度。其次，腰椎椎體、椎間盤及韌帶等結構材料的生物力學特性較為復雜，在材料賦值時往往需要簡便處理，因此使得材質特性與真實情況難以維系。最后，腰椎處于人體復雜多變的環境，不僅受到體重的影響，肌肉等軟組織在腰椎運動過程中的牽張作用也是不容忽視的。因此，完善腰椎的有限元模型進行深入的生物力學特性分析，是我們進一步探索的工作。
