CN104793211A - 一种风廓线雷达干扰抑制方法 - Google Patents

Abstract

本发明公开了一种风廓线雷达干扰抑制方法，包括以下步骤：步骤一：选定最优的分数阶次α_opt；步骤二：对经相干积累后的时域数据在最优分数阶α_opt下进行分数阶傅立叶变换，得到相应的分数阶域信号形式；步骤三：以接收信号在分数阶域内的各阶矩的统计特性为依据来确定杂波干扰检测阈值；步骤四：对抑制后的数据进行逆分数阶傅立叶变换，获得去除杂波后的时域风廓线雷达数据，完成杂波抑制处理；风廓线雷达接收信号经过此方法处理后，能够降低无线电信号对风廓线雷达的影响，且能够显著的改善信噪比，实用性强。

Description

一种风廓线雷达干扰抑制方法

技术领域

本发明涉及风廓线雷达技术领域，具体是一种风廓线雷达干扰抑制方法。

背景技术

风廓线雷达是通过向高空发射不同方向的电磁波束，接收并处理这些电磁波束因大气垂直结构不均匀而返回的信息进行高空风场探测的一种遥感设备。风廓线雷达利用多普勒效应能够探测其上空风向、风速等气象要素随高度的变化情况，具有探测时空分辨率高、自动化程度高等优点。在风廓线雷达基础上增加声发射装置构成无线电——声探测系统，可以遥感探测大气中温度的垂直廓线。风廓线雷达在气象、环保、军事等领域具有广阔的应用和研究前景。风廓线雷达回波信号功率较弱，而地杂波和间歇性杂波的幅度要高于风信号数个量级，造成风廓线雷达有效数据率偏低，另外，由于风廓线雷达与广播通信等设备的工作频段比较接近，再加上风廓线雷达回波信号功率较弱，所以容易受到无线电波的干扰，如何有效抑制杂波及无线电波对风廓线雷达的影响是亟待解决的难题之一。目前风廓线雷达主要采用经典的信号处理技术，例如脉冲压缩编码、去直流、谱积累、谱对消等，虽然对有源干扰有一定抑制作用，但信干比改善不大，效果不明显。

发明内容

本发明的目的在于提供一种风廓线雷达干扰抑制方法，以解决上述背景技术中提出的问题。

为实现上述目的，本发明提供如下技术方案：

一种风廓线雷达干扰抑制方法，包括以下步骤：

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

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

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

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

步骤五：在风廓线雷达发射端，同步产生脉间、脉内两种伪随机码序列，其中脉间伪随机序列采用M序列，序列长度根据实际工作参数选择，假设其码长为L，脉内编码序列采用巴克码，假设其码长为P，设脉间伪随机编码的码本为R：

R＝[R₁ R₂ … R_i … R_L]^T        (8)

式中R_i为M序列伪随机码码元，i＝1，2...L，T表示矩阵转置；脉内编码码本为巴克码，其位数为P，

BK＝[bk(1) bk(2) … bk(i) … bk(P)]^T     (9)

式中码元为bk(i)，i＝1，2...P，T表示矩阵转置；

步骤六：在风廓线雷达发射端，上述公式(8)和公式(9)码本序列经脉间脉内联合编码序列为：

$\mathrm{RB}=\left[\begin{array}{c}R{B}_1\\ R{B}_2\\ .\\ .\\ .\\ R{B}_i\\ .\\ .\\ .\\ R{B}_L\end{array}\right]---\left(10\right)$

式中RB_i为编码后第i个脉冲内的码本信号，i＝1，2...L，

$R{B}_i=(\mathrm{BK}+R_i)\mathrm{mod}\left(2\right)=\left[\begin{array}{c}\mathrm{bk}\left(1\right)\mathrm{xor}{R}_i\\ \mathrm{bk}\left(2\right)\mathrm{xor}{R}_i\\ .\\ .\\ .\\ \mathrm{bk}\left(P\right)\mathrm{xor}{R}_i\end{array}\right]---\left(11\right)$

式中mod()为求余运算，xor为异或运算符；

步骤七：将公式(10)产生的码本序列传送至风廓线雷达发射信号调制器并发射，同时将公式(8)和公式(9)产生的码本序列保存在接收信号处理机中；

步骤八：在风廓线雷达接收端对步骤七中的发射信号进行接收，设时域积累为J次，FFT点数为N，距离门数为G，则信号处理机需处理的信号样本空间I为：

$I={\left[\begin{array}{cccc}{A}_{11}& A_{12}& ...& A_{1,J}\\ A_{21}& A_{22}& ...& A_{2J}\\ .& .& & .\\ .& .& A_{\mathrm{ij}}& .\\ .& .& & .\\ A_{N1}& A_{N2}& ...& A_{\mathrm{NJ}}\end{array}\right]}_{N\×J}---\left(12\right)$

式中A_ij分别代表第i×j次采样信号a_ij(1)，a_ij(2)，...a_ij(G+P-1)构成的样本空间，其中i＝1，2...N，j＝1，2...J分别与信号处理中的FFT点数和相干积累次数有关，L＝N×J；

A_ij＝[a_ij(1)a_ij(2)…a_ij(G)a_ij(G+1)…a_ij(G+P-1)]^T  (13)

步骤九：在接收端信号处理机内，将接收到的信号通过采样样本公式(12)转化为信号空间，再将该信号空间内的信号进行脉间、脉内两次解码处理，即可获得抑制非相参无线电干扰方法处理后的结果；

第一次解码为公式(12)信号与公式(8)产生的脉间伪随机码本序列进行解码处理得：

$\mathrm{RI}=\left[\begin{array}{c}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{A}_{1j}{R}_j\\ \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{A}_{2j}{R}_{J+j}\\ .\\ .\\ .\\ \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{A}_{\mathrm{Nj}}{R}_{(N-1)J+j}\end{array}\right]={\left[\begin{array}{cccc}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}\left(1\right)R_j& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}\left(2\right)R_j& ...& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}(G+P-1)R_j\\ \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}\left(1\right)R_{J+j}& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}\left(2\right)R_{J+j}& ...& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}(G+P-1)R_{J+j}\\ .& .& & .\\ .& .& ...& .\\ .& .& & .\\ \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}\left(1\right)R_{(N-1)J+j}& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}\left(2\right)R_{(N-1)J+j}& ...& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}(G+P-1)R_{(N-1)J+j}\end{array}\right]}_{N\×G}---\left(14\right)$

第二次解码为公式(14)的计算结果与公式(9)产生的脉内码本序列进行解码处理得：

$\mathrm{BRI}={\left[\begin{array}{cccc}\underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}\left(p\right)R_j]& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}(p+1)R_j]& ...& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}(G+p-1)R_j]\\ \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}\left(p\right)R_{J+j}]& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}(p+1)R_{J+j}]& ...& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}(G+p-1)R_{J+j}]\\ .& .& & .\\ .& .& ...& .\\ .& .& & .\\ \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}\left(p\right)R_{(N-1)J+j}]& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}(p+1)R_{(N-1)J+j}]& ...& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}(G+p-1)R_{(N-1)J+j}]\end{array}\right]}_{N\×G}---\left(15\right)$

完成无线电干扰抑制处理。

作为本发明进一步的方案：步骤一中所述选定最优的分数阶次α_opt的方法为：

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

$\begin{array}{cc}{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\π}\end{array}---\left(1\right)$

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

$\begin{array}{cc}\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)\end{array}---\left(2\right)$

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

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

作为本发明再进一步的方案：在步骤二中应用FRFT的分解算法，在最优分数阶α_opt下进行分数阶傅立叶变换，即通过公式

$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)$

得到相应的分数阶域信号形式，其中Δx为信号x(t)的带宽，N＝(Δx)²，函数exp为指数函数，α为分数阶次。

作为本发明再进一步的方案：步骤三使用二阶矩阵法，通过公式(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)$

其中 $\mathrm{\μ}=\frac{1}{N}\·\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\left|F_{\mathrm{\α}}\left(\frac{m}{2\mathrm{\Δx}}\right)\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}$ 为幅度标准偏差的估计值，β为门限优化因子，N为FRFT的点数，α为分数阶次。

与现有技术相比，本发明的有益效果是：风廓线雷达接收信号经过此方法处理后，能够降低无线电信号对风廓线雷达的影响，且能够显著的改善信噪比，实用性强。

附图说明

图1为风廓线雷达干扰抑制方法的流程图。

图2为风廓线雷达干扰抑制方法的原理框图。

具体实施方式

下面结合具体实施方式对本专利的技术方案作进一步详细地说明。

图1为风廓线雷达干扰抑制方法处理流程图。图2为风廓线雷达干扰抑制方法的原理框图。风廓线雷达监控组合控制射频调制器、发射机和信号处理机，风廓线雷达监控组合内的伪随机码产生器输出伪随机码序列，然后送给射频调制器、发射机和信号处理机。接收信号与发射调制脉冲完全同步的脉间伪随机码进行第一次解码处理，其输出结果与脉内伪随机码本进行第二次解码处理，经以上两次解码运算后即可获得非相参无线电波干扰抑制的处理结果，该结果可以作为风廓线雷达后续信号处理的输入信号。

一种风廓线雷达干扰抑制方法，具体步骤为：

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

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

$\begin{array}{cc}{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\π}\end{array}---\left(1\right)$

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

则最优分数阶次α_opt的分数阶傅立叶变换峰值和带宽σ_αopt，u应满足(2)和(3)：

$\begin{array}{cc}\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)\end{array}---\left(2\right)$

$\begin{array}{cc}{\mathrm{\σ}}_{{\mathrm{\α}}_{\mathrm{opt}},u}\≤{\mathrm{\σ}}_{\mathrm{\α},u}& \∀\mathrm{\α}\∈\left(\mathrm{0,2}\mathrm{\π}\right)\end{array}---\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)$

其中，Δ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{\mathrm{j\π}\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个样本的列向量，则分解算法用矩阵形式表示为其中，F^α＝DK_pJ，D和J分别表示内插和抽取运算的矩阵，K_p(m，n)表示如(5)式所示：

$F_{\mathrm{\α}}(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{1}{\mathrm{\Σ}}}\left|F_{\mathrm{\α}}\left(\frac{m}{2\mathrm{\Δx}}\right)\right|---\left(6\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)对抑制后的数据进行逆分数阶傅立叶变换，即求得经过理后的时域数据，至此获得去除杂波后的时域风廓线雷达数据，完成杂波抑制处理。

步骤五：在风廓线雷达发射端，同步产生脉间、脉内两种伪随机码序列；其中脉间伪随机序列采用M序列，序列长度根据实际工作参数选择，假设其码长为L，脉内编码序列采用巴克码，假设其码长为P，

设脉间伪随机编码的码本为R：

R＝[R₁ R₂ … R_i … R_L]^T    (8)

式中R_i为M序列伪随机码码元，i＝1，2...L，T表示矩阵转置；脉内编码码本为巴克码，其位数为P，

BK＝[bk(1) bk(2) …b k(i) … bk(P)]^T     (9)

式中码元为bk(i)，i＝1，2...P，T表示矩阵转置；

步骤六：在风廓线雷达发射端，上述公式(8)和公式(9)码本序列经脉间脉内联合编码序列为：

$\mathrm{RB}=\left[\begin{array}{c}R{B}_1\\ R{B}_2\\ .\\ .\\ .\\ R{B}_i\\ .\\ .\\ .\\ R{B}_L\end{array}\right]---\left(10\right)$

式中RB_i为编码后第i个脉冲内的码本信号，i＝1，2...L，

$R{B}_i=(\mathrm{BK}+R_i)\mathrm{mod}\left(2\right)=\left[\begin{array}{c}\mathrm{bk}\left(1\right)\mathrm{xor}{R}_i\\ \mathrm{bk}\left(2\right)\mathrm{xor}{R}_i\\ .\\ .\\ .\\ \mathrm{bk}\left(P\right)\mathrm{xor}{R}_i\end{array}\right]---\left(11\right)$

式中mod()为求余运算，xor为异或运算符；

步骤七：将公式(10)产生的码本序列传送至风廓线雷达发射信号调制器并发射，同时将公式(8)和公式(9)产生的码本序列保存在接收信号处理机中；

步骤八：在风廓线雷达接收端对步骤七中的发射信号进行接收；

设时域积累为J次，FFT点数为N，距离门数为G，则信号处理机需处理的信号样本空间I为：

$I={\left[\begin{array}{cccc}{A}_{11}& A_{12}& ...& A_{1,J}\\ A_{21}& A_{22}& ...& A_{2J}\\ .& .& & .\\ .& .& A_{\mathrm{ij}}& .\\ .& .& & .\\ A_{N1}& A_{N2}& ...& A_{\mathrm{NJ}}\end{array}\right]}_{N\×J}---\left(12\right)$

式中A_ij分别代表第i×j次采样信号a_ij(1)，a_ij(2)，...a_ij(G+P-1)构成的样本空间，其中i＝1，2...N，j＝1，2...J分别与信号处理中的FFT点数和相干积累次数有关，L＝N×J；

A_ij＝[a_ij(1)a_ij(2)…a_ij(G)a_ij(G+1)…a_ij(G+P-1)]^T  (13)

步骤九：在接收端信号处理机内，将接收到的信号通过采样样本公式(12)转化为信号空间，再将该信号空间内的信号进行脉间、脉内两次解码处理，即可获得抑制非相参无线电干扰方法处理后的结果；

第一次解码为公式(12)信号与公式(8)产生的脉间伪随机码本序列进行解码处理得：

$\mathrm{RI}=\left[\begin{array}{c}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{A}_{1j}{R}_j\\ \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{A}_{2j}{R}_{J+j}\\ .\\ .\\ .\\ \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{A}_{\mathrm{Nj}}{R}_{(N-1)J+j}\end{array}\right]={\left[\begin{array}{cccc}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}\left(1\right)R_j& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}\left(2\right)R_j& ...& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}(G+P-1)R_j\\ \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}\left(1\right)R_{J+j}& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}\left(2\right)R_{J+j}& ...& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}(G+P-1)R_{J+j}\\ .& .& & .\\ .& .& ...& .\\ .& .& & .\\ \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}\left(1\right)R_{(N-1)J+j}& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}\left(2\right)R_{(N-1)J+j}& ...& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}(G+P-1)R_{(N-1)J+j}\end{array}\right]}_{N\×G}---\left(14\right)$

第二次解码为公式(14)的计算结果与公式(9)产生的脉内码本序列进行解码处理得：

$\mathrm{BRI}={\left[\begin{array}{cccc}\underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}\left(p\right)R_j]& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}(p+1)R_j]& ...& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}(G+p-1)R_j]\\ \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}\left(p\right)R_{J+j}]& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}(p+1)R_{J+j}]& ...& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}(G+p-1)R_{J+j}]\\ .& .& & .\\ .& .& ...& .\\ .& .& & .\\ \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}\left(p\right)R_{(N-1)J+j}]& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}(p+1)R_{(N-1)J+j}]& ...& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}(G+p-1)R_{(N-1)J+j}]\end{array}\right]}_{N\×G}---\left(15\right)$

至此，该风廓线雷达无线电干扰抑制处理完毕。

举例说明，设大气信号频率为10Hz，无线电干扰为30Hz，两者功率相同，脉间伪随机码选码长为255的8级结构M序列，级间反馈连接为8、6、5、4，脉内编码选择5位巴克码，则经过全相参脉间脉内联合伪随机编码解码处理后的无限电干扰抑制约20dB，风信号功率提高约13dB，改善信噪比达20dB以上。

上面对本专利的较佳实施方式作了详细说明，但是本专利并不限于上述实施方式，在本领域的普通技术人员所具备的知识范围内，还可以在不脱离本专利宗旨的前提下作出各种变化。

Claims (4)

1.一种风廓线雷达干扰抑制方法，其特征在于，包括以下步骤：

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

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

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

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

步骤五：在风廓线雷达发射端，同步产生脉间、脉内两种伪随机码序列，其中脉间伪随机序列采用M序列，序列长度根据实际工作参数选择，假设其码长为L，脉内编码序列采用巴克码，假设其码长为P，设脉间伪随机编码的码本为R：

R＝[R₁ R₂ … R_i … R_L]^T            (8)

式中R_i为M序列伪随机码码元，i＝1，2...L，T表示矩阵转置；

脉内编码码本为巴克码，其位数为P，

BK＝[bk(1) bk(2) … bk(i) … bk(P)]^T           (9)

式中码元为bk(i)，i＝1，2...P，T表示矩阵转置；

步骤六：在风廓线雷达发射端，上述公式(8)和公式(9)码本序列经脉间脉内联合编码序列为：

$\mathrm{RB}=\left[\begin{array}{c}{\mathrm{RB}}_1\\ {\mathrm{RB}}_2\\ .\\ .\\ .\\ {\mathrm{RB}}_i\\ .\\ .\\ .\\ {\mathrm{RB}}_L\end{array}\right]---\left(10\right)$

式中RB_i为编码后第i个脉冲内的码本信号，i＝1，2...L，

${\mathrm{RB}}_i=(\mathrm{BK}+R_i)\mathrm{mod}\left(2\right)=\left[\begin{array}{c}\mathrm{bk}\left(1\right)\mathrm{xor}{R}_i\\ \mathrm{bk}\left(2\right)\mathrm{xor}{R}_i\\ .\\ .\\ .\\ \mathrm{bk}\left(P\right)\mathrm{xor}{R}_i\end{array}\right]---\left(11\right)$

式中mod()为求余运算，xor为异或运算符；

步骤七：将公式(10)产生的码本序列传送至风廓线雷达发射信号调制器并发射，同时将公式(8)和公式(9)产生的码本序列保存在接收信号处理机中；

步骤八：在风廓线雷达接收端对步骤七中的发射信号进行接收，设时域积累为J次，FFT点数为N，距离门数为G，则信号处理机需处理的信号样本空间I为：

$I={\left[\begin{array}{cccc}{A}_{11}& A_{12}& ...& A_{1J}\\ A_{21}& A_{22}& ...& A_{2J}\\ .& .& & .\\ .& .& A_{\mathrm{ij}}& .\\ .& .& & .\\ A_{N1}& A_{N2}& ...& A_{\mathrm{NJ}}\end{array}\right]}_{N\×J}---\left(12\right)$

式中A_ij分别代表第i×j次采样信号a_ij(1)，a_ij(2)，...a_ij(G+P-1)构成的样本空间，其中i＝1，2...N，j＝1，2...J分别与信号处理中的FFT点数和相干积累次数有关，L＝N×J；

A_ij＝[a_ij(1)a_ij(2)…a_ij(G)a_ij(G+1)…a_ij(G+P-1)]^T   (13)

步骤九：在接收端信号处理机内，将接收到的信号通过采样样本公式(12)转化为信号空间，再将该信号空间内的信号进行脉间、脉内两次解码处理，即可获得抑制非相参无线电干扰方法处理后的结果；

第一次解码为公式(12)信号与公式(8)产生的脉间伪随机码本序列进行解码处理得：

$\mathrm{RI}=\left[\begin{array}{c}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{A}_{1j}{R}_j\\ \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{A}_{2j}{R}_{J+j}\\ .\\ .\\ .\\ \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{A}_{\mathrm{Nj}}{R}_{(N-1)J+j}\end{array}\right]={\left[\begin{array}{cccc}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}\left(1\right)R_j& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}\left(2\right)R_j& ...& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}(G+P-1)R_j\\ \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}\left(1\right)R_{J+j}& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}\left(2\right)R_{J+j}& ...& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}(G+P-1)R_{J+j}\\ .& .& & .\\ .& .& ...& .\\ .& .& & .\\ \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}\left(1\right)R_{(N-1)J+j}& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}\left(2\right)R_{(N-1)J+j}& ...& \underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}(G+P-1)R_{(N-1)J+j}\end{array}\right]}_{N\×G}---\left(14\right)$

第二次解码为公式(14)的计算结果与公式(9)产生的脉内码本序列进行解码处理得：

$\mathrm{BRI}={\left[\begin{array}{cccc}\underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}\left(p\right)R_j]& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}(p+1)R_j]& ...& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{1j}(G+p-1)R_j]\\ \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}\left(p\right)R_{J+j}]& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}(p+1)R_{J+j}]& ...& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{2j}(G+p-1)R_{J+j}]\\ .& .& & .\\ .& .& ...& .\\ .& .& & .\\ \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}\left(p\right)R_{(N-1)J+j}]& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}(p+1)R_{(N-1)J+j}]& ...& \underset{p=1}{\overset{P}{\mathrm{\Σ}}}[\mathrm{bk}\left(p\right)\·\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{a}_{\mathrm{Nj}}(G+p-1)R_{(N-1)J+j}]\end{array}\right]}_{N\×G}---\left(15\right)$

完成无线电干扰抑制处理。

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

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

$\begin{array}{cc}{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\π}\end{array}---\left(1\right)$

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

$\begin{array}{cc}\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)\end{array}---\left(2\right)$

$\begin{array}{cc}{\mathrm{\σ}}_{{\mathrm{\α}}_{\mathrm{opt}},u}\≤{\mathrm{\σ}}_{\mathrm{\α},u}& \∀\mathrm{\α}\∈\left(\mathrm{0,2}\mathrm{\π}\right)\end{array}---\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}}\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)$

得到相应的分数阶域信号形式，其中Δx为信号x(t)的带宽，N＝(Δx)²，函数exp为指数函数，α为分数阶次。

4.根据权利要求1或2或3所述的风廓线雷达干扰抑制方法，其特征在于，步骤三使用二阶矩阵法，通过公式(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)$

其中为幅度均值的估计值， $\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}$ 为幅度标准偏差的估计值，β为门限优化因子，N为FRFT的点数，α为分数阶次。