CN110109066B - 一种新的迭代stap优化方法 - Google Patents

Abstract

本发明公开了一种新的迭代STAP优化方法，属于雷达技术领域，具体涉及一种期望约束下的迭代STAP优化技术。本发明首先对优化问题的约束进行分解并定义理想滤波器的期望，重新构造约束优化的目标函数后并进一步加入对权向量的约束，最后对构造的目标函数通过迭代算法进行迭代来求解滤波器权重。本发明与现有直接求STAP权向量的方法相比，本发明所带来的技术效果为：将大的全局问题分解为多个较小、较容易求解的局部子问题，并通过协调子问题的解而得到大的全局问题的解，能确保得到全局的最优解；同时，其杂波的抑制性能更好。

Description

一种新的迭代STAP优化方法

技术领域

本发明属于雷达技术领域，具体涉及一种期望约束下的迭代STAP(空时自适应处理)优化技术。

背景技术

空时自适应处理(STAP)是一种重要的杂波和干扰抑制方法，并广泛应用于雷达，声纳，地震探测和通信领域中。在这些领域，通常由于干扰或杂波在时空域中与目标信号耦合，所以需要进行空时自适应处理。

在准确估计待检测单元(cell under test,CUT)的杂波协方差矩阵后，通常通过拉格朗日乘子法去求解滤波器权向量，例如在文献《Y.Wu,T.Wang,J.Wu,and J.Duan,“Training sample selection for space-time adaptive processing inheterogeneous environments,”IEEE Geoscience and Remote Sensing Letters,vol.12,no.4,pp.691–695,2015.][H.Li,W.Bao,J.Hu,J.Xie,and R.Liu,“A trainingsamples selection method based on system identification for stap,”SignalProcessing,vol.142,pp.119–124,2018.][Li X,Feng D,Liu H W,et al.Dimension-Reduced Space-Time Adaptive Clutter Suppression Algorithm Based on Lower-RankApproximation to Weight Matrix in Airborne Radar[J].Aerospace&ElectronicSystems IEEE Transactions on,2014,50(1):53-69.》所公开的一种的求解滤波器权向量的方式。但在该方案中，一方面需要协方差矩阵的求逆运算；另一方面，当杂波协方差矩阵的秩远小于系统自由度时，STAP滤波器权重将是稀疏的。

文献《H.Li,W.Bao,J.Hu,J.Xie,and R.Liu,“A training samples selectionmethod based on system identification for stap,”Signal Processing,vol.142,pp.119–124,2018.][Li X,Feng D,Liu H W,et al.Dimension-Reduced Space-TimeAdaptive Clutter Suppression Algorithm Based on Lower-Rank Approximation toWeight Matrix in Airborne Radar[J].Aerospace&Electronic Systems IEEETransactions on,2014,50(1):53-69.][Fa R,Lamare R C D.Reduced-Rank STAPAlgorithms using Joint Iterative Optimization of Filters[J].IEEE Transactionson Aerospace and Electronic Systems,2011,47(3):1668-1684》公开了另一种求解滤波器权向量的方式，在该方案中，其未将STAP滤波器权重的稀疏考虑在内。

随着压缩感知技术的发展，产生了很多稀疏的STAP算法，如文献《Yang Z,DeLamare R C,Li X.-Regularized STAP Algorithms With a Generalized SidelobeCanceler Architecture for Airborne Radar[J].2012,60(2):674-686.][Jichao Z,Shengqi Z,Zhiqi G,et al.L 1-regularised joint iterative optimisation space-time adaptive processing algorithm[J].Iet Radar Sonar&Navigation,2016,10(3):435-441》公开的方案中，其通过在广义旁瓣相消器中加入了L1范数正则化，但是该方案中仍然存在求逆运算。从而导致运算复杂度高。

发明内容

本发明的发明目的在于：为了避免对杂波协方差矩阵的求逆，同时将STAP滤波器权重的稀疏性考虑在内，提出了一种期望约束下的迭代求解的STAP优化方法。本发明针对目前STAP滤波器权向量的求解和优化处理方案，即首先对优化SCNR进行等价变形，然后将滤波器权重的稀疏性考虑在内，之后对构造的目标函数通过迭代算法进行迭代来求解滤波器权重。

本发明的一种新的迭代STAP优化方法，具体包括下列步骤：

步骤1：设置最优权向量的优化模型为：

其中，w_l表示对雷达接收数据的第l个距离做空时自适应处理的滤波器权向量，R_x表示待检测单元的杂波协方差矩阵，s(w_s,w_t)表示目标的空时导向向量，w_t表示归一化的多普勒频率，w_s表示空间频率；符号()^H表示共轭转置；

即本发明将SCNR(信杂噪比)优化条件变形为公式(1)所示的形式。

用P(m)来表示待检测距离单元的第m个频率通道的能量，针对约束条件w_l ^Hs(w_s,w_t)＝1，可以转换为：

P(m)＝|w_l ^Hs(w_s,w_t)|²＝w_l ^Hs(w_s,w_t)s(w_s,w_t)^Hw_l (2)

令：

β_m＝s(w_s,w_t)s(w_s,w_t)^H (3)

则有：

P(m)＝ω_l ^Hβ_mω_l (4)

步骤2：设置滤波器期望约束。

本发明通过参数d(m)来描述滤波器杂波抑制性能得好坏：

在上式中，表示在对某个距离单元进行STAP处理时，m＝0表示待检测的频率通道，m≠0表示其他频率通道。

步骤3：设置目标函数。

本发明中，将STAP滤波器权向量的优化求解转换为：最小化ω_l ^HR_xω_l同时最小化d(m)和P(m)之间的平方误差。故将针对STAP滤波器权向量的优化求解的目标函数(代价函数)定义为：

其中，J(α,w_l)表示目标函数，ω_m和δ表示加权系数，可以根据具体要求进行调整，即ω_m表示针对第m个频率通道的权重，δ表示|w_l ^HR_lw_l|²的权重，R_l表示待检测单元的杂波协方差矩阵，α为待优化的缩放参数，M表示雷达接收阵列在一个相干处理周期内的脉冲数。

为了便于求解，将上式进一步转换成双凸优化问题：

首先定义参量r、A(m)和

则目标函数可以转换成双凸优化问题：

对权向量加入正则化约束：

w_l ^Hw_l＝1 (11)

该条件可以转化为

其中，

是方阵，除了E_i(i,i)为1其他都是0，sqrt()表示开平方，K是阵元数目。

所以，最优权向量的优化模型可以转化为：

步骤4：用迭代优化算法求解步骤3给出的最优权向量的优化模型，基于求解结果得到STAP处理中的滤波器权重。从而用于生成所需的雷达波形。

综上所述，由于采用了上述技术方案，与现有直接求STAP权向量进而生成所需雷达波形的方法相比，本发明的有益效果是：

(1)将大的全局问题分解为多个较小、较容易求解的局部子问题，并通过协调子问题的解而得到大的全局问题的解；能确保得到全局的最优解。

(2)所生成的雷达波形对杂波的抑制性能更好。

附图说明

图1是实施例中，第380距离单元雷达回波信号的频谱示意图；

图2是实施例中，本发明与现有方案的杂波抑制结果，其中图2-a为现有方案的杂波抑制结果，图2-b为本发明的杂波抑制结果；

图3是实施例中，本发明与现有方案的在目标所在方位角的杂波抑制结果，其中图3-a为对应图2-a中目标所在方位角的处理结果，图3-b为对应图2-b中目标所在方位角的处理结果；

图4是实施例中，不同输入SCNR条件下的杂波抑制性能。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合实施方式和附图，对本发明作进一步地详细描述。

为了避免对杂波协方差矩阵的求逆，同时将STAP滤波器权重的稀疏性考虑在内。本发明首先对对优化问题条件进行分解变形；然后定义滤波器的理想期望，然后重写目标函数并对权向量进行约束，最后通过ADMM(交替方向乘子法)算法去迭代求解。

本发明的工作原理如下：

对于含有N个阵元的线型阵列的雷达天线，其阵元间距用d表示；每K个脉冲作为一个相干处理周期。对于雷达接收的数据，对于第l个距离环的数据，将其表示为x_l，其中x_l∈C^NK。

用c_l，n_l分别表示杂波和噪声；a_l表示目标信号幅值；s(w_s,w_t)表示目标的空时导向向量，则s(w_s,w_t)的可以用下式进行表示：

其中，s_t(w_t)表示时域导向矢量，w_t为归一化的多普勒频率；s_s(w_s)表示空域导向矢量，w_s为空间频率；

表示Kronnecker乘积。s_t(w_t)、s_s(w_s)的具体形式如下：

其中，e表示自然底数，j表示叙述单位。

信号检测可以写做两类检测问题：

其中，H₀表示待检测信号存在有用信号；H₁表示待检测信号仅存在杂波和噪声。

本发明基于上述两类信号检测问题的表示方式实现STAP优化处理，其具体实现步骤如下：

步骤1：设置最优权向量的优化模型，即公式(1)所示的初始优化模型；

步骤2：设置滤波器期望约束，如公式(5)所示；

步骤3：设置目标函数J(r)，如公式(10)所示，进而将初始优化模型转换为：

步骤4：用迭代优化算法求解步骤3转换后的最优权向量的优化模型。

由于，上述优化问题中目标函数和约束都是r的四阶非凸多项式。为了解决这个问题，本发明通过引入一个辅助原始变量

将公式(13)所示的最优权向量的优化模型可以改写如下：

s.t.h-r＝0

定义：

G₁(h,r)＝h-r＝0 (20)

G₂(h,r)＝T(h)r-j＝T'(r)h-j＝0 (21)

T(h)＝[h^HE₂,h^HE₃,...,h^HE_L+1] (22)

T'(r)＝[E₂r,E₃r,...,E_L+1r] (23)

由于F(h,r)是双凸的，因此对于每个r它是h的二阶多项式，对于每个h它是r的二阶多项式。将带有原始变量r的非凸问题转换为具有两个原始变量r和h的双凸问题。

首先将双变量引入公式(18)从而产生增广拉格朗日形式，然后通过交替迭代原始变量和双变量来获得结果，则公式(18)的增广拉格朗日形式为：

其中

是双变量,ρ₁、ρ₂是惩罚参数。

在第(k+1)次迭代中的原始变量和双变量的更新为：

y^k+1:＝y^k+ρ₁G₁(h^k+1,r^k+1) (28)

z^k+1:＝z^k+ρ₂G₂(h^k+1,r^k+1) (29)

通过组合增广拉格朗日量中的线性和二次项并定义缩放的双变量u＝(1/ρ₁)y和v＝(1/ρ₂)z，重新构造公式(25)如下：

所以有：

u^k+1:＝u^k+G₁(h^k+1,r^k+1) (33)

v^k+1:＝v^k+G₂(h^k+1,r^k+1) (34)

其中，h^k、r^k、u^k和v^k分别为第k次迭代时h,r,u,v的表达式(即计算结果)；h^k+1、r^k+1、u^k+1和v^k+1分别为第k+1次迭代时h,r,u,v的表达式。

引理1：更新h和r是凸问题，因此它们可以完全执行：

h^k+1＝Ψ^-1χ (35)

其中

I表示单位矩阵，其维度为(M+1)×(M+1)。

所以有：

r^k+1＝Φ^-1γ (38)

通常情况下，两种参数原始残差和双残差，可以用作收敛的标准。

在上式中，

是原始残差，_s ^k是双残差。当算法收敛时，原始残差和双残差都应该很小。因此，收敛条件可以表述如下：

在公式(41)中，

和ε^dual分别是上述三个残差(

和_s ^k)的合理误差，可以用公式(42)计算它们的值：

其中ε^abs＞0是绝对误差，ε^rel＞0是相对误差。

综上，步骤4中，迭代优化算法求解步骤3转换后的最优权向量的优化模型的具体处理过程可以描述为：

初始化变量h、r、u和v值，分别记为h⁰,r⁰,u⁰,v⁰；以及初始化惩罚参数ρ₁和ρ₂，绝对误差ε^abs和相对误差ε^rel的值；

在第k(本具体实施方式中，k的初始值为0)次迭代处理中，

步骤Step 1：判断是否满足公式(42)所示的收敛条件，若是，则终止迭代处理，得到求解结果；否则继续执行Step 2；

步骤Step 2：更新变量h^k+1、r^k+1、u^k+1和v^k+1；

即根据公式(35)更新变量h^k+1、根据公式(38)更新变量r^k+1；

根据公式(33)更新变量u^k+1，根据公式(34)更新变量v^k+1；

再更新k＝k+1，继续进行迭代处理，即继续执行步骤Step 1～Step 2。

实施例

为了验证所提方案，采用现有方案和所提方案对海杂波实测数据进行杂波抑制处理。

本实施例中，用到的实测数据的采集参数是：雷达工作频率f₀＝18.3MHZ，每512个脉冲作为一个相干处理周期，脉冲间隔T＝12ms。在第380个距离单元内，有一个目标处于多普勒频率1.139Hz处。其频谱如图1所示。

本发明先对优化SCNR进行等价变形，然后将滤波器权重的稀疏性考虑在内，之后对构造的目标函数通过迭代算法进行迭代来求解滤波器权重。

现有方案(参考文献《Yang Z,De Lamare R C,Li X.-Regularized STAPAlgorithms With a Generalized Sidelobe Canceler Architecture for AirborneRadar[J].2012,60(2):674-686.][Jichao Z,Shengqi Z,Zhiqi G,et al.L 1-regularised joint iterative optimisation space-time adaptive processingalgorithm[J].Iet Radar Sonar&Navigation,2016,10(3):435-441》、本发明所提方法对第380距离单元进行杂波抑制。

图2-a中的残留杂波依然很强，图2-b相对于图2-a中的残留杂波较弱，表明本文所提算法求得的滤波器权向量进行滤波后的杂波抑制性能相对于现有方案较好。下面分别从图2中截取目标所在方位角即0.3704rad的数据处理结果，如图3-a、图3-b所示。

图3-a中最大残留杂波的值为-6.421dB，图3-b最大残留杂波的值为-8.991dB，本文所提算法求得的滤波器权向量进行滤波后的杂波抑制性能相对于现有方案提高了2.57dB。

为了比较本文所提方法与上述文献中的方法在不同输入SCNR条件下的杂波抑制性能，本次仿真实验在输入SCNR的取值从-40dB到-15dB取值时，比较三种方法在不同输入SCNR条件下的杂波抑制性能，结果如图4。纵坐标为目标频率处输出SCNR与其余频率范围内最大输出SCNR的差值。该差值表征了目标多普勒频率处的输出SCNR的凸起程度，差值越大，说明杂波抑制效果越好，目标更容易被检测到。经计算，与现有方案文献《Yang Z,DeLamare R C,Li X.-Regularized STAP Algorithms With a Generalized SidelobeCanceler Architecture for Airborne Radar[J].2012,60(2):674-686.][Jichao Z,Shengqi Z,Zhiqi G,et al.L 1-regularised joint iterative optimisation space-time adaptive processing algorithm[J].Iet Radar Sonar&Navigation,2016,10(3):435-441》相比，所提方法的输出SCNR提高了2dB左右；

通过以上仿真，验证了本发明所公开的STAP优化方法的性能优势：杂波的抑制性能更好。

以上所述，仅为本发明的具体实施方式，本说明书中所公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换；所公开的所有特征、或所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以任何方式组合。

Claims (1)

1.一种新的迭代STAP优化方法，其特征在于，包括下列步骤：

设置最优权向量的优化模型为：

其中，M表示一个相干处理周期内的脉冲数，ω_m和δ表示加权系数；

参量r、A(m)和

分别为：

其中，α为待优化的缩放参数，w_l表示空时自适应滤波器权向量，R_l表示杂波协方差矩阵；

参数d(m)表示第m个频率通道的滤波器杂波抑制性能，其取值为：

表示待检测的频率通道，m≠0表示非待检测的频率通道；

参量β_m＝s(w_s,w_t)s(w_s,w_t)^H，其中，s(w_s,w_t)表示目标的空时导向向量，w_t表示归一化的多普勒频率，w_s表示空间频率；

矩阵E_i+1为除了元素E_i+1(i+1,i+1)为1，其余元素均为0的(MK+1)×(MK+1)的方阵，K表示雷达天线的阵元数；

sqrt( )表示开平方；

采用迭代优化方式求解最优权向量的优化模型，基于求解结果得到STAP处理中的滤波器权重；

所述迭代优化算法求解最优权向量的优化模型的具体过程为：

初始化变量h、r、u和v值，分别记为h⁰,r⁰,u⁰,v⁰；以及初始化惩罚参数ρ₁和ρ₂，绝对误差ε^abs和相对误差ε^rel的值；

在第k次迭代处理中，执行下列步骤，其中k的初始值为0；

步骤Step 1：判断是否满足迭代收敛条件，若是，则终止迭代处理，得到求解结果；否则继续执行Step 2；

所述迭代收敛条件为：双残差s^k和两个原始残差

是否均小于或等于各自的误差ε^dual、

和

其中，双残差s^k、原始残差

和

具体为：

其中，T(h^k)＝[(h^k)^HE₂,(h^k)^HE₃,...,(h^k)^HE_L+1]，

原始残差

的误差

为：

原始残差

的误差

为：

双残差s^k的误差为ε^dual为：

h^k、r^k、u^k和v^k分别为变量h,r,u,v在第k次迭代时的值；

步骤Step 2：更新变量h、r、u和v值在第k+1次迭代时的值h^k+1、r^k+1、u^k+1和v^k+1，再更新k＝k+1，继续进行迭代处理；

其中，h^k+1、r^k+1、u^k+1和v^k+1具体为：

h^k+1＝Ψ^-1χ；

r^k+1＝Φ^-1γ；

u^k+1＝u^k+(h^k+1-r^k+1)；

v^k+1＝v^k+(T(h^k+1)r^k+1-j)；

其中，

I表示单位矩阵；

T'(r^k)＝[E₂r^k,E₃r^k,...,E_L+1r^k]；

Images