一种机载雷达低空风切变风场回波风速的估计方法与流程

本发明涉及雷达信号处理技术领域，尤其是一种信号的估计方法的研究，特别是对于回波的风速估计。

背景技术：

风切变是一种大气现象，是风向和风速在水平和垂直方向突然变化的一种大气紊流，它主要是由大气运动本身的变化所造成的，同时地理和环境因素也会造成。600米以下空气层中风向和风速突然改变的现象称为低空风切变，它是航空界公认的飞机在起飞和着陆阶段时的“杀手”。低空风切变具有发生突然、时间短、尺度小、强度大等特点，为保证飞行安全，目前对付低空风切变最好的办法就是避开它，要避开风切变首先要能够探测或预测它出现的时间和地点，但是常规气象雷达很难完成此任务，从目前相关文献可见对于可以进行风切变探测和预警的机载前视风切变雷达的一系列研究，用以确保飞机在起飞和着陆时的安全。

目前已提出的用于低空风切变风场回波风速估计的算法，主要包括快速傅里叶变换法(fastfouriertransform，fft)、脉冲对处理法(pulsepairprocessing，ppp)，但是这两种方法均需要先对地杂波进行滤波，然后再进行风速估计，可能导致极点向地杂波方向的漂移，从而使得风场回波估计不够准确，不能够很好的反映出真实风速。

技术实现要素：

为了克服现有技术的不足，本发明提供一种机载雷达低空风切变风场回波风速的估计方法，用以克服现有方法存在的局限，使得无需先对地杂波进行滤波就可以进行低空风切变风场回波风速的估计，保证估计结果的准确，更好的反映真实风速。

本发明解决其技术问题所采用的技术方案包括下列步骤：

步骤1：定义高斯白噪声x(n)、距离门的数量n以及脉冲间隔ts三个变量，风切变雷达回波信号y(n)由白噪声通过p阶全极点滤波器建模产生，表示为：

式中，ak为滤波器的系数，k为系数向量；

步骤2：计算低空风切变雷达回波信号的自相关系数，若滤波器为m阶，则自相关系数r(m)的公式为：

其中，为x(n)的共轭，由此可得2阶ar模型里的第0、1、2个自相关系数r(0)、r(1)、r(2)分别为

步骤3：计算自回归(auto-regressive，ar)模型的系数

用levinson-durbin算法求解ar模型的系数，定义amk为m阶ar模型的第k个系数，ρm为m阶时的前向预测的最小误差功率，当ar模型为一阶时，有：

ρ1＝(1-|a11|²)r(0)(8)

定义第m阶时候的第m个系数为amm，取amm＝km，km称为反射系数，由toeplitz矩阵的性质，得到如下levinson-durbin递推算法：

amk＝a(m-1)k+kma(m-1)(m-k)(10)

其中a(m-1)k为第m-1阶时候的第k个系数，r(m-k)为第m-k个自相关系数；

则对于二阶的ar模型，系数估计由公式(9)-公式(11)得到，计算式如下：

σ²＝(1-|a11|²)r(0)(13)

其中，σ²是白噪声功率；

步骤4：计算风切变雷达回波功率谱

利用步骤3求出的二阶ar模型的系数计算风切变雷达回波功率谱，由于步骤1的y(n)是由白噪声通过滤波器得到的，因此噪声的谱函数就即滤波器的频率响应，滤波器的参数将完全决定p阶模型的功率谱sy(f)：

其中，σ²是白噪声功率，ai为滤波器的系数，f为频率自变量，i为系数向量，ts为脉冲采样间隔；

步骤5：计算方程的极点

滤波器的极点由下面二阶特征方程的根决定：

μ²+a21μ+a22＝0(17)

由公式(17)得到两个极点分别为：

步骤6：求出极点对应的速度

两个极点对应的速度v1，v2的计算公式为：

v1＝cμ1(20)

v2＝cμ2(21)

式中，为常量，用于多普勒频率f和频谱的速度v之间的转换，其中λ为波长，考虑到μ代表复指数exp(j2πfts)，因此有：

其中arg(μ1)、arg(μ2)分别为μ1、μ2的幅角；

步骤7：对极点的速度进行判断

若|v1|≥|v2|，则|v1|与|v2|值互换，且相应的|μ1|与|μ2|值也互换，并进入步骤8；否则，直接进入步骤8；

步骤8：利用分类曲线g(v)对速度进行判断

设计一个二阶多项式的判别曲线g(v)：

其中，v为多普勒速度，g(v)为多普勒速度对应的幅值，以公式(25)做判决：

其中g(v1)为速度v1对应的判别曲线，即当极点μ1幅值的平方小于等于判别曲线g(v1)时，该极点具有的多普勒速度v被置零；

步骤9：输出功率和平均风速

由步骤4得到回波信号功率谱p，由步骤8得到多普勒速度v，由此即可得到风切变雷达回波功率谱及风切变风场回波的风速估计值。

本发明的有益效果在于提出了一种机载雷达低空风切变风场回波风速估计方法，首先利用二阶自回归模型产生数据，再由二阶特征方程计算滤波器的极点，最后计算出极点对应的风速；该方法计算量小，克服了以往方案的局限，不需要对杂波进行滤波，保证了风场回波的准确性；不存在极点向地杂波方向的漂移，从而使结果更准确。从仿真结果容易验证，本发明提出的风速估计方法能够很好的反映出真实风速，方法有效、可靠。本发明方法同样适用于湍流情形和地杂波较强情况下的风切变信号处理。

附图说明

图1是本发明的流程图。

图2是顺风情况下本发明估计的回波风速与真实风速比较图。

图3是偏风情况下本发明估计的回波风速与真实风速比较图。

图4是逆风情况下本发明估计的回波风速与真实风速比较图。

图5是侧风情况下本发明估计的回波风速与真实风速比较图。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

图1是本发明的流程图，详细步骤如下：

步骤1：定义高斯白噪声x(n)、距离门的数量n以及脉冲间隔ts三个变量，风切变雷达回波信号y(n)由白噪声通过p阶全极点滤波器建模产生，表示为：

式中，ak为滤波器的系数，k为系数向量；

步骤2：计算低空风切变雷达回波信号的自相关系数，若滤波器为m阶，则自相关系数r(m)的公式为：

其中，为x(n)的共轭，由此可得2阶ar模型里的第0、1、2个自相关系数r(0)、r(1)、r(2)分别为

步骤3：计算自回归(auto-regressive，ar)模型的系数

用levinson-durbin算法求解ar模型的系数，定义amk为m阶ar模型的第k个系数，ρm为m阶时的前向预测的最小误差功率，当ar模型为一阶时，有：

ρ1＝(1-|a11|²)r(0)(8)

定义第m阶时候的第m个系数为amm，取amm＝km，km称为反射系数，由toeplitz矩阵的性质，得到如下levinson-durbin递推算法：

amk＝a(m-1)k+kma(m-1)(m-k)(10)

其中a(m-1)k为第m-1阶时候的第k个系数，r(m-k)为第m-k个自相关系数；

则对于二阶的ar模型，系数估计由公式(9)-公式(11)得到，计算式如下：

σ²＝(1-|a11|²)r(0)(13)

其中，σ²是白噪声功率；

步骤4：计算风切变雷达回波功率谱

利用步骤3求出的二阶ar模型的系数计算风切变雷达回波功率谱，由于步骤1的y(n)是由白噪声通过滤波器得到的，因此噪声的谱函数就即滤波器的频率响应，滤波器的参数将完全决定p阶模型的功率谱sy(f)：

其中，σ²是白噪声功率，ai为滤波器的系数，f为频率自变量，i为系数向量，ts为脉冲采样间隔；

步骤5：计算方程的极点

滤波器的极点由下面二阶特征方程的根决定：

μ²+a21μ+a22＝0(17)

由公式(17)得到两个极点分别为：

步骤6：求出极点对应的速度

两个极点对应的速度v1，v2的计算公式为：

v1＝cμ1(20)

v2＝cμ2(21)

式中，为常量，用于多普勒频率f和频谱的速度v之间的转换，其中λ为波长，考虑到μ代表复指数exp(j2πfts)，因此有：

其中arg(μ1)、arg(μ2)分别为μ1、μ2的幅角；

步骤7：对极点的速度进行判断

若|v1|≥|v2|，则|v1|与|v2|值互换，且相应的|μ1|与|μ2|值也互换，并进入步骤8；否则，直接进入步骤8；

步骤8：利用分类曲线g(v)对速度进行判断

设计一个二阶多项式的判别曲线g(v)：

其中，v为多普勒速度，g(v)为多普勒速度对应的幅值，以公式(25)做判决：

其中g(v1)为速度v1对应的判别曲线，即当极点μ1幅值的平方小于等于判别曲线g(v1)时，该极点具有的多普勒速度v被置零；

步骤9：输出功率和平均风速

由步骤4得到回波信号功率谱p，由步骤8得到多普勒速度v，由此即可得到风切变雷达回波功率谱及风切变风场回波的风速估计值，该方法可用于新型机载气象雷达的功能模块设计，结果可为雷达在飞机起飞和着陆阶段对低空风切变等的探测、预警提供依据。

由方案可知，本发明所提出的方法可以直接对风切变风场进行风速估计，无需事先对地杂波进行滤波等预处理，因此可避免因极点向地杂波方向漂移而使得风场回波估计不够准确、无法准确反映出真实风速的问题，弥补现有方法的不足。

参照实际装机使用的wxr-2100型机载气象雷达，基本的仿真参数设定如下：飞机高度600m，飞机速度75m/s，雷达发射功率100w，波束宽度3.5°，波束下视角2.75°，脉冲重复频率7000hz，发射脉宽1μs，天线增益30db，相干脉冲数128，杂噪比30db，信噪比5db。

下面以第30个距离门处的回波速度计算为例，以此类推可以得到全部50个距离门上的回波速度。需要说明的是，本方法可在任选适合的参数范围内实现，并不局限为以下的典型实现过程：

步骤1：设定n＝50，ts＝1μs,在matlab里产生一组回波信号，

步骤2：用matlab计算回波信号的自相关系数r0，r1，r2得

r(0)＝1.602

r(1)＝1.036

r(2)＝0.819

步骤3：

step1:用matlab产生均值为零、功率为0.01均匀分布的白噪声信号，过程为：

p＝0.01；n＝50000；

u＝rand(1,n),u＝u-mean(u)；

a＝sqrt(12*p)；u1＝u*a；

power_u1＝dot(u1,u1)/n

step2:计算出二阶ar模型的系数：

σ²＝(1-|a11|²)r(0)＝0.96

步骤4：代入n＝50，ts＝1μs，再结合之前得到的ar模型的系数，可计算任意阶数的风切变雷达回波功率谱

步骤5：解二阶特征方程：μ²+a21μ+a22＝0,

得到两个极点分别为：

μ1＝19.6,

μ2＝23.1。

步骤6：先计算再由v1＝cμ1，v2＝cμ2计算两个极点对应的速度：

得v1＝12.5

得v2＝10.3

步骤7：因为|v1|≥|v2|，将|v1|与|v2|值互换，同时|μ1|与|μ2|互换，得到v1＝10.3

步骤8：构造二阶多项式的判别曲线

代入v1＝10.3，可见g(v1)≤|μ1|²，依据得v＝10.3

步骤9：依据以上结果分别输出功率谱与平均风速值v＝10.3。

图2，图3，图4和图5是扫描线扫到风场中心，分别在顺风、逆风、偏风、侧风情况下第30根扫描线上测得的50个距离门上的回波风速跟采用本发明方法所估计得到的数据的对比，可见，在这四种情况下，本发明提出的风速估计方法能够很好的反映出真实风速，方法是有效可靠的。