线性调频连续波雷达低空小目标长时间相参积累方法

2022-08-06郭立民莫禹涵刘鲁涛郭晓冉

国防科技大学学报 2022年4期

郭立民，莫禹涵，刘鲁涛，郭晓冉

(1. 哈尔滨工程大学信息与通信工程学院，黑龙江哈尔滨 150001; 2. 中国人民解放军32181部队，河北石家庄 050000)

随着无人机技术的发展，操作小型无人机的技术门槛越来越低。为人们生活娱乐带来便利的同时也为重点区域的低空安全带来了极大威胁[1]。目前探测无人机的手段包括光学探测和雷达探测等。其中光学探测的精度较高，实现简单，但会受到光线条件的影响，且作用距离也受制于光学设备的分辨率，不利于全天候、远距离的探测；而雷达探测可以不受时间和气候的影响，作用距离较大，但实现起来较复杂。因此很多无人机反制系统选择复合型方案，选用光学系统和雷达系统结合的方式实现无人机反制系统。作为总体系统的一部分，雷达系统的成本必然不能过高。因此研制低成本、小型化雷达来探测无人机势在必行。线性调频连续波(linear frequency modulation continuous wave，LFMCW)雷达结构简单、成本低、易于实现，且无距离盲区，十分适合探测无人机目标。但受其成本限制，硬件参数一般较差。因此深入研究适用于LFMCW雷达探测小型无人机的方法十分具有现实意义。与常规的目标不同，无人机的雷达散射截面积(radar cross section，RCS)很小。在远距离探测目标时，无人机的回波很容易淹没在杂波和噪声当中，需要增加发射信号的周期和相参积累的时间以获取更大信噪比[2-3]。这就意味着雷达对无人机的一次积累时间会很长，目标会在积累周期内跨距离单元。直接对回波进行动目标检测(moving target detect，MTD)处理会导致目标能量无法有效聚集，即便增加积累周期也无法获得理想的信噪比增益[4]。

关于长时间相参积累的研究成果很多，核心内容都是校正目标在积累周期内的走动，实现对远距离微弱目标的长时间相参积累。其中包括霍夫(Hough)变换[5]、拉东(Radon)变换[6-9]等，这些方法引入了图像处理中的一些思想，取得了一定的效果。文献[6]中的拉东-傅里叶变换(Radon-Fourier transform，RFT)综合了图像和雷达信号处理中的思想，很大程度上将长时间接收的目标能量有效积累起来。但是目标轨迹是由时域移位或寻址运算确定的，偏移量为整数，这样会损失一部分目标能量。若要减少能量损失，可以对数据进行插值处理，又增加了算法的运算量。文献[7]针对这种缺点，提出了基于频域补偿实现RFT的方法及其快速实现方法。利用Chirp-Z变换(Chirp-Z transform，CZT)减小了文中直接计算式的运算量，具有一定的参考意义。但在无人机反制系统中，LFMCW雷达仅作为其中一部分，信号处理能力比不上文献[7]中的脉冲多普勒(pulse Doppler，PD)雷达。直接应用文中的算法运算量较大，不利于信号的实时处理。因此本文立足于低成本、小型化LFMCW雷达的工作特点，提出了基于时域差频信号CZT的RFT实现方法，规避了上述算法中的一些不足，与传统的频域RFT算法相比减少了快速傅里叶变换(fast Fourier transform，FFT)和逆快速傅里叶变换(invert fast Fourier transform，IFFT)的次数。在仿真数据和实测数据上进行长时间的相参积累验证，取得了理想的积累效果。

1 LFMCW雷达基本原理

LFMCW雷达的收发天线分置，发射线性调频连续波信号，能够以较小的发射功率获得较大的探测距离，且无距离盲区。假设雷达收发信号均为复指数信号，则发射信号可由式(1)表示。

(1)

式中:f0为载波频率，μ=B/T为调频斜率。B为发射的线性调频信号的带宽，T为调频周期。假设空中的一个匀速运动目标位于距离雷达R0处，速度为v(远离雷达为正)，则目标回波经过接收机放大，混频滤波后归一化处理的差频信号如式(2)所示。

(2)

此差频信号即为接收机输出到信号处理单元模数转换器(analog to digital converter, ADC)的信号，其中fd为多普勒频率。从式(2)中可以看出，差频信号主要由三部分构成。首先是与目标初始距离R0有关的相位项，代表着回波信号的时延。第二项是由目标运动引起的多普勒相位项，随着慢时间tm的变化，运动目标回波幅度会受到多普勒频率的调制。第三项是差频信号的主体，可视为中心频率为fd+2μ(R0+vtm)/c的单载频复指数信号。若仅关注某一个周期的回波信号，即将慢时间tm视为常量，则前两项为常数。差频信号就是单载频复指数信号，频率与目标距离相关。

若要得到直观的目标距离信息，需要对差频信号做一次快时间FFT运算，既能够解析距离信息，又实现了时域信号的脉冲压缩处理，提高了信噪比。对差频信号做快时间FFT运算后的表达式如式(3)所示。

(3)

由于FFT运算相当于为信号加了矩形窗，所以差频信号的频谱是辛克函数的形状，谱峰位于f=fd+2μ(R0+vtm)/c。对不同周期的回波进行MTD处理，可以实现相参积累，过程如式(4)所示。

(4)

式中的f与距离r等价，多普勒频率fd与速度v等价。故式中的Sb(f,fd)等价于Sb(r,v)。MTD处理可由慢时间FFT实现。至此经过两个维度的FFT运算，就可以得到目标的距离-多普勒(range-Doppler，RD)二维谱。并由此得到目标距离和速度信息。

2 RFT快速实现方法

2.1 频域实现RFT方法原理

由第1节中的分析可知，回波差频信号在经过快时间FFT处理后，为峰值位于f=fd+2μ(R0+vtm)/c处的辛克函数。设r=cf/(2μ)，式(3)可改写为式(5)。

(5)

式中，cfd/(2μ)项为距离-速度耦合项，可通过发射对称三角调频信号，对回波上下扫频段目标进行配对来抑制[10]。但式(5)中的R0+vtm项则会导致每个周期回波的峰值位置随着慢时间tm发生变化。若vtm大于一个距离单元，就会发生距离走动现象。因为目标匀速运动，所以相邻周期的距离走动ΔR=vT为定值。只要将Sb(r,tm)中的距离信号进行相应的移位，保证目标峰值在同一距离单元即可实现距离走动的校正。但对式(5)直接进行移位比较困难，可以对式(5)做傅里叶变换，再乘以相移因子exp(j2πf′vtm)，最后经过傅里叶反变换后即实现了距离校正。其中f′为距离变量r所对应的频率域变量。设Sb(r,tm)的傅里叶变换为Sb(f′,tm)，频域实现RFT的过程可由式(6)表示。

(6)

式(6)可以用FFT来实现。由于不知道目标实际的运动速度v，在用此方法进行RFT运算的时候需要对可能的速度v进行遍历，选择积累峰值最高的作为最终结果。式(6)的两次积分可以调换顺序得到如式(7)的结果。

(7)

综上，用FFT在频域实现RFT的框图如图1所示。

图1 频域RFT实现框图Fig.1 Frequency domain RFT implementation block diagram

2.2 时域RFT快速实现算法

频域实现RFT的基本思想是利用傅里叶变换的时移定理，对距离维信号做一次FFT，通过频域移相实现了原始序列的位移。在PD雷达中，脉冲压缩是通过匹配滤波实现的。而匹配滤波中的卷积运算多通过频域相乘实现。因此PD雷达实现频域RFT时可以在匹配滤波的过程中进行，不会额外增加FFT的次数。而LFMCW雷达的脉冲压缩则是通过对差频信号做一次快时间FFT实现的，用2.1节的方法需要额外做一次FFT和一次IFFT。增加了复数乘法次数，消耗了大量的时间和硬件资源。

观察1.1节中的式(3)，快时间维信号已经是频域信号。2.1节中分析的对Sb(r,tm)中不同周期的距离信号做移位对齐，实质上就是对Sb(f,tm)中的频率信号做频移。由傅里叶变换的时频对称定理可知，频域信号的频移可由时域移相实现。时域信号可通过对Sb(f,tm)进行傅里叶反变换得到，而差频信号sb(t,tm)就是Sb(f,tm)的傅里叶反变换。直接对差频信号进行时域移相即可实现目标距离走动的校正。

式(3)中，造成不同周期回波走动的项为2μvtm/c。所以时域差频信号相应的相移因子为exp(-j4πμvtmt/c)。每个周期的差频信号都乘以相移因子exp(-j4πμvtmt/c)后，再进行第1节所述的二维FFT，即可得到时域RFT结果。和2.1节一样，两次FFT的顺序可以调换。用二重积分表示此过程，如式(8)所示。

Sb(f,fd)=∬sb(t,tm)·exp(-j4πμvtmt/c)·

exp(-j2πfdtm)·exp(-j2πft)dtmdt

(8)

式(8)用FFT实现时，和式(7)一样，对需要对可能的速度v进行遍历，最后选取最大值作为积累结果。结构如图2所示。不同的是式(8)直接对差频信号sb(t,tm)进行运算，而式(7)是对差频信号快时间FFT后等价的Sb(r,tm)进行运算。因此式(8)比式(7)要少了一次FFT运算和一次IFFT运算。与频域实现RFT相比，时域实现更适用于LFMCW雷达。

图2 时域RFT实现框图Fig.2 Time domain RFT implementation block diagram

将式(8)进一步合并整理，可以得到式(9)。

Sb(f,fd)=∬sb(t,tm)·exp[-j2π(1+tμ/f0)fdtm]dtm·

exp(-j2πft)dt

(9)

将式(9)中的慢时间积分离散化处理，设a=1+tμ/f0，快时间采样点数为N，用n表示快时间离散变量，对应快时间t；积累周期数为M，用m表示慢时间离散变量，对应慢时间tm；k表示多普勒采样点，对应fd。从而得到式(10)。

(10)

直接计算式(10)运算量较大，可以将其做一定的变形得到式(11)。

exp(-jπam2/M)·exp[jπa(k-m)2/M]

(11)

通过对比可以发现，式(11)与式(12)CZT的表达式一致，即可以通过CZT来计算式(10)，达到沿着目标运动轨迹积累的目的。

(12)

CZT变换能够实现任意长度的傅里叶变换，多用来进行频谱细化。CZT可由FFT实现，相比直接计算式(10)计算量会小很多。利用CZT实现时域快速RFT的框图如图3所示。

图3 基于CZT的时域快速RFT实现框图Fig.3 Block diagram of fast RFT implementation in time domain based on CZT

图3中l表示快时间频域的离散变量，与f对应。图3中系数exp(jπam2/M)中的变量都是由雷达参数决定的，可事先计算其FFT。与图1频域RFT的方法相比，不需要对速度v遍历，极大提高了系统的效率。频域RFT也可以通过CZT来实现[7]，此处不对其做推导，仅在下一章对频域CZT-RFT做运算量的分析。

3 仿真实验对比分析

3.1 多普勒模糊情况分析

3.2 仿真结果分析

本节对LFMCW雷达回波进行仿真，比较了频域RFT、时域CZT-RFT和MTD三种积累方法的积累性能。雷达和目标的仿真参数如下：载频f0=16 GHz，调频带宽B=100 MHz，调频周期T=2 ms(分上下调频)，采样频率fs=2.5 MHz，快时间采样点数N=2 048，积累周期数M=512。目标距离R=600 m，速度v=15 m/s。积累后距离维仅截取前512个点进行显示。时域差频信号信噪比-20 dB时的仿真结果如下所示。仿真中LFMCW雷达发射信号为对称三角调频，即一个周期分为上调频和下调频两段，一个周期共2 ms，对应的重频fr=500 Hz。下面给出的均为上调频段的结果。

图4是对动目标显示(moving target indication, MTI)后的数据做快时间FFT的结果，可以看出目标从600 m处的位置远离雷达飞行，目标走动比较明显。直接对图4做相参积累得到的结果如图5所示，能够看出目标能量没有被完全聚集在一点上，能量分散严重。采用文献[7]中频域RFT对回波进行距离走动补偿和相参积累，得到的结果如图6所示。从图中能看到频域RFT能够很大限度地将能量聚集在RD谱上一点，校正了目标在积累周期内的距离走动。时域CZT-RFT的积累结果如图7所示，对比图6结果，时域RFT可以达到频域RFT的积累效果，幅度较频域RFT更高。这是因为频域处理时忽略了目标在快时间内运动带来的多普勒变化，造成目标多普勒维能量没有完全聚集。而时域RFT处理时则最大限度地保留了这些信息。为比较几种积累方法的性能，进行了3 100次蒙特卡洛实验，在虚警概率为10-6时得出了几种积累方法在不同信噪比下的性能曲线，如图8所示。

图4 目标回波包络走动Fig.4 Target echo envelope migration

图5 直接MTD结果Fig.5 MTD accumulation results

图6 频域RFT积累结果Fig.6 Frequency domain RFT accumulation results

图7 时域CZT-RFT积累结果Fig.7 Time domain CZT-RFT accumulation results

图8 虚警概率为10-6时的检测概率曲线Fig.8 Detection probability curve when false alarm probability is 10-6

图8中的信噪比为时域差频信号的信噪比，由于时域CZT-RFT是时域RFT的快速算法，本质上是相同的，因此图中仅给出时域CZT-RFT的检测概率曲线与另外两种积累方法对比。可以看出，频域RFT和时域CZT-RFT都比传统的MTD检测性能好得多。本文所提的时域CZT-RFT比文献[7]中的频域RFT的检测性能略好，检测概率相同的情况下对信噪比的要求小了1～2 dB。

3.3 算法计算量分析

由于MTD实现比较简单，因此计算量要小于频域RFT和时域CZT-RFT方法。但在长时间相参积累时MTD的性能远低于前两种方法，因此本小节仅讨论频域RFT和时域CZT-RFT的计算量。频域RFT在用图1的方法实现时，需要先通过对差频信号做快时间FFT得到Sb(r,tm)。设快时间采样点为N，积累周期为M。对sb(t,tm)做M次N点FFT需要的复数乘法次数为0.5MNlog2N。加上图1中的运算后，总的复数乘法次数为MN(1+1.5log2N+0.5log2M)。由2.1节分析可知，利用此方法实现RFT时要对速度v进行遍历。搜索步长Δv一般为积累周期内能够引起目标走动一个距离单元的速度，即Δv=ΔR/(MT)。对于本章的仿真参数，ΔR=1.831 m，M=512，T=2 ms。所以Δv=1.788 m/s。根据无人机手册，小型四旋翼无人机的速度一般不会超过20 m/s，因此速度搜索范围为[-20 m/s， 20 m/s]。共需要搜索23次。因此频域RFT总的复数乘法次数为MN(23+12.5log2N+11.5log2M)。

文献[7]中也推导了频域RFT通过CZT实现的方法，性能和图1的方法一致，计算量更小。不需要对v遍历，但需要对模糊数i遍历。根据本章的仿真参数可以计算出无模糊测速范围为[-2.34 m/s， 2.34 m/s]。需要遍历的模糊数为4。考虑速度方向，总遍历次数为9。频域CZT-RFT总的复数乘法次数为MN(27+5.5log2N+9log2M)。

由图3方式实现的时域CZT-RFT需要对模糊数进行9次遍历，总复数乘法次数为9MN(3+log2M+0.5log2N)。三种RFT实现方法的复数乘法次数随着积累周期M的变化趋势如图9所示。

图9 几种积累方式的运算量对比Fig.9 Comparison of calculation amount of several accumulation methods

从图9中的对比可以看出，时域CZT-RFT比频域RFT的计算量小得多。即使是文献[7]中快速实现频域RFT的频域CZT-RFT算法的计算量也要高于本文的时域CZT-RFT算法。以512周期积累为例，时域CZT-RFT比频域CZT-RFT复数乘法次数要少1.153 4×107次。

由图8和图9可知，本文的时域CZT-RFT算法在长时间相参积累时的检测性能要远强于传统的MTD，略强于文献[7]中的频域RFT。而在计算量上时域CZT-RFT也小于频域RFT及其快速算法频域CZT-RFT，更易于工程实现。

4 实测数据分析

本章在外场采集了大疆精灵4无人机飞行的雷达回波，采用MTD和时域CZT-RFT算法对回波进行相参积累并对比，观察本文算法是否能够提高目标信噪比。LFMCW雷达具体参数见表1。

表1 LFMCW雷达参数Tab.1 LFMCW radar parameters

所录取的实测数据为ADC采集到的差频信号。对其中的512周期回波截取中间2 048个采样点进行快时间FFT后，可观察到距离走动情况如图10所示。

图10 实测的目标走动现象Fig.10 Target moving phenomenon actually tested

从图10中能够看出目标在积累周期内从298.9 m处运动到了317.2 m处，走动了约3个距离单元。若直接对此结果进行MTD处理后得到图11的结果。

图11 实测数据直接MTD积累结果Fig.11 Measured data MTD accumulation result

由于目标较近，直接MTD处理的信噪比能够达到约39.32 dB，谱峰位于311.1 m处。下面给出时域CZT-RFT的积累结果，如图12所示。

图12 实测数据时域CZT-RFT积累结果Fig.12 Measured data time domain CZT-RFT accumulation result

从图12中可以观察到经过CZT-RFT对距离走动的校正后，谱峰位于298.9 m处。目标幅度增大了约2.77 dB，使得信噪比达到了42.09 dB。在多普勒维度上，可以看到目标在多普勒维度上还存在展宽现象。一方面目标相对雷达波束不是严格的匀速运动[11]，另一方面是由无人机旋翼高速旋转所导致的微多普勒效应[12-13]。因此也导致目标幅度没有增大到理论的4.77 dB。这些现象也只有在目标距离走动被完全校正后才能看到。数据录取时间共75 s，扇扫范围为[-5°，5°]，间隔2°，波束驻留时间为1.026 s。共录取了12次扇扫的数据，除两次目标运动超过雷达俯仰波束范围外，共检测到10次目标，恢复的航迹与无人机飞行记录的对比情况如图13所示，证明了本文算法的有效性。

图13(a)的无人机位置即为开始录取数据时的位置，白色轨迹即为开始录取数据后的无人机轨迹；图13(b)为还原的航迹，与飞行记录一致。从实测结果也可以看出，由于积累时间较长，实测数据的积累效果没有3.2节中的仿真结果理想。时域CZT-RFT算法能够补偿目标在长时间积累中的距离走动，但目标变速运动时会造成目标能量在多普勒上的扩散，同时长时间积累时可能存在目标跨波束的情况，影响积累效果。因此，在运用本文算法时需要根据所要探测的距离以及雷达自身的参数如波束宽度、脉冲重复周期等来预估相参积累的时间TCIT。保证在相参积累过程中目标的切向运动距离ΔRt不会跨越天线波束照射范围，即ΔRt=vtTCIT≤Rrθa，其中vt为目标切向速度分量，Rr为感兴趣的探测距离，θa为方位波束宽度。目标速度变化范围Δv也要尽可能小于一个速度分辨单元，即Δv≤λfr/(2M)=λ/(2TCIT)，以获取最佳的积累性能。