基于SWOMP算法的EPMA影像優(yōu)化重構(gòu)

金安安，丁雙雙，熊卿智，許婷婷，習(xí)淑萍，馬繼忠

(1.東華理工大學(xué)信息工程學(xué)院,江西 南昌 330013； 2.東華理工大學(xué)江西省核地學(xué)數(shù)據(jù)科學(xué)與系統(tǒng)工程技術(shù)研究中心,江西 南昌 330013； 3.東華理工大學(xué)軟件學(xué)院,江西 南昌 330013)

0 引 言

隨著互聯(lián)網(wǎng)技術(shù)越發(fā)成熟，5G時(shí)代已經(jīng)來臨，數(shù)據(jù)無時(shí)不在，人們面對(duì)海量的數(shù)據(jù)，使用傳統(tǒng)的香農(nóng)奈奎斯特定理(Shannon-Nyquist)[1]已經(jīng)顯得力不從心，因?yàn)樵摱ɡ砻鞔_提出，想要采樣信號(hào)不失真地精確重構(gòu)出來，必須滿足采樣速率達(dá)到信號(hào)帶寬的2倍或以上，這無疑給人們?cè)黾恿撕芏嘭?fù)擔(dān)，因?yàn)槭褂眠@種方法采集的數(shù)據(jù)有很多冗余，對(duì)人們后續(xù)的存儲(chǔ)及運(yùn)輸帶來極大的不利，特別是對(duì)電子探測(cè)顯微分析(Electron Probe Microanalysis， EPMA)[2]極為重要的電子探針影像來說。EPMA是一種微區(qū)分析，結(jié)合了顯微分析和成分分析，該分析對(duì)微小區(qū)域的化學(xué)元素分析極為有效，是一種用來分析材料內(nèi)部的元素分布狀態(tài)和組織結(jié)構(gòu)的分析方法，故對(duì)EPMA影像的存儲(chǔ)、運(yùn)輸及研究，都極為重要。而壓縮感知(Compressed Sensing， CS)[3-5]可以有效地解決這種缺陷，該方法由Tao等人于2004年提出，并在2006年由Candes[6]驗(yàn)證了其正確性，從此壓縮感知便受到學(xué)術(shù)界廣泛的關(guān)注和研究，經(jīng)過十多年的發(fā)展，其研究成果已經(jīng)越來越多，在眾多領(lǐng)域都有著廣泛的應(yīng)用，例如：醫(yī)學(xué)成像[7-9]、圖像壓縮[10]、雷達(dá)成像[11-13]、圖像重構(gòu)[14]、照相機(jī)設(shè)計(jì)[15]、激光雷達(dá)[16]、傳感器網(wǎng)絡(luò)[17]等。

CS主要由稀疏表示、測(cè)量矩陣、重構(gòu)算法[18]三者組成。這三者缺一不可，其中最重要的是重構(gòu)算法，因?yàn)檫@是壓縮感知的精髓所在，一個(gè)信號(hào)能否精準(zhǔn)并快速地恢復(fù)，都有賴于重構(gòu)算法。而經(jīng)過十多年的發(fā)展，重構(gòu)算法已經(jīng)越來越豐富，現(xiàn)在大致分成凸優(yōu)化重構(gòu)算法、貪婪重構(gòu)算法、組合重構(gòu)算法3大類。本文所研究的是分段弱正交匹配追蹤(Stagewise Weak Orthogonal Matching Pursuit， SWOMP)算法，它屬于第二類重構(gòu)算法。

SWOMP算法[19]是貪婪重構(gòu)類算法中最具代表性的，它是由Thomas Blumensath等人提出的，該算法的參數(shù)(迭代次數(shù)S和閾值alpha)是人為事先設(shè)置好的，對(duì)圖像的重構(gòu)效果并不是最佳。針對(duì)SWOMP算法重構(gòu)質(zhì)量不佳的問題，眾多學(xué)者紛紛出招，例如石曼曼等人[20]在選取原子過程中使用模糊閾值的回溯方法對(duì)錯(cuò)誤的原子進(jìn)行刪除；賀紹琪等人[21]將高斯矩陣替換為部分哈達(dá)瑪矩陣；王烈等人[22]在原子選擇過程中進(jìn)行自適應(yīng)設(shè)定閾值；王頂?shù)热薣23]將原子選擇過程中使用的內(nèi)積匹配準(zhǔn)則替換為Dice系數(shù)匹配準(zhǔn)則；侯成艷[24]同樣也是將Dice系數(shù)匹配準(zhǔn)則替換原來的內(nèi)積準(zhǔn)則進(jìn)行優(yōu)化支撐集，從而提高了重構(gòu)性能；江曉林等人[25]提出了在不同階段采取回溯方法；Zhang等人[26]采用算術(shù)閾值策略提高原子選擇準(zhǔn)確性和回溯方法刪除錯(cuò)誤原子等，這些方法都可以提高SWOMP算法的重構(gòu)質(zhì)量。本文針對(duì)這個(gè)問題，通過多次試驗(yàn)，對(duì)測(cè)量矩陣進(jìn)行修改，再修改迭代次數(shù)S和閾值alpha，找出最佳的重構(gòu)效果的參數(shù)，實(shí)驗(yàn)結(jié)果表明，本文所做出的優(yōu)化在重構(gòu)效果上有較大提升。

1 壓縮感知理論

壓縮感知理論指出：對(duì)于可以壓縮的信號(hào)(換句話來說是稀疏的)，可以使用測(cè)量矩陣對(duì)其進(jìn)行變換，然后將該變換結(jié)果映射到某個(gè)低維的空間上，但前提是該測(cè)量矩陣與變換基是不相關(guān)的，所以可以通過求解優(yōu)化問題來解決該問題，從而把原始信號(hào)重構(gòu)出來。

假設(shè)信號(hào)x是長(zhǎng)度為N的一維離散信號(hào)，并且在某個(gè)稀疏矩陣ψ上是稀疏的(即ψ包含K個(gè)非零值，且K?N)，則可表示為：

(1)

式中ψ為稀疏矩陣，ψ=[ψ1,ψ2,…,ψN]，s為稀疏系數(shù)，s=[s1,s2,…,sN]。由假設(shè)可知，該信號(hào)x是稀疏的，且稀疏度為K，所以可以使用一個(gè)測(cè)量矩陣Φ(M×N的二維矩陣，且M?N)，對(duì)信號(hào)x進(jìn)行線性測(cè)量，從而可以得到測(cè)量向量y(該向量為M×1維)，則可表示為：

y=Φx=Φψs=As

(2)

式中A=φψ，A為傳感矩陣，是一個(gè)M×N的矩陣。

如果傳感矩陣A滿足約束等距性(Restricted Isometry Property, RIP)條件[27]，即：

(3)

式中δk為常數(shù)，且δk∈(0,1)，k=1,2,…,K,K為信號(hào)x的稀疏度。當(dāng)滿足式(3)時(shí)，則s可以通過l0范數(shù)來進(jìn)行求解最優(yōu)化問題得到。

s′=arg min ‖s‖l0s.t.y=As

(4)

式中“‖?‖l0”表示0范數(shù)，也就是“?”中含有多少個(gè)非零元素。通過式(4)解出稀疏系數(shù)的逼近值s′。則原信號(hào)x=ψs′，就可以對(duì)信號(hào)x進(jìn)行準(zhǔn)確的恢復(fù)。但是式(4)是一個(gè)NP難問題，所以求解l0最小范數(shù)是非常困難的，由于l1最小范數(shù)下在一定條件下和l0最小范數(shù)具有等價(jià)性，可以得到相同的解，所以式(4)可以轉(zhuǎn)化為l1最小范數(shù)下的最優(yōu)化問題，即：

s′=arg min ‖s‖l1s.t.y=As

(5)

式(5)是一個(gè)凸最優(yōu)問題，眾所周知凸優(yōu)化問題是一種線性規(guī)劃問題，比較容易求解，所以求解公式(5)即可對(duì)原始信號(hào)x進(jìn)行重構(gòu)。

2 基于傅里葉矩陣的分段弱正交匹配追蹤算法

2.1 SWOMP算法改進(jìn)

SWOMP算法由4個(gè)參數(shù)決定，分別是：觀測(cè)向量(y)、傳感矩陣(A)、閾值(alpha)[28]、迭代次數(shù)(S)。本文對(duì)A、alpha、S這3個(gè)參數(shù)進(jìn)行修改，以尋找出最佳的參數(shù)匹配，得到最佳的重構(gòu)算法。SWOMP算法到目前為止已有不同的改進(jìn)算法[20-25]，由于不同的觀測(cè)矩陣[29]對(duì)算法的重構(gòu)質(zhì)量有不同的影響，本文在SWOMP算法的基礎(chǔ)上，將傅里葉矩陣引入其中，并對(duì)alpha進(jìn)行修正，提出一種新的算法，基于傅里葉矩陣的分段弱正交匹配追蹤(Fourier Stagewise Weak Orthogonal Matching Pursuit, FSWOMP)算法。

SWOMP算法在選擇原子的時(shí)候，這種選擇是將g[n]中最大元素的索引添加到索引集Γ[n]中，如式(6)所示：

(6)

StOMP[30]算法在選擇原子的時(shí)候引入閾值λ，選擇g[n]中最大的λ個(gè)，如式(7)所示：

(7)

(8)

式中，M為信號(hào)的維數(shù)M×N中的M，r[n]為當(dāng)前殘差。

SWOMP算法將StOMP算法的思想引入其中，選擇g[n]中最大的α個(gè)加入到索引集Γ[n]中，如式(9)所示：

(9)

由于α是事先設(shè)定的，α取值一般為0～1之間，通常情況下一般設(shè)定為α=0.5，但是在圖像重構(gòu)時(shí)，發(fā)現(xiàn)重構(gòu)效果并不是最好的，故對(duì)SWOMP算法的α進(jìn)行修改，然后更新到本算法的α中來，更新為αnew，更新過程如式(10)所示。

(10)

(11)

(12)

其中，σM(Φn)為測(cè)量矩陣Φn的第M個(gè)奇異值，N為信號(hào)的維數(shù)M×N中的N。

將式(8)中的λ代入式(12)中，得到式(13)：

(13)

再將式(13)化簡(jiǎn)得到式(14)：

(14)

因此在選擇原子的時(shí)候按照式(15)進(jìn)行更新。

(15)

2.2 FSWOMP算法流程

FSWOMP算法是在SWOMP算法的基礎(chǔ)上進(jìn)行一定的改進(jìn)。在迭代過程中，將迭代過程分成許多階段，因此在不同階段迭代的時(shí)候，對(duì)原子選擇的個(gè)數(shù)也不同，然后對(duì)索引集進(jìn)行更新，根據(jù)最后的索引集，最后對(duì)稀疏信號(hào)重構(gòu)。FSWOMP是在原來的SWOMP處引入了“弱”選擇，在對(duì)原子選擇時(shí)，對(duì)alpha進(jìn)行動(dòng)態(tài)修改，在迭代的時(shí)候，根據(jù)alpha挑選多個(gè)原子，而alpha是一個(gè)關(guān)于殘差和原子之間相關(guān)性的函數(shù)，這樣對(duì)測(cè)量矩陣沒有嚴(yán)格的要求。

FSWOMP算法的流程如圖1所示。

圖1 FSWOMP算法流程圖

Step1對(duì)各參數(shù)進(jìn)行初始化，令r0=y，Γ0=?，A0=?，t=1。其中r0代表的是第0次迭代時(shí)的殘差，t代表的是迭代次數(shù)，Γ0代表的是第0次迭代的索引集合，A0代表的是按索引Γ0所選出的A0的列集合，?代表的是空集。

Step2計(jì)算測(cè)量矩陣Φn(傅里葉正交矩陣)和殘差的內(nèi)積值，計(jì)算u=abs[ATr[t-1]]，即計(jì)算〈r[t-1],aj〉，其中1jN，選擇u中大于門限Th的值，Th=αnew·max [abs(u)]，然后將這些值對(duì)應(yīng)A的列序號(hào)為j的序列構(gòu)成集合J0。其中，abs[?]代表的是求絕對(duì)值，〈?,?〉代表的是對(duì)向量進(jìn)行求內(nèi)積，J0代表的是每次迭代所找到的索引，也就是列序號(hào)集合。

Step3對(duì)索引集和支撐集進(jìn)行更新，令A(yù)t=At-1∪aj，Γ[t]=Γ[t-1]∪J0，其中j∈J0；如果Γ[t]=Γ[t-1]，意味著新列中沒有被選到，因此立即停止迭代，從而跳轉(zhuǎn)到Step7。其中Γ[t]代表的是第t次迭代時(shí)所得到的索引集合，∪代表的是集合運(yùn)算中的并運(yùn)算，aj代表的是A的第j列構(gòu)成的集合，At所代表的是按索引Γ[t]所組成的A的列集合。

Step6令t=t+1，如果tS，則立即跳轉(zhuǎn)到Step2，再次進(jìn)行迭代，如果t>S或殘差rt=0，則終止迭代，立即跳轉(zhuǎn)到Step7。其中S代表的是迭代次數(shù)。

2.3 FSWOMP算法實(shí)驗(yàn)論證

本文對(duì)所提傅里葉矩陣的重構(gòu)率進(jìn)行驗(yàn)證，與其他常用的測(cè)量矩陣進(jìn)行對(duì)比測(cè)試，實(shí)驗(yàn)過程中使用信號(hào)長(zhǎng)度為256的高斯稀疏信號(hào)，測(cè)試間隔為5，重復(fù)進(jìn)行1000次，實(shí)驗(yàn)結(jié)果如圖2所示。

圖2 9種測(cè)量矩陣進(jìn)行信號(hào)重構(gòu)所得重構(gòu)比

由圖2可以看出，隨著K的增大，信號(hào)重構(gòu)概率越來越低，不同矩陣對(duì)SWOMP進(jìn)行重構(gòu)的影響不盡相同。傅里葉正交矩陣信號(hào)重構(gòu)率最高，無論在何稀疏度下，均能百分百對(duì)信號(hào)進(jìn)行重構(gòu)；而ToeplitzMtx的恢復(fù)概率最差，不能很好地對(duì)信號(hào)進(jìn)行恢復(fù)；PartFourierMtx的重構(gòu)質(zhì)量次之。故本文選擇傅里葉正交矩陣作為FSWOMP算法的測(cè)量矩陣是正確的。

眾所周知，迭代次數(shù)越多，則圖像重構(gòu)所需要花費(fèi)的時(shí)間越長(zhǎng)，當(dāng)然圖像所重構(gòu)的效果也越好；迭代次數(shù)越少，圖像重構(gòu)所需花費(fèi)的時(shí)間也就越短，因此圖像所重構(gòu)的質(zhì)量也就越差。但不能一味地追求質(zhì)量而舍棄時(shí)間，故圖像的好壞，需要兩者綜合來衡量，需要尋找一個(gè)合適的S，既可以使恢復(fù)的圖像質(zhì)量較好又費(fèi)時(shí)較少。故對(duì)alpha和S之間的關(guān)系進(jìn)行試驗(yàn)論證，實(shí)驗(yàn)過程中使用的信號(hào)長(zhǎng)度為256的高斯稀疏信號(hào)，而傳感矩陣為高斯隨機(jī)測(cè)量矩陣，測(cè)試間隔為5，重復(fù)進(jìn)行1000次，對(duì)重構(gòu)概率進(jìn)行繪圖顯示，實(shí)驗(yàn)結(jié)果如圖3所示。

(a) alpha=0.1

(b) alpha=0.2

(c) alpha=0.3

(d) alpha=0.4

(e) alpha=0.5

(f) alpha=0.6

(g) alpha=0.7

(h) alpha=0.8

(i) alpha=0.9

(j) alpha=1.0圖3 不同alpha下圖像重構(gòu)比

由圖3可以看出，隨著觀測(cè)數(shù)M的增大，重構(gòu)率也越來越高，在不同的alpha和S下，重構(gòu)的效果也不一樣，當(dāng)alpha=0.2和S=18時(shí)，圖像的重構(gòu)質(zhì)量較高，故綜合考慮，本文所提FSWOMP算法中alpha選擇為0.2，S選擇為18。

針對(duì)以上設(shè)置的參數(shù)進(jìn)行測(cè)試，論證其是否是最佳參數(shù)搭配。在仿真實(shí)驗(yàn)過程中，所用測(cè)試圖片BMP格式的Lena圖像，其大小為512×512，仿真實(shí)驗(yàn)結(jié)果如圖4和圖5所示。

圖4 不同閾值alpha下圖像重構(gòu)之后的圖像PSNR

此實(shí)驗(yàn)結(jié)果在測(cè)量矩陣為傅里葉矩陣和S=10時(shí)的情況下所得，實(shí)驗(yàn)結(jié)果表明，在對(duì)Lena圖像重構(gòu)時(shí)，alpha越小，圖像重構(gòu)的質(zhì)量也變得越好，相反，隨著alpha的增大，所恢復(fù)的圖像的質(zhì)量也變得越差。當(dāng)alpha=0.1和alpha=0.2時(shí)，圖像的重構(gòu)效果最好，圖像信噪比PSNR(Peak Signal to Noise Ratio)高達(dá)153.2478 dB；而alpha=1時(shí)，圖像的重構(gòu)效果相對(duì)較差，圖像信噪比PSNR為30.3176 dB。故本文將alpha選擇為0.2是正確的。

圖5 不同S下圖像重構(gòu)之后的圖像PSNR

此實(shí)驗(yàn)結(jié)果在測(cè)量矩陣為傅里葉矩陣和alpha=0.5時(shí)的情況下所得，根據(jù)圖5的實(shí)驗(yàn)結(jié)果可知，在對(duì)Lena圖像進(jìn)行恢復(fù)時(shí)，S越大，則圖像所恢復(fù)的質(zhì)量也越高，當(dāng)S=18時(shí)，圖像重構(gòu)效果最好，圖像的PSNR高達(dá)153.2478 dB，當(dāng)S大于18時(shí)，圖像的質(zhì)量沒有變化，且圖像重構(gòu)時(shí)間更長(zhǎng)，故綜合考慮，本文S選擇為18是正確的。綜上所述，本文所提FSWOMP算法中alpha選擇0.2，S選擇18。

3 EPMA影像重構(gòu)結(jié)果與分析

實(shí)驗(yàn)仿真環(huán)境為Windows 10，處理器為Intel(R) Core(TM) i5-6300HQ CPU，主頻為2.30 GHz，內(nèi)存為8 GB，64位操作系統(tǒng)，MATLAB 2016a。圖片選取的是EPMA的影像圖片，大小為1024×1024，格式為BMP。測(cè)試圖片分別是：EPMA的SEI、COMP、TOPO的3種影像，SEI、COMP、TOPO都是背散射電子所成的像，其中SEI反映表面，COMP反映成分，TOPO反映形貌。首先對(duì)EPMA影像做小波變換處理，然后獲得小波域稀疏信號(hào)X1，再利用傅里葉正交觀測(cè)矩陣，對(duì)X1進(jìn)行線性測(cè)量，用PSNR來衡量圖像重構(gòu)的效果，本文所提FSWOMP算法所得到的結(jié)果如圖6所示。

(a) Original

(b) Recovery

實(shí)驗(yàn)結(jié)果表明圖像經(jīng)過小波變換后，再利用傅里葉正交矩陣進(jìn)行觀測(cè)，所重構(gòu)圖像質(zhì)量較好，然后在alpha=0.2和迭代次數(shù)S=18的基礎(chǔ)上，進(jìn)行圖像重構(gòu)。在主觀上，紋理細(xì)節(jié)清晰可見，更為接近原始圖像，在客觀上，圖像信噪比PSNR高達(dá)130 dB，幾乎和原圖像一模一樣，甚至比原圖像效果更好，達(dá)到對(duì)EPMA影像分析的要求。

為了更加客觀地體現(xiàn)本文所重構(gòu)的圖像重構(gòu)質(zhì)量，本文對(duì)常用的幾種算法進(jìn)行仿真實(shí)驗(yàn)對(duì)比，本文選取了基追蹤算法[31](Basis Pursuit, BP)、壓縮采樣匹配追蹤[32](compressive Sampling Matching Pursuit， CoSaMP)、廣義正交匹配追蹤[33](generalized Orthogonal Matching Pursuit, gOMP)、正交匹配追蹤算法[34](Orthogonal Matching Pursuit， OMP)、正則化正交匹配追蹤[35](Regularized Orthogonal Matching Pursuit， ROMP)、子空間追蹤算法[36](Subspace Pursuit, SP)、分段正交匹配追蹤[37](StOMP)、分段弱正交匹配追蹤(SWOMP)算法等8種算法與本文所改進(jìn)的算法(FSWOMP)做一個(gè)簡(jiǎn)單的對(duì)比，對(duì)比的結(jié)果如表1所示。

表1 不同算法對(duì)3種EPMA影像重構(gòu)的PSNR 單位：dB

從表1可以清楚看出，8種不同的算法對(duì)3種測(cè)試影像進(jìn)行恢復(fù)，恢復(fù)的結(jié)果都普遍較低。其中CoSaMP算法恢復(fù)的效果最差，該算法所恢復(fù)的圖像PSNR大致穩(wěn)定在30 dB左右，而StOMP算法較好，每次重構(gòu)的圖像PSNR均在36 dB左右，至于本文研究的SWOMP算法效果一般，而本文所提出的改進(jìn)的算法，圖像重構(gòu)質(zhì)量最好，圖像的PSNR每次均高達(dá)130 dB，無論在哪種測(cè)試圖像上，重構(gòu)效果都遠(yuǎn)超其他算法，故該優(yōu)化算法在EPMA影像重構(gòu)上具有較大的提升。

圖7 不同算法對(duì)SEI影像的恢復(fù)結(jié)果

圖8 不同算法對(duì)COMP影像的恢復(fù)結(jié)果

圖9 不同算法對(duì)TOPO影像的恢復(fù)結(jié)果

從圖7～圖9中可以看出，各算法對(duì)EPMA影像進(jìn)行重構(gòu)，所得到的質(zhì)量都一般。其中CoSaMP算法和SP算法恢復(fù)的效果最差，從肉眼上都能明顯區(qū)分，因此對(duì)后續(xù)的EPMA圖像處理幾乎沒有作用，而gOMP算法和StOMP算法重構(gòu)效果較好，但相對(duì)于本文提出的優(yōu)化算法，還是有較大差距。本文所提出的算法優(yōu)化在EPMA影像上，具有很好的重構(gòu)效果，重構(gòu)的圖像極為接近原始圖像。因此本文所提出的改進(jìn)算法大大提高了圖像重構(gòu)的質(zhì)量，達(dá)到超分辨率恢復(fù)的要求。

4 結(jié)束語

本文對(duì)SWOMP算法進(jìn)行深入研究，對(duì)各種測(cè)量矩陣進(jìn)行試驗(yàn)論證，發(fā)現(xiàn)傅里葉矩陣比高斯隨機(jī)矩陣具有更高的圖像重構(gòu)質(zhì)量，能較好地達(dá)到超分辨率恢復(fù)[38-39]的要求；同時(shí)對(duì)alpha和S進(jìn)行深入研究，通過反復(fù)試驗(yàn)論證，對(duì)原有的參數(shù)進(jìn)行修改，取得了較好的圖像質(zhì)量，重構(gòu)的圖像質(zhì)量較高。針對(duì)EPMA影像，本文提出FSWOMP算法，該優(yōu)化算法相對(duì)于原算法可以精準(zhǔn)地對(duì)圖像進(jìn)行重構(gòu)，且圖像重構(gòu)的質(zhì)量遠(yuǎn)超原算法，大大提高了該算法的重構(gòu)質(zhì)量，并且達(dá)到EPMA影像超分辨率的恢復(fù)。