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

（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 中用以顯著提升成像速度，稱(chēng)其為CS?MRI。近年來(lái)，在磁共振成像領(lǐng)域越來(lái)越多學(xué)者開(kāi)始使用冗余變換。緊標(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 重建問(wèn)題。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ò)來(lái)解決并行磁共振成像重建問(wèn)題。但以上算法中采用的軟閾值函數(shù)收縮功能較差，導(dǎo)致成像速度慢、質(zhì)量差。

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

針對(duì)pFISTA 使用的軟閾值函數(shù)收縮功能較差的問(wèn)題，本文利用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)過(guò)某個(gè)變換后是稀疏的，因此也對(duì)應(yīng)有兩類(lèi)重建模型：綜合型和分解型［16］。

綜合型

分解型

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

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

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

1.2 迭代軟閾值投影算法

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

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

對(duì)式（4）中的最小化問(wèn)題應(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閾值算法是求解正則化逆問(wèn)題的有效算法。

文獻(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］，在圖像重建過(guò)程中，通過(guò)p閾值函數(shù)與小系數(shù)之間的比較來(lái)抑制小系數(shù)進(jìn)行去噪，從而獲得更好的重建圖像。當(dāng)p=1 時(shí)，變成了軟閾值算子［4］，當(dāng)p→∞時(shí)，為硬閾值算子。通過(guò)保留大系數(shù)來(lái)克服軟閾值函數(shù)容易失真的問(wèn)題，通過(guò)收縮中間系數(shù)來(lái)減少硬閾值函數(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ù)的近端映射可以解決正則化反問(wèn)題，因此軟閾值函數(shù)作為相應(yīng)的懲罰函數(shù)的近端映射十分有效。p閾值函數(shù)可以看作是具有稀疏約束更寬泛的懲罰函數(shù)的映射。針對(duì)軟閾值函數(shù)收縮功能較差的問(wèn)題，本文提出將收縮性能更好的p閾值函數(shù)替換pISTA 中的軟閾值，并設(shè)計(jì)了相關(guān)的算法。

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

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

通過(guò)靈活改變p值來(lái)設(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 的西門(mén)子 Trio Tim 磁共振成像儀掃描一個(gè)健康的志愿者頭部得到的圖像，成像儀使用32 線(xiàn)圈，脈沖序列為T(mén)2 加權(quán)的Turbo 自旋回波序列，兩幅圖像分別來(lái)自第7、10 層，其中TR/TE=6 100/99 ms，F(xiàn)OV=220 mm×220 mm，厚度為 3 mm。其中圖2（c）模擬三維成像中的二維相位編碼，圖2（d）模擬二維成像中的偽放射性線(xiàn)采樣，它的采樣點(diǎn)都是在笛卡爾坐標(biāo)的網(wǎng)格上最接近真實(shí)放射線(xiàn)采樣軌跡的點(diǎn)。k空間數(shù)據(jù)的中心區(qū)域信噪比較高，遠(yuǎn)離中心區(qū)域的噪聲較為凸顯，因此CS?MRI 中的采樣模板在k空間的中心采集更多的數(shù)據(jù)。以上所有數(shù)據(jù)來(lái)源于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 為T(mén)2 加權(quán)腦圖第7 層圖像重建結(jié)果，其中第1 行為原圖和4 種算法得到的重建圖像，第2 行第1幅為30%偽放射線(xiàn)采樣模板，剩余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)的收斂曲線(xiàn)如圖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)用圖像的局部邊緣方向性來(lái)進(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 層圖像收斂曲線(xiàn)Fig.5 Convergence curve of T2-weighted brain images with #7

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

圖7 Contourlets 下的收斂曲線(xiàn)Fig.7 Convergence curve using Con?tourlets

因?yàn)閜FIpTA 求解的是平衡型模型，而非精確的分解型模型，這部分通過(guò)對(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 種算法相比重建誤差更??；在50%采樣率下，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í)的收斂曲線(xiàn)，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 下所代表的曲線(xiàn)，也表明在0.7 時(shí)重建精度更好一點(diǎn)。綜合考慮，本文建議將p的值設(shè)置為0.7。

圖8 PFIpTA 不同p 值的收斂曲線(xiàn)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é)束語(yǔ)

針對(duì)迭代軟閾值投影算法中軟閾值函數(shù)收縮功能較差的問(wèn)題，本文提出了迭代p閾值投影算法和它的加速版快速迭代p閾值投影算法來(lái)求解緊標(biāo)架下的壓縮感知磁共振圖像重建問(wèn)題，將收縮性能更好的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)效果。