曹偉浩，姚直象，夏文杰，閆肅

(1.海軍工程大學(xué) 電子工程學(xué)院, 湖北 武漢 430033;2.國家海洋技術(shù)中心漳州基地籌建辦公室，福建 廈門 361000; 3.91388部隊(duì)，廣東 湛江 524000)

0 引言

線性調(diào)頻(LFM)信號是主動聲納應(yīng)用較廣泛的一種信號形式，由于水聲環(huán)境的復(fù)雜性及水中可利用的信號頻帶有限，對噪聲污染的LFM信號進(jìn)行瞬時頻率估計(jì)變得非常重要。LFM信號高精度的瞬時頻率估計(jì)，不僅在水下目標(biāo)探測、水聲通信多普勒補(bǔ)償?shù)确矫姘缪萘酥匾巧?，而且在雷達(dá)、通信等領(lǐng)域得到了廣泛應(yīng)用。此外，提高LFM信號參數(shù)估計(jì)的精度，對水下高精度目標(biāo)識別也具有重要意義。

短時傅里葉變換(STFT)[1-2]和維格納- 威利分布(WVD)[3-4]是兩種常用的對LFM信號進(jìn)行時頻分析的方法，但是STFT對于寬帶信號頻率估計(jì)效果并不理想，WVD是雙線性變換，時頻分析存在交叉項(xiàng)問題。分?jǐn)?shù)階傅里葉變換(FRFT)對LFM信號具有較好的能量聚集性，利用此變換對LFM信號進(jìn)行時頻分析的應(yīng)用前景廣闊[5-10]。

本文針對參數(shù)估計(jì)問題，利用聲納、雷達(dá)等領(lǐng)域常用的一種頻譜修正算法——Rife插值法，提出一種基于Rife插值短時FRFT(STFRFT)方法來估計(jì)LFM信號瞬時頻率，通過變權(quán)擬合(VWSF)估計(jì)信號的初始頻率和調(diào)頻率，進(jìn)一步推導(dǎo)了初始頻率和調(diào)頻率估計(jì)的克拉美- 羅限(CRLB)，分析了CRLB隨參數(shù)變化的規(guī)律。

1 時頻分析和參數(shù)估計(jì)方法原理

1.1 短時分?jǐn)?shù)階傅里葉變換

LFM信號作為水聲信號處理領(lǐng)域中一種常用的信號，其表達(dá)式如下：

x(t)=Aexp {-j(2πf0t+πμt2)}，

式中：A為信號幅度；f0和μ分別為信號初始頻率和調(diào)頻率；t為信號持續(xù)時間。實(shí)際處理中，信號采用離散形式，故給出離散FRFT的表達(dá)式[10]為

(1)

式中：Aα=exp {-jπsgn(sinα)/4+jα/2}/|sinα|1/2，α=(p-1)π/2，p為FRFT最優(yōu)階數(shù)，0<|p|<2；x為量綱歸一化后新的橫坐標(biāo)；Δx為量綱歸一化寬度；X1(·)為信號的傅里葉變換；N為采樣點(diǎn)數(shù)，|m0|≤N，|n0|≤N.

若對信號進(jìn)行加窗處理，則由(1)式可得第i個短時窗內(nèi)離散FRFT的表達(dá)式為

(2)

式中：i=1,2,…,K，K為窗移動的總次數(shù)，ΔK為移動步進(jìn)；m=0,1,2,…,M-1；h(·)為長度為M的短時窗；αi=(pi-1)π/2，pi為第i個窗內(nèi)FRFT的最優(yōu)階數(shù)。通?？紤]到估計(jì)精度與計(jì)算量等因素，移動步進(jìn)取窗長的0.1～0.3倍之間比較合適。

圖1 LFM信號歸一化功率譜密度

FRFT作為傅里葉變換的一種推廣，直觀上可以認(rèn)為是Chirp基分解，非常適用于處理Chirp信號。對一個窗內(nèi)的LFM信號分別作傅里葉變換和FRFT，得到歸一化功率譜密度如圖1所示。從圖1中可以明顯看出，經(jīng)過傅里葉變換后頻譜幅度在一定帶寬內(nèi)起伏變化，而經(jīng)過FRFT后頻譜幅度位于某一頻率，表明FRFT對LFM信號具有更好的能量聚集性。

1.2 基于STFRFT的LFM瞬時頻率估計(jì)

對信號進(jìn)行瞬時頻率估計(jì)時，在一定頻率范圍內(nèi)，短時窗內(nèi)的信號可近似認(rèn)為是單頻信號。STFT中通常取信號頻譜幅度的均值作為一個窗內(nèi)的頻譜幅度，而實(shí)際頻譜幅度高于此均值。因此，可以令第i個短時窗內(nèi)FRFT幅度最大值對應(yīng)的頻率為該窗內(nèi)LFM信號瞬時頻率fi的估計(jì)值i，

i=kiΔf，

由于FRFT滿足Parseval關(guān)系，可得

(3)

式中：Efi、Eti分別為信號頻域和時域能量；Ni為第i個短時窗內(nèi)時域采樣點(diǎn)數(shù)；Nαi為第i個短時窗內(nèi)分?jǐn)?shù)階域采樣點(diǎn)數(shù)，Nαi=Nfisinαi=μT2sinαi，Nfi為第i個短時窗內(nèi)頻域采樣點(diǎn)數(shù)，αi=piπ/2，pi為第i個短時窗內(nèi)FRFT的最優(yōu)階數(shù)，μ為調(diào)頻率，T為短時窗時長；Ai為第i個短時窗內(nèi)信號的頻譜幅度。由(3)式可以得到第i個短時窗內(nèi)信號經(jīng)過FRFT后的頻譜幅度為

式中：fs為采樣率。

對于離散頻譜，考慮柵欄效應(yīng)后的頻譜幅度為

(4)

當(dāng)信號頻率不是頻率分辨率的整數(shù)倍時，不可避免地會存在頻譜泄露情況，例如1個主瓣內(nèi)會存在2個峰值。采用較常用的Rife插值法進(jìn)行頻譜修正，該方法的頻差表示[11]為

δi=

分別采用STFT、插值STFT(ISTFT)、STFRFT和插值STFRFT(ISTFRFT)4種方法估計(jì)LFM信號的瞬時頻率，結(jié)果如圖2所示。由圖2可見，整體而言，STFRFT的方法較STFT方法更優(yōu)，對于STFRFT方法而言，ISTFRFT方法比STFRFT方法估計(jì)的頻率更加接近實(shí)際頻率。

圖2 LFM信號瞬時頻率估計(jì)

在得到LFM信號每一時刻的瞬時頻率后，利用曲線擬合方法估計(jì)其初始頻率和調(diào)頻率。最小二乘擬合(CLSF)是將觀測值賦予相同的權(quán)值，使觀測值與擬合值的誤差平方和達(dá)到最小。當(dāng)Rife插值估計(jì)的頻率存在異常點(diǎn)時，將每次觀測值的影響同等看待顯然不合理，應(yīng)該對誤差平方和較大的項(xiàng)賦予較小的權(quán)值，誤差平方和較小的項(xiàng)賦予較大的權(quán)值。因此采用VWSF方法估計(jì)信號的初始頻率和調(diào)頻率。

1.3 基于VWSF的LFM初始頻率和調(diào)頻率估計(jì)

式中：f′i=a1+a2t，a1、a2為2個待估計(jì)參量；Wi為第i個時刻估計(jì)瞬時頻率的加權(quán)值。分別求RW(f)關(guān)于a1、a2的偏導(dǎo)并令其偏導(dǎo)為0，可解得2個參量的值[12]如下：

以上兩個公式中的權(quán)值可由預(yù)先設(shè)定的精度控制實(shí)現(xiàn)。

分別利用CLSF和VWSF兩種方法估計(jì)的瞬時頻率結(jié)果如圖3所示。從圖3中可以看出，當(dāng)估計(jì)頻率出現(xiàn)異常點(diǎn)時，利用常規(guī)CLSF方法得到的曲線偏離實(shí)際頻率，而VWSF方法能夠很好地?cái)M合實(shí)際頻率。

圖3 LFM信號瞬時頻率擬合曲線

2 初始頻率和調(diào)頻率估計(jì)的CRLB研究

2.1 插值STFRFT頻率估計(jì)的CRLB

在一定條件下，任何估計(jì)量都存在一個均方誤差下限，即CRLB.對于未知確定參量f0和μ的均方誤差等于估計(jì)量的方差，該估計(jì)為無偏估計(jì)。因此，求得估計(jì)量的方差限即為其CRLB[13].

本文采用Rife插值法，參照文獻(xiàn)[14]中的方法，給出一個窗內(nèi)插值STFRFT方法估計(jì)頻率的均方誤差為

(5)

假設(shè)估計(jì)的LFM信號瞬時頻率為

f′i=f′0+μ′i+wi，

式中：f′0和μ′分別為估計(jì)的初始頻率和調(diào)頻率；wi為高斯白噪聲，其方差為σ2.則瞬時頻率的似然函數(shù)為

P(x|f′0,μ′)=

關(guān)于未知參量的Fisher信息矩陣為

式中：Pf′0f′0為函數(shù)P對f′0的2階偏導(dǎo)；Pf′0μ′、Pμ′f′0為函數(shù)P對f′0和μ′的2階偏導(dǎo)；Pμ′μ′為函數(shù)P對μ′的2階偏導(dǎo)。由此可得

(6)

2.2 初始頻率估計(jì)的CRLB分析

根據(jù)(6)式可以得到初始頻率估計(jì)的CRLB為

(7)

在輸入信噪比一定的情況下，取δ=0.5，初始頻率的CRLB與滑窗次數(shù)、FRFT的最優(yōu)階數(shù)有關(guān)。對(7)式的K求偏導(dǎo)，有

(8)

考慮實(shí)際情況，當(dāng)K大于1時(8)式恒小于0，表明初始頻率的估計(jì)精度隨滑窗次數(shù)增加而提高，即CRLB隨滑窗次數(shù)的增加而減小。同理，對α求偏導(dǎo)，有

(9)

由于階數(shù)p是以4為周期的參數(shù)，只需考察區(qū)間p∈[-2,2]即可，根據(jù)對稱性，取階數(shù)范圍為p∈[0,2]，則α∈[0,π]。當(dāng)α∈[0，π/2]時cosα>0，(9)式為正值；當(dāng)α∈[π/2,π]時cosα<0，(9)式為負(fù)值。因此初始頻率估計(jì)精度隨階數(shù)的增加先降低、后提高，即CRLB隨階數(shù)的增加先增大、后減小。當(dāng)階數(shù)為1時，初始頻率估計(jì)性能最低，此時CRLB達(dá)到最大。當(dāng)階數(shù)不變時，初始頻率估計(jì)性能與滑窗次數(shù)呈正比，但是滑窗次數(shù)的確定還需兼顧頻率分辨率、計(jì)算量等因素。因此，分析CRLB的變化規(guī)律對窗長等參數(shù)選取具有一定的指導(dǎo)意義。

初始頻率估計(jì)的CRLB仿真結(jié)果歸一化后如圖4所示。其中，圖4(a)的階數(shù)取值為[1，2]。由圖4(a)可以看出CRLB隨滑窗次數(shù)增加呈遞減趨勢。圖4(b)的階數(shù)取值為[0，2]，由圖4(b)可以看出CRLB隨階數(shù)的增加先增大、后減小，階數(shù)為1時估計(jì)性能最低。這一結(jié)論與理論分析結(jié)果一致。

2.3 調(diào)頻率估計(jì)的CRLB分析

根據(jù)(6)式可以得到調(diào)頻率估計(jì)的CRLB為

(10)

式中：τ=ΔK/fs，即短時窗移動步進(jìn)對應(yīng)時長。

在輸入信噪比一定的情況下，取δ=0.5，移動步進(jìn)取窗長的一半時，調(diào)頻率估計(jì)的CRLB與滑窗次數(shù)、信號時長和FRFT的最優(yōu)階數(shù)有關(guān)。調(diào)頻率的CRLB表達(dá)式為

(11)

式中：Tx為信號時長，為便于分析，取Tx=1 s，此時調(diào)頻率估計(jì)的CRLB與滑窗次數(shù)、FRFT的最優(yōu)階數(shù)有關(guān)，分析方法與初始頻率估計(jì)的CRLB分析相同。分別求(11)式關(guān)于K和α的偏導(dǎo)，有

(12)

(13)

當(dāng)K>1時，(12)式的值恒小于0，故調(diào)頻率的估計(jì)精度隨滑窗次數(shù)的增加而提高，即CRLB隨滑窗次數(shù)的增加而減小。(13)式的正負(fù)由cosα決定，變化規(guī)律與初始頻率估計(jì)的CRLB情況一致。

對調(diào)頻率估計(jì)的CRLB仿真結(jié)果歸一化后如圖5所示。同樣地，圖5(a)的階數(shù)取值為[1，2]。由圖5(a)可以明顯看出，調(diào)頻率估計(jì)的CRLB隨滑窗次數(shù)增加而減小。圖5(b)的階數(shù)取值為[0,2]。由圖5(b)可以看出，調(diào)頻率估計(jì)的CRLB隨階數(shù)的增加先增大、后減小，階數(shù)為1時估計(jì)性能最低。這一結(jié)論與理論分析結(jié)果一致。

圖5 歸一化調(diào)頻率估計(jì)的CRLB

3 實(shí)驗(yàn)驗(yàn)證

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

仿真信號采用LFM信號，初始頻率為600 Hz，調(diào)頻率為800 Hz/s，采樣率為4 096 Hz，脈寬為1 s，噪聲為高斯白噪聲，信噪比變化區(qū)間為-12～0 dB.經(jīng)過5 000次蒙特卡洛仿真，得到估計(jì)的調(diào)頻率均值及其均方根誤差如圖6所示；估計(jì)的初始頻率均值及其均方根誤差如圖7所示。圖6和圖7中，VWSF-ISTFT表示利用插值STFT估計(jì)信號瞬時頻率，然后通過VWSF估計(jì)參數(shù)；VWSF-ISTFRFT表示利用插值STFRFT估計(jì)信號瞬時頻率，然后通過VWSF估計(jì)參數(shù)。

圖6 調(diào)頻率估計(jì)

從圖6(a)可以看出，輸入信噪比為-11 dB時，VWSF-ISTFRFT方法估計(jì)的調(diào)頻率已經(jīng)接近真值，與VWSF-ISTFT方法相比，所要求的輸入信噪比降低了1 dB；相比于CLSF-ISTFRFT方法，VWSF-ISTFRFT方法所要求的信噪比降低了約3 dB.從圖6(b)可以看出，輸入信噪比為-11 dB時，VWSF-ISTFRFT方法估計(jì)的調(diào)頻率精度已經(jīng)接近CRLB;與VWSF-ISTFT方法相比，相同估計(jì)精度下，VWSF-ISTFRFT方法要求的輸入信噪比降低了1 dB；達(dá)到相同的估計(jì)精度時，VWSF方法同CLSF方法相比，要求的輸入信噪比降低了約3 dB.

圖7 初始頻率估計(jì)

從圖7(a)可以看出，輸入信噪比為-11 dB時，VWSF-ISTFRFT方法估計(jì)的初始頻率已經(jīng)接近真值，與VWSF-ISTFT方法相比信噪比降低了1 dB；相比于CLSF-ISTFRFT方法，VWSF-ISTFRFT方法所要求的信噪比降低了約3 dB.從圖7(b)可以看出，相同估計(jì)精度下，VWSF-ISTFRFT方法估計(jì)初始頻率所要求的輸入信噪比較VWSF-ISTFT方法降低了1 dB，VWSF方法所要求的輸入信噪比較CLSF方法降低了約3 dB.

針對被動聲納檢測到的脈沖信號參數(shù)估計(jì)問題，為了證明本文方法的普適性，在被動聲納工作頻率范圍內(nèi)隨機(jī)選取參數(shù)進(jìn)行仿真。根據(jù)文獻(xiàn)[15]可知，現(xiàn)有的被動聲納工作的大致頻率范圍為0.1～10 kHz，調(diào)頻信號通常作為窄帶信號應(yīng)用，脈寬由幾十毫秒到幾秒，在合理的參數(shù)取值范圍內(nèi)取3組參數(shù)進(jìn)行仿真。信號均采用LFM信號，噪聲為高斯白噪聲，信噪比變化范圍為-12～0 dB，其他仿真條件及得到的結(jié)果如下。

1)仿真條件1.采樣率為2 kHz，初始頻率為50 Hz，帶寬為80 Hz，脈寬為4 s，經(jīng)過5 000次蒙特卡洛仿真，得到調(diào)頻率和初始頻率的估計(jì)誤差如圖8和圖9所示。

圖8 調(diào)頻率估計(jì)誤差(仿真條件1)

圖9 初始頻率估計(jì)誤差(仿真條件1)

2)仿真條件2.采樣率為8 kHz，初始頻率為2 kHz，帶寬為100 Hz，脈寬為1 s，經(jīng)過5 000次蒙特卡洛仿真，得到調(diào)頻率和初始頻率的估計(jì)誤差如圖10和圖11所示。

圖10 調(diào)頻率估計(jì)誤差(仿真條件2)

圖11 初始頻率估計(jì)誤差(仿真條件2)

3)仿真條件3.采樣率為20 kHz，初始頻率為8 kHz，帶寬為400 Hz，脈寬為0.4 s，經(jīng)過5 000次蒙特卡洛仿真，得到調(diào)頻率和初始頻率的估計(jì)誤差如圖12和圖13所示。

圖12 調(diào)頻率估計(jì)誤差(仿真條件3)

圖13 初始頻率估計(jì)誤差(仿真條件3)

從上述3個不同條件的仿真結(jié)果中均可以看出，輸入信噪比為-12 dB時，VWSF-ISTFRFT方法仍有較高的估計(jì)精度，與基于STFT方法相比，VWSF-ISTFRFT方法的初始頻率和調(diào)頻率估計(jì)精度更高，估計(jì)性能更優(yōu)；達(dá)到相同的估計(jì)精度時，VWSF-ISTFRFT方法所要求的輸入信噪比更低。

3.2 海試驗(yàn)證

2017年12月在某地海試，發(fā)射信號為LFM信號，采樣率為64 kHz，脈寬為128 ms，歸一化后初始頻率為0.96，信噪比為-8 dB.分別采用VWSF-ISTFT、CLSF-ISTFRFT和VWSF-ISTFRFT方法3種方法估計(jì)的瞬時頻率如圖14所示。從圖14中可以看出，VWSF-ISTFRFT方法估計(jì)的瞬時頻率最接近實(shí)際頻率。將初始頻率和調(diào)頻率歸一化，估計(jì)的參數(shù)如表1所示。從表1中可以看出，VWSF-ISTFRFT方法估計(jì)的值最接近真實(shí)值，CLSF-ISTFRFT方法估計(jì)參數(shù)誤差最大。

圖14 LFM信號瞬時頻率估計(jì)

從海試數(shù)據(jù)的處理結(jié)果中可以看出，VWSF-ISTFRFT方法估計(jì)的瞬時頻率最接近實(shí)際頻率；從表1中數(shù)據(jù)可以看出，VWSF-ISTFRFT方法估計(jì)的參數(shù)精度最高，與CLSF方法相比，VWSF方法的參數(shù)估計(jì)精度由百分之一量級提高到千分之一量級，提高了一個數(shù)量級。

表1 估計(jì)參數(shù)及相對誤差

注：括號內(nèi)的值為參數(shù)估計(jì)相對誤差。

4 結(jié)論

本文針對低信噪比下LFM信號參數(shù)估計(jì)問題，利用Rife插值STFRFT方法估計(jì)LFM信號的瞬時頻率，通過VWSF估計(jì)初始頻率和調(diào)頻率，并推導(dǎo)分析了初始頻率和調(diào)頻率估計(jì)的CRLB影響因素。通過仿真分析和海試數(shù)據(jù)驗(yàn)證，得到主要結(jié)論如下：

1)插值STFRFT方法估計(jì)的LFM信號瞬時頻率較插值STFT精度更高，有著更優(yōu)的估計(jì)性能，估計(jì)精度相同時，對輸入信噪比的要求前者較后者更低。

2)利用瞬時頻率進(jìn)行參數(shù)估計(jì)時，VWSF方法可以消除插值異常點(diǎn)的影響，與CLSF方法相比，得到的參數(shù)估計(jì)精度更高。

3)初始頻率和調(diào)頻率估計(jì)的CRLB受滑窗次數(shù)和FRFT階數(shù)影響，隨滑窗次數(shù)的增加而減小，隨FRFT階數(shù)的增加先增大、后減小，階數(shù)為1時，CRLB達(dá)到最大值，此時參數(shù)的估計(jì)性能最低。