在越来越强调电子系统隐蔽攻击和硬杀伤功能的驱使下,无源探测技术[1]为空间目标隐蔽探测和精确定轨提供了十分重要的手段。基于非合作外辐射源的无源雷达系统通过估计外辐射源的直达波信号的调制参数,构建匹配滤波器,检测和分析目标反射辐射源发射的信号能量,从而实现对目标的定位和跟踪[2]。因此,准确快速地估计直达波的参数是无源雷达完成检测和跟踪任务的前提。
线性调频信号(chirp signal)是一种应用非常广泛的非平稳外辐射源信号,chirp信号的低截获概率特性也是其广泛应用在各个雷达体制中的重要原因。在以chirp雷达为外辐射源的无源雷达系统中,参考天线按一定的周期来截获直达波信号。在一个观测周期内,chirp信号往往以多分量的形式出现,分量间存在时频混叠,而且时间不同步。这种广义形式(占空比η≤100%)并不满足chirp信号的一般数学表达式(占空比η0=100%),这使得一些基于最大似然估计[3-4]和基于检测时频脊线[5-6]的方法无法直接应用于直达波信号的参数估计。尽管有学者提出在没有交叉干扰的时频图上检测时频线段的方法[7],但这些时频图大都因为平滑处理而无法表示信号的瞬时相位,因此估计精度较差,而且缺乏计算效率。
国内学者[8-9]提出了一种基于分数阶傅里叶域滤波的信号分离方法,其基本思想是根据chirp信号在某个分数阶傅里叶域的稀疏性和FrFT的无损可逆性,将chirp信号分离并恢复为单分量chirp信号,这一思想同样适用于占空比η≤100%的广义形式的chirp信号[10]。然而,即使工程上已存在分数阶傅里叶变换(Fractional Fourier Transform,FrFT)的快速算法[11],但在缺少关于调频斜率先验信息的条件下,则需要在所有变换阶内搜索,这对信号的实时处理带来了极大的压力。在工程背景下,非合作chirp雷达信号通常有多套固定的参数模板以完成跟踪和搜索任务,在不同的时间段内产生不同的频率参数组合。经过长期的工作积累,不难获得这些雷达信号的参数模板。因此,在该先验信息的帮助下,对非合作外辐射源雷达直达波的实时参数估计可以在有限的变换阶上实现。
本文将利用分数阶傅里叶域滤波方法实现无源雷达系统中直达波的参数估计。本文首先介绍了chirp雷达作为外辐射源时直达波信号的一般形式,然后介绍了基于分数阶傅里叶域滤波的多分量chirp信号的参数估计方法,最后通过仿真,分析了该方法处理非合作chirp雷达信号的性能及表现。
无源雷达[12-13]作为一种绿色环保、经济安全的被动探测手段近年来获得学术界的广泛关注。以chirp雷达为外辐射源的无源雷达系统实现目标的探测与跟踪的模型[14]如图 1所示,a,b,c分别表示非合作chirp雷达与目标的距离、目标与接收机之间的距离、非合作chirp雷达与接收机间的距离;φ表示目标与直达波间的方位角。一般情况下,c是已知的,φ可通过特定的DOA估计算法获得,因此,只需求得直达波与目标散射回波间的到达时间差ΔTOA=τa+τb-τc,即可求解a,b,c。
图1 无源雷达目标检测与跟踪模型
采样频率为Fs的离散chirp序列的数学表达式为
s(nΩ)=
n=1,…,N
(1)
式中,N为当前观测帧的长度,采样间隔观测周期W0=N/Fs,chirp序列的脉冲宽度Pw=W0,占空比η0=100%。
在工程应用中,以采样频率Fs采集到的chirp雷达信号可表示为
(2)
si(nΩ)=
n=1,…,N
(3)
式中,Nc表示当前观测帧内chirp分量的数目,v为背景噪声,si为广义形式的chirp序列,其脉冲宽度受rect(·)函数的约束,占空比信号参数Ai,τi,Pw,i,fi,0,ki分别表示si的幅度、时延、脉宽、中心频率与调频斜率。在工程应用中,只需估计参数Θ=,即可构造匹配滤波器对目标散射回波进行匹配滤波。
连续全脉冲chirp信号的p阶分数阶傅里叶变换Sp(μ)表达式为
Sp(μ)=FrFTp[s(t)]=
s(t)Ka(t,μ)dt
(4)
式中,FrFTp表示p阶FrFT变换算子,p为FrFT的变换阶,为信号在时频平面上的旋转角,Ka为p阶FrFT的变换核,表示为
Ka(t,μ)=
(5)
令将式(5)代入式(4)可得
Sp0(μ)=
exp(-jπμ2k)δ(csca0-f0)
(6)
当cota0=-k时,chirp信号在p0阶分数阶傅里叶域表现为冲击函数,谱峰位置为μ=μ0=f0sina0。对于占空比η<100%的chirp信号,谱峰估计仅能获得chirp信号在时频图上所在直线的调制参数,即
(7)
式中,表示时频脊线的中心频率。显然,当chirp分量的占空比ηi<100%时,只通过分数阶傅里叶变换无法获得多分量chirp信号的中心频率fi,0及脉冲宽度Pw,i,因此,我们可以借助如下所述的分数阶傅里叶域滤波的方法。
基于chirp分量在分数阶傅里叶域的稀疏特性及FrFT的可逆无损特性,信号分离是一种可行的方法。分数阶傅里叶域滤波实现信号分离的基本思想是在每个p0阶分数阶傅里叶域构造一个中心频率为μ0的带通滤波器,根据检测到的谱峰,对每个Sp0进行滤波,因为分数阶傅里叶变换核满足有
s(t)=Sp(μ)K-a(t,μ)dμ
(8)
因此,可通过逆分数阶傅里叶变换(inverse Fractional Fourier Transform,iFrFT)将滤波后的分数阶傅里叶域信号恢复为单分量时域信号,完成信号分离,该过程如图 2所示。
图2 分数阶傅里叶域滤波实现chirp信号分离
以步进Δp对chirp信号进行FrFT后得到(p,μ)的参数平面;通过谱峰检测得到峰值坐标(pi,μi),由式(7)可直接获得chirp分量的调频斜率抽取pi对应的Spi(μ),通过以μi为中心频率的带通滤波器,然后对其进行pi阶的FrFT,得到分离后的chirp分量通过门限检测获得Pw,i,最后以构建dechirp滤波器对解线性调频获得
显然,在缺少先验信息的条件下,如果以步进Δp无差别的在所有分数阶搜索,会给硬件系统带来大量额外的计算负担。目前国内外有许多搜索方法来简化这一过程,例如,先以较大的步进获得分辨率较低的参数矩阵,确定谱峰范围后再在该范围内以较小的步进获得高分辨的参数矩阵[15];也有学者通过延时自相关的方法粗略估计chirp分量的调频斜率再以小的步进Δp在相应的分数阶区间上搜索[16]。这些方法都是可行的,考虑到直达波信号有固定的参数模板,本文不再就搜索方法展开讨论,假设模板库是经过前期工作积累获得的,而且短时间内不会有新的参数模板加入。
本节将借助分数阶傅里叶域滤波的方法来实现非合作chirp雷达直达波信号的实时参数估计,根据模板信息,系统可以设置多个通道来并行处理一段直达波信号,如图 3所示,其中第i个通道对应pi阶的FrFT处理,只输出与pi阶FrFT匹配的chirp分量的参数。
图3 直达波信号实时参数估计系统
在实际工程应用中,雷达信号的采集通常伴随着大的采样速率,因此,应用中通常会分帧处理一段直达波信号,待整段信号处理完毕后,系统实时输出该段直达波的参数。图 4展示了某段直达波分量的处理过程。假设该段直达波分量的长度为L(L≫N),即一段直达波信号出现在多个连续观测帧,假设以第t帧为观测起点,直达波在分数阶傅里叶域的冲击函数出现在第t+1帧,消失于t+T+1帧,显然,t+2,…,t+T-1帧的直达波分量满足占空比因此,直接通过FrFT即可获得中间观测帧的参数,而无需进行FrFT的滤波;相应地,只需对占空比ηi<100%的t+1帧与t+T帧的直达波分量进行分数阶傅里叶域滤波操作即可获得整段直达波分量的估计参数。
图4 某通道直达波信号的处理过程
接下来以仿真chirp信号为例,来验证所提方法处理直达波信号的有效性。仿真数据中,采样频率Fs为60 MHz,chirp信号带宽Bi∈{5 MHz,10 MHz,15 MHz},脉宽Pw,i∈{100 μs,1 ms,10 ms},中心频率的选择满足奈奎斯特采样定理。假设这些参数已知,即均为模板库中可选的参数,则可以设计9路通道对接收的直达波信号进行参数估计,9路通道分别对应不同的变换阶p。
估计直达波参数的具体步骤如下:
Step 1: 确定观测帧的长度。分析模板库中所有chirp分量的Pw信息,选择合适尺寸的W0使具有不同Pw的chirp分量占有不同数目的观测帧T。在本例给出的仿真数据中,3种脉宽所占的采样点数分别为6 000、60 000和600 000,因此若取W0=8 192/Fs,则3种脉宽可能占有的连续观测帧数如表 1所示。
表1 脉宽匹配模板
PwT100μs1,21ms8,910ms73,74
W0选择不唯一,本例中取W0=8 192/Fs。
Step 2: 模板匹配。由式(7)知,ki与pi存在一一对应的关系,因此,记录pi阶通道下出现冲击特征的连续帧的长度Ti,可以完成(pi,Ti)与(ki,Pw,i)的匹配。以图5中的仿真信号为例,图5展示了某个时间段内观测帧在直达波信号si{Pw,i,fi,0,Bi=Pw,i×ki}上的滑动情况,信息如下:
·第1帧出现两个窄带chirp分量,记为s1{1 ms,1 MHz,10 MHz},s2{1 ms,2 MHz,5 MHz},都于第8帧消失。
图5 样例信号模板匹配过程
图6 直达波信号模板匹配流程
·第3帧出现一个宽带chirp分量,记为s3{100 μs,1.5 MHz,15 MHz},仅存在于第3帧。
本例中,在预先设计好的9路通道中,有3个通道在这段信号上表现出冲击特性。匹配流程如图 6所示:每个通道在对应的pi阶分数阶傅里叶域上进行局部极大值检测,我们用μ0,j来表示第j帧中分数阶傅里叶域出现的峰值的横坐标;若检测到满足阈值条件的极大值,则记录其直线中心频率并记录出现连续峰值的观测帧的长度T;若峰值消失,且脉宽匹配,T∈模板,则完成匹配,否则不匹配,继续寻找下一个满足条件的chirp分量;此外,对于窄带chirp分量,因为其调频斜率接近于0,分量间p的区别度有限(这是由余切函数的性质决定的),一个窄带分量可能在另一个窄带分量对应的FrFT域也表现出冲击特性,产生虚假的峰值,如图 5中的s1,s2,通道1和通道2中会出现两个峰值。因此,匹配准则中又引入了另一个判决准则,即带宽匹配,若峰值消失时,||μ0,j-μ0,1|-Bi|<ξ,则完全匹配。本例中,通道1中的两个峰值都满足T=8,符合1 ms脉冲匹配条件,但两个chirp分量的|μ0,j-μ0,1|分别等于9.555 MHz和5.12 MHz,显然,前者是完全匹配的,是一个10 MHz/1 ms脉冲,而后者更有可能是一个5 MHz/1 ms脉冲,它在通道1中被筛选掉,相反地,在通道2中被检测出来。而s3为宽带信号,其调频斜率与窄带信号区分度明显,因此可以在通道3中很轻松地完成匹配。在Step 2中,我们可以精确地估计同时借助T粗略地估计
Step 3: 对匹配的分量的起始帧和终止帧做iFrFT,恢复起始帧中和终止帧的chirp分量,精确估计图7以通道1为例对s1进行了分离,展示了带通滤波和iFrFT的过程。
Step 4: 根据各通道匹配检测到的和构造dechirp滤波器对恢复后的单分量chirp信号解线性调频,获得
Step 5: 输出当前分量的参数
Step 6: 重复Step 2~Step 4。
基于时频图像的线段检测方法是一种有效处理式(3)所示的广义chirp信号的方法,本文以文献[7]中SPWVD-线段检测方法为参考,对比本文所述方法展开蒙特卡洛仿真实验,仿真实验在2.6 GHz,Intel Core i7 CPU上行进行。图8展示了在分数阶傅里叶域滤波与SPWVD-线段检测[7]的方法下,上述样例信号中s1的估计精度随信噪比的变化曲线。表 2展示了两种方法处理100 MB仿真数据的计算效率。结果表明,相较于SPWVD-线段检测,本文所提方法所得估计值的估计精度提升了约与的估计精度提升了约0.25%。此外,所提方法在100 MB仿真数据上的运算速度是SPWVD-线段检测方法的10倍。
图7 s1的分离和恢复过程
图8 s1的参数估计误差曲线
表2 处理100 MB数据的计算效率
方法计算时间/s分数阶傅里叶域滤波32.7SPWVD-线段检测291.1
由此我们可以得出结论,基于线段检测的方法由于时频图像的模糊而导致低的估计精度,且不适用于低信噪比背景。基于分数阶傅里叶域滤波的方法可以在分数阶傅里叶域直接完成对k的估计,具有良好的抗噪性。虽然滤波和iFrFT损失了原始信号的部分信息和抗噪性,但对Pw和f0的估计依然可以达到较高的精度。而且,相较于求解全局能量分布的时频图像而言,在有限次的分数阶上进行信号的分离与恢复无疑大大节省了计算成本,因此,基于分数阶傅里叶域滤波的方法拥有更高的计算效率,结合FFT的快速算法,该方法可以在FPGA上实现对已知参数库直达波的实时估计。
本文通过分数阶傅里叶域滤波的方法实现非合作chirp雷达直达波信号的参数估计,本文基于已有的模板库,设计有限的通道对直达波信号进行实时的参数估计,利用分数阶傅里叶变换无损可逆的特性对时频混叠的chirp信号进行分离和恢复,将多分量chirp信号的参数估计问题转化为单分量问题,大大提高了复杂工程背景下chirp信号的参数估计精度和计算效率。实验结果表明,相比基于时频图的线段检测方法,基于分数阶傅里叶域滤波的方法能够达到更高的估计精度,而且易于实现,适用于高采样速率下的大数据背景。
[1] 燕改. 自适应滤波在无源探测中对杂波抑制的应用研究[J]. 信息系统工程,2019(2):152.
[2] 景桐,黄高明,田威. 一种联合目标和外辐射源的无源跟踪模型[J]. 现代雷达,2019,41(5):45-51.
[3] 索毅毅,鲍庆龙,王亚森,等. 基于捷变频非合作雷达辐射源的无源雷达时频同步方法[J]. 雷达科学与技术,2014,12(5):510-516.
SUO Yiyi,BAO Qinglong,WANG Yasen,et al. Time and Frequency Synchronization for Passive Radar Based on Frequency Agile Non-Cooperative Radar Illuminator[J]. Radar Science and Technology,2014,12(5):510-516. (in Chinese)
[4] 冯小平,李晨阳. 线性调频信号参数快速估计[J]. 电子对抗技术,2004,19(6):7-10.
[5] GULUM T O,ERDOGAN A Y,ATA L D,et al. Enhanced LPI Waveform Representation by Ambiguity-Domain Elliptical Gaussian Filtering[J]. IEEE Trans on Aerospace and Electronic Systems,2017,53(2):762-777.
[6] ERDOGAN A Y,GULUM T O,DURAK-ATA L,et al. FMCW Signal Detection and Parameter Extraction by Cross Wigner-Hough Transform[J]. IEEE Trans on Aerospace and Electronic Systems,2017,53(1):334-344.
[7] BOUCHIKHI A,BOUDRAA A,CEXUS J,et al. Analysis of Multicomponent LFM Signals by Teager Huang-Hough Transform[J]. IEEE Trans on Aerospace and Electronic Systems,2014,50(2):1222-1233.
[8] QI L,TAO R,ZHOU S,et al. Detection and Parame-ter Estimation of Multicomponent LFM Signal Based on the Fractional Fourier Transform[J]. Science in China (Series F: Information Sciences),2004,47(2):184.
[9] 段皓楠. 基于分数阶傅里叶变换的线性调频信号估计与分离研究[D]. 南京:南京理工大学,2017.
[10] COWELL D M J,FREEAR S. Separation of Overlapping Linear Frequency Modulated (LFM) Signals Using the Fractional Fourier Transform[J]. IEEE Trans on Ultrasonics,Ferroelectrics,and Frequency Control,2010,57(10):2324-2333.
[11] CANDAN C,KUTAY M A,OZAKTAS H M. The Discrete Fractional Fourier Transform[J]. IEEE Trans on Signal Processing,2000,48(5):1329-1337.
[12] 许道明,张宏伟.雷达低慢小目标检测技术综述[J]. 现代防御技术,2018,46(1):148-155.
[13] 夏建明,郑恒.基于CDR数字音频广播的外辐射源雷达实验研究[J]. 信息技术与信息化,2019 (7):124-127.
[14] 苏卫民,顾红,张先义,等. 基于外辐射源的雷达目标探测与跟踪技术研究[J]. 现代雷达,2005,27(4):19-22.
[15] SERBES A. On the Estimation of LFM Signal Parameters: Analytical Formulation[J]. IEEE Trans on Aerospace and Electronic Systems,2017,54(2):848-860.
[16] 祝俊,陈兵,唐斌. 快速多分量LFM信号的检测与参数估计方法[J]. 电子测量与仪器学报,2008(1):25-29.
苏汉宁 男,1995年生,山东潍坊人,硕士研究生,主要研究方向为阵列信号处理、深度学习和无源雷达。
E-mail:hanningsu18@163.com
鲍庆龙 男,1981年生,吉林河口人,博士,国防科技大学ATR国防科技重点实验室副教授,主要研究方向为雷达数据采集和信号处理。
王 森 男,1993年生,陕西西安人,博士研究生,主要研究方向为阵列信号处理、目标跟踪和无源雷达。
孙玉朋 男,1995年生,山东青岛人,硕士研究生,主要研究方向为雷达信号处理。