王悅斌，蔣景飛，張建秋，＊

1.復(fù)旦大學(xué) 智慧網(wǎng)絡(luò)系統(tǒng)研究中心和電子工程系，上海 200433 2.電子信息控制重點(diǎn)實(shí)驗(yàn)室，成都 610036

非平穩(wěn)信號(hào)的分析和參數(shù)估計(jì)，在雷達(dá)［1］、語(yǔ)音分析［2］和生物醫(yī)學(xué)［3］等領(lǐng)域中正發(fā)揮著越來(lái)越重要的作用。然而，傳統(tǒng)的譜估計(jì)算法只能分析出全局譜，無(wú)法揭示信號(hào)頻率的局部變化。為了獲取不同時(shí)刻信號(hào)分量的存在與否，以及各頻率成分隨時(shí)間變化情況，時(shí)頻分布（Time-Frequency Distribution，TFD）的分析方法應(yīng)運(yùn)而生。其中短時(shí) 傅 里 葉 變 換 （Short-Time Fourier Transform，STFT）、Wigner-Ville分布和S 變換［4］是3種典型的時(shí)頻分析算法。這些方法簡(jiǎn)單易懂，不需要額外的參數(shù)模型，應(yīng)用十分方便。

可是，對(duì)多分量非穩(wěn)態(tài)信號(hào)來(lái)說(shuō)，這些傳統(tǒng)算法在頻率分辨率，非線性調(diào)頻特性，以及動(dòng)態(tài)時(shí)頻譜（據(jù)文獻(xiàn)［5］的定義：若時(shí)頻譜中的信號(hào)分量隨著時(shí)間動(dòng)態(tài)地出現(xiàn)和／或消失，則稱其為“動(dòng)態(tài)時(shí)頻譜”）的分析等方面還存在局限性。例如，Wigner-Ville分布對(duì)于單分量信號(hào)具有較好的時(shí)頻分析能力，而對(duì)于多分量信號(hào)來(lái)說(shuō)，則會(huì)引入較強(qiáng)的交叉項(xiàng)［6］；S變換雖然具有可變的時(shí)頻分辨率，但由于窗寬度是頻率的函數(shù)，因而在不同的信號(hào)分析中，會(huì)表現(xiàn)出不同的能量集中度，從而限制了應(yīng)用范圍［7］。

為了提高對(duì)非平穩(wěn)信號(hào)的時(shí)頻表示（Time-Frequency Representations，TFRs）能力，近年來(lái)，人們研究出了許多新的算法。其中，廣義S變換（Generalized S-Transform，GST）［8］通過(guò)多參數(shù)的窗函數(shù)，為克服時(shí)頻分析方法所存在時(shí)頻分辨率單一的缺點(diǎn)提供了一條新途徑。然而，由于依賴于窗函數(shù)的選擇，GST無(wú)法與不同類型的信號(hào)相適應(yīng)。不同于GST，文獻(xiàn)［9］提出的自適應(yīng) 的 STFT （Adaptive Short-Time Fourier Transform，ASTFT）算法，可自適應(yīng)非平穩(wěn)信號(hào)特性，即據(jù)信號(hào)調(diào)頻斜率的變化來(lái)不斷調(diào)節(jié)自身STFT窗的長(zhǎng)度?？墒牵瑢?duì)于具有多個(gè)分量且調(diào)頻特性較復(fù)雜的時(shí)頻信號(hào)，該算法不能獲得較為可靠的調(diào)頻斜率變化信息。匹配解調(diào)變換（Matching Demodulation Transform，MDT）［10］可以獲得較高的瞬時(shí)頻率估計(jì)精度，并且能有效提高時(shí)頻表示的能量集中度，但該算法需要預(yù)先給定信號(hào)的全局參數(shù)模型，而大多數(shù)情況下瞬時(shí)頻率的模型是未知的，因此該算法實(shí)用性較低。文獻(xiàn)［11］提出了不需要信號(hào)全局參數(shù)模型的變分非線性調(diào)頻模式分解（Variational Nonlinear Chirp Mode Decomposition，VNCMD）算法。該算法對(duì)于非線性調(diào)頻信號(hào)的時(shí)頻分析，在具有較好性能的同時(shí)，也具有較高的時(shí)頻分辨率，并且可以解決信號(hào)頻率分量交叉問(wèn)題。然而，該算法需要預(yù)知信號(hào)分量數(shù)目，因而難以解決動(dòng)態(tài)出現(xiàn)和消失多分量時(shí)頻信號(hào)的分析問(wèn)題。為解決這一問(wèn)題，文獻(xiàn)［12］提出一種基于Capon譜估計(jì)的多分量信號(hào)檢測(cè)和追蹤算法。該算法首先利用Capon譜估計(jì)得到時(shí)頻譜的粗估計(jì)，進(jìn)而估計(jì)出噪聲譜，再進(jìn)行信號(hào)峰值檢測(cè)，最后進(jìn)行匹配追蹤，獲得了去噪效果較好的時(shí)頻分析結(jié)果。然而，該算法需要依靠全局譜判斷信號(hào)分量的存在與否，因而實(shí)時(shí)性較差。最近，文獻(xiàn)［5］提出一種利用狄利克雷混合模型來(lái)描述和估計(jì)動(dòng)態(tài)時(shí)頻譜的算法，進(jìn)而給出了一種基于自適應(yīng)狄利克雷混合模型的Rao-Blackwellised粒子濾波（Adaptive Dirichlet Process Mixture Model-based Rao-Blackwellised Particle Filter，ADPM-RBPF）算法。該算法實(shí)時(shí)性較好，提高了低信噪比下的時(shí)頻譜估計(jì)的精度?？墒?，由于該算法利用STFT得到的原始時(shí)頻譜作為觀測(cè)，時(shí)頻分辨率較差。

針對(duì)動(dòng)態(tài)出現(xiàn)和消失多分量非平穩(wěn)信號(hào)，本文提出了一種基于隨機(jī)有限集的時(shí)頻分析法。該算法利用時(shí)頻變換在每個(gè)時(shí)刻獲得的信號(hào)分量幅度的局部極值為測(cè)量的隨機(jī)有限集，然后將時(shí)頻信號(hào)各分量幅度和頻率的變化描述為多項(xiàng)式預(yù)測(cè)模型［13］，進(jìn)而為每個(gè)信號(hào)分量建立起了一種新的線性高斯?fàn)顟B(tài)空間模型。這樣，本文就將動(dòng)態(tài)出現(xiàn)和消失的時(shí)頻譜分析、探測(cè)和跟蹤問(wèn)題，歸納為可利用隨機(jī)有限集進(jìn)行多目標(biāo)追蹤的問(wèn)題。借助于高斯混合概率假設(shè)密度（Gaussian Mixture Probability Hypothesis Density，GM-PHD）濾波器［14］，本文算法可以分析動(dòng)態(tài)的時(shí)頻譜。對(duì)于GM-PHD濾波器，初始權(quán)重的賦值和權(quán)重的更新尤為關(guān)鍵。據(jù)時(shí)頻分布的特點(diǎn)，本文提出了一種初始權(quán)重賦值算法，結(jié)合提出的譜分量幅度和頻率的聯(lián)合似然函數(shù)，進(jìn)而解決了權(quán)重的更新問(wèn)題。分析和數(shù)值仿真結(jié)果均證明：所提算法可有效提升動(dòng)態(tài)時(shí)頻譜的跟蹤精度，對(duì)微弱時(shí)頻譜分量的探測(cè)能力，以及對(duì)載頻差異的分析能力。

1 時(shí)頻分析的隨機(jī)有限集模型

對(duì)于一個(gè)由多個(gè)調(diào)頻調(diào)幅分量構(gòu)成的時(shí)頻信號(hào)，為了描述其分量隨時(shí)間動(dòng)態(tài)出現(xiàn)和消失的情況，采用Kn來(lái)表示隨離散時(shí)間n變化的信號(hào)分量數(shù)目，這樣本文考慮的時(shí)頻信號(hào)的離散形式可描述為［5，12］

式中：yn為信號(hào)的觀測(cè)值；akn和fkn分別為第k個(gè)分量信號(hào)的幅度和頻率函數(shù)；θk0為第k個(gè)分量信號(hào)的初始相位；fs為采樣頻率；ηn為觀測(cè)噪聲。對(duì)時(shí)頻信號(hào)式（1）進(jìn)行STFT時(shí)頻變換，可以得到信號(hào)的時(shí)頻分布為［15］

式中：n和τ分別為離散的時(shí)間和頻率；N 為STFT短時(shí)窗的長(zhǎng)度；ωm為窗函數(shù)。由于時(shí)頻分布TFDN（n，τ）反映的是信號(hào)在（n，τ）時(shí)頻鄰域內(nèi)的能量，那么任意時(shí)刻信號(hào)分量的能量，就必定會(huì)集中在該時(shí)刻的瞬時(shí)頻率（Instantaneous Frequency，IF）附近。如果視每一時(shí)刻時(shí)頻分布TFDN（n，τ）的局部極大值即式（3）中的Zn為測(cè)量值，這意味著任意時(shí)刻的信號(hào)分量必定包含在該時(shí)刻的Zn中。這是因?yàn)槿我鈺r(shí)刻時(shí)頻分布TFDN（n，τ）的能量是信號(hào)能量和噪聲能量的疊加，對(duì)高斯白噪聲而言，其理想時(shí)頻分布的能量是某一常數(shù)，這樣當(dāng)一個(gè)常數(shù)加上任一小的信號(hào)能量時(shí)，都將為高斯白噪聲的理想時(shí)頻分布能量增加一個(gè)局部極值。

考慮到微弱信號(hào)分量的幅值可能很小，而觀測(cè)噪聲可能是非高斯的，以及高斯觀測(cè)噪聲的樣本可能不是足夠大，因此本文將所有檢測(cè)到的局部極值均視為測(cè)量，如圖1所示。由于觀測(cè)噪聲的存在，以及時(shí)頻信號(hào)的分量可能隨機(jī)出現(xiàn)或消失，那么每一時(shí)刻的局部極值也將是隨機(jī)的。這樣，可把Zn視為一個(gè)測(cè)量隨機(jī)有限集（RFS）集合［16－17］：

式中：Jn為n時(shí)刻局部極值（測(cè)量值）的個(gè)數(shù)；zjn＝［］T為n時(shí)刻第j個(gè)局部極值幅度和頻率的測(cè)量值集合，其中j＝1，2，…，Jn；Ez為多分量時(shí)頻信號(hào)的測(cè)量空間。因此，n時(shí)刻的測(cè)量隨機(jī)集可以進(jìn)一步表示為

式中：ηn為n時(shí)刻源于觀測(cè)噪聲和／或雜波的測(cè)量RFS；Θn（x）為n時(shí)刻源于時(shí)頻信號(hào)多分量狀態(tài)x的測(cè)量RFS；∪為并運(yùn)算。視信號(hào)分量為目標(biāo)，則多分量信號(hào)的目標(biāo)狀態(tài)RFS就可描述為［16-17］

式中：Kn為n時(shí)刻目標(biāo)個(gè)數(shù)（時(shí)頻信號(hào)分量的個(gè)數(shù)）；xkn為n時(shí)刻第k個(gè)目標(biāo)的狀態(tài)，其中k＝1，2，…，Kn；Ex為多目標(biāo)狀態(tài)空間。n時(shí)刻的多目標(biāo)狀態(tài)隨機(jī)集可以進(jìn)一步表示為

式中：Sn｜n－1（xn－1）為n－1時(shí)刻的目標(biāo)狀態(tài)xn－1在n時(shí)刻的存活RFS；Γn為n時(shí)刻新生分量目標(biāo)RFS。這樣，可推導(dǎo)出基于RFS的全體多目標(biāo)聯(lián)合后驗(yàn)概率密度的貝葉斯遞推公式［16－17］：

圖1 峰值檢測(cè)示意圖Fig.1 Diagram of peak detection

式中：Fn｜n－1（·｜·）為 多 目 標(biāo) 轉(zhuǎn) 移 概 率 密 度；gn（·｜·）為多目標(biāo)似然函數(shù)；Z1：n＝ ｛Z1，Z2，…，Zn｝為前n 個(gè)時(shí)刻所有的測(cè)量 RFS；pn｜n－1（·｜Z1：n－1）為多目標(biāo)預(yù)測(cè)概率密度；pn｜n（·｜Z1：n）為多目標(biāo)聯(lián)合后驗(yàn)概率密度。

式（8）和式（9）的遞推公式需要在多目標(biāo)狀態(tài)空間上計(jì)算集積分，計(jì)算量非常大，在實(shí)際計(jì)算中一般難以實(shí)現(xiàn)。在線性高斯條件下，基于RFS的GM-PHD，通過(guò)把目標(biāo)的概率假設(shè)密度（Probability Hypothesis Density，PHD）近似為混合高斯形式，可以得到遞推形式式（8）和式（9）的閉合解［14］。針對(duì)式（4）和式（6）的 RFS變量，本文將為其提出一種新的線性高斯?fàn)顟B(tài)空間模型，使得本文可以利用GM-PHD濾波計(jì)算目標(biāo)數(shù)量并估計(jì)各個(gè)目標(biāo)的狀態(tài)，進(jìn)而得到了一種可分析動(dòng)態(tài)出現(xiàn)和／或消失多分量時(shí)頻信號(hào)的新復(fù)算法。

2 多分量時(shí)頻信號(hào)分析

2.1 時(shí)頻分量信號(hào)的狀態(tài)空間模型

考慮到式（1）中信號(hào)分量的幅度akn和頻率fkn，既可能是線性，也可能是非線性的，為了使討論具有一般性，本文假設(shè)式（1）中的akn和fkn都是非線性函數(shù)，此時(shí)線性就可認(rèn)為是非線性的特例。對(duì)于任一非線性函數(shù)，Weierstrass逼近定理［18］可知：任一有界閉區(qū)間上的連續(xù)函數(shù)都可以由該區(qū)間內(nèi)的多項(xiàng)式以任意精度逼近。這意味著對(duì)任意連續(xù)函數(shù)，如果將其閉區(qū)間進(jìn)行分割，且讓每個(gè)分割出的小區(qū)間足夠小，以致在每一小區(qū)間內(nèi)，該函數(shù)都可以用一個(gè)低階多項(xiàng)式（例如一階多項(xiàng)式）以任一精度逼近。據(jù)這個(gè)原理，假設(shè)信號(hào)分量的瞬時(shí)頻率fkn，在每一分割的小區(qū)間內(nèi)可由L階多項(xiàng)式描述：

式中：p（l）為多項(xiàng)式的系數(shù)。如果想用前M個(gè)時(shí)刻的頻率值來(lái)預(yù)測(cè)當(dāng)前時(shí)刻的頻率，即

顯然，式（11）是以h（m）（m＝1，2，…，M）為系數(shù)的有限沖激響應(yīng)（FIR）濾波器。將式（10）代入式（11），消去多項(xiàng)式系數(shù)p（l），可得

由式（12）無(wú)法唯一確定系數(shù)h（m）?？紤]到信號(hào)的觀測(cè)通常都是存在噪聲的，這樣就需要濾波器的噪聲增益：

具有最小值。在式（12）的約束條件下，利用拉格朗日法就可求得式（13）的最優(yōu)解為［13］：

1）當(dāng)L＝1時(shí)

2）當(dāng)L＝2時(shí)

式中：m＝1，2，…，M，其中M 表示濾波器的長(zhǎng)度。當(dāng)L為其他值時(shí)，計(jì)算結(jié)果詳見(jiàn)文獻(xiàn)［13］。

上面的分析表明：多項(xiàng)式預(yù)測(cè)模型式（11）與FIR預(yù)測(cè)濾波器在數(shù)學(xué)上嚴(yán)格等效。據(jù)上述多項(xiàng)式的預(yù)測(cè)模型，就可對(duì)時(shí)頻信號(hào)分量的瞬時(shí)幅度和頻率函數(shù)進(jìn)行建模。定義第k個(gè)時(shí)頻信號(hào)分量在n時(shí)刻的狀態(tài)向量為

式（17）也可視作零階多項(xiàng)式的預(yù)測(cè)模型，即用前一時(shí)刻的幅度值來(lái)預(yù)測(cè)當(dāng)前時(shí)刻的幅度值。根據(jù)式（11）和式（17）的模型描述，式（16）的狀態(tài)方程就可以表示為

有了式（18）幅度和頻率的聯(lián)合狀態(tài)方程后，還需要建立信號(hào)分量幅度和頻率作為聯(lián)合參數(shù)的測(cè)量方程。假設(shè)第k個(gè)信號(hào)分量的測(cè)量為zn＝［za，n，zf，n］T，其中za，n和zf，n分別為幅度和頻率的測(cè)量值，則式（16）的測(cè)量方程可以表示為

式中：σn以0均值的加性噪聲來(lái)描述式（3）的不確定性，噪聲方差為R，觀測(cè)矩陣為H ＝

在式（18）和式（20）所描述多分量時(shí)頻信號(hào)的分析中，必須確定式（20）的測(cè)量zn屬于式（1）中哪一個(gè)信號(hào)分量或雜波，才能根據(jù)式（18）和式（20）的狀態(tài)空間模型進(jìn)行時(shí)頻分析，2.2節(jié)將討論。

2.2 GM-PHD濾波的時(shí)頻分析法

通過(guò)時(shí)頻變換，本文將非線性觀測(cè)方程式（1）轉(zhuǎn)換成了式（20）的線性觀測(cè)方程，而2.1節(jié)對(duì)線性／非線性調(diào)頻函數(shù)，通過(guò)多項(xiàng)式預(yù)測(cè)方法，則將其變換成了線性的狀態(tài)方程式（18）。這樣式（18）和式（20）構(gòu)成的時(shí)頻分析狀態(tài)空間模型是線性的。此外，由于本文假設(shè)噪聲是高斯的，因此采用了GM-PHD濾波器。當(dāng)式（1）觀測(cè)噪聲是非高斯時(shí)，則可以采用SMC-PHD［20］等濾波器，但建立的時(shí)頻分析狀態(tài)空間模型式（18）和式（20）有效，且將簡(jiǎn)化SMC-PHD等濾波器的計(jì)算復(fù)雜度。

2.2.1 算法預(yù)測(cè)

假設(shè)n－1時(shí)刻時(shí)頻分量的后驗(yàn)概率密度是高斯混合形式［14］

式中：Kn－1為n－1時(shí)刻的時(shí)頻分量個(gè)數(shù)；和分別為分量k在n－1時(shí)刻的狀態(tài)均值和方差；為分量k的權(quán)重（表示該分量為真實(shí)分量的概率）；Ν（x；a，b）為隨機(jī)矢量x服從均值為a、方差為b的正態(tài)分布。n時(shí)刻已存在分量的預(yù)測(cè)概率密度υS，n｜n－1（x）為

式中：pS，n為時(shí)頻分量存活下來(lái)的概率［14］，分量k的預(yù)測(cè)均值和預(yù)測(cè)方差可通過(guò) Kalman預(yù)測(cè)得到。注意，n－1時(shí)刻不屬于任何時(shí)頻分量的測(cè)量值將被假設(shè)為n時(shí)刻的新生時(shí)頻分量。假設(shè)時(shí)頻分量的預(yù)測(cè)概率密度也是高斯混合形式：

式中：Kγ，n為n時(shí)刻假設(shè)時(shí)頻分量的個(gè)數(shù)。假設(shè)時(shí)頻分量k的初始權(quán)重ωkγ，n的賦值方式需要據(jù)應(yīng)用場(chǎng)景來(lái)確定。在時(shí)頻分析中，考慮到測(cè)量值的瞬時(shí)幅度越大，其屬于真實(shí)時(shí)頻分量的概率就越大，而在信噪比高的環(huán)境下，雜波較少，那么測(cè)量值來(lái)自于真實(shí)時(shí)頻分量的概率也越大。綜合上述兩個(gè)因素，本文給出權(quán)重ωkγ，n的初始值為

2.2.2 算法更新

獲得式（1）的觀測(cè)后，對(duì)N 個(gè)觀測(cè)值進(jìn)行時(shí)頻變換，例如STFT，時(shí)頻變換后的結(jié)果通過(guò)式（3）獲得n時(shí)刻測(cè)量值后，需要對(duì)權(quán)重進(jìn)行相應(yīng)的更新。對(duì)于每個(gè)測(cè)量，權(quán)重）根據(jù)貝葉斯公式更新為［14］

式中：pD，n為檢測(cè)到時(shí)頻分量的概率；）為雜波在上的概率密度，一般取常數(shù)1／V，其中V是測(cè)量空間的大?。环謩e為來(lái)自時(shí)頻分量k測(cè)量的預(yù)測(cè)均值和方差；Kn｜n－1＝Kn－1＋Kγ，n為n 時(shí)刻預(yù)測(cè)的時(shí)頻分量個(gè)數(shù)。當(dāng)假設(shè)時(shí)頻分量的頻率和幅度相互獨(dú)立時(shí)，式（25）中的似然函數(shù)可以寫為

事實(shí)上，測(cè)量zjn只可能來(lái)自某一個(gè)時(shí)頻分量或者來(lái)自雜波。這樣，為了降低計(jì)算復(fù)雜度，本文根據(jù)權(quán)重的大小將測(cè)量與時(shí)頻分量進(jìn)行關(guān)聯(lián)式（27）表明：當(dāng)n時(shí)刻所有時(shí)頻分量的權(quán)重都更新后，其中沒(méi)有與目標(biāo)關(guān)聯(lián)過(guò)的測(cè)量值，將假設(shè)為下一時(shí)刻的新生時(shí)頻分量。然而由于雜波的存在，必須要有一個(gè)判別準(zhǔn)則來(lái)決定假設(shè)時(shí)頻分量是否來(lái)自真實(shí)的新時(shí)頻分量，還是來(lái)自雜波。由于雜波能量較為分散，而時(shí)頻分量的能量相對(duì)較集中，而據(jù)式（25）可知，真實(shí)時(shí)頻分量的權(quán)重將接近于1，而虛假時(shí)頻分量（雜波）的更新權(quán)重將接近于0，應(yīng)該舍去。這樣本文的判別準(zhǔn)則如下：

1）若ωkn＞0.5，該假設(shè)分量為真實(shí)分量。

2）若ωkn＞1／V，該分量為假設(shè)分量。

3）否則，該假設(shè)分量為雜波（該準(zhǔn)則也同樣適用于判斷時(shí)頻分量的消失）。

3 時(shí)頻分析的多目標(biāo)追蹤算法

步驟1 對(duì)存在的時(shí)頻分量k進(jìn)行Kalman預(yù)測(cè)，其中k＝1，2，…，Kn－1。

式中：Q為過(guò)程噪聲協(xié)方差矩陣。

步驟2 由式（27）可知，對(duì)在n－1時(shí)刻判別為來(lái)自雜波的Kγ，n個(gè)測(cè)量視為n時(shí)刻的假設(shè)時(shí)頻分量，并據(jù)式（24）賦予相應(yīng)的初始權(quán)重。

步驟3 對(duì)存在時(shí)頻分量或假設(shè)時(shí)頻分量k的測(cè)量進(jìn)行預(yù)測(cè)，其中k＝1，2，…，Kn－1＋ Kγ，n－1。

式中：R為測(cè)量噪聲方差。

步驟4 對(duì)于給定的zjn，其中j＝1，2，…，Jn，式（27）將該測(cè)量值與時(shí)頻分量相關(guān)聯(lián)，并獲得該分量權(quán)重的更新值；更新時(shí)頻分量的狀態(tài)方差。

1）若該測(cè)量由某一時(shí)頻分量k產(chǎn)生，更新該時(shí)頻分量的狀態(tài)均值mkn｜n

步驟5 根據(jù)判別準(zhǔn)則，刪除權(quán)重較小的分量，接受權(quán)重較大的假設(shè)時(shí)頻分量為真實(shí)分量。

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

為了使仿真實(shí)驗(yàn)更具說(shuō)服力，本文采用了文獻(xiàn)［5，12］的仿真實(shí)驗(yàn)思路來(lái)設(shè)計(jì)動(dòng)態(tài)時(shí)頻譜。仿真實(shí)驗(yàn)分為3部分：① 仿真中的多分量線性調(diào)頻信號(hào)來(lái)自文獻(xiàn)［5］；② 仿真中包含微弱信號(hào)和非線性調(diào)頻信號(hào)的情況參考文獻(xiàn)［12］；③ 主要考慮載頻差異較小信號(hào)的時(shí)頻分析，以借此來(lái)說(shuō)明所提算法的高分辨率性能。

4.1 多分量線性調(diào)頻信號(hào)

第1組仿真信號(hào)包含10個(gè)線性調(diào)頻分量，形式為

式中：αk和βk分別為第k個(gè)信號(hào)分量的初始頻率和調(diào)頻斜率；Ik（t）為第k個(gè)信號(hào)分量是否存在的指示函數(shù)，其形式為

加性復(fù)高斯白噪聲η（t）的方差為σ2η，信噪比（SNR）定義為［21］

仿真信號(hào)長(zhǎng)度為100s，采樣頻率為512Hz，信噪比為－8dB。圖2分別給出了信號(hào)進(jìn)行STFT、Capon 變 換［12］和迭 代 自 適 應(yīng) （Iterative Adaptive Approach，IAA）譜估計(jì)（見(jiàn)附錄 A）的時(shí)頻分布圖，其中STFT、Capon和IAA均采用矩形窗，窗長(zhǎng)N＝128，窗的移動(dòng)步長(zhǎng)為N。對(duì)于圖2的含噪時(shí)頻圖，分別采用ADPM-RBPF算法［5］、基于 Capon的峰值檢測(cè)法［12］和本文算法進(jìn)行分析。本文算法中狀態(tài)空間模型參數(shù)為L(zhǎng)＝2，M ＝3，σa＝10－2，σfm＝10－4，其中m＝1，2，3。時(shí)頻分量存活下來(lái)的概率pS，n＝0.9，時(shí)頻分量被檢測(cè)到的概率pD，n＝0.9，S型曲線的陡峭程度參數(shù)δ＝20，權(quán)重因子ρ1＝2ρ2。

圖2 線性調(diào)頻信號(hào)時(shí)頻圖Fig.2 Time-frequency spectra of chirp signal

從圖3（a）可以看到，在低信噪比下，ADPMRBPF算法不能很好地去除雜波；從圖3（b）可以看到，基于Capon的峰值檢測(cè)法可以更好地去除雜波，但是該算法必須依靠更多的譜來(lái)判斷真實(shí)的信號(hào)分量和雜波［12］，因此實(shí)時(shí)性差。在相同的窗長(zhǎng)下，本文算法的實(shí)時(shí)性與ADPM-RBPF相當(dāng)，均優(yōu)于基于Capon的峰值檢測(cè)法。而圖3（c）和圖3（d）則表明：本文算法明顯優(yōu)于ADPM-RBPF算法，當(dāng)短時(shí)譜是利用IAA估計(jì)得到時(shí)，本文算法效果最佳。

單時(shí)頻分量背景下的均方根誤差已不再適用于對(duì)本文算法的性能評(píng)價(jià)，因?yàn)樾枰紤]狀態(tài)估計(jì)值和真實(shí)值之間的對(duì)應(yīng)關(guān)系。為此，本文利用OSPA（Optimal SubPattern Assignment）距離［22］來(lái)評(píng)價(jià)多分量估計(jì)的誤差。對(duì)于瞬時(shí)頻率真實(shí)值珚F＝［珚f1，珚f2，…，珚f珡K］和估計(jì)值＾F＝［＾f1，＾f2，…，＾fK＾］，當(dāng)珡K≤＾K時(shí)，OSPA距離定義為［22］

式中：

式中：參數(shù)c為目標(biāo)狀態(tài)估計(jì)誤差的閾值，同時(shí)可調(diào)節(jié)集合勢(shì)的估計(jì)誤差比重。本文采用p＝2階，c＝0.1的評(píng)價(jià)指標(biāo)。

圖4給出了3種算法時(shí)頻譜估計(jì)的平均OSPA曲線，每種算法均進(jìn)行100次蒙特卡羅仿真。通過(guò)圖4可以看到，本文算法瞬時(shí)頻率的估計(jì)精度優(yōu)于ADPM-RBPF算法和基于Capon的峰值檢測(cè)法，在低信噪比（SNR＜－6dB）下，本文算法優(yōu)勢(shì)更加明顯。這是因?yàn)楸疚乃惴檎{(diào)頻函數(shù)建立一個(gè)可靠的狀態(tài)轉(zhuǎn)移模型（18），這個(gè)新模型使得在低信噪比下，本文算法能依據(jù)該模型魯棒地跟蹤時(shí)頻分量。需要強(qiáng)調(diào)的是，基于IAA短時(shí)譜的本文算法估計(jì)性能最佳，這是因?yàn)镮AA算法優(yōu)于STFT算法。

圖3 線性調(diào)頻信號(hào)的時(shí)頻重構(gòu)圖Fig.3 Time-frequency reconstruction spectra of chirp signal

圖4 IF估計(jì)誤差Fig.4 IF estimation error

4.2 復(fù)雜多分量時(shí)頻信號(hào)

實(shí)際中，信號(hào)分量幅度大小往往不同，調(diào)頻方式多樣。為了驗(yàn)證本文算法在各種復(fù)雜情況下的魯棒性和有效性，本文將非線性調(diào)頻特性、信號(hào)分量的動(dòng)態(tài)出現(xiàn)和消失以及微弱信號(hào)分量都結(jié)合起來(lái)進(jìn)行分析。第2組仿真信號(hào)形式為

仿真信號(hào)長(zhǎng)度為100s，采樣頻率為512Hz，信噪比為0dB。圖5給出了信號(hào)進(jìn)行STFT變換、Capon變換和IAA譜估計(jì)的時(shí)頻分布圖。從圖中可以看出，信號(hào)分量k＝4較強(qiáng)，而信號(hào)分量k＝3非常微弱，其局部信噪比相當(dāng)于－10.46dB，檢測(cè)難度變大。

圖6給出了本文算法與對(duì)比算法的時(shí)頻重構(gòu)圖，可以發(fā)現(xiàn)，ADPM-RBPF算法（見(jiàn)圖6（a））不能完整重構(gòu)出信號(hào)分量k＝3，這是因?yàn)樵撍惴ㄐ韪鶕?jù)功率譜大小來(lái)進(jìn)行隨機(jī)采樣獲得測(cè)量值，從而導(dǎo)致微弱信號(hào)分量被檢測(cè)到概率很小?；贑apon的峰值檢測(cè)法（見(jiàn)圖6（b））在微弱信號(hào)以及調(diào)頻斜率較大區(qū)域均出現(xiàn)漏檢情況，這是因?yàn)樵撍惴ú捎玫募僭O(shè)檢驗(yàn)在信號(hào)能量較低區(qū)域（低信噪比區(qū)域）漏報(bào)率高的緣故。本文算法則可利用建立的預(yù)測(cè)模型式（18）來(lái)預(yù)測(cè)非線性調(diào)頻分量在下一個(gè)采樣時(shí)刻可能的近似值，同時(shí)也考慮到了時(shí)頻分量未被檢測(cè)到的概率（見(jiàn)式（25）），因此不論基于STFT或IAA短時(shí)譜，均可以完整重構(gòu)出全部4個(gè)信號(hào)分量。

圖7給出了ADPM-RBPF算法、基于Capon的峰值檢測(cè)法，以及本文算法估計(jì)出的信號(hào)分量數(shù)目，每種算法均進(jìn)行100次蒙特卡羅仿真。當(dāng)微弱信號(hào)分量存在時(shí)，ADPM-RBPF算法估計(jì)誤差較大；而當(dāng)調(diào)頻斜率較大或微弱信號(hào)分量存在時(shí)，基于Capon的峰值檢測(cè)法對(duì)于信號(hào)分量數(shù)目的估計(jì)誤差較大。本文算法分析結(jié)果較為理想，這是因?yàn)楸疚乃惴紤]到了目標(biāo)未被檢測(cè)到的情況。

4.3 載頻差異較小的時(shí)頻信號(hào)

為了獲得較高的時(shí)間分辨率，必須減小時(shí)間窗的長(zhǎng)度。然而，對(duì)于載頻差異較小且變化劇烈的時(shí)頻信號(hào)，STFT很難同時(shí)滿足時(shí)間和頻率分辨率的要求，而IAA譜估計(jì)算法具有高時(shí)頻分辨率和魯棒性，因而適合載頻差異較小的時(shí)頻分析問(wèn)題。第3組仿真中，時(shí)頻信號(hào)包含2組載頻差異較小的線性和非線性調(diào)頻信號(hào)為

圖6 復(fù)雜信號(hào)時(shí)頻重構(gòu)圖Fig.6 Time-frequency reconstruction spectra of complicated signal

圖7 復(fù)雜信號(hào)分量個(gè)數(shù)估計(jì)Fig.7 Components estimation of complicated signal

仿真信號(hào)最小載頻差異約為0.008（歸一化），信號(hào)長(zhǎng)度為100s，采樣頻率為512Hz，信噪比為0dB。圖8給出了信號(hào)進(jìn)行STFT變換、Capon變換和IAA譜估計(jì)的時(shí)頻分布圖，3種算法均采用矩形窗，窗長(zhǎng)N＝128，窗的移動(dòng)步長(zhǎng)為N?？梢钥吹?，在信號(hào)載頻差異較小的區(qū)域，STFT和Capon下信號(hào)分量均存在嚴(yán)重的重疊。

從圖9（a）～圖9（c）中可以看出，當(dāng)信號(hào)分量間載頻差異較小時(shí)，ADPM-RBPF算法、基于Capon的峰值檢測(cè)法以及基于STFT的本文算法均無(wú)法分析出兩個(gè)載頻差異較小的信號(hào)分量。相反，如圖9（d）所示，基于IAA時(shí)頻譜的本文方法即使在載頻差異較小區(qū)域也可獲得較為理想的分析結(jié)果。圖10給出了每種算法都進(jìn)行100次蒙特卡羅仿真的信號(hào)分量數(shù)目的估計(jì)結(jié)果。可以看到，只有基于IAA時(shí)頻譜的本文算法具有較為理想的估計(jì)效果。

圖8 載頻差異較小信號(hào)時(shí)頻圖Fig.8 Time-frequency spectra in close modes

圖9 載頻差異較小信號(hào)時(shí)頻重構(gòu)圖Fig.9Time-frequency reconstructionspectrainclosemodes

圖10 載頻差異較小信號(hào)分量個(gè)數(shù)估計(jì)Fig.10 Components estimation of signal in close modes

5 結(jié) 論

針對(duì)由多個(gè)信號(hào)分量構(gòu)成的時(shí)頻信號(hào)且其分量動(dòng)態(tài)出現(xiàn)和消失的時(shí)頻譜分析問(wèn)題，本文提出了一種分析、探測(cè)和跟蹤動(dòng)態(tài)時(shí)頻譜的隨機(jī)有限集法。

1）該算法首先將觀測(cè)到的信號(hào)進(jìn)行時(shí)頻變換，進(jìn)而視每個(gè)時(shí)刻的時(shí)頻譜幅度的局部極值為測(cè)量。分析表明：每個(gè)時(shí)刻測(cè)量的集合是一個(gè)隨機(jī)有限集。

2）基于提出的狀態(tài)和測(cè)量隨機(jī)有限集模型，本文給出了利用高斯混合概率假設(shè)密度濾波器，并結(jié)合多項(xiàng)式預(yù)測(cè)模型以及提出的初始權(quán)重賦值算法，來(lái)進(jìn)行動(dòng)態(tài)時(shí)頻譜的分析、探測(cè)和跟蹤的算法。

3）仿真結(jié)果均證明：所提算法在分析精確度、魯棒性、對(duì)微弱時(shí)頻譜分量的探測(cè)能力、以及對(duì)載頻差異的分析能力等諸方面均優(yōu)于文獻(xiàn)中報(bào)道的算法。

式中：y＝為時(shí)間窗內(nèi)的測(cè)量數(shù)據(jù)，N為時(shí)間窗的寬度；I為離散化頻率 格 點(diǎn) 的 個(gè) 數(shù)；為離散化頻點(diǎn)τ對(duì)應(yīng)的頻率矢量，tn為采樣時(shí)刻；eτ為離散化頻點(diǎn)τ對(duì)應(yīng)的幅值；ε為均值為0且與信號(hào)相互獨(dú)立的高斯白噪聲。

對(duì)式（A1）所描述的信號(hào)模型，IAA算法的時(shí)頻譜估計(jì)結(jié)果為

式中：IAAy（n，τ）為n時(shí)刻IAA算法估計(jì)得到的譜；WJτ為第J次迭代中針對(duì)頻點(diǎn)τ的加權(quán)矩陣。IAA算法采用迭代技術(shù)，即每次迭代開(kāi)始時(shí)都利用了上一次譜估計(jì)結(jié)果來(lái)構(gòu)建加權(quán)矩陣

式中：為第j次迭代中針對(duì)頻點(diǎn)τ的加權(quán)矩陣；為第j次迭代開(kāi)始時(shí)針對(duì)頻點(diǎn)τ估計(jì)的信號(hào) 協(xié)方 差 矩 陣 ；為 第j－1次 迭 代 的 譜估計(jì)結(jié)果，其中初始譜估計(jì)結(jié)果由周期圖法獲得［23］。