基于迭代自适应处理的低空风切变风速估计方法与流程

本发明属于雷达信号处理技术领域，特别是涉及一种基于迭代自适应处理的低空风切变风速估计方法。

背景技术：

风切变是一种发生在很短一段距离内的微小气象现象，即在大气中相对较短的距离内的风速和(或)方向的急速变化。风切变尤其是低空风切变，对飞机的飞行(尤其是起降阶段)有重大影响。当飞机靠近风切变场时，首先进入逆风区域，此时飞机升力增加，飞行员会降低发动机功率，当飞机飞出逆风区域进入顺风区域时，局部逆风减弱，降低了飞机的空速并增加了下沉速率，这时飞机处于低功率快速下降状态。如果飞机太低而无法在接触地面之前增加升力，则可能导致事故。发生在低空的风切变是飞机起飞和着陆阶段的一个重要危险因素，被人们称为“无形杀手”。

机载气象雷达能够沿航路实时感知暴雨、湍流、风切变等有害自然现象并及时发出预警，可以称得上是飞机的“眼睛”。但当其处于“前视”模式时也使得雷达接收地杂波的空时分布与距离相耦合，此时地杂波产生了距离依赖性(即不同距离单元内杂波的空时分布不同)，这种非均匀特性使得可用于估计杂波协方差矩阵的iid样本数量大大减少。空时自适应处理技术主要用于机载气象雷达中抑制杂波及多种干扰，从而可以有效地检测空中与地(海)面的目标，其利用雷达回波中的时域信息和空域信息进行二维联合自适应处理，最大限度地提高了待检测目标的输出信噪比，具有很强大的灵活性。但由于系统维数较高，难以应用于实际工程中。减少空时自适应处理技术运算量的手段主要是通过设计次优的降维处理器以实现运算量级的降低，主要方法有辅助通道法、2维capon法、局域联合处理法等。目前，空时自适应处理技术主要被应用于点目标的参数估计上，而低空风切变属于分布式气象目标范畴。分布式信号伴随着源发生一定的角度扩展和频率扩展，往往比点信号源具有更复杂的空间分布特征。但目前尚缺少相应的方法。

技术实现要素：

为了解决上述问题，本发明的目的在于提供一种能够保证参数估计精度，同时降低运算复杂度的基于迭代自适应处理的低空风切变风速估计方法。

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

1)根据雷达接收数据，构造空时二维导向矢量，并获得空时导向矢量完备矩阵；

2)根据上述空时导向矢量完备矩阵对待检测距离单元的雷达接收数据进行迭代自适应处理，得到高分辨率空时二维功率谱；

3)利用邻近距离单元的空时二维谱功率确定地杂波功率水平，并设计中值滤波器，以检测风切变信号；

4)将上述经过中值滤波器滤波的信号利用频率质心法获得风切变场的中心风速估计值及分布情况。

在步骤1)中，所述的根据雷达接收数据，构造空时二维导向矢量，并获得空时导向矢量完备矩阵的方法是：首先，依据雷达接收数据设定角度域与多普勒频率域的扫描点数为雷达接收数据的慢时间-脉冲数的二倍，然后，依照角度域与多普勒频率域各个扫描点对应的角频率构造空时导向矢量，并用各个扫描点对应的导向矢量构造一个空时导向矢量完备矩阵。

在步骤2)中，所述的根据上述空时导向矢量完备矩阵对待检测距离单元的雷达接收数据进行迭代自适应处理，得到高分辨率空时二维功率谱的方法是：首先，使用待检距离单元的雷达接收数据与空时导向矢量以波束形成的方式初始化对应于空时导向矢量完备矩阵的功率对角矩阵的估计值，其次，通过空时导向矢量完备矩阵与功率对角矩阵来计算自适应协方差矩阵，然后，利用计算得到的协方差矩阵再重新估计功率对角矩阵，并反复进行迭代，最后，经过迭代多次得到一个收敛的自适应协方差矩阵。

在步骤3)中，所述的利用邻近距离单元的空时二维谱功率确定地杂波功率水平，并设计中值滤波器，以检测风切变信号的方法是：首先，选取待检测距离单元前后各3个邻近的距离单元作为保护单元，其次，选取待检测距离单元前后除保护单元外的各5个邻近的距离单元作为参考单元，并通过迭代自适应过程求得各自的功率对角矩阵，最后，使用各参考单元的功率对角矩阵的中值确定地杂波功率水平，并将地杂波功率水平作为滤波门限设计中值滤波器。

在步骤4)中，所述的将上述经过中值滤波器滤波的信号利用频率质心法获得风切变场的中心风速估计值及分布情况的方法是：首先，利用迭代自适应处理得到的自适应协方差矩阵估计待检测距离单元目标方向上的多普勒频谱，其次，利用频率质心法求出该距离单元内风切变场的中心风速估计值，最后，循环处理全部距离单元内所有多普勒通道的雷达接收数据，得到待检测区域内风切变场的中心风速分布情况。

本发明提供的基于迭代自适应处理的低空风切变风速估计方法无需抑制杂波，而是利用迭代算法估计协方差矩阵，从而重构出高分辨率的空时二维功率谱，由空时二维功率谱的高分辨率估计值直接经由中值滤波器后检测出风切变并估计其中心风速。本发明方法可以在低信噪比、强杂噪比条件下有效地估计出风切变中心风速，仿真实验验证了本方法的有效性，且由于无需抑制杂波，与最优处理器相比较，在运算量上有了明显减少。

附图说明

图1为本发明提供的基于迭代自适应处理的低空风切变风速估计方法流程图。

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

图3为本发明步骤2)中对待检测距离单元的雷达接收数据进行迭代自适应处理过程流程图。

图4为仿真所得的地杂波的最小方差谱图。

图5实际雷达接收信号谱图。

图6为直接使用第50号距离单元的雷达接收数据得出的空时二维功率谱图。

图7为利用本发明方法对第50号距离单元的雷达接收数据进行处理后得到的空时二维功率谱图。

图8为利用本发明方法对第50号距离单元的雷达接收数据进行处理后得到的多普勒谱图。

图9为本发明方法与固定凹口滤波器法、自适应频率置零法、广义capon方法及多通道方法的风速估计结果对比。

具体实施方法

下面结合附图和具体实例对本发明提供的基于迭代自适应处理的低空风切变风速估计方法进行详细说明。

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

1)根据雷达接收数据，构造空时二维导向矢量，并获得空时导向矢量完备矩阵；

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

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

其中，矩阵sl的第n行、第k列元素表示雷达第n(n＝1,2,…n)个阵元、第k(k＝1,2,…k)个脉冲对风切变场回波的采样，q表示该距离单元内雷达波束照射范围内气象散射粒子的个数。

式中：

和分别表示第q(q＝1,2,…,q)个气象散射粒子的空间角频率和时间角频率，θq，分别表示该气象散射粒子相对于雷达的方位角和俯仰角。aq表示第q个气象散射粒子的一次散射幅度，rq表示第q个散射粒子与载机的斜距。将矩阵sl展开成为nk×1维列向量，即为风切变场产生的雷达回波信号sl。

在本发明中，yl表示第l个距离单元的nk×1维雷达接收数据，其表达式如下：

yl＝sl+cl+nl(3)

其中，sl为第l个距离单元内风切变场产生的雷达回波信号，cl为地杂波信号，在此假设地杂波无起伏无模糊，nl为加性高斯白噪声。

将角度-多普勒轴均匀量化，假设角度和多普勒扫描(网格)点的数量分别为和则扫描的总点数为每次只处理一个距离单元，则空时导向矢量可以表示如下：

其中表示克罗内克积，表示空间导向矢量，表示时间导向矢量，(·)^t表示共轭转置。

定义是一个nk×m的空时导向矢量完备矩阵。

2)根据上述空时导向矢量完备矩阵对待检测距离单元的雷达接收数据进行迭代自适应处理，得到高分辨率空时二维功率谱；

根据步骤1)所得的空时导向矢量完备矩阵，可将式(3)所示的雷达接收数据写成如下形式：

yl＝axm+el(5)

其中xm＝[x0,x1,...,xm-1]，xm表示与频率相对应的复幅度，el为噪声项。

对于每一个待检测区域(regionofinterest,roi)，本发明只使用当前距离单元的数据，扫描角度和多普勒维，从而形成其角度-多普勒域的谱，例如要计算的目标以及地杂波和噪声的空时二维功率谱。

iaa算法(迭代自适应处理算法)的基本思想是优化如下所示的加权最小二乘问题：

其中表示除了网格点的信号分量以外的所有信号成分的协方差矩阵，具体可写成如下形式：

其中，pm＝|xm|²表示在网格点为的信号分量功率，r表示由离散化假设模型得到协方差矩阵，从而有iaa协方差矩阵的表达式:

riaa＝apma^h(8)

其中pm表示功率对角矩阵，其对角元素为[p0,p1,...,pm-1]，(·)^h表示共轭转置。应用矩阵求逆引理，我们获得复幅度xm的iaa估计值：

由iaa估计值xm^iaa进而求得功率对角矩阵pm的估计值：

并用重新估计所得的对角矩阵pm更新式(8)，并求得新的iaa协方差矩阵riaa。将此过程反复进行迭代直至收敛，最终获得迭代自适应处理结果，即iaa协方差矩阵：

riaa＝apma^h(11)

利用iaa协方差矩阵得到高分辨率的空时二维功率谱。本步骤的处理过程如图3所示。

3)利用邻近距离单元的空时二维谱功率确定地杂波功率水平，并设计中值滤波器，以检测风切变信号；

假设在固定的角度和给定的距离单元时，在与待检测区域相邻的几个距离单元内地杂波峰值几乎相同，这使得我们可以利用邻近距离单元的空时二维功率谱来确定地杂波功率水平。本步骤选取待检测距离单元前后各3个邻近的距离单元作为保护单元，选取待检测距离单元前后除保护单元外的各5个邻近的距离单元作为参考单元，并通过上述步骤2)中的迭代自适应处理过程求得各自的功率对角矩阵。在得到高分辨率的空时二维功率谱后，利用邻近距离单元的空时二维功率谱取中值来确定地杂波功率水平ηm(fd)，并将该地杂波功率水平作为滤波门限设计中值滤波器：

4)将上述经过中值滤波器滤波的信号利用频率质心法获得风切变场的中心风速估计值及分布情况。

进行风速估计时，由于我们已知了雷达主瓣的指向参数，所以式(13)中雷达接收数据yl与频率为已知,ωt(fd)为关于多普勒频率fd的函数。将上述步骤2)得到的最终达到收敛的iaa协方差矩阵riaa带入式(13)中,并遍历多普勒频率fd即可得到对应于各个多普勒频率的信号功率

依本步骤方法逐距离单元对雷达接收数据进行处理，选取该距离单元的信号功率的峰值所对应的风速即为中心风速，由各距离单元的中心风速即可得到风切变场速度随距离变化的分布情况，进而获得待检测区域内风切变场的中心风速分布情况。

本发明提供的基于迭代自适应处理的低空风切变风速估计方法的效果可以通过以下仿真结果进一步说明。

仿真参数设置：风切变分布于飞机前方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。

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

图6为直接使用第50号距离单元的雷达接收数据得出的空时二维功率谱图，图7为利用本发明方法对第50号距离单元的雷达接收数据进行处理后得到的空时二维功率谱图，从图中可明显看出地杂波的旁瓣和目标。由图6和图7可以看出，同样只使用单个距离单元的数据情况下，直接利用雷达接收数据得到的空时二维功率谱只能看到地杂波的主瓣，而经过本发明方法处理后所得的空时二维功率谱能明显看出地杂波的旁瓣和目标，即经过本发明方法处理后可得到高分辨率的空时二维功率谱。

图8为利用本发明方法对第50号距离单元的雷达接收数据进行处理后得到的多普勒谱图，从图中可明显看出经过本发明方法对雷达接收数据进行处理后所得的空时二维功率谱分辨率更高，更易进行目标检测和参数估计。

图9为本发明方法与固定凹口滤波器法、自适应频率置零法、广义capon方法及多通道方法的风速估计结果对比。

本发明方法与固定凹口滤波器法、自适应频率置零法、广义capon方法及多通道方法的风速估计结果对比如图9所示。可以看出在8.5～16.5km范围内，风速随距离呈反“s”变化；在相同杂噪比、信噪比下，本发明方法的风速估计结果与最优处理器结果性能相当，优于传统方法。