基于迭代p 閾值投影算法的壓縮感知磁共振成像

杜秀麗 ，劉晉廷 ，呂亞娜 ，邱少明

（1.大連大學(xué)通信與網(wǎng)絡(luò)重點(diǎn)實(shí)驗(yàn)室，大連，116622；2.大連大學(xué)信息工程學(xué)院，大連，116622）

引 言

磁共振成像（Magnetic resonance imaging，MRI）［1］在臨床診斷中扮演著重要角色。盡管MRI 可以提供具有優(yōu)異軟組織對(duì)比度信息的高質(zhì)量圖像，但其成像速度不盡如人意，主要是因?yàn)槭艿轿锢砩希ㄈ缣荻让}沖的幅度和變化率）以及生理上（神經(jīng)刺激）的約束。壓縮感知技術(shù)（Compressed sensing，CS）［2］被引入到MRI 中用以顯著提升成像速度，稱其為CS?MRI。近年來，在磁共振成像領(lǐng)域越來越多學(xué)者開始使用冗余變換。緊標(biāo)架［3］的冗余變換提高了對(duì)復(fù)雜圖像的稀疏表示能力。文獻(xiàn)［4］提出了一個(gè)簡(jiǎn)單而有效的基于緊標(biāo)架的壓縮感知磁共振圖像重建算法——快速迭代軟閾值投影算法（projected fast iterative soft?thresholding algorithm，pFISTA），它借鑒了快速迭代軟閾值算法（Fast iterative soft?thresh?olding algorithm，F(xiàn)ISTA），利用編程使得pFISTA 避免存儲(chǔ)整個(gè)冗余系統(tǒng)下的表示系數(shù)，因此pFISTA可以處理高冗余度變換下的大規(guī)模MRI 重建問題。pFISTA 的另一個(gè)優(yōu)點(diǎn)是重建誤差對(duì)步長(zhǎng)不敏感，因此可以在磁共振圖像重建中廣泛應(yīng)用于不同的緊標(biāo)架。與迭代軟閾值算法（Iterative soft?threshold?ing algorithm，ISTA）［5］、FISTA［6］、光滑化的快速迭代軟閾值算法（Smoothed FISTA，SFISTA）［7］、基于Moreau 包絡(luò)的近似平滑迭代算法（Proximal smoothing iterative algorithm，PSIA）［8］、交替方向乘子法（Alternating direction multiplier method，ADMM）［9］相比取得了比較好的重建效果。文獻(xiàn)［10］進(jìn)一步證明了并行成像版本的pFISTA 的收斂性，特別證明了由pFISTA 求解的兩個(gè)著名的并行成像重建模型SENSE 和SPIRiT 的收斂性。文獻(xiàn)［11］提出了pISTA-SENSE 殘差網(wǎng)絡(luò)來解決并行磁共振成像重建問題。但以上算法中采用的軟閾值函數(shù)收縮功能較差，導(dǎo)致成像速度慢、質(zhì)量差。

文獻(xiàn)［12］提出了一種有效的p閾值算法，使用p閾值替換ISTA 中軟閾值來解決稀疏信號(hào)恢復(fù)問題，并證明了所提出算法的全局收斂性。文獻(xiàn)［13］建立了用于恢復(fù)稀疏信號(hào)的確切條件。Li 等［14］提出用迭代奇異值p閾值算法來解決低秩矩陣恢復(fù)問題，并使用隨機(jī)奇異值分解給出了它的兩個(gè)加速版本。文獻(xiàn)［15］提出了一種基于交替方向法的算法——迭代p閾值算法，用于解決低秩張量恢復(fù)問題，進(jìn)一步證明了迭代p閾值算法的有效性；并得出：最小化p閾值算子的懲罰函數(shù)是經(jīng)典凸L1范數(shù)和非凸Schat?ten?p范數(shù)的良好替代。

針對(duì)pFISTA 使用的軟閾值函數(shù)收縮功能較差的問題，本文利用p閾值對(duì)小系數(shù)的懲罰更大、對(duì)大系數(shù)產(chǎn)生較小偏置的優(yōu)勢(shì)，提出迭代p閾值投影算法（projected iterativep?thresholding algorithm，pIp?TA），并利用文獻(xiàn)［6］中的Nesterov 加速方法，提出快速迭代p閾值投影算法（projected fast iterativep?thresholding algorithm，pFIpTA），以提高有效提高M(jìn)RI 的重構(gòu)精度和重構(gòu)速度。

1 迭代軟閾值投影算法

1.1 緊標(biāo)架下CS?MRI 的重建模型

CS?MRI 對(duì)全采樣k空間數(shù)據(jù)進(jìn)行欠采樣，數(shù)學(xué)模型可以表示為

式中：y∈ CM表示采樣得到的k空間數(shù)據(jù)，U∈ RM×N(M

磁共振圖像本身不具有稀疏性，磁共振圖像的稀疏性有兩種表現(xiàn)：一種是被某個(gè)字典中的原子信號(hào)稀疏表示，另一種是經(jīng)過某個(gè)變換后是稀疏的，因此也對(duì)應(yīng)有兩類重建模型：綜合型和分解型［16］。

綜合型

分解型

式中：Φ表示標(biāo)架，當(dāng)Ψ正交且滿足Φ=Ψ*=Ψ-1時(shí)，綜合型模型和分解型模型重建效果相同。α為標(biāo)架表示的圖像系數(shù)，λ為正則化參數(shù)，權(quán)衡稀疏項(xiàng)和數(shù)據(jù)保真項(xiàng)，標(biāo)架能夠有效地利用圖像的冗余信息實(shí)現(xiàn)圖像信息與噪聲的分離，通過求解‖Ψx‖1L1范數(shù)的凸優(yōu)化問題重構(gòu)出高質(zhì)量的圖像，通過求解范數(shù)的數(shù)據(jù)保真項(xiàng)來擬合噪聲。

為了聯(lián)合緊標(biāo)架下分解型和綜合型模型，文獻(xiàn)［16］同時(shí)給出了平衡型重建模型，即

式中：ΦH表示一個(gè)分解矩陣，當(dāng)ΦHΦ=Ι時(shí)，稱Φ為標(biāo)準(zhǔn)緊標(biāo)架，β為一個(gè)近似參數(shù)。

1.2 迭代軟閾值投影算法

文獻(xiàn)［4］提出，pISTA 先將分解型模型改寫成一個(gè)形似綜合型模型的等價(jià)模型，如式（4）所示，然后計(jì)算目標(biāo)函數(shù)中非光滑項(xiàng)的近鄰映射，并引進(jìn)Beck 和Teboulle 在文獻(xiàn)［6］中提出的加速策略，得到快速迭代軟閾值投影算法。

pISTA 和pFISTA 實(shí)際求解的是平衡型模型式（4），而不是精確的分解型模型。對(duì)于平衡型模型式（4），首先令

對(duì)式（4）中的最小化問題應(yīng)用ISTA 就需要求解如下的近鄰映射（Proximal map）［5］

式中：γ為步長(zhǎng)，αk+1表示第k+1 迭代的系數(shù)，F(xiàn)H為離散傅里葉反變換，UT為U的轉(zhuǎn)置。prox 表示近鄰映射，其定義為

式中：z表示重建圖像，x表示磁共振圖像，f代表一個(gè)基于L1范數(shù)的函數(shù)。對(duì)于綜合型模型，f(z) =‖z‖1，因此近鄰映射即為

式中：Tλ(?)代表對(duì)每一個(gè)數(shù)值計(jì)算的軟閾值函數(shù)，即

式中：sgn(xi)表示符號(hào)函數(shù)，xi表示x其中一個(gè)元素，即當(dāng)xi>0 時(shí)為 1，當(dāng)xi<0 時(shí)為-1；λ表示閾值。

將式（8）的左右兩邊同乘以Φ，同時(shí)令xk+1=ΦTλ(ΦH(xk+γFHUT(y-UFxk)))，可得

式中：xk+1表示第k+1 迭代的圖像。式（12）為pISTA 的迭代公式，再引入FSITA 的加速策略就可以得到pFISTA。

2 迭代p 閾值投影算法

2.1 迭代p 閾值

不論是ISTA 還是FISTA，求解懲罰函數(shù)的近端映射的迭代廣義p閾值算法是求解正則化逆問題的有效算法。

文獻(xiàn)［12］提出一種新的迭代p閾值函數(shù)

式中：p值表示為xi增加權(quán)重的值。p閾值操作函數(shù)如圖1 所示。

從圖1 可以看出，當(dāng)p<1 時(shí)，p閾值函數(shù)比軟閾值函數(shù)對(duì)小系數(shù)的懲罰更大，而對(duì)大系數(shù)產(chǎn)生更小的偏差。噪聲分量在所有頻帶中以小系數(shù)進(jìn)行散射［17］，在圖像重建過程中，通過p閾值函數(shù)與小系數(shù)之間的比較來抑制小系數(shù)進(jìn)行去噪，從而獲得更好的重建圖像。當(dāng)p=1 時(shí)，變成了軟閾值算子［4］，當(dāng)p→∞時(shí)，為硬閾值算子。通過保留大系數(shù)來克服軟閾值函數(shù)容易失真的問題，通過收縮中間系數(shù)來減少硬閾值函數(shù)不連續(xù)。文獻(xiàn)［13］已經(jīng)證明，當(dāng)p<1 時(shí)，稀疏恢復(fù)所需的測(cè)量比使用p=1的ISTA 情況下少，或者使用相同數(shù)量的測(cè)量時(shí)，使用p<1 可獲得更好的質(zhì)量。

圖1 p 閾值函數(shù)圖Fig.1 Plot of p-thresholding function

2.2 迭代p 閾值投影算法

懲罰函數(shù)的近端映射可以解決正則化反問題，因此軟閾值函數(shù)作為相應(yīng)的懲罰函數(shù)的近端映射十分有效。p閾值函數(shù)可以看作是具有稀疏約束更寬泛的懲罰函數(shù)的映射。針對(duì)軟閾值函數(shù)收縮功能較差的問題，本文提出將收縮性能更好的p閾值函數(shù)替換pISTA 中的軟閾值，并設(shè)計(jì)了相關(guān)的算法。

因此，針對(duì)緊標(biāo)架下模型重建問題設(shè)計(jì)迭代p閾值算法。迭代p閾值投影算法求解緊標(biāo)架下模型重建問題模型如式（4）所示，通過1.2 節(jié)了解到迭代軟閾值投影算法的迭代公式可以解決平衡型模型重建問題，將p閾值替換軟閾值，并融入到pISTA 的核心公式中，得到pIpTA 的迭代公式為

式中：Tλp(t)為迭代p閾值函數(shù)，λ表示閾值，γ為步長(zhǎng)。

通過靈活改變p值來設(shè)計(jì)新的稀疏目標(biāo)函數(shù)，以獲得更好的重建效果。同時(shí)，引入文獻(xiàn)［6］提出的Nesterov 加速策略，Nesterov 方法采用從上一步迭代點(diǎn)處朝前走一步處的梯度，以極少的額外計(jì)算量大幅度提高了迭代p閾值投影算法的收斂速度，從而得到pFIpTA。緊標(biāo)架下MRI 重建的pIpTA 和pFIp?TA 算法分別如算法1 和算法2 所示。

算法1pIpTA 算法

參數(shù)：λ，γ

初始化：x0

循環(huán)直到停止：

輸出：x

算法2pFIpTA 算法

參數(shù)：λ，γ

循環(huán)直到停止：

輸出：x

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

3.1 實(shí)驗(yàn)環(huán)境及數(shù)據(jù)集

圖2（a）和圖2（b）中的腦圖均是使用一臺(tái)3T 的西門子 Trio Tim 磁共振成像儀掃描一個(gè)健康的志愿者頭部得到的圖像，成像儀使用32 線圈，脈沖序列為T2 加權(quán)的Turbo 自旋回波序列，兩幅圖像分別來自第7、10 層，其中TR/TE=6 100/99 ms，F(xiàn)OV=220 mm×220 mm，厚度為 3 mm。其中圖2（c）模擬三維成像中的二維相位編碼，圖2（d）模擬二維成像中的偽放射性線采樣，它的采樣點(diǎn)都是在笛卡爾坐標(biāo)的網(wǎng)格上最接近真實(shí)放射線采樣軌跡的點(diǎn)。k空間數(shù)據(jù)的中心區(qū)域信噪比較高，遠(yuǎn)離中心區(qū)域的噪聲較為凸顯，因此CS?MRI 中的采樣模板在k空間的中心采集更多的數(shù)據(jù)。以上所有數(shù)據(jù)來源于http：//csrc.xmu.edu.cn/。

仿真實(shí)驗(yàn)的硬件配置為：因特爾酷睿i5?3360M 處理器、主頻3.4 GHz、內(nèi)存16 GB；軟件配置為Mat?lab 2081b。

利用式（15）定義的相對(duì)L2范數(shù)誤差作為圖像重建質(zhì)量的數(shù)字指標(biāo)。

圖2 實(shí)驗(yàn)數(shù)據(jù)集Fig.2 Experimental datasets

3.2 實(shí)驗(yàn)結(jié)果及討論

3.2.1 pFIpTA 算法性能分析

本文選用平移不變離散小波變換（Shift?invariant discrete wavelet transform，SIDWT）作為典型實(shí)驗(yàn)中的緊標(biāo)架，使用的是具有4 個(gè)分解尺度的Daubechies 小波。對(duì)于實(shí)驗(yàn)中所用到的算法，首先需要對(duì)其調(diào)節(jié)參數(shù)，這些參數(shù)包括：FIpTA 中的步長(zhǎng)γFp，SFISTA 中的步長(zhǎng)γS和光滑參數(shù)μ，pFISTA 中的步長(zhǎng)γp。按如下規(guī)則設(shè)置：γFp=1；根據(jù)文獻(xiàn)［6］，γS=1/(1+1/μ)，而μ也根據(jù)文獻(xiàn)取不同的推薦值；對(duì)于pFISTA，根據(jù)文獻(xiàn)［2］，將γp設(shè)置為1。pFIpTA 算法中p的值為0.7。

圖3 為T2 加權(quán)腦圖第7 層圖像重建結(jié)果，其中第1 行為原圖和4 種算法得到的重建圖像，第2 行第1幅為30%偽放射線采樣模板，剩余4 幅為4 種算法的殘差圖。圖中FIpTA，SFISTA，pFISTA 和pFIp?TA 的 RLNE 分別為 0.132 426，0.097 078，0.097 029 和 0.083 288。圖4 為加權(quán)腦圖第 10 層重建結(jié)果，使用30% 高斯采樣模板進(jìn)行采樣。圖中FIpTA，SFISTA，pFISTA 和pFIpTA 的RLNE 分別為0.104 166，0.091 703，0.091 573 和 0.069 123。從圖3，4 可以看出，F(xiàn)IpTA 的重建圖像含有明顯的偽影，而SFISTA 和pFISTA 的重建圖像中這些偽影被很好地壓制了，pFIpTA 中的偽影最少。pFIpTA 的重建誤差比其他3 種算法的重建誤差要低。對(duì)應(yīng)的收斂曲線如圖5，6 所示，所提出的pFIpTA 比pFIS?TA，SFISTA 提高了大約0.013 79 個(gè)單位。

圖3 T2 加權(quán)腦圖第7 層圖像重建結(jié)果Fig.3 Recconstructed T2-weighted brain images with #7

圖4 T2 加權(quán)腦圖第10 層圖像重建結(jié)果Fig.4 Recconstructed T2-weighted brain images with #10

不同的稀疏化磁共振圖像的緊標(biāo)架對(duì)CS?MRI 的重建誤差影響很大。為了全面評(píng)估pFIpTA，在本節(jié)實(shí)驗(yàn)中選用輪廓波（Contourlets）作為緊標(biāo)架，Contourlets 應(yīng)用圖像的局部邊緣方向性來進(jìn)一步稀疏化圖像，使得重建結(jié)果能夠盡可能保留邊緣信息，并和求解分解型模型的通用算法ADMM 進(jìn)行了對(duì)比，實(shí)驗(yàn)結(jié)果如圖7 所示。圖中結(jié)果顯示，和使用SIDWT 作為緊標(biāo)架效果一樣，pFIpTA 比其他的重建誤差都要低，而且pFIpTA 收斂更快。因此，pFIpTA 具備pFISTA 的優(yōu)勢(shì)不隨使用的緊標(biāo)架而變化，并且效果更好。

圖5 T2 加權(quán)腦圖第7 層圖像收斂曲線Fig.5 Convergence curve of T2-weighted brain images with #7

圖6 T2 加權(quán)腦圖第10 層圖像收斂曲線Fig.6 Convergence curve of T2-weighted brain images with#10

圖7 Contourlets 下的收斂曲線Fig.7 Convergence curve using Con?tourlets

因?yàn)閜FIpTA 求解的是平衡型模型，而非精確的分解型模型，這部分通過對(duì)比ADMM，可了解pFIpTA 與分解型重建結(jié)果的效果。從圖7 可知，SFISTA 和pFISTA 重建誤差比ADMM 的重建誤差更小，且本文所提pFIpTA 算法收斂速度更快，精度更高。

為了驗(yàn)證本文提出的重構(gòu)算法在不同采樣率下的磁共振圖像的重建性能，對(duì)不同采樣率下SFIS?TA、pFISTA、ADMM 和pFIpTA 的重建誤差進(jìn)行了分析。表1 是不同采樣率下各算法的重建誤差，實(shí)驗(yàn)分別采用1%，5%，10%，20%，30%，40%和50%的高斯采樣模板，緊標(biāo)架為輪廓波。從表1 可以看出在1%，5%，10%，20%，30%，40%和50%采樣率下，SFISTA 和pFISTA 重建誤差差距不大，pFIp?TA 與其他3 種算法相比重建誤差更?。辉?0%采樣率下，4 種算法重建誤差趨于平穩(wěn)。

表1 不同采樣率下重建誤差Table 1 Reconstruction error at different sampling ratios

3.2.2p值的選擇

本節(jié)將從數(shù)值實(shí)驗(yàn)驗(yàn)證p值如何影響重建速度和重建誤差。圖8 中，選用p值分別為1.0，0.8，0.7，0.6 時(shí)的收斂曲線，p=1 為pFISTA。圖9 為在緊標(biāo)架為SIDWT 和Contourlets 時(shí)，取不同p值對(duì)應(yīng)的RLNE。

從圖8 和圖9 可以看出，隨著p值減小，收斂速度變快，但是RLNE 并不遵循這個(gè)規(guī)律。從圖8 可以看出在0.8，0.7 出現(xiàn)拐點(diǎn)，0.8 和0.7 的收斂精度差不多，但是在0.7 時(shí)，速度有所提升。圖9 中在Contour?lets 下所代表的曲線，也表明在0.7 時(shí)重建精度更好一點(diǎn)。綜合考慮，本文建議將p的值設(shè)置為0.7。

圖8 PFIpTA 不同p 值的收斂曲線Fig.8 Convergence curve of PFIpTA at different p values

圖9 SIDWT 和Contourlets 條件下不同p 值對(duì)應(yīng)的RLNEFig.9 RLNE with different p values under SIDWT and Contourlets

4 結(jié)束語

針對(duì)迭代軟閾值投影算法中軟閾值函數(shù)收縮功能較差的問題，本文提出了迭代p閾值投影算法和它的加速版快速迭代p閾值投影算法來求解緊標(biāo)架下的壓縮感知磁共振圖像重建問題，將收縮性能更好的p閾值函數(shù)替換pISTA 中的軟閾值。數(shù)值實(shí)驗(yàn)顯示，pFIpTA 的重建結(jié)果比SFISTA 和pFISTA的好，而且pFIpTA 收斂速度比pFISTA 更快；與ADMM 相比收斂精度也有一定提高。p閾值具有能靈活的收縮功能，驗(yàn)證了pFIpTA 在p為0.7 時(shí)表現(xiàn)出色。pFIpTA 不僅繼承了pFISTA 不隨使用的緊標(biāo)架而變化的優(yōu)勢(shì)，甚至效果更好，而且繼承了p閾值函數(shù)的優(yōu)勢(shì)，獲得了更好的重構(gòu)效果。