董正山，林 耿

(閩江學(xué)院數(shù)學(xué)與數(shù)據(jù)科學(xué)學(xué)院(軟件學(xué)院)， 福建 福州 350108)

0 引言

圖像在一些變換(如小波變換、curvelet變換等)下是可以稀疏表示或壓縮的，從而可以用稀疏編碼的方法研究圖像處理中涉及到的問題[1-2]，如可以用稀疏優(yōu)化的方法分析圖像的超分辨率、圖像的去噪、以及圖像恢復(fù)等問題。

假設(shè)I表示圖像，Ψ表示稀疏變換，且圖像I在該變換下是稀疏或可壓縮的，即有I=Ψx，其中，x有較多的零分量(稀疏性)或較多較小的分量(可壓縮性)，Φ表示采樣矩陣，從而觀察到的數(shù)據(jù)可以表示為y=ΦΨx+σ，其中σ表示噪聲。為了從觀察數(shù)據(jù)y中恢復(fù)原始圖像I，根據(jù)問題的稀疏性特點，可建立如下求解模型：

其中,‖·‖0表示向量中非零元素的個數(shù)，其被稱為零?；颉傲惴稊?shù)”，s>0控制變量的稀疏性。通過模型(1)求解出變量x后，就可以通過I=Ψx恢復(fù)出原始圖像。但是問題(1)的最優(yōu)解主要任務(wù)是找準(zhǔn)非零分量的位置，這是一個組合優(yōu)化問題，且是NP-hard問題[3]。為了求解該稀疏優(yōu)化問題，目前主要的方法有貪心算法、松弛算法、混合網(wǎng)模型算法、基于零模的算法等。

主要的貪心算法有匹配追蹤(MP)算法[4-6]等，這些算法可以用于準(zhǔn)確快速地解決無噪聲的小規(guī)模問題，但不適于解決有噪的或大規(guī)模的問題。松弛算法主要是為解決難以處理的(離散的非凸的)“零范數(shù)”問題。松弛算法又分為凸松弛和非凸松弛，代表分別有(重加權(quán))l1范數(shù)松弛[7-13]和lp“范數(shù)”(0

非凸松弛算法主要是lp“范數(shù)”模型算法。與l1范數(shù)模型算法相比，lp“范數(shù)”模型可以得到更加稀疏的結(jié)果，所需采樣率比l1范數(shù)模型低。但由于是非凸優(yōu)化，求的解只是算法的一個鞍點或不動點，算法設(shè)計可能更復(fù)雜。

混合網(wǎng)模型算法主要是同時考慮多種范數(shù)模型的優(yōu)化算法，如混合范數(shù)模型同時考慮兩種范數(shù)(l1和l2范數(shù))模型等[16]。這些算法結(jié)合了不同范數(shù)算法的優(yōu)點，取得了較好的效果。但這些方法還是沒有直接求解“零范數(shù)”問題模型，得到的解可能并非稀疏解。

基于“零范數(shù)”的算法直接求解“零范數(shù)”優(yōu)化模型，可以得到更稀疏的解。這類算法的代表有基于硬閾值算子的迭代方法[17-22]。由于硬閾值算子可以得到更稀疏的解，本文將利用迭代硬閾值方法求解問題(1)的子問題。

1 模型分析

1.1 拉格朗日對偶問題

本節(jié)研究問題(1)的對偶問題。考慮到問題(1)的拉格朗日函數(shù)

(2)

以及對偶函數(shù)

(3)

(4)

根據(jù)這兩個性質(zhì)，本文將設(shè)計基于拉格朗日和硬閾值迭代的求解算法。

1.2 硬閾值迭代算法

硬閾值迭代算法可以用于求問題(3)中最小化問題的解，即給定λ值時求對偶函數(shù)的值。求對偶函數(shù)值的關(guān)鍵是求解“零范數(shù)”正則化問題，雖然該問題也是一個NP難問題，但可以使用迭代硬閾值算法求解。在迭代過程中，該方法為了使得問題可解，在當(dāng)前點近似f(x)，即

(5)

從而關(guān)于x的最小化問題可以變成變量可分離結(jié)構(gòu)，對偶函數(shù)中最小化問題即可近似用下面子問題求解

(6)

(7)

通常L取函數(shù)f(x)的李普希茲常數(shù)上界，但有時候這個上界不易求解，下面主要通過一個線搜索算法[19]找到一個合適的Lk，算法1中給出了求解對偶函數(shù)中最小化問題的迭代硬閾值算法。

算法1 線搜索的迭代硬閾值算法

輸入：x0,λ,L0,Lmin,Lmax

輸出：x*

1)初始化，k=0,γ>1；

2)通過線搜索找到Lk；

4)如果達(dá)到終止條件，終止算法，否則轉(zhuǎn)到2。

2 基于拉格朗日方法的圖像恢復(fù)算法

本節(jié)中介紹利用拉格朗日方法和迭代硬閾值算法求解問題(4)。根據(jù)性質(zhì)1，一個較大的λ值，可以得到更稀疏的解；相反，一個較小的λ得到一個相對稠密的解。所以在迭代過程中，發(fā)現(xiàn)當(dāng)前解的非零元素的個數(shù)大于s，即所得的解較稠密。根據(jù)性質(zhì)2，這時候適當(dāng)?shù)卦黾应说闹担沟酶潞蟮慕獗M量滿足原問題的約束；若在迭代過程中，當(dāng)前解的非零元素的個數(shù)小于s，即所得的解可能過于稀疏，此時需要減少λ的值。這樣通過這種搜索的方法，最終可以得到滿足原問題的一個稀疏解。求解問題(1)的算法框架流程如下。

算法2 基于拉格朗日方法的圖像恢復(fù)算法(LIHT)

輸入：x0,λ0,L0,Lmin,Lmax

輸出：x*

1)初始化，k=0,γ>1；

2)調(diào)用算法1獲得當(dāng)前解xk；

3)計算當(dāng)前解非零元素的個數(shù)nk，如果nk>s，增加λ的值，否則減小λ的值；

4)如果滿足終止條件，終止算法，否則轉(zhuǎn)到2。

文[23]中指出，在一定的條件下，算法2可以得到原問題的一個L-穩(wěn)定點。

3 實驗

本節(jié)通過在不同的采樣率下恢復(fù)CT圖像來說明本文提出算法的有效性。本文實驗都是在Windows10系統(tǒng)(Intel(R) Xeon(R) W-10885M CPU @ 2.40GHz，內(nèi)存32G)中使用工具M(jìn)atlab完成。

實驗采用標(biāo)準(zhǔn)的Shepp-Logan圖像來完成CT圖像的恢復(fù)，Shepp-Logan圖像的大小是256×256(圖1中左上角的第一張圖是原圖)。本實驗中用不同的徑向切片在星形域上通過二維的離散傅里葉變換采樣層析數(shù)據(jù)，并用逆Haar小波變換作為稀疏化變換Ψ,設(shè)定s=3 769。圖1、圖2給出了本文提出的算法LIHT與其它算法(文[20]中的AIHTCG, 文[21]中的DORE,文[22]中的NIHT，文[24]中基于l1范數(shù)的方法FPCAS)的實驗對比，包括算法運行時間以及PSNR值。

圖1 各比較算法恢復(fù)CT圖像時運行時間及恢復(fù)圖像的PSNR值Fig.1 Running times and PSNR value of all compared method on CT images

圖1展示了采樣率為0.176時的結(jié)果對比。從圖1中可以看出，本文提出的算法LIHT獲得了最大PSNR值，雖然算法時間上略微比AIHTCG和DORE多一些，但可以接受。LIHT比NIHT和FPCAS在時間和PSNR上都有優(yōu)勢。圖2展示了在不同的采樣率下各對比算法取得的PSNR值，相應(yīng)算法的運行時間在圖3中給出。對比發(fā)現(xiàn)，在采樣率較小時，F(xiàn)PCAS有些PSNR值的優(yōu)勢，但算法耗時較多。當(dāng)采樣率適當(dāng)增大時，本文給出的算法更有PSNR值的優(yōu)勢。

圖2 各比較算法在不同采樣率下恢復(fù)CT圖像時的PSNR值Fig.2 PSNR value of all compared method on CT images reconstruction under different sampling ratio

圖3 各比較算法在不同采樣率下恢復(fù)CT圖像時的運行時間Fig.3 Running time (second) of all compared method on CT images reconstruction under different sampling ratio

4 結(jié)語

首先，本文通過分析建立的稀疏約束問題的對偶問題，找到從對偶問題出發(fā)求解原約束問題的方法。該方法結(jié)合迭代硬閾值算法和對偶問題的特點，設(shè)計出有效的求解算法；最終通過數(shù)值實驗驗證了本文算法的有效性。后期將繼續(xù)研究使用增廣拉格朗日方法及鄰近點算法求解稀疏約束的CT圖像恢復(fù)模型，并將模型應(yīng)用于胸部CT識別新冠肺炎病變等場景。