具有噪聲魯棒性的三維磁性納米粒子成像快速重構(gòu)方法

史玥婷，任仕偉，王曉華

（1. 北京理工大學(xué) 集成電路與電子學(xué)院， 北京 100081；2. 北京理工大學(xué) 長(zhǎng)三角研究院(嘉興)， 浙江， 嘉興314019；3. 北京理工大學(xué) 重慶微電子中心， 重慶 401332）

磁性納米粒子成像（magnetic particle imaging, MPI）是一種全新的示蹤劑成像技術(shù)[1]，在血管造影[2]、腫瘤成像[3]、體細(xì)胞追蹤[4]和磁熱療[5]等領(lǐng)域具有重要的應(yīng)用前景. 該技術(shù)利用超順磁性納米粒子對(duì)變化磁場(chǎng)會(huì)產(chǎn)生非線性響應(yīng)的原理，完成對(duì)磁粒子濃度空間分布的成像. 如何快速且高質(zhì)量地從接收到磁性粒子電信號(hào)重構(gòu)出的MPI 圖像，一直是國(guó)內(nèi)外研究的熱點(diǎn)與難點(diǎn).

MPI 重構(gòu)的主要方法可以分為系統(tǒng)矩陣法和XSpace 重構(gòu)法. 系統(tǒng)矩陣重構(gòu)算法[1]是將MPI 接收線圈檢測(cè)到感應(yīng)信號(hào)與磁性納米粒子示蹤劑（super paramagnetic iron Oxide，SPIO）空間分布之間的映射關(guān)系描述為一個(gè)系統(tǒng)矩陣，然后再利用校驗(yàn)好的系統(tǒng)矩陣對(duì)待檢測(cè)的圖像進(jìn)行重構(gòu). 該方法雖然可以準(zhǔn)確地描述系統(tǒng)特性，但校準(zhǔn)數(shù)據(jù)采集時(shí)間較長(zhǎng)、計(jì)算復(fù)雜度高和對(duì)內(nèi)存要求高. X-Space 重構(gòu)算法的主要思想是使編碼直接在接收到信號(hào)的空間中進(jìn)行[6]，直接快速的在信號(hào)接收空間完成空間編碼并重建出MPI 圖像，可適用于一維至多維度MPI 成像和重構(gòu).在進(jìn)行三維MPI 成像時(shí)，通常需要獲取符合奈奎斯特采樣定律要求的投影，再通過(guò)濾波反投影等方法進(jìn)行重構(gòu).

根據(jù)MPI 基本原理，可以通過(guò)提高硬件性能的方法來(lái)提升成像速度，比如提高梯度磁場(chǎng)強(qiáng)度等等，但由于硬件的成本過(guò)高和磁場(chǎng)強(qiáng)度大小的限制，導(dǎo)致亟需在成像和重構(gòu)算法上進(jìn)行研究改進(jìn).DONOHO[7]提出的壓縮感知理論表明信號(hào)具有稀疏性，當(dāng)信號(hào)在某個(gè)變換域下表示為稀疏時(shí)，通過(guò)求解一個(gè)非線性最優(yōu)化問(wèn)題，可以以遠(yuǎn)低于奈奎斯特采樣定理要求的采樣數(shù)據(jù)來(lái)重建信號(hào). 隨后，國(guó)內(nèi)外研究人員相繼提出可用于CT、MRI 醫(yī)學(xué)圖像的壓縮感知重構(gòu)算法[8?9]. 就稀疏變換和重構(gòu)算法而言，現(xiàn)有算法在充分利用圖像的幾何正則性、選取最優(yōu)和最稀疏表達(dá)、重構(gòu)圖像計(jì)算復(fù)雜等方面還存在進(jìn)一步優(yōu)化空間. 此外，這些方法沒(méi)有考慮到MPI 的成像原理，無(wú)法直接用于MPI 圖像的稀疏采樣與重構(gòu). 2013年，KONKLE 等[10]提出了基于X-Space 理論的零磁場(chǎng)點(diǎn)磁性納米粒子成像（field free point MPI，F(xiàn)FP-MPI）直流信號(hào)恢復(fù)凸優(yōu)化的重構(gòu)算法，引入了MPI 圖像的連續(xù)性和非負(fù)性兩個(gè)先驗(yàn)條件，為進(jìn)一步研究零磁場(chǎng)線磁性納米粒子成像（field free line MPI ，F(xiàn)FLMPI）凸優(yōu)化重構(gòu)算法提供了新的方向.

本文通過(guò)設(shè)計(jì)基于X-Space 理論的磁性納米粒子成像稀疏采樣與重構(gòu)方法，提出具有噪聲魯棒性的MPI 成像模型正則優(yōu)化函數(shù)，實(shí)現(xiàn)3D TV 稀疏算子及低運(yùn)算復(fù)雜度重構(gòu)的軟件程序框架，充分利用MPI 先驗(yàn)信息解決傳統(tǒng)方法對(duì)成像投影數(shù)量完備性的要求，大幅提升MPI 三維成像和重構(gòu)速度，為今后MPI 三維成像在各種生命科學(xué)及應(yīng)用研究中發(fā)揮更好作用提供了重要基礎(chǔ).

1 磁性納米粒子成像模型

1.1 點(diǎn)擴(kuò)散函數(shù)及基礎(chǔ)成像模型

磁性納米粒子成像系統(tǒng)與常見(jiàn)的醫(yī)學(xué)成像方法同樣具有線性移不變（linear shift invariant, LSI）特性，MPI 圖像可表示為納米粒子濃度與MPI 成像點(diǎn)擴(kuò)散函數(shù)（point spread function, PSF）的卷積[6,11]. MPI 的點(diǎn)擴(kuò)散函數(shù)為L(zhǎng)angevin 響應(yīng)的導(dǎo)數(shù)形式，也可以稱為成像模糊方程，MPI 的分辨率由PSF 的寬度決定.

在MPI 系統(tǒng)中，零磁場(chǎng)區(qū)域（field free region, FFR）是由兩個(gè)相向的磁場(chǎng)創(chuàng)建產(chǎn)生的，由一個(gè)偏置磁場(chǎng)驅(qū)動(dòng)FFR 在視場(chǎng)內(nèi)按設(shè)定的軌跡進(jìn)行掃描，還有一個(gè)高頻正弦波激勵(lì)磁場(chǎng)用于激發(fā)SPIO 產(chǎn)生信號(hào). 當(dāng)FFR 經(jīng)過(guò)SPIO 位置時(shí)，檢測(cè)線圈即可以采集到SPIO 磁矩方向翻轉(zhuǎn)的信號(hào)，由于SPIO 信號(hào)在高頻諧波部分，因此需要濾除基頻信號(hào)，然后使用濾波后的磁粒子信號(hào)進(jìn)行圖像重建.

施加一個(gè)強(qiáng)磁梯度為 ?μ0G的磁場(chǎng)，那么待成像樣品中的SPIO 處于的磁場(chǎng)大小為

為了刺激粒子翻轉(zhuǎn)，需要對(duì)樣品再施加一個(gè)激勵(lì)磁場(chǎng)Hs，該激勵(lì)磁場(chǎng)一般為正弦激勵(lì)場(chǎng). 通過(guò)求解H(x)+Hs(t)=0，得到瞬時(shí)場(chǎng)FFP 信號(hào)為

如前所述，MPI 信號(hào)是通過(guò)接收線圈的磁通 ?變化的結(jié)果. 這種變化是由于納米粒子的總磁化強(qiáng)度M的變化引起的. 在一維情況下，可以得出[6]

其中，B1只是一個(gè)常數(shù)，它模擬了接收線圈的靈敏度.朗之萬(wàn)理論用來(lái)描述SPIO 的磁化強(qiáng)度是其所受外加磁場(chǎng)的函數(shù)，磁化強(qiáng)度的大小可以看作SPIO 粒子中與外加磁場(chǎng)對(duì)齊的原子偶極子的數(shù)量. 如果所有的原子偶極子都對(duì)齊，則稱之為磁飽和. 朗之萬(wàn)理論認(rèn)為磁場(chǎng)H中單點(diǎn)x的磁化強(qiáng)度為

由式(2)和式(3)可進(jìn)一步推導(dǎo)出

式中： γ ?B1mG/Hsat；m[Am2]為磁性納米粒子的磁矩；Hsat[A/m]為使納米顆粒示蹤物半飽和所需的磁場(chǎng)振幅；s(t) [m/s]為瞬時(shí)FFP 速度； L為表征SPIO 磁化特性的Langevin 函數(shù).

利用X-Space 重構(gòu)方法可以將時(shí)域信號(hào)方程直接解析地轉(zhuǎn)換成本地MPI 圖像. 這只需要將接收到的信號(hào)按瞬時(shí)速度歸一化，然后網(wǎng)格化到FFP 的瞬時(shí)位置.

1.2 零磁場(chǎng)線投影成像及重構(gòu)

在本文實(shí)驗(yàn)平臺(tái)上，F(xiàn)FL 是沿y軸方向驅(qū)動(dòng)場(chǎng)激勵(lì)矢量沿z軸方向，移位矢量沿x軸方向，掃描系統(tǒng)[x y z]T是 非旋轉(zhuǎn)框架，樣品成像平 臺(tái) [x′y′z′]T是旋轉(zhuǎn)平臺(tái). 為了獲得投影重建所需的多個(gè)投影圖像，將待成像樣本沿z軸軸心旋轉(zhuǎn)，這等效于FFL 磁場(chǎng)繞z軸旋轉(zhuǎn).

沿y軸線方向的FFL 磁場(chǎng)可以用梯度矩陣H(x)=Gx來(lái) 描述. FFL 由均勻移位磁場(chǎng)Hs=Hx控制在x軸方向移動(dòng). 粒子沿著z軸被激勵(lì)場(chǎng)He=Hz激發(fā). 這里同樣假設(shè)SPIO 粒子響應(yīng)是瞬時(shí)的，沒(méi)有弛豫效應(yīng)，疊加后的磁場(chǎng)可表示為[12]

式中：x為 掃描儀坐標(biāo)系上的位置向量；Gxx=?Gzz為滿足麥克斯韋方程組.

要確定FFL 在空間中的位置，需要計(jì)算磁場(chǎng)的幅值（為簡(jiǎn)化結(jié)果，公式中為該幅值的平方值）:

其中， θ為旋轉(zhuǎn)角度. 對(duì)式（8）進(jìn)行整理得到零磁場(chǎng)的條件是：z′=?Hz/Gzz和x′cosθ+y′sinθ=Hx/Gzz. 定義投影線:

由此，求解出FFL 處于

并由Hz激發(fā)場(chǎng)和梯度場(chǎng)決定其在z軸的位置. 注意，物理上l表示垂直于FFL 的移位位置. SPIO 信號(hào)沿著這條線進(jìn)行積分，從而可以任何旋轉(zhuǎn)角度獲取樣品的投影圖像. 通過(guò)調(diào)節(jié)變量Hx,Hz和 θ以精確控制并收集一組完整的投影數(shù)據(jù).

2 磁性納米粒子成像稀疏采樣及三維重構(gòu)算法

2.1 模型假設(shè)

不同于系統(tǒng)矩陣法，X-Space 重構(gòu)算法的主要思想是在接收到信號(hào)的空間中進(jìn)行直接編碼. MPI 需要強(qiáng)大的磁場(chǎng)梯度使FFP 以外所有位置的粒子處于飽和狀態(tài)，進(jìn)而進(jìn)行空間編碼. 換而言之，X-Space 重構(gòu)依賴于3 個(gè)主要假設(shè)：1) 磁性納米粒子絕熱地與所施加的磁場(chǎng)對(duì)齊；2) FFP 位置是唯一的；3) 由于對(duì)基波進(jìn)行濾波而導(dǎo)致的低頻信息丟失是可恢復(fù)的[1,6].除此之外，MPI 圖像還具有連續(xù)性和非負(fù)性兩個(gè)先驗(yàn)條件[9].

2.2 3D TV 稀疏約束算子

由于醫(yī)學(xué)圖像一般具有平滑性，即圖像多具有接近的灰度分布僅在某些結(jié)構(gòu)的邊緣產(chǎn)生灰度急劇變化，TV 稀疏的約束可以看作是一種自適應(yīng)的正則化形式，3D TV 準(zhǔn)則可表示為

其中， D1x、 D2x、 D3x分別為沿i方向、j方向以及k方向的差分算子，約束并在保障三維圖像邊緣的稀疏度，對(duì)于位于遠(yuǎn)離邊緣的平滑區(qū)域的體素，強(qiáng)正則化可以濾除振蕩和高頻分量. 但是對(duì)于位于邊緣的體素，TV 正則化濾除很弱. 因此，對(duì)于邊緣上的體素，高頻分量不會(huì)被過(guò)濾掉，在重建過(guò)程中保留了邊緣.與l2正則化相比，TV 正則化的這種邊緣保持特性實(shí)現(xiàn)了更高的空間分辨率[13]，且具有較好的噪聲抑制能力.

為了實(shí)現(xiàn)快速重構(gòu)，本文設(shè)計(jì)了稀疏矩陣形式的3D TV 變換函數(shù)，可以通過(guò)Matrix Free 方法和一階快速鄰近梯度重構(gòu)算法軟件實(shí)現(xiàn). 稀疏3D TV 算子的示意圖如圖1 所示，分別為x、y、z三個(gè)維度和完整3D 算子的可視化圖，黑色和白色位置為進(jìn)行全變分計(jì)算應(yīng)對(duì)的稀疏點(diǎn)位置.

圖1 稀疏3D TV 算子示意圖Fig. 1 Diagram of the sparse 3D TV operator

2.3 具有噪聲魯棒性的稀疏采樣與重構(gòu)算法

本文將MPI 投影重構(gòu)轉(zhuǎn)換為凸優(yōu)化重構(gòu)形式，3D NRSS-MPI 優(yōu)化問(wèn)題可以描述為

在MPI 中，函數(shù)g(x)可以是空操作，也可以選擇利用一些先驗(yàn)信息，如圖像連續(xù)性等. 在本文3D NRSS-MPI 算法中，為了實(shí)現(xiàn)對(duì)MPI 含噪數(shù)據(jù)進(jìn)行優(yōu)化重構(gòu)，設(shè)置了正則化函數(shù)為

其中，TV 范數(shù)具有n維度，D 表示一階有限差分，那么具體的目標(biāo)函數(shù)可定義為

進(jìn)一步通過(guò)任意鄰近梯度算法[14]求解重構(gòu)圖像:

換而言之，NRSS-MPI 方法解決了一個(gè)凸優(yōu)化問(wèn)題，試圖最小化兩個(gè)數(shù)據(jù)一致性項(xiàng)懲罰（模擬掃描和采集數(shù)據(jù)之間）和全變分稀疏域懲罰，其中MPI 成像模型用Matrix Free 線性變換表示. 在進(jìn)行三維MPI成像時(shí)，根據(jù)如圖2 所示，本文方法提供了基于廣義Radon 變換的MPI 正向模型和3D TV 稀疏算子，可通過(guò)鄰近梯度優(yōu)化算法快速實(shí)現(xiàn)三維圖像重構(gòu).

3 模型實(shí)驗(yàn)與性能分析

3.1 欠采樣下三維MPI 圖像重構(gòu)實(shí)驗(yàn)與分析

在磁性納米粒子成像系統(tǒng)中的重構(gòu)誤差通常表現(xiàn)為條狀偽影. 偽影會(huì)降低圖像的有效分辨率，影響圖像的清晰度. 為了避免條紋偽影和分辨率損失，傳統(tǒng)方法通常選擇符合Nyquist 標(biāo)準(zhǔn)的完備投影采樣集，但MPI 掃描和重構(gòu)所需時(shí)間卻大幅提升，這限制了許多MPI 應(yīng)用場(chǎng)景，如活體MPI 成像，在漫長(zhǎng)掃描過(guò)程中任何移動(dòng)誤差都將嚴(yán)重影響成像質(zhì)量.

為了分析本文3D NRSS-MPI 算法和傳統(tǒng)FBP 算法的性能，以及三維重構(gòu)結(jié)果對(duì)空間分辨率的影響，進(jìn)行了點(diǎn)陣模型實(shí)驗(yàn). 點(diǎn)陣模型一組豎條紋中間間隔為1.5 mm，每組條紋之間間隔為12 mm，兩組之間最鄰近距離為7.5 mm. 根據(jù)奈奎斯特采樣定律可以推導(dǎo)出避免混疊所需的投影數(shù)量[12]

式中：FOVxy為FOV 最小值，根據(jù)待成像實(shí)驗(yàn)設(shè)定；Kmax[m?1]為系統(tǒng)中最大的空間頻率，由系統(tǒng)PSF決定.

通過(guò)濾波反投影（filtered back projection，F(xiàn)BP）方法和NRSS-MPI 方法對(duì)點(diǎn)陣模型成像并欠采樣重構(gòu)，實(shí)驗(yàn)結(jié)果如圖3 所示，在1/3 完全投影數(shù)情況下，F(xiàn)BP 算法重構(gòu)圖像質(zhì)量較低，有星狀偽影，但依舊可以分清兩條細(xì)豎線. 在1/4 完全投影數(shù)情況下，F(xiàn)BP 方法重構(gòu)圖雖然能看清細(xì)豎線，但圖像質(zhì)量已經(jīng)不能用于診斷等后續(xù)應(yīng)用了，NRSS-MPI 算法重構(gòu)圖雖然沒(méi)有明顯星狀偽影，但細(xì)豎線間存在模糊.

圖3 不同采樣數(shù)據(jù)量下本文NRSS-MPI 方法與傳統(tǒng)FBP 方法重構(gòu)結(jié)果Fig. 3 The reconstruction results of NRSS-MPI and FBP with different sampling projections

在成像和重構(gòu)時(shí)間方面，由于NRSS-MPI 方法可以以1/4 欠采樣投影數(shù)進(jìn)行重構(gòu)，這意味著獲取圖像時(shí)間縮小了4 倍. 另外在重構(gòu)的計(jì)算上，NRSS-MPI方法花了10 s 重建了一個(gè)3D 版本的模型，在實(shí)際應(yīng)用中證明了該算法的快速性.

3.2 噪聲魯棒性實(shí)驗(yàn)與分析

為了定量評(píng)價(jià)重構(gòu)圖像的相似性和分析采樣數(shù)據(jù)含噪聲情況下CSMPI 算法的魯棒性質(zhì)，用SSIM指標(biāo)對(duì)RSN分別為5、10、20 下的含噪聲采樣數(shù)據(jù)通過(guò)CS-MPI 方法和FBP 方法重構(gòu)圖進(jìn)行評(píng)價(jià)測(cè)試.

式中： μx和 μy分 別代表x、y的平均值； σ和 σ分別代表x、y的方差； σxy為協(xié)方差；c1和c2為常數(shù)，對(duì)比結(jié)果如表1 所示.

表1 不同含噪聲數(shù)據(jù)在本文NRSS-MPI 方法和傳統(tǒng)FBP 方法重構(gòu)圖結(jié)構(gòu)相似度比較表Tab. 1 The comparison of SSIM between the proposed NRSSMPI and FBP with different noise level

從SSIM 水平而言，NRSS-MPI 優(yōu)于FBP 重建，尤其是在RSN較低的情況下. 當(dāng)噪聲水平增加時(shí)，可以觀察到NRSS-MPI 和FBP 的SSIM 收斂水平之間的差距增大，提出的CS-MPI 方法在重構(gòu)圖像的結(jié)構(gòu)相似性分析中表現(xiàn)優(yōu)越.

此外，對(duì)不同信噪比下FBP 算法和3D NRSSMPI 算法1D 重構(gòu)輪廓進(jìn)行了對(duì)比，分析本文NRSSMPI 算法對(duì)含噪聲采樣數(shù)據(jù)的重構(gòu)能力，仿真實(shí)驗(yàn)結(jié)果如圖4 所示. 對(duì)冠狀血管模型的水平位置（y=40）的最強(qiáng)一幀進(jìn)行1D 輪廓勾畫(huà)，圖中（x=19,x=40）位置是冠狀動(dòng)脈的兩條分支血管. 在RSN=20 情況下，3D NRSS-MPI 和FBP 算法總體上都能重構(gòu)出兩處極大值位置，但是在1/4 欠采樣時(shí)，F(xiàn)BP 方法出現(xiàn)了比較明顯的尾部波動(dòng)，即產(chǎn)生了背景偽影，而3D NRSSMPI 方法依然可以準(zhǔn)確重構(gòu). 在RSN=5 的采樣噪聲較大的情況下，F(xiàn)BP 在不同采樣投影數(shù)的實(shí)驗(yàn)中都存在比較明顯的波動(dòng)，3D NRSS-MPI 算法明顯比FBP穩(wěn)定、準(zhǔn)確. 1D 輪廓的重構(gòu)比較更直觀的看出FBP算法在背景處噪聲波動(dòng)大于3D SS-MPI 算法，尤其在含噪聲數(shù)據(jù)重構(gòu)實(shí)驗(yàn)時(shí)FBP 算法重構(gòu)數(shù)據(jù)與真值存在明顯差距. 再一次驗(yàn)證了3D NRSS-MPI 算法能夠快速、穩(wěn)健地重建欠采樣、噪聲較大的FFL-MPI數(shù)據(jù).

圖4 不同含噪數(shù)據(jù)下本文NRSS-MPI 方法與傳統(tǒng)FBP 方法重構(gòu)結(jié)果1D 輪廓圖Fig. 4 The comparison of 1D profile of the reconstruction results between the proposed NRSS-MPI and FBP with different noise level

3.3 三維成像與重構(gòu)實(shí)驗(yàn)

用未經(jīng)稀釋的Resovist 納米粒子示蹤劑制作了5 點(diǎn)源模型，然后將它們固定在設(shè)計(jì)定制的3D 支架中，使得在3D 點(diǎn)源模型在x、y、z視角方向具有不同間隔值. 此外，設(shè)計(jì)并用激光切割機(jī)（universal laser systems, Inc., Scottsdale, AZ, USA）制作雕刻了3 mm厚的點(diǎn)陣模型和冠狀動(dòng)脈血管樣待測(cè)模型，用膠帶進(jìn)行密封用針頭注入稀釋后的Resovist 納米粒子示蹤劑，針孔處用透明指甲油密封防止SPIO 蒸發(fā)或流出污染掃描平臺(tái). 用3D 打印的模型支架將血管模型呈60°角放置形成3D 模型，如圖5（a）（c）所示. 最終待測(cè)模型通過(guò)3D 打印模型底托準(zhǔn)確固定在成像孔，并加載到高分辨率FFL-MPI 掃描儀中，成像結(jié)果如圖5（b）（d）所示，實(shí)驗(yàn)證明了本文算法在實(shí)際應(yīng)用中是快速有效的.

圖5 三維成像模型及重構(gòu)圖Fig. 5 3D Phantoms and the reconstruction results

4 結(jié) 論

本文構(gòu)建了FFL-MPI 前向模型并提出了一種新的具有噪聲魯棒性的稀疏采樣MPI 三維重構(gòu)算法，該算法通過(guò)最小化模擬掃描和采集數(shù)據(jù)之間的差異性和3D 全變分域上的差異，然后通過(guò)鄰近優(yōu)化算法求解重構(gòu)出MPI 三維圖像. 所提出的MPI 成像及重構(gòu)模型和快速稀疏的3D TV 算子充分利用了MPI 成像的先驗(yàn)信息，具有良好的噪聲魯棒性. 通過(guò)成像實(shí)驗(yàn)和含噪數(shù)據(jù)魯棒性關(guān)鍵指標(biāo)測(cè)試分析了濾波反投影FBP 方法和本文NRSS-MPI 方法在不同投影數(shù)下的重構(gòu)圖像質(zhì)量，結(jié)果表明，當(dāng)FBP 對(duì)約1/4 欠采樣投影進(jìn)行圖像重構(gòu)時(shí)會(huì)產(chǎn)生明顯的偽影，但NRSSMPI 方法能正確地重建冠狀動(dòng)脈體模. 在不同含噪水平的數(shù)據(jù)重建中，本文NRSS-MPI 方法重構(gòu)結(jié)果的結(jié)構(gòu)相似性優(yōu)于傳統(tǒng)方法. 此外，還通過(guò)三維成像實(shí)驗(yàn)證明了本文算法在實(shí)際應(yīng)用中是快速有效的.