CN101334469B - 基于分数阶傅立叶变换的风廓线雷达杂波抑制方法 - Google Patents

Abstract

本发明涉及一种杂波抑制方法，具体涉及一种基于分数阶傅立叶变换的风廓线雷达杂波抑制方法，属于信号处理和检测技术领域。本发明的方法包括以下步骤：首先，选定最优的分数阶次α_opt；其次，对经相干积累后的时域数据在最优分数阶α_opt下进行分数阶傅立叶变换，得到相应的分数阶域信号形式；然后，以接收信号在分数阶域内的各阶矩的统计特性为依据来确定杂波干扰检测阈值；最后对抑制后的数据进行逆分数阶傅立叶变换，获得去除杂波后的时域风廓线雷达数据，完成基于分数阶傅立叶变换的杂波抑制处理。经过实测，经此方法处理的风廓线雷达接收信号能够显著的改善信噪比。

Description

基于分数阶傅立叶变换的风廓线雷达杂波抑制方法

技术领域

本发明是一种杂波抑制方法，是一种风廓线雷达在分数阶傅立叶域抑制杂波干扰的方法，属于信号处理和检测技术领域。

背景技术

风廓线雷达是对大气风向风速随高度分布连续观测的一种全天候型的Doppler遥感设备，其通过检测大气湍流反射的微弱信号来实现测风的目的。该技术的研究始于美国20世纪八十年代初，由于其具有较高的测风实时性和精度，是中小尺度大气探测的关键技术之一，风廓线雷达在气象、环保、军事等领域具有广阔的应用和研究前景。风廓线雷达回波信号功率较弱(甚至小于-140dBm)，而地杂波和间歇性杂波的幅度要高于风信号数个量级，是造成风廓线雷达有效数据率偏低的主要原因，更是限制其在机场、环保等领域应用的主要原因，如何有效抑制杂波而尽可能减少有用信号的损失是风廓线雷达信号处理亟待解决的难题之一。目前风廓线雷达主要采用经典的信号处理技术，例如脉冲压缩编码、去直流、谱积累、谱对消等，虽然对有源干扰有一定抑制作用，但信干比改善不大，效果不明显。

发明内容

本发明的目的在于将分数阶傅立叶变换与风廓线雷达信号处理结合，提供一个解决风廓线雷达杂波干扰问题的新方法，该方法能有效解决风廓线雷达实际应用的杂波干扰问题。

风廓线信号具有调频信号的特点，飞机等产生的间歇性杂波具有类似线性调频信号的特点，大气湍流回波、地杂波、间歇性杂波具有不同的时频相关性。基于此信号特点，本发明提出在寻找到最优变换阶次后，在相应的分数阶傅立叶域利用最佳杂波抑制准则，进行杂波干扰抑制的方法。

本发明的基于分数阶傅立叶变换的风廓线雷达杂波抑制方法，具体步骤为：

步骤一：选定最优的分数阶次α_opt；

步骤二：对经相干积累后的时域数据在最优分数阶α_opt下进行分数阶傅立叶变换，得到相应的分数阶域信号形式；

步骤三：以接收信号在分数阶域内的各阶矩的统计特性为依据来确定杂波干扰检测阈值，根据阈值设定在分数阶域内进行杂波抑制处理，将超过阈值的分数域数据识别为杂波并进行抑制；

步骤四：对抑制后的数据进行逆分数阶傅立叶变换，获得去除杂波后的时域风廓线雷达数据，完成基于分数阶傅立叶变换的杂波抑制处理。

经过实测，经此方法处理的风廓线雷达接收信号能够显著的改善信噪比。

附图说明

图1是本发明的基于分数阶傅立叶变换的风廓线雷达杂波抑制方法处理流程。

图2是本发明的基于分数阶傅立叶变换的风廓线雷达杂波抑制方法实现原理框图。

具体实施方式

图1是基于分数阶傅立叶变换的风廓线雷达杂波抑制方法处理流程，图2是基于分数阶傅立叶变换的风廓线雷达杂波抑制方法实现原理框图。风廓线雷达监控组合控制射频调制器、发射机和信号处理机。本实施例的具体实现过程为：

步骤一：选定最优的分数阶次α_opt(取值在0到π之间)，原则是使分数阶傅立叶变换的基函数与杂波信号能够很好的匹配，经过分数阶变换后杂波信号的特征可以很好的表现，使得有用信号分量和杂波分量在分数阶傅立叶域内可以得到很好的区分。

分数阶傅立叶变换的定义式如下：

$F_{\mathrm{\α}}\left(u\right)=\sqrt{\frac{1-j\cot \mathrm{\α}}{2\mathrm{\π}}}{e}^{\mathrm{j\π}\left(u^2\cot \mathrm{\α}\right)}\×{\∫}_{-\∞}^{+\∞}x\left(t\right)\·e^{\mathrm{j\π}\left(t^2\cot \mathrm{\α}\right)}\·e^{-\mathrm{j\π}\left(2\mathrm{ut}\csc \mathrm{\α}\right)}\mathrm{dt},\mathrm{\α}\≠\mathrm{n\π}---\left(1\right)$

其中，α为分数阶次；x(t)为经相干积累后的时域数据，即输入数据；F_α(u)为变换后的分数域信号。

则最优分数阶次α_opt的分数阶傅立叶变换峰值和带宽

应满足(2)和(3)：

$\underset{u}{\arg \max }\left|F_{{\mathrm{\α}}_{\mathrm{opt}}}\left(u\right)\right|\≥\left|F_{\mathrm{\α}}\left(u\right)\right|\∀\mathrm{\α}\∈\left(\mathrm{0,2}\mathrm{\π}\right)---\left(2\right)$

${\mathrm{\σ}}_{{\mathrm{\α}}_{\mathrm{opt},}u}\≤{\mathrm{\σ}}_{\mathrm{\α},u}\∀\mathrm{\α}\∈\left(\mathrm{0,2}\mathrm{\π}\right)---\left(3\right)$

通过代入α在0到π之间的不同取值，搜索满足式(2)和式(3)的取值，即为最优的分数阶次α_opt。

步骤二：对经相干积累后的时域数据进行分数阶傅立叶变换，由于风廓线雷达回波信号非常弱，信噪比太低，经过长时间的相干积累可以有效改善信噪比，同时降低分数阶傅立叶变换的数据率。应用H.M.Ozaktas等人提出的FRFT的分解算法，在最优分数阶α_opt下进行分数阶傅立叶变换，得到相应的分数阶域信号形式。

应用FRFT的分解算法，则公式(1)可以通过对连续分数阶傅立叶变换采样转换为(4)式：

$F_{\mathrm{\α}}\left(\frac{m}{2\mathrm{\Δx}}\right)=\frac{A_{\mathrm{\α}}}{2\mathrm{\Δx}}\underset{n=-N}{\overset{N}{\mathrm{\Σ}}}\exp [\frac{\mathrm{j\π}\left(\cot \mathrm{\α}\right)m^2}{{\left(2\mathrm{\Δx}\right)}^2}-\frac{j2\mathrm{\π}\left(\csc \mathrm{\α}\right)\mathrm{mn}}{{\left(2\mathrm{\Δx}\right)}^2}+\frac{\mathrm{j\π}\left(\cot \mathrm{\α}\right)n^2}{{\left(2\mathrm{\Δx}\right)}^2}]x\left(\frac{n}{2\mathrm{\Δx}}\right)---\left(4\right)$

其中， $A_{\mathrm{\α}}=\sqrt{\frac{1-j\cot \mathrm{\α}}{2\mathrm{\π}}},$ Δx为信号x(t)的带宽，N＝(Δx)²，函数exp为指数函数。按上式直接计算的计算复杂度为O(N²)，为了使计算复杂度减小为O(NlogN)。(4)式可以写作

$F_{\mathrm{\α}}\left(\frac{m}{2\mathrm{\Δx}}\right)=\frac{A_{\mathrm{\α}}}{2\mathrm{\Δx}}\exp \left[\frac{\mathrm{j\π}(\cot \mathrm{\α}-\csc \mathrm{\α})m^2}{{\left(2\mathrm{\Δx}\right)}^2}\right]\underset{n=-N}{\overset{N}{\mathrm{\Σ}}}\exp \left[\frac{j\mathrm{\π}\left(\cot \mathrm{\α}\right){(m-n)}^2}{{\left(2\mathrm{\Δx}\right)}^2}\right]\exp \left[\frac{\mathrm{j\π}(\cot \mathrm{\α}-\csc \mathrm{\α})n^2}{{\left(2\mathrm{\Δx}\right)}^2}\right]x\left(\frac{n}{2\mathrm{\Δx}}\right)$

此式的求和实际上是两个信号的卷积运算。该卷积运算可以用FFT计算，输出样本值然后可以通过最后的线性调频调制求得。因此，总的计算复杂度应为O(NlogN)。

实际工程应用中通常使用矩阵的方式来进行计算。用

表示x(t)的N个样本的列向量，表示分数阶傅立叶变换的N个样本的列向量，则分解算法用矩阵形式表示为 $\stackrel{\→}{X_{\mathrm{\α}}}=F^{\mathrm{\α}}\·\stackrel{\→}{x},$ 其中，F^α＝DK_pJ，D和J分别表示内插和抽取运算的矩阵，K_p(m，n)表示如(5)式所示：

$K_p(m,n)=\frac{A_{\mathrm{\α}}}{2\mathrm{\Δx}}\exp [\frac{\mathrm{j\π}\left(\cot \mathrm{\α}\right)m^2}{{\left(2\mathrm{\Δx}\right)}^2}-\frac{j2\mathrm{\π}\left(\csc \mathrm{\α}\right)\mathrm{mn}}{{\left(2\mathrm{\Δx}\right)}^2}+\frac{\mathrm{j\π}\left(\cot \mathrm{\α}\right)n^2}{{\left(2\mathrm{\Δx}\right)}^2}],\left|m\right|,\left|n\right|\≤N---\left(5\right)$

步骤三：以接收信号在分数阶域内的各阶矩的统计特性为依据来确定杂波干扰检测阈值。例如可以使用一阶矩法和二阶矩法实现此步骤，具体过程分别为：

一阶矩法

设检测阈值为Level，从雷达接收信号在分数域的均值特征判断何处为杂波干扰信号即为一阶矩法。

$\mathrm{Level}=\mathrm{\θ}\·\mathrm{\μ}=\mathrm{\θ}\·\frac{1}{N}\·\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\left|F_{\mathrm{\α}}\left(\frac{m}{2\mathrm{\Δx}}\right)\right|---\left(6\right)$

其中 $\mathrm{\μ}=\frac{1}{N}\·\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\left|F_{\mathrm{\α}}\left(\frac{m}{2\mathrm{\Δx}}\right)\right|$ 为幅度均值的估计值，N为FRFT的点数，

为谱线幅度值，θ为门限优化因子。

二阶矩法

$\mathrm{Level}=\mathrm{\μ}+\mathrm{\β}\·\mathrm{\σ}=\mathrm{\μ}+\mathrm{\β}\·\sqrt{\frac{1}{N-1}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left(\right|F_{\mathrm{\α}}\left(\frac{m}{2\mathrm{\Δx}}\right)|-\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}|F_{\mathrm{\α}}\left(\frac{m}{2\mathrm{\Δx}}\right)\left|\right)}^2}---\left(7\right)$

其中μ同上， $\mathrm{\σ}=\sqrt{\frac{1}{N-1}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left(\right|F_{\mathrm{\α}}\left(\frac{m}{2\mathrm{\Δx}}\right)|-\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}|F_{\mathrm{\α}}\left(\frac{m}{2\mathrm{\Δx}}\right)\left|\right)}^2}$ 为幅度标准偏差的估计值，β为门限优化因子。

根据阈值设定在分数阶域内进行杂波抑制处理，将超过阈值的分数域数据识别为杂波并进行抑制。

步骤四：根据公式(1)对抑制后的数据进行逆分数阶傅立叶变换，即求得经过处理后的时域数据，至此获得去除杂波后的时域风廓线雷达数据，完成基于分数阶傅立叶变换的杂波抑制处理。

至此，基于分数阶傅立叶变换的风廓线雷达杂波抑制方法处理完毕。

在实际工程中采用了二阶矩法来确定阈值，门限优化因子β取值为2，经以上两次运算后即可获得干扰抑制的处理结果，该结果可以作为风廓线雷达后续信号处理的输入信号。经此方法处理的风廓线雷达接收信号，经过实测数据显示能够改善信噪比达20dB以上。

Claims (6)

1.一种基于分数阶傅立叶变换的风廓线雷达杂波抑制方法，其特征在于包括以下步骤：

步骤一：选定最优的分数阶次α_opt；

步骤二：对经相干积累后的时域数据在最优分数阶α_opt下进行分数阶傅立叶变换，得到相应的分数阶域信号形式；

步骤三：以接收信号在分数阶域内的各阶矩的统计特性为依据来确定杂波干扰检测阈值，根据阈值设定在分数阶域内进行杂波抑制处理，将超过阈值的分数域数据识别为杂波并进行抑制；

步骤四：对抑制后的数据进行逆分数阶傅立叶变换，获得去除杂波后的时域风廓线雷达数据，完成基于分数阶傅立叶变换的杂波抑制处理。

2.根据权利要求1所述的风廓线雷达杂波抑制方法，其特征在于：步骤一中所述选定最优的分数阶次α_opt的方法为：

使用公式(1)进行分数阶傅立叶变换：

$F_{\mathrm{\α}}\left(u\right)=\sqrt{\frac{1-j\cot \mathrm{\α}}{2\mathrm{\π}}}{e}^{\mathrm{j\π}\left(u^2\cot \mathrm{\α}\right)}\×{\∫}_{-\∞}^{+\∞}x\left(t\right)\·e^{\mathrm{j\π}\left(t^2\cot \mathrm{\α}\right)}\·e^{-\mathrm{j\π}\left(2\mathrm{ut}\csc \mathrm{\α}\right)}\mathrm{dt},\mathrm{\α}\≠\mathrm{n\π}---\left(1\right)$

其中α为分数阶次；x(t)为经相干积累后的时域数据，即输入数据；F_α(u)为变换后的分数域信号；则最优分数阶次α_opt的分数阶傅立叶变换峰值和带宽应满足(2)和(3)：

$\underset{u}{\arg \max }\left|F_{{\mathrm{\α}}_{\mathrm{opt}}}\left(u\right)\right|\≥\left|F_{\mathrm{\α}}\left(u\right)\right|,\∀\mathrm{\α}\∈\left(\mathrm{0,2}\mathrm{\π}\right)---\left(2\right)$

${\mathrm{\σ}}_{{\mathrm{\α}}_{\mathrm{opt}},u}\≤{\mathrm{\σ}}_{\mathrm{\α},u}$ $\∀\mathrm{\α}\∈\left(\mathrm{0,2}\mathrm{\π}\right)---\left(3\right)$

通过代入α在0到π之间的不同取值，搜索满足式(2)和式(3)的取值，即为最优的分数阶次α_opt。

3.根据权利要求1所述的风廓线雷达杂波抑制方法，其特征在于：在步骤二中应用FRFT的分解算法，在最优分数阶α_opt下进行分数阶傅立叶变换，即通过公式

$F_{\mathrm{\α}}\left(\frac{m}{2\mathrm{\Δx}}\right)=\frac{A_{\mathrm{\α}}}{2\mathrm{\Δx}}\exp \left[\frac{\mathrm{j\π}(\cot \mathrm{\α}-\csc \mathrm{\α})m^2}{{\left(2\mathrm{\Δx}\right)}^2}\right]\underset{n=-N}{\overset{N}{\mathrm{\Σ}}}\exp \left[\frac{\mathrm{j\π}\left(\cot \mathrm{\α}\right){(m-n)}^2}{{\left(2\mathrm{\Δx}\right)}^2}\right]\exp [\frac{{\mathrm{j\π}(\cot \mathrm{\α}-\csc \mathrm{\α})n}^2}{{\left(2\mathrm{\Δx}\right)}^2}x\left(\frac{n}{2\mathrm{\Δx}}\right)$

得到相应的分数阶域信号形式，其中，

Δx为信号x(t)的带宽，N＝(Δx)²，函数exp为指数函数，α为分数阶次。

4.根据权利要求1所述的风廓线雷达杂波抑制方法，其特征在于：步骤三使用一阶矩法，通过公式(6)确定杂波干扰检测阈值：

$\mathrm{Level}=\mathrm{\θ}\·\mathrm{\μ}=\mathrm{\θ}\·\frac{1}{N}\·\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\left|F_{\mathrm{\α}}\left(\frac{m}{2\mathrm{\Δx}}\right)\right|---\left(6\right)$

其中

为幅度均值的估计值，N为FRFT的点数，

为谱线幅度值，θ为门限优化因子，α为分数阶次。

5.根据权利要求1所述的风廓线雷达杂波抑制方法，其特征在于：步骤三使用二阶矩法，通过公式(7)确定杂波干扰检测阈值：

$\mathrm{Level}=\mathrm{\μ}+\mathrm{\β}\·\mathrm{\σ}=\mathrm{\μ}+\mathrm{\β}\·\sqrt{\frac{1}{N-1}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left(\right|F_{\mathrm{\α}}\left(\frac{m}{2\mathrm{\Δx}}\right)|-\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}|F_{\mathrm{\α}}\left(\frac{m}{2\mathrm{\Δx}}\right)\left|\right)}^2}---\left(7\right)$

其中为幅度均值的估计值，

为幅度标准偏差的估计值，N为FRFT的点数，β为门限优化因子，α为分数阶次。

6.根据权利要求3所述的风廓线雷达杂波抑制方法，其特征在于：应用FRFT的分解算法，在最优分数阶α_opt下进行分数阶傅立叶变换使用矩阵的方式来进行计算：用

表示x(t)的N个样本的列向量，

表示分数阶傅立叶变换的N个样本的列向量，则分解算法用矩阵形式表示为

其中，F^α＝DK_pJ，D和J分别表示内插和抽取运算的矩阵，K_p(m，n)按照公式

$K_p(m,n)=\frac{A_{\mathrm{\α}}}{2\mathrm{\Δx}}\exp [\frac{{\mathrm{j\π}\left(\cot \mathrm{\α}\right)m}^2}{{\left(2\mathrm{\Δx}\right)}^2}-\frac{j2\mathrm{\π}\left(\csc \mathrm{\α}\right)\mathrm{mn}}{{\left(2\mathrm{\Δx}\right)}^2}+\frac{\mathrm{j\π}{\left(\cot \mathrm{\α}\right)n}^2}{{\left(2\mathrm{\Δx}\right)}^2}],\left|m\right|,\left|n\right|\≤N---\left(5\right)$

进行计算得到，其中

Δx为信号x(t)的带宽。

Images