聯(lián)合稀疏信號(hào)恢復(fù)的貪婪增強(qiáng)貝葉斯算法

王友華，張建秋

（1.復(fù)旦大學(xué)電子工程系，上海200433；2.模擬集成電路重點(diǎn)實(shí)驗(yàn)室，重慶400060）

聯(lián)合稀疏信號(hào)恢復(fù)的貪婪增強(qiáng)貝葉斯算法

王友華1，2，張建秋1

（1.復(fù)旦大學(xué)電子工程系，上海200433；2.模擬集成電路重點(diǎn)實(shí)驗(yàn)室，重慶400060）

本文針對(duì)聯(lián)合稀疏信號(hào)恢復(fù)問(wèn)題，提出了一種貪婪增強(qiáng)貝葉斯算法.算法首先利用聯(lián)合稀疏的特點(diǎn)對(duì)信號(hào)進(jìn)行建模，然后在貝葉斯框架下，提出一種貪婪推理方式對(duì)信號(hào)恢復(fù)問(wèn)題進(jìn)行迭代求解.在迭代過(guò)程中，提出算法利用貝葉斯估計(jì)的方差信息來(lái)增強(qiáng)支撐恢復(fù)的結(jié)果，極大地提高了算法對(duì)信號(hào)恢復(fù)性能.理論分析表明：提出算法與同步正交匹配追蹤算法具有相同的計(jì)算復(fù)雜度，遠(yuǎn)低于其他聯(lián)合稀疏信號(hào)恢復(fù)算法.提出方法在具有高恢復(fù)精度和較低計(jì)算復(fù)雜度的同時(shí)，兼具貝葉斯方法和貪婪算法的優(yōu)點(diǎn).數(shù)值仿真驗(yàn)證了理論分析的有效性.

聯(lián)合稀疏；信號(hào)恢復(fù)；貪婪算法；貪婪增強(qiáng)貝葉斯算法

1 引言

由于可以用遠(yuǎn)低于香農(nóng)采樣定理所需的測(cè)量來(lái)恢復(fù)原信號(hào)，壓縮感知理論一涌現(xiàn)就得到了廣泛的關(guān)注和研究.稀疏信號(hào)恢復(fù)是壓縮感知領(lǐng)域的一個(gè)重要研究方向［1～3］.稀疏信號(hào)恢復(fù)問(wèn)題可以分為單測(cè)量和多測(cè)量壓縮感知信號(hào)恢復(fù).多測(cè)量壓縮感知信號(hào)恢復(fù)問(wèn)題也稱(chēng)為聯(lián)合稀疏信號(hào)恢復(fù)問(wèn)題.當(dāng)前，聯(lián)合稀疏信號(hào)恢復(fù)問(wèn)題逐漸成為壓縮感知研究領(lǐng)域中的一個(gè)重要研究?jī)?nèi)容［4～8］.聯(lián)合稀疏信號(hào)恢復(fù)問(wèn)題在實(shí)際中有著廣泛的應(yīng)用；如：波達(dá)角估計(jì)［7］、磁共振成像［9］、雷達(dá)成像［10］、雷達(dá)波形設(shè)計(jì)［11］等.

與單測(cè)量壓縮感知信號(hào)恢復(fù)相比，聯(lián)合稀疏信號(hào)恢復(fù)算法可恢復(fù)的信號(hào)最大稀疏度，隨著測(cè)量數(shù)的增加而增加［4，5］.目前，針對(duì)聯(lián)合稀疏信號(hào)恢復(fù)問(wèn)題，存在諸如貪婪算法、基于lp（0＜p≤1）范數(shù)最小化方法、子空間方法等幾類(lèi)算法.文獻(xiàn)［4，12］提出了一種同步正交匹配追蹤算法（S-OMP）.文獻(xiàn)［7，13］擴(kuò)展了l1范數(shù)最小化方法，提出了混合范數(shù)方法（M-BP）.文獻(xiàn)［14］將單測(cè)量壓縮感知問(wèn)題中的稀疏貝葉斯學(xué)習(xí)方法推廣到聯(lián)合稀疏信號(hào)恢復(fù)問(wèn)題中，從而提出了M-SBL方法.文獻(xiàn)［15～17］將多測(cè)量壓縮感知問(wèn)題進(jìn)行“拉長(zhǎng)”處理，將問(wèn)題轉(zhuǎn)化為單測(cè)量壓縮感知問(wèn)題進(jìn)行求解，利用“塊稀疏”的先驗(yàn)，提出了B-OMP方法.文獻(xiàn)［6］以陣列信號(hào)處理中的MUSIC方法為基礎(chǔ)，提出了CS-MUSIC方法.最近，沿著單測(cè)量信號(hào)恢復(fù)的思路，文獻(xiàn)［18］將單測(cè)量壓縮感知中的子空間追蹤方法進(jìn)行擴(kuò)展，并應(yīng)用于多測(cè)量壓縮感知中，提出了GSP方法.文獻(xiàn)［19］將單測(cè)量壓縮感知文章中的貪婪算法統(tǒng)一推廣到多測(cè)量壓縮感知問(wèn)題中，并對(duì)多測(cè)量壓縮感知問(wèn)題的相變進(jìn)行了研究.

前述方法中，貪婪算法和子空間方法具有運(yùn)算速度快的優(yōu)點(diǎn)，但是信號(hào)恢復(fù)性能不如其他方法好.基于lp（0＜p≤1）范數(shù)最小化方法和基于貝葉斯推理的算法對(duì)信號(hào)恢復(fù)性能好，但是運(yùn)算量較大.為了克服目前算法存在運(yùn)算量和恢復(fù)性能折中的問(wèn)題，本文在貝葉斯框架下，提出了一種新的聯(lián)合稀疏信號(hào)恢復(fù)算法，稱(chēng)之為貪婪增強(qiáng)貝葉斯算法（GRBA）.GRBA首先采用稀疏的先驗(yàn)知識(shí)對(duì)信號(hào)的聯(lián)合稀疏特性進(jìn)行建模，以便獲得更為準(zhǔn)確的信號(hào)估計(jì).在所建模型基礎(chǔ)上，GRBA在貝葉斯框架下，自適應(yīng)地學(xué)習(xí)稀疏模型中的超參數(shù).超參數(shù)學(xué)習(xí)過(guò)程中，GRBA利用本文提出一種貪婪的策略對(duì)信號(hào)進(jìn)行恢復(fù)，從而極大的降低了算法的復(fù)雜度.同時(shí)，GRBA充分利用貝葉斯估計(jì)不僅能夠得到信號(hào)均值信息，還能得到信號(hào)方差信息的特點(diǎn)，首次將獲得的方差信息對(duì)每次估計(jì)均值進(jìn)行增強(qiáng)，極大地提高了算法對(duì)支撐恢復(fù)的準(zhǔn)確性.理論分析和數(shù)值仿真表明：GRBA既具有傳統(tǒng)貝葉斯方法對(duì)信號(hào)恢復(fù)精度高的優(yōu)點(diǎn)，又具有貪婪算法運(yùn)算量小的特點(diǎn).特別的，GRBA相對(duì)于常見(jiàn)的貝葉斯推理算法，由于更充分利用了貝葉斯推理過(guò)程得到的均值和方差信息，極大地降低了計(jì)算量、提高了信號(hào)恢復(fù)的精度.

本文的數(shù)學(xué)符號(hào)按照如下規(guī)則進(jìn)行處理：大寫(xiě)黑斜體字母表示矩陣，小寫(xiě)斜黑體字母表示向量.矩陣X的第i行表示為Xi·，第j列表示為 X·j，第 i行、第 j列元素表示為Xij.xi表示向量x的第i個(gè)元素.向量x的2-范數(shù)表示為‖x‖2.p（x|y）表示給定y時(shí)x的條件概率密度函數(shù)；p（x；α）表示x滿(mǎn)足參數(shù)為α的概率密度函數(shù)p.N（μ；σ）表示均值為μ，均方差為σ的正態(tài)分布.diag（x）表示一對(duì)角矩陣，其對(duì)角線元素為x.矩陣X的支撐記為supp（X），對(duì)于一M×N的矩陣X，其支撐定義為：supp（X）=｛1≤i≤N：‖Xi·‖2≠0｝，|I|表示集合I的非空元素個(gè)數(shù).IL表示L×L的單位矩陣，0表示零向量，其長(zhǎng)度由上下文確定.

2 問(wèn)題描述及提出算法

2.1問(wèn)題描述

聯(lián)合稀疏信號(hào)X∈RN×L可以描述為［6～8］：

其中，W∈RM×L，T∈RM×L和A∈RM×N分別為測(cè)量噪聲矩陣、觀測(cè)數(shù)據(jù)矩陣和測(cè)量矩陣.本文假設(shè)測(cè)量噪聲為零均值高斯白噪聲.根據(jù)觀測(cè)數(shù)據(jù)矩陣T和測(cè)量矩陣A，聯(lián)合稀疏信號(hào)恢復(fù)問(wèn)題中信號(hào)的估計(jì)可表示為［6，7］：

式（1）中信號(hào)X的最大后驗(yàn)估計(jì)為［14］：

當(dāng)測(cè)量噪聲W為零均值，方差與λ相關(guān)的高斯白噪聲，且信號(hào)的先驗(yàn)分布為exp（-‖X‖0）時(shí)；式（3）和式（2）的解具有相同形式.在信號(hào)空間中，先驗(yàn)分布exp （-‖X‖0）在全零處具有最大峰值，能夠很好地描述信號(hào)的稀疏特性.但是，與該信號(hào)先驗(yàn)相對(duì)應(yīng)的式（2）卻存在大量的局部最優(yōu)解［14］.因此，研究者希望找到一種稀疏的先驗(yàn)，既能對(duì)信號(hào)的稀疏性進(jìn)行描述，又能避開(kāi)局部最優(yōu)解.而層次化的先驗(yàn)具有上述的特征.

本文設(shè)信號(hào)X的第i行Xi·（1≤i≤N）滿(mǎn)足層次化的先驗(yàn)分布［20～22］；其中，第一層為零均值高斯分布：

其中γi是信號(hào)第i行高斯分布的方差.由式（4）知，當(dāng)方差γi為零時(shí)，信號(hào)X的第i行 Xi·為零的概率是1.進(jìn)一步，當(dāng)未知的方差γi滿(mǎn)足平移截?cái)噘ゑR分布時(shí)，得到先驗(yàn)分布的第二層為［22］：

其中γ=［γ1，…，γN］T，概率分布密度函數(shù)Γτ（γi；ε，η）定義為：

式中，形狀參數(shù)ε、變化率η和平移τ等參數(shù)均大于零；式（6）等號(hào)右側(cè)分母中伽馬函數(shù)Γητ（ε）定義為：

式（4）與式（5）所定義的層次化先驗(yàn)分布推廣了高斯-伽馬分布［20］和拉普拉斯先驗(yàn)分布［21］，它們能夠?qū)π盘?hào)稀疏性進(jìn)行更好的描述.在先驗(yàn)分布的第二層次，由于未知方差參數(shù)γi（i=1，…，N）滿(mǎn)足獨(dú)立平移截?cái)嗟馁ゑR分布，有：

由式（8）可以看出：當(dāng) ε＜1時(shí)，γi以極大的概率為零；由先驗(yàn)分布的第一層次知，此時(shí)Xi·以極大的概率為零.因此，這樣的層次化先驗(yàn)分布可以很好地描述信號(hào)X的聯(lián)合稀疏特性.進(jìn)一步，設(shè)變化率參數(shù)η滿(mǎn)足先驗(yàn)分布：

易證明：當(dāng)τ，ε→0時(shí)，X的先驗(yàn)分布在零處具有無(wú)窮大概率.上述分析表明：利用上述先驗(yàn)分布得到的最大后驗(yàn)估計(jì)，能夠非常好的表示信號(hào)的稀疏解，即能夠很好的逼近式（2）.

2.2貝葉斯推理算法

由式（4）定義的信號(hào)X每行先驗(yàn)分布；可得：當(dāng)給定參數(shù)γ時(shí)，信號(hào)X的每列都為零均值的高斯分布.設(shè)信號(hào)X全部列的協(xié)方差矩陣記為Γ，則有：Γ=diag（γ1，…，γN）.設(shè)測(cè)量噪聲W的方差為σ2，由測(cè)量方程式（1）知，當(dāng)給定觀測(cè)數(shù)據(jù)矩陣T時(shí)，X的第j列X·j（1≤j≤L）后驗(yàn)分布為：

經(jīng)過(guò)計(jì)算，X第j列X·j（1≤j≤L）的后驗(yàn)分布為：

其中：

由式（11）知，該后驗(yàn)分布為高斯分布.將 X的所有列寫(xiě)成矩陣形式，得到X的后驗(yàn)分布為：

其中，均值矩陣Π和方差矩陣Σ分別為：

矩陣C定義為：

當(dāng)超參數(shù)γ已知時(shí)，均值矩陣Π可作為信號(hào)X的估計(jì)值 ^X.通過(guò)式（15），將信號(hào)X的估計(jì)問(wèn)題轉(zhuǎn)化成為對(duì)超參數(shù)γ的估計(jì)問(wèn)題.

（1）超參數(shù)估計(jì)

由于不同的超參數(shù)γ對(duì)應(yīng)不同的先驗(yàn)分布，不同的先驗(yàn)分布又決定了信號(hào)非零行的位置，從而決定測(cè)量矩陣A中對(duì)測(cè)量數(shù)據(jù)T起作用的列.因此：確定合適的超參數(shù)γ等效于根據(jù)已知測(cè)量數(shù)據(jù)矩陣T來(lái)選擇測(cè)量矩陣A中不同的列（基）.在貝葉斯框架下，可將X作為隱變量，將其積分掉［23］，從而可完成對(duì)γ的估計(jì).當(dāng)c= d=0時(shí)，將log（η）作為隱變量，可以得到下面關(guān)于γ和log（η）的對(duì)數(shù)似然函數(shù)：

其中c1為常數(shù).將式（17）代入式（18）中，可得到γ的似然函數(shù)；對(duì)其關(guān)于γi求導(dǎo)，可得：

其中：

令式（21）為零，即可求得γi.由式（20）知，f（t）的系數(shù)與方差矩Σ和均值矩陣Π相關(guān).因此，γi與方差矩陣Σ和均值矩陣Π相關(guān).而由式（15）和（16）知：方差矩陣Σ和均值矩陣Π又與超參數(shù)γ相關(guān).因此，可以采用迭代方式對(duì)γ、Σ、Π進(jìn)行求解.式（21）中一元三次方程的三次項(xiàng)系數(shù)和二次項(xiàng)系數(shù)均大于零，常數(shù)項(xiàng)小于零.由附錄中一元三次方程求根引理一可知：f（γi）=0存在一正實(shí)根，記為則似然函數(shù) L（γ，log（η））在處取得最大值.同理，可求得參數(shù)η的估計(jì)值.當(dāng)求得γ估計(jì)值后，由式（15）和式（16）可以求得矩陣Σ和Π；從而可以進(jìn)行下次迭代，直至滿(mǎn)足終止條件.

（2）噪聲方差估計(jì)

在上述分析中，假設(shè)測(cè)量噪聲W的方差σ2是已知的.當(dāng)σ2未知時(shí)，將噪聲方差σ2作為參數(shù)，利用EM算法［24］，對(duì)σ2進(jìn)行迭代求解，在每次迭代中，σ2的估計(jì)為：

3 貪婪推理方法

在前述貝葉斯推理過(guò)程中，每次迭代需對(duì)式（15）中的矩陣C求逆.當(dāng)信號(hào)維度較大時(shí)，矩陣求逆將帶來(lái)很大的運(yùn)算量.針對(duì)該問(wèn)題，本節(jié)在前述框架下，提出一種貪婪的推理方法來(lái)進(jìn)行貝葉斯估計(jì).由于信號(hào)X是聯(lián)合稀疏的，即參數(shù)γ只有較少的非零元素.由式（15）知，均值矩陣Π也只有較少非零行.這表明在計(jì)算Σ和Π時(shí)，可以只針對(duì)非零的γi進(jìn)行計(jì)算，以達(dá)到降低計(jì)算中Σ和Π維度的目的.基于這個(gè)思路，本文首先從γ= 0開(kāi)始，通過(guò)貪婪方式逐漸將非零γi加入到模型中.因此在整個(gè)計(jì)算過(guò)程中，Σ和Π的實(shí)際維度都較小.為了避免復(fù)雜的記號(hào)，設(shè)測(cè)量矩陣A按列分塊為A=［a1，a2，…，aN］.

對(duì)式（17）中C作如下變換：

其中C-i為C中不含ai和γi的部分.將式（23）其代入式（18），可得：

式中c1和c3為常數(shù).通過(guò)式（24），將對(duì)數(shù)似然函數(shù)L（γ）分成兩部分：與ai及γi不相關(guān)的部分L（γ-i），與ai及γi相關(guān)的部分l（γi）.l（γi）定義為：

其中，si，qij分別為：

由于矩陣C-i和C-1-i均為正定矩陣，信號(hào)采樣快拍數(shù)L ≥1，si＞0（i=1，…，N）.因此式（29）中系數(shù)c1，c2＞0.利用附錄所述的一元三次方程求根基本知識(shí)，下面分三種情況來(lái)討論 l（γi）在何處取得最大值.

（1）c4＜0.由引理一知h（γi）=0在（0，+∞）上有唯一的解，設(shè)為容易證明l（γi）在上單調(diào)上升，在上單調(diào)下降.因此處取得最大值，即

（2）c4≥0，c3＜0，Δ＞0.因?yàn)棣ぃ?，由引理二知，存在三個(gè)不同的實(shí)根.由于表示函數(shù)h（γi）關(guān)于γi的導(dǎo)數(shù)）且c3＜0，所以h（γi）在y軸的兩邊各有一駐點(diǎn).由于h（0）= c4≥0，所以三個(gè)不同實(shí)根包含一個(gè)負(fù)實(shí)根和兩個(gè)正實(shí)根.令為最大實(shí)根，則l（γi）從零點(diǎn)處到某點(diǎn)單調(diào)下降，然后從某點(diǎn)到l（γi）處單調(diào)上升，接著單調(diào)下降.故：l（γi）要么在零點(diǎn)處獲得最大值，要么在處獲得最大值.因此：如果，則處獲得最大值，則；如 果，則l（γi）在零處獲得最大值，則

上述分析表明：通過(guò)對(duì)一元三次方程求解和邏輯判斷可保證每次迭代后似然函數(shù)都在不斷的增加；通過(guò)迭代更新γi，以確定測(cè)量矩陣A中第i列ai是否對(duì)應(yīng)信號(hào)的非零行.在迭代求解過(guò)程中，設(shè)前次迭代的γi值記為取得最大值處對(duì)應(yīng)的自變量為如果則表明基矢量ai對(duì)增加似然函數(shù)是有貢獻(xiàn)的，應(yīng)將其納入到模型中，同時(shí)更新如果，則表明基矢量ai已存在于模型中，只需要更新γi的值為當(dāng)前迭代值即可.如果則表明基矢量ai對(duì)增加似然函數(shù)無(wú)貢獻(xiàn)，當(dāng)前的迭代應(yīng)該將其從模型中刪除，同時(shí)令完成γ的一次迭代后，更新參數(shù)η.

上述迭代過(guò)程可以從空集開(kāi)始，逐漸將基矢量ai添加到模型中，或者從模型中刪除，直到收斂.迭代過(guò)程僅需求解一元三次方程和作簡(jiǎn)單的邏輯判斷，不需要進(jìn)行高維矩陣求逆，從而極大地降低了算法的復(fù)雜度.與文獻(xiàn)［25］類(lèi)似，為了更高效地更新式（26）和式（27）中的參數(shù)si、qij，令：

其中

上式中Σactive和Aactive僅包含γi≠0對(duì)應(yīng)的ai所構(gòu)成的方差矩陣和測(cè)量矩陣.由于信號(hào)具有聯(lián)合稀疏特性，所以γ中非零元素非常少.因此，Σactive和Aactive的維度是很小.故對(duì)式（32）和式（33）進(jìn)行計(jì)算只需要較小的計(jì)算量.詳細(xì)的計(jì)算量分析見(jiàn)算法復(fù)雜度部分.

4 貝葉斯增強(qiáng)

傳統(tǒng)基于貝葉斯方法的聯(lián)合稀疏信號(hào)恢復(fù)算法［20～22］，僅僅利用了后驗(yàn)分布的均值信息，并將其作為信號(hào)的估計(jì)值.但是貝葉斯估計(jì)方法不僅能得到估計(jì)的均值信息，也能夠得到方差信息.但是，除文獻(xiàn)［26］以外，尚未有其他的文獻(xiàn)利用該信息.文獻(xiàn)［26］將方差信息作為自適應(yīng)設(shè)計(jì)測(cè)量矩陣的一種度量，并未討論如何利用方差信息來(lái)提高恢復(fù)算法本身的性能.本文將提出一種增強(qiáng)方法，它將利用方差信息來(lái)提高算法對(duì)支撐估計(jì)的正確率.由測(cè)量方程式（1）可知，X·j，T·j構(gòu)成的評(píng)價(jià)函數(shù)為：

因此，其Fisher信息矩陣為［24］：

設(shè)信號(hào)和測(cè)量噪聲不相關(guān)，有：

將式（36）至式（38）代入式（35），可得到：

不失一般性，設(shè)測(cè)量矩陣A為列歸一化的矩陣，可得：

因此：

由于Δ中元素較小，因此式（42）表示：對(duì)任給的列信號(hào)，其估計(jì)的克拉美勞下限約為噪聲功率.基于前述推導(dǎo)，算法在每次迭代過(guò)程中，根據(jù)當(dāng)前的方差矩陣和噪聲估計(jì)值進(jìn)行判斷.如果方差矩陣的某一（些）對(duì)角元素Σii，遠(yuǎn)遠(yuǎn)比噪聲方差?。患纯烧J(rèn)為，與Σii相對(duì)應(yīng)的γi值較小，其估計(jì)的可信度不高.與之對(duì)應(yīng)的基向量ai不應(yīng)在信號(hào)生成模型Aactive中.如果在迭代過(guò)程中，該基向量已經(jīng)在Aactive中；那么當(dāng)前迭代中應(yīng)該將其刪除.

5 算法復(fù)雜度分析

本文提出算法計(jì)算量主要在對(duì)式（29）至（33）的計(jì)算及做相應(yīng)的邏輯判斷上.由前述分析知：Σactive列數(shù)在每次迭代過(guò)程中動(dòng)態(tài)增加，但是最多為K列（K為X中非零行數(shù)）.因此，最壞情況下計(jì)算Σactive需要O（K2M）+O（K3）次操作.由于K＜M，Σactive計(jì)算僅需要O（K2M）次操作.在計(jì)算Qij時(shí)，需要O（MK）操作.因此，在每次迭代過(guò)程中，更新完所有的Qij時(shí)，需要O（MKNL）次操作.因此本文提出算法的復(fù)雜度與S-OMP方法的算法復(fù)雜度相當(dāng)［12］.由文獻(xiàn)［14］知，M-SBL每次迭代需要的算法復(fù)雜度為O（M2NL），用標(biāo)準(zhǔn)二階錐實(shí)現(xiàn)的M-BP算法的復(fù)雜度為O（N3L3）.由于在聯(lián)合稀疏問(wèn)題中K＜＜M ＜＜N，因此，本文提出算法的復(fù)雜度遠(yuǎn)遠(yuǎn)小于其他方法的算法復(fù)雜度.

6 仿真結(jié)果

本節(jié)在仿真中，除GRBA的性能外，還列出了S-OMP［12］、B-OMP［15］、M-BP［7，13］、M-SBL［14］及CS-MUSIC［6］的性能.在仿真中，信號(hào)相對(duì)均方誤差（RMSE）定義為：

本節(jié)通過(guò)仿真，將驗(yàn)證貝葉斯推理過(guò)程中方差矩陣對(duì)角元素分布情況.在仿真中，噪聲方差 σ2=0.01. 圖1表示了沒(méi)有增強(qiáng)時(shí)，方差矩陣Σ對(duì)角元素分布情況.圖中x-軸為信號(hào)X的行號(hào)，y-軸為估計(jì)的方差值.

圖中星號(hào)點(diǎn)為非零行對(duì)應(yīng)的方差值.從圖1中可以看出，當(dāng)無(wú)增強(qiáng)時(shí)，算法共恢復(fù)出了31個(gè)非零行；其中第183行被錯(cuò)誤地恢復(fù)成信號(hào)的支撐.從圖中可以看到正確支撐所對(duì)應(yīng)的方差值集中在σ2附近.與第183行對(duì)應(yīng)的方差明顯小于其他行所對(duì)應(yīng)的方差值.因此，可以通過(guò)前述的貝葉斯增強(qiáng)方法，提高算法對(duì)支撐正確恢復(fù)的概率.

圖2為本文提出算法有增強(qiáng)和無(wú)增強(qiáng)兩種情況下，支撐恢復(fù)率對(duì)比.在仿真中，信號(hào)非零行數(shù)從26到45均勻變化，其余仿真參數(shù)設(shè)置同上一仿真設(shè)置.

圖2中，菱形線為有方差增強(qiáng)時(shí)的支撐回復(fù)率，星號(hào)線為算法無(wú)方差修正時(shí)的結(jié)果.從圖中可知，當(dāng)算法采用方差增強(qiáng)時(shí)，可以大幅度的提高算法對(duì)支撐的恢復(fù)率.特別是在低稀疏度時(shí)，兩者具有較大的差距.在下文中，GRBA方法的仿真結(jié)果均是指具有方差增強(qiáng)時(shí)算法的性能.本小節(jié)給出了算法性能隨著信噪比變化的趨勢(shì).在仿真中，信噪比從12dB到30dB均勻變化.

圖3所示為算法的支撐恢復(fù)率與信噪比關(guān)系.從圖中可以看到，當(dāng)信噪比大于一定“閾值”后，GRBA方法和M-SBL方法對(duì)支撐恢復(fù)性能迅速提高.對(duì)于GRBA方法，該“閾值”約為13dB；對(duì)于M-SBL方法，該“閾值”約為25dB.因此，GRBA對(duì)測(cè)量噪聲敏感性遠(yuǎn)低于MSBL方法.同時(shí)，由圖3知：CS-MUSIC方法的支撐恢復(fù)率隨著信噪比的增加而線性增加；但是增加的速度較慢.兩種貪婪的恢復(fù)算法：S-OMP和B-OMP方法在仿真的信噪比區(qū)間均不能正確的恢復(fù)信號(hào)支撐.

圖4所示為信號(hào)的均方誤差與信噪比關(guān)系.從圖中可以看到，基于統(tǒng)計(jì)推理的算法對(duì)信號(hào)幅度恢復(fù)明顯優(yōu)于確定性算法.在所有算法中，GRBA方法幅度恢復(fù)性能最優(yōu)，CS-MUSIC方法對(duì)幅度恢復(fù)性能最差.在低信噪比區(qū)域，由于GRBA和M-SBL對(duì)支撐恢復(fù)率均較低，因此兩者間的差距較小.隨著信噪比的增加，GRBA對(duì)支撐成功恢復(fù)率迅速增加，而M-SBL對(duì)支撐恢復(fù)率仍然較低；因此，GRBA和M-SBL之間的差距拉大.隨著信噪比的進(jìn)一步增加，GRBA和M-SBL均能夠以極高的概率恢復(fù)信號(hào)支撐，此時(shí)GRBA和 M-SBL之間的差距縮寫(xiě).從圖中可以看到，兩種貪婪的恢復(fù)算法：S-OMP方法和B-OMP方法對(duì)信號(hào)幅度恢復(fù)差別較小.

圖5所示為算法運(yùn)行時(shí)間隨著信噪比變化曲線.圖5上圖為所有算法運(yùn)行時(shí)間的對(duì)比圖.從上圖可以看到，M-BP算法的運(yùn)算時(shí)間最長(zhǎng)，遠(yuǎn)遠(yuǎn)高于其他的算法運(yùn)行時(shí)間；其次為M-SBL方法，算法S-OMP方法的運(yùn)行時(shí)間最短.圖5下圖顯示了GRBA方法、B-OMP方法、SOMP方法和CS-MUSIC方法的運(yùn)行時(shí)間.從圖中可以看到，雖然本文提出的GRBA方法為基于貝葉斯推理的恢復(fù)方法；但是采用了貪婪推理方法，運(yùn)算時(shí)間不僅僅遠(yuǎn)低于同樣為統(tǒng)計(jì)推理的M-SBL方法；還低于貪婪算法B-OMP的運(yùn)行時(shí)間，和貪婪方法S-OMP的運(yùn)行時(shí)間相似.且隨著信噪比的增加，GRBA運(yùn)行時(shí)間和S-OMP方法運(yùn)行時(shí)間的差距越來(lái)越小.從圖中看到，隨著信噪比的減少，噪聲功率增加；M-SBL方法迭代次數(shù)增加，運(yùn)行時(shí)間變長(zhǎng).而GRBA方法采用估計(jì)方差信息增強(qiáng)，將估計(jì)方差和估計(jì)噪聲功率進(jìn)行比較，能自適應(yīng)地估計(jì)信號(hào)的支撐.因此，與M-SBL方法相比，其運(yùn)行時(shí)間隨著信噪比變化波動(dòng)較小.另外，從圖5中可以看到，同為貪婪算法，B-OMP方法運(yùn)行時(shí)間較S-OMP方法運(yùn)行時(shí)間長(zhǎng).造成這種情況是因?yàn)椋築-OMP在求解過(guò)程中，將多測(cè)量壓縮感知問(wèn)題轉(zhuǎn)化成單測(cè)量壓縮感知問(wèn)題，將信號(hào)X和測(cè)量數(shù)據(jù)T進(jìn)行“拉長(zhǎng)”處理，極大的增加了求解問(wèn)題的維度，增加了算法運(yùn)行時(shí)間.

當(dāng)改變信號(hào)快拍數(shù)、信號(hào)稀疏度時(shí)，算法取得了與前述類(lèi)似的性能，由于篇幅限制，本文沒(méi)有將其性能一一列出.

7 結(jié)論

本文針對(duì)聯(lián)合稀疏信號(hào)恢復(fù)問(wèn)題，提出了一種增強(qiáng)貝葉斯信號(hào)恢復(fù)算法.算法在貝葉斯框架下，采用了一種貪婪的方法和一種增強(qiáng)的策略對(duì)信號(hào)稀疏模型進(jìn)行求解.從而使本文提出方法既具有貝葉斯方法恢復(fù)精度高的優(yōu)點(diǎn)又具有貪婪算法運(yùn)算量小的特點(diǎn).理論分析和數(shù)值仿真結(jié)果證明了提出算法的良好特性.

附錄

引理二［28］對(duì)一元三次方程令其判別式為：

如果Δ＞0，則g（t）=0有三個(gè)不同的實(shí)根；如果Δ =0，則方程有三個(gè)實(shí)根，其中包含重根；如果Δ＜0，則方程包含一實(shí)根和一對(duì)共軛復(fù)根.

［1］D L Donoho.Compressed sensing［J］.IEEE Trans on Information Theory，2006，52（4）：1289-1306.

［2］E Candes，T Tao.Decoding by linear programming［J］. IEEE Trans on Information Theory，2005，51（12）：4203 -4215.

［3］焦李成，楊淑媛，劉芳，侯彪.壓縮感知回顧與展望［J］.電子學(xué)報(bào)，2011，20（7）：1651-1662. JIAO Li-cheng，YANG Shu-yuan，LIU Fang，HOU Biao. Development and prospect of compressive sensing［J］.Acta Electronica Sinica，2011，20（7）：1651-1662.（in Chinese）

［4］J Chen，X Huo.Theoretical results on sparse representations of multiple measurement vectors［J］.IEEE Trans on Signal Processing，2006，54（12）：4634-4643.

［5］P Feng.Universal Minimum-Rate Sampling and Spectrum-Blind Reconstruction for Multiband Signals［D］.Champaign：University of Illinois，1997.

［6］J M Kim，et al.Compressive MUSIC：A missing link between compressive sensing and array signal processing［J］.IEEE Trans on Information Theory，2012，58（1）：278 -301.

［7］D Malioutov，et al.A sparse signal reconstruction perspective for source localization with sensor arrays［J］.IEEE Trans on Signal Processing，2005，53（8）：3010-3022.

［8］M F Duarte，et al.Distributed compressed sensing of jointly sparse signals［A］.Asilomar Conference on Signals，Systems and Computers［C］.USA：Asilomar，2005.1537 -1541.

［9］O Lee，et al.Compressive diffuse optical tomography：noniterative exact reconstruction using joint sparsity［J］.IEEE Trans on Medical Imaging，2011，30（5）：1129-1142.

［10］周漢飛，李禹，粟毅.基于壓縮感知的多角度SAR特征提?。跩］.電子學(xué)報(bào)，2013，41（3）：543-548. ZHOU Han-fei，LI Yu，SU Yi.Multi-aspect SAR feature extraction based on compressive sensing［J］.Acta Electronica Sinica，2013，41（3）：543-548.（in Chinese）

［11］賀亞鵬，朱曉華，等.噪聲干擾背景下壓縮感知雷達(dá)波形優(yōu)化設(shè)計(jì)［J］.電子學(xué)報(bào)，2014，42（3）：469-476. HE Ya-peng，ZHU Xiao-hua，et al.Waveform design for compressive sensing radar in the presence of interference and noise［J］.Acta Electronica Sinica，2014，42（3）：469-476.（in Chinese）

［12］J A Tropp，A C Gilbert，M J Strauss.Algorithms for simultaneous sparse approximation，Part I：Greedy pursuit［J］. Signal Processing，2006，86（3）：572-588.

［13］J A Tropp.Algorithms for simultaneous sparse approximation，Part II：Convex relaxation［J］.Signal Processing，2006，86（3）：589-602.

［14］D P Wipf.Bayesian Methods for Finding Sparse Representations［D］.San Diego：University of California，2006.

［15］Y C Eldar，P Kuppinger，et al.Block-sparse signals：uncertainty relations and efficient recovery［J］.IEEE Transa on Signal Processing，2010，58（6）：3042-3054.

［16］R G Baraniuk，V Cevher，et al.Model-based compressive sensing［J］.IEEE Trans on Information Theory，2010，56 （4）：1982-2001.

［17］Y C Eldar，M Mishali.Robust recovery of signals from a structured union of subspaces［J］.IEEE Trans on Information Theory，2009，55（11）：5302-5316.

［18］J M Feng.Generalized subspace pursuit for signal recovery from multiple-measurement vectors［A］.Proceedings of IEEE Wireless Communications and Networking Conference［C］.Shanghai：IEEE，2013.2874-2878.

［19］J D Blanchard，M Cermak，et al.Greedy algorithms for joint sparse recovery［J］.IEEE Trans on Signal Processing，2014，62（7）：1694-1704.

［20］N Pedersen，D Shutin，et al.Sparse Estimation Using Bayesian Hierarchical Prior Modeling for Real and ComplexModels［OL］.http：//arxiv.org/pdf/1108. 4324v2，2011.

［21］S Babacan，R Molina，et al.Bayesian compressive sensing using Laplace priors［J］.IEEE Trans on Image Processing，2010，19（1）：53-63.

［22］Z Yang，L H Xie，C S Zhang.Bayesian Compressed Sensing With New Sparsity-Inducing Prior［OL］.http：//arxiv.org/pdf/1208.6464v1，2012.

［23］D J C MacKay.Bayesian interpolation［J］.Neural Computation，1992，4（3）：415-447.

［24］M R Gupta，Y H Chen.Theory and use of the EM algorithm［J］.Foundations and Trends in Signal Processing，2011，4（3）：223-296.

［25］M Tipping，A Faul.Fast marginal likelihood maximisation for sparse Bayesian models［A］.Proceedings of the 9th International Workshop Artificial Intelligence and Statistics ［C］.Florida：Key West，2003.1-13.

［26］S Ji，Y Xue，L Carin.Bayesian compressive sensing［J］. IEEE Trans on Signal Processing，2008，56（6）：2346 -2356.

［27］Y C Eldar，G Kutynoik.Compressed Sensing［M］.Cambridge：Cambridge Press，2012.

［28］R Irving.Integers，Polynomials，and Rings：A Course in Algebra［M］.Berlin：Springer Verlag，2004.

王友華 男，1981年生，重慶人.2011年進(jìn)入復(fù)旦大學(xué)電子工程系，現(xiàn)為博士生，從事壓縮感知信號(hào)恢復(fù)，混合信號(hào)集成電路設(shè)計(jì)等工作.

E-mail：wangyouhua@fudan.edu.cn

張建秋 男，1962年生于湖南隆回縣，現(xiàn)任復(fù)旦大學(xué)電子工程系教授、博士生導(dǎo)師、IEEE高級(jí)會(huì)員，主要研究領(lǐng)域有信息處理理論及其在測(cè)量和儀器、新型傳感器、控制和通信中的應(yīng)用.

E-mail：jqzhang01@fudan.edu.cn

A Greedy Refinement Bayesian Approach to Joint Sparse Signal Recovery

WANG You-hua1，2，ZHANG Jian-qiu1

（1.Department of Electronic Engineering，F(xiàn)udan University，Shanghai 200433，China；2.Science and Technology on Analog Integrated Circuits Laboratory，Chongqing 400060，China）

In this paper，a new greedy refinement bayesian approach（GRBA），used to solve the joint sparse signal recovery problem，is proposed.The joint sparse property of signals is first used to model the signals.Based on the model，a greedy Bayesian inference method used to estimate the signals is then presented.In order to enhance the performance of the recovery，the covariance matrix got by the Bayesian inference is utilized to refine the support recovery results in our inference process.The analytical results show that GRBA outperforms the reported algorithms in the literature in terms of both the signal recovery accuracy and computational complexity.It keeps both the advantages of Bayesian methods and greedy methods. Numerical simulations verify the effectiveness of the analytical results.

joint sparsity；signal recovery；greedy algorithm；greedy refinement Bayesian approach

TN919.72

A

0372-2112（2016）04-0780-08

電子學(xué)報(bào)URL：http：//www.ejournal.org.cn 10.3969/j.issn.0372-2112.2016.04.005

2014-08-18；

2014-12-14；責(zé)任編輯：孫瑤

國(guó)家自然科學(xué)基金（No.61171127，No.61571131）；模擬集成電路重點(diǎn)實(shí)驗(yàn)室基金（No.9140C090110130C09003）