在手術仿真、手術規劃、醫學診斷等眾多醫學應用中，人體軟組織力學特性建模都是一個關鍵問題。為了開發手術仿真中人體軟組織大變形生物力學模型，本文介紹了自適應準線性粘彈性（AQLV）模型，對人體前臂軟組織進行了印壓試驗，采用增量斜坡保持試驗得到模型參數，通過增量斜坡保持試驗、大應變斜坡保持試驗和大應變正弦循環加載試驗來檢驗該模型的預測能力。本文研究結果表明，AQLV 模型能夠較好地預測這三種加載工況下的試驗數據。因此，本研究的結論認為，該 AQLV 模型對在體軟組織在大應變時的非線性粘彈性的描述具有一定的實用性和可參照性，有望作為手術仿真或醫學診斷中軟組織模型的選擇之一。
引用本文： 王恒, 桑元俊. 在體軟組織非線性粘彈性的自適應準線性粘彈性模型研究. 生物醫學工程學雜志, 2017, 34(5): 702-712. doi: 10.7507/1001-5515.201610028 復制
引言
在很多醫學應用領域中，比如手術仿真（surgery simulation）、手術規劃（surgery planning）以及醫學診斷（diagnosis）等，人體軟組織力學特性建模都是一個關鍵問題[1]。在實際手術過程中，“剛性”的手術器械常常需要與“柔性”的軟組織進行接觸操作，即對軟組織進行探查、夾持、切割等操作。在虛擬手術仿真系統（virtual surgery simulation system）中，手術器械與人體軟組織之間的交互過程需要通過建立軟組織形變模型來描述，軟組織形變模型的精度直接決定了手術仿真系統中視覺反饋和觸覺反饋的精度、速度和仿真效果[2-3]；同時，軟組織形變模型精度越高，計算實時性會降低。也就是說，在虛擬手術仿真設計中，人體軟組織形變模型的精度和實時性形成最主要的一對矛盾。所以說，在虛擬手術仿真或者診斷中，很有必要開展對人體軟組織力學特性建模的研究。
前人對手術仿真中軟組織形變模型進行了大量的研究。d’Aulignac 等[4]采用彈簧振子模型來模擬人體軟組織特性，采用隱式積分方法提高仿真實時性。Okamura 等[5]建立了針穿刺入牛肝臟軟組織的力學模型，軟組織總穿刺力由接觸剛度力、摩擦力和切割力疊加得到，采用非線性彈簧模擬接觸剛度力，改進的 Karnopp 模型模擬摩擦力，并以恒定常數模擬切割力。Hu 等[6]比較了超彈性模型 Mooney-Rivlin 模型與 Ogeden 模型的差異性，認為 Ogeden 超彈性模型可以較好地模擬肝臟組織大變形情況下的力學特性。Moreira 等[7]采用線性粘彈性 Kelvin-Boltzmann 模型模擬軟組織來提高實時性。Kobayashi 等[8]采用分數階導數模型描述表面接觸形變，但存在計算量大、耗時長等問題。Pourhosseini 等[9]采用線性粘彈性 Kelvin-Voigt 模型模擬軟組織在印壓試驗下的變形。
從上述文獻來看，在大部分的虛擬手術仿真設計中，由于受到實時性的制約，采用的模型主要為線彈性或超彈性，以及線性粘彈性等。但是，當軟組織的形變較大時，軟組織實際上會表現出明顯的非線性粘彈性等，此時如果繼續采用線彈性或超彈性模型，則形變仿真結果與實際形變結果之間就會出現較大誤差，導致軟組織在大變形時仿真精度降低。這是因為軟組織本身是一種具有非線性、不可壓縮性、各向異性和粘彈性等復雜性質的材料[10]，所以為了提高虛擬手術仿真設計中軟組織大變形時的仿真精度，有必要對人體軟組織的大變形非線性粘彈性屬性進行研究。
Fung 等[11]提出的準線性粘彈性（quasi-linear viscoelastic，QLV）模型較好地描述了各種軟組織大變形時的特性。但是該模型存在兩個局限性[12]，一是該模型假設在任何瞬時應變條件下應力響應都必須與同一個減松弛函數成一定函數關系，但是實際上很多軟組織不滿足這條假設；二是該模型參數的擬合由于存在卷積部分，擬合步驟多，計算耗時長。Nekouzadeh 等[13]對 QLV 模型進行了簡化，其減松弛函數可以根據不同應變水平自適應地采取不同的參數，稱之為自適應準線性粘彈性（adaptive quasi-linear viscoelastic model，AQLV）模型，并以膠原凝膠為實驗材料進行了增量斜坡保持實驗，經過分析后認為該 AQLV 模型對增量斜坡保持實驗數據具有較好的預測能力。AQLV 模型在線性粘彈性函數里加入了非線性屬性，所以屬于一種非線性粘彈性模型，該模型假設在不同的應變水平時，松弛函數的各參數與應變水平成一定的函數關系，從而使得應力變化可以根據不同的應變水平自適應地調用松弛函數的相應參數。Quaia 等[14]利用 AQLV 模型對成年恒河猴眼肌組織麻醉狀態時的增量斜坡保持印壓特性進行了研究，但是他們并沒有對大應變斜坡保持加載和循環加載狀況下的 AQLV 模型的預測能力進行研究。
基于現有的文獻研究基礎，考慮到在體試驗的操作方便性，本文對正常人體前臂軟組織進行了大變形印壓試驗。為了避免 Fung 等[11]提出的 QLV 模型的局限性，本文采用增量斜坡保持試驗分析不同應變條件下的組織力學特性，并采用 AQLV 模型進行參數擬合辨識。為了進一步驗證 AQLV 模型的預測能力，本文在增量斜坡保持試驗、大應變斜坡保持試驗以及大應變正弦循環加載試驗共 3 種典型工況下，對正常人前臂軟組織的生物力學特性進行預測分析。具體步驟為：首先，對 AQLV 模型的理論進行了詳細分析，提出了不同工況加載下的應力預測函數；其次，對該模型的增量斜坡保持加載進行了試驗以及擬合分析，并進一步分析了模型各參數與應變水平的函數關系；最后，對增量斜坡保持加載、大應變斜坡保持加載和大應變正弦循環加載進行了模型的預測能力分析。本文研究成果有望作為手術仿真或醫學診斷中軟組織模型的選擇之一。
1 材料與方法
1.1 AQLV 模型的理論分析
AQLV 模型假設應力
與時間 t、應變
存在非線性函數關系。在一定應變
下，軟組織的應力由長期的穩態應力
以及多個不同時間松弛常數的松弛應力
累加而成，即
 |
在 AQLV 模型中，對第 i 個松弛應力
而言，其與相應的應變通過在一個線性卷積積分的中間變量
（稱之為粘彈性應變）形成函數關系，即
 |
其中，
為應變的純非線性函數，
表示為應力基于應變時間歷程的關系。
可以假設由并聯的 Maxwell 單元（注：Maxwell 單元是由一個理想彈簧和一個理想黏壺串聯組成的一種線性粘彈性模型）組成，設定時間常數與應變無關。即
 |
根據式（1）、式（2）、式（3），AQLV 模型可以寫成一般形式，即
 |
假設在 t ≤ 0 時，
為 0，同時，假設 N=3，則加載時間區間為
時，式（4）可以寫為
 |
當采用階躍加載或斜坡保持加載時，若 t ≥ T，則
將保持不變，即維持為
，式（5）轉換為
 |
在式（5）或式（6）中，
、k1、k2、k3、
、
、
共 7 個材料參數需通過辨識得到。
、k1、k2、k3 等 4 個參數假設分別與應變水平存在某種函數關系，需要試驗擬合。
、
、
等 3 個時間常數假設與應變水平變化無關，僅通過某一個應變水平下的試驗數據進行擬合。
1.2 在不同加載工況下 AQLV 模型的理論推導
軟組織的粘彈性屬性包括應力松弛、應變蠕變、回滯等，AQLV 模型可根據軟組織的應力松弛試驗進行分析。理論上，應力松弛試驗應該通過應變階躍加載到某應變水平，然后應變值保持一定時間不變，測量應力值隨時間的衰減情況。但是，實際操作時應變加載過程需要一定的時間，過快的加載速度也將導致測量數據的震蕩，并有可能導致組織的破壞。所以，應變的階躍加載方法實際上并不可行。通常的試驗方法是通過應變在一定速度下進行加載，然后應變值保持一定時間不變，測量組織應力松弛過程的應力數據。該試驗方法是實際操作中軟組織應力松弛試驗的一種方法，稱之為斜坡保持加載方法。在斜坡保持加載方法中，當應變從初始值加載至設定值時，該過程稱之為上升階段；當應變達到設定值并保持不變時，應力隨時間發生松弛至穩定狀態，該過程稱之為保持階段。
軟組織在大變形情形下表現出非線性粘彈性特性，已有理論模型（如 QLV 模型），并假設模型各參數值不隨應變水平變化而變化。為了考察 AQLV 模型各參數值在不同應變水平下的變化情況，可通過等增量斜坡保持試驗方法進行試驗，具體為：設定最終大變形時應變水平數值為
，然后平均分為 N 段，以應變 0 時為初始值開始斜坡加載保持試驗，各段應變加載值為
/N，各段加載后應變保持不變，直到應力松弛至穩定狀態，然后再繼續下一段斜坡保持加載。該試驗方法稱之為增量斜坡保持加載方法。該方法中各應變增量并不是一定要相等，但建議為等增量加載，方便后續參數擬合。本文采用增量斜坡保持試驗來辨識得到各個材料參數，然后再擬合各個材料參數與應變的函數關系。
當采用增量斜坡保持加載方法試驗時，第 n 段斜坡保持試驗的應變為
 |
式（7）中，
為第 n 段松弛試驗的應變歷程，
為第 n 段松弛試驗的應變幅度，T 為第 n 段松弛試驗的斜坡加載時間。則根據式（5）、式（6）可以得到第 n 段斜坡保持試驗的應力為
 |
式（8）中，
、
分別為第 n 段斜坡保持試驗上升階段和保持階段的應力預測。采用
來擬合第 n 段保持階段試驗數據以得到各個材料參數，擬合得到的材料參數代入
中來檢驗第 n 段上升階段模型的預測能力。
為了檢驗 AQLV 模型的大應變粘彈性的預測能力，通過典型加載工況，如大應變斜坡保持加載、大應變正弦循環加載來檢驗模型的預測能力。大應變斜坡保持加載方法為應變在一定速度下直接加載至設定數值，然后應變保持一定時間不變直到應力松弛至穩定狀態。因為該應變數值較大，所以稱之為大應變斜坡保持加載方法。在大應變情形下，還可以設定應變為正弦循環加載方式，該正弦信號的峰峰值即為應變設定最大值，該加載方式稱之為大應變正弦循環加載方法。
當加載方式為大應變斜坡保持加載時，
，代入公式（5）得到上升階段的應力預測公式，代入式（6）得到保持階段預測公式，即
 |
式（9）中，v 為加載速度，H 為軟組織厚度，
為應變率，
、
分別為大應變斜坡保持試驗上升階段和保持階段的應力預測。
當加載方式為大應變正弦循環加載時，應變加載公式為
，其中 A 為加載應變幅值，f 為加載頻率。把該應變函數公式代入式（5），則得到具體應力預測公式。
1.3 試驗方法
對一名正常男性志愿者（年齡：35 歲；身高：172 cm；體重：79 kg）的左前臂進行測試，志愿者前臂皮膚等軟組織良好。試驗前明確告知其試驗目的、流程等，并得到志愿者知情同意書。為降低化妝品對軟組織力學性能的影響，試驗前 3 天不可涂防曬霜等任何化妝品，試驗前清洗手臂且自然晾干。試驗位置為從肘窩中心向前臂遠端軸向移動 4 cm，再轉向橈骨側移動 1 cm 處，該區域軟組織的底層橈骨面較為平緩。
試驗裝置如圖 1 所示。采用長度為 40 mm，直徑為 4 mm 的光滑平頭不銹鋼圓柱為印壓頭。采用萬能試驗機（3330，Bose co. inc，美國）實現位移驅動控制。采用量程為 3 N，誤差精度為 0.02% 的力傳感器（FGA-3 N，Forsentek co. inc，中國）測量軟組織的力學反應。采用數據采集平臺（PCD-300A，KEYENCE co.inc，日本）進行數據采集與上位機處理，數據采樣頻率為 100 Hz。前臂印壓位置的軟組織厚度由電子計算機斷層掃描（computed tomography，CT）采集得到。
圖1
試驗裝置圖
Figure1.
Figure of the indentation equipment
前臂軟組織非線性粘彈性特性的模型參數辨識試驗為增量斜坡保持試驗。模型的預測能力檢驗的試驗包括三種加載工況，即增量斜坡保持試驗、大應變斜坡保持試驗、大應變正弦循環試驗。考慮到正弦循環加載工況時應變施加值較為簡單，故不單獨列圖說明。而增量斜坡保持試驗與大應變斜坡保持試驗時應變施加值分別如圖 2 所示。在增量斜坡保持試驗中：設定試驗加載速率為 10 mm/s，每次加載位移為 2 mm，則每階段的上升時間為 0.2 s；印壓區域軟組織厚度根據 CT 掃描得到，為 30.72 mm，則加載應變分別為 0.065 1、0.130 2、0.195 3、0.260 4，保持時間分別為 30 s、30 s、150 s、150 s。由于整個試驗過程中受試者需要盡量保持靜止以實現數據測量的準確性，而長時間保持靜止容易導致受試者不舒適，并且受試者輕微動作容易導致數據測量誤差，所以根據初期試驗結果分析后發現，在較小應變（即應變分別為 0.065 1、0.130 2）時，由于軟組織的應力松弛可以較快達到平衡狀態（約 10 s 內即可達到平衡狀態），由此設定小應變保持時間為 30 s，這樣可使得整個試驗時間縮短，從而一定程度地減免受試者輕微動作改變帶來的數據測量誤差。在較大應變（即應變分別為 0.195 3、0.260 4）時，根據初期試驗結果分析后發現，軟組織的應力松弛可以在約 50 s 內達到平衡狀態，由此設定大應變保持時間為 150 s，這樣可以觀測到應力松弛至穩定狀態的應力數據。根據式（8）可知，將增量斜坡保持試驗的 4 個保持階段試驗數據用來做擬合分析，4 個上升階段試驗數據可分別用來檢驗模型的預測能力。
圖2
前臂軟組織非線性粘彈性特性測試方法
Figure2.
Test protocols for the nonlinear viscoelastic properties of forearm soft tissues
為了進一步檢測模型的大應變加載預測能力，如圖 2 大應變斜坡保持試驗所示，加載速率為 10 mm/s，應變為 0.325 5（即位移幅值為 10 mm），檢驗模型在應變高速加載下的預測能力；同時，為了檢驗模型在應變低速加載下的預測能力，改變加載速率為 0.5 mm/s，應變設定值仍為 0.325 5，對試驗數據和預測數據進行比較分析。另外，還采用位移幅值為 6 mm，頻率為 0.2 Hz 的正弦循環加載試驗來檢驗模型對循環加載的預測能力。
所有試驗均進行 10 次三角波預調制，以實現軟組織達到穩定狀態。所有試驗均測量 3 次，取其平均值，從而提高數據測量的準確性。
2 結果
2.1 AQLV 模型的參數擬合
根據如圖 2 所示測試方法，可以分別得到如圖 3 所示的增量斜坡保持試驗結果和大應變斜坡保持試驗結果。結果表明，經過 7 次三角波預調制后軟組織測量結果趨于穩定，因此試驗時采用 10 次預調制是可行的。
如圖 3 所示的增量斜坡保持試驗結果，可以用來辨識 AQLV 模型的各個材料參數，辨識方法為根據式（8）中
來擬合 4 段保持階段試驗數據。首先對應變為 0.260 4 時的保持階段數據進行擬合來優化選擇時間常數，擬合方法為非線性最小二乘法，擬合算法為信賴域算法，得到時間常數
、
、
分別為 0.4、4、40，此時擬合精度最高，擬合優度為 R2 = 0.946 1。因為本文假設時間常數與應變水平變化無關，所以根據確定后的時間常數，再分別對其他各段應變水平時保持階段松弛數據進行擬合，從而分別得到各段應變水平時模型的材料參數。各應變水平保持階段試驗數據的擬合結果如圖 4 所示。辨識得到的 AQLV 模型中各個材料參數數值，如表 1 所示，
為與應變水平存在函數關系的穩態應力，k1、k2、k3 分別為應變水平的純非線性函數，
、
、
分別為 3 個與應變水平變化不相關的時間常數。根據表 1 得到的材料參數的離散值，可以進一步擬合各個參數與應變水平的函數關系，擬合方法均采用三次多項式擬合，即擬合公式為 f（x）= p1 · x3 + p2 · x2 + p3 · x + p4，需要辨識得到多項式中各個系數的數值，擬合結果分別如圖 5 所示，多項式中各個系數數值及擬合優度如表 2 所示，擬合優度均大于 0.98，說明采用三次多項式擬合各參數與應變水平的函數關系具有一定的合理性。
圖3
加載試驗結果
Figure3.
Experimental results
圖4
增量斜坡保持試驗模型參數擬合結果
Figure4.
Incremental ramp-and-hold test fitting results
表1
擬合得到的 AQLV 模型的材料參數（τ = 0.4，τ = 4，τ = 40）
Table1.
Fitted mechanical parameters of the AQLV model (τ = 0.4, τ = 4, τ = 40)
圖5
AQLV 模型各個材料參數與應變水平的函數關系
Figure5.
The functional relationship between each material parameter of the AQLV model and the strain level
表2
AQLV 模型各個材料參數與應變水平的函數關系
Table2.
The functional relationship between each material parame-ter of the AQLV model and the strain level
2.2 AQLV 模型的預測能力
根據增量斜坡保持試驗的各段保持階段數據可以擬合得到 AQLV 模型的各個材料參數。在此基礎上，采用式（8）來預測增量斜坡保持試驗的試驗數據。如圖 6 所示為增量斜坡保持試驗的預測結果和試驗結果的對比，為了更清晰地觀測數據的對比效果，對試驗數據與預測數據的對比圖部分區域進行局部放大，得到對比數據的局部放大圖，可以看到，得到的 AQLV 模型可以較好地預測增量斜坡保持試驗結果。進一步采用大應變斜坡保持試驗、大應變正弦循環加載試驗來檢驗 AQLV 模型的預測能力。
圖6
增量斜坡保持試驗的模型預測結果
Figure6.
Predictions of the AQLV model for incremental ramp-and-hold test
當加載方式為大應變斜坡保持加載時，如圖 2 所示，應變設置為 0.325 5，加載速度為 10 mm/s，檢驗模型在應變高速加載下的預測能力。根據式（9）可以得到預測結果如圖 7 所示，可以看到在加載上升階段，預測結果要略高于試驗結果，上升到 1 s 時應力預測值達到最大值，此時誤差為 13.5%；而在加載保持階段，預測結果在約 40 s 就達到平衡值，并和試驗結果保持一致。
圖7
大應變斜坡保持試驗預測結果
Figure7.
Predictions of the AQLV model for the single large amplitude ramp-and-hold test
為了進一步檢驗 AQLV 模型對大應變斜坡保持加載的預測性能，仍設置應變為 0.325 5，加載速度設置為 0.5 mm/s，檢驗模型在應變低速加載下的預測能力。根據式（9）可以得到預測結果如圖 7 所示，在加載上升階段，上升到 20 s 時應力預測值達到最大值，此時誤差為 7.8%；而在保持階段，預測結果同樣在約 40 s 就達到平衡值，并和試驗結果保持一致。比較圖 7 的對比結果可知，實際軟組織的力學響應受到應變加載速率的影響，而 AQLV 模型對這兩種不同加載速率的大應變加載情形都能較好地預測。
雖然 AQLV 模型對大應變斜坡保持加載的上升加載階段存在一定的跟隨誤差，但誤差范圍較小；此外，該模型可以快速達到穩態且與試驗結果保持一致。所以，可以認為 AQLV 模型能夠較好地預測大應變斜坡保持試驗結果。
當加載方式為大應變正弦循環加載時，加載頻率為 0.2 Hz，加載應變幅值為 0.097 7，則應變變化范圍為 0~0.195 3。根據式（5）可以得到預測結果。大應變正弦循環加載的試驗與預測結果如圖 8 所示。因為 AQLV 模型預測結果需要經過一段時間才能達到平衡穩定值，而正弦循環加載試驗結果為經過預調制后測量得到的穩態值。因此如圖 8 所示，為 AQLV 模型穩定狀態的預測數據和正弦循環試驗數據對比，可以認為 AQLV 模型能夠較好地預測大應變正弦循環加載試驗結果。
圖8
大應變正弦循環加載預測結果
Figure8.
Predictions of the AQLV model for the sinusoidal cyclic loading under large strain
3 討論與分析
手術器械與軟組織交互作用的模擬可以通過計算機模擬仿真來實現。其中對軟組織的大變形模擬分析時，現有較高精度的模型由于存在卷積積分步驟容易導致計算響應時間長，操控響應存在滯后遲鈍等缺點；而較低精度的模型容易導致仿真誤差變大。特別地，現有仿真模型對軟組織大變形粘彈性的生物力學特性考慮不足，為了提高大變形粘彈性模型的精度，同時簡化計算過程，本文開展了對人體軟組織大變形的非線性粘彈性特性的研究探討。軟組織大變形時盡管可以采用超彈性模型來模擬加載過程的非線性變化，但是忽略了軟組織的粘彈性屬性。在經典線性粘彈性理論中，松弛函數的參數是獨立于應變的，僅為時間的函數。在 AQLV 模型中，假設模型松弛函數的參數可以根據應變水平自適應地選擇不同的參數，即假設松弛函數的參數不僅是時間的非線性函數，還是應變的非線性函數，則 AQLV 模型應該屬于一種非線性粘彈性模型，符合軟組織自身的生物力學特性。
QLV 模型假設松弛函數的參數與應變水平無關，然而 AQLV 模型假設松弛函數參數與應變水平存在函數關系。從圖 4 和表 1 試驗結果可以看到，松弛函數的參數實際上不是完全獨立于應變水平的，所以 AQLV 模型假設松弛函數模型參數不獨立于應變水平是合理的。因為 AQLV 模型假設松弛特性由并聯 Maxwell 單元來模擬，各 Maxwell 單元的時間常數設定為不變，這樣設定是為了擬合計算方便，所以假設時間常數與應變水平無關，即
、
、
分別為常數。而在實際軟組織內部，其實是不存在 Maxwell 單元這種理想化的物質。時間常數的選擇并無嚴格的標準，一般設定為相互之間存在一定的比例，比如采用對數間隔的方法[14]，然后采用設定的優化算法來比較選擇。
如圖 5 和表 2 所示的模型各個參數與應變的擬合結果來看，各個參數與應變水平采用三次多項式擬合可以得到較高精度。當應變為 0 時，假設軟組織的狀態處于零應力狀態，所以各參數均為 0。這樣假設具有一定的計算合理性，在實際試驗時候可以通過多次預調制實現軟組織的狀態穩定。著名生物力學學家馮元楨先生也提出：“要認識生物組織的改造，必須在零應力狀態來比較”[15]。
AQLV 模型的預測能力根據增量斜坡保持試驗、大應變斜坡保持試驗及大應變正弦循環加載試驗來檢驗。從圖 6、圖 7、圖 8 的預測結果來看，AQLV 模型可以較好地預測這三種加載工況的試驗結果。其中需要說明是在圖 7 中，在相同加載應變水平但不同應變加載速率下，軟組織實際穩態應力與 AQLV 模型預測的穩態應力均較好地吻合，并且應變加載速率對軟組織穩態應力值的影響很小。但是，應變加載速率對瞬態特性的影響較大，較大的應變加載速率導致軟組織實際瞬態應力變大，且 AQLV 模型對大應變斜坡保持加載的瞬態特性的擬合存在一定的誤差。這是因為 AQLV 模型不滿足疊加原理（superposition principle），Quaia 等[14, 16]的研究也證明了該模型不滿足疊加原理。
本文研究存在一定程度的局限和不足之處。AQLV 模型在不同應變速率下的預測能力還有待于后續試驗，但是計算公式可以直接代入本文所得到的公式中來擬合和預測。本文設定 AQLV 模型由 3 項松弛函數疊加組成，在手術仿真應用時還可以根據實時性的要求減少為 1 項松弛函數，即僅需擬合
和
。另外，計算機仿真的精度不僅受到軟組織非線性粘彈性模型精度影響，仿真精度和實時性還受到網格密度的影響[9]，這些工作有待于后續研究完善。
本文對手術仿真或醫學診斷中的人體軟組織特性建模進行了探討，著重對在體軟組織大變形非線性粘彈性特性進行了理論和試驗分析。以人體前臂軟組織生物力學特性為例進行印壓試驗，認為 AQLV 模型對前臂軟組織增量斜坡保持、大應變斜坡保持以及大應變正弦循環等三種加載工況下的試驗結果都能較好地預測。后續研究工作將著重對 AQLV 模型在手術仿真中的軟件設計進行研究。
引言
在很多醫學應用領域中，比如手術仿真（surgery simulation）、手術規劃（surgery planning）以及醫學診斷（diagnosis）等，人體軟組織力學特性建模都是一個關鍵問題[1]。在實際手術過程中，“剛性”的手術器械常常需要與“柔性”的軟組織進行接觸操作，即對軟組織進行探查、夾持、切割等操作。在虛擬手術仿真系統（virtual surgery simulation system）中，手術器械與人體軟組織之間的交互過程需要通過建立軟組織形變模型來描述，軟組織形變模型的精度直接決定了手術仿真系統中視覺反饋和觸覺反饋的精度、速度和仿真效果[2-3]；同時，軟組織形變模型精度越高，計算實時性會降低。也就是說，在虛擬手術仿真設計中，人體軟組織形變模型的精度和實時性形成最主要的一對矛盾。所以說，在虛擬手術仿真或者診斷中，很有必要開展對人體軟組織力學特性建模的研究。
前人對手術仿真中軟組織形變模型進行了大量的研究。d’Aulignac 等[4]采用彈簧振子模型來模擬人體軟組織特性，采用隱式積分方法提高仿真實時性。Okamura 等[5]建立了針穿刺入牛肝臟軟組織的力學模型，軟組織總穿刺力由接觸剛度力、摩擦力和切割力疊加得到，采用非線性彈簧模擬接觸剛度力，改進的 Karnopp 模型模擬摩擦力，并以恒定常數模擬切割力。Hu 等[6]比較了超彈性模型 Mooney-Rivlin 模型與 Ogeden 模型的差異性，認為 Ogeden 超彈性模型可以較好地模擬肝臟組織大變形情況下的力學特性。Moreira 等[7]采用線性粘彈性 Kelvin-Boltzmann 模型模擬軟組織來提高實時性。Kobayashi 等[8]采用分數階導數模型描述表面接觸形變，但存在計算量大、耗時長等問題。Pourhosseini 等[9]采用線性粘彈性 Kelvin-Voigt 模型模擬軟組織在印壓試驗下的變形。
從上述文獻來看，在大部分的虛擬手術仿真設計中，由于受到實時性的制約，采用的模型主要為線彈性或超彈性，以及線性粘彈性等。但是，當軟組織的形變較大時，軟組織實際上會表現出明顯的非線性粘彈性等，此時如果繼續采用線彈性或超彈性模型，則形變仿真結果與實際形變結果之間就會出現較大誤差，導致軟組織在大變形時仿真精度降低。這是因為軟組織本身是一種具有非線性、不可壓縮性、各向異性和粘彈性等復雜性質的材料[10]，所以為了提高虛擬手術仿真設計中軟組織大變形時的仿真精度，有必要對人體軟組織的大變形非線性粘彈性屬性進行研究。
Fung 等[11]提出的準線性粘彈性（quasi-linear viscoelastic，QLV）模型較好地描述了各種軟組織大變形時的特性。但是該模型存在兩個局限性[12]，一是該模型假設在任何瞬時應變條件下應力響應都必須與同一個減松弛函數成一定函數關系，但是實際上很多軟組織不滿足這條假設；二是該模型參數的擬合由于存在卷積部分，擬合步驟多，計算耗時長。Nekouzadeh 等[13]對 QLV 模型進行了簡化，其減松弛函數可以根據不同應變水平自適應地采取不同的參數，稱之為自適應準線性粘彈性（adaptive quasi-linear viscoelastic model，AQLV）模型，并以膠原凝膠為實驗材料進行了增量斜坡保持實驗，經過分析后認為該 AQLV 模型對增量斜坡保持實驗數據具有較好的預測能力。AQLV 模型在線性粘彈性函數里加入了非線性屬性，所以屬于一種非線性粘彈性模型，該模型假設在不同的應變水平時，松弛函數的各參數與應變水平成一定的函數關系，從而使得應力變化可以根據不同的應變水平自適應地調用松弛函數的相應參數。Quaia 等[14]利用 AQLV 模型對成年恒河猴眼肌組織麻醉狀態時的增量斜坡保持印壓特性進行了研究，但是他們并沒有對大應變斜坡保持加載和循環加載狀況下的 AQLV 模型的預測能力進行研究。
基于現有的文獻研究基礎，考慮到在體試驗的操作方便性，本文對正常人體前臂軟組織進行了大變形印壓試驗。為了避免 Fung 等[11]提出的 QLV 模型的局限性，本文采用增量斜坡保持試驗分析不同應變條件下的組織力學特性，并采用 AQLV 模型進行參數擬合辨識。為了進一步驗證 AQLV 模型的預測能力，本文在增量斜坡保持試驗、大應變斜坡保持試驗以及大應變正弦循環加載試驗共 3 種典型工況下，對正常人前臂軟組織的生物力學特性進行預測分析。具體步驟為：首先，對 AQLV 模型的理論進行了詳細分析，提出了不同工況加載下的應力預測函數；其次，對該模型的增量斜坡保持加載進行了試驗以及擬合分析，并進一步分析了模型各參數與應變水平的函數關系；最后，對增量斜坡保持加載、大應變斜坡保持加載和大應變正弦循環加載進行了模型的預測能力分析。本文研究成果有望作為手術仿真或醫學診斷中軟組織模型的選擇之一。
1 材料與方法
1.1 AQLV 模型的理論分析
AQLV 模型假設應力
與時間 t、應變
存在非線性函數關系。在一定應變
下，軟組織的應力由長期的穩態應力
以及多個不同時間松弛常數的松弛應力
累加而成，即
|
在 AQLV 模型中，對第 i 個松弛應力
而言，其與相應的應變通過在一個線性卷積積分的中間變量
（稱之為粘彈性應變）形成函數關系，即
|
其中，
為應變的純非線性函數，
表示為應力基于應變時間歷程的關系。
可以假設由并聯的 Maxwell 單元（注：Maxwell 單元是由一個理想彈簧和一個理想黏壺串聯組成的一種線性粘彈性模型）組成，設定時間常數與應變無關。即
|
根據式（1）、式（2）、式（3），AQLV 模型可以寫成一般形式，即
|
假設在 t ≤ 0 時，
為 0，同時，假設 N=3，則加載時間區間為
時，式（4）可以寫為
|
當采用階躍加載或斜坡保持加載時，若 t ≥ T，則
將保持不變，即維持為
，式（5）轉換為
|
在式（5）或式（6）中，
、k1、k2、k3、
、
、
共 7 個材料參數需通過辨識得到。
、k1、k2、k3 等 4 個參數假設分別與應變水平存在某種函數關系，需要試驗擬合。
、
、
等 3 個時間常數假設與應變水平變化無關，僅通過某一個應變水平下的試驗數據進行擬合。
1.2 在不同加載工況下 AQLV 模型的理論推導
軟組織的粘彈性屬性包括應力松弛、應變蠕變、回滯等，AQLV 模型可根據軟組織的應力松弛試驗進行分析。理論上，應力松弛試驗應該通過應變階躍加載到某應變水平，然后應變值保持一定時間不變，測量應力值隨時間的衰減情況。但是，實際操作時應變加載過程需要一定的時間，過快的加載速度也將導致測量數據的震蕩，并有可能導致組織的破壞。所以，應變的階躍加載方法實際上并不可行。通常的試驗方法是通過應變在一定速度下進行加載，然后應變值保持一定時間不變，測量組織應力松弛過程的應力數據。該試驗方法是實際操作中軟組織應力松弛試驗的一種方法，稱之為斜坡保持加載方法。在斜坡保持加載方法中，當應變從初始值加載至設定值時，該過程稱之為上升階段；當應變達到設定值并保持不變時，應力隨時間發生松弛至穩定狀態，該過程稱之為保持階段。
軟組織在大變形情形下表現出非線性粘彈性特性，已有理論模型（如 QLV 模型），并假設模型各參數值不隨應變水平變化而變化。為了考察 AQLV 模型各參數值在不同應變水平下的變化情況，可通過等增量斜坡保持試驗方法進行試驗，具體為：設定最終大變形時應變水平數值為
，然后平均分為 N 段，以應變 0 時為初始值開始斜坡加載保持試驗，各段應變加載值為
/N，各段加載后應變保持不變，直到應力松弛至穩定狀態，然后再繼續下一段斜坡保持加載。該試驗方法稱之為增量斜坡保持加載方法。該方法中各應變增量并不是一定要相等，但建議為等增量加載，方便后續參數擬合。本文采用增量斜坡保持試驗來辨識得到各個材料參數，然后再擬合各個材料參數與應變的函數關系。
當采用增量斜坡保持加載方法試驗時，第 n 段斜坡保持試驗的應變為
|
式（7）中，
為第 n 段松弛試驗的應變歷程，
為第 n 段松弛試驗的應變幅度，T 為第 n 段松弛試驗的斜坡加載時間。則根據式（5）、式（6）可以得到第 n 段斜坡保持試驗的應力為
|
式（8）中，
、
分別為第 n 段斜坡保持試驗上升階段和保持階段的應力預測。采用
來擬合第 n 段保持階段試驗數據以得到各個材料參數，擬合得到的材料參數代入
中來檢驗第 n 段上升階段模型的預測能力。
為了檢驗 AQLV 模型的大應變粘彈性的預測能力，通過典型加載工況，如大應變斜坡保持加載、大應變正弦循環加載來檢驗模型的預測能力。大應變斜坡保持加載方法為應變在一定速度下直接加載至設定數值，然后應變保持一定時間不變直到應力松弛至穩定狀態。因為該應變數值較大，所以稱之為大應變斜坡保持加載方法。在大應變情形下，還可以設定應變為正弦循環加載方式，該正弦信號的峰峰值即為應變設定最大值，該加載方式稱之為大應變正弦循環加載方法。
當加載方式為大應變斜坡保持加載時，
，代入公式（5）得到上升階段的應力預測公式，代入式（6）得到保持階段預測公式，即
|
式（9）中，v 為加載速度，H 為軟組織厚度，
為應變率，
、
分別為大應變斜坡保持試驗上升階段和保持階段的應力預測。
當加載方式為大應變正弦循環加載時，應變加載公式為
，其中 A 為加載應變幅值，f 為加載頻率。把該應變函數公式代入式（5），則得到具體應力預測公式。
1.3 試驗方法
對一名正常男性志愿者（年齡：35 歲；身高：172 cm；體重：79 kg）的左前臂進行測試，志愿者前臂皮膚等軟組織良好。試驗前明確告知其試驗目的、流程等，并得到志愿者知情同意書。為降低化妝品對軟組織力學性能的影響，試驗前 3 天不可涂防曬霜等任何化妝品，試驗前清洗手臂且自然晾干。試驗位置為從肘窩中心向前臂遠端軸向移動 4 cm，再轉向橈骨側移動 1 cm 處，該區域軟組織的底層橈骨面較為平緩。
試驗裝置如圖 1 所示。采用長度為 40 mm，直徑為 4 mm 的光滑平頭不銹鋼圓柱為印壓頭。采用萬能試驗機（3330，Bose co. inc，美國）實現位移驅動控制。采用量程為 3 N，誤差精度為 0.02% 的力傳感器（FGA-3 N，Forsentek co. inc，中國）測量軟組織的力學反應。采用數據采集平臺（PCD-300A，KEYENCE co.inc，日本）進行數據采集與上位機處理，數據采樣頻率為 100 Hz。前臂印壓位置的軟組織厚度由電子計算機斷層掃描（computed tomography，CT）采集得到。
圖1
試驗裝置圖
Figure1.
Figure of the indentation equipment
前臂軟組織非線性粘彈性特性的模型參數辨識試驗為增量斜坡保持試驗。模型的預測能力檢驗的試驗包括三種加載工況，即增量斜坡保持試驗、大應變斜坡保持試驗、大應變正弦循環試驗。考慮到正弦循環加載工況時應變施加值較為簡單，故不單獨列圖說明。而增量斜坡保持試驗與大應變斜坡保持試驗時應變施加值分別如圖 2 所示。在增量斜坡保持試驗中：設定試驗加載速率為 10 mm/s，每次加載位移為 2 mm，則每階段的上升時間為 0.2 s；印壓區域軟組織厚度根據 CT 掃描得到，為 30.72 mm，則加載應變分別為 0.065 1、0.130 2、0.195 3、0.260 4，保持時間分別為 30 s、30 s、150 s、150 s。由于整個試驗過程中受試者需要盡量保持靜止以實現數據測量的準確性，而長時間保持靜止容易導致受試者不舒適，并且受試者輕微動作容易導致數據測量誤差，所以根據初期試驗結果分析后發現，在較小應變（即應變分別為 0.065 1、0.130 2）時，由于軟組織的應力松弛可以較快達到平衡狀態（約 10 s 內即可達到平衡狀態），由此設定小應變保持時間為 30 s，這樣可使得整個試驗時間縮短，從而一定程度地減免受試者輕微動作改變帶來的數據測量誤差。在較大應變（即應變分別為 0.195 3、0.260 4）時，根據初期試驗結果分析后發現，軟組織的應力松弛可以在約 50 s 內達到平衡狀態，由此設定大應變保持時間為 150 s，這樣可以觀測到應力松弛至穩定狀態的應力數據。根據式（8）可知，將增量斜坡保持試驗的 4 個保持階段試驗數據用來做擬合分析，4 個上升階段試驗數據可分別用來檢驗模型的預測能力。
圖2
前臂軟組織非線性粘彈性特性測試方法
Figure2.
Test protocols for the nonlinear viscoelastic properties of forearm soft tissues
為了進一步檢測模型的大應變加載預測能力，如圖 2 大應變斜坡保持試驗所示，加載速率為 10 mm/s，應變為 0.325 5（即位移幅值為 10 mm），檢驗模型在應變高速加載下的預測能力；同時，為了檢驗模型在應變低速加載下的預測能力，改變加載速率為 0.5 mm/s，應變設定值仍為 0.325 5，對試驗數據和預測數據進行比較分析。另外，還采用位移幅值為 6 mm，頻率為 0.2 Hz 的正弦循環加載試驗來檢驗模型對循環加載的預測能力。
所有試驗均進行 10 次三角波預調制，以實現軟組織達到穩定狀態。所有試驗均測量 3 次，取其平均值，從而提高數據測量的準確性。
2 結果
2.1 AQLV 模型的參數擬合
根據如圖 2 所示測試方法，可以分別得到如圖 3 所示的增量斜坡保持試驗結果和大應變斜坡保持試驗結果。結果表明，經過 7 次三角波預調制后軟組織測量結果趨于穩定，因此試驗時采用 10 次預調制是可行的。
如圖 3 所示的增量斜坡保持試驗結果，可以用來辨識 AQLV 模型的各個材料參數，辨識方法為根據式（8）中
來擬合 4 段保持階段試驗數據。首先對應變為 0.260 4 時的保持階段數據進行擬合來優化選擇時間常數，擬合方法為非線性最小二乘法，擬合算法為信賴域算法，得到時間常數
、
、
分別為 0.4、4、40，此時擬合精度最高，擬合優度為 R2 = 0.946 1。因為本文假設時間常數與應變水平變化無關，所以根據確定后的時間常數，再分別對其他各段應變水平時保持階段松弛數據進行擬合，從而分別得到各段應變水平時模型的材料參數。各應變水平保持階段試驗數據的擬合結果如圖 4 所示。辨識得到的 AQLV 模型中各個材料參數數值，如表 1 所示，
為與應變水平存在函數關系的穩態應力，k1、k2、k3 分別為應變水平的純非線性函數，
、
、
分別為 3 個與應變水平變化不相關的時間常數。根據表 1 得到的材料參數的離散值，可以進一步擬合各個參數與應變水平的函數關系，擬合方法均采用三次多項式擬合，即擬合公式為 f（x）= p1 · x3 + p2 · x2 + p3 · x + p4，需要辨識得到多項式中各個系數的數值，擬合結果分別如圖 5 所示，多項式中各個系數數值及擬合優度如表 2 所示，擬合優度均大于 0.98，說明采用三次多項式擬合各參數與應變水平的函數關系具有一定的合理性。
圖3
加載試驗結果
Figure3.
Experimental results
圖4
增量斜坡保持試驗模型參數擬合結果
Figure4.
Incremental ramp-and-hold test fitting results
表1
擬合得到的 AQLV 模型的材料參數（τ = 0.4，τ = 4，τ = 40）
Table1.
Fitted mechanical parameters of the AQLV model (τ = 0.4, τ = 4, τ = 40)
圖5
AQLV 模型各個材料參數與應變水平的函數關系
Figure5.
The functional relationship between each material parameter of the AQLV model and the strain level
表2
AQLV 模型各個材料參數與應變水平的函數關系
Table2.
The functional relationship between each material parame-ter of the AQLV model and the strain level
2.2 AQLV 模型的預測能力
根據增量斜坡保持試驗的各段保持階段數據可以擬合得到 AQLV 模型的各個材料參數。在此基礎上，采用式（8）來預測增量斜坡保持試驗的試驗數據。如圖 6 所示為增量斜坡保持試驗的預測結果和試驗結果的對比，為了更清晰地觀測數據的對比效果，對試驗數據與預測數據的對比圖部分區域進行局部放大，得到對比數據的局部放大圖，可以看到，得到的 AQLV 模型可以較好地預測增量斜坡保持試驗結果。進一步采用大應變斜坡保持試驗、大應變正弦循環加載試驗來檢驗 AQLV 模型的預測能力。
圖6
增量斜坡保持試驗的模型預測結果
Figure6.
Predictions of the AQLV model for incremental ramp-and-hold test
當加載方式為大應變斜坡保持加載時，如圖 2 所示，應變設置為 0.325 5，加載速度為 10 mm/s，檢驗模型在應變高速加載下的預測能力。根據式（9）可以得到預測結果如圖 7 所示，可以看到在加載上升階段，預測結果要略高于試驗結果，上升到 1 s 時應力預測值達到最大值，此時誤差為 13.5%；而在加載保持階段，預測結果在約 40 s 就達到平衡值，并和試驗結果保持一致。
圖7
大應變斜坡保持試驗預測結果
Figure7.
Predictions of the AQLV model for the single large amplitude ramp-and-hold test
為了進一步檢驗 AQLV 模型對大應變斜坡保持加載的預測性能，仍設置應變為 0.325 5，加載速度設置為 0.5 mm/s，檢驗模型在應變低速加載下的預測能力。根據式（9）可以得到預測結果如圖 7 所示，在加載上升階段，上升到 20 s 時應力預測值達到最大值，此時誤差為 7.8%；而在保持階段，預測結果同樣在約 40 s 就達到平衡值，并和試驗結果保持一致。比較圖 7 的對比結果可知，實際軟組織的力學響應受到應變加載速率的影響，而 AQLV 模型對這兩種不同加載速率的大應變加載情形都能較好地預測。
雖然 AQLV 模型對大應變斜坡保持加載的上升加載階段存在一定的跟隨誤差，但誤差范圍較小；此外，該模型可以快速達到穩態且與試驗結果保持一致。所以，可以認為 AQLV 模型能夠較好地預測大應變斜坡保持試驗結果。
當加載方式為大應變正弦循環加載時，加載頻率為 0.2 Hz，加載應變幅值為 0.097 7，則應變變化范圍為 0~0.195 3。根據式（5）可以得到預測結果。大應變正弦循環加載的試驗與預測結果如圖 8 所示。因為 AQLV 模型預測結果需要經過一段時間才能達到平衡穩定值，而正弦循環加載試驗結果為經過預調制后測量得到的穩態值。因此如圖 8 所示，為 AQLV 模型穩定狀態的預測數據和正弦循環試驗數據對比，可以認為 AQLV 模型能夠較好地預測大應變正弦循環加載試驗結果。
圖8
大應變正弦循環加載預測結果
Figure8.
Predictions of the AQLV model for the sinusoidal cyclic loading under large strain
3 討論與分析
手術器械與軟組織交互作用的模擬可以通過計算機模擬仿真來實現。其中對軟組織的大變形模擬分析時，現有較高精度的模型由于存在卷積積分步驟容易導致計算響應時間長，操控響應存在滯后遲鈍等缺點；而較低精度的模型容易導致仿真誤差變大。特別地，現有仿真模型對軟組織大變形粘彈性的生物力學特性考慮不足，為了提高大變形粘彈性模型的精度，同時簡化計算過程，本文開展了對人體軟組織大變形的非線性粘彈性特性的研究探討。軟組織大變形時盡管可以采用超彈性模型來模擬加載過程的非線性變化，但是忽略了軟組織的粘彈性屬性。在經典線性粘彈性理論中，松弛函數的參數是獨立于應變的，僅為時間的函數。在 AQLV 模型中，假設模型松弛函數的參數可以根據應變水平自適應地選擇不同的參數，即假設松弛函數的參數不僅是時間的非線性函數，還是應變的非線性函數，則 AQLV 模型應該屬于一種非線性粘彈性模型，符合軟組織自身的生物力學特性。
QLV 模型假設松弛函數的參數與應變水平無關，然而 AQLV 模型假設松弛函數參數與應變水平存在函數關系。從圖 4 和表 1 試驗結果可以看到，松弛函數的參數實際上不是完全獨立于應變水平的，所以 AQLV 模型假設松弛函數模型參數不獨立于應變水平是合理的。因為 AQLV 模型假設松弛特性由并聯 Maxwell 單元來模擬，各 Maxwell 單元的時間常數設定為不變，這樣設定是為了擬合計算方便，所以假設時間常數與應變水平無關，即
、
、
分別為常數。而在實際軟組織內部，其實是不存在 Maxwell 單元這種理想化的物質。時間常數的選擇并無嚴格的標準，一般設定為相互之間存在一定的比例，比如采用對數間隔的方法[14]，然后采用設定的優化算法來比較選擇。
如圖 5 和表 2 所示的模型各個參數與應變的擬合結果來看，各個參數與應變水平采用三次多項式擬合可以得到較高精度。當應變為 0 時，假設軟組織的狀態處于零應力狀態，所以各參數均為 0。這樣假設具有一定的計算合理性，在實際試驗時候可以通過多次預調制實現軟組織的狀態穩定。著名生物力學學家馮元楨先生也提出：“要認識生物組織的改造，必須在零應力狀態來比較”[15]。
AQLV 模型的預測能力根據增量斜坡保持試驗、大應變斜坡保持試驗及大應變正弦循環加載試驗來檢驗。從圖 6、圖 7、圖 8 的預測結果來看，AQLV 模型可以較好地預測這三種加載工況的試驗結果。其中需要說明是在圖 7 中，在相同加載應變水平但不同應變加載速率下，軟組織實際穩態應力與 AQLV 模型預測的穩態應力均較好地吻合，并且應變加載速率對軟組織穩態應力值的影響很小。但是，應變加載速率對瞬態特性的影響較大，較大的應變加載速率導致軟組織實際瞬態應力變大，且 AQLV 模型對大應變斜坡保持加載的瞬態特性的擬合存在一定的誤差。這是因為 AQLV 模型不滿足疊加原理（superposition principle），Quaia 等[14, 16]的研究也證明了該模型不滿足疊加原理。
本文研究存在一定程度的局限和不足之處。AQLV 模型在不同應變速率下的預測能力還有待于后續試驗，但是計算公式可以直接代入本文所得到的公式中來擬合和預測。本文設定 AQLV 模型由 3 項松弛函數疊加組成，在手術仿真應用時還可以根據實時性的要求減少為 1 項松弛函數，即僅需擬合
和
。另外，計算機仿真的精度不僅受到軟組織非線性粘彈性模型精度影響，仿真精度和實時性還受到網格密度的影響[9]，這些工作有待于后續研究完善。
本文對手術仿真或醫學診斷中的人體軟組織特性建模進行了探討，著重對在體軟組織大變形非線性粘彈性特性進行了理論和試驗分析。以人體前臂軟組織生物力學特性為例進行印壓試驗，認為 AQLV 模型對前臂軟組織增量斜坡保持、大應變斜坡保持以及大應變正弦循環等三種加載工況下的試驗結果都能較好地預測。后續研究工作將著重對 AQLV 模型在手術仿真中的軟件設計進行研究。
