基于CT影像的骨質疏松性胸腰椎壓縮骨折PKP術后三維有限元力學模型的建立與驗證

權禎，張曉剛，，秦大平，，宋敏，，張華，張宏偉，趙希云，王志鵬，馬濤，陳缽

1.甘肅中醫藥大學中醫臨床學院，甘肅蘭州730000；2.甘肅中醫藥大學附屬醫院脊柱外科，甘肅蘭州730020

前言

骨質疏松性胸腰椎壓縮骨折（Osteoporotic Thoracolumbar Compression Fracture,OVCF）是在骨質疏松的前提下，暴力釋放后胸腰椎局部應力集中，導致椎體發生骨質疏松性骨折，輕者椎體塌陷，高度丟失，脊柱后凸畸形，重者椎體爆裂骨折，壓迫脊髓神經，甚至截癱［1］。胸腰段在整個脊柱解剖結構中特殊，處于活動度較小的胸椎和活動范圍較大的腰椎的應力移行區域，不同模載下的生物力學環境復雜，在脊柱損傷中發生率高。胸腰椎骨折發生生物力學機制核心為“力學失衡”，應用有限元生物力學途徑研究壓縮骨折發生的力學機制及特點，完善骨質疏松性椎體壓縮骨折的防治策略至關重要，有限元分析作為生物力學研究的一種新興方法，在脊柱骨折生物力學研究中逐漸被重視，本文就圍繞OVCF經過經皮椎體后凸成形術（PKP）治療后三維有限元模型的建立展開，為OVCF內置物椎體強化治療的生物力學研究提供模型支撐。

有限元分析在20 世紀70年代初逐漸被應用于骨科生物力學研究［2］。1974年，Belytschko 等［3］首次應用有限元方法進行脊柱生物力學分析，使得有限元分析在脊柱生物力學研究中得到應用與拓展［4］。有限元分析是一種離散化的數理求解方法，建立高仿真的生物力學模型，對模型有限單元、自由度、節點添加設定限制，應用計算機軟件輔助模擬人體真實受力模式［5］。其核心是將無限個質點構成的自由連續體分割成有限的單元及節點，明確模型不同結構的材料屬性和受力特性，計算出每個單元的剛度矩陣，每個小單元的剛度矩陣最后集合成總剛度矩陣，運用數字表達出來［6-7］。人體的生理環境和承載特征極其復雜，傳統尸體標本建模很難適應脊柱生物力學研究，對椎體的應力測試僅局限于椎體表面應力，無法測量椎體內部的生物力學變化，而有限元分析彌補了這一缺陷，對椎體內部的力學分布特性進行測試。近年來隨著計算機軟件的開發及生物力學數理研究的不斷深入，有限元分析作為生物力學研究的新興模式而被廣泛應用，其具備模型逼真、精準測量、重復性強等優勢，彌補了傳統骨生物力學研究的諸多干擾因素和不足，在骨質疏松性椎體壓縮骨折、內固定置入治療及PKP 椎體強化治療的生物力學研究領域逐漸被深入應用［8］。

1 研究對象

選取甘肅中醫藥大學附屬醫院脊柱外科收治的T12 骨質疏松性椎體壓縮骨折患者1 例，女性，72 歲，身高165 cm，體質量52 kg，骨密度測量T＜-2.5，且并發T12 椎體壓縮骨折，確診為骨質疏松性T12 椎體壓縮骨折。排除T12 椎體壓縮骨折以外的脊柱其他病變。本研究通過醫院醫學倫理委員會審核，患者表示知情，并自愿參加實驗。

2 實驗設備及軟件

本實驗選用西門子64 排高分辨率螺旋CT（甘肅中醫藥大學附屬醫院放射科提供）；高配置臺式計算機輔助建模平臺（電腦配置：處理器：銳龍AMD R5 2600X；硬盤：金士頓480 G；8 GB 運行內存；8 GB 顯卡）。Mimics 19.0 建模軟件（Materialise 公司，魯汶，比利時）進行原始模型的提取與轉化；Geomagic Warp 2017軟件（Raindrop公司，馬布希爾希，美國）進行特征去除、光順、曲面實體擬合；Solidworks 2017軟件（達索公司，馬薩諸塞州，美國）進行零件結構的組裝與附屬結構生成；Ansys Workbench 17.0 軟件（Ansys 公司，賓州匹茲堡，美國）添加材料屬性、邊界條件、坐標及載荷設定及進行生物力學分析。

3 建模方法

3.1 基于Mimics軟件進行原始三維模型的提取與轉化

3.1.1 骨質疏松性胸腰段CT影像數據的重構與導入對胸腰段T11椎體上緣至L2椎體下緣進行多序列掃描，層厚0.625 mm，電壓140 kV，電流200 mA，獲取512×512像素的高清連續性圖層數據，提取骨組織窗層，并以DICOM格式保存［9］。將DICOM格式的骨質疏松性胸腰椎（T11-L2）的二維CT數據導入Mimics 19.0軟件，利用導入（Import Image）命令，對圖像的前后、左右、上下不同方位進行調整與重置。圖層位置調整完畢后，轉化完畢的影像資料會自動重組顯示為斷面影像，影像斷面序列在Mimics操作界面中分別以橫斷面、矢狀面及冠狀面的形式顯示。

3.1.2 灰度閾值設定及“骨骼－軟組織”邊界分割二維層面數據導入成功后，進行閾值的選擇與界定，點擊灰度閾值（Thresholding）命令，選擇并調整CT的灰度閾值，設定合理的窗寬，選擇閾值區間為239～1 809 HU。只有選擇恰當的閾值區間，才能將骨骼和軟組織進行精確分離，最大限度保留模型的原始形態（圖1）。實際操作中，原始模型提取過程中由于分割對象主要為骨骼組織，則常規選擇系統默認的骨骼（Bone）的灰度閾值。再利用編輯蒙版（Edit Masks）命令對圖層進行逐層分離，分別對T11、T12、L1、L2、T12 骨水泥與毗鄰軟組織進行準確分割。選擇恰當的分離工具圓形（Circle）、拉索（Lasso）及方形（Square）等，提取過程中不同的解剖部位選用不同的分離工具，應用上述分離工具準確擦除椎體與相鄰組織的銜接點［10］。要求精準細致，盡量保留所需結構，刪除冗雜結構，若擦除與保留操作銜接不準，則易導致下一步模型重建的失真或表面缺失。

圖1 T12椎體的閾值選擇與蒙版分離界面Fig.1 Threshold selection of T12 vertebral body and mask separation interface

3.1.3 漸進化區域增長及三維原始單結構模型提取轉化單結構模型提取與轉化的前提是漸進化的區域增長（Region Growing），具體而言就是將椎體與毗鄰組織結構精確分割，最大限度地保持單結構模型的原始形態（圖2）。Mimics軟件的3D 視圖模型預覽（Toggle Mask 3D Preview）命令可對模型椎的分割情況進行即時追蹤與預覽（圖3），分析模型椎是否被準確分離提取。若模型椎與毗鄰組織有限接觸均可導致區域增長失敗，蒙版分離不足則結構粘連，分割過度則模型結構丟失，內部缺損需用腔隙填充（Cavity Fill）命令進行編輯，邊緣缺損需用邊緣繪制（Draw）命令進行修補［11］。待模型椎區域增長成功后，利用3D 模型計算（Caulatate 3D）命令進行原始模型提取與轉化，分別構建T11、T12、L1、L2 及骨水泥的三維原始模型，最后將提取完畢的模型椎以STL 格式儲存（圖4）。

圖2 模型椎的區域增長與轉化Fig.2 Regional growth and transformation of model vertebrae

圖3 模型椎二維CT平面中三維視圖Fig.3 Three-dimensional view of model vertebrae in two-dimensional CT plane

圖4 提取完畢后的T12高質量原始模型Fig.4 High-quality original model of T12 after extraction

3.2 基于Geomagic軟件的T11-L2節段三維曲面模型優化與構建

從Mimics 軟件提取的三維模型表面粗糙，內部冗雜，無法進行曲面片劃分、擬合實體等操作，不滿足脊柱生物力學分析條件，因此，需對原始三維模型的結構缺陷進行優化處理。

3.2.1 原始三維模型的導入及“多邊形”命令下模型的優化處理Mimics 軟件中提取的三維原始模型存在結構缺陷，不符合生物力學分析特性，需將原始模型導入Geomagic Warp 2017 軟件進行前處理。開始/導入命令將STL 格式的原始模型導入優化程序，開始命令僅可導入一個文件，進行單椎體優化操作時，可以使用此命令，導入命令則可同時導入多個模型椎，便于多模型同時優化［12］。模型導入后，軟件會自動推薦網格醫生命令，由于模型椎的內部冗雜及邊緣缺失，暫不推薦使用網格醫生命令進行模型檢測（圖5）。

圖5 原始三維模型椎導入與“多邊形”修飾界面Fig.5 Import of original 3D model vertebrae and"polygon"modification interface

模型優化處理中椎體內部冗雜刪減、邊緣孔彌合是關鍵，由于在CT解剖斷面提取層面的邊緣空洞未被完全彌合，部分邊緣空洞所形成的表面邊界向椎體內部發生“逃逸”，造成椎體內部多余的冗雜結構（圖6），限制了曲面實體的擬合。多邊形/填充單個孔命令將對模型邊緣孔和內部孔進行修補，并通過曲面/斜率追蹤命令，進行邊緣孔“搭橋”處理（圖7），降低模型彌合后的變形系數［13］。椎體內部冗雜及邊緣孔修復完畢后，模型椎表面結構仍凹凸不平，呈現為釘狀物、凹陷及粗糙特性，可利用多邊形/去除特征/刪除釘狀物/光順命令進行模型椎表面優化，待模型表面光滑、輪廓清晰、結構完整后進入精確曲面程序擬合實體。

圖6 模型椎內部冗雜剔除Fig.6 Elimination of miscellaneous inside model vertebrae

圖7 “搭橋”輔助修復邊緣孔Fig.7 "Bridging"assisted repair of the edge hole

3.2.2 對優化模型進行精確曲面、構造格柵及擬合實體建模曲面實體擬合的前提是繪制合理的輪廓線，輪廓線引導曲面片的繪制至關重要，規整的輪廓線有利于有限元分析前處理階段接觸關系的設定和求解（圖8）。待模型外觀修正與真實狀態完全一致后進行曲面片分割，使用精確曲面命令中的探測輪廓線工具，此時軟件本身會自動根據模型的邊界輪廓繪制輪廓線，但是軟件本身繪制的輪廓線比較粗糙扭曲，不規整的輪廓線對曲面片的生成造成影響，此處則需要手動編輯輪廓線（圖9）。將曲率敏感度降低至20，軟件自動探測的輪廓線將消失，需要手動根據椎體表面邊界繪制輪廓線，為曲面片的生成提供線性引導（圖10）。至生成規整、正確的曲面片后，重畫網格（圖11）。構造格柵可以檢查到曲面分割的缺陷，通過編輯曲面命令可以進行修改，格柵越均勻則曲面表面越光滑，格柵構造的均勻程度將會影響到后期的曲面提取和網格劃分等諸多操作（圖12）［14］。生成的皮質骨的曲面模型分別保存為STP格式（圖13）。松質骨位于皮質骨的內部，通過多邊形/偏移命令進行建模，將皮質骨向內部偏移1.00 mm，到達松質骨層面，此時松質骨模型的外觀形態及邊緣較為銳利，仍需用多邊形命令中的特征去除、磨砂、光順等工具對模型進行處理，待模型前處理接近于椎體松質骨真實形態時，進行網格重畫、繪制輪廓線、構造曲面片、編造格柵及擬合實體等順序性建模（圖14～圖16），最終構建T11-L2節段的椎體及骨水泥（圖17～圖18）的擬合實體模型。

圖8 探測輪廓線Fig.8 Detecting profile

圖9 編輯輪廓線Fig.9 Editing profile

圖10 構造曲面片Fig.10 Constructing curved surface

圖11 重畫網格Fig.11 Redrawing grid

圖12 構造格柵Fig.12 Constructing grid

圖13 擬合曲面實體Fig.13 Fitting curved surface entity

圖14 松質骨輪廓線Fig.14 Contour of cancellous bone

3.3 基于Solidworks 軟件的T11-L2 節段模型組裝與結構重建

圖15 松質骨曲面片構造Fig.15 Construction of cancellous bone curved plate

圖16 松質骨構造格柵Fig.16 Structural grid of cancellous bone

圖17 松質骨擬合實體Fig.17 Cancellous bone fitting entity

圖18 彌散型骨水泥Fig.18 Dispersive bone cement

3.3.1 單椎體零件結構的生成將STP 格式優化修飾后的曲面實體導入Solidworks軟件中，進行模型組裝及單元結構的添加。利用菜單工具欄中的“文件打開”命令，對已經生成的皮質骨擬合曲面和松質骨擬合曲面及骨水泥逐一打開，分別導入后將生成單一椎體皮質骨、松質骨、骨水泥的零件文件，將生成文件以零件（SLDPRT）格式進行保存，如此反復，將所有整體三維模型的各個部分均以零件（SLDPRT）格式保存，生成T11-L2 節段的皮質骨、松質骨、骨水泥的單椎體零件。

3.3.2 基于“原點重合”理論的裝配體零件合成及干涉檢查生成的模型結構零件是松散的，具有若干單獨“坐標”，需按照CT 三維坐標“原點重合”理論將單零件進行裝配，然后將其組裝為裝配體零件，方便下一步整體操作。點擊“新建裝配體”命令，以“裝配體原點重合原理”對離散的T11、T12、L1、L2 及T12 椎體內部的骨水泥進行結構整體化裝配（圖19），建立T11-L1 節段包含皮質骨、松質骨、骨水泥的整體模型，并以SLDPRT 格式保存。待整體模型組裝完畢后，需要對模型局部接觸關系進行干涉檢查，尤其是關節突關節的局部接觸構建，避免模型曲面生成過程中“邊界逃逸”后使得模型局部變形，如果整體模型干涉檢查陽性，則會導致模型失真和力學測試失敗，需要進一步返回Geomagic 軟件中對失真模型局部進行重建，以確保模型的高仿真效應和精準性（圖20）。

圖19 “離散”的T11-L2單結構零件Fig.19 "Discrete"single-structure part of T11-L2

圖20 組裝后T11-L2節段整體模型Fig.20 Overall model of T11-L2 after assembly

3.3.3 椎間盤、終板、關節突關節軟骨的立體模型構建椎間盤、終板、關節突關節軟骨由于其結構特性的限制，柔性結構在CT 骨窗層面的分辨率較低，提取較為困難，要依靠高質量的MRI 數據進行提取。本實驗是基于胸腰椎CT 數據的仿真模擬，無法直接提取，需要在Solidworks 實體建模軟件中進行重構［15］。

椎間盤的生成以椎體表面網格結構基準面為基礎，“三點一面”理論結合草圖框制，在椎體表面的網絡結構上選擇同一平面的不同三點，構建平行于椎體上下表面的基準面，借助于測量命令，繪制椎間盤的平面輪廓，利用“拉伸-凸臺”命令，初步構建椎間盤的立體形態，此時椎間盤與上椎體的下表面和下椎體的上表面立體結構重合，不符合脊柱的生理特性，需要對模型局部結構進行微調。“等距曲面”命令可以刪除重合的冗雜結構，并且等距曲面刀具距離須設置為零，至此T11-12、T12-L1、L1-2 胸腰段椎間盤的立體結構構建完畢（圖21）。

圖21 T11-L2節段椎間盤組織Fig.21 Intervertebral disc tissue at T11-L2

關節突關節的構建類似于椎間盤建模，同樣需要借助于“基準面-草圖-等距曲面”命令，在上位椎體的下關節突和下位椎體的上關節突結構重合部位選擇基準面并繪制關節軟骨平面結構，再利用“拉伸-凸臺”命令，構建關節軟骨的立體結構，將等距曲面厚度設置為零，對冗雜結構進行刪減，最終建立T11-12、T12-L1、L1-2 胸腰段左右關節突關節的立體結構（圖22）。

圖22 T11-L2節段關節突關節Fig.22 Facet joint at T11-L2

髓核結構外于椎間盤內部，利用測量工具測量髓核體積，髓核大小占椎間盤體積的50%，重建草圖，將草圖大小牽拉為接近于髓核體積，利用等距曲面進行切割，生成纖維環和髓核［16］圖23。上下軟骨終板從立體椎間盤凸臺結構中進行分離，軟骨終板的厚度約為1.5 mm，利用“等距曲面”命令，調整等距曲面厚度為1.5 mm，以接近于真實解剖結構中軟骨終板厚度（圖24）。皮質骨外殼與松質骨內核的分離借助于“移動復制實體/組合”命令，對T11、T12、L1、L2松質骨實體及骨水泥團塊進行復制，“組合”命令實質是將皮質骨與內部的松質骨、骨水泥重合部分進行刪減，得到單純的皮質骨外殼和內部的松質骨、骨水泥團塊，以提高模型仿真度。

圖23 T11-L2節段髓核、纖維環Fig.23 Nucleus pulposus and annulus fibrosus at T11-L2

圖24 T11-L2節段軟骨終板、骨水泥結構Fig.24 Cartilage endplate and bone cement structure at T11-L2

3.4 基于Ansys軟件構建有限元分析前處理模型

3.4.1 彈性模量和泊松比的設定彈性模量和泊松比是表示在縱向壓縮和橫向拉伸下材料結構強度和剛度的物理量（表1），就椎體而言，彈性模量和泊松比越大，則材料強度及剛度越強，相反則越低（圖25）。骨質疏松模型的材料屬性均低于正常骨結構，因此將骨質疏松模型松質骨的所有骨結構降低67%，將皮質骨、終板和后部的彈性模量降低34%［17］。

表1 三維有限元模型材料屬性Tab.1 Material properties of three-dimensional finite element model

圖25 模型結構的材料參數添加視圖Fig.25 Structural material properties addition view

3.4.2 模型坐標系的重置與接觸關系設置坐標系為力的加載與釋放提供了約束，通過坐標系（construc-tion geometry）命令，可以對模型的坐標進行修正與重置，以適應人體真實靜力測試的要求，Z 軸設定為垂直應力加載，模擬垂直壓縮受力，同時模型圍繞垂直力模擬左右旋轉受力，Y軸設定為水平向后的應力加載，模擬前屈、后伸受力，后伸時為正向力，前屈時為反向力，X 軸設定為水平向左力的加載，模擬左右側彎時受力載荷，左側彎曲時為正向力，右側彎曲時為反向力［18］。

合理的結構間接觸是應力傳遞的必要前提，其接觸關系主要有固定接觸和有限接觸，“面-面接觸”、“點-面接觸”是模型接觸的主要形式，可以通過接觸（connection）命令對接觸關系進行設定，將椎體-終板、皮質骨-松質骨、終板-椎間盤設定為綁定（bonded）模式，關節軟骨具有一定的活動度但又不會脫離活動軌跡，其接觸為一定無摩擦位移滑動范圍內的有限接觸，故將其接觸模式設定為不綁定的（no sperarete）模式。

3.4.3 模型網格化處理模型網格化將結構精確化，劃分為有限個利于求解和控制的單元和節點（表2），不同結構的差異性與復雜性導致模型網格化也存在差異，適當的網格大小影響著模型求解精細度，利用Generate Mesh（網格生成）命令進行網格重置，皮質骨、松質骨結構網格大小設定為2.0 mm，椎間盤、關節軟骨、終板、骨水泥等附屬結構網格大小設定為1.5 mm。

表2 結構的節點與單元Tab.2 Nodes and units of different structures

3.5 彈簧韌帶的添加及剛度系數設定

腰椎韌帶結構復雜，在CT 骨窗層面難以提取，需要在Ansys17.0 軟件中用彈簧韌帶對生理韌帶結構進行替代，利用彈簧韌帶（spring elements）命令進行添加，并將其生理特性設置為僅拉伸（tensiononly）模式，依次添加前縱韌帶、后縱韌帶、黃韌帶、棘間韌帶、棘上韌帶、橫突間韌帶等，并將其模擬特性設置為線性彈性材料。

3.6 力學加載及邊界條件設定

人體脊柱的整體運動狀態及機制是復雜的，我們無法在軟件中進行真實還原，高仿真下模擬不同載荷求解工況，在圍繞Z 軸在T11 椎體表面施加垂直向下500 N 的壓力，以模擬壓縮載荷（圖26），-7.5 N·m的反向扭力，模擬前屈載荷（圖27），圍繞Y軸施加水平向后的7.5 N·m的正向扭力，模擬后伸載荷（圖28）；X 軸設定為水平向左的7.5 N·m 扭力和向右的-7.5 N·m扭力，模擬左右側彎時受力載荷（圖29～圖30）；圍繞Z軸施加7.5 N·m左旋力和-7.5 N·m的反向右旋力，模擬左右側彎載荷（圖31～圖32）。為了使其仿真度達到最優，需要對T11-L2節段的運動軌跡進行約束，實驗中固定了L2椎體下表面及下關節區域（圖33），限制其自由度［19］。

圖26 軸向受力加載Fig.26 Axial force-loading

圖27 前屈受力加載Fig.27 Anterior flexion force-loading

圖28 后伸受力加載Fig.28 Posterior extension force-loading

圖29 左側彎受力加載Fig.29 Left flexion force-loading

圖30 右側彎受力加載Fig.30 Right flexion force-loading

圖31 右旋轉受力加載Fig.31 Right rotation force-loading

圖32 左旋轉受力加載Fig.32 Left rotation force-loading

圖33 L2椎體約束自由度Fig.33 Constraint freedom of L2 vertebral body

4 OVCF有限元模型的有效性驗證

模型建立成功后，利用有限元軟件進行計算，在進行生物力學實驗之前，對三維有限元模型進行生物力學驗證是必要前提，以驗證模型的精確性、可靠性。脊柱有限元模型有效性驗證的方法主要有直接驗證和間接驗證，直接驗證需要與離體模型進行驗證，可靠性及難度較大。間接驗證是將有限元模型的角位移與既往生物力學實驗進行對比驗證，本模型采用了精確化建模路徑，精細劃分網格及單元，模型仿真度高，準確性強，有效性驗證合格，可進行骨質疏松性胸腰椎壓縮骨折PKP術后生物力學測試［20］。記錄模型在各個運動方向上的最大角位移即活動范圍（Range of Motion,ROM）。之后參照標準脊柱有限元實驗和生物力學實驗的結果，對比所構建模型的屈曲度變化判斷所構建模型的有效性，評判其是否可以適用于進一步的模型處理和有限元分析（表3）。

表3 T11-L2節段模型ROM與既往數據T11-L2胸腰段ROM活動度均值（°）Tab.3 Average range of motion(°)of T11-L2 between model and previous data

5 結語

隨著計算機應用技術和有限元分析軟件的開發與升級，有限元作為一種生物力學分析的有效方法被深入應用于骨科生物力學研究當中，通過仿真模擬途徑揭示骨科諸多現象的生物力學機制，就骨質疏松性椎體壓縮骨折而言，其在椎體壓縮骨折及慢性疼痛的生物力學機制應用廣泛［23］。以往脊柱骨折的生物力學研究以尸體建模為主，通過機械壓力機測試材料生物力學性能，模型單一、可重復性差，與人體真實復雜的生物力學受力情況差異性大，而有限元分析作為一種數據仿真模擬，具有可重復性、仿真度高、精細性強等優勢。

骨質疏松性椎體壓縮骨折的生物力學模型仍需要優化升級，目前建模數據是基于DICOM格式的CT數據的建模，對CT 數據質量要求比較高，高精度的三維有限元模型必須以多序列、薄層、高分辨率的CT數據為基礎，但只能對骨窗組織進行提取與轉化，柔性結構很難被提取，椎間盤、軟骨終板、關節軟骨及韌帶等組織在原始模型中依然缺失，需要在Solidworks軟件中進行附屬結構的重建及Ansys軟件中進行韌帶重構。隨著生物力學測試的精準度不斷提高，基于3.0T 高精度MRI 數據的模型提取與建立將是有限元建模的熱點與難點，由于核磁對軟組織的識別度比較高，不同組織在MRI 序列中均有不同信號的顯影成像，依托于腰椎的MRI 數據有望直接提取椎間盤、韌帶、關節軟骨、軟骨終板等柔性結構，部分學者對基于MRI 數據建模進行了探索［24］，但腰椎韌帶、肌肉結構混雜，單純利用腰椎MRI 數據提取腰椎包括韌帶在內的整體模型案例較少，后期隨著有限元技術的開發與應用，為提高模型仿真度，利用核磁構建腰椎三維整體模型已成為必然趨勢。