基于多重高阶同步压缩的高动态信号捕获技术

程凌峰，张英健，倪淑燕

(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方法,这也是后续要继续改进的方向。