張譽(yù)馨,張建秋

(復(fù)旦大學(xué) 電子工程系 智慧網(wǎng)絡(luò)與系統(tǒng)研究中心，上海 200433)

線譜估計(jì)在雷達(dá)、聲吶、地震學(xué)等領(lǐng)域有著廣泛的應(yīng)用,因此它一直是信號(hào)處理領(lǐng)域中研究的重要問(wèn)題[1].經(jīng)典譜估計(jì)算法包括周期圖法、非線性最小二乘法(Nonlinear Least Squares, NLS)及多重信號(hào)分類法(Multiple Signal Classification, MUSIC)等算法.但這些算法都有各自的限制和缺陷[2].

近年來(lái),隨著稀疏表示和壓縮感知理論[3]的發(fā)展,稀疏譜估計(jì)算法[4-5]得到了廣泛關(guān)注.稀疏譜估計(jì)算法將連續(xù)的頻譜離散化為均勻分布的頻率格點(diǎn),并假設(shè)待估計(jì)譜信號(hào)的真實(shí)頻率分量均落在預(yù)先離散化的格點(diǎn)上.這樣,對(duì)譜信號(hào)的頻率和幅度的估計(jì),就轉(zhuǎn)化成了求解稀疏支撐向量及其對(duì)應(yīng)的系數(shù).稀疏恢復(fù)算法通常假設(shè)觀測(cè)噪聲是高斯的.但是,在一些應(yīng)用場(chǎng)景中,由突發(fā)干擾、儀器突然失常或測(cè)量者操作不當(dāng)?shù)纫蛩囟鴮?dǎo)致的觀測(cè)異常值,使得觀測(cè)噪聲不再能用高斯噪聲作為其模型[6].異常值的存在,往往意味著觀測(cè)數(shù)據(jù)具有重尾分布[7].文獻(xiàn)[7]的研究表明： 重尾數(shù)據(jù)與高斯數(shù)據(jù)相比,其產(chǎn)生與均值差異較大的觀測(cè)值可能性較大.當(dāng)觀測(cè)存在異常值時(shí),傳統(tǒng)稀疏恢復(fù)算法性能將下降甚至失效[8].文獻(xiàn)[8]研究表明,使用lp(0

可是,待估計(jì)的頻率未必落在預(yù)先離散化的格點(diǎn)上,此時(shí)稀疏譜估計(jì)的性能會(huì)下降甚至失效,文獻(xiàn)[10]稱此問(wèn)題為格點(diǎn)失配問(wèn)題.近年來(lái),研究者針對(duì)格點(diǎn)失配問(wèn)題展開(kāi)了廣泛的研究,提出一系列無(wú)格點(diǎn)稀疏譜估計(jì)算法.此類算法可直接求解連續(xù)頻域的譜,因而有效避免了格點(diǎn)失配問(wèn)題.原子范數(shù)最小化(Atomic Norm Minimization, ANM)[10]和基于協(xié)方差的無(wú)格點(diǎn)稀疏迭代估計(jì)(Gridless SPICE, GLS)算法[12]就是它們的典型代表.然而,ANM算法對(duì)觀測(cè)噪聲中可能存在的異常值是不魯棒的.GLS算法雖然對(duì)沖激噪聲魯棒,但是存在頻率分裂問(wèn)題[12].頻率分裂意味著對(duì)一個(gè)幅度較大的待估計(jì)線頻,算法將在其兩側(cè)估計(jì)得到兩個(gè)幅度較小的線頻.為了解決GLS算法的頻率分裂問(wèn)題,文獻(xiàn)[12]進(jìn)一步提出了一個(gè)譜估計(jì)方法,該方法首先使用GLS算法,來(lái)估計(jì)待估計(jì)信號(hào)的協(xié)方差矩陣,然后采用統(tǒng)計(jì)方法[13]來(lái)估計(jì)模型的階,最后再使用MUSIC算法來(lái)進(jìn)行頻率估計(jì).盡管上述無(wú)格點(diǎn)稀疏算法完全消除了頻域離散化方法存在的問(wèn)題,然而它們始終需要待估計(jì)線頻之間具有足夠的間隔才能保證精確恢復(fù).文獻(xiàn)[14]的超分辨率譜估計(jì)的理論分析表明： 無(wú)格點(diǎn)稀疏譜估計(jì)算法,只對(duì)線頻間隔大于2倍采樣數(shù)據(jù)窗主瓣寬度(以后將簡(jiǎn)稱主瓣寬度)以上信號(hào)才能精確恢復(fù).文獻(xiàn)[15]的研究則表明： 上述線頻間隔僅為充分條件,實(shí)際的限制則與主瓣寬度和模型階有關(guān).

為了提高線頻間隔位于c倍主瓣寬度內(nèi)的譜估計(jì)精度,本文提出了一種對(duì)異常值魯棒的無(wú)格點(diǎn)譜估計(jì)(Outlier Robust and Grid Free Spectral Estimation, ORGFSE)算法.該算法首先將GLS算法的頻率估計(jì)誤差描述成誤差向量,再利用l1范數(shù)對(duì)異常值的魯棒性,來(lái)分別約束信號(hào)幅度和誤差向量的擬合誤差,進(jìn)而通過(guò)交替迭代來(lái)使它們同時(shí)到達(dá)最小,以便同時(shí)對(duì)信號(hào)幅度和誤差向量的聯(lián)合魯棒進(jìn)行估計(jì).分析表明： 提出的算法在保證收斂的同時(shí),還能提高在c倍主瓣寬度內(nèi)的譜估計(jì)精度.仿真實(shí)驗(yàn)結(jié)果則在驗(yàn)證分析結(jié)果有效性的同時(shí),也表明當(dāng)信號(hào)的實(shí)際頻率間隔大于1倍主瓣寬度時(shí),大多數(shù)情況下,由GLS得出的結(jié)果是準(zhǔn)確的,只有在極少數(shù)情況下(概率通常小于3%),GLS算法才可能出現(xiàn)頻率分裂.因此,對(duì)于GLS算法估計(jì)出的信號(hào),當(dāng)其線頻間距在1倍至c倍主瓣寬度之間時(shí),可利用本文的聯(lián)合魯棒估計(jì)來(lái)進(jìn)行線譜估計(jì)；而一旦發(fā)現(xiàn)其間距小于1倍主瓣寬度時(shí),則采用MUSIC算法來(lái)進(jìn)行線譜估計(jì).這樣本文提出的算法就能在提高譜估計(jì)精度和對(duì)重尾分布噪聲魯棒的同時(shí),減小算法的運(yùn)算量.

本文中的數(shù)學(xué)符號(hào)遵照如下規(guī)定： 大寫斜黑體字母表示矩陣,小寫斜黑體字母表示向量.矩陣diag(x)表示對(duì)角線元素為x的對(duì)角線矩陣.矩陣T(x)表示向量x的Toeplitz陣.tr(R)表示矩陣R的跡.xi表示向量x的第i個(gè)元素.‖x‖p(p≥1)表示向量x的p范數(shù).表示取模.

1 信號(hào)模型

在線譜估計(jì)問(wèn)題中,觀測(cè)信號(hào)是含噪聲的K個(gè)復(fù)正弦信號(hào)的線性組合[1],它可描述為：

(1)

定義頻率導(dǎo)向矢量a(w)=[e-jw,e-j2w,…,e-jMw]T,則模型式(1)可以重寫為如下矩陣形式：

(2)

(3)

(4)

文獻(xiàn)[14]和[15]研究表明： 盡管GLS算法不存在格點(diǎn)失配問(wèn)題,但當(dāng)實(shí)際信號(hào)頻率位于c(1

(5)

(6)

結(jié)合式(2)和式(6)可得一個(gè)新的觀測(cè)模型如下：

(7)

2 算法描述

針對(duì)GLS頻率估計(jì)誤差和觀測(cè)異常值同時(shí)存在的稀疏信號(hào)恢復(fù)問(wèn)題,本文將優(yōu)化問(wèn)題的目標(biāo)函數(shù)描述成如下形式：

(8)

由于式(8)所示的目標(biāo)函數(shù)存在兩個(gè)未知量s和δ,難以直接求解其全局最優(yōu)解.可是,當(dāng)求解待恢復(fù)信號(hào)的幅度矢量s時(shí),如果假設(shè)誤差矢量δ已知,此時(shí)式(8)可轉(zhuǎn)化為：

(9)

同理,當(dāng)求解誤差矢量δ時(shí),如果假設(shè)待恢復(fù)信號(hào)的幅度矢量s已知,那么式(8)可重寫為：

(10)

(11)

(12)

本文稱在利用GLS算法估計(jì)值的基礎(chǔ)上,再由式(11)和(12)構(gòu)成的交替迭代譜估計(jì)算法,為對(duì)異常值魯棒的無(wú)格點(diǎn)譜估計(jì)(ORGFSE)算法.由文獻(xiàn)[16]的分析可知,本文提出算法ORGFSE的運(yùn)算復(fù)雜度為O(KM2),而MUSIC算法[17]的運(yùn)算復(fù)雜度為O(M3),由于K?M,因而本文提出算法的運(yùn)算復(fù)雜度更低.

此外,當(dāng)信號(hào)的實(shí)際頻率間隔大于1倍主瓣寬度時(shí),大多數(shù)情況下,由GLS得出的結(jié)果是準(zhǔn)確的,只有在極少數(shù)情況下(概率通常小于3%),GLS算法才可能出現(xiàn)頻率分裂.因此,為了減小算法的運(yùn)算量,本文在GLS譜估計(jì)方法的基礎(chǔ)上提出了一種對(duì)異常值魯棒新穎的無(wú)格點(diǎn)譜估計(jì)算法,其算法流程總結(jié)如下：

1) 采用GLS算法估計(jì)信號(hào)的頻率及協(xié)方差矩陣；

2) 采用統(tǒng)計(jì)方法估計(jì)模型階K,即線譜的個(gè)數(shù)；

3) 對(duì)于GLS算法估計(jì)出的K個(gè)線頻,若其中任意兩個(gè)之間的線頻間距在c(1

3 收斂性分析

本文提出算法在文獻(xiàn)[18]中稱為循環(huán)最小算法.由于文獻(xiàn)[18]證明了循環(huán)最小算法與最大最小算法等價(jià),因此只需證明本文提出算法中,所構(gòu)造目標(biāo)函數(shù)的替代函數(shù)滿足文獻(xiàn)[18]中定理1所示的上界特性,便可以由文獻(xiàn)[18]的結(jié)果,直接斷定本文提出算法是收斂的.

定理1[18]上界特性： 在點(diǎn)xt處,最大最小算法構(gòu)造了目標(biāo)函數(shù)f(x)的一個(gè)連續(xù)的替代函數(shù)g(·|xt).此替代函數(shù)滿足如下的上界特性：

g(x|xt)≥f(x)+ct, ?x,

其中ct=g(xt|xt)-f(xt),進(jìn)而x可以通過(guò)迭代算法更新如下：

在本文算法中,可以將問(wèn)題式(7)描述為：

本文目的是通過(guò)交替迭代算法,來(lái)使關(guān)于s和δ的函數(shù)同時(shí)達(dá)到最小,以實(shí)現(xiàn)了目標(biāo)函數(shù)f(s)的最小化.也就是說(shuō),(s,δ)可以通過(guò)迭代算法更新如下：

現(xiàn)在,證明本文構(gòu)造目標(biāo)函數(shù)f(s)的替代函數(shù)滿足上界特性.

g(s,δ*(s(k)))≥g(s,δ*(s))=f(s),

其中c(k)=g(s(k)|s(k))-f(s(k))=g(s(k),δ*(s(k)))-f(s(k))=0.此外,由于

4 仿真實(shí)驗(yàn)與分析

4.1 實(shí)驗(yàn)條件

在仿真中,采樣點(diǎn)數(shù)M為50,信號(hào)包含實(shí)際線頻數(shù)K為3,IAA,SPICE算法的頻域離散化格點(diǎn)數(shù)為500,迭代次數(shù)為15次.定義ORGFSE算法中的第k次迭代的殘差為：

r(k)=y-(A+Bdiag(δ(k)))s(k).

算法中含異常值的觀測(cè)噪聲e分別使用高斯混合模型(Gaussian Mixture Model, GMM)噪聲和α-穩(wěn)定噪聲來(lái)作為其模型.

GMM噪聲的概率密度函數(shù)為：

其中SNR為信噪比.

對(duì)稱α-穩(wěn)定噪聲的特征函數(shù)為

φ(en)=exp(γα|en|α),

其中參數(shù)α稱為特征參數(shù),α∈(0,2],仿真中取α=1.4；分散系數(shù)γ由SNR和信號(hào)功率決定：

4.2 仿真實(shí)驗(yàn)結(jié)果

本文進(jìn)行了6組仿真實(shí)驗(yàn).第1,3,5組分別研究算法恢復(fù)精度與GMM、α-穩(wěn)定噪聲和高斯噪聲及其信噪比的關(guān)系.實(shí)驗(yàn)中信號(hào)的3個(gè)線頻頻率分別在各自的區(qū)間內(nèi)隨機(jī)分布,且f1∈[0.102,0.104],f2=f1+Δf,f3∈[0.499,0.501],其中Δf=c/M.3個(gè)線頻的幅度分別為2,2,1.第2,4,6組實(shí)驗(yàn)分別研究了3種噪聲背景下算法恢復(fù)精度與頻率間距Δf的關(guān)系,實(shí)驗(yàn)中3個(gè)線頻的幅度在區(qū)間[1,2]內(nèi)隨機(jī)分布.

第1組仿真實(shí)驗(yàn),研究了頻率估計(jì)誤差與α-穩(wěn)定噪聲和信噪比的關(guān)系.α-穩(wěn)定噪聲的信噪比范圍為0dB 到30dB.仿真中取c=1.1.圖1為信號(hào)的頻率估計(jì)誤差f_MSE隨著α-穩(wěn)定噪聲信噪比增加而變化的仿真結(jié)果.仿真實(shí)驗(yàn)表明,隨著α-穩(wěn)定噪聲信噪比的增加,所有算法的頻率估計(jì)誤差呈現(xiàn)減小的趨勢(shì).當(dāng)信噪比小于8dB時(shí),所有算法均存在較大的頻率估計(jì)誤差.當(dāng)信噪比大于8dB時(shí),本文提出算法的頻率估計(jì)誤差最小.此外,SPICE算法將頻域離散化,存在格點(diǎn)失配問(wèn)題.當(dāng)信噪比為10dB以上時(shí),SPICE的估計(jì)誤差基本不隨信噪比增大而變化,此時(shí)的誤差主要是格點(diǎn)失配誤差.由于以上算法均采用l1范數(shù)約束擬合誤差,因而對(duì)α-穩(wěn)定噪聲魯棒IAA, ANM算法的估計(jì)性能較差,原因是這兩種算法對(duì)α-穩(wěn)定噪聲非常敏感.實(shí)驗(yàn)表明,由于采用l1范數(shù)約束擬合誤差,本文提出算法對(duì)α-穩(wěn)定噪聲具有魯棒性,而且在所有算法中恢復(fù)精度最高.當(dāng)f_MSE為-65dB時(shí),與GLS-MUSIC和SPICE-MUSIC算法相比,本文算法獲得了5dB的增益.

圖1 頻率估計(jì)誤差與α-穩(wěn)定噪聲和信噪比的關(guān)系Fig.1 Frequency estimation error analysis with α-stable noise SNR varying

圖2 α-穩(wěn)定噪聲下頻率估計(jì)誤差與頻率間距的關(guān)系Fig.2 Frequency estimation error analysis with frequency intervals varying under α-stable noise

第2組仿真實(shí)驗(yàn)研究了信噪比為20dB的α-穩(wěn)定噪聲背景下頻率估計(jì)誤差與頻率間距的關(guān)系.圖2為信號(hào)的頻率估計(jì)誤差f_MSE隨著頻率間距Δf=c/M的系數(shù)c增加而變化的仿真結(jié)果.仿真實(shí)驗(yàn)表明,隨著c的增加,對(duì)異常值魯棒的所有算法的頻率估計(jì)誤差均呈現(xiàn)減小的趨勢(shì).當(dāng)c小于1.4時(shí),本文提出算法的頻率估計(jì)誤差最小.當(dāng)c大于1.4時(shí),這3種算法的頻率估計(jì)誤差基本相同.此外,在c大于1.2時(shí),SPICE算法的估計(jì)誤差基本不隨信噪比增大而變化,此時(shí)的誤差主要是格點(diǎn)失配誤差；由于IAA和ANM算法對(duì)α-穩(wěn)定噪聲非常敏感,因而估計(jì)性能較差.仿真實(shí)驗(yàn)表明,在α-穩(wěn)定噪聲背景下,當(dāng)信號(hào)存在頻率間距c倍主瓣寬度以內(nèi)的線頻時(shí),本文提出算法在所有算法中恢復(fù)精度最高.當(dāng)f_MSE為-65dB時(shí),與GLS-MUSIC和SPICE-MUSIC算法相比,本文算法可分辨的頻率間距系數(shù)c縮小了0.06.

第5組仿真實(shí)驗(yàn),研究了頻率估計(jì)誤差與高斯噪聲和信噪比的關(guān)系.高斯噪聲信噪比的范圍為-5dB到25dB.仿真中取c=1.1.圖5(看98頁(yè))為信號(hào)的頻率估計(jì)誤差f_MSE隨著高斯噪聲信噪比增加而變化的仿真結(jié)果.仿真實(shí)驗(yàn)表明,隨著高斯噪聲信噪比的增加,所有算法的頻率估計(jì)誤差呈現(xiàn)減小的趨勢(shì).當(dāng)信噪比小于0dB時(shí),所有算法均存在較大的頻率估計(jì)誤差.當(dāng)信噪比大于0dB時(shí),本文提出算法的頻率估計(jì)誤差最小.當(dāng)信噪比大于2dB時(shí),IAA和SPICE算法的估計(jì)誤差基本不隨信噪比增大而變化,此時(shí)的誤差主要是格點(diǎn)失配誤差.本文仿真實(shí)驗(yàn)表明,即使在高斯噪聲背景下,本文提出算法也是恢復(fù)精度最高.當(dāng)f_MSE為-65dB時(shí),與ANM,GLS-MUSIC和SPICE-MUSIC算法相比,本文算法獲得了6dB的增益.

圖3 頻率估計(jì)誤差與GMM噪聲及信噪比的關(guān)系Fig.3 Frequency estimation error analysis with GMM noise SNR varying

圖4 GMM噪聲下頻率估計(jì)誤差與頻率間距的關(guān)系Fig.4 Frequency estimation error analysis with frequency intervals varying under GMM noise

第6組仿真實(shí)驗(yàn)研究了信噪比為15dB的高斯噪聲背景下頻率估計(jì)誤差與頻率間距的關(guān)系.圖6為信號(hào)的頻率估計(jì)誤差f_MSE隨著頻率間距Δf=c/M的系數(shù)c增加而變化的仿真結(jié)果.仿真實(shí)驗(yàn)表明,隨著c的增加,所有算法的頻率估計(jì)誤差呈現(xiàn)減小的趨勢(shì).當(dāng)c小于1.4時(shí),本文算法的恢復(fù)精度在所有基于l1范數(shù)的算法中仍是最高.當(dāng)c大于1.4時(shí),這3種算法的頻率估計(jì)誤差基本相同.此外,當(dāng)c大于1.15時(shí),SPICE和IAA算法的估計(jì)誤差基本不隨信噪比增大而變化,此時(shí)的誤差主要是格點(diǎn)失配誤差.仿真實(shí)驗(yàn)表明,在高斯噪聲背景下,當(dāng)信號(hào)存在頻率間距c倍主瓣寬度以內(nèi)的線頻時(shí),本文提出算法在所有算法中恢復(fù)精度最高.當(dāng)f_MSE為-60dB時(shí),與GLS-MUSIC和SPICE-MUSIC算法相比,本文算法可分辨的頻率間距系數(shù)c縮小了0.08.

圖5 頻率估計(jì)誤差與高斯噪聲及信噪比的關(guān)系Fig.5 Frequency estimation error analysis with Gaussian noise SNR varying

圖6 高斯噪聲背景下頻率估計(jì)誤差與頻率間距的關(guān)系Fig.6 Frequency estimation error analysis with frequency intervals varying under Gaussian noise

5 結(jié) 論

本文針對(duì)稀疏譜估計(jì)中,GLS算法對(duì)頻率間隔位于c(1

參考文獻(xiàn)：

[1] LI H. Spectral analysis of signals [J].IEEESignalProcessingMagazine, 2007,24(1)： 148-150.

[2] YARDIBI T, LI J, STOICA P, et al. Source localization and sensing： Anonparametric iterative adaptive approach based on weighted least squares [J].IEEETAeroElecSys, 2010,46(1)： 425-443.

[3] WANG Y H, ZHANG J Q. Robust signal recovery algorithm for structured perturbation compressive sensing [J].JournalofSystemsEngineeringandElectronics, 2016,27(02)： 319-325.

[4] DUARTE M F, BARANIU KR G. Spectral compressive sensing [J].ApplComputHarmonicAnal, 2013,35(1)： 111-129.

[5] STOICA P, BABU P, LI J. SPICE： A sparse covariance-based estimation method for array processing [J].IEEETransactionsonSignalProcessing, 2011,59(2)： 629-638.

[6] BEN-HAIM Z, ELDAR Y C, ELAD M. Coherence-based performance guarantees for estimating a sparse vector under random noise [J].IEEETransactionsonSignalProcessing, 2010,58(10)： 5030-5043.

[7] MIDDLETON D. Non-Gaussian noise models in signal processing for telecommunications： New methods and results for class A and class B noise models [J].IEEETransactionsonInformationTheory, 1999,45(4)： 1129-1149.

[8] ZENG W J, SO H C, JIANG X. Outlier-robust greedy pursuit algorithms inlp-space for sparse approximation [J].IEEETransactionsonSignalProcessing, 2016,64(1)： 60-75.

[9] ROJAS C R, KATSELIS D,HJALMARSSON H. A note on the SPICE method [J].IEEETransactionsonSignalProcessing, 2013,61(18)： 4545-4551.

[10] TANG G, BHASKAR B N, SHAH P,et al. Compressed sensing off the grid [J].IEEETransactionsonInformationTheory, 2013,59(11)： 7465-7490.

[11] BHASKAR B N, TANG G, RECHT B. Atomic norm denoising with applications to line spectral estimation [J].IEEETransactionsonSignalProcessing, 2013,61(23)： 5987-5999.

[12] YANG Z, XIE L. On gridless sparse methods for line spectral estimation from complete and incomplete data [J].IEEETransactionsonSignalProcessing, 2015,63(12)： 3139-3153.

[13] HAN K, NEHORAI A. Improved source number detection and direction estimation with nested arrays and ULAs using jackknifing [J].IEEETransactionsonSignalProcessing, 2013,61(23)： 6118-6128.

[15] YANG Z, XIE L. Exact joint sparse frequency recovery via optimization methods [J].IEEETransactionsonSignalProcessing, 2016,64(19)： 5145-5157.

[16] BOYD S, PARIKH N, CHU E, et al. Distributed optimization and statistical learning via the alternating direction method of multipliers [J].Foundations&TrendsinMachineLearning, 2011,3(1)： 1-122.

[17] SCHMIDT R. Multiple emitter location and signal parameter estimation [J].IEEETransactionsAntennPropag, 1986,34(3)： 276-280.

[18] SUN Y, BABU P, PALOMAR D. Majorization-Minimization algorithms in signal processing, communications, and machine learning [J].IEEETransactionsonSignalProcessing, 2017,65(3)： 794-816.