基于空时最优处理器的低空风切变风速估计方法

摘要

本发明涉及一种基于空时最优处理器的低空风切变风速估计方法。本发明属于机载气象雷达信号处理领域。一种基于空时最优处理器的低空风切变风速估计方法，包括以下步骤：1)从关联的机载设备中读取飞机状态参数、雷达天线扫描参数等作为先验信息；2)选取待检测距离单元，利用空时插值法校正地杂波的距离依赖性；3)利用先验信息补偿飞机速度；4)估计待检测距离单元的地杂波协方差矩阵；5)利用空时最优处理器抑制地杂波并进行风速估计；6)重复步骤2)-5)，依次对范围内所有距离单元的回波数据进行风速估计。本发明针对强杂波背景下风速估计精度变差的问题，在机载气象雷达中引入空时信息，提出了一种基于空时最优处理器的低空风切变风速估计方法。仿真结果证明了该方法的有效性。

说明书

技术领域

本发明属于机载气象雷达低空风切变检测技术领域，特别是涉及一种基于空时最优处 理器的低空风切变风速估计方法。

背景技术

机载气象雷达可对雷雨、风切变、湍流等灾害天气进行探测与预警，是飞机实时感知 航路气象的重要设备，其对危险天气的探测能力对保证飞机行驶的安全性至关重要。低空 风切变常指高度600米以下风向、风速突然变化的气象现象，通常是强气流由空中冲击到 地面后向四周扩散的过程。低空风切变具有发生突然、时间短、尺度小、强度大和不易检 测等特点，当飞机起飞和着陆进入强低空风切变区域时，会面临由风向、风速变化带来的 升力损失，若处理不得当或缺乏足够的调节空间就可能遭遇飞行事故。因此，低空风切变 检测技术的研究已是航空运输领域的一项重要课题。

机载气象雷达下视检测低空风切变时，有用信号会被强地杂波所覆盖。在低空风切变 风速估计的工作流程中，地杂波抑制一直是最大的难点。传统的在地杂波背景下估计风场 风速的方法包括杂波图法、基于参数化模型的谱估计方法、利用模式分析的扩展Prony方 法、零陷滤波器法等，韩雁飞等提出了基于压缩感知的低空风切变风速估计方法，能够在 脉冲数较少时实现风速的精确估计，但该方法在应用时需要确定杂波谱宽，且并未涉及风 速估计过程中的强杂波抑制问题。这些方法的本质在于寻找合适的凹口在抑制地杂波的同 时保留风切变信号，但在强杂波背景下，地杂波难以完全消除，残余杂波仍会覆盖掉风切 变信号的多普勒信息，影响风速估计结果的准确性。

相对于传统机载气象雷达的单天线体制，相控阵天线体制具有灵活度高、扫描速度快、 易于波束赋形的优点，其回波信号中包含目标空间采样信息。目前国外已经对机载相控阵 气象雷达开展了相关的研究工作。2013年9月，美国国家大气研究中心(National Center for  Atmospheric Research,NCAR)表示，美国正在研制的下一代机载气象雷达采用双极化相 控阵体制。空时自适应处理(Space Time Adaptive Processing,STAP)是机载相控阵雷达杂 波抑制与目标检测的关键技术，空时最优处理器可在对有用信号输出功率不变的条件下， 尽可能地抑制杂波、干扰及噪声。在我国，STAP技术在机载预警雷达上的应用已经颇为 成熟，但此种应用多针对点目标，不能直接用于低空风切变等分布式目标的检测与估计， 且机载气象雷达工作在前视状态，前视阵的杂波谱在距离上不平稳，因此传统的空时最优 处理器不能直接应用于机载气象雷达中。

发明内容

为了解决上述问题，本发明的目的在于提供一种在高杂噪比、低信噪比时，仍然可以 精确估计风场风速的基于空时最优处理器的低空风切变风速估计方法。

为了达到上述目的，本发明提供的基于空时最优处理器的低空风切变风速估计方法， 包括按顺序进行的下列步骤：

1)从关联的机载设备中读取飞机状态参数、雷达天线扫描参数等作为先验信息；

2)选取待检测距离单元，利用空时插值法校正地杂波的距离依赖性；

3)利用先验信息补偿飞机速度；

4)估计待检测距离单元的地杂波协方差矩阵；

5)利用空时最优处理器抑制地杂波并进行风速估计；

6)重复步骤2)-5)，依次对范围内所有距离单元的回波数据进行风速估计。

在步骤1)中，所述的从关联的机载设备中读取飞机状态参数、雷达天线扫描参数等 作为先验信息是读取飞机速度、飞行高度、天线扫描的主瓣方位角和俯仰角等作为先验信 息。

在步骤2)中，所述的选取待检测距离单元，利用空时插值法校正地杂波的距离依赖 性是利用空时插值法，构造某一距离单元与待检测距离单元的插值变换矩阵T，使用该距 离单元数据与插值变换矩阵T相乘，使该距离单元内地杂波的空时二维分布与待检测距离 单元趋于一致。依次处理所有距离单元数据，可得到与待检测距离单元地杂波分布特性相 同的多个独立同分布(Independent and Identically Distributed,IID)样本。

在步骤3)中，所述的利用先验信息补偿飞机速度是利用上述先验信息对接收信号进 行相位补偿，剔除机载平台运动带来的多普勒偏移，以消除机载平台的运动性对参数估计 的影响。

在步骤4)中，所述的估计待检测距离单元的地杂波协方差矩阵是利用步骤2)中得 到的多个独立同分布样本求取其协方差矩阵的统计平均值。

在步骤5)中，所述的利用空时最优处理器抑制地杂波并进行风速估计是指，构造可 自适应抑制机载前视阵地杂波并积累待检测单元内风切变信号的最优处理器；根据低空风 切变的分布特性，设定多普勒谱宽与中心频率，构造风切变信号的协方差矩阵，并联合步 骤4)中计算得到的待检测距离单元内地杂波协方差矩阵求解最优处理器的权矢量；按多 普勒中心频率做循环，更新最优处理器权矢量，当输出信号功率最大时，此时设定的多普 勒中心频率即为待检测距离单元内风切变信号的多普勒频率估计值，进而能够计算待检测 距离单元内的风场速度。

在步骤6)中，所述的重复步骤2)-5)，依次对范围内所有距离单元的回波数据进行 风速估计是指更新待检测距离单元，再次利用步骤2)-5)估计其风切变信号多普勒中心 频率。将范围内所有距离单元处理完毕后，得到低空风切变场风速随距离的变化曲线。

本发明提供的基于空时最优处理器的低空风切变风速估计方法是将相控阵体系引入 到机载气象雷达中，利用空时插值法校正前视阵地杂波的距离依赖性，利用最优处理器自 适应地抑制地杂波并匹配低空风切变信号，估计风切变信号的多普勒频率后得到风场速度 估计。在强杂波条件下，本发明方法仍能有效地抑制地杂波并精确地估计风场风速。

附图说明

图1为本发明提供的基于空时最优处理器的低空风切变风速估计方法的工作原理图；

图2为机载前视阵回波仿真示意图；

图3为本发明提供的基于空时最优处理器的低空风切变风速估计方法的基本流程图

图4为补偿飞机速度后的仿真信号空时二维谱；

图5为空时插值法补偿地杂波距离依赖性效果图；

图6为空时处理器频响特性；

图7为本发明方法多普勒估计结果；

图8为本发明方法风场速度估计结果。

具体实施方式

下面结合附图和具体实施例对本发明提供的基于空时最优处理器的低空风切变风速 估计方法进行详细说明。

如图1所示，本发明提供的基于空时最优处理器的低空风切变风速估计方法包括按顺 序进行的下列步骤：

1)从关联的机载设备中读取飞机状态参数、雷达天线扫描参数等作为先验信息：

从机载导航设备中读取飞机速度、飞行高度、天线扫描的主瓣方位角和俯仰角等作为 信号处理的先验信息。

2)选取待检测距离单元，利用空时插值法校正地杂波的距离依赖性：

本发明假设天线体制为均匀线阵。如图2所示为机载前视阵回波仿真示意图，设定载 机速度为V；脉冲重复频率为f_r；阵元间隔为d；天线主瓣波束水平方位角为θ₀，俯仰角 为杂波单元水平方位角和俯仰角表示为θ和杂波单元空间角频率和时间角频率分 别表示为和且有

雷达回波信号x＝s+c+n，其中x为低空风切变信号，c为地杂波，n为高斯白噪声， 矩阵形式可写为

x＝[x₁  x₂  …  x_L]^T                 (2)

其中$x_l={\left(\begin{array}{cccccccccccccccc}{x}_{\mathrm{1,1}}^{\left(l\right)}& x_{\mathrm{1,2}}^{\left(l\right)}& ...& x_{1,N}^{\left(l\right)}& ,& x_{\mathrm{2,1}}^{\left(l\right)}& x_{\mathrm{2,2}}^{\left(l\right)}& ...& x_{2,N}^{\left(l\right)}& ,& ...& ,& x_{K,1}^{\left(l\right)}& x_{K,2}^{\left(l\right)}& ...& x_{K,N}^{\left(l\right)}\end{array}\right)}^T,$表 示将每个距离单元的接收数据排成NK×1的列向量N为阵元数，K为脉冲数，L为最大距 离单元数。

空时插值法通过一个插值变换矩阵使得训练距离单元的地杂波空时分布与待检测距 离单元一致，从而达到校正地杂波距离相关性的目的，具体实现步骤如下：

a)选取第i号距离单元为待检测距离单元，构造待检测距离单元的空时导向矢量矩阵 V_i。首先将水平方位角均匀分成M份，得到对应{θ_m}|_{m＝1,2,…,M}。然后构造每个水平方位角 θ_m和距离单元俯仰角对应的空时导向矢量其中

$\left(\begin{array}{c}{v}_t\left({\theta }_m\right)={[1,e^{j{\omega }_t\left({\theta }_m\right)},...,e^{j(K-1){\omega }_t\left({\theta }_m\right)}]}^T\\ v_s\left({\theta }_m\right)={[1,e^{j{\omega }_s\left({\theta }_m\right)},...,e^{j(N-1){\omega }_s\left({\theta }_m\right)}]}^T\end{array}\right)---\left(3\right)$

式中j表示虚数单位，可以得到待检测距离单元的空时导向矢量矩阵为

V_i＝[v(θ₁)  v(θ₂)  …  v(θ_M)]             (4)

b)选取第j号距离单元为参考距离单元，构造该参考距离单元的空时导向矢量矩阵 V_j，可得到第j号与第i号距离单元的插值变换矩阵为T_j,i＝V_i(V_j)⁺，其中(·)⁺表示求伪逆 运算。处理后的参考距离单元数据为y_j＝T_j,ix_j，其与待检测距离单元内地杂波的空时二 维分布一致。

c)更新参考距离单元，同样利用空时插值法进行处理，可以得到地杂波分布特性与 待检测距离单元相同的多个独立同分布样本。

3)利用先验信息补偿飞机速度：

利用上述先验信息对接收信号进行相位补偿，剔除机载平台运动带来的多普勒频移， 以消除机载平台的运动性对参数估计的影响。

4)估计待检测距离单元的地杂波协方差矩阵：

插值变换处理后，估计待检测距离单元地杂波数据的统计协方差矩阵为

${\hat{R}}_i=\frac{1}{L}(\underset{j=1,j\ne i}{\overset{L}{\Sigma }}{y}_j{y}_j^H+x_i{x}_i^H)---\left(5\right)$

5)利用空时最优处理器抑制地杂波并进行风速估计：

可在抑制地杂波的同时匹配低空风切变信号的空时最优处理器可以描述为如下的数 学优化问题：

$\left(\begin{array}{c}\min {\omega }^H{R}_{c+n}\omega \\ \mathrm{st}.{\omega }^H{R}_s\omega =1\end{array}\right)---\left(6\right)$

其中ω为最优处理器的权矢量，R_c+n为杂波噪声协方差矩阵，R_s为低空风切变信号 协方差矩阵。

构造代价函数

H(ω,λ)＝ω^HR_c+nω+λ(1-ω^HR_sω)             (7)

其中，λ为拉格朗日算子，当时，式(6)可等效为求解下式广义特征值问题：

R_c+nω＝λR_sω                   (8)

此时，拉格朗日算子λ可以看做是相应的广义特征值。由于矩阵R_c+n和R_s是半正定 的，所以上式得到的所有广义特征值都是非负实数。当式(8)成立时，代价函数取得最小值 λ。将左乘到式(8)两边，得到

$R_{c+n}^{-1}{R}_s\omega =\frac{1}{\lambda }\omega ---\left(9\right)$

则有其中p{·}表示求解矩阵最大特征值对应的特征向量。由于地杂波 功率远高于噪声功率，可近似认为得到

$\omega =p\left\{{\hat{R}}_i^{-1}{R}_s\right\}---\left(10\right)$

由上式可知，求解最优处理器权矢量ω需要杂波协方差矩阵和信号协方差矩阵R_s两个关键信息，其中杂波协方差矩阵可由式(5)得到，而在实际检测过程中，风切变信 号协方差矩阵R_s是未知的，上式仅代表对目标R_s进行积累。

风切变目标R_s的参数化模型为

${\hat{R}}_s(f_0,f_s,{\sigma }_f)=\mathrm{PA}(f_0,f_s)B\left({\sigma }_f\right)A{(f_0,f_s)}^H---\left(11\right)$

其中，P为回波信号功率，f_s为信号空间频率，f₀为信号多普勒中心频率，σ_f为信号多 普勒谱宽，且有

$A=\mathrm{diag}(A_t\otimes A_s)$

$\left(\begin{array}{c}{A}_t={\left(\begin{array}{cccc}1& e^{-\mathrm{j\pi }{f}_0}& ...& e^{-j(K-1)\pi f_0}\end{array}\right)}^T\\ A_s={\left(\begin{array}{cccc}1& e^{-\mathrm{j\pi }{f}_s}& ...& e^{-j(N-1)\pi f_s}\end{array}\right)}^T\end{array}\right)---\left(12\right)$

$B=B_t\otimes B_s$

上式中，B_t为K×K维矩阵，且B_s为N×N维全1矩阵。

风切变目标R_s的参数化模型可由其回波模型推导得到。低空风切变属于分布式目标， 其多普勒频率随距离的分布呈反“S”形，且单距离单元内的风切变信号在多普勒域内存 在连续性展宽。基于此种特性，设定待检测距离单元内多普勒中心频率和多普勒谱宽，可 以构造信号协方差矩阵某一距离单元内低空风切变回波模型可写为

$s_{n,k}=z_{n,k}{e}^{\mathrm{j\pi }(k-1)f+\mathrm{j\pi }(n-1)f_s},k=\mathrm{1,2},...,K,n=\mathrm{1,2},...,N---\left(13\right)$

其中z_n,k为接收信号幅度；f为信号归一化多普勒频率，且通常认为此 时，回波信号相关函数为

$r_{n,k}=E\left(s_{n_0,k_0}{s}_{n_0+n,k_0+k}^{*}\right)=\mathrm{PE}\left(e^{-\mathrm{j\pi kf}-\mathrm{j\pi n}{f}_s}\right)={\mathrm{Pe}}^{-\mathrm{j\pi n}{f}_s}E\left(e^{-\mathrm{j\pi kf}}\right)={\mathrm{Pe}}^{-\mathrm{j\pi n}{f}_s}{e}^{-\mathrm{j\pi k}{f}_0}{e}^{-\frac{1}{2}{\pi }^2{k}^2{\sigma }_f^2}---\left(14\right)$

由信号相关函数即可推得风切变目标的参数化模型如式(11)所示。

将信号协方差矩阵带入式(10)中，可以求得最优处理器权矢量ω，此处理器可以自 适应地抑制地杂波，并积累多普勒频率处于内的目标信号。地杂波抑制 与信号匹配的结果为

y＝ω^Hx_i                     (15)

更新求解权矢量ω时设定的多普勒中心频率f₀，得到对不同多普勒频段目标信号进行 匹配滤波的最优处理器权矢量，并分别求解处理器输出信号功率。当输出功率最大时，最 优处理器可以有效地抑制地杂波并积累低空风切变信号，此时设定的多普勒中心频率即为 待检测距离单元内风切变信号的多普勒频率估计值进而得到该距离单元的风速估计值 为

${\hat{v}}_i=\frac{\lambda {\hat{f}}_i{f}_r}{4}---\left(16\right)$

谱宽σ_f表示了低空风切变场的风速变化率，可取为实际观测中风速谱宽的统计平均 值，本步骤本质是固定空间锥角ψ₀和谱宽σ_f，利用空时最优处理器对待检测距离单元的 多普勒中心频率做一维搜索。

6)重复步骤2)-5)，依次对范围内所有距离单元的回波数据进行风速估计：

利用步骤2)-5)所述算法估计出每个距离单元的风场速度，即完成低空风切变场的 检测过程。

本发明提供的基于空时最优处理器的低空风切变风速估计方法的效果可以通过以下 仿真结果进一步说明。

仿真条件描述：低空风切变场分布于飞机前方约8.5-16.5km处，天线阵为阵元数为 N＝8的均匀线阵，阵元间距d＝λ/2，主瓣波束水平方位角为60°，俯仰角为0°，机载气 象雷达工作波长为0.05m，脉冲重复频率为7000Hz，最小可分辨距离150m，相干处理脉 冲数K＝64，信噪比5dB，杂噪比40dB，飞机速度为75m/s，飞行高度为600m，归一化 多普勒谱宽σ_f＝0.05。

图3为本发明提供的基于空时最优处理器的低空风切变风速估计方法的基本流程图。 其关键步骤在于：利用空时插值法校正地杂波的距离依赖性，构造符合待检测距离单元地 杂波二维分布特征的多个IID样本；设计自适应空时处理器，在抑制地杂波的同时对风切 变信号进行积累，运用目标检测的思想，搜索风切变信号的多普勒频率。

图4为补偿飞机速度后的仿真信号空时二维谱示意图，为便于观察，此时设置风切变 信号与地杂波功率相同。如图所示，地杂波的空时二维谱呈现椭圆形分布，低空风切变信 号的空时二维谱表现为主瓣波束空间锥角方向上的一条“窄带”。本发明提供方法利用空 时最优处理器自适应地在地杂波方向形成凹口，在风切变信号方向形成增益，进而得到风 场多普勒频率的有效估计结果。

图5(a)给出了利用空时插值法校正地杂波距离依赖性前的杂波功率谱，由于机载前视 阵雷达的地杂波分布具有距离依赖性，其空时二维谱有明显发散现象，直接使用回波信号 估计地杂波协方差矩阵，会致使最优处理器的地杂波抑制性能严重下降。图5(b)给出了利 用空时插值法向第20号距离单元做杂波谱补偿的结果，可以看出，补偿后的地杂波空时 二维谱明显变窄。

图6以20号距离单元为例，显示了f₀＝0.5,σ_f＝0.05时的最优处理器频响特性。可 以看出，最优处理器沿地杂波分布方向形成自适应凹口，同时对中心频率为f₀的宽带信号 进行功率积累。按f₀循环，更新最优处理器的权矢量，得到最大输出功率对应的即可 计算此距离单元的风场风速。

图7为本发明方法多普勒估计结果。如图所示，在8.5-16.5km范围内，信号的多普勒 频率随距离呈现反“S”形分布，进一步计算风速梯度和告警因子后可完成低空风切变的 检测过程。

图8为本发明方法与传统方法的风速估计结果对比。在同等杂噪比与信噪比条件下， 本发明方法处理效果明显优于自适应频域置零法与固定凹口法等传统方法。结果证明，本 发明提供的基于空时最优处理器的低空风切变风速估计方法可以在强杂波条件下，有效地 抑制地杂波，并精确估计风场速度。