基于仿真和智能算法骨骼肌超彈性本構(gòu)參數(shù)的反演方法研究1)

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

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

引言

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

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

1 有限元模型的建立

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

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

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

圖1 模擬實(shí)驗(yàn)設(shè)置的骨骼肌有限元模型Fig.1 Finite element model of skeletal muscles for simulating experiment

2 K 近鄰與支持向量機(jī)預(yù)測(cè)模型構(gòu)建

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

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

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

3 預(yù)測(cè)結(jié)果與討論

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

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

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

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

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

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

圖3 KNN 模型和SVR 模型測(cè)試集預(yù)測(cè)準(zhǔn)確率比較Fig.3 Comparison of test set prediction accuracy between KNN model and SVR model

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

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

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

圖5 利用SVR 模型獲得的本構(gòu)參數(shù)的計(jì)算響應(yīng)與實(shí)驗(yàn)數(shù)據(jù)以及與利用KNN 模型獲得的計(jì)算響應(yīng)對(duì)比Fig.5 Comparison of calculated response between KNN model and SVR model

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

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

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

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

表2 利用KNN 模型獲得的骨骼肌組織本構(gòu)參數(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 結(jié)論

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

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

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