基于非凸全變差正則的核磁共振圖像重構(gòu)算法

沈馬銳，李金城，張 亞，鄒 健

（長江大學(xué)信息與數(shù)學(xué)學(xué)院，湖北荊州434023）

0 引言

核磁共振成像（Magnetic Resonance Imaging，MRI）由于其有效的描繪軟組織變化的能力，被廣泛地運用于醫(yī)學(xué)領(lǐng)域［1］。MRI 應(yīng)用的主要限制因素是其成像的采集時間相對較長，為了加快掃描速度、減少數(shù)據(jù)采集量和采集時間而不降低重建質(zhì)量始終是關(guān)注的焦點。稀疏重構(gòu)模型能夠從高度欠采樣的信號中精確地重構(gòu)圖像，解決由于采集數(shù)據(jù)的限制和采集時間長的問題，因此在核磁共振成像中發(fā)揮著重要的作用［2-3］。

MRI 系統(tǒng)可以抽象為一個欠定線性方程組Ax = y，其中：x ∈RN為待重構(gòu)的核磁共振圖像；y ∈RM是磁共振線圈采集到的k空間信號；A ∈ RM×N（M ＜ N）為欠采樣傅里葉變換。一般來說 A = P × F，P ∈ RM×N為欠采樣模板，如徑向采樣模板、變密度隨機采樣模板、笛卡爾采樣模板等；F ∈ RN×N為稀疏變換算子，如傅里葉變換、小波變換等。一般情況下，上述欠定線性方程組有無窮多解，但可以通過求解如下優(yōu)化問題來唯一重構(gòu)x，求解模型如下：

常用的稀疏正則包括L0正則、L1正則、Lp正則和全變差（Total Variation，TV）正則等。L0正則的一般形式為φ0(x) =，它表示稀疏解x 中非零元素的個數(shù)；但求解L0范數(shù)是個NP 難問題。L1正則是L0正則的凸松弛，目前使用最為廣泛；但L1范數(shù)不可導(dǎo)且是一個有偏估計，對稀疏解中較大的值欠估計［5］。Lp正則雖然能更有效地估計稀疏解x 中的較大值，但由于目標(biāo)函數(shù)本身是非凸的，求解過程往往非常困難，例如會陷入局部最優(yōu)值。為了克服L1和Lp正則的缺點，Selesnick 等［6-9］提出了一類L1非凸正則構(gòu)造方法，其基本思想是利用最小最大凹（Minmax-Concave，MC）罰函數(shù)，在 L1范數(shù)上減去其Moreau 包絡(luò)。該正則項是非凸的，能有效估計稀疏解x中的較大值，同時該非凸正則能使目標(biāo)函數(shù)整體保持凸性，可以利用已有的凸優(yōu)化算法進(jìn)行高效求解。該類非凸正則被廣泛地應(yīng)用于圖像的去噪和重構(gòu)上［10-11］。

TV 正則是另一種常見的稀疏正則，可有效地解決圖像在噪聲情況下的光滑問題［12］。它不僅能抑制圖像的噪聲，而且能保留圖像的紋理、邊緣等細(xì)節(jié)，因此被廣泛地用于信號、圖像的去噪或圖像的重構(gòu)［13-15］。TV 正則可分為各向異性和各向同性兩種情況，通過引入差分矩陣，可以分別轉(zhuǎn)化為L1范數(shù)和 L2范數(shù)［16］。目前對 TV 正則的研究主要以可轉(zhuǎn)化為 L1正則的各向異性情況為主，如：文獻(xiàn)［17］中基于Selesnick 的工作，利用基于L1范數(shù)的最小最大凹全變差正則（Minmax-Concave Total Variation of L1norm，MCTV-L1）來重構(gòu)MR 腦圖像。但正如文獻(xiàn)［18］中所述，各向異性TV模型適用于那些呈現(xiàn)“塊狀”結(jié)構(gòu)的圖像，而對于自然圖像，各向同性TV 模型要更加適合。因此，對可轉(zhuǎn)化L2范數(shù)的各向同性TV模型的研究是很有必要的。

在稀疏重構(gòu)問題中，高效的稀疏重構(gòu)算法是保證稀疏重構(gòu)得到廣泛應(yīng)用的關(guān)鍵，常用的稀疏重構(gòu)算法大致可以分為貪婪算法、凸優(yōu)化算法和神經(jīng)網(wǎng)絡(luò)算法三類［19-20］。交替方向乘子法（Alternating Direction Method of Multipilers，ADMM）是一種簡單而有效的凸優(yōu)化算法，其主要思想是通過分裂變量將無約束優(yōu)化問題轉(zhuǎn)化為約束問題，然后交替迭代求解。該算法被廣泛應(yīng)用于求解各種稀疏優(yōu)化模型［21］，本文也考慮利用ADMM求解新提出的模型。

本文對核磁共振圖像的重構(gòu)模型和算法進(jìn)行研究，提出了一種基于L2正則的非凸全變差正則重構(gòu)模型。首先將Selesnick 構(gòu)造非凸正則的思想應(yīng)用到L2范數(shù)上，利用Moreau包絡(luò)函數(shù)和MC 罰函數(shù)構(gòu)造L2范數(shù)的非凸正則，提出基于L2范數(shù)的非凸TV 正則模型——（Minmax-Concave Total Variation of L2norm，MCTV-L2），并證明了該模型的目標(biāo)函數(shù)在一定條件下是凸的；接著使用ADMM 優(yōu)化算法來有效地求解新模型；最后利用仿真實驗驗證了所提模型的有效性。實驗結(jié)果表明，該模型能更好地重構(gòu)圖像的紋理、邊緣等細(xì)節(jié)，重構(gòu)性能有顯著提升。

1 L2范數(shù)的非凸正則

1.1 L2范數(shù)的Moreau包絡(luò)和MC罰函數(shù)

Moreau 包絡(luò)函數(shù)和MC 罰函數(shù)是構(gòu)造L2范數(shù)的非凸正則常用的構(gòu)造工具。其定義如下：

定義1 對于一個半連續(xù)函數(shù)f，定義其Moreau 包絡(luò)函數(shù)為：

定義2 對于一個半連續(xù)函數(shù)f，定義其MC罰函數(shù)為：

對L2范數(shù)，其非凸正則構(gòu)造過程［22］如下。

圖1（a）為二維L2范數(shù)的圖像，圖1（b）為對L2范數(shù)進(jìn)行平滑化處理后的Moreau 包絡(luò)圖像，平滑后的函數(shù)是原函數(shù)的可微緊下界，具有迭代形式的解析解，并且也是原函數(shù)的解［23］。

圖1（c）為二維L2范數(shù)的MC 罰函數(shù)圖像，從圖像上可以看出MC罰函數(shù)是非凸的。

圖1 L2范數(shù)及其Moreau包絡(luò)和MC罰函數(shù)Fig. 1 L2 norm and its Moreau envelope and MC penalty function

1.2 L2范數(shù)的尺度Moreau包絡(luò)和尺度MC罰函數(shù)

式（5）中MC 罰函數(shù)的非凸性是固定的，為了更好地控制其非凸性，可以引入一個參數(shù)b(b ＞0)來構(gòu)造L2范數(shù)的尺度Moreau包絡(luò)和尺度MC罰函數(shù)。

定義 5 令 b ＞ 0，對于函數(shù) f(x) =‖x ‖2，定義其尺度Moreau包絡(luò)函數(shù)為：

圖2 L2范數(shù)的尺度Moreau包絡(luò)和尺度MC罰函數(shù)Fig. 2 Scaled Moreau envelope and scaled MC penalty function of L2 norm

2 TV范數(shù)的非凸正則

L1正則化模型在圖像重構(gòu)過程中未能考慮圖像中像素點與其鄰近點水平和豎直方向上的關(guān)系，常導(dǎo)致圖像邊緣特征嚴(yán)重缺失和重構(gòu)圖像不夠光滑、包含噪聲等問題，TV 正則充分考慮到圖像的結(jié)構(gòu)特征，能有效抑制圖像的噪聲而且能保留圖像的紋理、邊緣等細(xì)節(jié)，克服階梯效應(yīng)，實現(xiàn)空間分段光滑，使得空間信息得到更加充分的利用，極大提升圖像的重構(gòu)質(zhì)量。TV 正則有各向異性和各向同性兩種情況，可以分別轉(zhuǎn)化為L1范數(shù)和L2范數(shù)。

假設(shè)x ∈RN為一幅圖像，各向異性的情況可以轉(zhuǎn)化為L1范數(shù)，即：

各向同性的情況可以轉(zhuǎn)化為L2范數(shù)，即：

其中：X ∈ Rn×n表示與圖像x相對應(yīng)的矩陣；D ∈ R2×n2為有限差分算子；此外，通篇使用循環(huán)邊界條件，記為

Selesnick 在文獻(xiàn)［6］中提出了式（8）的非凸正則模型——MCTV-L1，文獻(xiàn)［17］將該模型應(yīng)用到MRI 重構(gòu)上。本文將討論式（9）的非凸正則模型——MCTV-L2。

定義7 令b ＞0，對于各向同性的TV 正則，定義其MC罰函數(shù)（MCTV-L2）為：

3 核磁共振稀疏重構(gòu)模型

結(jié)合稀疏重構(gòu)模型式（1）和MCTV-L2正則式（10），可得到基于非凸全變差正則的核磁共振圖像重構(gòu)模型如下：

定理1 令λ ＞ 0，b ＞ 0，x ∈ RN，定義Gb(x)為：

由于式（11）在一定條件下是一個凸優(yōu)化問題，因此可以利用已有的凸優(yōu)化算法進(jìn)行求解，下一章將考慮利用ADMM算法來求解該問題。

4 ADMM算法

利用ADMM求解式（11）的主要步驟如下：

將式（10）代入式（11），可得：

令z = Dx，可得到式（14）的增廣拉格朗日形式如下：

其中：u 是拉格朗日算子（對偶變量）；ρ 是懲罰參數(shù)項。ADMM 算法通過交替對 x、z、u 進(jìn)行更新，找到式（11）的最優(yōu)解。

步驟1 更新xk+1：

對式（16）求偏導(dǎo)，令其偏導(dǎo)等于零，可得：

步驟2 更新zk+1：

與參考文獻(xiàn)［8，17］類似，式（18）可通過式（19）、（20）迭代步驟求解：

其中：

步驟3 更新uk+1：

綜上可得求解模型式（11）的ADMM算法如算法1所示。

算法1 基于各向同性的非凸TV 正則化的圖像重構(gòu)算法。

輸入：y

1） 初始化：x0= 0，z0= 0，u0= 0，λ ＞ 0，ρ ＞ 0，0 ＜ b ＜ ρ，

8） End Do

9） 采用式（23）計算uk+1

10） End Do

輸出：x

5 仿真實驗與結(jié)果分析

5.1 實驗設(shè)置

為了驗證本文MCTV-L2模型的有效性，將該模型與經(jīng)典的各向異性TV 正則和各向同性TV 正則（分別記為TV-L1和TV-L2）以及文獻(xiàn)［17］中的MCTV-L1模型對核磁共振圖像重構(gòu)的結(jié)果進(jìn)行了對比。實驗操作都是在Inter Core i5-4200M、2.50 GHz CPU、4.00 GB 的計算機上進(jìn)行，編程環(huán)境使用Matlab 2018a。實驗使用三種欠采樣模板進(jìn)行測試：徑向采樣模板、變密度隨機采樣模板以及笛卡爾采樣模板，分別對Shepp Logan、大腦以及腦部血管共3 幅核磁共振圖像進(jìn)行重構(gòu)結(jié)果對比，采樣模板和實驗圖像如圖3 所示。為方便比較，所有圖像尺寸都設(shè)為256 × 256 像素大小，灰度取值范圍為0～255，稀疏變換算子F 采用傅里葉變換。本文的實驗代碼和實驗結(jié)果可以從網(wǎng)站鏈接https：//github. com/zj15001/MCTV_L2.git下載。

為保證比較的公平性，算法1 中參數(shù)的設(shè)置參考文獻(xiàn)［17］中提到的最好的參數(shù)設(shè)置條件，即λ= 0.01，b= 0.05/λ，δ1= δ2= 0.000 1，ρ = 150。

5.2 實驗評估標(biāo)準(zhǔn)

本文采用相對誤差（Relative Error）和峰值信噪比（Peak Signal-to-Noise Ratio，PSNR）來評估圖像重構(gòu)的質(zhì)量與準(zhǔn)確性。相對誤差VRE和峰值信噪比VPSNR的計算公式為：

其中：x0為觀測圖像；xr為重構(gòu)圖像；x^ 為x0的均值。在這兩種定量評估標(biāo)準(zhǔn)中，相對誤差的值越小模型重構(gòu)效果越好；而對于峰值信噪比來說，則是值越大，重構(gòu)效果越好。在視覺效果評估中，對于每幅核磁共振觀測圖像分別給出了不同模型的重構(gòu)效果以及重構(gòu)圖像與觀測圖像的差值圖像，以便于更好地進(jìn)行結(jié)果比較分析。

圖3 采樣模板與測試圖像Fig.3 Sampling templates and test images

5.3 結(jié)果分析

本文實驗分為兩個部分，第一部分展示了算法針對不同的核磁共振圖像和三種采樣模板下的重構(gòu)效果，視覺效果如圖4～6所示，數(shù)值比較結(jié)果如表1所示。

首先，采用10 條線的徑向采樣模板，選擇Shepp Logan 原始圖像來對所提出的MCTV-L2重構(gòu)模型進(jìn)行測量，四種模型的重構(gòu)效果與差值圖像如圖4所示?？梢钥吹?，使用TV-L1重構(gòu)后的圖像邊緣比較模糊且伴有重影，而MCTV-L2模型重構(gòu)后的圖像更接近原始圖像，在處理圖像邊緣信息時效果明顯最好。如表1所示，使用TV-L1重構(gòu)和MCTV-L2重構(gòu)后的相對誤差分別為0.36 和0.25，峰值信噪比分別為20.8841 dB 和24.0202 dB，顯然MCTV-L2模型重構(gòu)的實驗數(shù)據(jù)最優(yōu)。

為了驗證MCTV-L2模型在變密度隨機采樣模板下重構(gòu)圖像的質(zhì)量，選擇采樣率為30%的隨機采樣模板對像素大小為256 × 256 的大腦圖像進(jìn)行測試。如圖5 所示：本文提出MCTV-L2模型在重構(gòu)視覺效果上比前幾種模型都好，重構(gòu)后的大腦圖像在視覺效果上與原始觀測圖像最相近，具體體現(xiàn)在大腦圖像中邊緣細(xì)節(jié)的重建以及對差值圖像的對比。而TV-L1模型重構(gòu)效果最差，不僅沒有很好地重構(gòu)圖像，而且還丟失圖像的很多邊緣信息，看起來很模糊。

使用笛卡爾采樣模板對腦部血管圖像進(jìn)行測量，沿相位編碼方向僅用100條線對k空間數(shù)據(jù)進(jìn)行欠采樣，同時保持對頻率編碼方向的完全采樣。四種模型重構(gòu)后的圖像效果與差值圖像如圖6 所示。顯然，和其他測量圖像一樣，從視覺效果上看，新模型MCTV-L2重構(gòu)后的圖像效果優(yōu)于其他三種重構(gòu)方法。實驗結(jié)果顯示，MCTV-L2重構(gòu)后的圖像誤差明顯降低，峰值信噪比提升了大約4 dB。

由表1 可知，無論使用徑向采樣模板、隨機采樣模板還是笛卡爾采樣模板，所有實驗結(jié)果均表明，相較于TV-L1重構(gòu)模型、MCTV-L1重構(gòu)模型以及TV-L2重構(gòu)模型，在核磁共振圖像重構(gòu)評估參數(shù)方面，MCTV-L2重構(gòu)模型效果都是最好的，能得到更低的相對誤差和更高的峰值信噪比，更好地重建平滑區(qū)域。

實驗第二部分展示了采樣率對核磁共振圖像重構(gòu)結(jié)果的影響。在變密度隨機采樣模板下，通過改變采樣率，得到四種模型重構(gòu)的圖像與原始大腦圖像的相對誤差值的對比結(jié)果，如圖7所示。

表1 不同采樣模板下四種模型重構(gòu)結(jié)果Tab. 1 Reconstruction results of four models under different sampling templates

圖4 Shepp Logan重構(gòu)結(jié)果與差值圖像Fig. 4 Reconstruction results and difference images of Shepp Logan

圖5 核磁共振大腦重構(gòu)結(jié)果與差值圖像Fig. 5 Reconstruction results and difference images of MR brain

圖6 腦部血管重構(gòu)結(jié)果與差值圖像Fig. 6 Reconstruction results and difference images of brain vessels

圖7 大腦圖像在四種重構(gòu)模型下的相對誤差與采樣率的比較Fig. 7 Comparison of relative error and sampling rates of brain image under four reconstruction models

在不同采樣率下，MCTV-L2模型重構(gòu)后的圖像與原始圖像的相對誤差值均低于另外幾種模型重構(gòu)后的圖像與原始圖像的相對誤差值，且采樣率越高，MCTV-L2重構(gòu)后的圖像與原始圖像的相對誤差值越低，在視覺效果上與原始圖像更接近。實驗結(jié)果表明，MCTV-L2模型具有非常不錯的重構(gòu)效果，與理論推斷相一致。

6 結(jié)語

本文針對核磁共振圖像重構(gòu)中由于欠采樣導(dǎo)致重構(gòu)圖像不夠完整、邊緣模糊以及噪聲殘留等問題，以Moreau 包絡(luò)和最小最大凹罰函數(shù)為工具提出了一種基于L2正則的非凸全變差正則重構(gòu)模型。相較于已有的重構(gòu)模型，所提模型不僅可以有效地避免凸正則中對較大非零元欠估計現(xiàn)象，能夠更有效地處理圖像的紋理、邊緣等細(xì)節(jié)，同時在一定條件下可以保證目標(biāo)函數(shù)的整體凸性。在仿真實驗中，利用高效的ADMM 算法對實驗?zāi)Ｐ瓦M(jìn)行求解，將其運用于核磁共振圖像的重構(gòu)中。實驗結(jié)果表明，相較于其他幾種經(jīng)典模型，所提模型MCTV-L2在重構(gòu)效果上有明顯的提升，能獲得更高的峰值信噪比，視覺效果改善顯著，能有效地消除圖像偽影并在保留圖像邊緣輪廓信息基礎(chǔ)上使圖像更加光滑，具有非常不錯的效果。由于基于各向同性的非凸全變差正則對于自然圖像的處理效果要優(yōu)于各向異性的方法，在進(jìn)一步研究中，我們會考慮將該模型應(yīng)用到其他圖像處理領(lǐng)域。