CN103412294B - 基于双重直积分解的机载雷达空时三维杂波抑制方法 - Google Patents

Abstract

本发明公开了一种基于双重直积分解的机载雷达空时三维杂波抑制方法，主要解决机载雷达空时三维杂波抑制计算复杂度高，样本需求量大，阵元误差导致检测性能明显下降的问题。其实现过程是：1）将空时三维滤波器权矢量进行一次直积分解；2）将空域权矢量组进行二次直积分解；3）通过截断处理，得到空时三维滤波器权矢量的双重直积分解形式；4）用双重双迭代算法，求解时域权矢量组、空域方位维权矢量组和空域俯仰维权矢量组；5）根据所求的时域权矢量组、空域方位维权矢量组和空域俯仰维权矢量组，恢复出空时三维滤波器权矢量，进行杂波抑制。本发明计算复杂度低，样本需求量小，具有良好的误差稳健性，能提高动目标检测性能。

Description

基于双重直积分解的机载雷达空时三维杂波抑制方法

技术领域

本发明属于雷达信号处理技术领域，具体涉及一种基于双重直积分解的机载雷达空时三维降维自适应杂波抑制方法，可用于提高对运动目标的检测性能。

背景技术

机载雷达被用来识别强杂波环境中的运动目标，杂波散射单元由于方向不同，对应的杂波多普率频率也不同，这一现象称之为空时耦合特性。寻找简单高效的杂波抑制方法，一直是机载雷达领域的研究重点之一。三维空时自适应处理（3D-STAP）是在空域二维和时域上自适应处理，能够充分利用俯仰维自由度，对阵元误差具有更好的稳健性。同时，在远射程和一定的脉冲重复频率（PRF）下，俯仰维自适应波束可以抑制距离模糊杂波，因此3D-STAP可以减轻距离模糊问题。尽管3D-STAP具有一些更好的性能，但随着加入俯仰维进行处理，全维自适应处理器维数更高，计算复杂度更大，难以实时处理，且随着处理器维数的增加，三维处理需要更多独立同分布的样本。为了降低计算复杂度和样本需求量，同时保证自适应损益最小，3D-STAP降维技术得到深入研究。

目前，在机载雷达3D-STAP中，经常用的方法有扩展因子化方法（EFA），该方法先在时域进行多普勒滤波，然后在空域俯仰维和方位维进行联合自适应处理，当存在阵元误差时，杂波谱在空域发生扩散，有更多的杂波进入时域多普勒通道，算法性能下降，且该方法还需对空域全维处理器求逆，计算复杂度高，不利于实时处理。此外，已经证明，为保证自适应损益小于3dB，EFA方法所需独立同分布样本数应大于空域全维处理器维数的2倍以上，而在实际快时变的杂波环境中，通常难以得到如此多的独立同分布样本。同样，现有的其他一些机载雷达3D-STAP方法，也会出现类似的问题。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于双重直积分解的机载雷达杂波抑制方法，降低计算复杂度和样本需求量，提高误差稳健性，改善动目标检测性能。

实现本发明目的的技术方案可以概括为：首先，将空时三维滤波器权矢量进行一次直积分解，得到空时三维滤波器权矢量关于时域权矢量组、空域权矢量组的直积构成形式；其次，将空域权矢量组进行二次直积分解，得到空域权矢量组关于空域方位维权矢量组、空域俯仰维权矢量组的直积构成形式；然后，截取多对时域权矢量组、空域方位维权矢量组及空域俯仰维权矢量组逼近空时三维滤波器权矢量；接着，利用双重双迭代算法（BBIA），求解截取的多对时域权矢量组、空域方位维权矢量组及空域俯仰维权矢量组；最后，由所求得的多对权矢量组恢复出空时三维滤波器权矢量，用来抑制杂波。

假设机载雷达阵列为M×N矩形面阵，其中M表示面阵行数，N表示面阵列数，一个相干处理时间（CPI）内的脉冲数为K。

具体实现过程如下：

（1）将空时三维滤波器权矢量w进行一次直积分解，将空时三维滤波器权矢量w=[w_1,1,1…w_M,1,1…w_1,2,1…w_M,2,1…w_1,N,K…w_M,N,K]^T排列成一个K×MN维的权矩阵W₁，其中符号T表示转置，利用满秩分解方法将权矩阵W₁分解为其中r₁表示权矩阵W₁的秩，u_n表示空域权矢量组，d_n表示时域权矢量组，符号H表示共轭转置；将等式左右两边列矢量化，得到一次直积分解其中符号*表示取共轭，符号表示直积；

（2）将空域权矢量组u_n进行二次直积分解，将空域权矢量u_n排成M×N维矩阵U_n，利用满秩分解方法将权矩阵U_n分解为其中r₂表示U_n的秩，θ_n表示空域俯仰维权矢量组，表示空域方位维权矢量组；将等式左右两边列矢量化，得到二次直积分解

（3）将步骤（1）中权矢量w的一次直积分解和步骤（2）中空域权矢量组u_n的二次直积分解进行截断处理，得到空时三维滤波器权矢量w的双重直积分解形式：D＜＜r₁，符号<<表示远小于；

（4）分别将d_n、和θ_n（n=1,…,D）列矢量化，形成列矢量 $p_d={[d_1^T,...,d_D^T]}^T,$ $p_{\mathrm{\&theta;}}={[{\mathrm{\&theta;}}_1^T,...,{\mathrm{\&theta;}}_D^T]}^T,$ 并利用双重双迭代算法，求解列矢量p_d、和p_θ；

（5）根据求得的列矢量p_d、和p_θ恢复出空时三维滤波器权矢量w，用来抑制杂波。

本发明与现有技术相比具有以下特点：

1、本发明与现有技术相比，样本需求量低，在小样本情况下，即可获得优良的性能。理论分析已证明，为保证自适应损益小于3dB，独立同分布的样本数应大于杂波协方差矩阵维数2倍。以扩展因子化方法（EFA）为例，该方法在空域进行联合自适应处理时，需要对维数为3MN的协方差矩阵进行求逆，样本需求量应大于6MN。本发明用D(D≤3)组空域和时域权矢量逼近空时三维滤波器权矢量，利用双重双迭代算法求解D组空域和时域权矢量，其中涉及的协方差矩阵维数分别为DM，DN和DK，为保证协方差矩阵非奇异，本发明方法需要的样本数只需不小于max{DM,DN,DK},max表示取括号中数的最大值。一般情况下，max{DM,DN,DK}<<6MN，本发明大大降低了样本需求量。

2、本发明与现有技术相比，降低了计算复杂度，有利于实时处理。以典型的扩展因子化方法（EFA）为例，对3MN维的协方差矩阵进行求逆，计算复杂度为O(27(MN)³)。本发明中涉及的协方差矩阵维数分别为DM，DN和DK，求逆运算的复杂度分别为O((DM)³)，O((DN)³)及O((DK)³)，本发明中外层双迭代算法2步即可收敛，内层双迭代算法5步即可收敛，则本发明总的计算复杂度为O(10(DM)³+10(DN)³+2(DK)³)。由于D取值很小（D<<3），通常O(10(DM)³+10(DN)³+2(DK)³)<<O(27(MN)³)，因此本发明大大降低了计算复杂度，有利于实时处理。

3、本发明与现有技术相比，在存在阵元误差时，具有更好的动目标多普勒频率检测性能，即具有更好的误差稳健性。以典型的EFA算法为例，由于存在阵元误差，杂波谱在空域发生扩散，在使用多普勒滤波器进行滤波时，有更多的杂波进入多普勒主瓣，造成空域自适应处理所形成的凹口宽度和深度难以满足杂波抑制要求，使得该算法性能下降。本发明在时域维、空域方位维和空域俯仰维进行双重迭代联合自适应处理，使得波束方向图在空域形成深凹口，进一步提高了杂波抑制能力。

附图说明

图1是机载雷达阵列模型图

图2是本发明的实现流程图

图3是本发明中双重双迭代算法的内层迭代收敛曲线

图4是本发明中双重双迭代算法的外层迭代收敛曲线

图5是本发明与EFA方法输出信噪比改善比较图

图6是本发明与EFA方法处理仿真数据时，剩余杂波归一化输出功率比较图

图7是本发明与EFA方法的输出信噪比改善随样本数的变化曲线图

具体实施方式

参照图2，本发明的具体实现步骤如下：

1.将空时三维滤波器权矢量w进行一次直积分解,将空时三维滤波器权矢量w=[w_1,1,1…w_M,1,1…w_1,2,1…w_M,2,1…w_1,N,K…w_M,N,K]^T排列成一个K×MN维的权矩阵W₁，利用满秩分解方法将权矩阵W₁分解为空域权矢量组和时域权矢量组的乘积形式，得到其中r₁为权矩阵W₁的秩，u_n表示空域权矢量组，d_n表示时域权矢量组；将权矩阵W₁列矢量化可恢复出w，则将等式左右两边列矢量化，可得一次直积分解

2.将空域权矢量组u_n进行二次直积分解，将空域权矢量组u_n排成M×N维空域权矩阵U_n，利用满秩分解方法将空域权矩阵U_n分解为空域方位维权矢量组和空域俯仰维权矢量组的乘积形式，得到其中r₂为权矩阵U_n的秩，θ_n表示空域俯仰维权矢量组，表示空域方位维权矢量组；将空域权矩阵U_n列矢量化可恢复出空域权矢量组u_n，则将等式左右两边都列矢量化，可得到二次直积分解形式

本发明应用的满秩分解方法是：对于秩为r的矩阵Z，可以分解为一个列满秩且秩为r的矩阵F与一个行满秩且秩为r的矩阵G的乘积形式，即有其中F=[f₁,...,f_r]，G=[g₁,...,g_r]^H，列矢量f_i(i=1,...,r)彼此不相关，列矢量g_i(i=1,...,r)彼此不相关。

3.对步骤1中得到的权矢量w的一次直积分解和步骤2中得到的空域权矢量组u_n的二次直积分解进行截断处理，先截取D个时域权矢量组和D个空域权矢量组u_n来表示权矢量w，得到再截取1个空域俯仰维权矢量组和1个空域方位维权矢量组θ_n来表示空域权矢量组u_n，得到将其带入得到权矢量w的双重直积分解形式：

4.分别将d_n、和θ_n（n=1,…,D）列矢量化，形成列矢量 $p_{\mathrm{\&theta;}}={[{\mathrm{\&theta;}}_1^T,...,{\mathrm{\&theta;}}_D^T]}^T,$ 然后利用双重双迭代算法求解列矢量p_d、和p_θ。

下面具体介绍求解过程。

假设x=[x_1,1,1…x_M,1,1…x_1,2,1…x_M,2,1…x_1,N,K…x_M,N,K]^T为空时三维接收采样数据矢量，目标时域导向矢量为a_d，目标空域方位维导向矢量为目标空域俯仰维导向矢量为a_θ，则目标导向矢量为

1）固定p_d，求p_θ的最优解。给定初值建立如下代价函数：

其中 $X=\mathrm{diag}(\mathrm{reshape}((d_1^H\&CircleTimes;I_{\mathrm{MN}})x,M,N),...,\mathrm{reshape}((d_D^H\&CircleTimes;I_{\mathrm{MN}})x,M,N)),$ 符号E表示求期望， $\mathrm{reshape}((d_i^H\&CircleTimes;I_{\mathrm{MN}})x,M,N),i=1,...,D$ 表示将矢量排成M×N矩阵，I_MN表示MN维单位矩阵，符号diag表示对角化；

2）利用双迭代算法求解式（2），得到和p_θ1。求解过程如下：

2.1）给定初始值对其进行归一化符号||||表示取模。

2.2）计算p_θ0

其中L表示样本数，X_l表示第l个空时三维采样数据。

2.3）计算

2.4）判断迭代停止公式是否成立，其中ε表示给定的迭代停止参数。若成立，则迭代终止；若不成立，将更新为然后重复执行步骤2.2）和2.3），直到本内层迭代终止，最后，输出最优权和p_θ1=p_θ0；

3）固定权矢量和p_θ1，求解权矢量p_d1；

令 建立如下代价函数：

$\left\{\begin{array}{c}\min E{\left|\right|p_{\mathrm{\&theta;}1}^H{\mathrm{Yp}}_{d1}\left|\right|}^2\\ s.t.p_{\mathrm{\&theta;}1}^H{\mathrm{bp}}_{d1}=1\end{array}\right.---\left(5\right)$

其中Y=diag(Y(1),…,Y(D))，利用拉格朗日乘子法求得：

$p_{d1}=\frac{{\left(\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}\left(Y_l^H{p}_{\mathrm{\&theta;}1}\right){\left(Y_l^H{p}_{\mathrm{\&theta;}1}\right)}^H\right)}^{-1}{b}^H{p}_{\mathrm{\&theta;}1}}{\left|\right|{\left(\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}\left(Y_l^H{p}_{\mathrm{\&theta;}1}\right){\left(Y_l^H{p}_{\mathrm{\&theta;}1}\right)}^H\right)}^{-1}{b}^H{p}_{\mathrm{\&theta;}1}\left|\right|}---\left(6\right)$

其中Y_i表示第l个训练样本；

4）判断迭代停止公式||p_d1-p_d0||<z,0<z<<1是否成立，其中z表示给定的迭代停止参数。若成立则迭代终止；若不成立，将p_d0更新为p_d1，重复执行步骤1）、2）、3）直到外层迭代终止。最后，输出最优解p_d、和p_θ。

5.由求得的变量p_d、和p_θ的最优解p_d1、和p_θ1，根据步骤3恢复出空时三维滤波器权矢量，用以抑制杂波。

本发明的效果通过以下仿真试验进一步说明：

1.仿真条件

实验采用矩形平面阵如图1所示，其中M=4表示阵元行数，N=8表示阵元列数，V=90m/s表示载机速度，R=500km表示射程，θ_t=90°表示目标方位角，表示目标俯仰角，ψ_t=90°表示目标锥面角，v_t=36m/s表示目标相对于载机的径向速度。载机高度8km，波长0.3m，阵元间距为0.15m，PRF为1200Hz，一个CPI内脉冲数为16，天线波束指向阵面法线方向，发射俯仰维和方位维均加30dB固定权，单个阵元输入杂噪比（CNR）为60dB,信噪比（SNR）为0dB，阵元幅相误差2%，仿真中采用2组双重直积分解逼近空时三维滤波器权矢量，即令D=2，其中D≤r₁，r₁表示权矩阵W₁的秩。

2.仿真内容

仿真1，在小样本时（样本数为100）时，本发明中双重双迭代算法中，内层迭代时输出信噪比改善随迭代次数的变化曲线图，如图3。

仿真2，在小样本时（样本数为100）时，本发明中双重双迭代算法中，外层迭代时输出信噪比改善随迭代次数的变化曲线图，如图4。

仿真3，在小样本时（样本数为100）时，本发明与EFA方法输出信噪比改善性能比较图，结果如图5。

仿真4，在小样本时（样本数为100）时，本发明与EFA方法处理仿真数据时，剩余杂波归一化输出功率比较图，如图6。

仿真5，本发明与EFA方法的输出信噪比改善随样本数的变化曲线比较图，如图7。

3.仿真分析

从图3可以看出，小样本情况下，在存在阵元误差时，本发明中双重双迭代算法的内层迭代经过5步即可收敛，从图4可以看出，本发明中双重双迭代算法的外层迭代经过2步即可收敛，则双重双迭代算法总共经过10步即可收敛，总的计算复杂度为O(10(8)³+10(16)³+2(32)³)。而EFA方法需要对96维的协方差矩阵进行求逆，计算复杂度为O(96³)。本发明的计算复杂度O(10(8)³+10(16)³+2(32)³)<<O(96³)，与EFA方法相比大大降低了计算复杂度，节约了计算量。

从图5可以看出，在小样本时（样本数为100），本发明与EFA方法相比，有5dB左右的性能改善。在100个样本的情况下，由于样本数小于EFA方法中协方差矩阵维数的2倍，即100<<192，此时EFA方法自适应损益较大，而本发明中协方差矩阵最高维数为32，(32×2)<<100，因此本发明自适应损益很小。此外，由于幅相误差的存在，EFA方法在时域进行多普勒滤波时，有较多杂波进入主瓣，使得算法性能下降，而本发明在时域、空域方位维、空域俯仰维进行迭代联合自适应处理，对误差具有更好的稳健性。因此，本发明具有更好的动目标检测性能。

从图6可以看出，在归一化多普率频率为0.4时，本发明与EFA方法相比，在各距离门有约5dB的性能改善，本发明具有更好的动目标检测性能。

从图7可以看出，本发明在64个样本时开始接近收敛，而EFA方法在192个样本时开始接近收敛，说明本发明在小样本情况下即可收敛，而EFA方法在大样本数情况下才会收敛。而在真实机载雷达快时变的杂波环境中，符合独立同分布的样本较少，本发明具有更好的实用性。

以上是本发明的优选实例，并不构成对本发明的任何限制，显然在不脱离本发明原理的情况下，均可用不同组双重直积分解逼近空时三维滤波器权矢量，但这些均在本发明的保护之列。

Claims (1)

1.一种基于双重直积分解的机载雷达空时三维杂波抑制方法，首先，将空时三维滤波器权矢量进行一次直积分解，得到空时三维滤波器权矢量关于时域权矢量组、空域权矢量组的直积构成形式；其次，将空域权矢量组进行二次直积分解，得到空域权矢量组关于空域方位维权矢量组、空域俯仰维权矢量组的直积构成形式；然后，截取多对时域权矢量组、空域方位维权矢量组及空域俯仰维权矢量组逼近空时三维滤波器权矢量；接着，利用双重双迭代算法求解截取的多对时域权矢量组、空域方位维权矢量组及空域俯仰维权矢量组；最后，由所求得的多对权矢量组恢复出空时三维滤波器权矢量，用来抑制杂波，具体实现过程如下：

假设机载雷达阵列为M×N矩形面阵，其中M表示面阵行数，N表示面阵列数，一个相干处理时间CPI内的脉冲数为K；

(1)对空时三维滤波器权矢量w进行一次直积分解

将空时三维滤波器权矢量w＝[w_1，1，1…w_M，1，1…w_1，2，1…w_M，2，1…w_1，N，K…w_M，N，K]^T排列成一个K×MN维的权矩阵W₁，其中符号T表示转置，利用满秩分解方法将权矩阵W₁分解为其中r₁表示权矩阵W₁的秩，u_n表示空域权矢量组，d_n表示时域权矢量组，符号H表示共轭转置；将等式左右两边列矢量化，得到一次直积分解其中符号*表示取共轭，符号表示直积；

(2)对空域权矢量组u_n进行二次直积分解

将空域权矢量组u_n排成M×N维矩阵U_n，利用满秩分解方法将权矩阵U_n分解为其中r₂表示U_n的秩，θ_n表示空域俯仰维权矢量组，表示空域方位维权矢量组；将等式左右两边列矢量化，得到二次直积分解

(3)将步骤(1)中得到的权矢量w的一次直积分解和步骤(2)中得到的空域权矢量组u_n的二次直积分解进行截断处理，得到空时三维滤波器权矢量w的双重直积分解形式：符号＜＜表示远小于；

(4)分别将d_n，和θ_n列矢量化，其中n＝1，…，D，形成列矢量 $p_d={[d_1^T,\&CenterDot;\&CenterDot;\&CenterDot;,d_D^T]}^T,$ $p_{\mathrm{\&theta;}}={[{\mathrm{\&theta;}}_1^T,\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{\&theta;}}_D^T]}^T,$ 并利用双重双迭代算法，求解列矢量p_d、和p_θ，所谓的双重双迭代算法其具体按如下步骤进行：

假设空时三维接收采样数据矢量x＝[x_1，1，1…x_M，1，1…x_1，2，1…x_M，2，1…x_1，N，K…x_M，N，K]^T，目标时域导向矢量为a_d，目标空域方位维导向矢量为目标空域俯仰维导向矢量为a_θ，则目标导向矢量为

①固定p_d，求p_θ的最优解，给定初值建立如下代价函数：

其中 $X=\mathrm{diag}(\mathrm{reshape}((d_1^H\&CircleTimes;I_{\mathrm{MN}})x,M,N),\&CenterDot;\&CenterDot;\&CenterDot;,\mathrm{reshape}((d_D^H\&CircleTimes;I_{\mathrm{MN}})x,M,N)),$ 符号E表示求期望，表示将矢量排成M×N矩阵，I_MN表示MN维单位矩阵，符号diag表示对角化；

②利用双迭代算法求解式(1)，得到p_θ1，求解过程如下：

[2.1]给定初始值对其进行归一化符号|| ||表示取模；

[2.2]计算p_θ0

其中L表示样本数，X_l表示第l个空时三维采样数据；

[2.3]计算

[2.4]判断迭代停止公式是否成立，其中ε表示给定的迭代停止参数，若成立，则迭代终止；若不成立，将更新为然后重复执行步骤[2.2]和[2.3]，直到本内层迭代终止；输出最优权和p_θ1＝p_θ0；

③固定权矢量和p_θ1，求解权矢量p_d1，

令 建立如下代价函数：

$\left\{\begin{array}{c}\min E{\left|\right|p_{\mathrm{\&theta;}1}^H{\mathrm{Yp}}_{d1}\left|\right|}^2\\ s.t.p_{\mathrm{\&theta;}1}^H{\mathrm{bp}}_{d1}=1\end{array}\right.---\left(4\right)$

其中Y＝diag(Y(1)，…，Y(D))，利用拉格朗日乘子法求得：

$p_{d1}=\frac{{\left(\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}\left(Y_l^H{p}_{\mathrm{\&theta;}1}\right){\left(Y_l^H{p}_{\mathrm{\&theta;}1}\right)}^H\right)}^{-1}{b}^H{p}_{\mathrm{\&theta;}1}}{\left|\right|{\left(\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}\left(Y_l^H{p}_{\mathrm{\&theta;}1}\right){\left(Y_l^H{p}_{\mathrm{\&theta;}1}\right)}^H\right)}^{-1}{b}^H{p}_{\mathrm{\&theta;}1}\left|\right|}---\left(5\right)$

其中Y_l表示第l个训练样本；

④判断迭代停止公式||p_d1-p_d0||＜z，0＜z＜＜1是否成立，其中z表示给定的迭代停止参数，若成立则迭代终止；若不成立，将p_d0更新为p_d1，重复执行步骤①、②、③直到外层迭代终止，最后，输出最优解p_d、和p_θ；

(5)根据求得的列矢量p_d、和p_θ恢复出空时三维滤波器权矢量w，进行杂波抑制。