金杰鋒劉華鋒*胡紅杰

1(浙江大學現(xiàn)代光學儀器國家重點實驗室，杭州 310027)2(浙江大學醫(yī)學院附屬邵逸夫醫(yī)院，杭州 310016)

基于MRI圖像的心肌三維彈性成像

金杰鋒1劉華鋒1*胡紅杰2

1(浙江大學現(xiàn)代光學儀器國家重點實驗室，杭州 310027)2(浙江大學醫(yī)學院附屬邵逸夫醫(yī)院，杭州 310016)

材料參數(shù)比起運動參數(shù)能夠更顯著地反映物體的內(nèi)在屬性變化，定量獲得病變組織的材料參數(shù)對心臟疾病診斷具有非常重要的意義。本研究提出一種在利用MRI圖像獲得部分形變位移等運動信息和約束條件的情況下，依據(jù)連續(xù)力學模型，使用有限元分析結(jié)合H∞濾波估計方法實現(xiàn)三維彈性成像的算法。模擬數(shù)據(jù)仿真給出了與最小二乘法結(jié)果的對比，驗證了方法的有效性。最后給出了對病人心臟MRI圖像三維彈性成像結(jié)果。

左心室;三維彈性成像;有限元方法;H∞濾波

Abstract:Quantification of organic tissue’s material parameters has significant implication in clinic diagnosis of cardiac diseases，since material parameters are far more sensitive measures of the object intrinsic property changes.With known nodal kinematic data and boundary constrains，a new 3D elastic modulus estimation algorithm was based on finite-element method and H∞ filter.Simulation results proved the efficiency of the method compared with the least square method.A myocardiac 3D elastography using cardiac MRI images was presented as well.

Key words:left ventricle;3D elastography;finite elements method;H∞filter

引言

心臟是人體最重要的器官之一，它運作的方式和機理復(fù)雜而精巧，對心臟的研究一直是醫(yī)學上一個非常重要的領(lǐng)域。早期認識心臟主要是根據(jù)解剖學的知識，由于技術(shù)手段的限制人們無法了解其在人體內(nèi)部工作的細節(jié)。隨著現(xiàn)代醫(yī)學影像技術(shù)的發(fā)展，人們已經(jīng)可以無創(chuàng)地觀察人體內(nèi)部組織器官，利用CT、MRI等手段，不僅可以獲得二維的人體斷層圖像，而且可以利用先進的計算機技術(shù)生成心臟的三維模型，更加有效地對心臟的運作規(guī)律進行研究，從而將對心臟疾病的研究和診斷提升到一個新的高度。

很多研究表明，生物組織的病變會引起彈性系數(shù)的變化，尤其是對于像心臟這樣由肌肉組織構(gòu)成的器官。因此在臨床上，對心肌組織的材料系數(shù)進行量化的測定有著非常重要的實際意義。目前，超聲應(yīng)變成像與彈性成像［1］以及 MR彈性成像(MRE)［2］都已經(jīng)在心臟運動研究領(lǐng)域得到了應(yīng)用。但超聲應(yīng)變(應(yīng)變率)成像只能提供一維信息，而且受到操作手法、成像角度等因素影響。MRE成像需要外部激發(fā)裝置，目前只能成平面圖像，而且現(xiàn)有的圖像處理方法也存在分辨率不高或者易受噪聲影響等缺點。

連續(xù)生物力學模型將通過實驗得到的心臟材料屬性先驗知識引入到心臟圖像分析中，取得了比較好的效果，很多研究者在研究中采用了這種方法［3-5］。利用連續(xù)生物力學模型研究心臟器官的材料屬性的方法是假設(shè)心臟的運動信息已知(比如依靠植入種子獲得或者通過圖像處理的方法測量運動信息，如位移值或瞬時速度值)，然后從這些運動信息中來估計材料參數(shù)如楊氏模量、泊松比等，由此可以實現(xiàn)對材料的成像。

以往基于醫(yī)學圖像信息求解組織的彈性模量分布信息的工作，通常只對某一斷層圖像進行單獨處理［6］，忽略了垂直于圖像平面方向上的運動信息，且所采用的最小二乘法重建無法應(yīng)用于更復(fù)雜的三維數(shù)據(jù)，容易受到噪聲的影響。隨著三維成像和基于圖像的三維輪廓運動追蹤方法的不斷發(fā)展［7-8］，已經(jīng)可以從多層斷層圖像中提取三維運動信息，利用三維運動信息估計整個心肌的彈性模量，從而更準確地反映實際情況。因此，本研究提出了一種基于有限元與H∞濾波估計的方法，利用MRI心臟圖像信息對心肌實現(xiàn)三維彈性成像，并通過實驗證明了其有效性。

1 方法

1.1 心肌的有限元表達

使用有限元方法來表達心肌求解域。有限元方法是一種在工程問題中使用非常廣泛的近似數(shù)值解法，其基本思想將求解域分解為較小單元來處理，在三維問題中，常見的單元類型是四面體單元和六面體單元。由于四面體單元對于復(fù)雜形狀的網(wǎng)格化有更好的適應(yīng)性，網(wǎng)格劃分算法也比較成熟，而且計算量小，因此本研究采用這種單元形狀。對左心室劃分網(wǎng)格的方法是，先對MRI等斷層成像圖像使用圖像分割方法獲得左心室心肌表層的內(nèi)外輪廓邊界，并在內(nèi)外輪廓線界定的區(qū)域均勻布點，然后使用三維Delaunay四面體剖分算法生成所需的四面體網(wǎng)格。

當四面體單元形變時，四個頂點都有沿 x，y，z三個方向的位移，結(jié)合起來用一個矩陣表示:

四面體單元內(nèi)的位移u可以僅用其四個頂點的位移ue來表示

其中 N為四面體單元形狀函數(shù)矩陣，具體形式:

其中 Ni，Nj，Nk，Nm分別是四面體單元四個頂點對應(yīng)的形狀函數(shù)，可以由下式計算得到:

記由四面體單元四個頂點坐標構(gòu)成的矩陣為

則ap，n為矩陣 Λ 的第 p行，第 n列元素的代數(shù)余子式，單元體積Ve=|Λ|/6。

1.2 心肌三維連續(xù)生物力學模型

實際心肌的力學屬性非常復(fù)雜，從計算的角度考慮，假設(shè)心臟是各向同性的線性彈性體，其應(yīng)力向量σ與應(yīng)變向量ε成線性關(guān)系:

式中，材料矩陣

式中，E是彈性模量，表示物體產(chǎn)生形變的難易度;ν是泊松比，表示橫向應(yīng)變與縱向應(yīng)變之比。這兩項都屬于物體的材料性質(zhì)，實際中泊松比的變化較小，可以認為是常量。

空間中的應(yīng)變與位移的關(guān)系由三維幾何方程表示式中，ρ是心肌組織平均密度，I是12×12的單位矩陣。

1.3 從圖像獲取運動信息

建立序列圖像不同時間幀之間像素點的對應(yīng)關(guān)系一直是計算機視覺領(lǐng)域的重要研究課題。在心臟圖像研究領(lǐng)域中，對于普通 MRI圖像，所能依據(jù)的僅僅是心臟內(nèi)外邊界的匹配。一種方法是基于幾何形態(tài)標記匹配和定位的心肌壁輪廓運動追蹤方法［7］。該方法基于左心室輪廓形變在連續(xù)的兩個時間幀之間盡可能小的假設(shè)，利用曲率能量作為匹配原則，來得到節(jié)點在前后兩個時間幀中的匹配關(guān)系。如果利用MRI標記成像和MRI相對比成像技術(shù)，那就可以匹配更多像素點，得到更豐富的運動信息，使彈性成像的結(jié)果更好。

在獲取左心室內(nèi)外壁三維空間位移場時，本研究利用文獻［12］中提出的心臟左心室形狀和運動聯(lián)合估計框架，其主要思想是通過使內(nèi)能與外能達到平衡獲得幀與幀之間的邊界點的對應(yīng)關(guān)系，其中內(nèi)能由采用點的力學特性所決定，外能為邊界力、先驗力、時間相關(guān)力和保形力的組合。

1.4 重建算法

基于有限元方法原理得到心臟運動的系統(tǒng)動

由式(2)和式(8)可以得到:

由此可見四面體單元內(nèi)各點的應(yīng)變是一樣的，是常應(yīng)變單元，應(yīng)變大小由四個節(jié)點的位移確定，B是單元應(yīng)變矩陣。四面體單元剛度矩陣的計算如下式:

四面體單元的質(zhì)量矩陣采用集中質(zhì)量矩陣的形式，即認為質(zhì)量集中在單元的節(jié)點上:態(tài)方程的矩陣形式如下［9］

式中，R是系統(tǒng)所有節(jié)點的載荷向量，反映物體所受作用力的大小和分布;U是系統(tǒng)所有節(jié)點的位移向量，其一階微分和二階微分分別表示速度向量和加速度向量。M和K分別是系統(tǒng)的質(zhì)量矩陣和剛度矩陣，均利用相應(yīng)的單元矩陣(10)和(11)通過有限元方法構(gòu)造而成。C為阻尼矩陣，假設(shè)為瑞利阻尼，滿足公式

α，β均取 0.01。

本研究的目標是重建彈性模量，即通過重組系統(tǒng)動態(tài)方程(式(12))，利用已有的位移信息和邊界力信息來求得彈性模量E的分布。因此我們需要把式(12)的左邊變換成 E向量的函數(shù)，即 KU=G1E，K=G2E，由方程(7)可知這樣的變換是可行的0。Ke可以表示成含Ee項的形式:

由此，根據(jù)剛度矩陣構(gòu)造的方法，可以得到構(gòu)造G1和 G2的方法:

1)初始化兩個N×Ne的空矩陣 G1和 G2，N為系統(tǒng)節(jié)點變量數(shù)，在三維空間中，每個節(jié)點有三個方向的自由度，即系統(tǒng)總節(jié)點數(shù)×3，Ne為系統(tǒng)總單元數(shù)。

3)對于一個編號為j的單元，利用方程(14)構(gòu)造不包含該單元彈性模量Ej的剛度矩陣。

4)按照系統(tǒng)生成網(wǎng)格時單元和節(jié)點的編號規(guī)則，把K′j中的元素插入系統(tǒng)臨時剛度矩陣中相應(yīng)的位置中，并將其單元內(nèi)下標改成相應(yīng)的系統(tǒng)下標。

6)返回到第2)步直到所有節(jié)點遍歷。

由此構(gòu)造出 N×Ne的矩陣 G1，G2。式(12)經(jīng)過變形:

為了求解E，將式(15)整理成如下

令 G=βG2+G1，R′=R-Mü- αM，得到簡化形式:

這樣利用一定的邊界條件，系統(tǒng)彈性模量向量E就能直接利用最小二乘法得到近似的解:

在三維空間中，單元的數(shù)量常常大于節(jié)點自由度。

1.5 H∞濾波估計

由成像手段獲得的運動信息必然含有一定噪聲，最小二乘法不能很好的消除噪聲帶來的影響，本研究采用一種基于H∞的噪聲擾動全狀態(tài)信息方法(NPFSI)來實現(xiàn)估計過程。

為了應(yīng)用H∞濾波估計方法，首先必須建立系統(tǒng)狀態(tài)方程。設(shè)狀態(tài)向量為x(t)=［U(t)(t)］T，控制向量為w(t)=［0 R］T，將式(11)變形，建立連續(xù)時間線性隨機系統(tǒng)狀態(tài)空間表達式:

式中，D為測量矩陣，D的行數(shù)為三維運動信息可以測量得到的有限元節(jié)點自由度，D的列數(shù)為系統(tǒng)總節(jié)點自由度。D的每行中，索引位置對應(yīng)的有限元節(jié)點的三維運動信息如果可以測量得到，該索引位置的元素設(shè)為1，否則設(shè)為0。v(t)代表過程噪聲，e(t)代表觀測噪聲，本研究均假設(shè)為0均值的白噪聲(E(v(t))=0，E(v(t)v(s)′)=Qv，E(e(t))=0，E(e(t)e(s)’)=Re)。

由于彈性模量E的是包含在參數(shù)矩陣中的，為了實現(xiàn)估計的目的，設(shè)材料向量為θ=［E］，再將原運動狀態(tài)量x和材料參數(shù)θ組成一個新的狀態(tài)向量，根據(jù)式(15)將原狀態(tài)方程變形，得到如下一組新的狀態(tài)方程，這里假設(shè)材料參數(shù)θ不隨時間變化

根據(jù)文獻［11］由狀態(tài)方程得到如下一組連續(xù)狀態(tài)估計方程:

為了使估計算法達到預(yù)期的效果，在開始估計之前，必須先設(shè)定合適的噪聲干擾衰減水平，權(quán)重矩陣 Q，Q0，Q1和 γ 值。通常很難決定最優(yōu)的 γ*，它可能是無窮的，意味著無法達到期望的噪聲控制水平。然而，實際上如果Q被選定了，γ*能夠通過分析獲得［12］。

在計算過程中，彈性模量E的初值可以根據(jù)對心肌的經(jīng)驗知識，設(shè)定與所成像對象的真實彈性模量接近的某一個先驗的值，然后利用式(25)～(28)對預(yù)設(shè)的彈性模量分布E進行遞歸估計修正，直到獲得最優(yōu)解。

2 實驗與討論

為了驗證本方法，首先采用了類似心臟左心室結(jié)構(gòu)的仿真模型進行了實驗。仿真模型高100 mm，內(nèi)徑60 mm，外徑100 mm，深灰色部分表示正常組織，彈性模量取正常心肌組織的先驗值75 kPa，淺灰色代表較硬組織，彈性模量為105 kPa，黑色代表較軟組織，彈性模量為45 kPa，其形變狀態(tài)與位移場分布如圖2所示。

將模型的底面固定，頂面施加1 000 N垂直向下的壓力，使用 ABAQUS軟件生成節(jié)點位移數(shù)據(jù)，作為用于重建的理想數(shù)據(jù)。

從圖2中可以看到不同材料之間的位移差別非常小。這一方面說明從位移數(shù)據(jù)中重建物體的材料參數(shù)非常困難，尤其當數(shù)據(jù)包含噪聲時。另一方面，如果在這樣的條件下可靠的重建出了物體的彈性模量分布，那么這說明物體的材料參數(shù)能夠比運動參數(shù)更顯著地反映物體屬性細微變化。

圖3是在不同噪聲程度影響下，使用H∞濾波估計方法得到三維可視化結(jié)果。

圖1 三維仿真模型Fig.1 3D synthetic model

圖2 施加力后物體的形變狀態(tài)和位移場。(a)形變狀態(tài);(b)單元節(jié)點位移向量Fig.2 Deformed state and displacement field after applying forces.(a)Deformedstate;(b)Nodal displacement field.

圖3 不同噪聲情況H∞濾波估計三維可視化結(jié)果(左列為外側(cè)正面視圖，右列為內(nèi)側(cè)反向視圖)。(a)和(b)不含噪聲;(c)和(d)SNR=40 dB;(e)和(f)SNR=30 dBFig.3 Visual results of reconstruction from the synthetic data corrupted by noise using H∞filter algorithm.(The left column is front view，the right column is semi-sectional view from opposite direction)(a)and(b)noise free;(c)and(d)SNR=40 dB;(e)and(f)SNR=30 dB

表1是加入不同程度和種類的噪聲后重建的結(jié)果，可以發(fā)現(xiàn)信噪比在30 dB的水平時，最小二乘法已經(jīng)不能重建出真實分布，但使用H∞濾波估計仍然能得到可以接受的接近真實彈性模量分布結(jié)果，而且對不同類型的噪聲都有很好的適應(yīng)性。

真實心臟三維彈性成像實驗采用來自醫(yī)院病人心臟數(shù)據(jù)，整個心臟分10個斷層，每層原始圖像大小512像素(512像素，層內(nèi)分辨率0.704 4 mm/像素，層間隔為8 mm，整個心臟運動周期內(nèi)包含20幀，從心臟舒張末期開始。

為了獲得心肌三維運動信息量，本研究利用文獻［12］中提出的心臟左心室形狀和運動聯(lián)合估計框架，其主要思想是利用連續(xù)生物力學心臟有限元模型(如圖4所示)，每個采樣點的受力表示為邊界力、先驗力、時間相關(guān)力和保形力的組合，從而同時求得左心室心肌的運動軌跡和序列圖像分割的結(jié)果。使用該框架得到的分割結(jié)果如圖5所示。圖6是得到左心室舒張末期到收縮末期的位移值。

表1 H∞濾波重建和LS重建在不同噪聲情況下的結(jié)果對比 (均值±標準差，kPa)Tab.1 Results of reconstruction using H∞and LS from the synthetic data corrupted by different kinds and degrees of noise(mean±SD，kPa)

圖4 心臟左心室有限元模型Fig.4 Left ventricular finite elements model

圖5 左心室MRI圖像序列的分割結(jié)果。(a)舒張末期 (b)收縮末期Fig.5 Segmentation of LV from cardiac MR image sequence(a)end-diastole(ED);(b)end-systole(ES).

圖6 左心室從舒張末期到收縮末期的三維位移信息Fig.6 The displacement field of LV model between ED and ES

該框架中加入了幾何形態(tài)匹配的約束條件，因此對于左心室圖像的分割和心肌壁上節(jié)點的運動跟蹤能夠得到比較好的效果，但由于其框架假設(shè)心肌材料是均勻一致的，而且所采用的作用力僅是來自圖像信息的虛擬力，所以對于壁間節(jié)點的運動重建與實際有一定的差異。因此僅采用其框架得到的內(nèi)外心肌壁上的節(jié)點舒張末期至收縮末期的圖像幀的三維空間位移值，作為觀測信息即式(19)的觀測值y。由于沒有相應(yīng)的相對比MRI圖像，而心臟舒張末期和收縮末期速度值比較小，y中的速度觀測量可以設(shè)為0。采用病人血壓的收縮壓和舒張壓之差作為心肌內(nèi)壁上的壓力載荷R。利用本算法進行估計，最終得到彈性成像的結(jié)果圖7所示。從圖6可以看出，上部分心肌在同樣壓力作用下產(chǎn)生的位移較小，其彈性模量應(yīng)該比較大，圖7的結(jié)果正好反映了這一點，說明本算法得到的結(jié)果是符合實際意義的。

圖7 左心室三維彈性成像結(jié)果Fig.7 Visual results of LV elastic modulus distribution

3 結(jié)論

本研究提出了利用MRI圖像信息，使用有限元結(jié)合H∞的濾波估計方法來進行心肌組織的三維彈性成像。仿真實驗表明算法在有噪聲干擾的情況下，仍可以估計出接近真實的結(jié)果，驗證了本方法的有效性。進一步的實驗利用真實病人心臟MRI圖像數(shù)據(jù)得到了彈性成像結(jié)果。未來在對本方法加以一定的完善和實用化之后，可以利用MRI圖像數(shù)據(jù)給醫(yī)生提供輔助診斷所需的組織材料信息，對于心臟疾病的診斷具有重要參考意義。

［1］Konofagou E，D'hooge J，Ophir J.Cardiac elastography-A feasibility study［J］.2000 IEEE Ultrasonics Symposium Proceedings，2000，1-2:1273-1276.

［2］Elgeti T.Cardiac magnetic resonance elastography initial results［J］.Investigative Radiology，2008;43(11):762-772.

［3］Frangi AF，Niessen WJ，Viergever MA.Three-dimensional modeling for functional analysis of cardiac images:A review［J］.IEEE Transactions on Medical Imaging，2001，20(1):2-25.

［4］Papademetris X，Sinusas AJ，Dione DP，Duncan JS.Estimation of 3D left ventricular deformation from echocardiography［J］.Medical Image Analysis，2001，5(1):17-28.

［5］Shi Pengcheng，Sinusas AJ，Constable RT，et al.Volumetric deformation analysis using mechanics-based data fusion:Applications in cardiac motion recovery［J］.International Journal of Computer Vision，1999，35(1):87-107.

［6］Zhu Yanning，HallTJ，Jiang Jingfeng.A finite-element approach forYoung's modulus reconstruction［J］.IEEE Transactions on Medical Imaging.2003，22(7):890-901.

［7］Shi Pengcheng，Sinusas AJ，Constable RT，et al.Point-tracked quantitative analysis of left ventricular surface motion from 3-D image sequences［J］.IEEE Transactions on Medical Imaging，2000，19(1):36-50.

［8］莊玲，劉華鋒，鮑虎軍.基于圖像的心臟左心室形狀和運動參數(shù)聯(lián)合估計框架的研究［J］.自然科學進展，2007，17(3):396-400.

［9］ShiPengcheng，Liu Huafeng.Stochastic finite element framework forsimultaneousestimation ofcardiackinematic functions and material parameters［J］.Medical Image Analysis，2003，7(4):445-464.

［10］Shi Pengcheng，Liu Huafeng.Robust identification of object elasticity［J］.Computer Vision and Mathematical Methods in Medical and Biomedical Image Analysis，2004，3117:423-435.

［11］Didinsky G，Pan ZG，Basar T.Parameter-identification for uncertain plants using H-infinity methods［J］.Automatica，1995，31(9):1227-1250.

［12］Didinsky G，Basar T，Bernhard P.Structural-properties of minimax controllers for a class of differential-games arising in nonlinear H-infinity control［J］.Systems ＆ Control Letters，1993，21(6):433-441.

Myocardiac 3D Elastography Using MRI Images

JIN Jie-Feng1LIU Hua-Feng1＊HU Hong-Jie2
1(State Key Laboratory of Modern Optical Instrumentation，Zhejiang University，Hangzhou 310027，China)2(The Sir Run Run Shaw Hospital，College of Medical Sciences Zhejiang University，Hangzhou 310016，China)

TP391.41

A

0258-8021(2010)04-0564-07

10.3969/j.issn.0258-8021.2010.04.014

2009-11-20，

2010-02-20

國家自然科學基金資助項目(60872068);浙江省自然基金項目(R207119)

*通訊作者。 E-mail:liuhf@zju.edu.cn