基于三次相位变换和虚拟相位函数的空中机动目标检测方法

摘要

一种基于三次相位变换和虚拟相位函数的高机动目标检测和参数估计方法。其包括对机载预警雷达接收到的总回波数据进行杂波抑制，得到杂波抑制后信号；对杂波抑制后信号进行空间相位校正，得到相位校正后的接收信号矩阵；选取多个时刻对各阵元相位校正后的信号进行双线性变换，得到多个时刻的相位函数矩阵；分别对不同时刻各阵元的相位函数补偿空时相位，完成虚拟相位函数构造；构造代价函数，对不同时刻重构的虚拟相位函数进行参数搜索，得到相应时刻的瞬时频率；对得到的一簇瞬时频率进行线性拟合，求得运动参数估值。本发明利用重构虚拟相位函数，等效于增加了相干积累时间，提高了参数估计的精度，能在脉冲点数有限的情况下得到精确估计结果。

说明书

技术领域

本发明属于雷达信号处理技术领域，特别是涉及一种基于三次相位变换和虚拟相位函数的空中机动目标检测方法。

技术背景

机载相控阵雷达又称机载预警雷达，其是以飞机作为载体，具有机动性高、生存力强的特点，与地基雷达相比，其可视范围远，因此可担任警戒、指挥等重要任务，所以在现代战争中的作用越来越重要。但是，由于机载相控阵雷达处于下视状态，其杂波具有分布范围广、强度大等特点，因此在很大程度上淹没了目标信号，从而严重影响了雷达对目标信号的检测和估计性能。传统空时自适应处理(Space-Time Adaptive Processing，STAP)技术可在假设目标回波平稳的条件下，通过空时二维滤波有效地抑制地杂波，从而减小杂波对目标检测及参数估计性能的影响。随着飞行器机动性能的不断提高，目标在运动过程中随时会出现转弯、闪避或其他特殊攻击姿态等高机动现象，当目标运动参数包含加速度及加速度的变化率—急动度(Jerk)时，将这种目标称为高机动目标。高机动目标的多普勒带宽因其高机动性而被严重展宽，使得传统对其的检测及参数估计方法性能下降，因此对高机动目标运动检测和参数估计方法的研究具有重要意义。

当目标做匀急动度运动(即加速度变化率不变)时，目标回波信号为三次相位信号(Cubic Phase Signal，CPS)。近年来，各种CPS信号检测与估计方法不断出现，其中比较典型的方法有离散多项式相位变换(Discrete Polynomial Phase Transform，DPT)、三次相位变换(Cubic Phase Transform，CPT)等。虽然DPT方法能够提供较高的估计精度，但是由于需要对信号进行四阶非线性变换，因此对采样脉冲数和信噪比门限都有较高的要求。CPT方法只对信号进行二阶非线性变换，与DPT方法相比，在低信噪比下能够保持较高的参数估计性能，但是，由于机载预警雷达采样脉冲数较少，因此CPT方法依然难以直接应用。

发明内容

为了解决上述问题，本发明的目的在于提供一种能够提高参数估计精度的基于三次相位变换和虚拟相位函数的空中高机动目标检测和参数估计方法。

为了达到上述目的，本发明提供了基于三次相位变换和虚拟相位函数的空中高机动目标检测和参数估计方法包括按顺序进行的下列步骤：

1)对机载预警雷达接收到的总回波数据进行杂波抑制，得到杂波抑制后的信号；

2)对上述杂波抑制后的信号进行空间相位校正，得到相位校正后的接收信号矩阵；

3)选取多个时刻对上述各阵元相位校正后的信号进行双线性变换，得到多个时刻的相位函数矩阵；

4)分别对不同时刻上述各阵元的相位函数补偿空时相位，完成虚拟相位函数构造；

5)构造代价函数，对不同时刻重构的虚拟相位函数进行参数搜索，得到相应时刻的瞬时频率；

6)对得到的一簇瞬时频率进行线性拟合，求得运动参数估值。

在步骤2)中，所述的对上述杂波抑制后的信号进行空间相位校正，得到相位校正后的接收信号矩阵的方法是：利用阵列流行补偿各阵元杂波抑制后信号的空间相位，得到只含时间相位的接收信号矩阵。

在步骤3)中，所述的选取多个时刻对上述各阵元相位校正后的信号进行双线性变换，得到多个时刻的相位函数矩阵的方法是：利用杂波抑制后数据构造原数据解调数据，利用解调数据实现由杂波抑制后信号向信号相位函数的二次解调，其中杂波抑制后信号为系数是高机动目标运动参数，自变量为采样时间的三次相位函数，信号相位函数为系数是信号瞬时频率，自变量为变换延迟的二次相位函数。

在步骤4)中，所述的分别对不同时刻上述各阵元的相位函数补偿空时相位，完成虚拟相位函数构造的方法是：对各阵元的相位函数进行时间相位补偿，利用各阵元短相干积累的相位函数重构参考阵元长相干积累的相位函数，达到增加相干积累时间的目的。

在步骤5)中，所述的构造代价函数，对不同时刻重构的虚拟相位函数进行参数搜索，得到不同时刻的瞬时频率的方法是：构造代价函数，通过对不同重构参数下代价函数峰值的搜索，确定原信号关于双线性变换时刻的瞬时频率。

在步骤6)中，所述的对得到的一簇瞬时频率进行线性拟合，求得运动参数估值的方法是：通过对多个瞬时频率进行线性拟合，求得高机动目标二阶及三阶参数估值，通过所得的二阶及三阶参数解调原接收信号，实现对高机动目标一阶参数的估计，进而完成对高机动目标运动参数的估计。

本发明提供的基于三次相位变换和虚拟相位函数的空中高机动目标检测和参数估计方法是对多阵元杂波抑制后所得相位函数进行虚拟重构，等效于增加了阵元相位函数相干积累时间，提高了参数估计的精度。本发明能在脉冲点数有限的情况下得到精确的估计结果。

附图说明

图1为本发明提供基于三次相位变换和虚拟相位函数的空中高机动目标检测和参数估计方法流程图。

图2为相控阵体制下天线接收信号示意图。

图3为N个阵元拼接前和拼接后的时频图(N＝8)。

图4为杂波抑制前的总功率谱。

图5为杂波抑制后的总功率谱。

图6为未进行拼接前的时频域能量积累图。

图7为进行拼接前的时频域能量积累图。

图8为本发明方法估计的初始速度均方根误差与CRB界的比较结果图。

图9为本发明方法估计的初始加速度均方根误差与CRB界的比较结果图。

图10为本发明方法估计的急动度均方根误差与CRB界的比较结果图。

具体实施方法

下面结合附图和具体实施例对本发明提供的基于三次相位变换和虚拟相位函数的空中高机动目标检测和参数估计方法进行详细说明。

如图1所示，本发明提供的基于三次相位变换和虚拟相位函数的空中高机动目标检测和参数估计方法包括按顺序进行的下列步骤：

1)对机载预警雷达接收到的总回波数据进行杂波抑制，得到杂波抵制后的信号；

在此阶段中，本发明利用子空间投影技术对机载预警雷达接收到的总回波数据进行杂波抑制，将待检测单元数据投影到杂波子空间的正交子空间中，得到投影后的无杂波数据。杂波加噪声协方差矩阵R为：

R＝E{(x_c+x_n)(x_c+x_n)^H}＝R_c+R_n    (1)

在上式中，x_c表示杂波数据，x_n表示噪声数据，R_c表示杂波协方差矩阵，R_n是噪声项协方差矩阵,E{·}表示统计平均，(·)^H表示转置操作。对杂波加噪声协方差矩阵R进行特征值分解，可以得到：

$>R=\underset{m=1}{\overset{M}{\Sigma }}{\lambda }_m{u}_m{u}_m^H\approx \underset{m=1}{\overset{Q}{\Sigma }}{\lambda }_m{u}_m{u}_m^H+{\sigma }_n^2\underset{m=Q+1}{\overset{M}{\Sigma }}{u}_m{u}_m^H---\left(2\right)$>

其中，λ_m(m＝1,…Q)为Q个大特征值，Q是杂波特征值数目，其余的M-Q个特征值相等，均为u_m(m＝1,…Q)是第m个大特征值对应的特征向量张成的杂波子空间，记为U_C＝span{u₁,…,u_Q}。则其正交补空间的投影矩阵为：

通过上述推导可知，投影后的无杂波数据x'_proj为：

其中，为杂波子空间的正交补空间的投影矩阵，x为机载预警雷达接收到的总回波数据，由于通常为估计得到的噪声方差，λ_i为协方差矩阵的大特征值；而确定杂波子空间维数较复杂，所以我们用杂波加噪声的协方差的逆矩阵R^-1代替抑制杂波，实际中杂波协方差矩阵由参考距离单元数据估计出，因此，杂波抑制后的数据x_proj可写为：

x_proj＝R^-1x＝x_i+x_nc    (5)

其中x_i为理想无噪声目标回波，x_nc为杂波抑制后的噪声项。

2)对上述杂波抑制后的信号进行空间相位校正，得到相位校正后的接收信号矩阵：

如图2所示，机载平台沿航向方向放置均匀线阵，阵元数为N，相干处理脉冲数为K，阵元间距d＝0.5λ，其中λ为雷达发射脉冲波长。当目标做匀急动度运动时，设目标初始速度为v，初始加速度为a，急动度为j_a，则对于机载预警雷达的每个阵元，其接收到的信号均为三次相位信号：

$>\left(\begin{array}{c}{x}_n\left(k\right)={\mathrm{Ae}}^{\phi \left(k\right)}={\mathrm{Ae}}^{j2\pi \frac{2\mathrm{vk}}{\lambda }}·e^{j2\pi \frac{{\mathrm{ak}}^2}{\lambda }}·e^{j2\pi \frac{{\mathrm{ak}}^2}{\lambda }}·e^{j2\pi \frac{1}{3}\frac{j_a{k}^3}{\lambda }},k=\mathrm{0,1}···,K,n=\mathrm{1,2}···,N\\ ={\mathrm{Ae}}^{j({\beta }_{1}k+{\beta }_2{k}^2{\beta }_3{k}^3)}\end{array}\right)---\left(6\right)$>

其中A为目标回波的复包络，β₁＝4πv/λ表示目标初始速度引起的多普勒一次项，β₂＝2πa/λ表示目标加速度引起的多普勒二次项，β₃＝2πj_a/3λ表示目标急动度引起的多普勒三次项。则目标的时域导向矢量可以表示为：

$>a\left({\omega }_t\right)={\left(\begin{array}{cccc}1& e^{j({\beta }_1·1+{\beta }_2·1^2+{\beta }_3·1^3)}& ···& e^{j({\beta }_1(K-1)+{\beta }_2{(K-1)}^2+{\beta }_3{(K-1)}^3)}\end{array}\right)}_{K\times 1}^T---\left(7\right)$>

空域导向矢量为：

$>a\left({\omega }_s\right)={\left(\begin{array}{cccc}1& e^{j2\pi \frac{d\cos \theta }{\lambda }}& ···& e^{j2\pi (N-1)\frac{d\cos \theta }{\lambda }}\end{array}\right)}_{N\times 1}^T---\left(8\right)$>

在式(7)及式(8)中，令ω_tk＝β₁(k-1)+β₂(k-1)²+β₃(k-1)³，(k＝0,1,…,K-1)表示第k个脉冲的时间相位，(n＝1,2,…,N)表示第n个阵元的空间相位，则将理想无噪声目标回波向量x_i可以表示为：

$>\left(\begin{array}{c}{x}_i=\tilde{A}a\left({\omega }_t\right)\otimes a\left({\omega }_s\right)=\tilde{A}{[1,e^{{\mathrm{j\omega }}_{t2}},···,e^{{\mathrm{j\omega }}_{\mathrm{tK}}}]}_{K\times 1}^T\otimes {[{1,e}^{{\mathrm{j\omega }}_{s2}},···,e^{{\mathrm{j\omega }}_{\mathrm{sN}}}]}_{N\times 1}^T\\ =\tilde{A}{[e^{{\mathrm{j\omega }}_{s1}+{\mathrm{j\omega }}_{t1}},e^{{\mathrm{j\omega }}_{s1}+{\mathrm{j\omega }}_{t2}},···,e^{{\mathrm{j\omega }}_{\mathrm{sN}}{+\mathrm{j\omega }}_{t1}},,e^{{\mathrm{j\omega }}_{\mathrm{sN}}+{\mathrm{j\omega }}_{t2}},···,e^{{\mathrm{j\omega }}_{\mathrm{sN}}{\mathrm{j\omega }}_{\mathrm{tK}}}]}_{\mathrm{NK}\times 1}^T\end{array}\right)---\left(9\right)$>

将x_i展开成为K×N阶接收信号矩阵X_i的形式：

其中，矩阵每一列代表一个阵元接收数据，可以发现，各阵元相应数据间仅相差一空间相位。以1号阵元为参考阵元，定义各阵元到参考阵元空间相位矫正矩阵：

则对原接收信号进行空间相位矫正后的信号可以表示为：

3)选取多个时刻对上述各阵元相位校正后的信号进行双线性变换，得到多个时刻的相位函数矩阵；

式(12)中，X的每个列向量均具有三次相位信号形式，并且依式(6)可以表示为：

$>s_n\left(k\right){\tilde{A}e}^{{\omega }_{\mathrm{sn}}}{e}^{{\mathrm{j\beta }}_{1}k+{\mathrm{j\beta }}_2{k}^2+{\mathrm{j\beta }}_3{k}^3}={\tilde{A}e}^{{\omega }_{\mathrm{sn}}}{e}^{{\omega }_{\mathrm{tk}}},(k=\mathrm{0,1}···,K-1,n=\mathrm{1,2},···,N)---\left(13\right)$>

定义k₀时刻回波信号的双线性变换：

$>\left(\begin{array}{c}{q}_{\mathrm{nk}0}\left(\tau \right)=s_n(k_0+\tau )s_n(k_0-\tau )\\ ={\tilde{A}}^2{e}^{j({\beta }_1{k}_0+{\beta }_2{k}_0^2+{\beta }_2{k}_0^3)e^{j2({\beta }_2+3{\beta }_3{k}_0){\tau }^2}}\end{array}\right),\left(\begin{array}{c}\tau =\mathrm{0,1},···,t_0-1,s.t.k_0<\frac{K-1}{2}\\ \tau =\mathrm{0,1},···,K{-t}_0,s.t.k_0\ge \frac{K-1}{2}\end{array}\right),K=\mathrm{odd}---\left(14\right)$>

其中，称q_nk0(τ)为k₀时刻第n个阵元的接收信号s_n(k)的相位函数，τ≥0是引入的变换延迟，为变换后信号的复包络。

Ω(k₀)＝2(a₂+3a₃k₀)    (15)

在k₀时刻的瞬时频率(instantaneous frequency rate，IFR)。

由式(14)可知，经过k₀的双线性变换，将关于时间k的三次相位信号转化为了关于延迟τ的二次相位信号。在k_i(1<i<K)时刻对每个阵元接收的信号进行双线性变换，得到k_i时刻总回波信号的相位函数矩阵：

$>D_i=\tilde{A}{\left(\begin{array}{cccc}{q}_{i1}& q_{i2}& ···& q_{\mathrm{iN}}\end{array}\right)}_{\frac{K-1}{2}\times N}---\left(16\right)$>

其中，表示第n个阵元的相位函数，Ω(k_i)为k_i时刻回波信号的瞬时频率。

4)分别对不同时刻各阵元的相位函数补偿空时相位，完成虚拟相位函数构造；

由式(16)可知，各阵元的相位函数具有二次相位信号形式，每个阵元相位函数只相差一定的空时相位，如图3中实线所示。以某一阵元接收信号的相位函数为参考，补偿其它各阵元相位函数的空时相位，然后进行拼接，等效为增加了参考阵元相位函数相干积累时间的效果，如图3中虚线所示。不失一般性，当阵元数为2时，根据两阵元数据重构的相位函数可以表示为：

其中q_i1、q_i2分别表示k_i时刻第一个阵元及第二个阵元接收信号的相位函数，为Hadamard积，为q_i1与q_i2之间的时间相位差。

可以看出，对q_i2补偿时间相位后再跟q_i1拼接，等效为增加了q_i1时间样本点数。当阵元数增加到N时，k_i时刻重构相位函数的表达式为：

其中：

5)构造代价函数，对不同时刻重构的虚拟相位函数进行参数搜索，得到相应时刻的瞬时频率；

由上述分析可知，上述方法中，补偿的空时相位中包含了频率Ω，这个参数由双线性变换时刻k以及未知的目标运动参数决定，即：

Ω(k)＝2(β₂+3β₃k)    (21)

瞬时频率未知导致无法直接对相位函数矩阵多行数据进行拼接，本发明对重构的相位函数进行信号瞬时频率Ω的搜索，以此构造代价函数，实现对某个双线性变换时刻k_i信号瞬时频率Ω(k_i)的估计，则对瞬时频率Ω进行参数估计的代价函数可以表示为：

$>\hat{\Omega }\left(k_i\right)\arg \underset{\Omega }{\min }{\left|\right|q_{\mathrm{ircp}}-q_{\mathrm{iz}}\left|\right|}_2---\left(22\right)$>

其中

表示对接收信号k_i时刻的相位函数进行重构相位函数，q_iz表示当补偿的瞬时频率为Ω_i时，理想无噪声的相位函数。

6)对得到的一簇瞬时频率进行线性拟合，求得运动参数估值。

式(13)表明，信号在变换时刻k处得到的瞬时频率Ω为待估参数β₂,β₃关于变换时刻k的一次函数，为了实现对运动参数的估计，对一簇时刻瞬时频率的估计为：

$>\hat{\Omega }\left(k_i\right)\arg \underset{\Omega }{\min }{\left|\right|q_{\mathrm{ircp}}{-q}_{\mathrm{iz}}\left|\right|}_2,0<k_i=k_1,k_2,···,k_M<K-1---\left(24\right)$>

其中，M表示进行估计的瞬时频率的个数，然后对所得瞬时频率进行线性拟合即可求得β₂,β₃的估计值：

$>\left(\begin{array}{c}{\hat{\beta }}_3=\frac{1}{6}\frac{M\underset{i}{\overset{M}{\Sigma }}{k}_i\Omega \left(k_i\right)-\underset{i}{\overset{M}{\Sigma }}{k}_i\underset{i}{\overset{M}{\Sigma }}\Omega \left(k_i\right)}{M\underset{i}{\overset{M}{\Sigma }}{k}_i^2-\underset{i}{\overset{M}{\Sigma }}{k}_i\underset{i}{\overset{M}{\Sigma }}{k}_i}\\ {\hat{\beta }}_2=\frac{1}{2}\frac{\underset{i}{\overset{M}{\Sigma }}\Omega \left(k_i\right)-{\hat{\beta }}_3\underset{i}{\overset{M}{\Sigma }}\Omega \left(k_i\right)}{M}\end{array}\right)---\left(25\right)$>

利用解调式(6)，可得

$>{\hat{\beta }}_1=\arg \underset{a_1}{\max }\left|\underset{t=-(N-1)/2}{\overset{(N-1)/2}{\Sigma }}s\left(t\right)e^{{\mathrm{j\beta }}_{1}t{\beta }_2{t}^2-{\beta }_3{t}^3}\right|---\left(26\right)$>

再由式(6)求解运动参数估计值，即可完成对高级动目标的参数估计：

$>\left(\begin{array}{c}\hat{v}=\frac{\lambda {\hat{\beta }}_1}{4\pi }\\ \hat{a}=\frac{\lambda {\hat{\beta }}_2}{2\pi }\\ {\hat{j}}_a=\frac{3\lambda {\hat{\beta }}_3}{2\pi }\end{array}\right)---\left(27\right)$>

本发明提供的基于三次相位变换和虚拟相位函数的空中高机动目标检测方法的效果可以通过以下仿真结果进一步说明。

仿真参数设置：天线为阵元数N＝8的正侧阵理想均匀线阵，阵元间距d＝0.5λ。载机速度为120m/s，雷达工作波长为0.32m，平台高度为10km，雷达的距离分辨率为20m，脉冲重复频率为1500Hz，相干处理脉冲数K＝128，输入信噪比为SNR＝0dB，杂噪比为CNR＝50dB。机动目标处于检测单元内，处于方位角90°处，初始速度为20.1m/s，加速度为99.1m/s²，急动度为40.1m/s³，实验中假设目标方位已知，该实验中估计参数均方误差均进行了400次蒙特卡罗实验。

图4为杂波抑制前功率谱，从图中可以看出，由于信杂比很低，因此信号完全被淹没在杂波中。图5为杂波抑制后的功率谱，从图中可以看出，杂波得到了有效的抑制，目标凸显出来，但是由于存在二阶和三阶参量，因此其在多普勒域仍存在一定的展宽。

图6、图7为不同方法下进行能量积累的功率谱，可以看到，与图6重构相位函数相比，图7中通过重构相位函数，对接收数据进行相干积累，可使目标能量大大提高，能够很好地进行目标检测和估计。

图8-图10为不同方法得到的参数均方根误差与CRB界的比较结果图，其中图8为初始速度均方根误差与CRB界的比较结果图，图9为初始加速度均方根误差与CRB界的比较结果图，图10为急动度均方根误差与CRB界的比较结果图，可以看出，较传统CPT方法，本发明方法估计性能更接近CRB界，因此估计效果较好。