基于小波稀疏的磁性納米粒子成像算法研究

張玉錄 柯麗 杜強 趙宇楠 祖婉妮

0 引言

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

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

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

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

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

1 磁性納米粒子成像基礎研究

1.1 磁性納米粒子成像原理

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

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

隨后，對成像區域進行空間編碼[13]，通過成像方法建立接收線圈感應產生的MPI信號與各個編碼點粒子濃度值間的關系，是實現磁性納米粒子成像的關鍵。

1.2 系統矩陣成像方法

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

(1)

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

(2)

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

G×c=u

(3)

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

(4)

展開表示為：

(5)

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

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

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

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

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

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

2 基于小波稀疏的MPI算法

2.1 小波稀疏處理

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

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

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

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

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

2.2 基于小波稀疏的MPI算法

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

(6)

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

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

(7)

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

通過以上的迭代過程，可以實現對待測區域中各編碼點的粒子濃度求解。需要注意的是，粒子濃度值是一個非負的值，這要求對迭代后的粒子濃度數值進行非負運算。在提出的基于小波稀疏MPI算法中，算法的收斂性與系統噪聲水平有關，噪聲水平越大，算法收斂性越差。

3 仿真實驗與結果分析

3.1 仿真實驗

采用圖4所示開放式線性零磁場MPI電磁系統模型進行仿真，獲取MPI信號。

圖4 線性零磁場MPI電磁系統Figure 4 MPI electromagnetic system of linear zero magnetic field

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

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

系統搭建完畢后，對兩個接收線圈中間面積為3.1 cm×3.1 cm大小的成像區域按照1 mm的上下間距進行空間編碼，得到31×31共961個編碼點。通過模型與計算相結合的方法確定各編碼點單獨存在粒子時在接收線圈中感應產生的MPI信號，提取各編碼點處粒子產生MPI信號頻譜的1 000個頻率分量得到大小為961×1 000的系統矩陣。

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

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

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

本文仿真實驗分別采用圖6所示的幾何形粒子分布模型和血管模型進行成像實驗，其中點狀編碼點對應區域的粒子濃度為1 mol/L，“x”狀編碼點對應區域的粒子濃度為0.5 mol/L。

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

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

(10)

在相同的終止條件下，繪制出PSNR隨閾值λ的變化曲線，當PSNR達到最大值時，選擇對應的閾值λ為最終的閾值參數值。算法實現過程中采用基于Haar小波的2級SWT對迭代變量進行處理。

3.2 結果分析

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

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

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

對代表兩種算法所成圖像的質量差異參數值進行顯示，如圖9所示。

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

橫向對比發現，兩種成像模型下，基于小波稀疏的MPI算法所成圖像均表現出更好的成像質量。由于閾值算子統一將閾值以下的SWT各分量置零，結果在每次迭代結束后都對得到的粒子濃度信息進行稀疏處理，不含粒子編碼點對應的粒子濃度值被有效降低，無粒子區域與粒子均勻分布區域間的濃度梯度相應地增加，實現了對MPI所成圖像中粒子分布邊界非平穩特征的良好成像，所以無論采用哪種閾值算子，基于小波稀疏的MPI算法所成圖像都表現出更好的成像質量。

縱向對比發現，兩種成像模型下，當系統信噪比為30 dB、20 dB、10 dB時，基于小波稀疏的MPI算法與經典ART算法相比所成圖像的PSNR參數平均提升了67.83%、18.66%、8.05%。在基于小波稀疏的MPI算法中，小波稀疏處理對圖像質量的提升幅度取決于閾值算子對SWT各分量的稀疏程度，隨著噪聲水平的增加，部分不含粒子編碼點處對應的分量值接近于存在粒子編碼點處對應的分量，這時小波稀疏處理也無法將這些過大的干擾信號有效除去，所以圖像質量提升的幅度會隨之下降。

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

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

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

4 討論與總結

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

在提出的基于小波稀疏的MPI算法中，將小波稀疏處理融入到經典ART迭代過程中，可在每次迭代后耗費極少時間強化所成圖像中的粒子分布邊緣非平穩特征，當算法達到收斂時，計算得到的各編碼點粒子濃度值就更接近真實粒子濃度值。算法收斂速度與系統中的噪聲水平有關，系統噪聲水平越低，算法達到收斂所需的時間就越少，所成圖像的質量也越高。綜上，在較低的噪聲水平下，基于小波稀疏的MPI算法可以在短時間內實現對粒子分布狀況的高質量成像。

在今后的工作中，將著重研究基于小波稀疏的MPI算法與MPI成像系統的兼容性，從而推動本文方法在MPI實時精準成像方面的實際應用。