基于Hankel矩陣的磁共振波譜成像量化方法

黃敏 陳晨 陳軍波

430074武漢，中南民族大學(xué)生物醫(yī)學(xué)工程學(xué)院

基于Hankel矩陣的磁共振波譜成像量化方法

黃敏 陳晨 陳軍波

430074武漢，中南民族大學(xué)生物醫(yī)學(xué)工程學(xué)院

目的 針對磁共振波譜成像（MRSI）數(shù)據(jù)，研究基于Hankel矩陣的量化方法、抑水方法及代謝物信息圖像成像方法。 方法 分析不同Hankel矩陣結(jié)構(gòu)對幾種MRSI數(shù)據(jù)方法的影響，得到量化效率較高的Hankel矩陣結(jié)構(gòu)。應(yīng)用水成分信號強(qiáng)度最大的特點，提出基于最大幅值抑水處理MRSI方法。從MRSI數(shù)據(jù)中提取感興趣代謝物信息，再通過雙線性插值進(jìn)行代謝物信息成像。 結(jié)果 Hankel矩陣列數(shù)為量化信號點數(shù)3/4時MRSI數(shù)據(jù)幅值及頻率參數(shù)誤差達(dá)到最小，應(yīng)用3/4信號長度構(gòu)成Hankel矩陣的基于Hankel矩陣及部分重正交Lanczos算法的奇異值分解（HLSVDPRO）方法，仿真數(shù)據(jù)的幅值、頻率、衰減系數(shù)準(zhǔn)確度分別為96.94%、99.72%、95.55%。數(shù)據(jù)量化速度隨采樣點數(shù)的增加而減小，采樣點數(shù)為512點時，量化參數(shù)誤差達(dá)到最小?；谧畲蠓狄炙姆椒ㄔ诜抡鏀?shù)據(jù)的抑水程度為99.55%。 結(jié)論 將優(yōu)化的Hankel矩陣結(jié)構(gòu)應(yīng)用于基于Hankel矩陣的MRSI量化方法中對參數(shù)準(zhǔn)確度和速度均有提升，最佳采樣點數(shù)為512，基于最大幅值抑水方法對MRSI數(shù)據(jù)抑水徹底，通過對代謝物信息成像，可在（超）早期對疾病進(jìn)行診斷。

磁共振波譜成像； Hankel矩陣； 奇異值分解； 抑水； 代謝物圖

0 引 言

核磁共振現(xiàn)象以無創(chuàng)性、無輻射性等優(yōu)點被廣泛應(yīng)用于醫(yī)學(xué)檢測中，其包括磁共振成像、磁共振功能成像[1]及磁共振波譜成像（magnetic resonance spectrum imaging,MRSI）等。MRSI是一種獨特的分子成像方式，它同時提供了空間與波譜定位信息，正在成為強(qiáng)大的臨床研究工具。MRSI可非侵入性檢測代謝物濃度，能在組織或器官結(jié)構(gòu)變化之前檢測疾病。MRSI的信噪比和波譜分辨率可隨場強(qiáng)的增加而增加，因此為了獲得高信噪比和波譜分辨率，MRSI技術(shù)一般采用高場強(qiáng)，但由于其他方面技術(shù)的限制，目前MRSI可應(yīng)用1.5、3、5、7 T[2]，而在磁共振波譜（MRS）研究中，已應(yīng)用到9.4 T[3]。隨著MRSI技術(shù)研究的深入，其成像元素除H核以外，P核、C核成像也逐漸發(fā)展起來。MRSI技術(shù)應(yīng)用包括：①對腦部疾病的診斷[4]，通過N-乙酰天門冬氨酸（N-acetyl-aspartate,NAA）、肌酸（creatine,Cre）、乳酸（lactate,Lac）、膽堿（choline,Cho）、谷氨酰胺（Gln）、谷氨酸（Glu）及NAA/Cre等代謝物數(shù)值指標(biāo)以判斷患者病情，如對膠質(zhì)瘤、阿爾茲海默癥、帕金森病及腦缺血癥的診斷有很好的效果。②對前列腺癌的診斷，通過檸檬酸鹽（citrate,Cit）、Cho及Cre 3種代謝物含量，采用（Cho+Cre）/Cit值作為衡量指標(biāo)。③對乳腺癌的診斷，將組織Cho含量作為檢測標(biāo)準(zhǔn)。此外，MRSI還應(yīng)用到了宮頸癌、慢性丙型肝炎及肝硬化的診斷中。波譜量化一直是MRSI技術(shù)最耗時、并直接影響其準(zhǔn)確性的關(guān)鍵問題，除傳統(tǒng)的時域、頻域量化方法外，當(dāng)前基于主成分分析、獨立成分分析[5]、基于圓形軌跡量化[6]及時頻域結(jié)合量化[7]等新方法的研究也取得了較好的效果。

1 基于Hankel矩陣的MRSI數(shù)據(jù)時域量化方法模型

基于Hankel矩陣的MRSI時域量化方法具有數(shù)學(xué)模型完整及模型參數(shù)可完全估計的優(yōu)點。

MRSI時域信號可由Lorentzian模型進(jìn)行描述

式中：K為代謝物種類，ak、φk、dk、fk分別為代謝物的幅值、相位、衰減系數(shù)及頻率（是量化的目標(biāo)），en為噪聲。

算法實現(xiàn)過程如下。

第1步：排列yn數(shù)據(jù)點成為一個L×M的Hankel矩陣H，且N=M+L-1

式中：HK為用K截短的矩陣H，UK、VK分別為HK的左右奇異矩陣，ΣK為包含HK奇異值的對角陣。

第2步：進(jìn)行奇異值分解

式中：UK↑和UK↑分別由UK去掉最后一行和第一行所得。這一過程中，應(yīng)用最小二乘法。

第4步：ê的K個特征值即為數(shù)據(jù)的極點

第3步：估計矩陣ê

本研究討論3種基于Hankel矩陣的量化方法：基于Hankel矩陣的奇異值分解（Hankel singular value decomposition，HSVD）法[8]，其實現(xiàn)步驟與上述過程相同；基于Hankel矩陣及Lanczos算法的奇異值分解（Hankel Lanczos singular value decomposition，HLSVD）法[9]，其應(yīng)用Lanczos迭代算法進(jìn)行部分奇異值分解；基于Hankel矩陣及部分重正交Lanczos算法的奇異值分解（Hankel Lanczos with partial reorthogonalization singular value decomposition，HLSVDPRO）法[10]，即在HLSVD方法基礎(chǔ)上，將部分重正交應(yīng)用到Lanczos方法中。

2 分析與實現(xiàn)

2.1 Hankel矩陣大小和采樣點數(shù)對量化結(jié)果及耗時的影響

研究采用配置為Intel（R）Core（TM）i3-2120CPU處理器（美國Intel公司），2 G內(nèi)存的H61H2-CM PC機(jī)（清華同方股份有限公司），在MATLAB R2012a開發(fā)平臺下進(jìn)行編程。

因為HSVD、HLSVD及HLSVDPRO 3種方法必須要將信號排列成Hankel矩陣，所以首先對Hankel矩陣的結(jié)構(gòu)進(jìn)行研究，分析其對MRSI數(shù)據(jù)量化參數(shù)誤差及量化時間的影響。

根據(jù)式（1），筆者對單體素MRS的時域數(shù)據(jù)進(jìn)行仿真，仿真數(shù)據(jù)包含了3種代謝物數(shù)據(jù)及1個水信號，參數(shù)信息如表1所示。

表1 仿真數(shù)據(jù)

在場強(qiáng)為1.5 T的情況下，水信號位于4.7處，頻率為0 Hz，感興趣代謝物位于0～4，頻率為-299.28～-44.57Hz。人體MRSI數(shù)據(jù)中水信號的強(qiáng)度一般為代謝物信號強(qiáng)度的幾十甚至幾百倍，為達(dá)到該仿真目的，筆者仿真的水信號幅值為仿真代謝物幅值最大值的10倍。當(dāng)采樣頻率為1.25kHz，噪聲標(biāo)準(zhǔn)差為3dB，采樣點數(shù)為512點時，時域仿真數(shù)據(jù)如圖1所示。

圖1 磁共振波譜成像時域仿真數(shù)據(jù)

改變Hankel矩陣的大小，m為Hankel矩陣的列數(shù)，N為數(shù)據(jù)采樣點數(shù)，依次設(shè)置m為N/2、5N/8、3N/4、7N/8。分別采用HSVD、HLSVD、HLSVDPRO 3種方法對MRSI仿真數(shù)據(jù)進(jìn)行量化，得到幅值、頻率平均誤差及平均耗時（圖2）。其中，信號量化平均次數(shù)為500。

如圖2所示，Hankel矩陣大小的選取對量化結(jié)果的準(zhǔn)確度及效率均有影響。①Hankel矩陣列數(shù)為量化信號點數(shù)3/4時幅值及頻率的量化誤差達(dá)到最小。②Hankel矩陣大小對HSVD方法的量化時間影響較大，時間與Hankel矩陣列數(shù)成反比關(guān)系；但其對HLSVD和HLSVDPRO的量化時間影響很小。

由于采樣數(shù)據(jù)的點數(shù)直接影響代謝物量化的速度，因此也對仿真數(shù)據(jù)采樣點數(shù)為256點、512點、1024點的3種情況進(jìn)行了對比分析。結(jié)果表明，隨著采樣點數(shù)的增加，量化速度下降，在512點時誤差值為最?。▓D2）。綜合考慮量化時間與量化誤差，最佳采樣點數(shù)為512。

從圖2中可看出：3種方法中，HLSVDPRO方法的綜合性能最優(yōu)。利用表1中仿真數(shù)據(jù)，在采樣點數(shù)為512，噪聲標(biāo)準(zhǔn)差為3 dB時，采用HLSVDPRO方法，取矩陣長度為3/4信號長度，得到的幅值、頻率及衰減系數(shù)量化結(jié)果如表2所示。當(dāng)取50次量化平均值時，得到的幅值準(zhǔn)確度為96.94%，頻率準(zhǔn)確度為99.72%，衰減系數(shù)準(zhǔn)確度為95.55%，算法耗時為0.007 1 s。應(yīng)用式（6）計算各參數(shù)準(zhǔn)確度δ為

式中：xk為參數(shù)平均誤差值，Xk為參數(shù)準(zhǔn)確值，K為量化代謝物種類。

2.2 基于最大幅值抑水處理

在采集人體MRSI數(shù)據(jù)過程中，即便序列采用了抑水脈沖序列，仍包含非常突出的水成份信息，它們將待量化的代謝物信息淹沒，因此快速有效的抑水處理是必不可少的。

基于最大幅值抑水的處理方法，應(yīng)用水信號強(qiáng)度最大的特點，結(jié)合量化參數(shù)及水峰位于4.7附近的先驗知識移除水信號，達(dá)到抑水的目的。

圖2 Hankel矩陣大小與量化誤差及時間關(guān)系

表2 基于Hankel矩陣及部分重正交Lanczos算法的奇異值分解量化結(jié)果

由式（1）的信號模型，信號強(qiáng)度與模型中ak參數(shù)直接相關(guān)，抑水處理即轉(zhuǎn)化為抑制ak值最大的信號分量。人體MRSI數(shù)據(jù)中，水成分峰≥1個，筆者的方法是將量化得到的參數(shù)中取出最大幅值，即為，將它們?nèi)コ笾貥?gòu)得到抑水信號成分的代謝物量化結(jié)果，方法流程如圖3所示。

圖3 基于最大幅值抑水處理方法流程圖

應(yīng)用圖1中的仿真數(shù)據(jù)，采用3/4信號長度Hankel矩陣大小的HLSVDPRO方法進(jìn)行基于最大幅值抑水處理，結(jié)果如圖4所示。圖4A中藍(lán)色虛線表示仿真數(shù)據(jù)經(jīng)傅里葉變換得到的仿真頻域數(shù)據(jù)，紅線表示經(jīng)最大幅值抑水方法處理后的數(shù)據(jù)；圖4B為單獨顯示經(jīng)最大幅值抑水方法處理后的數(shù)據(jù)；圖4C為應(yīng)用式（1）及表1中水信號數(shù)據(jù)，構(gòu)造出水成分。為了計算水信號的抑制程度，建立公式（7），在水信號（即0 Hz）處，計算水信號抑制程度b為

式中：ys（nw）為經(jīng)抑水處理后信號0 Hz處值，yw（nw）為水信號0 Hz處值。圖4B、4C所示信號經(jīng)計算后，得到水信號抑制程度為99.55%。

基于最大幅值抑水方法是與MRSI量化同時進(jìn)行的，應(yīng)用量化獲得的代謝物參數(shù)進(jìn)行抑水處理，前期量化方法可變，即可移植性高；該方法的耗時只發(fā)生在重構(gòu)過程，大大減少了MRSI數(shù)據(jù)處理的總體時間；以往抑水方法只針對單個共振峰進(jìn)行，而真實人體MRSI數(shù)據(jù)中水成份往往不止一個峰，該基于最大幅值抑水方法解決了其他方法抑水不徹底的問題。

2.3 人體磁共振波譜成像數(shù)據(jù)

研究中應(yīng)用的人體前列腺MRSI數(shù)據(jù)信息如下：場強(qiáng)為1.5 T，采集體素大小為5.25 mm×4.5 mm× 4.125 mm，層間隙為4.125 mm，每一層矩陣大小為16×16，層數(shù)為16，采樣點數(shù)為512，采樣頻率為1.25 kHz，MRI數(shù)據(jù)的體素大小為 0.273 4 mm× 0.273 4 mm×4 mm，層間隙為4 mm，每一層矩陣大小為512×512，層數(shù)為20。

采集得到MRSI數(shù)據(jù)是k空間數(shù)據(jù)，首先對每一層進(jìn)行二維傅里葉變換，得到各像素處時域波譜信號，選取15種代謝物，采用HLSVDPRO方法，對該MRSI數(shù)據(jù)進(jìn)行量化處理，量化16層數(shù)據(jù)總用時386.94 s。

應(yīng)用基于最大幅值的方法對單個體素MRSI數(shù)據(jù)進(jìn)行抑水處理，抑制了4.7附近兩個水成分信號，結(jié)果如圖5所示。其中，圖5A中藍(lán)色虛線表示頻域數(shù)據(jù)，紅色線表示經(jīng)最大幅值抑水方法處理后的頻域數(shù)據(jù)；圖5B為單獨顯示經(jīng)最大幅值抑水方法處理后的頻域數(shù)據(jù)。對比圖5A與5B，可見基于最大幅值方法抑水徹底，同時對其他代謝物未產(chǎn)生影響。

MRSI技術(shù)檢測患者是否患病及定位患病部位，本研究將代謝物含量以圖像顯示，可使檢測結(jié)果直觀、高效。人體前列腺部位研究的代謝物包括檸檬酸（Cit，化學(xué)位移為2．62）、膽堿（Cho，化學(xué)位移為3．25）及肌酸（Cre，化學(xué)位移為3．05）。人體前列腺組織在病變情況下，Cho含量相對增加，大于Cit含量。

將MRSI數(shù)據(jù)在2個空間維進(jìn)行傅里葉變換，得到頻域數(shù)據(jù)點數(shù)為16×16×16×512，在每層中提取各代謝物信息，得到數(shù)據(jù)點數(shù)為16×16，再將各包含代謝物信息的數(shù)據(jù)進(jìn)行雙線性插值，得到每層數(shù)據(jù)點數(shù)為256×256，最后將Cit含量設(shè)置為紅色，Cho含量設(shè)置為綠色，Cre含量設(shè)置為藍(lán)色，進(jìn)行RGB偽彩色圖像顯示，得到圖6B，該過程共耗時0.432 6 s。

筆者應(yīng)用MRI及MRSI數(shù)據(jù)的空間采集位置信息，設(shè)置MRI與MRSI視野相同，同時，進(jìn)行MRI與MRSI層面的對應(yīng)，最后進(jìn)行雙線性插值，得到每層數(shù)據(jù)點數(shù)為256×256（圖6A）。綜合圖6A與6B中對應(yīng)層面圖像，在圖6B中具有綠色相對明顯的區(qū)域，即為病變區(qū)域，再結(jié)合圖6A中MRI圖像，可得到具有解剖學(xué)意義的診斷信息，即病變發(fā)生在前列腺外周帶。

將代謝物含量圖像與傳統(tǒng)MRI圖像結(jié)合，再加上MRI圖像分割技術(shù)[11]，可得到人體組織結(jié)構(gòu)發(fā)生變化的（超）早期病變情況。本研究對前列腺MRSI數(shù)據(jù)代謝物的成像方法簡單易行，可推廣到臨床應(yīng)用中。

3 結(jié) 論

首先從Hankel矩陣大小及采樣點數(shù)兩方面對MRSI量化方法進(jìn)行了研究，研究表明：①對基于Hankel矩陣的奇異值分解MRSI信號方法，當(dāng)Hankel矩陣構(gòu)成列數(shù)取信號長度3/4時，其量化效果最佳。②采樣點數(shù)為512時，量化效果最佳。然后提出了基于最大幅值抑水處理方法，實驗數(shù)據(jù)表明，該方法抑水徹底，其與量化處理同步進(jìn)行，不受量化方法約束，可移植性強(qiáng)。最后針對人體前列腺MRSI數(shù)據(jù)中感興趣代謝物進(jìn)行成像處理，提出了偽彩色顯示代謝物比值圖像的方法，可清晰直觀地得到病變區(qū)域。

在研究Hankel矩陣大小及信號采樣點數(shù)對MRSI量化效果影響過程中，筆者只對一種噪聲進(jìn)行了模擬，后續(xù)研究中，將模擬多種噪聲，分析其對Hankel矩陣結(jié)構(gòu)及采樣點數(shù)選取的影響；同時，將MRI與代謝物含量圖進(jìn)行融合顯示及病變位置準(zhǔn)確定位。

利益沖突 無

（圖4～6見插頁2-3）

[1]賈婧,李鏘,白煜,等.一種基于靜息態(tài)功能磁共振成像的快速丘腦分割算法[J].國際生物醫(yī)學(xué)工程雜志,2015,38(1):1-4.DOI: 10.3760/cma.j.issn.1673-4181.2015.01.001. Jia J,Li Q,Bai Y,et al.Thalamus segmentation based on restingstate functional magnetic resonance imaging[J].Int J Biomed Eng, 2015,38(1):1-4.DOI:10.3760/cma.j.issn.1673-4181.2015.01.001.

[2]Posse S,Otazo R,Dager SR,et al.MR spectroscopic imaging: principlesandrecentadvances[J].JMagnResonImaging,2013,37(6): 1301-1325.DOI:10.1002/jmri.23945.

[3]Deelchand DK,van de Moortele PF,Adriany G,et al.In vivo1H NMR spectroscopy of the human brain at 9.4 T:initial results[J].J Magn Reson,2010,206(1):74-80.DOI:10.1016/j.jmr.2010.06.006.

[4]霍珊珊,沈智威,陳耀文,等.多體素1H磁共振頻譜絕對定量腦內(nèi)代謝物的方法[J].國際生物醫(yī)學(xué)工程雜志,2011,34(3):174-178. DOI:10.3760/cma.j.issn.1673-4181.2011.03.011. Huo SS,Shen ZW,ChenYW,etal.Review ofabsolute quantification of brain metabolites by using multi-voxel1H-MRS[J]. Int J Biomed Eng,2011,34(3):174-178.DOI:10.3760/cma.j.issn. 1673-4181.2011.03.011.

[5]Ladroue C,Howe FA,Griffiths JR,et al.Independent component analysis for automated decomposition of in vivo magnetic resonance spectra[J].Magn Reson Med,2003,50(4):697-703.DOI:10.1002/ mrm.10595.

[6]Gabr RE,Ouwerkerk R,Bottomley PA.Quantifying in vivo MR spectra with circles[J].J Magn Reson,2006,179(1):152-163.DOI: 10.1016/j.jmr.2005.11.004.

[7]門闊,全紅,楊沛沛,等.基于奇異值分解和頻率選擇的磁共振波譜頻域量化方法[J].生物醫(yī)學(xué)工程學(xué)雜志,2010,27(2):249-252. DOI:1001-5515(2010)02-0249-04. Men K,Quan H,Yang PP,et al.Frequency-domain quantification based on the singular value decomposition and frequency-selection for magnetic resonance spectra[J].J BioMed Eng,2010,27(2):249-252.DOI:1001-5515(2010)02-0249-04.

[8]Vanhamme L,van Huffel S,van Hecke P,et al.Time-domain quantificationofseriesofbiomedicalmagneticresonancespectroscopy signals[J].J Magn Reson,1999,140(1):120-130.DOI:10.1006/jmre. 1999.1835.

[9]Barkhuijsen H,de Beer R,van Ormondt D.Improved algorithm for noniterative time-domain model fitting to exponentially damped magnetic resonance signals[J].J Magn Reson,1987,73(3):553-557. DOI:10.1016/0022-2364(87)90023-0.

[10]Laudadio T,Mastronardi N,Vanhamme L,et al.Improved Lanczos algorithms for blackbox MRS data quantitation[J].J Magn Reson, 2002,157(2):292-297.DOI:10.1016/jmre.2002.2593.

[11]羅清,秦文健,辜嘉,等.核磁共振圖像分割算法及應(yīng)用進(jìn)展[J].國際生物醫(yī)學(xué)工程雜志,2013,36(3):165-171.DOI:10.3760/cma. j.issn.1673-4181.2013.03.009. Luo Q,Qin WJ,Gu J,et al.Progress of the segmentation methods of magnetic resonance image and its application[J].Int J Biomed Eng, 2013,36(3):165-171.DOI:10.3760/cma.j.issn.1673-4181.2013.03. 009.

Quantification method for MRSI based on Hankel matrix

Huang Min,Chen Chen,Chen Junbo
School of Biomedical Engineering,South-Central University for Nationalities,Wuhan 430074,China Corresponding author:Huang Min,Email:minhuang@mail.scuec.edu.cn

Objective To study the quantification method based on Hankel matrix,the water suppression method and the metabolite imaging method for magnetic resonance spectroscopic imaging(MRSI)data.Methods Impact of Hankel data matrix on quantification MRSI methods were analyzed to obtain the most efficient Hankel matrix structure.The maximum amplitude method of the water signal peak was proposed for MRSI data water suppression.The interested metabolites information was extracted from MRSI data,and then metabolite image was obtained through bilinear interpolation algorithm.Results The minimum amplitude error and the minimum frequency error were acquired when columns number was 3/4 sampling points.The amplitude,frequency and the damping factor of the simulation data accuracy was 96.94%,99.72%and 95.55%respectively.Hankel lanczos with partial reorthogonalization singular value decomposition(HLSVDPRO)method with 3/4 sampling points was used to form Hankel matrices.The speed of quantification decreased with the increase of sampling points.The error of quantification parameter reached minimum when the number of sampling points was 512.The water suppression degree of simulation data was 99.55%with the maximum amplitude water suppression method.Conclusions The accuracy and the speed of the quantification are promoted with an optimized Hankel matrix structure for the MRSI quantization method.The optimal length of sampling points is 512.The maximum amplitude method can suppress water perfectly.Doctors can detect the presence of tumor regions in human body at the (super)early stage by metabolite information images.

Magnetic resonance spectrum imaging; Hankel matrix; Singular value decomposition; Water suppressed;Metabolite image

National Natural Science Foundation of China(30970782);Natural Science Foundation of Hubei Province of China(2014CFB918)

黃敏，Email:minhuang@mail.scuec.edu.cn

10．3760/cma．j．issn．1673-4181．2016．02．007

國家自然科學(xué)基金（30970782）；湖北省自然科學(xué)基金（2014CFB918）

2016-01-25）