基于BP神經網絡的左心室心肌組織參數反演方法的研究*

張其帥，楊克虎，李錄賢△

(1.西安交通大學航天航空學院 機械結構強度與振動國家重點實驗室 陜西省先進飛行器服役環境與控制重點實驗室，西安 710049；2.西安電子科技大學通信工程學院 綜合業務網理論與關鍵技術國家重點實驗室，西安 710071)

1 引 言

近年來，隨著心臟力學研究的不斷深入，人體心肌的組織特性逐漸成為表征心臟健康狀態的一個重要指標[1-2]。研究表明，相比于健康者，舒張型心力衰竭患者或心肌梗死患者的心肌硬度要高很多[3-5]，因此，開展心肌組織特性研究對心臟疾病的臨床診斷具有重要意義。

限于人體心肌組織標本很難獲得[6-7]，研究者對心肌組織特性的研究主要依賴于動物實驗模型，即假設心肌組織特性符合一種特定本構關系，通過對標本的單軸拉伸、雙軸拉伸或剪切試驗獲取應力應變曲線，確定本構方程中的材料參數[8-10]。此方法雖可定量分析心肌的組織特性，但離體失活的心肌測量結果與具有血流特性的真實心肌特性之間存在明顯的差異[3]，加之動物心肌組織與人體心肌組織存在差異，其結果的可參考性也一直存在爭議[5]。

針對獲取心肌組織特性存在的困難，本研究提出一種基于智能優化算法的反演方法。采用左心室動態響應數據建立參數反演的BP(back propagation,BP)神經網絡結構，映射左心室動態響應與心肌組織參數之間的非線性關系；將個體患者的臨床非侵入式左心室診斷數據輸入到已訓練好的BP神經網絡，反演出測試者的心肌組織參數。該基于BP神經網絡的反演方法在獲得心肌組織特性的同時，還能滿足臨床心臟疾病檢測的快速時間響應需求。

2 左心室的有限元建模

2.1 左心室的幾何模型

本研究中的幾何模型來自人體左心室(left ventricle，LV)的內外膜CT掃描數據。首先沿左心室軸向進行等間距(5 mm)掃描，獲得12個不同截面的左心室CT影像，見圖1(a)—(l)。

圖1 左心室的CT截面影像Fig.1 CT Image of LV sections

再利用Matlab語言編寫圖像識別程序，提取CT影像中的左心室內外膜位置點，分別見圖2(a)、(b)。

圖2 Matlab圖像識別的左心室CT數據Fig.2 The LV CT Data recognized by Matlab

進一步，通過SolidWorks軟件，建立左心室的幾何模型，見圖3。

2.2 左心室的有限元模型

經計算圖1中12個不同截面的CT影像，左心室心肌的平均厚度為5.65 mm，第5～12層共8個截面的平均厚度在6 mm左右。基于圖3的左心室幾何模型，抽取其中面，形成左心室的殼結構，取6 mm作為殼的厚度。利用Abaqus軟件自動劃分網格功能，控制左心室模型邊界布種密度和網格屬性，得到圖4所示的殼結構有限元網格，共計4 084個S4R殼單元和60個S3R殼單元，節點總數4 150個。

圖3 基于CT影像的左心室幾何模型Fig.3 CT image based geometric model of the LV

圖4 左心室有限元模型Fig.4 Finite element model of the LV

3 BP反演數據的建立方法

BP神經網絡可以映射任意復雜的非線性函數關系[11-12]，本研究運用BP神經網絡技術反演心肌組織的性能參數，為此，我們對左心室先進行動態有限元分析，以獲得BP神經網絡的訓練數據。

3.1 心肌組織特性

人體心肌是一種具有特異性質的超彈性不可壓天然生物材料。從心肌性能的不同特點出發，不少學者建立了心肌組織性能的描述方法，但至今仍沒有一種被廣泛認可的表達形式。本研究聚焦于反演方法研究，將心肌視為各向同性超彈性材料，對于實際的各向異性心肌組織，在增加參數數目和反演運算量后，本研究方法仍將適用。

對于各向同性超彈性材料，我們選取最常用的不可壓Mooney-Rivlin模型[13]，其能量密度函數的表達式為：

(1)

3.2 受載分析

心臟的周期性收縮與擴張來源于兩個因素，一個是心肌受到心電刺激所產生的主動力，另一個是心臟變形受到內部血液壓力反作用而產生的被動力[14]。主動力可通過心肌興奮力σ隨興奮時間t的變化表示為[14]：

(2)

其中τ為不同時刻的傳播時序，Te為興奮周期，σmax為肌小節興奮力。本研究取與文獻[15]相同的參數，即σmax=53.3 kPa和Te=0.353 s。興奮時序τ一般較難確定[15]，本研究取為0。左心室被動力通過動脈壓力予以間接獲得[7, 15]。一個心動周期內左心室的受載變化見圖5。

3.3 BP輸入數據的獲取

給定心肌組織參數(也就是BP神經網絡的目標向量)，對左心室有限元模型心底位置施加固支約束，對內外表面分別施加圖5中均布的徑向被動力和周向主動力，利用Abaqus軟件的Explicit求解器，對左心室一個心動周期0.8 s內的變形過程進行分析[15]，得到圖6所示的變形過程。可以看出，在一個心動周期內的不同時刻，左心室發生了明顯的收縮變形，與真實的左心室運動過程十分類似。

圖5 左心室受載隨時間變化關系[15-16]Fig.5 Variation of loading with time for the LV[15-16]

圖6 收縮過程中左心室幾何形態變化Fig.6 Geometric variation of the LV during a period of contraction

通過左心室的有限元分析，還可得到左心室任意位置的變形。本研究選取左心室Z=8、20、32、44 mm的四個典型截面位置(見圖7)，將心動周期內某些時刻的徑向位移值作為BP神經網格的輸入向量。

圖7 取樣點在左心室結構中的位置分布Fig.7 Distribution of sampling points in LV structure

4 基于BP神經網絡的心肌組織參數反演

4.1 BP神經網格反演的基本思想

BP神經網絡是一種按誤差逆傳播算法訓練的多層前饋網絡，是目前應用最廣泛的神經網絡模型之一[17]。BP神經網絡由輸入層、隱含層和輸出層組成，在不限制隱含層節點的情況下，三層的神經網絡理論上可逼近任意非線性函數關系式[18]，見圖8。

圖8 BP神經網絡拓撲結構Fig.8 The topologic structure of BP neural network

基于BP神經網絡的參數反演就是建立響應與參數之間的非線性關系，包括以下基本步驟[19]：

(1)確定待反演參數的取值范圍；

(2)在參數取值范圍內均勻分散取值獲得目標向量；

(3)通過響應獲得輸入向量；

(4)將輸入、目標向量輸入BP神經網絡進行訓練；

(5)利用BP神經網絡輸入和輸出的映射關系，將人體左心室診斷數據輸入神經網絡反演參數。

上述5個步驟又可概括為BP神經網絡的建立(包括步驟(1)—(4))和參數反演(即步驟(5))兩部分，下面予以分別研究。

4.2 BP神經網絡的建立

4.2.1輸入向量和目標向量的建立 由文獻[20]可知，以Mooney-Rivlin模型表示的正常左心室心肌組織參數分別為C10=0.178 MPa和C01=3.989 MPa，密度為1.370×103kg/m3。據此，兩個參數的可能取值范圍分別確定為0.100≤C10≤0.300和3.000≤C01≤5.000；基于正交試驗設計方法，對兩個參數進行組合，取值C10=0.100、0.150、0.200、0.250、0.300，C01=3.000、3.250、3.500、3.750、4.000、4.250、4.500、4.750、5.000，共形成45組不同的心肌組織參數組合，將它們作為BP神經網絡訓練的目標向量。

取心動周期內t=0.2 s(收縮初期)、t=0.4 s(收縮中期)和t=0.6 s(收縮末期)共3個特征時刻左心室4個典型截面上的徑向位移作為BP神經網絡的輸入向量。每個截面上等間隔地取12個點，共計48個點。不同心肌組織參數下的左心室模型對應不同的徑向位移，形成45組訓練樣本。在3個特征時刻下，得到對應時刻下的45組訓練樣本，分別見表1、表2、表3。

表1 不同心肌組織參數對應的左心室徑向位移(t=0.2 s)Table 1 Radial displacements at the LV from different myocardial tissue parameters (at t=0.2 s)

表2 不同心肌組織參數對應的左心室徑向位移(t=0.4 s)Table 2 Radial displacements at the LV from different myocardial tissue parameters (at t=0.4 s)

4.2.2隱含層和激活函數 隱含層的選擇根據經驗和多次實驗來確定，本研究中輸入單元數為48，輸出單元數為2，根據經驗公式計算隱含層單元數的范圍為[7, 18]，再通過BP神經算法訓練樣本，綜合考慮收斂時的誤差和迭代次數，確定最佳的隱含層單元數為15。

輸入層和隱含層的激活函數使用雙曲型正切tansig函數，輸出層使用purelin線性函數，訓練函數用負梯度下降動量BP算法。接下來，使用Matlab軟件編制BP神經網絡程序并運行，得到心肌參數反演的BP神經網絡結構，見圖9。

圖9 Matlab中建立的BP神經網絡結構Fig.9 BP neural network model developed via Matlab software

表3 不同心肌組織參數對應的左心室徑向位移(t=0.6 s)Table 3 Radial displacements at the LV from different myocardial tissue parameters (at t=0.6 s)

4.3 心肌參數反演結果及討論

臨床上通過左心室標記點方法可獲得心肌在不同位置點的徑向位移。因此，為驗證本研究方法，假定心肌組織性能參數為已知，據此通過有限元分析，獲得左心室結構的響應數據，再通過BP神經網絡對響應數據的反演，獲得心肌組織參數的大小，并與假定的心肌組織參數進行誤差比較。通過對4個假定實例的心肌組織參數反演得到在3個特征時刻下的反演參數的平均值，反演參數的平均值與真實值間的誤差見表4。

由表4可以看出，4個實例的反演值平均值與真實值間的誤差在3.58%以內，整體L2范數相對誤差為0.72%，與臨床醫學診斷過程中的常見誤差[21]相比，本研究參數反演的精度較高，可適合臨床診斷的要求，同時說明本研究建立的BP神經網絡具有高可靠性。

表4 真實心肌組織參數和反演心肌組織參數對比Table 4 The inversed results and their errors

5 結論

基于臨床CT圖像，建立了左心室的真實三維幾何模型；根據左心室的動態有限元分析結果，提出了基于BP神經網絡的左心室心肌組織參數的反演方法，獲得了心肌組織Mooney-Rivlin超彈性模型下的心肌組織參數，反演值平均值與真實值具有很好的一致性，證明了本研究心肌組織參數反演方法的正確性，可為人體心肌組織性能的準確獲得提供一種有效方法。

采用臨床上獲得的心臟患者左心室標記點變形數據，將此輸入到本研究訓練好的BP神經網絡，反演得到心肌組織參數，可作為臨床心臟疾病診斷的依據。本研究BP神經網絡算法為左心室表象和心肌組織本征特性之間建立了定量關系，可定量描述心臟的健康程度，從而為臨床心臟疾病的高效高精確診斷提供了一條可行途徑。