一种空时自适应处理中的自适应权值迭代计算方法

技术领域

本发明属于机载相控阵雷达杂波抑制技术领域，涉及一种空时自适应处理中的自适应权值迭代计算方法。 

背景技术

机载相控阵雷达能实现对地面运动目标的有效检测，不过处于下视工作状态的机载相控阵雷达将面临比地基雷达更为严重的地/海杂波问题。地/海杂波不仅分布范围广、强度大，而且呈现出很强的空时耦合特性。空时自适应处理（STAP）技术能充分利用空域和时域信息，对目标信号进行相干积累的同时，将自适应空域处理和自适应多普勒处理两者优势结合起来，在空时域联合自适应滤除杂波，能获得更好的主瓣杂波抑制性能，改善对慢速目标的检测；同时可以对受旁瓣杂波干扰的小目标进行有效的检测。 

常规的全维STAP算法需要大量满足独立同分布（I.I.D）条件的训练样本对协方差矩阵进行估计，在非均匀杂波环境中此条件尤其得不到满足；而且当系统自由度很高时，全维协方差矩阵直接求逆运算（SMI）在现有计算水平下几乎无法实现。虽然降维和非均匀STAP算法能够降低STAP自适应权值计算的运算量和改善非均匀杂波环境中STAP算法杂波抑制性能，但是常规降维STAP算法在求解自适应权值时依然面临着对协方差矩阵直接求逆的运算，这将耗费系统很大的运算量和设备量，使得STAP技术难以满足系统实时性的要求。另外，虽然基于协方差矩阵逆更新的SMI算法无需估计采样协方差矩阵，其迭代计算的次数等于训练样本的数目，但是因为零矩阵不存在逆矩阵，该算法将难以设置初始逆矩阵来等价求解SMI算法的自适应权矢量，只能获得近似的解。 

发明内容

本发明针对常规STAP算法自适应权值计算需对空时协方差矩阵直接求逆，耗费系统很大运算量和设备量，使得STAP技术难以满足实时性要求的问题，根据协方差矩阵为Hermitian矩阵的特性，利用脉冲数据阶数分块递推，提出一种空时自适应处理中的自适应权值迭代计算方法。 

步骤一、建立空时接收数据模型； 

假设雷达天线阵元数目为N，发射脉冲数目为M，阵元间距为d，载机飞行速度为v，高度为h，脉冲重复频率PRF为f_r，T_r＝1/f_r为脉冲重复时间；若将斜距为R_c处的杂波距离环在 方位角度θ上分成N_c个间隔为Δθ＝2π/N_c的杂波散射单元，θ和 是杂波散射单元的方位角和俯仰角， 和f_t＝βf_s分别为归一化空间频率和多普勒频率，β＝2vT_r/d表示杂波谱的斜率；那么第i个杂波散射单元的N×1维空间导向矢量c(f_s,i)和M×1维时间导向矢量c(f_t,i)即表示为 

$\left(\begin{array}{c}c\left(f_{s,i}\right)={[1,\exp \left(j2\pi f_{s,i}\right),...,\exp \left(j2\pi (N-1)f_{s,i}\right)]}^T\\ c\left(f_{t,i}\right)={[1,\exp \left(j2\pi f_{t,i}\right),...,\exp \left(j2\pi (M-1)f_{t,i}\right)]}^T\end{array}\right)---\left(5\right)$

斜距R_c处的杂波回波为N_c个在空间上相互独立的杂波散射源响应之和 

$x_c={\Sigma }_{i=1}^{N_c}{a}_{i}c(f_{s,i},f_{t,i})---\left(6\right)$

其中 为第i个杂波散射单元的空时导向矢量，其中α_i(θ)为杂波散射单元的回波信号复幅度，其中α_i(θ)不仅取决于发射天线方向图，还跟杂波散射特性有关，建模为广义平稳随机过程 

$E\left\{a_i{a}_j^{*}\right\}=0,\forall i,j:i\ne j---\left(7\right)$

而且第i个杂波散射单元的平均强度假设为和发射天线的增益成正比 

E{|α_i|²}＝G_i,for i＝1,...,N_c                (8) 

其中G_i为正的实常数，正比与发射天线增益，那么第l个距离门内的空时接收数据表示为如下的向量形式 

x_k＝x_c，k+x_n，k

＝[x_1，k，x_2，k，...，x_M，k]^T             (9) 

其中x_n,k表示零均值的高斯白噪声，而x_m，k＝[x_1，m，k,x_2，m，k，...,x_N，m，k]^T表示第m个脉冲接收的N×1维的阵列数据； 

步骤二：估计STAP协方差矩阵；利用训练样本数据对NM×NM维协方差矩阵进行估计 

$\hat{R}=\frac{1}{L}{\Sigma }_{l=1}^L{x}_l{x}_l^H---\left(6\right)$

其中L是满足I.I.D条件的训练样本数目；协方差矩阵估计 为非负定Hermitian矩阵，假设存在足够数量的I.I.D训练样本，则 满秩为正定的Hermitian矩阵； 

步骤三：建立基于脉冲阶数的协方差矩阵迭代求逆模型； 

待检测距离单元的协方差矩阵表示为 而 按照脉冲阶数分解如下 

$R_l\left(M\right)=E\left(\left(\begin{array}{c}{x}_{1,l}\\ .\\ .\\ .\\ x_{m,l}\\ .\\ .\\ .\\ x_{M,l}\end{array}\right){\left(\begin{array}{c}{x}_{1,l}\\ .\\ .\\ .\\ x_{m,l}\\ .\\ .\\ .\\ x_{M,l}\end{array}\right)}^H\right)=\left(\begin{array}{cc}{R}_l(M-1)& F_l(M-1)\\ F_l^H(M-1)& G_l\left(M\right)\end{array}\right)---\left(7\right)$

其中矩阵 

$R_l(M-1)=E\left(\left(\begin{array}{c}{x}_{1,l}\\ .\\ .\\ .\\ x_{m,l}\\ .\\ .\\ .\\ x_{M-1,l}\end{array}\right)\left(\begin{array}{ccccccccc}{x}_{1,l}^H& .& .& .& x_{m,l}^H& .& .& .& x_{M-1,l}^H\end{array}\right)\right)---\left(8\right)$

表示前M-1个接收脉冲数据形成的N(M-1)×N(M-1)维的协方差矩阵， 

$F_l(M-1)=E\left(\left(\begin{array}{c}{x}_{1,l}{x}_{M,l}^H\\ x_{2,l}{x}_{M,l}^H\\ .\\ .\\ .\\ x_{m,l}{x}_{M,l}^H\\ .\\ .\\ .\\ x_{M-1,l}{x}_{M,l}^H\end{array}\right)\right)---\left(9\right)$

表示第M个脉冲接收数据和前M-1个脉冲接收数据的N(M-1)×N维的互相关矩阵及 

$G_l\left(M\right)=E\left(x_{M,l}{x}_{M,l}^H\right)---\left(10\right)$

表示第M个脉冲接收数据的N×N维协方差矩阵； 

待检测距离门的M个脉冲接收数据形成的空时协方差矩阵R由前M-1个脉冲接收数据形成的空时协方差矩阵表示，即只要得到第1个脉冲的阵列接收数据，即按照脉冲阶数进行递归计算空时协方差矩阵的逆，进而得到空时自适应权值；利用Hermitian矩阵分块和脉冲递推的特性，计算得到前m个脉冲接收数据的协方差矩阵R_l(m)的逆与前m-1个脉冲接收数据的协方差矩阵R_l(m-1)的逆之间的迭代关系 

$R^{-1}\left(m\right)={\left(\begin{array}{cc}R(m-1)& F(m-1)\\ F^H(m-1)& G\left(m\right)\end{array}\right)}^{-1}$(11) 

$=\left(\begin{array}{cc}{R}^{-1}(m-1)+B\left(m\right)P^{-1}\left(m\right)B^H\left(m\right)& B\left(m\right)P^{-1}\left(m\right)\\ P^{-1}\left(m\right)B^H\left(m\right)& P^{-1}\left(m\right)\end{array}\right)$

其中矩阵B(m)＝-R^-1(m)F(m-1)，矩阵P(m)＝G(m)-F^H(m-1)R^-1(m)F(m-1)； 

步骤三中所述的按照脉冲阶数进行递归计算空时协方差矩阵的逆采用以下方法： 

在接收到第1个脉冲数据时，利用距离训练样本计算协方差矩阵R(1)，因为该矩阵是Hermitian矩阵，利用Hermitian矩阵的顺序主子矩阵迭代完成第1个脉冲协方差矩阵的求逆；接着在接收第2个直至第M个脉冲数据时，便利用以上迭代方法计算前2个直至前M个脉冲接收数据协方差矩阵的逆；递推过程中，中间变量矩阵P(2)直至矩阵P(M)的逆，同样根据Hermitian矩阵的性质，利用Hermitian矩阵的顺序主子矩阵迭代完成上述N×N维的变量矩阵逆的计算。 

步骤四、计算STAP自适应权矢量；STAP处理自适应权矢量通过下述线性约束优化问题得到 

w＝R^-1a(f_s0,f_t0)             (12) 

其中a(f_s0,f_t0)表示目标空时导向矢量，而第l个距离单元的滤波输出为 

y_l＝w^Hx_l           (13) 

其中H表示共轭转置运算，x_l表示待检测距离单元数据；自此，就完成了一种用于空时自适应处理自适应权值迭代的计算方法。 

本发明提出的一种用于空时自适应处理自适应权值迭代的计算方法，对比已有技术，避免了协方差矩阵直接求逆的运算，因此大大降低了计算STAP自适应权值的运算量，更利于工程实现，其效果具体如下： 

①本发明避免了协方差矩阵直接求逆的运算，能大大降低了计算STAP自适应权值的运算量； 

②本发明中基于脉冲阶数迭代计算的方法，可以任意选择进行STAP自适应权矢量计算的脉冲接收数据，能够在保持性能的基础上进一步降低计算复杂度； 

③本发明不需要设置初始逆矩阵来等价求解自适应权矢量，因此能够获得STAP自适应权矢量的精确解。 

附图说明

图1为机载相控阵雷达几何结构示意图； 

图2为计算复杂度比较； 

图3为空时自适应方向图比较； 

其中图(a)为直接矩阵求逆，(b)为Hermitian矩阵求逆，(c)为基于脉冲阶数迭代求逆； 

图4为实测MCARM数据距离-多普勒输出结果比较。 

其中图(a)为直接矩阵求逆，(b)为Hermitian矩阵求逆，(c)为基于脉冲阶数迭代求逆。 

具体实施方式

下面结合附图对本发明方法的实施方式做详细说明。 

为了更方便地描述本发明的内容，首先对利用Hermitian矩阵的顺序主子矩阵进行迭代计算协方差矩阵的逆作如下解释： 

假设矩阵R为D×D维的Hermitian矩阵，且矩阵R_d+1表示R的第d+1阶顺序主子矩阵，即R_d+1＝R(1:d+1,1:d+1)。根据分块Hermitian矩阵特性，矩阵R_d+1的逆可以利用R_d的逆来计算。由于矩阵R_d+1的逆矩阵Q_d+1也是Hermitian矩阵，则 

$Q_{d+1}=\left(\begin{array}{cc}{Q}_d& q_{d+1}\\ q_{d+1}^H& q_{d+1}\end{array}\right)---\left(1\right)$

其中q_d+1表示Q_d+1的第d+1个对角元素，即q_d+1＝Q_d+1(d+1,d+1)；q_d+1表示矩阵Q_d+1的第d+1列前d个元素组成的列矢量，即q_d+1＝Q_d+1(1:d,d+1)；Q_d是Q_d+1的第d阶顺序主子矩阵，即Q_d＝Q_d+1(1:d,1:d)。根据互为逆矩阵乘积为单位矩阵可得 

$R_{d+1}{Q}_{d+1}=\left(\begin{array}{cc}{R}_d& r_{d+1}\\ r_{d+1}^H& {\rho }_{d+1}^b\end{array}\right)\left(\begin{array}{cc}{Q}_d& q_{d+1}\\ q_{d+1}^H& q_{d+1}\end{array}\right)=\left(\begin{array}{cc}{I}_d& 0_{d+1}\\ 0_{d+1}^H& 1\end{array}\right)---\left(2\right)$

其中0_d+1是d×1维的零向量。通过计算，可以得到如下迭代公式 

$R_{d+1}^{-1}={\left(\begin{array}{cc}{R}_d& r_{d+1}\\ r_{d+1}^H& {\rho }_{d+1}\end{array}\right)}^{-1}=\left(\begin{array}{cc}{R}_d^{-1}& 0_{d+1}\\ 0_{d+1}^H& 0\end{array}\right)+\frac{1}{{\alpha }_{d+1}}\left(\begin{array}{cc}{b}_{d+1}{b}_{d+1}^H& b_{d+1}\\ b_{d+1}^H& 1\end{array}\right)---\left(3\right)$

其中向量b_d+1和系数α_d+1定义如下 

$b_{d+1}=-R_d^{-1}{r}_{d+1}$(4) 

${\alpha }_{d+1}={\rho }_{d+1}-r_{d+1}^H{R}_d^{-1}{r}_{d+1}={\rho }_{d+1}+r_{d+1}^H{b}_{d+1}$

其中ρ_d+1＝R_d+1(d+1,d+1)和r_d+1＝R_d+1(1:d,d+1)。 

利用Hermitian矩阵的顺序主子矩阵迭代计算协方差矩阵逆的具体步骤如下： 

*说明：Re(·)为取实部算子。通过取实部运算可以有效避免计算误差造成算法失效的影响，得到稳健的递推计算过程。 

一种用于空时自适应处理自适应权值迭代的计算方法，该方法实现的过程如下： 

步骤一、接收空时杂波数据模型 

①空时接收数据建模 

根据如图1所示的机载相控阵雷达几何结构，假设雷达天线阵元数目为N，发射脉冲数目为M，阵元间距为d。载机飞行速度为v，高度为h。脉冲重复频率（PRF）为f_r，T_r＝1/f_r为脉冲重复时间。若将斜距为R_c处的杂波距离环在方位角度θ上分成N_c个间隔为Δθ＝2π/N_c的杂波散射单元。θ和 是杂波散射单元的方位角和俯仰角。 和f_t＝βf_s分别为归一化空间频率和多普勒频率，β＝2vT_r/d表示杂波谱的斜率。那么第i个杂波散射单元的N×1维空间导向矢量c(f_s,i)和M×1维时间导向矢量c(f_t,i)可以表示为 

$\left(\begin{array}{c}c\left(f_{s,i}\right)={[1,\exp \left(j2\pi f_{s,i}\right),...,\exp \left(j2\pi (N-1)f_{s,i}\right)]}^T\\ c\left(f_{t,i}\right)={[1,\exp \left(j2\pi f_{t,i}\right),...,\exp \left(j2\pi (M-1)f_{t,i}\right)]}^T\end{array}\right)---\left(14\right)$

斜距R_c处的杂波回波为N_c个在空间上相互独立的杂波散射源响应之和 

$x_c={\Sigma }_{i=1}^{N_c}{a}_{i}c(f_{s,i},f_{t,i})---\left(15\right)$

其中 为第i个杂波散射单元的空时导向矢量，其中α_i(θ)为杂波散射 单元的回波信号复幅度。那么第l个距离门内的空时接收数据可以表示为如下的向量形式 

x_k＝x_c,k+x_n,k

                               (16) 

＝[x_1,k,x_2,k,...,x_M,k]^T

其中x_n,k表示零均值的高斯白噪声，而x_m,k＝[x_1，m，k,x_2，m，k，...,x_N，m，k]^T表示第m个脉冲接收的N×1维的阵列数据。 

②STAP协方差矩阵估计 

STAP处理中，待检测距离门的杂波特性是未知的，需要利用训练样本数据对NM×NM维协方差矩阵进行估计 

$\hat{R}=\frac{1}{L}{\Sigma }_{l=1}^L{x}_l{x}_l^H---\left(17\right)$

其中L是满足I.I.D条件的训练样本数目。显然，协方差矩阵估计 为非负定Hermitian矩阵。假设存在足够数量的I.I.D训练样本，则 满秩为正定的Hermitian矩阵。 

步骤二、基于脉冲阶数的STAP自适应权值迭代计算 

①基于脉冲阶数的协方差矩阵迭代求逆模型 

由式（7）可以看出，待检测距离门的M个脉冲接收数据形成的空时协方差矩阵R可以由前M-1个脉冲接收数据形成的空时协方差矩阵表示，即只要得到第1个脉冲的阵列接收数据，就可以不断的按照脉冲阶数进行递归计算空时协方差矩阵的逆，进而得到空时自适应权值。通过利用Hermitian矩阵分块和脉冲递推的特性，可以计算得到前m个脉冲接收数据的协方差矩阵R_l(m)的逆与前m-1个脉冲接收数据的协方差矩阵R_l(m-1)的逆之间的迭代关系 

$R^{-1}\left(m\right)={\left(\begin{array}{cc}R(m-1)& F(m-1)\\ F^H(m-1)& G\left(m\right)\end{array}\right)}^{-1}$(18) 

$=\left(\begin{array}{cc}{R}^{-1}(m-1)+B\left(m\right)P^{-1}\left(m\right)B^H\left(m\right)& B\left(m\right)P^{-1}\left(m\right)\\ P^{-1}\left(m\right)B^H\left(m\right)& P^{-1}\left(m\right)\end{array}\right)$

其中矩阵B(m)＝-R^-1(m)F(m-1)，矩阵P(m)＝(G(m)-F^H(m-1)R^-1(m)F(m-1))。 

②基于脉冲阶数的协方差矩阵迭代求逆具体步骤 

由式（11）可以看出，本发明提出的方法在接收到第1个脉冲数据时，利用距离训练样本计算协方差矩阵R(1)，因为该矩阵是Hermitian矩阵，可以利用Hermitian矩阵的顺序主子矩阵迭代完成第1个脉冲协方差矩阵的求逆。接着，在接收第2个直至第M个脉冲数据时，便可以利用以上迭代方法计算前2个直至前M个脉冲接收数据协方差矩阵的逆。不过在递推 的时候，需要求中间变量矩阵P(2)直至矩阵P(M)的逆，但是同样根据Hermitian矩阵的性质，可以利用Hermitian矩阵的顺序主子矩阵迭代完成上述N×N维的变量矩阵逆的计算。 

基于脉冲阶数的协方差矩阵迭代求逆的具体步骤如下： 

步骤三、计算STAP自适应权矢量 

STAP处理自适应权矢量可以通过下述线性约束优化问题得到 

w＝R^-1a(f_s0,f_t0)              (19) 

其中a(f_s0,f_t0)表示目标空时导向矢量，而第l个距离单元的滤波输出为 

y_l＝w^Hx_l               (20) 

自此，就完成了一种用于空时自适应处理自适应权值迭代的计算方法。本发明提出的方法避免了协方差矩阵直接求逆的运算，能够大大降低计算STAP自适应权值的运算量，因此更利于工程实现。 

为了验证本发明给出的STAP自适应权值迭代计算方法的性能，进行了如下仿真验证。首先比较本发明提出方法和常规STAP算法的计算复杂度，其次分别利用仿真数据和实测MCARM数据对其性能进行比较分析。 

实验一、计算复杂度比较 

假设STAP算法处理的是N×M×K维的空时数据立方体，其中N表示阵元数目，M表示CPI内的脉冲数目，K表示不模糊的距离门数目，直接协方差矩阵求逆SMI算法估计空时协方差矩阵的训练样本数目为K0。在全维STAP算法中，而言，估计NM×NM维的协方差矩阵需要的运算量为K₀(NM)²，滤波操作需要的运算量为KNM。直接协方差矩阵求逆的运算量是最大的，为(NM)³，其次是Hermitian协方差矩阵直接迭代求逆，需要的运算量大约为2(NM)³/3+NM/3，本发明提出的自适应权值迭代计算方法中协方差矩阵求逆需要的运算量大约为(NM)³/3+K₀(NM)²/2+K₀N²+2N³/3。将上述三种算法的运算量总结如下表所示。 

为了更直观的比较上述算法的运算量，假设阵元数目为N＝14，CPI内的脉冲数目为M＝12，用于协方差矩阵估计的训练数目为K₀=2NM，那么直接协方差矩阵求逆SMI算法、Hermitian协方差矩阵求逆SMI算法和基于脉冲阶数迭代协方差矩阵求逆SMI算法的运算量比较如图2所示。 

由图可以明显观察出，直接协方差矩阵求逆SMI算法的运算量最大，而Hermitian协方差矩阵求逆SMI算法的运算量仅次于直接矩阵求逆，利用Hermitian矩阵特性直接迭代计算并没有降低多少运算量。相比较而言，本发明提出的自适应权值迭代计算方法比前两种算法都低的多，递归的脉冲阶数越少，运算量越低，不过即使所有的脉冲数都利用，运算量也只是直接协方差矩阵求逆SMI算法的50%，这对于实际工程应用是十分有利的。 

实验二、仿真数据性能比较 

该实验对机载相控阵雷达杂波回波仿真，仿真的参数如下表所示。 

分别使用直接协方差矩阵求逆SMI算法、Hermitian协方差矩阵求逆SMI算法和基于脉冲阶数迭代协方差矩阵求逆SMI算法对上述产生的杂波数据进行STAP处理，各种算法得到的空时自适应方向图结果如图3所示。 

由图可以明显观察出，几种协方差矩阵求逆算法的空时自适应方向图在杂波方向都形成深的凹口，因此能有效抑制杂波检测出目标信号。另外，由于这些算法都能得到相同的自适应权矢量，因此空时自适应方向图都彼此相同。因此本发明提出的自适应权值迭代计算方法可以在保持性能的同时进一步降低协方差矩阵求逆的运算量，具有很好的工程实际意义。 

实验三、实测MCARM数据性能比较 

下面上述三种协方差矩阵求逆算法对实测MCARM数据进行STAP处理，得到的距离-多普勒输出结果如图4所示。 

由图可以明显观察出，上述三种协方差矩阵求逆算法对MCARM数据进行处理之后，距离-多普勒输出的目标峰值和背景噪声功率都基本一致，并且都能在归一化多普勒频率为-0.15，距离门为299号的位置处检测到目标信号。但是，本发明提出的自适应权值迭代计算方法能在保持目标检测性能的基础上降低计算复杂度，只有直接矩阵求逆算法的50%左右，因此更利于工程实现。