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

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

0 引言

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

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

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

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

在稀疏重構(gòu)問(wèn)題中，高效的稀疏重構(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）是一種簡(jiǎn)單而有效的凸優(yōu)化算法，其主要思想是通過(guò)分裂變量將無(wú)約束優(yōu)化問(wèn)題轉(zhuǎn)化為約束問(wèn)題，然后交替迭代求解。該算法被廣泛應(yīng)用于求解各種稀疏優(yōu)化模型［21］，本文也考慮利用ADMM求解新提出的模型。

本文對(duì)核磁共振圖像的重構(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)化算法來(lái)有效地求解新模型；最后利用仿真實(shí)驗(yàn)驗(yàn)證了所提模型的有效性。實(shí)驗(yàn)結(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 對(duì)于一個(gè)半連續(xù)函數(shù)f，定義其Moreau 包絡(luò)函數(shù)為：

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

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

圖1（a）為二維L2范數(shù)的圖像，圖1（b）為對(duì)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ù)的非凸性是固定的，為了更好地控制其非凸性，可以引入一個(gè)參數(shù)b(b ＞0)來(lái)構(gòu)造L2范數(shù)的尺度Moreau包絡(luò)和尺度MC罰函數(shù)。

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

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

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

其中：X ∈ Rn×n表示與圖像x相對(duì)應(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，對(duì)于各向同性的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）在一定條件下是一個(gè)凸優(yōu)化問(wèn)題，因此可以利用已有的凸優(yōu)化算法進(jìn)行求解，下一章將考慮利用ADMM算法來(lái)求解該問(wèn)題。

4 ADMM算法

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

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

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

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

步驟1 更新xk+1：

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

步驟2 更新zk+1：

與參考文獻(xiàn)［8，17］類似，式（18）可通過(guò)式（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）計(jì)算uk+1

10） End Do

輸出：x

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

5.1 實(shí)驗(yàn)設(shè)置

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

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

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

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

5.3 結(jié)果分析

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

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

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

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

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

實(shí)驗(yàn)第二部分展示了采樣率對(duì)核磁共振圖像重構(gòu)結(jié)果的影響。在變密度隨機(jī)采樣模板下，通過(guò)改變采樣率，得到四種模型重構(gòu)的圖像與原始大腦圖像的相對(duì)誤差值的對(duì)比結(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)模型下的相對(duì)誤差與采樣率的比較Fig. 7 Comparison of relative error and sampling rates of brain image under four reconstruction models

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

6 結(jié)語(yǔ)

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