李 洋 桑建兵 敖日汗 馬 鈺 魏新宇

(河北工業(yè)大學機械工程學院,天津 300130)

引言

骨骼肌約占生物體重的40%[1],是維持生物基本運動最重要的生理軟組織,而骨骼肌的損傷會導致人類行為嚴重受限,如常見的臨床病例:中風、截癱等[2].在對骨骼肌研究中,軟組織的本構關系可能是最基本和最關鍵的,因為其可以建立可靠的數(shù)學或計算模型.由于骨骼肌組織的特殊特性(各向異性[3]、非線性[4]、超彈性[5]、時間依賴性[6]) 的存在,骨骼肌組織的力學特性的表征具有很大的挑戰(zhàn)性.全面了解骨骼肌的組成行為,對于更好地了解骨骼肌的生物學特性,提高肌肉疾病和損傷的治療水平具有重要意義.許多學者進行了大量的實驗來研究骨骼肌的生物力學特性,包括壓縮[7]、拉伸[8]、應力松弛[9]等.王寶珍等[10]對豬后腿肌肉進行沖擊壓縮實驗,并利用改進后的分離式Hopkinson 壓桿技術獲取了不同加載方向和不同應變率下的應力應變數(shù)據(jù).Mohammadkhah 等[11]對新鮮雞胸大肌進行了加載方向平行于纖維方向和垂直于纖維方向以及加載方向與纖維方向成45°的單軸準靜態(tài)壓縮和拉伸實驗,證明雞肌肉的力學行為是非線性和各向異性的.B¨ol 等[12]使用3 種不同的加載模式對肌肉組織進行了一系列實驗.結果表明,纖維方向對壓縮應力有各向異性貢獻.如果有可靠的正問題模型,反問題分析是系統(tǒng)的識別材料性質的有效手段[13].Takaza 等[14]證明了反向有限元方法在確定材料性能方面的有效性.Silva 等[15]通過逆有限元方法并結合優(yōu)化算法獲得了女性恥骨聯(lián)合肌肉的主動和被動力學特性.雖然現(xiàn)有的研究在識別本構參數(shù)方面取得了進展,但現(xiàn)有預測的準確性還有待進一步提高.

隨著人工智能技術(AI)的發(fā)展,機器學習(ML)作為AI 的一個分支已被廣泛應用于生活中的各個領域,如數(shù)據(jù)融合和識別[16-17]、醫(yī)學檢測[18-19]、流體力學[20-22]和固體力學中[23-24].ML 在模擬輸入和輸出變量之間的非線性相關性方面具有很高的效率,并且已經(jīng)被用于預測生物軟組織的力學特性.Liang等[25]利用支持向量機模型和有限元方法相結合來預測主動脈瘤的風險,并取得了良好的效果.Lu 等[26]提出了一種基于動態(tài)多體方法和神經(jīng)網(wǎng)絡分析的創(chuàng)新方法,建立了基于有限元的軟骨神經(jīng)網(wǎng)絡模型.Liu等[27]對牛的肌肉組織進行了無約束壓縮實驗,并利用了改進后的小生境遺傳算法預測了肌肉組織的本構參數(shù).通常AI 的成功歸結于大規(guī)模的數(shù)據(jù)集,然而許多現(xiàn)實應用的場景例如醫(yī)學,軍事等領域,因為隱私或者安全等因素不允許收集足夠多的數(shù)據(jù)集,所以對于小樣本的學習便成了當今AI 領域的熱點問題.Koch 等[28]將深度學習與小樣本學習結合起來,提出了一種新的卷積神經(jīng)網(wǎng)絡模式來學習成對的樣本中與類別無關的變量.Liu 等[29]基于BP 神經(jīng)網(wǎng)絡,建立了一個可靠的正問題神經(jīng)網(wǎng)絡來代替有限元方法收集樣本,來縮短收集樣本的時間,也為解決小樣本學習問題提供了新的想法.而簡單的機器學習模型對小樣本學習也會有良好的效果.因此,本研究考慮將機器學習的方法應用于骨骼肌的本構參數(shù)的預測,分別提出了兩種反問題求解方法,一種是將K 近鄰模型(KNN)與有限元方法(FEM)相結合,另一種是將支持向量機回歸模型(SVR)與有限元方法相結合.在本研究中,通過訓練和測試Liu 等[27]收集的數(shù)據(jù)集,證明了所提出的兩種反問題求解方法的有效性.本研究還進一步對比了所提出的兩種方法,對比結果表明,利用K 近鄰模型結合有限元方法是預測骨骼肌本構參數(shù)的更精確有效的方法.

1 有限元模型的建立

骨骼肌的在生理學上的數(shù)學表述是一項復雜的任務,考慮到其大變形行為,骨骼肌的材料特性常被模擬為各向異性,非線性和超彈性.各向異性超彈性模型可以很好模擬骨骼肌的力學行為,其本構方程如下[30]

各向異性骨骼肌中本構參數(shù)的確定是一項具有挑戰(zhàn)性的任務,由于骨骼肌組織的高含水量,骨骼肌通常被認為是不可壓縮的物質[32].所以本文提出的反演程序預測的本構參數(shù)包括:C10,k1,k2和κ.

圖1 給出了為模擬骨骼肌單軸壓縮試驗的有限元模型.骨骼肌被建模為超彈性各向異性材料.有限元網(wǎng)格由3375 個八節(jié)點六面體單元和4096 個節(jié)點組成.八節(jié)點六面體單元相對于四面體單元,增加了節(jié)點的數(shù)目,可以顯著的提高計算精度.因為拉伸試驗機頂部和底部的壓板是由鋼制成的,相比于新鮮的骨骼肌組織要硬的多,所以在有限元仿真中他們被定義為剛體.在壓板與骨骼肌組織樣本之間設置了摩擦系數(shù)以更準確的模擬實驗.

圖1 模擬實驗設置的骨骼肌有限元模型Fig.1 Finite element model of skeletal muscles for simulating experiment

2 K 近鄰與支持向量機預測模型構建

KNN 是一種常用的有監(jiān)督學習的機器學習算法,因為其具有從高維特征空間中恢復相似實例的簡單性和直觀性,所以具有廣泛的應用前景.在KNN模型中,KNN 對于給定測試樣本,基于距離度量找出訓練集中與其最靠近的K個訓練樣本,然后基于這K個“鄰居”的信息來進行預測[33].KNN 模型的核心問題是K值的選取以及距離函數(shù)的選擇.在本研究中,通過大量的數(shù)值實驗,最終確定K=3.待求樣本的狀態(tài)向量與已知樣本的狀態(tài)向量之間的距離量度表示二者之間的相似程度.常見的KNN 模型的距離函數(shù)有歐式距離、曼哈頓距離、馬氏距離、切比雪夫距離等.在本研究中使用的歐氏距離公式如下:假設待求樣本的狀態(tài)向量為x=(x1,x2,···,xn),已知樣本的狀態(tài)向量為y=(y1,y2,···,yn),則其之間的距離d為

支持向量機(SVM)也是一種常用的機器學習算法,具有很強的非線性數(shù)據(jù)處理能力.而SVR 是支持向量機算法通過非線性回歸分析的延伸[34],是解決回歸問題的一種有效方法.支持向量機回歸模型利用核函數(shù)將非線性問題轉化為高維空間的線性問題,是替代高維空間復雜計算的一種非常有效的方法[35].常見的核函數(shù)有線性核函數(shù)、多項式核函數(shù)、徑向基核函數(shù)(radical basis function,RBF).其中應用最廣泛的是徑向基核函數(shù),其已廣泛應用在機械設備故障識別預測[36].RBF 定義為

這是一個從0 變化到1 的鐘形函數(shù).其中a和b代表樣本向量,γ 為超參數(shù).

3 預測結果與討論

在本研究中,Liu 等[27]對骨骼肌進行單軸壓縮實驗中獲得的名義應力和主伸長數(shù)據(jù)集被分成兩個子數(shù)據(jù)集,分別對應于加載方向平行于纖維方向和加載方向垂直于纖維方向.參數(shù)空間的必要參數(shù)為4 個,分別為C10(取值范圍0.2～1.5 kPa),k1(取值范圍0.05～0.4 MPa),k2(取值范圍10～60),κ(取值范圍0.04～0.32).本研究在參數(shù)空間中分別以固定增量采集了250 組骨骼肌超彈性本構參數(shù)以使得采集到的本構參數(shù)能夠充分代表整個參數(shù)空間,其中C10增量為5.0×10-3kPa,k1增量為1.4×10-3MPa,k2增量為0.2,κ 增量為1.12×10-3.將收集到的骨骼肌組織的本構參數(shù)通過建立的有限元模型得到了相應的名義應力和主伸長.之后以有限元模型收集的名義應力作為輸入,以相應的本構參數(shù)作為輸出分別構建KNN模型和SVR 模型.為了更準確,快速地進行回歸,本研究分別對名義應力和骨骼肌組織的本構參數(shù)進行歸一化處理

其中,X′是歸一化后的取值,X是實際取值;Xmax,Xmin分別對應于名義應力和骨骼肌組織的本構參數(shù)的最大值和最小值.當訓練完成后,將骨骼肌壓縮實驗中獲得的名義應力數(shù)據(jù)集分別輸入兩個模型,最終得到實驗樣本的本構參數(shù).

為了進行監(jiān)督學習,需要選擇合適的訓練和測試數(shù)據(jù)進行預處理.訓練集用于構建KNN 和SVR模型,測試集用于驗證模型的泛化能力.訓練集的大小必須代表整個數(shù)據(jù)集的特征.訓練集數(shù)目太少可能導致模型的學習能力差,而訓練集數(shù)目太多很可能導致驗證困難并導致計算時間過長.在本研究中,80%的有限元收集到的數(shù)據(jù)被分給訓練集,其余20%被分配給測試集[37].圖2 給出了本文中提出的兩個反向識別本構參數(shù)反演方法的流程圖.

圖2 本研究中提出的本構參數(shù)識別方法流程圖Fig.2 Flow chart of the parameter identification method proposed in this study

基于Liu 等[27]得到的骨骼肌壓縮實驗數(shù)據(jù),本研究對利用提出的參數(shù)識別程序得到的本構參數(shù)進行處理,以獲得一組本構參數(shù),這些參數(shù)既可以表示加載方向平行于纖維方向,又可以很好地表示加載方向垂直于纖維方向.在訓練集完成之后,測試集的數(shù)據(jù)用來評估KNN 模型和SVR 模型的準確性,并比較本構參數(shù)的真實值和預測值.在本研究中,選用了決定系數(shù)(R2)作為預測值準確率的指標.R2反映了因變量的所有變化可以通過回歸關系由自變量所解釋的比例,即利用數(shù)據(jù)的平均值作為誤差基準,觀察預測誤差與平均參考誤差的偏離程度.R2為

式中,N是樣本的數(shù)量,yi代表真實值,ypred代表預測值,代表真實值的平均值.圖3 描繪了測試集中本構參數(shù)的真實值與預測值準確率的直方圖.從直方圖中可以看出,無論是KNN 模型還是SVR 模型測試集中各本構參數(shù)的準確率都很好.

圖3 KNN 模型和SVR 模型測試集預測準確率比較Fig.3 Comparison of test set prediction accuracy between KNN model and SVR model

圖4 給出了當加載方向平行于纖維方向和加載方向垂直于纖維方向時,利用KNN 模型獲得的本構參數(shù)的計算響應曲線和實驗曲線的比較以及仿真加載過程的名義應力云圖.從圖4 可以看出名義應力和主伸長之間存在非線性關系,并且名義應力隨著主伸長的增大而增大,無論加載方向平行于纖維方向還是垂直于纖維方向,利用KNN 模型預測的骨骼肌本構參數(shù)的響應曲線與實驗曲線吻合良好.通過比較圖4(a) 和圖4(b) 可以發(fā)現(xiàn),當加載方向垂直于纖維方向時,骨骼肌組織對壓縮變形表現(xiàn)出更強的抵抗力.

圖4 利用KNN 模型獲得的計算響應以及仿真加載過程名義應力云圖Fig.4 Calculated response of KNN model and the nominal stress diagram of simulation loading process

圖5 給出了當加載方向平行于纖維方向和加載方向垂直于纖維方向時,利用SVR 模型獲得的本構參數(shù)的計算響應曲線與利用KNN 模型獲得的本構參數(shù)的計算響應曲線和實驗曲線的比較.從圖5 中可以看出,無論加載方向平行于纖維方向還是垂直于纖維方向,利用SVR 模型和KNN 模型預測的骨骼肌本構參數(shù)的響應曲線與實驗曲線吻合良好.但在加載方向垂直于纖維方向時,利用KNN 模型預測的骨骼肌本構參數(shù)的響應曲線比利用SVR 模型與實驗曲線吻合更好.

圖5 利用SVR 模型獲得的本構參數(shù)的計算響應與實驗數(shù)據(jù)以及與利用KNN 模型獲得的計算響應對比Fig.5 Comparison of calculated response between KNN model and SVR model

圖5 利用SVR 模型獲得的本構參數(shù)的計算響應與實驗數(shù)據(jù)以及與利用KNN 模型獲得的計算響應對比(續(xù))Fig.5 Comparison of calculated response between KNN model and SVR model(continued)

為了更直觀地觀察真實值與預測值的擬合程度,引入了2 種統(tǒng)計指標,包括相關系數(shù)(R)和決定系數(shù)(R2).R是一個相關性指標,可以用來衡量真實值與預測值之間的關系有多強,R2是一個效率指標如前所述.R為代表預測值的平均值.表1 給出了從統(tǒng)計指標方面比較了利用SVR 模型參數(shù)識別方法和利用KNN 模型參數(shù)識別方法.

從表1 可以更直觀看出,無論是加載方向平行與纖維方向還是垂直于纖維方向,KNN 模型和SVR模型的相關系數(shù)R和決定系數(shù)R2相差不大,在加載方向垂直于纖維方向時,KNN 模型的相關系數(shù)R和決定系數(shù)R2比SVR 模型好一些.所以KNN 模型和SVR 模型都是預測本研究中骨骼肌組織的本構參數(shù)的準確有效手段.本研究對利用KNN 模型得到的骨骼肌組織的本構參數(shù)進行處理,表2 給出了骨骼肌組織本構參數(shù)的均值和標準差.

表1 KNN 模型與SVR 模型統(tǒng)計指標比較Table 1 Comparison of statistical indicators between KNN model and SVR model

表2 利用KNN 模型獲得的骨骼肌組織本構參數(shù)范圍Table 2 Range of hyperelastic material parameters of skeletal muscles by KNN model

從表2 可以看出C10在0.508～0.684 kPa 之間,k1在373.361～422.881 kPa 之間,k2在31.235～36.394之間,κ 在0.261～0.299 之間.

4 結論

在本研究中提出了兩種有效的反演方法來識別骨骼肌的本構參數(shù).一種是利用K 近鄰模型并結合有限元方法來識別骨骼肌的本構參數(shù),另一種是利用支持向量機回歸模型并結合有限元方法來識別骨骼肌的本構參數(shù).本研究中的數(shù)據(jù)集由有限元方法來生成,從骨骼肌組織的單軸壓縮實驗獲得的實驗數(shù)據(jù)用來驗證本文提出的反演方法的準確性.最后,用相關系數(shù)R和決定系數(shù)R2來直觀地衡量本文提出的兩種模型的準確性,得出如下結論:

(1)本文利用K 近鄰模型和SVR 模型確定了可靠的骨骼肌組織的本構參數(shù)的范圍,所得到的本構參數(shù)的計算響應與骨骼肌組織的實驗曲線吻合良好.

(2)利用K 近鄰模型和SVR 模型在預測骨骼肌組織的本構參數(shù)都具有很大的潛力.本文提出的兩種模型和有限元相結合的反求方法也可以應用于其他具有超彈性結構的模型的軟組織.