本發(fā)明涉及生物醫(yī)學工程領域，特別涉及一種基于組織分化的長骨骨折愈合仿真系統(tǒng)。

背景技術：

骨折是一種常見的創(chuàng)傷，骨折的高發(fā)性使得對骨折機理及促進愈合的研究尤為迫切，一旦骨折發(fā)生，與其它組織損傷修復不同的是，骨折不是靠纖維結締組織連接，而是骨組織的完全再生。然而，并不是所有的骨折都可以完成愈合，有時會發(fā)生延遲愈合甚至是不愈合。骨折延遲愈合或者不愈合會引起患肢疼痛，功能障礙，導致患者失業(yè)，由此造成很大的社會經(jīng)濟負擔。因此，盡管關于骨折愈合的研究一直備受關注，但仍有5％～10％的骨折因各種原因發(fā)生延遲愈合甚至是不愈合。

骨折愈合受到特定的幾何因素、力學因素、生物學因素影響，良好的幾何因素、力學因素、生物學因素得到良好愈合效果。反之則會導致骨折的延遲愈合甚至是不愈合。目前缺少能夠精確表達骨折愈合這一復雜過程的計算機仿真系統(tǒng)。在骨折愈合仿真系統(tǒng)中存在以下缺陷：

1.沒有建立專門針對患者的個體化模型；

2.力學因素與骨折愈合過程沒有一個確定性關系；

3.骨折愈合區(qū)域的模型和生物力學材料設置過于簡化；

4.沒有在同一個仿真系統(tǒng)中綜合考慮力學因素和生物學因素對骨折愈合的影響。

技術實現(xiàn)要素：

本發(fā)明的目的是為了解決現(xiàn)有的骨折愈合仿真中不能綜合模擬力學因素和生物學因素對骨折愈合過程的影響，骨折愈合生物力學模型材料設置過于簡化的缺點，而提出的一種基于組織分化的長骨骨折愈合仿真系統(tǒng)。

本發(fā)明的目的通過下述技術方案實現(xiàn)：一種基于組織分化的長骨骨折愈合仿真系統(tǒng)，其特征在于，所述系統(tǒng)包括：

骨折區(qū)域幾何建模模塊、骨折區(qū)域生物力學有限元分析模塊、骨痂單元組織分化模塊和程序終止判斷模塊；

骨折區(qū)域幾何建模模塊用于根據(jù)導入的二維斷層掃描圖像數(shù)據(jù)，經(jīng)過圖像預處理后進行骨折部位的三維表面幾何模型的建立；

骨折區(qū)域生物力學有限元分析模塊用于對建立好的骨折區(qū)域模型進行網(wǎng)格劃分，施加外部載荷和設置邊界條件；

骨折區(qū)域生物力學有限元分析模塊還用于初始骨折區(qū)域環(huán)境的設置；

骨折區(qū)域生物力學有限元分析模塊還用于計算單元力學刺激；

骨痂單元組織分化模塊用于對單元內(nèi)組織分化進行仿真，使單元內(nèi)各組織含量得到更新，從而使單元材料屬性得到更新，從而得到下一迭代步中所需要的單元力學刺激；

程序終止判斷模塊用于判斷程序是否終止，若不滿足終止條件，程序進行下一迭代步；若滿足程序終止條件，則程序結束并輸出愈合時間。

本發(fā)明的有益效果為：

1.本發(fā)明提出的一種基于組織分化的長骨骨折愈合仿真系統(tǒng)是基于windows開發(fā)語言平臺來開發(fā)軟件，通過自主編程實現(xiàn)骨折愈合過程的動態(tài)模擬，基于對話框的形式，易于操作，培訓周期短；

2.將骨折區(qū)域設置為雙相多孔彈性模型，相比單相模型，更加符合骨折愈合區(qū)域的生物特性，使仿真結果更加精確；

3.將血供作為變量引入骨折愈合仿真系統(tǒng)中，能夠更加準確的模擬骨折愈合的過程以及血液對骨折愈合的影響；

4.利用模糊邏輯的方法對組織分化進行仿真，相較于傳統(tǒng)的采用偏微分方程組進行建模的方法，更加便于理解。減少了在建模過程中偏微分方程組的數(shù)量，便于建模且減少了仿真時間；

5.通過構建骨折愈合仿真系統(tǒng)，可以對醫(yī)生制定最優(yōu)的手術方案提供指導，進而提高手術成功率、提高骨折愈質(zhì)量，減少骨折不愈合和延遲愈合的情況；

6.通過構建骨折愈合仿真系統(tǒng)，可以對建立的仿真模型進行多次重復實驗研究，減少真實的生物實驗，節(jié)省時間，提高效率，節(jié)省費用，避免人道主義的爭議。

綜上，本發(fā)明的仿真平臺克服了現(xiàn)有技術的缺點與不足。

附圖說明

圖1為基于組織分化的骨折愈合仿真系統(tǒng)流程圖；

圖2為骨折區(qū)域幾何模型建立流程圖；

圖3為骨折區(qū)域生物力學有限元分析流程圖；

圖4為骨痂單元組織分化模糊控制示意圖；

圖5為組織分化與骨痂單元力學刺激之間的關系示意圖；

圖6為骨痂單元力學刺激隸屬度函數(shù)；

圖7為軟骨組織含量隸屬度函數(shù)；

圖8為骨組織含量隸屬函數(shù)；

圖9為血供隸屬度函數(shù)。

具體實施方式

具體實施方式一：如圖1所示，本實施方式所述的一種基于組織分化的長骨骨折愈合仿真系統(tǒng)包括：

骨折區(qū)域幾何建模模塊1、骨折區(qū)域生物力學有限元分析模塊2、骨痂單元組織分化模塊3、程序終止判斷模塊4；

骨折區(qū)域幾何建模模塊1用于根據(jù)導入的二維斷層掃描圖像數(shù)據(jù)，經(jīng)過圖像預處理后進行骨折部位的三維表面幾何模型的建立；

骨折區(qū)域生物力學有限元分析模塊2用于對建立好的骨折區(qū)域模型進行網(wǎng)格劃分，施加外部載荷和設置邊界條件；

骨折區(qū)域生物力學有限元分析模塊2還用于初始骨折區(qū)域環(huán)境的設置；

骨折區(qū)域生物力學有限元分析模塊2還用于計算單元力學刺激；

骨痂單元組織分化模塊3用于對單元內(nèi)組織分化進行仿真，使單元內(nèi)各組織含量得到更新，從而使單元材料屬性得到更新，從而得到下一迭代步中所需要的單元力學刺激；

程序終止判斷模塊4用于判斷程序是否終止，若不滿足終止條件，程序進行下一迭代步；若滿足程序終止條件，則程序結束并輸出愈合時間。

具體實施方式二：如圖1～9所示，本實施方式中，所述的骨折區(qū)域幾何建模模塊1實現(xiàn)其功能的具體過程為：

采用基于分割的三維醫(yī)學影像表面重建算法對圖像進行三維表面重構，通過閾值篩選、交互式分割和三維重建過程得到三維表面幾何模型；

所述的影像由影像設備CT得到，數(shù)據(jù)存儲格式為DICOM。

其他組成及連接與具體實施方式一相同。

具體實施方式三：如圖1～9所示，本實施方式中，所述的骨折區(qū)域生物力學有限元分析模塊2實現(xiàn)其功能的具體過程為：

1)將骨折區(qū)域三維表面幾何模型進行網(wǎng)格化分，使連續(xù)的幾何模型離散化，得到骨折區(qū)域有限元模型；

所述的網(wǎng)格劃分包括面網(wǎng)格劃分和體網(wǎng)格劃分兩個步驟；面網(wǎng)格劃分過程用于將三維表面模型進行優(yōu)化，包括：表面模型優(yōu)化，平滑處理，修補漏洞；表面模型的優(yōu)化通過減小表面模型的三角面片來實現(xiàn)，該過程只需將相鄰的兩個頂點合并到一個新的頂點上，并延續(xù)原有的拓撲關系；平滑處理的過程中，對三維的面網(wǎng)格模型進行去噪；修補漏洞的過程中，通過將模型當中的空洞提取成空間多邊形，然后對空洞多邊形進行三角化的方法實現(xiàn)；體網(wǎng)格劃分的過程是將面網(wǎng)格模型進行拉伸、旋轉步驟來實現(xiàn)的；

通過網(wǎng)格劃分得到的骨折區(qū)域有限元模型包括單元編號和節(jié)點坐標兩部分；

節(jié)點坐標包含三列數(shù)據(jù)，三列數(shù)據(jù)分別代表每個節(jié)點的空間坐標值；

單元編號包含四列數(shù)據(jù)，四列數(shù)據(jù)分別為每個單元的四個節(jié)點的節(jié)點序號。

2)在有限元模型上施加外加載荷，并設置邊界條件。載荷的大小由骨所承受的力的大小決定，實驗對象不同，所受的力也不同；

3)對骨折區(qū)域有限元模型進行骨折區(qū)域初始環(huán)境設置。骨折區(qū)域由皮質(zhì)骨和骨痂區(qū)域兩部分組成。初始骨折區(qū)域環(huán)境設置包括皮質(zhì)骨材料屬性賦值，皮質(zhì)骨血供賦值，初始骨痂材料屬性賦值，初始骨痂區(qū)域血供賦值；

骨折之后，骨折斷端處的皮質(zhì)骨收到損害，血供也遭到破壞，所以將距離骨折斷端5mm內(nèi)的皮質(zhì)骨血供設置為0％，其余部分皮質(zhì)骨結構完整，血供條件良好，所以將其余部分的皮質(zhì)骨血供設置為100％；

骨痂外周區(qū)域可以接受其周圍組織提供的血供，所以將骨痂外周3mm內(nèi)血供設置為30％，骨痂內(nèi)部血供設置為0％。

4)將骨折區(qū)域看作雙相多孔彈性模型，由多孔彈性理論得到骨痂單元的本構方程，平衡方程和幾何方程，并通過有限單元法計算骨痂單元應力刺激S，具體過程為：

a.本構方程

式中，σ_rr，σ_θθ，σ_zz為正應力，τ_rθ，τ_θz，τ_rz為剪應力；ε_rr，ε_θθ，ε_zz為正應變，γ_rθ，γ_θz，γ_rz為主應變；α，α'分別為各向同性彈性面的Biot系數(shù)和軸向Biot系數(shù)；p為骨痂單元中的流體壓力；M₁₁，M₁₂，M₁₃，M₃₃，M₄₄，M₅₅為脫水的彈性模量矩陣分量；

其中，M₁₁，M₁₂，M₁₃，M₃₃，M₄₄，M₅₅脫水的彈性模量矩陣分量表達式如下所示：

M₄₄＝E_r/2(1+ν_r) (6)

M₅₅＝G' (7)

式中，E_r，ν_r分別是各同性彈性層的彈性模量和泊松比；E_z，ν_z分別是軸向彈性模量和泊松比；G'為剪切模量；

b.平衡方程

式中，σ_rr，σ_θθ，σ_zz為正應力，τ_rθ，τ_θz，τ_rz為剪應力；r為徑向半徑；

c.幾何方程

式中，ε_rr，ε_θθ，ε_zz為正應變，γ_rθ，γ_θz，γ_rz為主應變；u_r，u_θ，u_z分別為三個方向上的位移；r為徑向半徑；

通過上述方程的求解得到骨痂單元的正應變σ_rr，σ_θθ，σ_zz，由正應變可得到骨痂單元受到的畸變應變：

式中，D為骨痂單元受到的畸變應變；σ_rr，σ_θθ，σ_zz分別為各個方向上的正應變；

骨痂單元中液體的流速V為：

其中，k為骨痂中液體的達西滲透系數(shù)；u液體粘度；p為液體壓力；

由此可得到骨痂單元所受的力學刺激S為：

其中，D為骨痂單元受到的畸變應變；V為骨痂單元中液體流速；a，b分別為經(jīng)驗常數(shù)。

其他組成及連接與具體實施方式一或二之一相同。

具體實施方式四：如圖1～9所示，本實施方式中，所述的骨痂單元組織分化模塊3實現(xiàn)其功能的具體過程為：

1)設置輸入變量隸屬度函數(shù)

骨痂單元組織分化過程由六個輸入變量協(xié)同決定，分別為單元力學刺激，骨痂單元骨含量，骨痂單元軟骨含量，骨痂單元血供，周圍骨痂單元骨含量，周圍骨痂單元血供。輸入變量的精確值通過隸屬度函數(shù)轉變?yōu)橄鄳斎胱兞康碾`屬度。輸入變量隸屬度函數(shù)由梯形函數(shù)表示；

2)輸入隸屬度函數(shù)模糊化

將輸入變量的隸屬度函數(shù)用模糊邏輯語言值來表示：

單元力學刺激可分為五個等級，分別為：低，中，高；

單元血供可分為三個等級，分別為：低，中，高；

單元骨含量可分為三個等級，分別為：低，中，高；

單元軟骨含量可分為三個等級，分別為：低，中，高。

3)設置輸出變量隸屬度函數(shù)

輸出變量共有三個，分別為骨痂單元血供改變量，骨痂單元骨含量改變量和骨痂單元軟骨含量改變量。設置輸出變量隸屬度函數(shù)，輸出隸屬度函數(shù)由梯形函數(shù)表示，且將輸出隸屬度函數(shù)用如下的模糊邏輯語言值表示：

骨痂單元血供改變量分為三個等級，分別為：低，中，高；

骨痂單元骨含量改變量分為三個等級，分別為：低，中，高；

骨痂單元軟骨含量改變量分為三個等級，分別為：低，中，高。

4)模糊規(guī)則的編寫

骨痂內(nèi)組織分化主要包括：血管再生，膜內(nèi)骨化，軟骨生成，軟骨鈣化，軟骨骨化五個過程。通過模糊規(guī)則分別表示這五個組織分化過程，由此可得到組織分化后，骨痂單元中骨的模糊邏輯語言值，軟骨的模糊邏輯語言值以及血供的模糊邏輯語言值。具體過程如下：

過程一、血管再生

如果單元應力刺激為非高且骨痂單元血供為低且周圍骨痂單元血供為非低，那么骨痂單元血供增加；

如果單元應力刺激為非高且骨痂單元血供為中且周圍骨痂單元血供為高，那么骨痂單元血供增加；

如果單元應力刺激為非高且骨痂單元血供為高，那么骨痂單元血供增加；

過程二、膜內(nèi)骨化

如果骨痂單元軟骨含量為低且單元應力刺激為中且骨痂單元血供為高且周圍骨痂單元血供為高，那么骨痂單元骨含量增加；

過程三、軟骨生成

如果骨痂單元骨含量為非高且骨痂單元軟骨含量為低且單元應力刺激為中，那么骨痂單元軟骨含量增加；

如果骨痂單元軟骨含量為非低且單元應力刺激為高，那么骨痂單元軟骨含量增加；

如果骨痂單元軟骨含量為高且單元應力刺激為中，那么骨痂單元軟骨含量增加；

過程四、軟骨鈣化

如果骨痂單元軟骨含量為非低且單元應力刺激為中且骨痂單元血供為非低且周圍骨痂單元骨含量為非低，那么骨痂單元骨含量增加，骨痂單元軟骨含量減少；

如果骨痂單元軟骨含量為非低且單元應力刺激為高且骨痂單元血供為非低且周圍骨痂單元骨含量為非低，那么骨痂單元骨含量增加，骨痂單元軟骨含量減少；

過程五、軟骨骨化

如果骨痂單元骨含量為高且骨痂單元軟骨含量為低且單元應力刺激為中且骨痂單元血供為非低且周圍骨痂單元骨含量為高，那么骨痂單元骨含量增加，骨痂單元軟骨含量減少；

如果骨痂單元骨含量為高且骨痂單元軟骨含量為低且單元應力刺激為低且骨痂單元血供為非低且周圍骨痂單元骨含量為高，那么骨痂單元骨含量增加，骨痂單元軟骨含量減少；

5)輸出變量反模糊化

通過模糊規(guī)則的作用，得到輸出變量的模糊邏輯語言值，利用面積中心法得到輸出變量的隸屬度；

6)更新輸出變量含量

將得到的輸出變量隸屬度經(jīng)過轉化得到輸出變量的改變量，從而得到經(jīng)過組織分化后，骨痂單元骨組織，軟骨組織，結締組織含量以及骨痂單元血供。具體過程如下：

a.計算經(jīng)組織分化后骨痂單元骨組織含量

式中，為第n個骨痂單元在第k+1迭代步中骨組織的含量；為第n個骨痂單元在第k迭代步中骨的含量；為單位時間第n個骨痂單元在第k迭代步中骨的生成量；Δt為時間步長；

其中單位時間第n個骨痂單元在第k迭代步中骨的改變量由下式得到：

式中，為單位時間第n個骨痂單元在第k迭代步中骨的生成量；為骨痂單元中骨含量的隸屬度；T_Bone在骨痂單元骨的轉化率；

b.計算經(jīng)組織分化后骨痂單元軟骨骨組織含量

式中，為第n個骨痂單元在第k+1迭代步中軟骨的含量；為第n個骨痂單元在第k迭代步中軟骨的含量；為單位時間第n個骨痂單元在第k迭代步中軟骨的生成量；Δt為時間步長；

其中單位時間第n個骨痂單元在第k迭代步中軟骨的改變量由下式得到：

式中，為單位時間第n個骨痂單元在第k迭代步中軟骨的生成量；為骨痂單元中軟骨含量的隸屬度；T_Cartilage在骨痂單元軟骨的轉化率；

c.計算經(jīng)組織分化后骨痂單元血供

式中，為第n個骨痂單元在第k+1迭代步中結締組織的含量；為第n個骨痂單元在第k迭代步中結締組織的含量；為單位時間第n個骨痂單元在第k迭代步中結締組織的生成量；Δt為時間步長；

其中單位時間第n個骨痂單元在第k迭代步中軟骨的改變量由下式得到：

式中，為單位時間第n個骨痂單元在第k迭代步中結締組織的生成量；為骨痂單元中結締組織含量的隸屬度；T_Perfusio_n在骨痂單元結締組織的轉化率；

d.計算經(jīng)組織分化后骨痂單元結締組織含量

骨痂單元由骨組織，軟骨組織，結締組織三部分組成，三者之間的關系如下所示：

μ_Bone+μ_Cartilage+μ_ConnTissue＝1 (26)

式中，μ_Bone為骨痂單元骨組織含量；μ_Cartilage為骨痂單元軟骨組織含量；μ_ConnTissue為骨痂單元結締組織含量；

由此可得組織分化后骨痂單元結締組織含量。

其他組成及連接關系與具體實施方式一至三之一相同。

具體實施方式五：如圖1～9所示，本實施方式中，所述的程序終止判斷模塊4實現(xiàn)其功能的具體過程為：

1)判斷骨痂單元材料屬性

判斷當前骨痂單元材料屬性是否與骨的材料屬性相同，若不同，則程序執(zhí)行下面的2)步驟；若不等，程序執(zhí)行3)步驟；

2)更新骨痂單元材料屬性

若骨痂單元材料屬性不等于骨的材料屬性，則骨痂單元材料屬性進行更新，進入下一個迭代步，骨痂材料屬性更新公式如下：

式中，為第n個骨痂單元在第k+1迭代步中的彈性模量；E_Bone為骨的彈性模量；為第n個骨痂單元在第k+1迭代步中骨的含量；E_Cartilage為軟骨的彈性模量；為第n個骨痂單元在k+1迭代步中的軟骨的含量；E_ConnTissue為結締組織的彈性模量；為第n個骨痂單元在k+1迭代步中結締組織的含量；

式中，為第n個骨痂單元在第k+1迭代步中的泊松比；ν_Bone為骨的泊松比；為第n個骨痂單元在第k+1迭代步中骨的含量；ν_Cartilage為軟骨的泊松比；為第n個骨痂單元在k+1迭代步中的軟骨的含量；ν_ConnTissue為結締組織的泊松比；為第n個骨痂單元在k+1迭代步中結締組織的含量；

3)程序結束

若所有骨痂單元材料屬性等于骨的材料屬性，程序終止并輸出愈合時間。

其他組成及連接關系與具體實施方式一至四之一相同。

實施例：

為了說明本系統(tǒng)的使用方法，下面具體舉一個例子闡述本發(fā)明的操作過程。

模擬羊跖骨骨折愈合過程

1.羊跖骨骨折區(qū)域有限元模型的建立

1)建立幾何模型

把CT圖像數(shù)據(jù)利用三維醫(yī)學圖像表面重建算法對圖像進行三維表面重構，得到羊跖骨骨折區(qū)域三維幾何模型。

2)網(wǎng)格劃分

上述的三維幾何模型導入網(wǎng)格劃分軟件中進行網(wǎng)格化分，將得到的有限元模型導入到MTLAB中進行預處理，只提取目標數(shù)據(jù)，根據(jù)目標數(shù)據(jù)生成后續(xù)有限元計算所需要的單元編號和節(jié)點坐標文件。單元編號和節(jié)點坐標兩個文件為txt文本格式的文件。單元編號文件包含四列數(shù)據(jù)，四列數(shù)據(jù)分別為每個單元的四個節(jié)點序號，節(jié)點坐標文件包含散三列數(shù)據(jù)，三列數(shù)據(jù)分別為每個節(jié)點的空間坐標值。

2.骨折區(qū)域初始環(huán)境設置

骨折區(qū)域初始環(huán)境設置包括皮質(zhì)骨材料屬性賦值，初始骨痂單元材料屬性賦值。將骨折區(qū)域看成是基于混合物的雙相多孔彈性模型，皮質(zhì)骨和骨痂單元材料屬性如表1所示。

骨折區(qū)域初始環(huán)境設置還包括皮質(zhì)骨血供賦值和骨痂區(qū)域血供賦值。骨折之后，骨折斷端處的皮質(zhì)骨收到損害，血供也遭到破壞，所以將距離骨折斷端5mm內(nèi)的皮質(zhì)骨血供設置為0％，其余部分皮質(zhì)骨結構完整，血供條件良好，所以將其余部分的皮質(zhì)骨血供設置為100％；骨痂外周區(qū)域可以接受其周圍組織提供的血供，所以將骨痂外周3mm內(nèi)血供設置為30％，骨痂內(nèi)部血供設置為0％。

表1各組織材料屬性

3.骨折區(qū)域生物力學有限元分析

對建立好的有限元模型施加邊界條件和載荷。設定骨折區(qū)域下端固定約束長度，即將羊跖骨下端約束的那段所有節(jié)點的自由度均賦值為0，包括3個方向的位移和3個方向的旋轉。由于羊在正常行走時，跖骨所受力的大小為500N，所以在骨折區(qū)域施加大小為500N的邊界載荷。與求解羊跖骨相關的經(jīng)驗常數(shù)a和b分別為a＝0.0375μm/s，b＝3μm/s。

接下來對有限元模型進行生物力學有限元分析，得到骨折區(qū)域單元的畸變應變和流體速度，從而得到骨折區(qū)域單元收到的單元力學刺激。

4.組織分化過程

將由上一步得到的單元力學刺激，以及骨痂單元骨組織含量、周圍骨痂單元骨組織含量、骨痂單元軟骨組織含量、骨痂單元血供和周圍骨痂單元血供作為輸入變量，導入到由模糊規(guī)則中，得到經(jīng)組織分化后骨痂單元骨組織含量、骨痂單元軟骨組織含量和骨痂單元血供，如圖2所示。其中，單元力學刺激和組織分化之間的關系如圖3所示，骨組織、軟骨組織和結締組織之間的組織分化關系如圖4所示，輸入變量及輸出變量的隸屬度函數(shù)如圖5、6、7、8所示。

5.判斷程序終止條件

判斷當前骨痂單元材料屬性是否等于骨的材料屬性，若不等，則利用更新當前骨痂單元材料屬性，程序進入下一個迭代步；若相等，則程序結束，并輸出骨折愈合時間。