基于多重高階同步壓縮的高動態信號捕獲技術

程凌峰，張英健，倪淑燕

(1.航天工程大學研究生院，北京 101416；2.航天工程大學 電子與光學工程系，北京 101416)

0 引言

低軌衛星軌道高度低,具有傳輸時延小、路徑損耗小的優點,近年來,Starlink等低軌衛星星座快速發展,已經成為一種重要的通信方式。但是對于飛機、導彈等用戶,一方面衛星軌道低,處于高速運動;另一方面自身速度很快,雙方處于高速的相對運動。在這種情況下接收信號會受到較強的多普勒效應的影響[1],體現為較大的多普勒頻偏及高階頻偏變化率[2-3]。當接收信號時間內頻偏變化超出FFT分辨率的情況下,傳統捕獲方法難以對信號進行有效的估計[4-6]。

高動態接收機在設計幀結構時,通常在每一幀的最前端添加導頻序列用于捕獲,該序列在經過高動態信道之后變為調頻信號。目前針對調頻信號的參數估計方法基本可分為2類:參數域類方法和時頻分析法。參數域類方法比較常用的有調頻率逼近法[7]、降階類參數估計算法[8]。調頻率逼近法對線性調頻信號有較好的估計性能,但是存在計算精度與計算復雜度之間的矛盾且無法對高階調頻信號進行參數估計。降階類參數估計算法比較適用于非線性調頻信號,但是存在高階參數估計量對低階參數估計時誤差積累的問題,同時可能存在偽峰,抗噪性較差。時頻分析法將時間表示和頻率表示聯合起來[9],是分析非平穩信號的理想工具,思路是對信號進行時頻分析,獲得時頻分布,提取時頻圖脊線即瞬時頻率(Instantaneous Frequency,IF),實現對頻率的粗估。文獻[10]提出用分數階傅里葉變換(FRFT)這種特殊的時頻分析方法對線性調頻信號進行估計,對于高階的調頻信號估計,過去通常使用傳統的短時傅里葉變換(STFT)和小波變換,但是受到海森堡測不準原理的限制[11],其時頻分布是模糊的,不利于脊線的提取。為了提高傳統方法的能力,最近幾十年提出了時頻重排方法(Reassignment Method,RM)[12]、同步壓縮(Synchrosqueezing Transform,SST)方法[13]等方法。RM雖然時頻聚集能力較好,但是無法實現信號重構;SST方法對快時變信號的聚集能力仍然較差,存在頻率估計誤差逐步增大的問題。針對SST的誤差逐漸增大和時域能量溢出現象[14],基于更精確的瞬時頻率估計的二階SST[15]和擴展形式的高階同步壓縮(High-Order Synchrosqueezing Transform,HOSST/SSTN)[16]陸續被提出。相比于傳統的幾種時頻分析工具,HOSST/SSTN具有更集中的時頻能量聚集能力和更高的分辨率,已經應用于機械系統的檢測等領域,但是存在對噪聲魯棒性較差、低信噪比下時頻圖發散嚴重的問題。根據本文的信號模型,后續對比實驗時,采用三階同步壓縮變換(Third-order Synchrosqueezing Transform,SST3)進行比較。

針對以上問題,本文提出一種基于多重高階同步壓縮變換(Multi-High-Order Synchrosqueezing Transform,MHOSST/MSSTN)的時頻分析法對信號進行載波捕獲,將HOSST/SSTN和多重同步壓縮(Multisynchrosqueezing Transform,MSST)相結合,對接收信號進行高分辨率的時頻分析,并針對常用的脊線提取能量泛函數法由于噪聲導致錯誤的提取點,提出將時頻面分成多個部分進行提取的方法,向前向后分別提取脊線,提高低信噪比下信號估計能力。

1 信號模型

進入載波同步過程的起始信號是通過傳輸信道傳輸,經過接收機自動增益控制、下變頻和匹配濾波后的含有噪聲的數字信號。

接收機的信號模型可表示為:

(1)

式中:s(k)為調制數據,k為整數;gT(t)為信號成形濾波器的沖激響應,T為符號周期,Δf、Δθ為接收機接收到的信號載波頻率、載波相位與本地載波頻率、相位的差值,w(nT)是均值為0的高斯白噪聲。

導頻序列通常為無數據變化的常數序列:

s(k)=1,k=0,1,…,P-1,

(2)

式中:P為導頻長。

暫不考慮成形濾波器的拖尾,經過下變頻和濾波后的基帶信號可表示為:

(3)

對于低軌道衛星,受徑向相對速度和相對速度變化率等影響,可將時變基帶信號相位近似高階泰勒級數展開:

(4)

式中:f0、f1、f2分別為頻偏、一階頻偏變化率和二階頻偏變化率。考慮到導頻數據較短的情況,忽略了更高階頻偏變化率分量造成的影響,從而將載波捕獲的問題轉換為非線性調頻信號的參數估計問題。

2 基于MSSTN的高動態信號時頻分析

非平穩信號的瞬時頻率為瞬時相位對時間的一階導數,理想情況下時頻能量只分布在瞬時頻率脊線上。由于快時變信號變化規律的復雜性,目前的時頻分析方法存在瞬時頻率估計精度有限、計算量大、噪聲魯棒性差的缺點,針對時頻分析后處理方法的局限性,本文提出了MSSTN這樣一種新的時頻分析方法,提出思路如圖1所示,后續仿真實驗對比足以證明所提方法的優越性。

圖1 MSSTN方法提出思路Fig.1 Source of ideas for MSSTN

2.1 SST原理

SST是基于STFT進行后處理得到的。對于信號r(t),STFT可表示為:

(5)

式中:g(t)為窗函數,r(t)為輸入信號。受海森堡測不準原理的制約,STFT無法同時實現時頻的高分辨率。SST方法實際上是對STFT在頻率方向進行能量重新分配的一種方法,如圖2所示,主要用于純諧波信號分析。其過程可以表示為[17]:

(6)

圖2 SST方法時頻能量分配示意Fig.2 Time-frequency energy distribution in SST method

接收信號的SST表達式為:

(7)

2.2 高動態信號的SST

對于研究的高動態信號,由于STFT的窗函數較短,二階頻偏變化率影響較小,因此可以將其所加窗內信號近似為線性調頻信號[18],即:

r(τ)=Aej[φ(t)+2πf0(τ-t)+πf1(τ-t)2]。

(10)

定義窗函數為高斯窗g(t)=e-0.5t2,則:

(11)

(12)

因此將式(12)兩邊對時間求偏導可得:

(13)

因為式(13)是復數形式,不能直接用于計算,取實部可得:

(14)

2.3 HOSST

對于SST方法的缺點,基于更為準確的瞬時頻率估計,二階同步壓縮變換(SST2)[17]被提出,首先定義一個二階局部調制因子,然后用來計算新的瞬時頻率估計。

對于給定的信號,二階局部調制系數為:

(15)

式中:

(16)

可得:

(17)

式中:Gg″(t,f)、Gg′(t,f)、Gtg(t,f)、Gtg′(t,f)分別表示窗函數為g″(t)、g′(t)、tg(t)、tg′(t)時的STFT的計算結果[18]。

該調制系數為SST瞬時頻率估計對時間的一階導數,二階瞬時頻率估計可表示為:

(18)

SST2已經被證明能夠提高線性調頻信號的時頻分布的效果[19-21]。針對更為復雜的快速變化頻率信號,根據SST2表達式遞推進一步提出了具有更為通用形式的HOSST方法。

(19)

定義y1(t,f)和xk,1(t,f)為:

(20)

式中:Gtk-1g(t,f)表示窗函數為tk-1g(t)的STFT。

從而可以將yj(t,f)和xk,j(t,f)用迭代的形式表示為:

(22)

HOSST的可表示為:

(23)

利用HOSST的方法,可以獲得更為準確的瞬時頻率估計,從而有效抑制時頻模糊現象。根據本文實際需求,研究信號為具有多普勒頻偏及一二階變化率的高動態信號,將其建模為三階非線性調頻信號。為了估計模型匹配,本文采用三階同步壓縮方法,即N=3。

2.4 MSST

區別于HOSST,對于快時變信號,MSST采用迭代的方法逐步接近瞬時頻率。

(24)

進一步可得:

(25)

圖3 MSST頻率壓縮過程示意Fig.3 MSST frequency compression process

(26)

(27)

雖然MSST中采用的固定點迭代原則能夠增強時頻分析能力,但理論分析已經證明了MSST是一種有偏瞬時頻率估計器[22],MSST計算的時頻表示結果保留了未重新分配能量的時頻點,時頻分布向理想分布收斂的速度有待提升。

2.5 MSSTN

類比MSST,式(23)中的SSTN表達式可寫為:

(28)

經過二重壓縮后的時頻分析結果可表示為:

(29)

以此類推,經過N次壓縮后的MSSTN的時頻分析結果可表示為:

本文考慮的場景主要涉及多重三階同步壓縮變換(MSST3)方法,具體實現的方法如算法1所示。

算法1 MSST3輸入:x為二階非線性調頻信號Gg=STFT(x),Gg′=STFTg′(x),Gg″(x)=STFTg″(x)Gtg=STFTtg(x),Gtg′=STFTtg′(x),Gt2g=STFTt2g(x)f～=f-1j2πGg′Gg()τ2=GtgGgτ3=Gt2gGgXk,j=GgGtkg-Gtj-1gGtk-j+1gW2=1j2π((Gg)2+GgGtg′-GtgGg′)W3=?fW2=1j2π(2GgGtg+GgGt2g-Gt2gGg′)q～[3,3]=W3X2,2-W2X3,3X4,3X2,2-X3,2X3,3q～[2,3]=W2X2,2-q～[3,3]X3,2X2,2f～[3]=f～+q～[2,3]τ2+q～[3,3]τ3

引入2個矩陣f和t,f和t分別對應于時頻分析結果的時間頻率集合。在進行多重壓縮時,頻率點集合表示為Sf,時間點集合表示為St,Sf和St的頻率時間點個數分別為Nf和Nt,重分配的次數為N。此外,對于W3=?fW2的求解需要用到:

?fGtk-1g=-j2πGtkg。

(31)

通過對上述變量和公式的確定,MSST3的具體實現算法見算法1。

3 仿真驗證

仿真條件:接收信號為二階非線性調頻信號,信號長度N=10 000,采樣率為20 kHz,頻率范圍為±10 kHz,調頻率范圍為±2 kHz/s,二階調頻率范圍為±5 kHz/s2,頻偏、一階調頻率和二階調頻率在取值范圍內隨機取值,在計算頻率估計均方根誤差和捕獲概率時采用500次蒙特卡洛仿真,比較STFT、RM、SST、SST3、MSST和MSST3這幾種時頻分析方法在高動態非線性調頻信號參數估計中的性能。

3.1 迭代次數的確定

作為一種迭代算法,需要確定MSST和MSST3迭代終止的條件。本文采用時頻分析瑞利熵和重構信號的信噪比變化來確定合適的迭代次數。

時頻聚集性是時頻分析方法性能優劣的重要指標之一,通過將時頻分析和信息熵結合,提出利用瑞利熵[23]來表征時頻分析方法性能,α階瑞利熵的定義為:

(32)

Ryα越小,時頻聚集程度越高。瑞利熵可以確保MSST時頻聚集到合適的水平,當迭代過程中連續時頻分析結果的瑞利熵沒有明顯的變化時,迭代過程可以中止。

對時頻分析結果進行重構,比較不同迭代次數重構后信號的信噪比,相鄰2次信噪比不再明顯變化時,可中止迭代[22]。圖4為迭代次數對瑞利熵和重構信號信噪比的影響,輸入信噪比為5 dB。可以看出,迭代次數越多,瑞利熵越小,重構信號的信噪比越高,瑞利熵在迭代次數為4后趨于穩定,信噪比在迭代3次后提升微弱。考慮到瑞利熵在迭代3次和4次差別不大和出于減少運算量的目的,本文擬采用迭代3次的多重壓縮方法。同時由圖4可以比較STFT、SST、MSST、SST3和MSST3這5種具有信號重構能力方法的重構性能,STFT、SST這2種方法重構時會造成信號信噪比的損失,SST3重構后對信號的信噪比有一定提升,經過多重壓縮后的MSST、MSST3在信噪比穩定后分別可提升約9、12 dB,這種能力可以考慮用于信號的增強。

(a)瑞利熵

3.2 時頻聚集能力

時頻聚集能力是準確估計瞬時頻率、提高捕獲精度的重要保證。假設信號信噪比為5 dB,比較STFT、SST、RM、SST3、MSST和MSST3這6種方法的時頻二維能量分布圖。將圖5中所有圖的能量刻度設置相同,通過比較可以看出,STFT的能量聚集性最差,局部放大后幾乎難以分辨;由于噪聲的加入和快時變信號的影響,SST方法的分析結果變得模糊;RM方法是從時間和頻率2個方向重新分配時頻能量,得到比SST更為集中的時頻圖;SST3方法雖然能夠應對二階非線性調頻信號的快時變,但是自身的噪聲魯棒性較差,時頻分析圖也出現了發散的現象;同時可以發現多重壓縮方法確實可以提高時頻聚集能力,MSST和MSST3這2種方法都獲得了清晰的瞬時頻率脊線,但是受到SST有偏估計的影響,MSST方法在一些時刻有2個能量低的時頻點,在圖中時頻能量顯示為藍色,并且無法隨迭代次數的增加而匯成一個時頻點,在提取瞬時頻率脊線時容易導致提取錯誤,這就是MSST方法存在的未分配時頻能量點的缺點,不利于頻率的準確提取。

(a)STFT時頻圖

用瑞利熵對4種方法的性能進行定量分析,結果如圖6所示。可以看出,輸入信噪比為-10～15 dB時,多重壓縮算法的瑞利熵遠小于其他幾種方法,在信噪比-2 dB以下時,多重壓縮方法不足以彌補HDSST對噪聲魯棒性較差的問題,出現了MSST3瑞利熵大于MSST的現象,這對在相應的場景選擇合適的方法有指導意義。

圖6 不同時頻分析方法的瑞利熵隨信噪比變化情況Fig.6 Variation of Rayleigh entropy with signal to noise ratio in different time-frequency analysis methods

3.3 頻率捕獲精度

圖7為不同方法的頻率估計均方根誤差隨信噪比變化情況。多重壓縮方法的頻率估計誤差始終是最小的,在較高信噪比(8 dB)條件下,各種方法的均方根誤差變化趨于平穩,MSST3的估計誤差降到3 Hz以下,約為2.4 Hz,MSST的估計誤差穩定在3.3 Hz左右,降至其他方法誤差的1/3以下。在較低信噪比(-2 dB)條件下,MSST3的頻率估計性能還是略弱于MSST方法。同時可以發現其他4種方法的均方根誤差不是完全按照時頻聚集能力的順序排列的,這說明作為STFT的后處理方法,不僅將有用信號的瞬時頻率時頻能量聚集,同時也會將噪聲聚集,在瞬時頻率估計的時候造成過大的誤差。很明顯,MSST和MSST3對信號的時頻聚集能力強于噪聲聚集能力,而SST、RM和SST3在低信噪比時瞬時頻率的時頻能量嚴重發散,無法準確地實現瞬時頻率估計,甚至相較于STFT的估計精度都出現下滑,這也是研究時頻分析方法改進的方向,將時頻聚集能力更多地發揮在對有用信號能量的積聚上。

圖7 頻率估計均方根誤差隨信噪比變化情況Fig.7 Variation of frequency estimation root mean square error with signal to noise ratio

3.4 捕獲概率與時間

捕獲概率通常定義為信號頻率估計誤差小于20 Hz的概率。圖8比較了不同算法的捕獲概率,在較高信噪比下,MSST3方法的捕獲概率穩定在99.5%左右,優于其他方法;在較低信噪比下MSST方法捕獲性能更優。同時也發現SST、RM和SST3這3種方法在高信噪比條件下捕獲概率均收斂到95%以上,而STFT的捕獲概率收斂到80%左右,證明了對STFT后處理的必要性。

圖8 不同方法捕獲概率隨信噪比變化情況Fig.8 Variation of acquisition probability with signal to noise ratio in different methods

在實際應用時,時頻分析方法的運行效率非常關鍵,決定了能否應用于實際工程中。假設時間點數和頻率點數分別為Nt和Nf,STFT的時間計算復雜度為ο(NtNflog(Nf),后處理方法的多重壓縮、求導及重新分配時頻能量等步驟的計算復雜度為低階項,可忽略。采用測試的計算機配置為i7-10750H 2.6 GHz,Windows 10系統,Matlab版本為R2020a。所需時間如表1所示,所有方法都在20 s內完成處理。可以看出,MSST3和SST3相對于其他方法,由于求導高階項步驟的增加,算法復雜度提高,運行時間增長。

表1 不同時頻分析方法所需計算時間Tab.1 Computing time required by different time-frequency analysis methods 單位:s

4 結論

本文主要針對大頻偏、快時變條件下的載波頻率捕獲問題進行研究,對基于時頻分析的捕獲方法進行研究,通過仿真對比,得到以下結論:

本文提出一種MSST3的時頻分析方法,可用于高動態載波捕獲,這種方法既考慮了SST3方法對二階非線性調頻信號這種快時變信號能夠準確估計的優點,又考慮了多重壓縮方法提高噪聲魯棒性的能力。仿真實驗表明,這種方法能夠有效應對多重同步壓縮方法多次迭代無法消除未分配能量點的問題,而且能夠提供更為集中的時頻分析結果;同時在對時頻分析后的信號重構時發現,信噪比提升約12 dB,有較大的的應用價值。在較高信噪比下,頻率估計誤差降到2.4 Hz左右,比MSST方法降低了0.9 Hz左右,約為STFT誤差的1/4,約為其他后處理方法的1/3。在頻率捕獲概率方面,捕獲概率為99.5%左右,比MSST方法提高約0.4%左右,也明顯優于其他方法。

作為SST3的多重壓縮方法,通過多重壓縮,極大地提高了噪聲魯棒性,但是在仿真時發現,較低信噪比下捕獲性能略遜于MSST方法,這也是后續要繼續改進的方向。