一種基于組織分化的長骨骨折愈合仿真系統的制作方法
本發明涉及生物醫學工程領域，特別涉及一種基于組織分化的長骨骨折愈合仿真系統。
背景技術：
骨折是一種常見的創傷，骨折的高發性使得對骨折機理及促進愈合的研究尤為迫切，一旦骨折發生，與其它組織損傷修復不同的是，骨折不是靠纖維結締組織連接，而是骨組織的完全再生。然而，并不是所有的骨折都可以完成愈合，有時會發生延遲愈合甚至是不愈合。骨折延遲愈合或者不愈合會引起患肢疼痛，功能障礙，導致患者失業，由此造成很大的社會經濟負擔。因此，盡管關于骨折愈合的研究一直備受關注，但仍有5％～10％的骨折因各種原因發生延遲愈合甚至是不愈合。
骨折愈合受到特定的幾何因素、力學因素、生物學因素影響，良好的幾何因素、力學因素、生物學因素得到良好愈合效果。反之則會導致骨折的延遲愈合甚至是不愈合。目前缺少能夠精確表達骨折愈合這一復雜過程的計算機仿真系統。在骨折愈合仿真系統中存在以下缺陷：
1.沒有建立專門針對患者的個體化模型；
2.力學因素與骨折愈合過程沒有一個確定性關系；
3.骨折愈合區域的模型和生物力學材料設置過于簡化；
4.沒有在同一個仿真系統中綜合考慮力學因素和生物學因素對骨折愈合的影響。
技術實現要素：
本發明的目的是為了解決現有的骨折愈合仿真中不能綜合模擬力學因素和生物學因素對骨折愈合過程的影響，骨折愈合生物力學模型材料設置過于簡化的缺點，而提出的一種基于組織分化的長骨骨折愈合仿真系統。
本發明的目的通過下述技術方案實現：一種基于組織分化的長骨骨折愈合仿真系統，其特征在于，所述系統包括：
骨折區域幾何建模模塊、骨折區域生物力學有限元分析模塊、骨痂單元組織分化模塊和程序終止判斷模塊；
骨折區域幾何建模模塊用于根據導入的二維斷層掃描圖像數據，經過圖像預處理后進行骨折部位的三維表面幾何模型的建立；
骨折區域生物力學有限元分析模塊用于對建立好的骨折區域模型進行網格劃分，施加外部載荷和設置邊界條件；
骨折區域生物力學有限元分析模塊還用于初始骨折區域環境的設置；
骨折區域生物力學有限元分析模塊還用于計算單元力學刺激；
骨痂單元組織分化模塊用于對單元內組織分化進行仿真，使單元內各組織含量得到更新，從而使單元材料屬性得到更新，從而得到下一迭代步中所需要的單元力學刺激；
程序終止判斷模塊用于判斷程序是否終止，若不滿足終止條件，程序進行下一迭代步；若滿足程序終止條件，則程序結束并輸出愈合時間。
本發明的有益效果為：
1.本發明提出的一種基于組織分化的長骨骨折愈合仿真系統是基于windows開發語言平臺來開發軟件，通過自主編程實現骨折愈合過程的動態模擬，基于對話框的形式，易于操作，培訓周期短；
2.將骨折區域設置為雙相多孔彈性模型，相比單相模型，更加符合骨折愈合區域的生物特性，使仿真結果更加精確；
3.將血供作為變量引入骨折愈合仿真系統中，能夠更加準確的模擬骨折愈合的過程以及血液對骨折愈合的影響；
4.利用模糊邏輯的方法對組織分化進行仿真，相較于傳統的采用偏微分方程組進行建模的方法，更加便于理解。減少了在建模過程中偏微分方程組的數量，便于建模且減少了仿真時間；
5.通過構建骨折愈合仿真系統，可以對醫生制定最優的手術方案提供指導，進而提高手術成功率、提高骨折愈質量，減少骨折不愈合和延遲愈合的情況；
6.通過構建骨折愈合仿真系統，可以對建立的仿真模型進行多次重復實驗研究，減少真實的生物實驗，節省時間，提高效率，節省費用，避免人道主義的爭議。
綜上，本發明的仿真平臺克服了現有技術的缺點與不足。
附圖說明
圖1為基于組織分化的骨折愈合仿真系統流程圖；
圖2為骨折區域幾何模型建立流程圖；
圖3為骨折區域生物力學有限元分析流程圖；
圖4為骨痂單元組織分化模糊控制示意圖；
圖5為組織分化與骨痂單元力學刺激之間的關系示意圖；
圖6為骨痂單元力學刺激隸屬度函數；
圖7為軟骨組織含量隸屬度函數；
圖8為骨組織含量隸屬函數；
圖9為血供隸屬度函數。
具體實施方式
具體實施方式一：如圖1所示，本實施方式所述的一種基于組織分化的長骨骨折愈合仿真系統包括：
骨折區域幾何建模模塊1、骨折區域生物力學有限元分析模塊2、骨痂單元組織分化模塊3、程序終止判斷模塊4；
骨折區域幾何建模模塊1用于根據導入的二維斷層掃描圖像數據，經過圖像預處理后進行骨折部位的三維表面幾何模型的建立；
骨折區域生物力學有限元分析模塊2用于對建立好的骨折區域模型進行網格劃分，施加外部載荷和設置邊界條件；
骨折區域生物力學有限元分析模塊2還用于初始骨折區域環境的設置；
骨折區域生物力學有限元分析模塊2還用于計算單元力學刺激；
骨痂單元組織分化模塊3用于對單元內組織分化進行仿真，使單元內各組織含量得到更新，從而使單元材料屬性得到更新，從而得到下一迭代步中所需要的單元力學刺激；
程序終止判斷模塊4用于判斷程序是否終止，若不滿足終止條件，程序進行下一迭代步；若滿足程序終止條件，則程序結束并輸出愈合時間。
具體實施方式二：如圖1～9所示，本實施方式中，所述的骨折區域幾何建模模塊1實現其功能的具體過程為：
采用基于分割的三維醫學影像表面重建算法對圖像進行三維表面重構，通過閾值篩選、交互式分割和三維重建過程得到三維表面幾何模型；
所述的影像由影像設備CT得到，數據存儲格式為DICOM。
其他組成及連接與具體實施方式一相同。
具體實施方式三：如圖1～9所示，本實施方式中，所述的骨折區域生物力學有限元分析模塊2實現其功能的具體過程為：
1)將骨折區域三維表面幾何模型進行網格化分，使連續的幾何模型離散化，得到骨折區域有限元模型；
所述的網格劃分包括面網格劃分和體網格劃分兩個步驟；面網格劃分過程用于將三維表面模型進行優化，包括：表面模型優化，平滑處理，修補漏洞；表面模型的優化通過減小表面模型的三角面片來實現，該過程只需將相鄰的兩個頂點合并到一個新的頂點上，并延續原有的拓撲關系；平滑處理的過程中，對三維的面網格模型進行去噪；修補漏洞的過程中，通過將模型當中的空洞提取成空間多邊形，然后對空洞多邊形進行三角化的方法實現；體網格劃分的過程是將面網格模型進行拉伸、旋轉步驟來實現的；
通過網格劃分得到的骨折區域有限元模型包括單元編號和節點坐標兩部分；
節點坐標包含三列數據，三列數據分別代表每個節點的空間坐標值；
單元編號包含四列數據，四列數據分別為每個單元的四個節點的節點序號。
2)在有限元模型上施加外加載荷，并設置邊界條件。載荷的大小由骨所承受的力的大小決定，實驗對象不同，所受的力也不同；
3)對骨折區域有限元模型進行骨折區域初始環境設置。骨折區域由皮質骨和骨痂區域兩部分組成。初始骨折區域環境設置包括皮質骨材料屬性賦值，皮質骨血供賦值，初始骨痂材料屬性賦值，初始骨痂區域血供賦值；
骨折之后，骨折斷端處的皮質骨收到損害，血供也遭到破壞，所以將距離骨折斷端5mm內的皮質骨血供設置為0％，其余部分皮質骨結構完整，血供條件良好，所以將其余部分的皮質骨血供設置為100％；
骨痂外周區域可以接受其周圍組織提供的血供，所以將骨痂外周3mm內血供設置為30％，骨痂內部血供設置為0％。
4)將骨折區域看作雙相多孔彈性模型，由多孔彈性理論得到骨痂單元的本構方程，平衡方程和幾何方程，并通過有限單元法計算骨痂單元應力刺激S，具體過程為：
a.本構方程
式中，σrr，σθθ，σzz為正應力，τrθ，τθz，τrz為剪應力；εrr，εθθ，εzz為正應變，γrθ，γθz，γrz為主應變；α，α'分別為各向同性彈性面的Biot系數和軸向Biot系數；p為骨痂單元中的流體壓力；M11，M12，M13，M33，M44，M55為脫水的彈性模量矩陣分量；
其中，M11，M12，M13，M33，M44，M55脫水的彈性模量矩陣分量表達式如下所示：
M44＝Er/2(1+νr) (6)
M55＝G' (7)
式中，Er，νr分別是各同性彈性層的彈性模量和泊松比；Ez，νz分別是軸向彈性模量和泊松比；G'為剪切模量；
b.平衡方程
式中，σrr，σθθ，σzz為正應力，τrθ，τθz，τrz為剪應力；r為徑向半徑；
c.幾何方程
式中，εrr，εθθ，εzz為正應變，γrθ，γθz，γrz為主應變；ur，uθ，uz分別為三個方向上的位移；r為徑向半徑；
通過上述方程的求解得到骨痂單元的正應變σrr，σθθ，σzz，由正應變可得到骨痂單元受到的畸變應變：
式中，D為骨痂單元受到的畸變應變；σrr，σθθ，σzz分別為各個方向上的正應變；
骨痂單元中液體的流速V為：
其中，k為骨痂中液體的達西滲透系數；u液體粘度；p為液體壓力；
由此可得到骨痂單元所受的力學刺激S為：
其中，D為骨痂單元受到的畸變應變；V為骨痂單元中液體流速；a，b分別為經驗常數。
其他組成及連接與具體實施方式一或二之一相同。
具體實施方式四：如圖1～9所示，本實施方式中，所述的骨痂單元組織分化模塊3實現其功能的具體過程為：
1)設置輸入變量隸屬度函數
骨痂單元組織分化過程由六個輸入變量協同決定，分別為單元力學刺激，骨痂單元骨含量，骨痂單元軟骨含量，骨痂單元血供，周圍骨痂單元骨含量，周圍骨痂單元血供。輸入變量的精確值通過隸屬度函數轉變為相應輸入變量的隸屬度。輸入變量隸屬度函數由梯形函數表示；
2)輸入隸屬度函數模糊化
將輸入變量的隸屬度函數用模糊邏輯語言值來表示：
單元力學刺激可分為五個等級，分別為：低，中，高；
單元血供可分為三個等級，分別為：低，中，高；
單元骨含量可分為三個等級，分別為：低，中，高；
單元軟骨含量可分為三個等級，分別為：低，中，高。
3)設置輸出變量隸屬度函數
輸出變量共有三個，分別為骨痂單元血供改變量，骨痂單元骨含量改變量和骨痂單元軟骨含量改變量。設置輸出變量隸屬度函數，輸出隸屬度函數由梯形函數表示，且將輸出隸屬度函數用如下的模糊邏輯語言值表示：
骨痂單元血供改變量分為三個等級，分別為：低，中，高；
骨痂單元骨含量改變量分為三個等級，分別為：低，中，高；
骨痂單元軟骨含量改變量分為三個等級，分別為：低，中，高。
4)模糊規則的編寫
骨痂內組織分化主要包括：血管再生，膜內骨化，軟骨生成，軟骨鈣化，軟骨骨化五個過程。通過模糊規則分別表示這五個組織分化過程，由此可得到組織分化后，骨痂單元中骨的模糊邏輯語言值，軟骨的模糊邏輯語言值以及血供的模糊邏輯語言值。具體過程如下：
過程一、血管再生
如果單元應力刺激為非高且骨痂單元血供為低且周圍骨痂單元血供為非低，那么骨痂單元血供增加；
如果單元應力刺激為非高且骨痂單元血供為中且周圍骨痂單元血供為高，那么骨痂單元血供增加；
如果單元應力刺激為非高且骨痂單元血供為高，那么骨痂單元血供增加；
過程二、膜內骨化
如果骨痂單元軟骨含量為低且單元應力刺激為中且骨痂單元血供為高且周圍骨痂單元血供為高，那么骨痂單元骨含量增加；
過程三、軟骨生成
如果骨痂單元骨含量為非高且骨痂單元軟骨含量為低且單元應力刺激為中，那么骨痂單元軟骨含量增加；
如果骨痂單元軟骨含量為非低且單元應力刺激為高，那么骨痂單元軟骨含量增加；
如果骨痂單元軟骨含量為高且單元應力刺激為中，那么骨痂單元軟骨含量增加；
過程四、軟骨鈣化
如果骨痂單元軟骨含量為非低且單元應力刺激為中且骨痂單元血供為非低且周圍骨痂單元骨含量為非低，那么骨痂單元骨含量增加，骨痂單元軟骨含量減少；
如果骨痂單元軟骨含量為非低且單元應力刺激為高且骨痂單元血供為非低且周圍骨痂單元骨含量為非低，那么骨痂單元骨含量增加，骨痂單元軟骨含量減少；
過程五、軟骨骨化
如果骨痂單元骨含量為高且骨痂單元軟骨含量為低且單元應力刺激為中且骨痂單元血供為非低且周圍骨痂單元骨含量為高，那么骨痂單元骨含量增加，骨痂單元軟骨含量減少；
如果骨痂單元骨含量為高且骨痂單元軟骨含量為低且單元應力刺激為低且骨痂單元血供為非低且周圍骨痂單元骨含量為高，那么骨痂單元骨含量增加，骨痂單元軟骨含量減少；
5)輸出變量反模糊化
通過模糊規則的作用，得到輸出變量的模糊邏輯語言值，利用面積中心法得到輸出變量的隸屬度；
6)更新輸出變量含量
將得到的輸出變量隸屬度經過轉化得到輸出變量的改變量，從而得到經過組織分化后，骨痂單元骨組織，軟骨組織，結締組織含量以及骨痂單元血供。具體過程如下：
a.計算經組織分化后骨痂單元骨組織含量
式中，為第n個骨痂單元在第k+1迭代步中骨組織的含量；為第n個骨痂單元在第k迭代步中骨的含量；為單位時間第n個骨痂單元在第k迭代步中骨的生成量；Δt為時間步長；
其中單位時間第n個骨痂單元在第k迭代步中骨的改變量由下式得到：
式中，為單位時間第n個骨痂單元在第k迭代步中骨的生成量；為骨痂單元中骨含量的隸屬度；TBone在骨痂單元骨的轉化率；
b.計算經組織分化后骨痂單元軟骨骨組織含量
式中，為第n個骨痂單元在第k+1迭代步中軟骨的含量；為第n個骨痂單元在第k迭代步中軟骨的含量；為單位時間第n個骨痂單元在第k迭代步中軟骨的生成量；Δt為時間步長；
其中單位時間第n個骨痂單元在第k迭代步中軟骨的改變量由下式得到：
式中，為單位時間第n個骨痂單元在第k迭代步中軟骨的生成量；為骨痂單元中軟骨含量的隸屬度；TCartilage在骨痂單元軟骨的轉化率；
c.計算經組織分化后骨痂單元血供
式中，為第n個骨痂單元在第k+1迭代步中結締組織的含量；為第n個骨痂單元在第k迭代步中結締組織的含量；為單位時間第n個骨痂單元在第k迭代步中結締組織的生成量；Δt為時間步長；
其中單位時間第n個骨痂單元在第k迭代步中軟骨的改變量由下式得到：
式中，為單位時間第n個骨痂單元在第k迭代步中結締組織的生成量；為骨痂單元中結締組織含量的隸屬度；TPerfusion在骨痂單元結締組織的轉化率；
d.計算經組織分化后骨痂單元結締組織含量
骨痂單元由骨組織，軟骨組織，結締組織三部分組成，三者之間的關系如下所示：
μBone+μCartilage+μConnTissue＝1 (26)
式中，μBone為骨痂單元骨組織含量；μCartilage為骨痂單元軟骨組織含量；μConnTissue為骨痂單元結締組織含量；
由此可得組織分化后骨痂單元結締組織含量。
其他組成及連接關系與具體實施方式一至三之一相同。
具體實施方式五：如圖1～9所示，本實施方式中，所述的程序終止判斷模塊4實現其功能的具體過程為：
1)判斷骨痂單元材料屬性
判斷當前骨痂單元材料屬性是否與骨的材料屬性相同，若不同，則程序執行下面的2)步驟；若不等，程序執行3)步驟；
2)更新骨痂單元材料屬性
若骨痂單元材料屬性不等于骨的材料屬性，則骨痂單元材料屬性進行更新，進入下一個迭代步，骨痂材料屬性更新公式如下：
式中，為第n個骨痂單元在第k+1迭代步中的彈性模量；EBone為骨的彈性模量；為第n個骨痂單元在第k+1迭代步中骨的含量；ECartilage為軟骨的彈性模量；為第n個骨痂單元在k+1迭代步中的軟骨的含量；EConnTissue為結締組織的彈性模量；為第n個骨痂單元在k+1迭代步中結締組織的含量；
式中，為第n個骨痂單元在第k+1迭代步中的泊松比；νBone為骨的泊松比；為第n個骨痂單元在第k+1迭代步中骨的含量；νCartilage為軟骨的泊松比；為第n個骨痂單元在k+1迭代步中的軟骨的含量；νConnTissue為結締組織的泊松比；為第n個骨痂單元在k+1迭代步中結締組織的含量；
3)程序結束
若所有骨痂單元材料屬性等于骨的材料屬性，程序終止并輸出愈合時間。
其他組成及連接關系與具體實施方式一至四之一相同。
實施例：
為了說明本系統的使用方法，下面具體舉一個例子闡述本發明的操作過程。
模擬羊跖骨骨折愈合過程
1.羊跖骨骨折區域有限元模型的建立
1)建立幾何模型
把CT圖像數據利用三維醫學圖像表面重建算法對圖像進行三維表面重構，得到羊跖骨骨折區域三維幾何模型。
2)網格劃分
上述的三維幾何模型導入網格劃分軟件中進行網格化分，將得到的有限元模型導入到MTLAB中進行預處理，只提取目標數據，根據目標數據生成后續有限元計算所需要的單元編號和節點坐標文件。單元編號和節點坐標兩個文件為txt文本格式的文件。單元編號文件包含四列數據，四列數據分別為每個單元的四個節點序號，節點坐標文件包含散三列數據，三列數據分別為每個節點的空間坐標值。
2.骨折區域初始環境設置
骨折區域初始環境設置包括皮質骨材料屬性賦值，初始骨痂單元材料屬性賦值。將骨折區域看成是基于混合物的雙相多孔彈性模型，皮質骨和骨痂單元材料屬性如表1所示。
骨折區域初始環境設置還包括皮質骨血供賦值和骨痂區域血供賦值。骨折之后，骨折斷端處的皮質骨收到損害，血供也遭到破壞，所以將距離骨折斷端5mm內的皮質骨血供設置為0％，其余部分皮質骨結構完整，血供條件良好，所以將其余部分的皮質骨血供設置為100％；骨痂外周區域可以接受其周圍組織提供的血供，所以將骨痂外周3mm內血供設置為30％，骨痂內部血供設置為0％。
表1各組織材料屬性
3.骨折區域生物力學有限元分析
對建立好的有限元模型施加邊界條件和載荷。設定骨折區域下端固定約束長度，即將羊跖骨下端約束的那段所有節點的自由度均賦值為0，包括3個方向的位移和3個方向的旋轉。由于羊在正常行走時，跖骨所受力的大小為500N，所以在骨折區域施加大小為500N的邊界載荷。與求解羊跖骨相關的經驗常數a和b分別為a＝0.0375μm/s，b＝3μm/s。
接下來對有限元模型進行生物力學有限元分析，得到骨折區域單元的畸變應變和流體速度，從而得到骨折區域單元收到的單元力學刺激。
4.組織分化過程
將由上一步得到的單元力學刺激，以及骨痂單元骨組織含量、周圍骨痂單元骨組織含量、骨痂單元軟骨組織含量、骨痂單元血供和周圍骨痂單元血供作為輸入變量，導入到由模糊規則中，得到經組織分化后骨痂單元骨組織含量、骨痂單元軟骨組織含量和骨痂單元血供，如圖2所示。其中，單元力學刺激和組織分化之間的關系如圖3所示，骨組織、軟骨組織和結締組織之間的組織分化關系如圖4所示，輸入變量及輸出變量的隸屬度函數如圖5、6、7、8所示。
5.判斷程序終止條件
判斷當前骨痂單元材料屬性是否等于骨的材料屬性，若不等，則利用更新當前骨痂單元材料屬性，程序進入下一個迭代步；若相等，則程序結束，并輸出骨折愈合時間。
- 還沒有人留言評論。精彩留言會獲得點贊！