基于奇異值分解方法的磁感應(yīng)磁聲斷層重建算法*

馬任，楊銘，張順起，周曉青，殷濤，劉志朋

（北京協(xié)和醫(yī)學院中國醫(yī)學科學院生物醫(yī)學工程研究所，天津300192）

1 引 言

醫(yī)學相關(guān)研究表明，癌癥的早期診斷可以提高治愈率[1]。通過醫(yī)學成像方式，例如X射線斷層成像、核磁共振成像（MRI）及正電子發(fā)射斷層成像（PET），可以進行癌癥的診斷，目前X射線斷層成像及核磁共振成像的早期診斷精確度不高，PET雖然可以進行癌癥的早期診斷，但其具有放射性，且檢查價格昂貴，并不適用于癌癥早期篩查及診斷。當人體組織發(fā)生癌變時，癌變部位的電特性會與其他部位不同，癌變組織和正常組織之間電特性參量也會存在明顯差異，科研人員通過離體組織電特性測量實驗證明在癌變的早期，癌變組織的電導(dǎo)率約為正常組織的6.4倍，而介電常數(shù)約為正常組織的3.8倍，最高可相差高達一個數(shù)量級[2]。電特性的差異可以作為成像對比度，從而實現(xiàn)對癌癥的早期診斷。目前，可以實現(xiàn)檢測這種電特性差異的方法有采用電磁場直接檢測方法，包括電阻抗成像[3]及磁感應(yīng)成像[4]。同時，為了提高成像分辨率，研究人員提出采用電磁場及聲場互相耦合的新型電阻抗成像方式，這些方法主要有霍爾效應(yīng)成像[5-6]，磁聲成像[7-8]以及磁感應(yīng)磁聲成像[9-10]等。

磁感應(yīng)磁聲成像（magnetoacoustic tomography with magnetic induction，MAT-MI）是將阻抗成像與超聲檢測相結(jié)合，利用超聲檢測分辨率高的優(yōu)點，將組織內(nèi)部電導(dǎo)率通過磁聲耦合效應(yīng)轉(zhuǎn)化為超聲信號傳出體外并重建組織電導(dǎo)率分布。近年來，由于激勵源特性、信號處理方法與重建算法及理論的改進，偶極子源[11]，矢量源重建算法[12]及基于真實聲換能器特性的重建算法[13-14]相繼被應(yīng)用于磁聲成像中，磁聲成像的實驗仿體也逐漸從高電導(dǎo)率的銅環(huán)、石墨等金屬擴展到低電導(dǎo)率的類生物組織明膠瓊脂仿體，隨后離體豬肉，肝臟也可進行粗略成像，電導(dǎo)率圖像質(zhì)量已經(jīng)能分辨出脂肪和肌肉的邊界輪廓，但是目前實驗結(jié)果尚不能利用檢測聲信號重建非邊界位置的內(nèi)部電導(dǎo)率信息。為了研究MAT-MI中電導(dǎo)率分布與檢測超聲信號之間的關(guān)聯(lián)，實現(xiàn)電導(dǎo)率定量重建，本研究基于MAT-MI原理，建立了一個系統(tǒng)矩陣，通過截斷奇異值方法（TSVD）對MATMI算法進行仿真研究。

2 理論分析

MAT-MI是從待成像組織的電阻抗分布到檢測到的超聲信號的轉(zhuǎn)換過程。其基本原理是將待檢測組織置于穩(wěn)恒磁場B中，同時外加相同方向的脈沖變化磁場B1，脈沖磁場在組織內(nèi)部產(chǎn)生感應(yīng)電流J，與靜磁場B相互作用會產(chǎn)生洛侖茲力，從而激發(fā)組織產(chǎn)生同頻振動，該振動會向組織外傳播包含成像物體電特性的超聲信號，通過重建算法即可重建出組織的電導(dǎo)率分布圖像。在MAT-MI中，聲源振動是由時變洛倫茲力產(chǎn)生，而洛倫茲力與感應(yīng)電流、靜態(tài)磁場相關(guān)，感應(yīng)電流又與時變磁場、成像物體的電導(dǎo)率相關(guān)，因而，通過聲源重建，并建立與成像物體電導(dǎo)率之間的關(guān)系，可以得到電阻抗斷層圖像。MAT-MI聲源的對比度和空間分辨率決定了電阻抗圖像的質(zhì)量。

根據(jù)上述成像原理，可知MAT-MI是通過采集邊界超聲信號重建組織內(nèi)部的電導(dǎo)率分布，其系統(tǒng)模型見圖1，包含從電導(dǎo)率分布到產(chǎn)生超聲信號過程中的各種物理耦合效應(yīng)，超聲換能器特性及磁感應(yīng)磁聲實驗系統(tǒng)的設(shè)置等參數(shù)。對于固定的實驗系統(tǒng)，該系統(tǒng)模型方法可以有效的抑制噪聲，提高成像分辨率。因此，MAT-MI的研究就轉(zhuǎn)化為對系統(tǒng)模型特征分析及其逆系統(tǒng)求解的研究。

圖1 磁感應(yīng)磁聲成像的系統(tǒng)模型Fig 1 The system model of MAT-MI

3 方法

3.1 系統(tǒng)矩陣及其偽逆矩陣

在MAT-MI中，磁場、感應(yīng)渦電流及聲壓是時間和空間的函數(shù)，根據(jù)生物組織中電場磁場聲場的機電耦合機制[15]，聲壓分布由以下波動方程表示：

其中cs是聲在組織中傳播的速度，p（r，t）是聲壓場的時空分布，J（r，t）是感應(yīng)渦流密度，r為無界空間中的任一點，?·[J（r，t）×B0]是聲振源。根據(jù)矢量分解公式，該聲振源可以表示為：

由聲場在自由空間中的格林函數(shù)公式，求解聲壓波動方程可得：

對于電導(dǎo)率均勻變化的組織，上述公式說明電導(dǎo)率與電導(dǎo)率的梯度經(jīng)過一個系統(tǒng)矩陣可以得到相應(yīng)的超聲信號，式（3）可以抽象為一個積分方程的形式，如下：

對于n×n的電導(dǎo)率成像區(qū)域，檢測換能器個數(shù)為m，每個換能器的采樣點數(shù)為k。磁感應(yīng)磁聲成像的系統(tǒng)矩陣模型建立方法見圖2，抽象為矩陣形式見式（5），其中A為磁感應(yīng)磁聲成像系統(tǒng)矩陣，其矩陣大小為mk×n2，x反映了待成像區(qū)域的電導(dǎo)率分布，b為檢測超聲信號。

圖2 磁感應(yīng)磁聲成像系統(tǒng)矩陣模型的建立Fig 2 The system matrix of MAT-MI

對于該不適定矩陣方程有很多正則化方法求解，本研究采用截斷奇異值分解方法（TSVD）[16]進行求解，即將矩陣A分解成特征值與相應(yīng)的特征向量，具體如下：

其中，U和V為正交矩陣，∑為對角矩陣，其對角線值為系統(tǒng)矩陣的特征值。這些特征值反映了電導(dǎo)率分布與超聲聲壓的線性關(guān)系，采用TSVD方法求A的逆矩陣為：

其重建電導(dǎo)率分布可由式（8）給出：

3.2 數(shù)值仿真方法

為驗證上述算法的有效性，建立一個原始電導(dǎo)率分布模型，模型內(nèi)部采用不同形狀及不同電導(dǎo)率進行區(qū)分，外部橢圓結(jié)構(gòu)的電導(dǎo)率設(shè)定為0.1 S/m，內(nèi)部包括一個電導(dǎo)率為0.25 S/m的小橢圓結(jié)構(gòu)，還有電導(dǎo)率為0.35 S/m的一個正方形結(jié)構(gòu)和一個三角形結(jié)構(gòu)，該原始電導(dǎo)率分布見圖3。采用的靜態(tài)磁場強度B0為1 T，為保證脈沖激勵的均勻性，激勵磁場采用赫姆霍茲線圈產(chǎn)生，該線圈注入電流密度設(shè)定為107A。

圖3 原始電導(dǎo)率分布仿體模型Fig 3 The phantom of original conductivity distribution

設(shè)定求解域為一個50 mm×50 mm的求解空間，將該求解空間分為201×201的正方體網(wǎng)格，使用MATLAB進行有限差分計算，設(shè)定采集換能器的個數(shù)為60個，采樣501個點，中心頻率為1 MHz，掃描半徑為45 mm，利用上述條件設(shè)置可以求解系統(tǒng)矩陣A，可知該系統(tǒng)矩陣A為一欠定矩陣。利用該系統(tǒng)矩陣經(jīng)過正問題計算可得到檢測聲壓正弦圖，通過這些檢測到的聲壓，對A進行SVD分解，按照其噪聲級別，便可重建仿體模型的電導(dǎo)率的分布，該數(shù)值仿真方法流程圖見圖4。

圖4 基于系統(tǒng)矩陣重建方法數(shù)值仿真示意圖Fig 4 Image of simulation based on system matrix method

4 結(jié)果

根據(jù)上述仿真過程，采用無噪聲及信噪比為40 dB的含噪聲磁聲信號分別重建，圖5給出了磁感應(yīng)磁聲系統(tǒng)矩陣A的特征值曲線，可以看出，特征值是不斷下降的，且下降的速率很快，最大特征值大約為1.5，最小特征值大約為10-16，因此，在當前實驗系統(tǒng)下，我們只能選取部分特征值進行重建，并根據(jù)噪聲水平可選擇截斷的位置，隨后采用式（8）進行電導(dǎo)率分布重建。

重建結(jié)果見圖6，當信噪比為40 dB時，圖6（a）給出了使用特征值為1 500時的電導(dǎo)率分布重建結(jié)果，圖6（b）給出了使用特征值為3 500時的電導(dǎo)率分布重建結(jié)果，圖6（c）給出了使用特征值為7 500的電導(dǎo)率分布重建結(jié)果，為無噪聲條件下的重建結(jié)果，可以得出，使用系統(tǒng)矩陣的特征值越多，重建結(jié)果越接近真實電導(dǎo)率分布，由于受到噪聲的影響，使用大于7 500個特征值的重建會被噪聲完全覆蓋。因此，在40 dB的噪聲水平下，重建中只能利用7 500個特征值進行電導(dǎo)率重建。為了進一步給出噪聲水平與磁感應(yīng)磁聲重建中使用特征值數(shù)目之間的關(guān)系，算法也對無噪聲理想條件下的重建進行了重建算法驗證，圖6（d）給出了無噪聲條件下使用特征值為12 000的電導(dǎo)率重建結(jié)果，可以看出在欠定條件下，電導(dǎo)率重建結(jié)果中尚存在一些紋波噪聲，但電導(dǎo)率重建結(jié)果與真實電導(dǎo)率分布基本吻合。

圖5 系統(tǒng)矩陣特征值曲線Fig 5 Eigenvalues curve of system matrix

圖6 不同條件下重建電導(dǎo)率分布結(jié)果Fig 6 Distribution of reconstructive conductivity results with different conditions

由上述結(jié)果可以得出，使用系統(tǒng)矩陣方法可以根據(jù)實驗條件最大限度的利用磁聲信號所提供的信息，使重建電導(dǎo)率分布結(jié)果最大限度的接近原始電導(dǎo)率分布。為了體現(xiàn)該算法的優(yōu)勢，我們使用相同仿真實驗條件，采用磁感應(yīng)磁聲的常用重建算法時間反轉(zhuǎn)方法進行重建，。圖7給出了信噪比為40 dB時，采用時間反轉(zhuǎn)方法進行電導(dǎo)率重建的結(jié)果。結(jié)果表明，使用時間反轉(zhuǎn)算法重建的電導(dǎo)率分布僅有邊界，僅相當于使用特征值約為1 500左右的系統(tǒng)矩陣算法的重建結(jié)果，其并未充分的利用磁聲信號的所包含的信息。因此，通過對比可知，本研究算法的重建結(jié)果要優(yōu)于時間反轉(zhuǎn)方法重建結(jié)果。

圖7 時間反轉(zhuǎn)算法電導(dǎo)率重建結(jié)果（SNR =40 d B）Fig 7 Reconstructive conductivity results by using time reversal algorithm（SNR =40 dB）

5 總結(jié)

本研究采用系統(tǒng)矩陣方法，對磁感應(yīng)磁聲電導(dǎo)率重建進行了仿真研究，仿真結(jié)果證明該方法可以在欠定條件及信噪比大約為40 dB的情況下重建出電導(dǎo)率的分布，在更低信噪比條件下，該算法也可實現(xiàn)電導(dǎo)率邊界分布的重建，證明了該算法在MATMI中的適用性。相較其他方法，系統(tǒng)矩陣方法將從電導(dǎo)率分布，到磁感應(yīng)磁聲聲源產(chǎn)生，傳播和接收等過程精確的集成于一個系統(tǒng)矩陣中，通過分析該矩陣的特征值，給出了MAT-MI電導(dǎo)率分布與磁聲信號之間的對應(yīng)關(guān)系，從而最大限度的利用磁聲信號中包含的信息進行電導(dǎo)率重建，提高了成像分辨率和精度，為磁聲成像檢測系統(tǒng)建立和應(yīng)用提供基礎(chǔ)。