CN109212500A - 一种基于稀疏重构的ka-stap杂噪协方差矩阵高精度估计方法 - Google Patents

Abstract

本发明公布了一种基于稀疏重构的KA‑STAP杂噪协方差矩阵高精度估计方法，包括如下步骤：S1、分析机载雷达回波数据，利用稀疏重构获取高分辨率二维空时谱；S2、筛选二维空时谱上的像素点，并计算像素点对应的加权权值；S3、利用加权最小二乘法对杂波轨迹进行拟合；S4、根据稀疏重构空时谱估计噪声功率，并构造先验杂波加噪声协方差矩阵；S5、利用先验杂噪协方差矩阵和降维STAP对检测单元进行自适应滤波与目标检测。本发明方法相比较常规STAP算法，运算量较小，而且可以有效提升STAP系统在非平稳杂波环境下的的杂波抑制与目标检测性能，易于工程实施。

Description

一种基于稀疏重构的KA-STAP杂噪协方差矩阵高精度估计 方法

技术领域

本发明属于知识辅助空时自适应处理(KA-STAP)技术领域，尤其涉及一种基于稀疏重构的KA-STAP杂噪协方差矩阵高精度估计方法。

背景技术

STAP(Space-Time Adaptive Processing,空时自适应处理)是当前机载雷达抑制地面杂波、实现地面动目标探测的一种重要技术手段。常规STAP方法是基于相邻距离单元的杂波样本与待检测单元内的杂波满足统计平稳性这一假设，利用杂波协方差矩阵的最大似然估计求解自适应权值。为了保证相对最优STAP处理的输出信杂噪比(SCNR)损失限制在3dB范围内，用以估计杂波协方差的训练样本需满足独立同分布(IID)的条件且数量应当超过两倍以上的自适应处理器维数。

然而在实际应用中,机载雷达所处的环境较为复杂，无法满足统计平稳性的假设，难以获取足够数量的IID样本，这使得常规STAP方法性能下降。针对这一问题，国内外学者就KA-STAP展开研究，KA-STAP利用雷达参数、地形信息、数字地图及地表覆盖数据等先验信息辅助滤波器的设计，从而有效提高STAP技术在非均匀杂波环境下的杂波抑制性能。例如，利用数字地形高程数据(DTED)、地面覆盖/地面使用数据(LCLU)等先验信息选取均匀的样本，进而估计检测单元的干扰协方差矩阵，将先验知识和观测数据融合在参数化模型中以捕捉待测距离单元的瞬时特征，并借助先验信息构造先验的杂波加噪声协方差矩阵，通过线性加权的方式将样本估计协方差矩阵与先验协方差矩阵同时用于滤波器的设计，这种方法可以有效的提高非均匀杂波环境下的杂波抑制性能。

发明内容

针对机载非正侧视阵雷达近程杂波非平稳性引起的STAP性能下降的问题，本发明提出了一种基于稀疏重构的KA-STAP杂噪协方差矩阵高精度估计方法。

为解决上述技术问题，本发明采用了如下技术手段：

一种基于稀疏重构的KA-STAP杂噪协方差矩阵高精度估计方法，具体包括如下步骤：

S1，分析机载雷达回波数据，利用稀疏重构获取高分辨率二维空时谱；

S2，筛选二维空时谱上的像素点，并计算像素点对应的加权权值；

S3，利用加权最小二乘法对杂波轨迹进行拟合；

S4，根据稀疏重构空时谱估计噪声功率，并构造先验杂波加噪声协方差矩阵；

S5，利用先验杂噪协方差矩阵和降维STAP对检测单元进行自适应滤波与目标检测。

进一步的，所述的步骤S1具体通过以下步骤实现：

S11，对机载雷达接受到的回波信号进行加窗处理，处理后的接收信号T_X_l满足公式：

T_X_l＝X_l·T_w＝[S_{Tl_1}…S_{Tl_i}…S_{Tl_K}]_N×K (1)

其中，X_l是第l个距离单元的接收信号，T_w是以窗函数T_f＝[w(1),w(2),…,w(i),…,w(K)]为对角元素的K×K维加窗矩阵，w(i)是窗函数系数，S_{Tl_i}是第l个距离单元第i个脉冲各阵元加窗处理后的输出信号，l＝1,2,…,M，i＝1,2,…,K，M是距离单元总数，N是雷达天线的阵元个数，K是一次相干处理间隔内的时域脉冲数。

S12，利用FFT域变换矩阵F_D获得阵元-多普勒域输出信号：

D_X_l＝T_X_l·F_D＝[S_{Dl_1}…S_{Dl_i}…S_{Dl_K}]_N×K (2)

其中，S_{Dl_i}表示第l个距离单元第i个多普勒单元各阵元输出信号。

S13，利用凸优化算法对S_{Dl_i}进行稀疏重构，稀疏重构的约束方程如下：

其中，σ_i是第l个距离单元第i个多普勒单元接收信号的空域幅度分布，是函数g(x)取最小值时所对应的变量x的值，观测矩阵ψ_i是空域导向矢量构成的一组超完备基，大小是N×N_s，N_s是空域量化单元数，||·||₁表示L₁范数运算，ε_i是允许误差。

S14，对阵元-多普勒域数据逐多普勒单元进行空域稀疏重构，获得高分辨率的二维空时谱。

进一步的，所述的步骤S2具体通过以下步骤实现：

S21，在二维空时谱上挑选出每个归一化空间频率上幅值最大的像素点并记录其坐标。

S22，将S21挑选出来的像素点中偏离杂波轨迹的噪声像素点剔除，得到筛选后的点集{(x_j,y_j)}，x_j是第j个像素点的横坐标，y_j是第j个像素点的纵坐标，j＝1,2,3...N_I，N_I是筛选出的符合要求的像素点的个数。

S23，对空时谱幅值向量A进行归一化处理，得到归一化后的空时谱幅值向量a：

其中，是第j个像素点对应的真实幅度值，用归一化后的空时谱幅值向量a＝[a₁,a₂,...,a_NI]^T中a_j作为第j个像素点对应的加权权值。

进一步的，所述的步骤S3利用加权最小二乘法对杂波轨迹进行拟合，对应的目标函数如下：

其中，表示求解使函数f(p₁,p₂)取得最小值时ψ的值， λ是机载雷达波长，f_r是脉冲重复频率，v是载机飞行速度，ψ是载机偏航角，是杂波散射体的俯仰角。

进一步的，所述的步骤S4中的先验杂波加噪声协方差矩阵满足下列公式：

其中，

δ_n ²是根据稀疏重构空时谱估计得到噪声功率，ψ_opt是根据加权最小二乘法解得ψ的最优解，I表示NK×NK维的单位阵。

进一步的，所述的步骤S5具体通过以下步骤实现：

S51，利用NK×r维列满秩降维变换矩阵T对采样数据X_l进行降维处理，r是降维后的自适应处理自由度。

S52，计算降维后的目标导引矢量S_T、数据矢量X_lT和杂噪协方差矩阵R_0T：

S_T＝T^HS (8)

X_lT＝T^HX_l (9)

R_0T＝T^HR₀T (10)

其中，S是没有降维的目标导引矢量，X_l是第l个距离单元的接收信号，R₀是先验杂波加噪声协方差矩阵。

S53，计算降维处理后的最优自适应权值w_T：

其中，μ_T是归一化复常数，

S54，根据最优自适应权值和回波数据矢量计算二维滤波输出Z_out：

采用以上技术手段后可以获得以下优势：

本发明公开的一种基于稀疏重构的KA-STAP杂噪协方差矩阵高精度估计方法，由于机载非正侧视阵雷达的近程回波具有严重的距离依赖性，若采用距离采样统计平均的方法估计杂波协方差矩阵，会使得杂波协方差矩阵估计的精确度急剧下降，本发明利用稀疏重构获取高分辨率二维空时谱作为先验知识，对杂波轨迹进行拟合，构造杂波加噪声协方差矩阵用于STAP权值计算及自适应滤波，相比较常规STAP算法，运算量较小，而且可以有效提升STAP系统在非平稳杂波环境下的的杂波抑制与目标检测性能，易于工程实施。

附图说明

图1是本发明一种基于稀疏重构的KA-STAP杂噪协方差矩阵高精度估计方法的流程图。

图2是本发明方法中机载雷达构型和杂波几何关系图。

图3是本发明方法稀疏重构二维空时平均谱图。

图4是本发明方法杂波轨迹拟合曲线图。

图5是利用先验杂波加噪声协方差矩阵计算权值后的滤波输出图。

图6是利用训练样本估计协方差矩阵计算权值后的滤波输出图。

其中，图2中H是载机飞行高度，ψ是偏航角，α是散射体P的杂波入射锥角，β是散射体P与天线轴向的夹角，θ是散射体P的方位角，散射体P的俯仰角，v是载机飞行速度。

具体实施方式

下面结合附图对本发明的技术方案作进一步说明：

如图1所示，一种基于稀疏重构的KA-STAP杂噪协方差矩阵高精度估计方法，具体包括如下步骤：

S1，分析机载雷达回波数据，利用稀疏重构获取高分辨率二维空时谱；

S2，筛选二维空时谱上的像素点，并计算像素点对应的加权权值；

S3，利用加权最小二乘法对杂波轨迹进行拟合；

S4，根据稀疏重构空时谱估计噪声功率，并构造先验杂波加噪声协方差矩阵；

S5，利用先验杂噪协方差矩阵和降维STAP对检测单元进行自适应滤波与目标检测。

对机载非正侧视阵雷达杂波特性进行分析，如图2所示，机载雷达天线是均匀线阵，天线阵元个数为N，一次相干处理间隔内时域脉冲数为K，阵元间距为d，雷达波长为λ，阵元间距为d，雷达波长为λ，d＝λ/2。载机沿X轴以速度v匀速飞行，偏航角为ψ，飞行高度为H，α、β分别为散射体P的杂波入射锥角和散射体P与天线轴向的夹角，θ、分别为散射体P的方位角及俯仰角。

设定机载雷达某一距离环内杂波可按方位角等分为N_c个杂波源，是第i个杂波源的归一化空间频率，是归一化多普勒频率，分别满足以下公式：

机载雷达平台的运动会导致地杂波的空时二维耦合特性，使得散射体P的杂波入射锥角余弦cosα和散射体P与天线轴向的夹角余弦cosβ存在相互联系，杂波空时分布轨迹满足公式：

将上式化为关于归一化多普勒频率f_d与归一化空间频率f_s的方程：

当0°＜ψ＜90°时，根据上式可以得出其杂波轨迹是一簇斜椭圆，当然，在实际机载雷达天线中，由于天线后向辐射的影响较小，所以实际杂波轨迹为椭圆的一半。

将第l个距离单元的接收信号X_l逐脉冲排列可得：

X_l＝[S_{l_1}S_{l_2}…S_{l_K}]_N×K (16)

其中，S_{l_i}是第l个距离单元第i个脉冲各阵元接收信号矢量，l＝1,2,…,M，i＝1,2,…,K，M是距离单元总数。

选取窗函数T_f＝[w(1),w(2),…,w(i),…,w(K)]来进行信号加窗处理，w(i)是窗函数系数，构建加窗矩阵T_w如下：

加窗处理后的接受信号T_X_l满足公式：

T_X_l＝X_l·T_w＝[S_{Tl_1}…S_{Tl_i}…S_{Tl_K}]_N×K (18)

其中，S_{Tl_i}是第l个距离单元第i个脉冲各阵元加窗处理后的输出信号。

利用FFT域变换矩阵F_D获得阵元-多普勒域输出信号：

D_X_l＝T_X_l·F_D＝[S_{Dl_1}…S_{Dl_i}…S_{Dl_K}]_N×K (19)

其中，S_{Dl_i}表示第l个距离单元第i个多普勒单元各阵元输出信号。

利用凸优化算法对S_Dl__i进行稀疏重构，稀疏重构的约束方程如下：

其中，σ_i是第l个距离单元第i个多普勒单元接收信号的空域幅度分布，是函数g(x)取最小值时所对应的变量x的值，观测矩阵ψ_i是空域导向矢量构成的一组超完备基，大小是N×N_s，N_s是空域量化单元数，||·|₁表示L₁范数运算，ε_i是允许误差。

利用公式(20)对阵元-多普勒域数据逐多普勒单元进行空域稀疏重构，可以获得高分辨率的二维空时谱，记作A(f_d,f_s)，其中f_d、f_s分别是杂波散射单元的归一化多普勒频率和归一化空间频率。

在二维空时谱上挑选出每个归一化空间频率上幅值最大的像素点并记录其坐标，将挑选出来的像素点中偏离杂波轨迹的噪声像素点剔除，得到筛选后的点集{(x_j,y_j)}，x_j是第j个像素点的横坐标，y_j是第j个像素点的纵坐标，j＝1,2,3...N_I，N_I是筛选出的符合要求的像素点的个数。

将筛选出的像素点集合中每一个像素点对应空时谱上的杂波幅值记为A_j,利用恢复系数校正加窗效应，可以得到校正后的幅值A'_j为：

A_j'＝A_j/w (21)

由于FFT变换会导致信号幅值在频域产生相干积累，所以像素点对应的真实幅度值为：

则有时谱幅值向量

对空时谱幅值向量A进行归一化处理，可以得到归一化后的空时谱幅值向量a：

用归一化后的空时谱幅值向量中的a_j作为第j个像素点对应的加权权值。

利用加权最小二乘法对杂波轨迹进行拟合，杂波轨迹最小二乘拟合对应的目标函数如下：

其中，表示求解使函数f(p₁,p₂)取得最小值时ψ的值， f_r是脉冲重复频率，v是载机飞行速度，ψ是载机偏航角，是杂波散射体的俯仰角，可以通过雷达测高仪精确测量得到。

为了使函数f(p₁,p₂)取极小值，必须有即：

其中，分别表示函数f(p₁,p₂)相对于p₁、p₂的偏导数。

联立公式(25)和(26)，可以得到p₁的最优解：

进而得到ψ的最优解ψ_opt为：

根据稀疏重构空时谱估计出噪声功率δ_n ²，则先验杂波加噪声协方差矩阵满足公式：

其中，

I表示NK×NK维的单位阵。

利用NK×r维列满秩降维变换矩阵T对采样数据进行降维处理，r是降维后的自适应处理自由度。降维后的目标导引矢量S_T、数据矢量X_lT和杂噪协方差矩阵R_0T分别为：

S_T＝T^HS (31)

X_lT＝T^HX_l (32)

R_0T＝T^HR₀T (33)

其中，S是没有降维的目标导引矢量，X_l是第l个距离单元的接收信号，R₀是先验杂波加噪声协方差矩阵。

进而计算出降维处理后的最优自适应权值w_T：

其中，μ_T是归一化复常数，

根据最优自适应权值和回波数据矢量可以得到二维滤波输出Z_out：

下面通过仿真实验进一步验证本发明方法的有效性，机载雷达系统仿真实验的参数如表1所示。

X_l为第360个距离单元接收信号，其斜距为10.8km。在检测单元注入待检测目标，信杂噪比为-10dB，归一化多普勒频率与归一化空间频率分别为0.0313与0.2734。对待检测单元及其相邻的两个保护单元进行稀疏重构，并利用稀疏重构二维空时平均谱作为先验知识，根据杂波协方差矩阵估计算法进行杂波轨迹拟合，图3、图4分别为稀疏重构二维空时平均谱及杂波轨迹拟合曲线图。

利用先验杂波加噪声协方差矩阵计算STAP权值后进行目标检测的滤波输出如图5所示，利用训练样本估计协方差矩阵计算STAP权值后进行目标检测的滤波输出如图6所示，训练样本是与待检测单元相邻的120个距离单元，目标检测方法采用的多普勒三通道联合自适应处理(3DT-STAP)算法。对比图5和图6，由于近程杂波的距离非平稳性，各距离单元数据不再满足IID条件，导致直接由训练样本估计得到的协方差矩阵误差急剧加大，杂波抑制性能严重下降，目标无法被有效检测，而利用先验杂波协方差矩阵进行权值计算并滤波输出后，目标被有效探测，本发明方法在近程杂波区的杂波抑制性能非常明显。

表1雷达系统参数

参数名称	参数数值
		阵元个数	16
距离单元总数	1000
		一次相干积累脉冲数	128
脉冲重复频率	5000HZ
		采样带宽	5MHZ
载机速度	130m/s
		载机飞行高度	8000m
阵元间距与波长比值	1/2

上面结合附图对本发明的实施方式作了详细地说明，但是本发明并不局限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。

Claims (6)

1.一种基于稀疏重构的KA-STAP杂噪协方差矩阵高精度估计方法，其特征在于，包括如下步骤：

S1、分析机载雷达回波数据，利用稀疏重构获取高分辨率二维空时谱；

S2、筛选二维空时谱上的像素点，并计算像素点对应的加权权值；

S3、利用加权最小二乘法对杂波轨迹进行拟合；

S4、根据稀疏重构空时谱估计噪声功率，并构造先验杂波加噪声协方差矩阵；

S5、利用先验杂噪协方差矩阵和降维STAP对检测单元进行自适应滤波与目标检测。

2.根据权利要求1所述的一种基于稀疏重构的KA-STAP杂噪协方差矩阵高精度估计方法，其特征在于，所述的步骤S1具体通过以下步骤实现：

S11，对机载雷达接受到的回波信号进行加窗处理，处理后的接收信号T_X_l满足公式：

T_X_l＝X_l·T_w＝[S_{Tl_1}…S_{Tl_i}…S_{Tl_K}]_N×K

其中，X_l是第l个距离单元的接收信号，T_w是以窗函数T_f＝[w(1),w(2),…,w(i),…,w(K)]为对角元素的K×K维加窗矩阵，w(i)是窗函数系数，S_{Tl_i}是第l个距离单元第i个脉冲各阵元加窗处理后的输出信号，l＝1,2,…,M，i＝1,2,…,K，M是距离单元总数，N是雷达天线的阵元个数，K是一次相干处理间隔内的时域脉冲数；

S12，利用FFT域变换矩阵F_D获得阵元-多普勒域输出信号：

D_X_l＝T_X_l·F_D＝[S_{Dl_1}…S_{Dl_i}…S_{Dl_K}]_N×K

其中，S_{Dl_i}表示第l个距离单元第i个多普勒单元各阵元输出信号；

S13，利用凸优化算法对S_{Dl_i}进行稀疏重构，稀疏重构的约束方程如下：

其中，σ_i是第l个距离单元第i个多普勒单元接收信号的空域幅度分布，是函数g(x)取最小值时所对应的变量x的值，观测矩阵ψ_i是空域导向矢量构成的一组超完备基，大小是N×N_s，N_s是空域量化单元数，||·||₁表示L₁范数运算，ε_i是允许误差；

S14，对阵元-多普勒域数据逐多普勒单元进行空域稀疏重构，获得高分辨率的二维空时谱。

3.根据权利要求1所述的一种基于稀疏重构的KA-STAP杂噪协方差矩阵高精度估计方法，其特征在于，所述的步骤S2具体通过以下步骤实现：

S21，在二维空时谱上挑选出每个归一化空间频率上幅值最大的像素点并记录其坐标；

S22，将S21挑选出来的像素点中偏离杂波轨迹的噪声像素点剔除，得到筛选后的点集{(x_j,y_j)}，x_j是第j个像素点的横坐标，y_j是第j个像素点的纵坐标，j＝1,2,3...N_I，N_I是筛选出的符合要求的像素点的个数；

S23，对空时谱幅值向量A进行归一化处理，得到归一化后的空时谱幅值向量a：

其中， 是第j个像素点对应的真实幅度值，用归一化后的空时谱幅值向量中a_j作为第j个像素点对应的加权权值。

4.根据权利要求3所述的一种基于稀疏重构的KA-STAP杂噪协方差矩阵高精度估计方法，其特征在于，步骤S3利用加权最小二乘法对杂波轨迹进行拟合，对应的目标函数如下：

其中，表示求解使函数f(p₁,p₂)取得最小值时ψ的值， λ是机载雷达波长，f_r是脉冲重复频率，v是载机飞行速度，ψ是载机偏航角，是杂波散射体的俯仰角。

5.根据权利要求3所述的一种基于稀疏重构的KA-STAP杂噪协方差矩阵高精度估计方法，其特征在于，步骤S4中的先验杂波加噪声协方差矩阵满足下列公式：

其中，

δ_n ²是根据稀疏重构空时谱估计得到噪声功率，ψ_opt是根据加权最小二乘法解得ψ的最优解，I表示NK×NK维的单位阵。

6.根据权利要求1所述的一种基于稀疏重构的KA-STAP杂噪协方差矩阵高精度估计方法，其特征在于，所述的步骤S5具体通过以下步骤实现：

S51，利用NK×r维列满秩降维变换矩阵T对采样数据X_l进行降维处理，r是降维后的自适应处理自由度；

S52，计算降维后的目标导引矢量S_T、数据矢量X_lT和杂噪协方差矩阵R_0T：

S_T＝T^HS

X_lT＝T^HX_l

R_0T＝T^HR₀T

其中，S是没有降维的目标导引矢量，X_l是第l个距离单元的接收信号，l＝1,2,…,M，M是距离单元总数，R₀是先验杂波加噪声协方差矩阵；

S53，计算降维处理后的最优自适应权值w_T：

其中，μ_T是归一化复常数，

S54，根据最优自适应权值和回波数据矢量计算二维滤波输出Z_out：