基于有限元仿真的數(shù)字化脊柱側(cè)彎矯形研究

盛然，蔣國璋，劉融

1.武漢科技大學(xué)機(jī)械自動化學(xué)院，湖北武漢430081；2.武漢市普仁醫(yī)院骨科，湖北武漢430081

前言

脊柱側(cè)彎是一種嚴(yán)重影響患者身體正常發(fā)育的疾病，通常會隨著患者的生長發(fā)育而持續(xù)加重，并且生長發(fā)育速度越快加重的速度越快。據(jù)統(tǒng)計，目前我國脊柱側(cè)彎患者超過300 萬人，并以每年30 萬人的速度遞增，其中超過半數(shù)為青少年［1］。中輕度的脊柱側(cè)彎多帶來身體形態(tài)與外觀異樣，重度則會影響患者的心肺功能，壓迫器官，對青少年身心帶來巨大的壓力與負(fù)擔(dān)，不利于他們的正常生活與學(xué)習(xí)。因此，關(guān)于脊柱側(cè)彎的康復(fù)治療越來越受到家長、醫(yī)生以及全社會的關(guān)注［2］。目前，國際上公認(rèn)的脊柱側(cè)彎治療手法主要有兩種，手術(shù)治療與保守治療。手術(shù)治療創(chuàng)傷大、風(fēng)險高、治療費(fèi)用高，一般來說Cobb 角大于40°時建議進(jìn)行手術(shù)，大多數(shù)青少年達(dá)不到手術(shù)指征，也不愿意采用手術(shù)治療。因此保守治療就成為了患者首選，保守治療包括脊柱側(cè)彎矯形器和其他輔助運(yùn)動鍛煉［3-6］。然而當(dāng)前脊柱側(cè)彎的矯形方案中，僅依據(jù)三點(diǎn)力矯形原理與臨床治療經(jīng)驗(yàn)對患者進(jìn)行矯形治療［7］。對矯形力大小、施力點(diǎn)、施力方向等并沒有明確其具體數(shù)值，在矯形器的制備與使用中，往往會因此導(dǎo)致矯形過度或不足，影響患者治療進(jìn)展［8］。為保證患者安全高效地完成矯形治療，提出一種更加科學(xué)合理的系統(tǒng)化矯形方案迫在眉睫。

1 數(shù)字化矯形方案分析

1.1 脊柱側(cè)彎數(shù)字化信息提取

為了實(shí)現(xiàn)患者矯形信息數(shù)字化，利用Matlab 軟件進(jìn)行編程計算，自動獲得患者Cobb 角與最大偏移量數(shù)值，并且通過數(shù)學(xué)建模計算對應(yīng)的矯形力與位置。

安排經(jīng)驗(yàn)豐富的脊柱外科醫(yī)生，利用Mimics 19.0軟件從患者的CT斷層掃描數(shù)據(jù)中，以人工標(biāo)定的方式，提取出冠狀面每個節(jié)段椎體上下終板的坐標(biāo)，導(dǎo)入Matlab數(shù)據(jù)欄中。利用cftools工具箱對離散的數(shù)據(jù)點(diǎn)進(jìn)行曲線擬合，綜合考慮擬合效果、均方根誤差（RMSE）、殘差誤差值（adjusted R-square）等影響，最終選定三階高斯函數(shù)進(jìn)行擬合［9］。如圖1所示。

圖1 三階高斯函數(shù)擬合曲線圖Fig.1 Third-order Gaussian function fitting curve

擬合曲線方程確定后，找到此彎的上下端椎坐標(biāo)，利用編程算法偽代碼（圖2），得到編程測算結(jié)果（圖3），計算獲得患者的Cobb 角為36.54°，最大偏移位置在距原點(diǎn)170.10 mm處，最大偏移量為40.76 mm。

圖2 編程算法偽代碼Fig.2 Programming algorithm pseudocode

圖3 編程測算結(jié)果圖Fig.3 Calculation results of programming

1.2 矯形力與矯形位置分析

脊柱側(cè)彎矯形支具設(shè)計中的矯形力大小與位置往往通過臨床經(jīng)驗(yàn)設(shè)定，尚未出現(xiàn)可靠的理論依據(jù)支撐［10-11］?？紤]到脊柱在矯形力作用下的形變類似于平面曲桿的受力方式，可以用彎曲梁模型對主彎部分進(jìn)行等效代替［12］。將近骨盆一端設(shè)為原點(diǎn)視作固定鉸支座，另一端為可動鉸支座。

由于脊柱側(cè)彎矯形過程為矯形力在空間中始終垂直于椎體部分，在模型等效計算中，可將脊柱的彎曲變形視作純彎曲。而截面形狀近似為圓形截面進(jìn)行計算，本例中截面直徑為43.82 mm（圖4）。

圖4 等效橫截面直徑示意圖Fig.4 Diagram of equivalent cross section diameter

假設(shè)在脊柱不同位置受到n個集中力作用，分別為F～Fn，集中力對應(yīng)位置分別為x1～xn，如圖5所示。其中，RA為固定鉸支座受矯形力時的支反力，RB為可動鉸支座支反力，l為主彎部分上下端椎之間的距離，本例中l(wèi)=152.16 mm。

圖5 等效彎曲梁模型受力分析圖Fig.5 Force analysis diagram of equivalent curved beam model

以x左端為研究對象，推導(dǎo)彎矩M與截面位置的函數(shù)為：

為在彎曲梁等效模型中計算出相關(guān)應(yīng)力與應(yīng)變，設(shè)r為中性層半徑；E綜為脊柱綜合彈性模量；Φ為支反力與x位置的夾角；△Φ為施力位置x處的矯形恢復(fù)角。由胡克定律σ=Eε以及彎矩的定義，推導(dǎo)可得到在橫截面x處彎矩：

其中，R為側(cè)彎偏移量，A為等效截面面積。在彎矩M的作用下，主彎部分在橫截面上以中性軸為基準(zhǔn)，一側(cè)受拉，一側(cè)受壓，則在脊柱上應(yīng)力總和為0。而在計算曲桿彎曲應(yīng)力時，圓形截面的中性層半徑的計算公式為：

最大偏移量R0=40.76，截面直徑d=43.82 帶入計算可得r=37.56。由上式可得到在x處的矯形力與矯形角度的函數(shù)關(guān)系式：

最終確定E綜=0.285 MPa。

本例患者脊柱在冠狀面未見明顯畸形，故僅考慮其矢狀面方向的矯形方案。本例中患者側(cè)彎為單彎，主彎部分主要在胸腰椎段，上下端椎長度l=152.16 mm，最大偏移量R0=40.76 mm，Cobb 角為36.54°。目前臨床治療經(jīng)驗(yàn)表明，矯形角度不足可能導(dǎo)致矯形效果不明顯，甚至引起側(cè)彎狀況惡化；而過大的矯形角度將對患者的康復(fù)過程造成不適，嚴(yán)重者可導(dǎo)致不可逆轉(zhuǎn)的脊柱損傷［13］。臨床上，人體可承受范圍內(nèi)最佳矯形角度為脊柱側(cè)彎角度的，即

再根據(jù)脊柱側(cè)彎三點(diǎn)力矯形原理計算得到最優(yōu)矯形位置x與最優(yōu)矯形力大小F，即在距離原點(diǎn)x=170.10 mm 的位置施加65.2 N 的矯形力，為最優(yōu)矯形方案。

2 脊柱側(cè)彎數(shù)字化仿真

本研究以患者的CT 斷層掃描數(shù)據(jù)為基礎(chǔ)，對患者脊柱進(jìn)行數(shù)字化三維重建。通過Mimics19.0 軟件CT 骨分割功能，提取患者全脊柱骨骼信息。使用3D蒙版重建功能對患者脊柱全長進(jìn)行初步重建，如圖6a所示。

圖6 側(cè)彎脊柱的數(shù)字化仿真Fig.6 Digital simulation of the scoliosis

本例脊柱側(cè)彎患者為C型單彎，主彎部位為腰彎，為了研究矯形力與矯形位置對患者矯形效果的影響，應(yīng)對主彎部分進(jìn)行精確建模。利用Mimics19.0軟件的remesh功能，將3D蒙版導(dǎo)入3-Matic中對其進(jìn)行實(shí)體化。為保證重建模型的精度，將每個椎體導(dǎo)入Geomagic Studio 2013軟件中進(jìn)行模型孔洞修復(fù)與壞片修復(fù)處理，并對實(shí)體化的椎體進(jìn)行體網(wǎng)格劃分，簡化三角面片數(shù)量，以確保計算速度。

將優(yōu)化后的椎體導(dǎo)入3-Matic 軟件中，準(zhǔn)備進(jìn)行椎間盤重建。由于在脊柱側(cè)彎矯形過程中，在矯形力作用下椎體形變可忽略不計，發(fā)生形變的主體實(shí)質(zhì)上是連接每段椎骨間的椎間盤，因此重建出符合人體生物結(jié)構(gòu)的椎間盤十分重要。當(dāng)前，在CT 斷層掃描圖像中無法通過灰度值劃分出椎間盤此類軟組織，故需要在椎體模型間對椎間盤進(jìn)行手動重建。

椎間盤結(jié)構(gòu)由外圍的纖維環(huán)與內(nèi)部的髓核構(gòu)成。在3-Matic中，使用Mark功能標(biāo)記兩塊椎骨的上下表面，利用解剖重建功能進(jìn)行椎間盤的重建，并以向內(nèi)抽殼方式構(gòu)建髓核，剖面圖如圖6b所示

不同于以往重建案例中通過workbench 對脊柱的單元類型、彈性模量與泊松比進(jìn)行均一化材料屬性賦值［15］（表1），本例經(jīng)過實(shí)體網(wǎng)格劃分后，將模型導(dǎo)入Mimics19.0 后對骨骼進(jìn)行基于灰度值劃分的材料屬性賦值。由于骨骼本身的密度層次復(fù)雜，基于灰度值賦值相比均一化賦值而言更能反映患者骨骼結(jié)構(gòu)的真實(shí)情況［16］（表2）。椎間盤的纖維環(huán)與髓核部分，按照表1中彈性模量與泊松比進(jìn)行賦值。

表1 ANSYS材料屬性賦值Tab.1 ANSYS material attribute assignment

表2 Mimics骨骼灰度值賦值經(jīng)驗(yàn)公式Tab.2 Empirical formula for the assignment of Mimics bone gray value

除了椎體與椎間盤外，各類韌帶也是脊柱的重要組成部分。脊柱的連接依賴前縱韌帶、后縱韌帶、棘上韌帶等多種韌帶。為了更好地模擬患者在各類矯形方案下脊柱的形變狀況，本例中通過添加彈簧連接模擬韌帶連接，如圖6c所示。至此，脊柱數(shù)字化仿真建模完成。

3 結(jié)果

利用前文建立的脊柱側(cè)彎三維模型，對計算出的最優(yōu)矯形方案進(jìn)行有限元分析驗(yàn)證。三維模型底端設(shè)為固定端，空間移動設(shè)置為0，在距離原點(diǎn)170.10 mm 的位置對脊柱側(cè)彎三維模型施加65.2 N的矯形力，應(yīng)力云圖如圖7所示。

圖7 側(cè)彎脊柱模型有限元仿真結(jié)果Fig.7 Finite element simulation results of scoliosis model

等效彈性應(yīng)變云圖反映出應(yīng)變主體部分為椎間盤，椎體部分不受形變影響，證明了脊柱側(cè)彎三維模型的有效性?？傋冃卧茍D顯示，在此矯形方案下，最大位移15.435 mm，形變后Cobb 角19.86°，相比矯形前的36.54°減少了16.68°，基本滿足的矯形要求，證明了本矯形方案的有效性。

4 討論

本研究基于患者全脊柱CT 斷層掃描數(shù)據(jù)，利用Matlab軟件編程，計算出患者的Cobb角、最大偏移量與偏移位置。能夠在治療前期，高效準(zhǔn)確地為醫(yī)生提供脊柱側(cè)彎患者的數(shù)字化信息，減少了人為測量導(dǎo)致的誤差與時間。輔助醫(yī)生快速判斷患者是否需要進(jìn)行手術(shù)治療，能否以佩戴矯形器的方式進(jìn)行矯形治療。

將脊柱側(cè)彎部分以圓截面彎曲梁模型進(jìn)行等效替代，進(jìn)而分析出脊柱側(cè)彎矯形力與矯形角度的關(guān)系，為患者提供了具體的矯形方案。區(qū)別于依賴臨床經(jīng)驗(yàn)的傳統(tǒng)矯形方案，本研究提出的矯形方案有相關(guān)力學(xué)理論支持，并且經(jīng)過患者脊柱側(cè)彎三維模型的有限元驗(yàn)證其有效性，更具科學(xué)性。

需要引起注意的是，若患者脊柱在矢狀面存在畸形時，需要將矢狀面脊柱以相同的方式在Matlab中進(jìn)行計算分析，并且在確定矯形力與矯形位置時在空間中計算合力的作用點(diǎn)與力的大小，以保證脊柱側(cè)彎矯形的全面性。考慮到本例中患者為C 型脊柱側(cè)彎，且主彎位置在腰部，近骨盆側(cè)才能等效為固定端進(jìn)行力學(xué)分析。當(dāng)患者的側(cè)彎類型與主彎位置發(fā)生變化時，需結(jié)合具體情況進(jìn)行等效建模與相關(guān)分析［17-19］。

基于數(shù)字化三維建模的脊柱側(cè)彎矯形研究旨在為患者提供最優(yōu)化矯形方案。不僅能在治療前期準(zhǔn)確地為醫(yī)生提供判斷患者側(cè)彎狀況的依據(jù)，也為后續(xù)脊柱側(cè)彎矯形器的設(shè)計打下了基礎(chǔ)。

當(dāng)前我國正在大力推進(jìn)智慧醫(yī)療進(jìn)程，在脊柱側(cè)彎矯形研究方面，智能化程度遠(yuǎn)遠(yuǎn)不足。未來智能化脊柱側(cè)彎矯形的發(fā)展，將集中在結(jié)合人工智能算法，自動生成個性化矯形方案。從CT 斷層掃描數(shù)據(jù)，自動提取患者身體信息，生成最優(yōu)矯形方案，并結(jié)合虛擬現(xiàn)實(shí)技術(shù)，直觀反映脊柱側(cè)彎矯形過程與個體化矯形器定制方案。筆者也將以此課題作為今后的關(guān)注重點(diǎn)，繼續(xù)進(jìn)行深入研究。