CN102360074A - 风廓线雷达谱分析方法 - Google Patents

Abstract

本发明涉及一种风廓线雷达谱分析方法，其包括：通过信号处理器输入时域采样信号；对该时域采样信号进行时域平均；对信号分别进行前向预测误差滤波以及后向预测误差滤波，并分别输出滤波器系数；计算预测误差能量；当预测误差能量开始增加时，输出当前滤波器系数，否则，将所得的滤波器系数代入进行递推，再次进行前向预测误差滤波以及后向预测误差滤波，并再次输出滤波器系数，直到预测误差能量开始增加；利用输出的滤波器系数将时域信号转化为频域信号，从而得到功率谱。本发明能够同时满足去噪性能和自适应性能方面的要求，提高了风廓线雷达谱分析的质量。

Description

风廓线雷达谱分析方法

技术领域

本发明涉及一种风廓线雷达谱分析方法，用于风廓线雷达的谱分析。

背景技术

风廓线雷达(wind profiling radar/wind profiler radar)是一种新型的测风雷达，能够无人值守24小时连续提供大气水平风场、垂直气流、大气折射率结构常数等气象要素随高度的分布，具有时空分辨率高、连续性和实时性好的特点，是进行高空气象探测的重要设备，是当前常规气球测风体制的重要补充，是开展天气预报和气象保障的新手段。

但是，由于风廓线雷达主要探测的是晴空湍流造成的回波，这种信号非常微弱，信噪比很低，雷达回波信号极易受到噪声的污染。噪声信号掺杂在微弱的雷达回波信号中会严重影响气象信息提取的精度。因此，开展风廓线雷达回波信号功率谱分析方法的研究是十分必要的。

风廓线雷达探测的目标物是弥散目标，回波信号是由它们对雷达波共同“散射”的结果，而每个湍流微团的移动速度未必一样，造成的多普勒频率也不一样，因此需要对回波信号进行多普勒频谱分析。只有得到具有一定信噪比的频谱，才有利于开展对目标的检测和矩数据的计算，也才有可能进一步计算风廓线。

目前风廓线雷达主要采用的就是经典信号谱分析技术，在用FFT法分析出功率谱之后，采用去直流、谱积累、谱对消和滤波等处理方法，对雷达噪声干扰进行抑制，但在很多情况下作用效果很有限。例如，消除零频附近的对称谱成分，有时会将重叠的信号部分也去除；用均值代替直流谱的方法会产生频率偏移问题。利用FFT变换将信号由时域变换到频域，一旦噪声频谱与雷达回波信号的频谱混杂在一起时，谱分析的结果就会是它们的集体呈现。当噪声较强时就难以得到有用的雷达回波信息，影响了谱的分析结果。因此，为提高数据质量，需要研究如何有效去除或抑制风廓线雷达信号中的噪声和地杂波，提高回波信号的信噪比。

目前多数研究都集中在如何将小波技术应用于风廓线雷达信号的杂波去除，并取得了一定的效果。但是，应用小波分析的主要困难在于需要找到一组合适的小波基，能够适应风廓线雷达的连续运行和无人值守式的工作，也就是要对环境具有自适应性，如果针对性太强，那就会限制其推广使用。

因此，需要研究一种新的谱分析方法，能够同时满足去噪性能和自适应性能方面的要求，以提高风廓线雷达谱分析的质量。

发明内容

本发明的目的在于提供一种风廓线雷达谱分析方法，该方法改善了雷达的功率谱计算，提高了信噪比，对随机白噪声有一定的抑制作用，另外，适用于各种气候条件下风廓线雷达信号谱的分析。

根据本发明的主要方面，该风廓线雷达谱分析方法包括以下步骤：

a、通过信号处理器输入时域采样信号；

b、对该时域采样信号进行时域平均；

c、对信号分别进行前向预测误差滤波以及后向预测误差滤波，并分别输出滤波器系数；

d、计算预测误差能量；

e、当预测误差能量开始增加时，输出当前滤波器系数，否则，将所得的滤波器系数代入进行递推，再次回到步骤c，重新依次执行后续步骤直到预测误差能量开始增加；

f、利用所输出的滤波器系数将时域信号转化为频域信号，得到功率谱。

根据本发明的一个方面，误差滤波的各阶系数的确定采用Burg递推算法。

根据本发明的一个方面，误差滤波的阶数的预测采用最终预测误差准则或信息论准则或自回归传输函数准则。

根据本发明的一个方面，对输入的时域采样信号进行A/D转换。

根据本发明的一个方面，对A/D转换后的信号还进行时域滤波。

应当认识到，本发明以上各方面中的特征可以在本发明的范围内自由组合，而并不受其顺序的限制——只要组合后的技术方案落在本发明的实质精神内。

附图说明

为了更清楚地说明本发明中的技术方案，下面将对本发明的附图作简单地介绍，其中：

图1显示了根据本发明的风廓线雷达谱分析方法的流程图。

具体实施方式

下文将结合本发明的优选实施例对本发明的技术方案进行详细说明。

需要明白，下文的描述(包括附图)仅仅是示例性的，而非对本发明的限制性描述。在以下描述中会涉及到部件的具体数量，然而也需要明白的是，这些数量也仅仅是示例性的，本领域技术人员可以参照本发明任意选取适当数量的部件。

根据本发明的一个实施例，所检测的风廓线雷达信号通过信号处理器输入，优选地，信号经过A/D转换处理；更优选的，还经过时域滤波处理；其中时域滤波器设计成线性相位滤波器。

然后，对信号处理器输入的信号进行时域平均，也即相干累积。由于噪声是随机信号，相位是随机的，因此噪声在累积中得不到额外的功率累加，因此，经相干累积后，输入信号的信噪比得到提高。

接着，对信号分别进行前向预测误差滤波以及后向预测误差滤波，并分别输出滤波器系数。“预测误差”是指用前次一个或多个观测值，采用某种预测模型对现时值作出的估计值

错误！未找到引用源。与现时实际观测值X(n)之差。所谓“预测误差滤波器”就是将现时观测值X(n)作为输入，而输出预测误差

的滤波器。预测误差滤波器一般由预测滤波器、延时器和减法器组成，预测滤波器采用历史观测值产生现时估计值

经延时一个采样周期后，与输入的现时观测值X(n)相减，然后输出预测误差

在本发明的一个实施例中，预测误差滤波器是具有最小相位的预测误差滤波器，该滤波器的各阶系数a(1，M)，a(2，M)，......，a(M，M)和滤波器的输出功率P_M，共同确定了输入数据序列{X(n)}的功率谱密度函数S(f)

$S\left(f\right)=\frac{P_M\mathrm{\Δt}}{{|1-\underset{m=1}{\overset{M}{\mathrm{\Σ}}}a(m,M)\exp (-j2\mathrm{\πfm\Δt})|}^2}---\left(1\right)$

优选地，使用Burg递推算法来给出各阶滤波系数；该算法在使滤波器的前向与后向预测误差能量之和为最小的约束下，给出了各阶滤波系数的递推方法；其不仅具有简单的解析求算形式，而且直接利用观测数据序列{X(n)}来计算滤波器系数。

对一阶预测误差滤波器，如使用前向的数据序列，且不超出观测数据序列，则预测误差输出为

f₁(n)＝X(n+1)-a(1，1)X(n)    (2)

其中n＝0，1，2，……，N-2。

若使用后向的数据序列，则预测误差输出为

b₁(n)＝X(n)-a^*(1，1)X(n+1)   (3)

则一阶预测误差滤波器的前后向输出能量之和为

$e_1=\frac{1}{2}\underset{n=0}{\overset{N-2}{\mathrm{\Σ}}}({\left|f_1\left(n\right)\right|}^2+{\left|b_1\left(n\right)\right|}^2)---\left(4\right)$

本发明的方法要求滤波器的输出能量为最小，则e₁对a(1，1)求最小，得

$a\left(\mathrm{1,1}\right)=\frac{2\underset{n=0}{\overset{N-2}{\mathrm{\Σ}}}{X}^{*}\left(n\right)X(n+1)}{\underset{n=0}{\overset{N-2}{\mathrm{\Σ}}}({\left|X\left(n\right)\right|}^2+{\left|X(n+1)\right|}^2)}---\left(5\right)$

对于二阶预测误差滤波器情况，此时预测误差输出分别是

f₂(n)＝X(n+2)-a(1，2)X(n+1)-a(2，2)X(n)   (6)

b₂(n)＝X(n)-a^*(1，2)X(n+1)-a^*(2，2)X(n+2)   (7)

其中n＝0，1，2，……，N-3。前、后向预测误差滤波器输出能量之和为

$e_2=\frac{1}{2}\underset{n=0}{\overset{N-3}{\mathrm{\Σ}}}({\left|f_2\left(n\right)\right|}^2+{\left|b_2\left(n\right)\right|}^2)---\left(8\right)$

为了能够确定出两个系数a(1，2)和a(2，2)，Burg应用了莱文森(Levinson)约束关系，即

a(1，2)＝a(1，1)-a(2，2)a^*(1，1)            (9)

可以解得

$a(2,2)=\frac{2\underset{n=0}{\overset{N-3}{\mathrm{\Σ}}}{b}_1^{*}\left(n\right)f_1(n+1)}{\underset{n=0}{\overset{N-3}{\mathrm{\Σ}}}({\left|b_1\left(n\right)\right|}^2+{\left|f_1(n+1)\right|}^2)}---\left(10\right)$

若令f₀(n)＝b₀(n)＝X(n)，代入公式(5)后，可以发现与(10)式很相似。于是可以将其推广到一般情况，得到如下递推公式

f_K(n)＝f_K-1(n+1)-a(K，K)b_K-1(n)        (11)

b_K(n)＝b_K-1(n)-a^*(K，K)f_K-1(n+1)       (12)

a(k，K)＝a(k，K-1)-a(K，K)a^*(K-1，K-1) (13)

$a(K,K)=\frac{2\underset{n=0}{\overset{N-1-K}{\mathrm{\Σ}}}{b}_{K-1}^{*}\left(n\right)f_{k-1}(n+1)}{\underset{n=0}{\overset{N-1-K}{\mathrm{\Σ}}}({\left|b_{K-1}\left(n\right)\right|}^2+{\left|f_{K-1}(n+1)\right|}^2)}---\left(14\right)$

其中K＝1，2，……，M。k＝1，2，……，K-1。K阶预测误差滤波器的输出能量的递推公式为

$P_0=\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{\left|X\left(n\right)\right|}^2---\left(15\right)$

P_K＝P_K-1(1-|a(K，K)|²)     (16)

公式(11)-(16)组成了Burg递推算法的计算公式，随着K从1逐渐递推到某个合适的最高阶M，用上述公式便可以计算出所有的预测滤波器系数。

在首次对信号分别进行前向预测误差滤波以及后向预测误差滤波，并分别输出滤波器系数后，需计算预测误差能量；如果预测误差能量开始增加，则输出当前滤波器系数，否则，将所得的滤波器系数代入进行递推，再次分别进行前向预测误差滤波以及后向预测误差滤波，并分别输出滤波器系数，直到预测误差能量开始增加。

最后，将时域信号转化为频域信号，其中根据所得到的滤波器系数利用公式(1)计算出功率谱。

另外，在本发明的方法中，误差滤波器阶数的预测是一个很重要的问题，它不仅影响计算时间，还关系到谱估计的质量。阶数小，计算时间短，但是谱会被严重平滑，降低了谱的分辨率；而阶数过大，谱分辨率会提高，但是计算时间长，而且可能会带来虚假的谱线分裂。所以有个最佳确定的问题。根据本发明的实施例，优选选择以下三个准则之一来进行误差滤波器阶数预测：

其一、最终预测误差准则(FPE准则)：其定义式为

$\mathrm{FPE}\left(K\right)=\frac{N+K+1}{N-K-1}{P}_K---\left(17\right)$

式中P_K为预测误差滤波器的输出功率，由公式(16)可知P_K随K的增加而减少，而公式(17)中右边前半部分随K的增加而增加，则FPE(K)必在某个K值处取得最小值，此时的K值就为所要确定的阶数。

其二、信息论准则(ALC准则)：其定义式为

ALC(K)＝ln(P_K)+2K/N    (18)

同样取使ALC(K)为最小的K值作为最后的阶数。

其三、自回归传输函数准则(CAT准则)：其定义式为

$\mathrm{CAT}\left(K\right)=\frac{1}{N}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\frac{N-k}{N{P}_k}-\frac{N-K}{N{P}_K}---\left(19\right)$

同样取使CAT(K)为最小的K值作为最后的阶数。

以上根据优选的实施方式对本发明作了详细的描述，不过需要理解的是，本发明的范围并不局限于这些具体的实施方式，而是包括本领域技术人员根据本发明的公开能够做出的任何改动和变更。

Claims (5)

1.一种风廓线雷达谱分析方法，所述方法包括以下步骤：

a、通过信号处理器输入时域采样信号；

b、对该时域采样信号进行时域平均；

c、对信号分别进行前向预测误差滤波以及后向预测误差滤波，并分别输出滤波器系数；

d、计算预测误差能量；

e、当预测误差能量开始增加时，输出当前滤波器系数，否则，将所得的滤波器系数代入进行递推，再次回到步骤c，重新依次执行后续步骤直到预测误差能量开始增加；

f、利用所输出的滤波器系数将时域信号转化为频域信号，得到功率谱。

2.根据权利要求1所述的风廓线雷达谱分析方法，其特征在于，所述递推采用Burg递推算法。

3.根据权利要求1所述的风廓线雷达谱分析方法，其特征在于，所述前向预测误差滤波和所述后向预测误差滤波的阶数的预测采用最终预测误差准则或信息论准则或自回归传输函数准则。

4.根据权利要求1所述的风廓线雷达谱分析方法，其特征在于，在所述信号处理器中，对时域采样信号进行A/D转换。

5.根据权利要求4所述的风廓线雷达谱分析方法，其特征在于，对A/D转换后的信号还进行时域滤波。