(重慶郵電大學(xué) 信號(hào)與信息處理重慶市重點(diǎn)實(shí)驗(yàn)室，重慶 400065)

1 引 言

寬帶Chirp信號(hào)作為一種大的時(shí)間帶寬積信號(hào)，廣泛應(yīng)用于各種信息系統(tǒng)，如通信、雷達(dá)、聲納和地震勘探等[1]，它也是一種特殊的非平穩(wěn)信號(hào)。在這類系統(tǒng)中，Chirp信號(hào)的檢測(cè)、參數(shù)和波達(dá)方向角的估計(jì)是一個(gè)重要的研究課題。

關(guān)于Chirp信號(hào)的參數(shù)估計(jì)，比較有代表性的是基于空間時(shí)頻分布(STFD)的方法[2-3]。該類方法通過(guò)計(jì)算陣元間的互Wigner-Ville分布(WVD)構(gòu)造時(shí)頻域數(shù)據(jù)向量，然后利用Chirp信號(hào)的局部窄帶特性，截取時(shí)頻分布譜峰附近的時(shí)頻點(diǎn)，對(duì)截取出的時(shí)頻點(diǎn)數(shù)據(jù)進(jìn)行特征值分解，再利用子空間類方法[4-5]實(shí)現(xiàn)信號(hào)參數(shù)估計(jì)。但是，該類算法只在信噪比較高的情況下才能有較好的性能，而且在多目標(biāo)的情況下，信號(hào)之間的交叉項(xiàng)將嚴(yán)重影響Chirp信號(hào)的參數(shù)估計(jì)性能，弱目標(biāo)很容易被交叉項(xiàng)所掩蓋。

對(duì)于稀疏分解在陣列信號(hào)處理領(lǐng)域中的應(yīng)用，已經(jīng)有人把稀疏分解中常用的匹配追蹤(MP)算法[6]引入到陣列信號(hào)處理領(lǐng)域中，實(shí)現(xiàn)了對(duì)信號(hào)參數(shù)和波達(dá)方向角的估計(jì)。而本文提出的基于正交匹配追蹤(Orthogonal Matching Pursuit,OMP)[7]分解的寬帶Chirp信號(hào)多參數(shù)估計(jì)方法與MP算法相比有更好的收斂性，在一定的過(guò)完備原子庫(kù)密度條件下，本文所述算法能很好地實(shí)現(xiàn)對(duì)Chirp信號(hào)起始頻率和調(diào)頻斜率的估計(jì)，且波達(dá)方向角(DOA)估計(jì)的性能高于WVD算法[8]，尤其在多信源和低信噪比情況下，算法仍然有良好的性能。

2 信號(hào)數(shù)學(xué)模型

假設(shè)M元線陣沿x軸均勻放置，陣元間距為d。第一個(gè)陣元位于坐標(biāo)原點(diǎn)，稱為參考陣元。假定有P個(gè)寬帶Chirp信號(hào)從不同方向入射到M元陣列上，則第m個(gè)陣元的輸出為

(1)

(2)

將式(2)代入式(1)可得：

(3)

(4)

將其寫(xiě)為向量形式可以得到陣列瞬時(shí)窄帶模型：

Aθ,tSt+Nt

(5)

式中，A(θ,t)=[a1(θ,t),a2(θ,t),…,aP(θ,t)]為M×P的陣列方向矩陣，包含了信源的方位信息，ai(θ,t)=[1,e-j2π(fi,0+kit)τ2,i,e-j2π(fi,0+kit)τ3,i,…,e-j2π(fi,0+kit)τM,i]T；P×1的入射信號(hào)向量為S(t)=[s1(t),s2(t),…,sP(t)]T;N(t)=[n1(t),n2(t),…,nM(t)]T為M×1的加性噪聲向量。

假設(shè)取N次快拍，問(wèn)題就變?yōu)樵谝阎獂(n),n=1,2,3，…,N的情況下，估計(jì)信號(hào)參數(shù)(fi,0,ki,θi),i=1,2,3，…,P。

3 基于OMP分解的寬帶Chirp信號(hào)多參數(shù)估計(jì)

給定一個(gè)集合D={gk，k=1,2,3,…,K}，其元素是張成整個(gè)Hilbert空間H=RN的單位矢量，K>>>N,我們稱集合D為原子庫(kù)，其元素為原子。集合Γ={γk，k=1,2,3,…,K}為原子參數(shù)組成的集合，其元素γk為每個(gè)原子的參數(shù)。

假設(shè)在D中所有原子的范數(shù)都為一，算法在初始化時(shí)，設(shè)R0x=x，可以首先對(duì)信號(hào)x進(jìn)行第一次OMP分解，從過(guò)完備庫(kù)中選出與信號(hào)x最為匹配的原子gγ0，即：

R0x=〈R0x,gγ0〉gγ0+R1x

(6)

式中，〈R0x,gγ0〉gγ0是R0x對(duì)gγ0的投影，R1x是剩余量。顯然，gγ0是與R1x正交的，所以可得到下式：

‖R0x‖2=〈R0x,gγ0〉2+‖R1x‖2

(7)

初始化u0=gγ0，對(duì)m≥0，OMP挑選gγm，使得：

(8)

利用Gram-Schmidt算法[8]將gγm關(guān)于gγp0≤p≤m正交化，定義：

(9)

對(duì)um歸一化：

(10)

將余項(xiàng)Rmf投影到um(而不是gγm)上，得到：

Rmx=〈Rmx,um〉um+Rm+1x

(11)

將此方程對(duì)0≤m

(12)

式中，Pvk是在um0≤m

〈Rmx,um〉=〈Rmx,gγm〉

(13)

因Vk的維數(shù)為k,故存在M使得x∈VM，從而Rmx=0。將式(13)代入式(12)，并令k=M,得：

(14)

做有限M次迭代可以得到收斂性。它是x在一個(gè)正交向量族上的分解，故：

(15)

此過(guò)程即為OMP的基本步驟，因?yàn)樾盘?hào)殘差和前面每一步中從原子庫(kù)中選擇的原子相正交，在有限維空間中，OMP能在有限次迭代后收斂。

3.1 寬帶Chirp信號(hào)起始頻率和調(diào)頻斜率的估計(jì)

假設(shè)有P個(gè)寬帶Chirp信號(hào)從不同方向入射到M元陣列上，則參考陣元接收數(shù)據(jù)為

(16)

對(duì)P個(gè)寬帶Chirp信號(hào)進(jìn)行離散化，得：

(17)

式中，i=1,2,3,…,P；n=1,2,3,…,N；fs為采樣頻率。

根據(jù)寬帶Chirp信號(hào)形式，建立過(guò)完備原子庫(kù)，其基本原子為

(18)

原子gγ的長(zhǎng)度與信號(hào)x1n本身長(zhǎng)度相同，原子參數(shù)組為γ=f,k，f、k分別對(duì)應(yīng)Chirp信號(hào)的起始頻率和調(diào)頻斜率參量。假設(shè)f的取值范圍為f∈0,0.6，k的取值范圍為k∈0,0.1，f、k按照需要的精度取值(本文設(shè)原子庫(kù)的搜索精度為0.001)，構(gòu)造出過(guò)完備原子庫(kù)Df=gγγ∈Γ：

Df=gγγ∈Γ=

(19)

式中，Γ為參數(shù)組γ所屬的集合，其值視具體情況而定。

由于OMP分解的特性，對(duì)Chirp信號(hào)進(jìn)行分解時(shí)，將在某一分量上達(dá)到最大匹配，而噪聲不具有Chirp信號(hào)的特征，因此在這一分量上投影接近零，從而可以達(dá)到濾除噪聲的目的。將淹沒(méi)在噪聲中的多個(gè)Chirp信號(hào)進(jìn)行OMP分解，將會(huì)在不同分量上得到最大匹配，根據(jù)這些與信號(hào)最大匹配的原子參數(shù)，即可估計(jì)出每個(gè)Chirp信號(hào)的起始頻率和調(diào)頻斜率。

經(jīng)過(guò)P步OMP分解后，可以得到P個(gè)原子gγi和原子系數(shù)αi=〈Ri-1xin,gγi〉。根據(jù)原子的參數(shù)γi=fi,ki就可以得到P個(gè)Chirp信號(hào)的起始頻率和調(diào)頻斜率參數(shù)估計(jì)：

(20)

(21)

由此，我們可以重新構(gòu)造出P個(gè)寬帶Chirp信號(hào)：

(22)

在得到Chirp信號(hào)具體形式的基礎(chǔ)上，我們?cè)龠M(jìn)行DOA估計(jì)。

3.2 寬帶Chirp信號(hào)的DOA估計(jì)

首先，針對(duì)各個(gè)不同的Chirp信號(hào)分別建立過(guò)完備庫(kù)的原子：

(23)

Gθgl=g1,g2,g3,…,gL

(24)

整個(gè)原子庫(kù)的原子數(shù)為L(zhǎng)個(gè)。

通過(guò)在原子庫(kù)Gθgl上對(duì)陣列信號(hào)用OMP算法做稀疏分解，由所獲得的最佳原子便可以得到信源DOA估計(jì)參數(shù)，我們通過(guò)計(jì)算陣列接收數(shù)據(jù)與原子庫(kù)Gθgl中原子之間的互相關(guān)矩陣的跡，搜索跡的最大峰值來(lái)獲得。假定互相關(guān)矩陣Rxy的m個(gè)特征值為λ1,λ2,λ3,…,λm。由矩陣的跡的概念可知：

(25)

即互相關(guān)矩陣Rxy的特征值之和trRxy反映了兩信號(hào)xt和yt之間的相關(guān)程度。由此，分別計(jì)算陣列接收數(shù)據(jù)矩陣xn與原子庫(kù)Gθgl中原子gl的互相關(guān)矩陣：

(26)

并求解互相關(guān)矩陣Rxg的跡，通過(guò)比較跡的大小獲得最佳匹配原子gi的參數(shù)，進(jìn)一步得到DOA估計(jì)：

(27)

綜合上述分析，將基于OMP分解的寬帶Chirp信號(hào)參數(shù)估計(jì)算法步驟總結(jié)如下：

(1)建立過(guò)完備庫(kù)的原子gγ，如式(18)所示，進(jìn)一步建立過(guò)完備原子庫(kù)Df=gγγ∈Γ，如式(19)所示；

(3)建立過(guò)完備庫(kù)的原子gl，如式(23)所示；進(jìn)一步建立原子庫(kù)Gθgl，如式(24)所示；

(5)多信號(hào)時(shí)重復(fù)第3步和第4步，得到各個(gè)信號(hào)的DOA估計(jì)。

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

為體現(xiàn)OMP算法中信號(hào)原子之間的正交性，考慮5個(gè)寬帶Chirp信號(hào)從遠(yuǎn)場(chǎng)入射到八陣元均勻線陣，信號(hào)的起始頻率和調(diào)頻斜率分別如下：

信號(hào)1：f10=0.150 000 1,K1=0.033 000 000 01;

信號(hào)2：f20=0.250 000 1,K2=0.066 000 000 01;

信號(hào)3：f30=0.350 000 1,K3=0.099 000 000 01;

信號(hào)4：f40=0.450 000 1,K4=0.055 000 000 01;

信號(hào)5：f50=0.550 000 1,K5=0.077 000 000 01。

5個(gè)信號(hào)的波達(dá)方向角分別為10°、20°、30°、40°、50°，陣元間距為信號(hào)最高頻率對(duì)應(yīng)的半波長(zhǎng)。

(1)實(shí)驗(yàn)1：用OMP算法和WVD算法進(jìn)行DOA估計(jì)的空間譜

(a)OMP的空間譜曲線

(b)WVD的空間譜曲線圖1 五信源的DOA估計(jì)的空間譜Fig.1 Spatial spectrum curve of DOA estimation of five source signals

從圖1中可見(jiàn)，OMP算法準(zhǔn)確有效地估計(jì)出了5個(gè)信源的波達(dá)方向，在波達(dá)方向角為10°、20°、30°、40°、50°處，空間譜出現(xiàn)了譜峰，而WVD算法在40°處出現(xiàn)偏差。仿真結(jié)果表明，OMP算法的估計(jì)性能優(yōu)于WVD算法。

(2)實(shí)驗(yàn)2：用OMP算法對(duì)寬帶Chirp信號(hào)起始頻率和調(diào)頻斜率的估計(jì)

信號(hào)采樣長(zhǎng)度取N=512。用OMP算法估計(jì)信號(hào)1、信號(hào)2、信號(hào)3的起始頻率和調(diào)頻斜率。信噪比從-20～20 dB,每隔1 dB進(jìn)行100次Monte-Carlo實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果如圖2和圖3所示。

(a)信號(hào)1

(b)信號(hào)2

(c)信號(hào)3圖2 3個(gè)信號(hào)的起始頻率估計(jì)Fig.2 The starting frequency estimation of three signals

(a)信號(hào)1

(b)信號(hào)2

(c)信號(hào)3圖3 3個(gè)信號(hào)的調(diào)頻斜率估計(jì)Fig.3 Modulation frequency estimation of three signals

由圖2和圖3可見(jiàn)，OMP算法很好地實(shí)現(xiàn)了對(duì)Chirp信號(hào)起始頻率和調(diào)頻斜率的估計(jì)。在信噪比為-20 dB時(shí)，OMP算法對(duì)起始頻率和調(diào)頻斜率的估計(jì)性能很好，估計(jì)誤差較小。3個(gè)信號(hào)的起始頻率估計(jì)誤差分別為0.11、0.15、0.14，調(diào)頻斜率估計(jì)誤差分別為0.03、0.03、0.053；在信噪比門(mén)限值-2 dB時(shí)，起始頻率和調(diào)頻斜率估計(jì)誤差(RMSE)值達(dá)到了克拉美羅限，并且當(dāng)信噪比進(jìn)一步增加時(shí)，由于本文所設(shè)的原子庫(kù)的搜索精度為0.001，OMP算法只能估計(jì)到起始頻率和調(diào)頻斜率真實(shí)數(shù)值的小數(shù)點(diǎn)后3位，且估計(jì)值不再隨信噪比的增加而變化，故估計(jì)誤差(RMSE)不隨信噪比的增加而增加，在仿真結(jié)果中表現(xiàn)為一條水平的線段；當(dāng)信噪比在信噪比門(mén)限值-2 dB以下時(shí)，OMP算法對(duì)起始頻率和調(diào)頻斜率的估計(jì)誤差隨信噪比的增加越來(lái)越小。

(3)實(shí)驗(yàn)3:用OMP算法和WVD算法對(duì)寬帶Chirp信號(hào)的DOA的估計(jì)性能與快拍數(shù)的關(guān)系

(a)信號(hào)1

(b)信號(hào)2

(c)信號(hào)3圖4 3個(gè)信號(hào)DOA估計(jì)的RMSEFig.4 RMSE of three signals DOA estimation

固定信噪比為10 dB，分別采用OMP算法和WVD算法，信號(hào)采樣長(zhǎng)度N(取值從100到1 000)每隔100進(jìn)行200次Monte-Carlo仿真實(shí)驗(yàn)。圖4給出了這兩種算法在信噪比為10 dB時(shí)，在不同快拍數(shù)下的3個(gè)信號(hào)DOA估計(jì)的RMSE曲線。

從圖4可見(jiàn)，在固定信噪比為10 dB條件下，隨著信號(hào)采樣長(zhǎng)度的增加，OMP算法的DOA估計(jì)的均方誤差減小，且估計(jì)性能明顯優(yōu)于WVD方法。

(4)實(shí)驗(yàn)4:用OMP算法和WVD算法對(duì)寬帶Chirp信號(hào)的DOA估計(jì)性能與信噪比的關(guān)系

實(shí)驗(yàn)中分別采用OMP算法(信噪比從-20～20 dB)和WVD算法(信噪比從-10～20 dB)，信噪比每隔2 dB進(jìn)行200次的Monte-Carlo實(shí)驗(yàn)。其中，陣元數(shù)M=8,信號(hào)采樣長(zhǎng)度N=512。結(jié)果如圖5所示。

(a)信號(hào)1

(b)信號(hào)2

(c)信號(hào)3圖5 3個(gè)信號(hào)DOA估計(jì)的RMSEFig.5 RMSE of three signals DOA estimation

由圖5明顯可見(jiàn)，OMP算法對(duì)信號(hào)DOA估計(jì)效果非常好。隨著信噪比的提高，DOA估計(jì)的均方誤差越來(lái)越小，估計(jì)精度明顯優(yōu)于WVD方法。在信噪比為-14 dB時(shí)，該算法仍具有很好的性能，而WVD方法在-10 dB幾乎失效。

5 結(jié) 論

本文通過(guò)大量的仿真實(shí)驗(yàn)證明了在多信源低信噪比的情況下，OMP對(duì)寬帶陣列信號(hào)參數(shù)估計(jì)的有效性。論文的研究尤其在以下方面得到了很大改進(jìn)：在信噪比為-20 dB時(shí)，OMP算法對(duì)起始頻率和調(diào)頻斜率的估計(jì)性能仍很好；在信噪比門(mén)限值為-2 dB時(shí)，該方法對(duì)起始頻率和調(diào)頻斜率的估計(jì)的RMSE值達(dá)到了克拉美羅限；在信噪比為-14 dB時(shí)，該方法對(duì)波達(dá)方向角度的估計(jì)仍然有良好的性能，且其估計(jì)精度明顯優(yōu)于傳統(tǒng)的WVD方法。

參考文獻(xiàn)：

[1] 張賢達(dá)，保錚.非平穩(wěn)信號(hào)分析與處理[M].北京：國(guó)防工業(yè)出版社，1998.

ZHANG Xian-da, BAO Zheng.Non-stationary signal analysis and process[M].Beijing:National Defense Industry Press, 1998.(in Chinese)

[2] Belouchrani A, Amin M G. Blind source separation based on time-frequency signal representations [J]. IEEE Transactions on Signal Processing, 1998,46(11):2888-2897.

[3] Belouchrani A, Amin M G. Time-frequency MUSIC [J].IEEE Signal Processing Letters,1999, 6(5):109-110.

[4] Schmidt R O. Multiple emitter location and signal parameter estimation [J].IEEE Transactions on Antennas and Propagation, 1986,34(3):276-280.

[5] Roy R,Kailath T.ESPRIT-estimation of signal parameters via rotational invariance technique[J].IEEE Transactions on Acoust, Speech and Signal Processing,1989,37(7):984-995.

[6] Mallat S，Zhang Z. Matching Pursuit with Time-frequency Dictionaries[J].IEEE Transactions on Signal Processing, 1993, 41(12):3397-3415.

[7] 楊愚.圖像稀疏分解快速實(shí)現(xiàn)與初步應(yīng)用[D].成都：西南交通大學(xué)出版社，2007.

YANG Yu.Fast implement and preliminary application of picture sparse decomposition[D].Chengdu: Southwest Jiaotong University Press, 2007. (in Chinese)

[8] 張賢達(dá).現(xiàn)代信號(hào)處理[M].北京：清華大學(xué)出版社，2002.

ZHANG Xian-da.Modern Signal Processing [M].Beijing:Tsinghua University Press, 2002. (in Chinese)