基于多普勒预滤波的降维STAP微下击瀑流中心风速估计方法与流程

本发明属于雷达信号处理技术领域，特别是涉及一种基于多普勒预滤波的降维STAP微下击瀑流中心风速估计方法。

背景技术：

下击瀑流是能在地面产生17.9m/s以上辐散风的一种强烈的下沉气流。强烈的空气对流或者较为闭塞的地理环境都有可能诱发下击瀑流。随着城市的不断发展，终端机场附近的地形环境愈发复杂，更易诱发微下击瀑流而威胁飞机起降安全。其中平辐散尺度小于4km，持续时间不超过10min的微下击瀑流对民用航空飞行安全威胁尤为巨大。当飞机在起降阶段进入微下击瀑流作用区域时，飞行员往往没有足够的时间和空间来调整飞机姿态以应对升力的突然增加或丢失对飞行高度、飞行速度以及飞机受力载荷的影响，因此飞机很容易迅速失速进而有可能产生坠毁。

机载气象雷达可以对雷雨、风切变、湍流等气象情况进行探测与预警，是飞机实时检测航路气象信息的重要设备，但机载气象雷达在前下视检测低空风切变时，有用信号会淹没在强杂波背景中。抑制地杂波的基本思想是选用合适的滤波器在抑制地杂波的同时保留风切变信号。传统方法包括杂波图法、基于参数化模型的谱估计方法、利用模式分析的扩展Prony方法、零陷滤波器法等，但在强杂波条件下，这些方法难以完全消除地杂波，因此残余杂波仍会影响风速估计结果的准确性。

相较于传统单天线雷达，相控阵雷达由于在接收回波中增加了目标信号的空域信息，使其在强杂波背景下的杂波抑制和信号检测更具优势，能够更好地实现对目标的精确检测，目前已有机构开始展开对新一代机载相控阵气象雷达的研究。在相控阵雷达中应用的空时自适应处理(Space-Time Adaptive Processing，STAP)技术可利用空时二维信息抑制地杂波，有效提高了相控阵雷达的目标检测能力。但是，传统STAP技术主要针对点目标的检测与估计，无法直接运用在风切变这样的分布式目标检测上。并且机载气象雷达工作在前视状态，前视阵的杂波谱在距离上不平稳，从而导致了传统的空时最优处理器不能直接应用于机载气象雷达中。同时由于构造空时最优处理器所需的协方差矩阵维数很高，对其进行估计和求逆运算量巨大，因此实现实时处理较为困难。综上，上述原因制约了STAP技术在机载气象雷达中的应用与推广。

技术实现要素：

为了解决上述问题，本发明的目的在于提供一种能够保证参数估计精度，同时降低运算复杂度的基于多普勒预滤波的降维STAP微下击瀑流中心风速估计方法。

为了达到上述目的，本发明提供的基于多普勒预滤波的降维STAP微下击瀑流中心风速估计方法包括按顺序进行的下列步骤：

1)利用杂波谱配准法逐距离单元处理雷达接收数据，求得待检测距离单元的重构杂波协方差矩阵估计所需的IID(独立同分布)样本；

2)利用多普勒滤波器选取某一多普勒频率作为中心多普勒频率，构造该频率下的时域降维矩阵，用于处理上述IID样本中的数据；

3)利用时域降维矩阵中的数据估计上述选取的中心多普勒频率下的空域二次协方差矩阵，求取降维的最优权矢量；

4)利用上述降维的最优权矢量处理降维数据，获得目标区域内风场的中心风速估计值及中心风速分布情况。

在步骤1)中，利用杂波谱配准法逐距离单元处理雷达接收数据，求得待检测距离单元的重构杂波协方差矩阵估计所需的IID样本的方法是：首先对雷达接收到的数据逐距离单元进行杂波配准，以消除不同距离单元的杂波由于距离依赖导致的非均匀性，然后将待检测距离单元相邻的距离单元作为参考单元，求得估计待检测距离单元的重构杂波协方差矩阵所需的IID样本。

在步骤2)中，利用多普勒滤波器选取某一多普勒频率作为中心多普勒频率，构造该频率下的时域降维矩阵，用于处理IID样本中的数据的方法是：首先，将杂波配准后的IID样本逐距离单元传给滤波器，其次，用设计好的多普勒滤波器对各个距离单元的雷达接收数据进行过滤处理，使得过滤后的信号频率分布在中心频率附近，最后，构造各个多普勒中心频率下的时域降维矩阵，用来对多普勒滤波后的数据进行时域降维。

在步骤3)中，利用时域降维矩阵中的数据估计上述选取的中心多普勒频率下的空域二次协方差矩阵，求取降维的最优权矢量的方法是：首先，取相邻的多普勒通道作为辅助通道，并取辅助通道的降维数据共同组成一个新的矢量，其次，用新的矢量构造该多普勒通道的空域二次协方差矩阵，最后，根据线性约束最小方差准则求解数学优化问题，得到匹配该多普勒通道的空域二次协方差矩阵的滤波器最优权矢量。

在步骤4)中，利用上述降维的最优权矢量处理降维数据，获得目标区域内风场的中 心风速估计值及中心风速分布情况的方法是：首先，利用降维的最优权矢量自适应抑制对应的多普勒通道空域上的有源定向窄带干扰，即通过滤波器的杂波，其次，积累主瓣方向上的目标信号，并利用频率质心法求出该距离单元内风场的中心风速估计值，最后，循环处理全部距离单元所有多普勒通道的接收数据，得到目标区域内风场的中心风速分布情况。

本发明提供的基于多普勒预滤波的降维STAP微下击瀑流中心风速估计方法是针对相控阵体制的机载气象雷达，利用微下击瀑流分布式气象目标特点，利用多通道联合自适应处理方法构造自适应处理器，估计微下击瀑流中心风速。本发明方法可以在低信噪比、强杂噪比条件下有效地估计出微下击瀑流中心风速，仿真实验验证了本方法的有效性，且由于采用了降维结构，较最优处理器在运算量上有了明显减少。

附图说明

图1为本发明提供的基于多普勒预滤波的降维STAP微下击瀑流中心风速估计方法流程图。

图2为机载前视阵雷达阵模型图。

图3为雷达接收信号的最小方差谱。

图4为估计所得的杂波的最小方差谱。

图5为第80号距离单元改善因子图。

图6为滤波器频响特性图。

图7为特定距离单元滤波器频响特性图。

图8为风速估计结果图.

具体实施方法

下面结合附图和具体实例对本发明提供的基于多普勒预滤波的降维STAP微下击瀑流中心风速估计方法进行详细说明。

如图1所示，本发明提供的基于多普勒预滤波的降维STAP微下击瀑流中心风速估计方法包括按顺序进行的下列步骤：

1)利用杂波谱配准法逐距离单元处理雷达接收数据，求得待检测距离单元的重构杂波协方差矩阵估计所需的IID样本；

机载前视阵雷达阵模型如图2所示，设载机速度为V_R，飞行高度为H，机载气象雷达(以下简称雷达)天线系统由N阵元均匀线阵(也可是由面阵经过微波合成的等效线 阵结构)组成，天线阵元间距d＝0.5λ，其中λ为雷达发射脉冲的波长。前视阵指天线阵面与载机速度方向夹角为90°。设雷达系统脉冲重复频率(Pulse Repetition Frequency，PRF)为f_r，相干处理脉冲数为K。图中，θ为方位角，为俯仰角，ψ为空间锥角，且满足

在本发明中，x_l表示第l个距离单元的NK×1维空时快拍数据，其表达式如下：

x_l＝s_l+c_l+n_l (1)

其中，s_l为第l个距离单元内微下击瀑流场产生的雷达回波信号，c_l为地杂波，在此假设杂波无起伏无模糊，n_l为加性高斯白噪声。

对于第l个距离单元内的风场，雷达对其的采样数据可以写成一个N×K的矩阵S_l。其中，矩阵S_l的第n行、第k列元素表示雷达第n(n＝1,2,…N)个阵元、第k(k＝1,2,…K)个脉冲对风场回波的采样，当该距离单元内雷达波束照射范围内共有Q个气象散射粒子时，其具体表达式如下：

式中

和分别表示第q(q＝1,2,…,Q)个气象散射粒子的空间角频率和时间角频率，θ_q、分别表示该气象散射粒子相对于雷达的方位角和俯仰角，A_q为第q个气象散射粒子的一次散射幅度，R_q为第q个气象散射粒子与载机的斜距。将矩阵S_l展开成为NK×1维列向量，即为微下击瀑流场空时快拍s_l。

假设第i个距离单元为待检测距离单元，接收信号为x(i)，第j个距离单元为参考距离单元，接收数据为x(j)。存在补偿矩阵T_j,i使得处理后的参考数据y(j)与接收信号x(i)服从同样的分布特性，其中y(j)＝T_j,ix(j)。

补偿矩阵T_j,i可由待检测距离单元与参考距离单元间杂波协方差矩阵的逻辑关系得到。待检测距离单元的重构杂波协方差矩阵为：

式中{θ_m}|_m＝1,2,…M由水平方位角离散化得到，a_m为对应方位内的回波幅度，表示空时导向矢量，且有：

上式中j为虚数单位，和a_t(v)分别表示空间与时间导向矢量，为待检测距离单元相对雷达的俯仰角。同样写出参考距离单元的重构杂波协方差矩阵为R_j，并且待检测距离单元的重构杂波协方差矩阵R_i与参考距离单元的重构杂波协方差矩阵R_j满足如下关系：

分别对对称的待检测距离单元的重构杂波协方差矩阵R_i和参考距离单元的重构杂波协方差矩阵R_j进行特征分解，可以得到：

其中U_i、U_j为正交矩阵，G_i、G_j为实对角矩阵，将式(7)带入式(6)中得到

处理后的参考数据y(j)与待检测距离单元内杂波的空时二维分布一致。更新参考距离单元，同样利用上述杂波谱配准法进行处理，可以得到与待检测距离单元杂波分布特性相同的多个IID样本。

2)利用多普勒滤波器选取某一多普勒频率作为中心多普勒频率，构造该频率下的时域降维矩阵，用于处理上述IID样本中的数据；

在本发明中，将雷达主瓣宽度作为先验信息，建立微下击瀑流场的空间导向矢量。

当雷达主瓣方向俯仰角为中心方位角为θ_i时，设其照射范围内微下击瀑流场的广义空间导向矢量为其表达式如下：

其中表示中心方位角为θ_i，中心俯仰角为方位上，点目标的空间导向矢量；为确定性角信号密度函数，本发明中表示微下击瀑流场在中心方位角θ_i和中心俯仰角上的扩展，通常是为以分布源中心DOA为对称中心的单峰对称函数。由于雷达主瓣范围照射内气象散射粒子的数量很大，并且没有一个处于主导地位，则根据中心极限定理，微下击瀑流场在中心方位角θ_i和中心俯仰角上的扩展可以表示为：

其中，σ_θ表示θ_i方向上的角度扩展，表示方向上的角度扩展。

本发明以三通道联合自适应处理器(3DT)为例进行说明。三通道联合自适应处理方法是通过前置加权多普勒滤波将全空时分布的杂波局域化为某个多普勒通道中心多普勒频率下的定向有源干扰，由此可实现将目标信号由空时二维降至空域一维。设第m个多普勒通道的中心多普勒频率为f_dm，则降维后的第m个多普勒通道的空域接收数据Y_m可由式(11)求得：

其中，I_N表示空域单位矩阵，x_l表示接收的第l个待检测距离单元的信号，W_tm表示第m个多普勒通道的多普勒权值，组成了实际的时域降维矩阵。本发明中：

3)利用时域降维矩阵中的数据估计上述选取的中心多普勒频率下的空域二次协方差矩阵，求取降维的最优权矢量；

在空域自适应处理中，取第(m1)和第(m+1)个多普勒通道作为辅助通道，其阵列数据矢量分别用Y_m+1和Y_m-1表示。定义一新的矢量：

利用新的矢量Z_m构造第m个多普勒通道的空域二次协方差矩阵R_Zm：

则根据LCMV准则可知，能够同时抑制杂波和实现信号匹配的降维处理器可以描述为如下数学优化问题：

其中，S_Zm为分布式目标的二维空时导向矢量，可得：

g_m+i(i＝±1)为常数，表示选定的辅助通道与主通道的归一化增益：

则可解得降维的最优权矢量w_Zm为：

4)利用上述降维的最优权矢量处理降维数据，获得目标区域内风场的中心风速估计值及中心风速分布情况。

式(18)示出的降维的最优权矢量可以自适应抑制第m个多普勒通道空域上的有源定向窄带干扰，并积累主瓣方向上的目标信号。干扰抑制和信号匹配的结果为：

更新多普勒主通道，求得不同中心多普勒频率下的最优权值，并分别求取此时的滤波器输出功率，进而求得该距离单元内风场的中心风速估计值为：

依上述步骤逐距离单元对雷达接收的回波数据进行处理，可得风场速度随距离变化的分布情况，进而获得目标区域范围内风场的中心风速分布情况。

本发明提供的基于多普勒预滤波的降维STAP微下击瀑流中心风速估计方法的效果可以通过以下仿真结果进一步说明。

仿真参数设置：微下击瀑流场分布于飞机前方8.5-16.5km处。天线阵为阵元数N＝8 的前视理想均匀线阵，阵元间距为d＝λ/2，主瓣波束水平方位角为60°，俯仰角为0°，波束宽度为3.5°，雷达波长为0.05m，脉冲重复频率为7000Hz，雷达最小分辨距离为150m，相干处理脉冲数K＝64，杂噪比为40dB，信噪比为5dB；载机速度为75m/s，飞行高度为600m，归一化多普勒谱宽σ_f＝0.05。

仿真所得的地杂波的最小方差谱如图3(a)所示，可以看出，机载前视阵雷达地杂波最小方差谱呈椭圆形分布，频率扩展比较严重。实际雷达接收信号为微下击瀑流信号、地杂波以及噪声的叠加，如图3(b)所示。由于雷达回波中微下击瀑流回波信号功率要远小于地杂波功率，地杂波在很大程度上淹没了微下击瀑流信号，造成了对微下击瀑流场检测的困难。

对选取的杂波训练样本利用杂波谱配准法消除距离依赖性后估计的第80号和第10号距离单元的杂波协方差矩阵的最小方差谱如图4所示，对比图4和图3(a)可知地杂波的多普勒展宽明显变窄，不同距离单元的杂波所呈现的“圆弧”的半径各不相同，由此可以看出雷达回波的非均匀性得到了较好地消除。

图5以80号距离单元为例，对比了可用于微下击瀑流风速估计的3DT处理器与STAP最优处理器的改善因子，可以看出，本发明方法在性能上可以较好地逼近STAP最优处理器的上界。

图6和图7为第100号距离单元3DT处理器的频响特性。其中图6为全多普勒域的频响特性图，图7显示了中心多普勒频率分别为-20、0、20、40、60、80Hz的六个多普勒通道的处理器空域响应，由频响特性图可以看出，各多普勒通道下的空域滤波器均在目标处形成了明显的增益。

本发明方法与最优处理器方法、传统地杂波抑制及风速估计方法的风速估计结果对比如图8所示。可以看出在8.5～16.5km范围内，风速随距离呈反“S”变化；在相同杂噪比、信噪比下，本发明方法的风速估计结果与最优处理器结果性能相当，优于传统方法。

表1分析了最优处理器和本发明方法的运算量，其中仅考虑了对消单一距离门接收数据，即N_t＝1。以文本仿真进行计算，最优处理器的运算量为4.19×108，本发明方法的运算量仅为2.6×106，由此可见，实时性得到了极大的提升。

表1