張玉錄 柯麗 杜強(qiáng) 趙宇楠 祖婉妮

0 引言

磁性納米粒子成像(magnetic particle imaging，MPI)是一種可實(shí)現(xiàn)快速、精準(zhǔn)成像的新興臨床成像技術(shù)。該技術(shù)通過檢測磁性納米粒子在變化磁場中的非線性磁化特性來重構(gòu)粒子在待測區(qū)域中的濃度分布狀況。在該過程中，變化的磁場產(chǎn)生移動的零場區(qū)域，只有零場區(qū)域中粒子磁矩會發(fā)生變化，進(jìn)而在接收裝置中產(chǎn)生MPI信號[1]，這不僅有利于MPI的高分辨率成像，也可實(shí)現(xiàn)MPI的快速動態(tài)成像，所以MPI在臨床上的快速診斷和實(shí)時(shí)監(jiān)控中具有廣泛的應(yīng)用前景，已被證實(shí)可以應(yīng)用于細(xì)胞成像[2]、血管成像[3]、腫瘤成像[4]、干細(xì)胞追蹤[5]、熱療[6]和藥物靶向[7]等多個(gè)領(lǐng)域。

MPI成像方法是實(shí)現(xiàn)快速、精準(zhǔn)成像的關(guān)鍵。目前MPI常用的成像方法主要有系統(tǒng)矩陣(system matrix，SM)方法和X-space方法[1]。由于X-space方法對零場區(qū)域的磁場和粒子在磁場中的弛豫時(shí)間提出更高要求[8]，SM方法成為MPI成像方案的首要選擇。在SM方法中，由SM和測量向量構(gòu)成矩陣方程的求解算法決定了成像速度和質(zhì)量。

在SM成像方案的早期研究過程中，由于噪聲的存在，矩陣方程的求解問題被轉(zhuǎn)化為最小二乘問題，采用正則化和矩陣求逆技術(shù)求解該最小二乘問題是求解粒子濃度信息的常用方法[1]，但是大型MPI矩陣的逆計(jì)算復(fù)雜。隨著MPI技術(shù)的發(fā)展，Lübeck大學(xué)的Panagiotopoulos等[1]針對大型SM在矩陣求逆中面臨的大計(jì)算量問題，提出一種采用收斂迅速的代數(shù)重建算法(algebraic reconstruction technique，ART)求解吉洪諾夫正則化最小二乘模型成像的方法，成像結(jié)果表明ART算法大大縮短了重建時(shí)間,但是所成圖像邊緣較為模糊。Heidelberg大學(xué)的Storath等[9]針對吉洪諾夫正則化最小二乘模型所成圖像邊緣模糊的問題，提出用交替方向乘子算法(alternating direction method of multipliers，ADMM)求解非負(fù)融合Lasso模型，成像結(jié)果表明該算法可以實(shí)現(xiàn)粒子分布的精準(zhǔn)成像，但是成像耗時(shí)是ART算法的7倍。近年來，研究人員相繼提出隨機(jī)奇異值分解加速的ART算法[10]、逐級背景分割的圖像重建算法[11]、編碼校準(zhǔn)場景框架重構(gòu)SM方法[12]等多種成像方法，雖然在求解速度或者成像質(zhì)量上有一定的優(yōu)勢，但難以滿足MPI的快速精準(zhǔn)成像要求。

在ART快速成像算法基礎(chǔ)上通過簡單的圖像處理改善圖像邊緣模糊問題，是實(shí)現(xiàn)MPI快速精準(zhǔn)成像的一種有效方法。小波變換具有良好的表征圖像局部特征的能力，且運(yùn)算耗時(shí)較少。就邊緣模糊的圖像而言，采用簡單數(shù)學(xué)工具對小波變換提取的圖像邊緣特征進(jìn)行處理，可以快速去除圖像中的干擾信息，提升邊緣區(qū)域的圖像質(zhì)量。

本文為了實(shí)現(xiàn)MPI中粒子濃度空間分布的快速精準(zhǔn)成像，針對系統(tǒng)矩陣成像方法構(gòu)建的矩陣方程求解問題，提出一種基于小波變換和閾值算子的迭代MPI算法。在ART算法每次迭代結(jié)束后，首先采用小波變換提取圖像中粒子分布邊界的非平穩(wěn)特征，然后通過閾值算子稀疏運(yùn)算去除圖像中的干擾信號，增加少量計(jì)算時(shí)間的同時(shí)通過強(qiáng)化圖像非平穩(wěn)特征提高了粒子分布邊緣圖像的清晰度，結(jié)合采用的線性零磁場MPI電磁系統(tǒng)可在短時(shí)間內(nèi)實(shí)現(xiàn)對粒子分布狀況的高質(zhì)量成像。

1 磁性納米粒子成像基礎(chǔ)研究

1.1 磁性納米粒子成像原理

MPI通過檢測磁性納米粒子在變化磁場中的非線性磁化特性[圖1(a)]來重構(gòu)粒子在待測區(qū)域中的濃度分布狀況。在MPI中，選擇場產(chǎn)生零場區(qū)域，聚焦場控制零場區(qū)域的自由移動，當(dāng)對成像區(qū)域施加一個(gè)時(shí)變驅(qū)動磁場[圖1(c)]時(shí)，只有零場區(qū)域中的粒子磁化強(qiáng)度改變，從而產(chǎn)生磁化響應(yīng)信號[圖1(b)]。粒子磁化響應(yīng)信號會在接收線圈中感應(yīng)產(chǎn)生電壓信號[圖1(d)]，該電壓信號即為MPI信號。

圖1 磁性納米粒子對外部磁場的響應(yīng)Figure 1 The response of magnetic particles to an external magnetic field

隨后，對成像區(qū)域進(jìn)行空間編碼[13]，通過成像方法建立接收線圈感應(yīng)產(chǎn)生的MPI信號與各個(gè)編碼點(diǎn)粒子濃度值間的關(guān)系，是實(shí)現(xiàn)磁性納米粒子成像的關(guān)鍵。

1.2 系統(tǒng)矩陣成像方法

為了實(shí)現(xiàn)對粒子分布狀況的精準(zhǔn)成像，采用SM成像方法[14]。在SM成像方法中，接收線圈感應(yīng)產(chǎn)生的MPI信號u(t)與各編碼點(diǎn)rn(n=1,…,N)處的粒子濃度信息在時(shí)域的關(guān)系可用式(1)表示：

(1)

式中：M(rn,t)為位于位置rn處的粒子產(chǎn)生的磁化響應(yīng)信號；s(rn)為線圈靈敏度；Δv代表rn處粒子對應(yīng)的體積，則：

(2)

用矩陣-向量形式表示線性系統(tǒng)方程組(1)如式(3)所示：

G×c=u

(3)

式中：u是接收線圈感應(yīng)產(chǎn)生的MPI信號；c是待測的磁性納米粒子的濃度；G是SM。對等式兩側(cè)的數(shù)據(jù)進(jìn)行傅里葉變換，就可以得到方程組在頻域的矩陣-向量形式：

(4)

展開表示為：

(5)

式中：m代表頻率分量的個(gè)數(shù)；n代表成像區(qū)域中編碼點(diǎn)的個(gè)數(shù)。通過求解式(5)所示的矩陣方程組即可獲得各個(gè)編碼點(diǎn)的粒子濃度信息，進(jìn)而成像。

1.3 基于ART迭代的矩陣方程求解算法

在實(shí)際的MPI中，通常檢測不同斷層血管中粒子的濃度分布狀況進(jìn)行成像。對于單一斷層而言，首先，總是存在粒子濃度為零的組織部分；其次，血管邊界使得粒子在血管中的分布相對均勻，這樣斷層上的粒子分布圖像就具有分段恒定區(qū)域。在對圖像恒定區(qū)域精確成像前提下實(shí)現(xiàn)粒子分布邊緣區(qū)域的良好成像，是實(shí)現(xiàn)高質(zhì)量成像的關(guān)鍵，這也對矩陣方程(5)的求解算法提出了更高的要求。

在理想的MPI電磁系統(tǒng)中，基于定點(diǎn)迭代的經(jīng)典ART算法是求解矩陣方程(5)的快速有效方法。而實(shí)際的MPI系統(tǒng)中，通常結(jié)合正則化運(yùn)算對矩陣方程(5)進(jìn)行擴(kuò)展來減小外界噪聲對系統(tǒng)的影響，從而實(shí)現(xiàn)粒子濃度分布的良好成像[15]。在含噪聲的一維MPI電磁系統(tǒng)中，對圖2(a)所示的粒子分布模型進(jìn)行成像，采用經(jīng)典ART算法求解擴(kuò)展后的矩陣方程所成圖像如圖2(b)所示。

圖2 一維MPI粒子濃度分布對比圖Figure 2 Contrast diagram of one-dimensional MPI particle concentration distribution

盡管已采用正則化運(yùn)算降低噪聲的影響，但經(jīng)典ART算法求得的粒子濃度值相較真實(shí)粒子濃度值仍然存在一定誤差，特別是粒子分布邊界附近的粒子濃度差異，弱化了粒子分布邊緣處的非平穩(wěn)特征，圖像相應(yīng)地變得模糊。

2 基于小波稀疏的MPI算法

2.1 小波稀疏處理

通過簡單運(yùn)算增強(qiáng)MPI中被弱化的非平穩(wěn)特征是實(shí)現(xiàn)MPI快速、精準(zhǔn)成像的關(guān)鍵。為確保MPI信號的完整性，采用平穩(wěn)小波變換(stationary wavelet transform，SWT)對MPI信號進(jìn)行處理，MPI信號中的非平穩(wěn)特征在SWT后的小波域表現(xiàn)為數(shù)值較大的小波系數(shù)，而不含粒子區(qū)域表現(xiàn)為數(shù)值較小的小波系數(shù)，對變換后的小波系數(shù)做稀疏處理可以強(qiáng)化信號的非平穩(wěn)特征，提升成像質(zhì)量。

閾值算子可實(shí)現(xiàn)對信號的稀疏處理。軟閾值(soft thresholding)、硬閾值(hard thresholding)、NNG(nonnegative garotte thresholding)閾值是3種最常見的閾值算子，3種閾值算子處理后的數(shù)據(jù)除了在閾值處及高于閾值部分產(chǎn)生一定差異性外，均可實(shí)現(xiàn)對信號的稀疏處理。SWT和閾值算子結(jié)合的小波稀疏處理可以降低MPI信號中干擾信號的強(qiáng)度，從而增強(qiáng)粒子分布圖像中被弱化的非平穩(wěn)特征。

以圖2(b)中采用ART算法求得的一維MPI信號為例，采用不同閾值算子下的小波稀疏處理對該信號進(jìn)行處理，得到不同狀況下的重構(gòu)信號如圖3所示。

圖 3 小波稀疏處理前后效果對比Figure 3 Comparison of effects before and after wavelet sparse processing

3種閾值算子下的小波稀疏處理均可實(shí)現(xiàn)對MPI信號的稀疏處理。不同的是，基于軟閾值的小波稀疏處理重構(gòu)信號中粒子濃度數(shù)值明顯小于真實(shí)值，信號失真明顯；基于NNG閾值的小波稀疏處理重構(gòu)信號在粒子濃度數(shù)值稍微減小的基礎(chǔ)上，實(shí)現(xiàn)對非平穩(wěn)特征的良好重構(gòu)；基于硬閾值的小波稀疏處理雖然可良好實(shí)現(xiàn)對粒子分布狀況成像，但是在粒子分布邊界處會出現(xiàn)振蕩，從而影響重構(gòu)信號的精度，據(jù)此本文采用NNG閾值算子進(jìn)行稀疏處理。綜上，基于NNG閾值的小波稀疏處理可以降低ART解中干擾信息的強(qiáng)度，進(jìn)而實(shí)現(xiàn)對存在非平穩(wěn)特征的MPI粒子信號的優(yōu)化。

2.2 基于小波稀疏的MPI算法

將小波稀疏處理應(yīng)用于經(jīng)典ART算法的迭代過程中，可在每次迭代后強(qiáng)化所成圖像中的粒子分布邊界非平穩(wěn)特征，從而有效解決粒子分布圖像中的邊界模糊問題。在這種情況下：令φ：Rn→RJ×n表示分解級數(shù)為J的平穩(wěn)小波變換，MPI中粒子濃度的重構(gòu)問題[式(5)]可以被轉(zhuǎn)化為式(6)的最小化問題：

(6)

在該最小化問題中，函數(shù)f代表通過閾值算子促進(jìn)平穩(wěn)小波變換后小波系數(shù)φc的稀疏處理過程。

綜上，對于式(6)的最小化問題，可應(yīng)用基于小波稀疏的MPI算法求解，算法實(shí)現(xiàn)過程為：首先用正則化運(yùn)算對線性系統(tǒng)[式(5)]進(jìn)行擴(kuò)展，然后采用經(jīng)典ART算法對擴(kuò)展后的線性系統(tǒng)進(jìn)行求解，之后對求解得到的粒子濃度c執(zhí)行SWT，從而得到粒子濃度c的近似分量和細(xì)節(jié)分量，然后用閾值算子對小波變換后的各分量進(jìn)行稀疏處理，再將閾值處理后的各分量通過小波反變換得到新的粒子濃度c，將所得到的粒子濃度重新帶入到算法中進(jìn)行下次迭代過程，重復(fù)迭代直到所得到的粒子濃度滿足某種迭代終止條件算法終止。單次迭代過程可用式(7)表示：

(7)

式中：gi表示系統(tǒng)矩陣G的行向量；α為正則化參數(shù)?；谛〔ㄏ∈璧腗PI算法偽代碼如下所示：

通過以上的迭代過程，可以實(shí)現(xiàn)對待測區(qū)域中各編碼點(diǎn)的粒子濃度求解。需要注意的是，粒子濃度值是一個(gè)非負(fù)的值，這要求對迭代后的粒子濃度數(shù)值進(jìn)行非負(fù)運(yùn)算。在提出的基于小波稀疏MPI算法中，算法的收斂性與系統(tǒng)噪聲水平有關(guān)，噪聲水平越大，算法收斂性越差。

3 仿真實(shí)驗(yàn)與結(jié)果分析

3.1 仿真實(shí)驗(yàn)

采用圖4所示開放式線性零磁場MPI電磁系統(tǒng)模型進(jìn)行仿真，獲取MPI信號。

圖4 線性零磁場MPI電磁系統(tǒng)Figure 4 MPI electromagnetic system of linear zero magnetic field

在線性零磁場MPI電磁系統(tǒng)中：采用上下各2個(gè)矩形線圈組成的雙平面梯度線圈結(jié)構(gòu)組成一組選擇場線圈，兩組選擇場線圈共同作用就可實(shí)現(xiàn)零場線在待成像平面的自由旋轉(zhuǎn)；采用同樣由矩形線圈組成的雙平面梯度線圈組成聚焦場線圈，主要作用是控制選擇場產(chǎn)生的零場線在待測平面平行移動。對選擇場和聚焦場同時(shí)施加不同的電流激勵，就可使零場線按照不同的掃描軌跡對待測區(qū)域進(jìn)行掃描，本文采用成像效果良好的徑向軌跡(圖5)對待測區(qū)域掃描。亥姆霍茲結(jié)構(gòu)的驅(qū)動場線圈在25 kHz激勵頻率的電流作用下可在成像區(qū)域z方向上產(chǎn)生磁通密度振幅為2.6 T的均勻變化磁場，該變化磁場可使所成圖像的分辨率達(dá)到1 mm以下[16]。

圖5 零場線徑向掃描軌跡Figure 5 Radial scan trace of free field line

系統(tǒng)搭建完畢后，對兩個(gè)接收線圈中間面積為3.1 cm×3.1 cm大小的成像區(qū)域按照1 mm的上下間距進(jìn)行空間編碼，得到31×31共961個(gè)編碼點(diǎn)。通過模型與計(jì)算相結(jié)合的方法確定各編碼點(diǎn)單獨(dú)存在粒子時(shí)在接收線圈中感應(yīng)產(chǎn)生的MPI信號，提取各編碼點(diǎn)處粒子產(chǎn)生MPI信號頻譜的1 000個(gè)頻率分量得到大小為961×1 000的系統(tǒng)矩陣。

實(shí)際MPI成像過程中，受外界噪聲的影響，接收線圈感應(yīng)到的MPI信號與真實(shí)MPI信號間存在隨機(jī)性的偏差，該偏差可用高斯白噪聲表示。添加高斯白噪聲的MPI信號與SM組成矩陣方程，求解矩陣方程即可得到各編碼點(diǎn)的粒子濃度信息，進(jìn)而成像。

所成圖像中的粒子濃度值與真實(shí)粒子濃度值間的差異可用峰值信噪比參數(shù)(peak signal-to-noise ratio，PSNR)、圖形相似性參數(shù)(structural similarity，SSIM)來衡量，所成圖像中的粒子濃度值與真實(shí)粒子濃度值間的差異性越小，峰值信噪比參數(shù)值越大、圖形相似性參數(shù)值也越大。兩個(gè)參數(shù)的計(jì)算方法由式(8)和式(9)決定：

式中：x代表理想的粒子分布圖像；y代表通過算法所成圖像；xmax是x中最大像素值；μx、μy分別表示x、y中各像素值的均值；δx、δy分別表示x、y中各像素值的方差；δxy表示x和y像素值間的協(xié)方差；c1和c2是兩個(gè)常數(shù)，主要的作用是避免式子中的分母為零。

本文仿真實(shí)驗(yàn)分別采用圖6所示的幾何形粒子分布模型和血管模型進(jìn)行成像實(shí)驗(yàn)，其中點(diǎn)狀編碼點(diǎn)對應(yīng)區(qū)域的粒子濃度為1 mol/L，“x”狀編碼點(diǎn)對應(yīng)區(qū)域的粒子濃度為0.5 mol/L。

圖6 粒子分布模型示意圖Figure 6 Schematic diagram of particle distribution model

實(shí)驗(yàn)過程中設(shè)置算法的最大迭代次數(shù)為5 000，當(dāng)兩個(gè)連續(xù)迭代之間的粒子濃度間的差異參數(shù)er[式(10)]小于10-5時(shí)，迭代終止。

(10)

在相同的終止條件下，繪制出PSNR隨閾值λ的變化曲線，當(dāng)PSNR達(dá)到最大值時(shí)，選擇對應(yīng)的閾值λ為最終的閾值參數(shù)值。算法實(shí)現(xiàn)過程中采用基于Haar小波的2級SWT對迭代變量進(jìn)行處理。

3.2 結(jié)果分析

在信噪比分別為30 dB、20 dB、10 dB的噪聲下，采用基于小波稀疏的MPI算法和經(jīng)典ART算法分別對圖6所示幾何形粒子分布模型和血管模型進(jìn)行成像，得到的成像結(jié)果如圖7和圖8所示。

圖7 基于小波稀疏的MPI算法和經(jīng)典ART算法成像結(jié)果(幾何模型)對比Figure 7 Comparison of imaging results (geometric model) between wavelet sparse MPI algorithm and classic ART algorithm

圖8 基于小波稀疏的MPI算法和經(jīng)典ART算法成像結(jié)果(血管模型)對比Figure 8 Comparison of imaging results (blood vessel model) between wavelet sparse MPI algorithm and classic ART algorithm

對代表兩種算法所成圖像的質(zhì)量差異參數(shù)值進(jìn)行顯示，如圖9所示。

圖9 基于小波稀疏的MPI算法和經(jīng)典ART算法成像質(zhì)量對比Figure 9 Comparison of imaging quality between MPI algorithm based on wavelet sparseness and classic ART algorithm

橫向?qū)Ρ劝l(fā)現(xiàn)，兩種成像模型下，基于小波稀疏的MPI算法所成圖像均表現(xiàn)出更好的成像質(zhì)量。由于閾值算子統(tǒng)一將閾值以下的SWT各分量置零，結(jié)果在每次迭代結(jié)束后都對得到的粒子濃度信息進(jìn)行稀疏處理，不含粒子編碼點(diǎn)對應(yīng)的粒子濃度值被有效降低，無粒子區(qū)域與粒子均勻分布區(qū)域間的濃度梯度相應(yīng)地增加，實(shí)現(xiàn)了對MPI所成圖像中粒子分布邊界非平穩(wěn)特征的良好成像，所以無論采用哪種閾值算子，基于小波稀疏的MPI算法所成圖像都表現(xiàn)出更好的成像質(zhì)量。

縱向?qū)Ρ劝l(fā)現(xiàn)，兩種成像模型下，當(dāng)系統(tǒng)信噪比為30 dB、20 dB、10 dB時(shí)，基于小波稀疏的MPI算法與經(jīng)典ART算法相比所成圖像的PSNR參數(shù)平均提升了67.83%、18.66%、8.05%。在基于小波稀疏的MPI算法中，小波稀疏處理對圖像質(zhì)量的提升幅度取決于閾值算子對SWT各分量的稀疏程度，隨著噪聲水平的增加，部分不含粒子編碼點(diǎn)處對應(yīng)的分量值接近于存在粒子編碼點(diǎn)處對應(yīng)的分量，這時(shí)小波稀疏處理也無法將這些過大的干擾信號有效除去，所以圖像質(zhì)量提升的幅度會隨之下降。

分別以幾何形粒子分布模型和血管模型為研究對象，在不同噪聲水平下對基于小波稀疏的MPI算法和經(jīng)典ART算法收斂過程進(jìn)行分析，得到的圖像如圖10所示。

圖10 基于小波稀疏的MPI算法和經(jīng)典ART算法收斂過程Figure 10 The convergence process diagram of the MPI algorithm based on wavelet sparseness and classic ART algorithm

對比發(fā)現(xiàn)，基于小波稀疏的MPI算法收斂速度相較經(jīng)典ART算法較慢，成像時(shí)間會有所增加。在基于小波稀疏的MPI算法迭代過程中，單次迭代所得的粒子濃度信號是由上次迭代所得粒子濃度信號經(jīng)過ART迭代與小波稀疏處理后產(chǎn)生，相鄰的兩次迭代過程所得粒子濃度信號之間的差異參數(shù)相較經(jīng)典ART算法更大，所以當(dāng)設(shè)置的終止條件一致時(shí)，基于小波稀疏的MPI算法達(dá)到穩(wěn)定就需要更多次的迭代，但是小波稀疏處理對粒子信號的優(yōu)化作用是有限的，所以兩種算法迭代次數(shù)的差異不是特別明顯。這樣，基于小波稀疏的MPI算法就可在稍微增加成像時(shí)間的前提下實(shí)現(xiàn)更高質(zhì)量的成像。

4 討論與總結(jié)

本文為了實(shí)現(xiàn)MPI中粒子濃度空間分布的快速精準(zhǔn)成像，提出在耗時(shí)較短的經(jīng)典ART算法每次迭代結(jié)束后對當(dāng)前解進(jìn)行小波稀疏處理的MPI算法。在小波稀疏處理過程中，SWT在保證信號完整性的前提下可實(shí)現(xiàn)對粒子分布邊緣非平穩(wěn)特征的良好提取，閾值算子對SWT后各分量的稀疏處理在確保含粒子區(qū)域粒子濃度值變化不大的前提下，有效降低了不含粒子區(qū)域的干擾信息強(qiáng)度，無粒子區(qū)域與粒子均勻分布區(qū)域間的濃度梯度相應(yīng)增加，實(shí)現(xiàn)了對MPI所成圖像中粒子分布邊界非平穩(wěn)特征的良好成像。

在提出的基于小波稀疏的MPI算法中，將小波稀疏處理融入到經(jīng)典ART迭代過程中，可在每次迭代后耗費(fèi)極少時(shí)間強(qiáng)化所成圖像中的粒子分布邊緣非平穩(wěn)特征，當(dāng)算法達(dá)到收斂時(shí)，計(jì)算得到的各編碼點(diǎn)粒子濃度值就更接近真實(shí)粒子濃度值。算法收斂速度與系統(tǒng)中的噪聲水平有關(guān)，系統(tǒng)噪聲水平越低，算法達(dá)到收斂所需的時(shí)間就越少，所成圖像的質(zhì)量也越高。綜上，在較低的噪聲水平下，基于小波稀疏的MPI算法可以在短時(shí)間內(nèi)實(shí)現(xiàn)對粒子分布狀況的高質(zhì)量成像。

在今后的工作中，將著重研究基于小波稀疏的MPI算法與MPI成像系統(tǒng)的兼容性，從而推動本文方法在MPI實(shí)時(shí)精準(zhǔn)成像方面的實(shí)際應(yīng)用。