CN102393525B - 子空间投影的导航干扰抑制与信号增强方法 - Google Patents

Abstract

本发明公开了一种子空间投影的导航干扰抑制与信号增强方法，主要解决现有方法抗干扰输出性能不佳甚至失效，且无法获得空间增益的问题。其实现过程为：(1)对阵列接收数据估计协方差矩阵，进行特征分解；(2)估计干扰个数，利用特征向量构造干扰正交补空间；(3)将阵列接收数据向干扰正交补空间投影；(4)设置卫星序号p为循环变量；(5)将抗干扰后数据频移，与本地第p颗卫星信号做相关运算，记录相关峰及相关峰主副比；(6)对相关峰主副比大于门限的卫星，用相关峰矢量作波束形成，增强卫星p的强度。本方法能稳健抑制干扰，对卫星信号增强效果显著，输出信号可获得接近理论值的空间增益，可用于导航卫星的干扰抑制及信号增强。

Description

子空间投影的导航干扰抑制与信号增强方法

技术领域

本发明属于信号处理技术领域，涉及阵列信号处理技术，可以用于导航接收机强干扰抑制和弱信号增强，提高导航接收机输出性能，实测数据验证了方法的有效性。

背景技术

卫星导航定位系统可全天候、全天时提供位置、速度以及时间等信息，在民用和国防领域有着广阔的应用需求。但是，由于电磁环境日益恶劣并且到达接收机的导航信号非常微弱，使得导航接收机对干扰极其敏感，这种对干扰敏感的特性严重影响导航接收机的输出性能，因此在干扰背景下提高导航接收机的输出性能是十分重要的。

传统的抗干扰算法如采样协方差矩阵求逆SMI（Sample Matrix Inverse）方法、最小均方误差LMS（Least Mean Square）迭代方法在理想条件下具有良好的抗干扰性能，但在实际中由于存在系统误差，如阵列方向图不一致、多通道幅相误差、指向误差以及训练样本数少等，导致常规抗干扰方法输出性能不佳甚至失效。还有一种比较通用的抗干扰算法是功率倒置PI（Power Inverse）算法，干扰功率越强，抑制效果越好，因此可以有效的抑制强干扰，但是该算法对干扰抑制的同时，对导航信号没有增强效果，甚至有可能将其削弱，无法获得信号的空间增益。还有采用自适应调零天线进行抗干扰的技术，该方法也能通过调节阵列指向，在干扰处形成零陷实现干扰抑制，但也存在着无法获得信号空间增益的问题。卢丹、吴仁彪等人在文献基于自相干MUSIC算法的全球定位系统干扰抑制方法中提出了一种基于导航干扰抑制技术，该技术在均匀线阵下进行投影抗干扰后，采用MUSIC谱峰搜索进行卫星信号方向估计，而实际场景中导航信号是完全淹没在噪声中的，进行MUSIC进行谱峰搜索具有很大难度，从而很难再通过波束形成获得导航信号空间增益。

通过阵列优化布阵结合上述抗干扰算法的研究也是抗干扰技术的一个方面，但实际系统中考虑到系统复杂度，一般阵元数有限的情况下阵列优化布阵也就受到限制，主要还是通过传统抗干扰技术实现干扰抑制，而这些干扰技术往往无法获得导航信号空间增益。

发明内容

本发明的目的在于克服上述已有技术的不足，提出了一种子空间投影的导航干扰抑制与信号增强方法，以在有系统误差和强干扰的情况下，实现干扰稳健抑制和导航卫星信号空间增益的获得，以便于后续导航电文的准确提取。

为实现上述目的，本发明方法包括如下步骤：

（1）对阵列接收数据x(t)进行子空间投影，得到抑制干扰后的数据x₁(t)；

（2）利用导航信号自相关特性，寻找x₁(t)中视野内卫星的最优自相关矢量，利用该相关矢量对数据x₁(t)进行波束形成，获得该卫星信号空间增益：

（2.1）设置外循环变量p=1，p表示卫星序号；

（2.2）产生本地第p颗卫星2ms的C/A码序列作为本地参考信号

并设多普勒频率为：

f_d(q)＝[-5+(q-1)×0.5]KHz，

其中上标ref表示该信号为本地参考信号，q为内循环变量，表示卫星信号多普勒频率序号，初始化q=1，多普勒频率f_d从-5KHz变化到5KHz，以0.5KHz为步长；

（2.3）将步骤（1）中产生的数据x₁(t)做频移后与步骤（2.2）产生的做相关运算，记录相关结果中最大峰值矢量和次大峰值矢量至矩阵

的第q列和矩阵

的第q列：

$Z_1^p(:,q)=\mathrm{peak}\left\{\right[x_1\left(t\right)e^{j2\mathrm{\π}(f_d\left(q\right)+f_0)t}]*\mathrm{conj}[s_p^{\mathrm{ref}}\left(t\right)],1\},$

$Z_2^p(:,q)=\mathrm{peak}\left\{\right[x_1\left(t\right)e^{j2\mathrm{\π}(f_d\left(q\right)+f_0)t}]*\mathrm{conj}[s_p^{\mathrm{ref}}\left(t\right)],2\},$

其中：

表示将x₁(t)各行数据乘以f_d(q)加上中频f₀＝7.5MHz的复单音信号实现频移，*表示相关运算，conj表示取共轭，peak{·,k}表示取第k高峰值矢量，和

分别表示最大峰值矩阵

的第q列和次大峰值矩阵

的第q列；

（2.4）更新内循环变量q＝q+1，当q≤21时，返回步骤（2.3），否则执行步骤（2.5）；

（2.5）在步骤（2.3）中得到的矩阵

的第1行中寻找最大值的列标q₀，计算卫星p在多普勒频率点q₀处的相关峰主副比矢量r_p：

$r_p=Z_1^p(:,q_0)./Z_2^p(:,q_0),$

其中./表示对应元素相除，

表示取矩阵

的第1行中最大值的列标；

（2.6）根据步骤（2.5）得到的相关峰主副比矢量r_p判断是否需要对第p颗卫星进行信号增强，如果r_p中所有元素都大于设定门限T_h说明视野中存在卫星p，需要对其增强，执行步骤（2.7），否则跳至步骤（2.8）；

（2.7）由步骤（2.5）中的列标q₀得到第p颗卫星的归一化最佳相关峰矢量b_p，用b_p对步骤（1）中得到的干扰抑制后数据x₁(t)进行波束形成，获得第p颗卫星的空间增益：

$y_p\left(t\right)=b_p^H{x}_1\left(t\right),$

其中：

它表示第p颗卫星的归一化最佳相关峰矢量，表示通道1在q₀频率点处的相关峰值，上标H表示对矢量取共轭转置；

（2.8）更新循环变量p＝p+1，如果p≤32，返回步骤（2.2），否则结束。

本发明与现有技术相比，具有以下优点：

（1）由于本发明的干扰子空间完全从数据中估计得到，抗干扰性能对系统误差具有稳健性；

（2）本发明通过利用视野内卫星的归一化最优相关峰矢量对抗干扰后数据进行波束形成，能够提高卫星的输出信噪比，获得接近理论值的空间增益10×log₁₀M；

（3）本发明由于对接收数据的协方差矩阵进行特征值分解，构造干扰的补空间，将接收数据直接向该干扰的补空间投影，故无需估计干扰或目标的来波方向，保证了抗干扰的有效性，简化了运算。

实测数据处理结果表明：与常规矩阵求逆SMI抗干扰和最小均方误差LMS迭代方法相比，本发明通过使用子空间投影抗干扰与信号增强技术，有效抑制干扰，使卫星信号相关峰的主副比MSR得到显著提高，从而获得了导航信号的空间增益。

附图说明

图1是本发明流程图

图2是本发明采集数据的天线阵列排布图；

图3是本发明中接收数据x(t)中存在的点频干扰频谱；

图4是用常规方法与本发明方法抑制点频干扰后的频谱图；

图5是用常规方法与本发明方法抑制点频干扰后信号相关峰值；

图6是本发明中接收数据x(t)中存在的宽带干扰频谱；

图7是用常规方法与本发明方法抑制宽带干扰后的频谱图；

图8是用常规方法与本发明方法抑制宽带干扰后信号相关峰值。

具体实施方式

参照图1，本发明的实现具体步骤如下：

步骤1．对阵列接收数据进行子空间投影，得到抑制干扰后的数据。

（1.1）使用如图2所示系统，利用M个阵元组成的阵列天线，此处M＝4，本发明不局限于4阵元情况，可适用于任意阵元数组成的阵列天线，得到t时刻一批包含L颗卫星信号和J个干扰的阵列接收数据x(t)：

其中

表示第l颗卫星的导向矢量，

分别表示第l颗卫星的方位角和俯仰角，s_l(t)表示t时刻第l颗卫星信号，此处l为卫星编号，

表示第k个干扰的导向矢量，

分别表示第k个干扰的方位角和俯仰角，j_k(t)表示t时刻的第k个干扰，此处k为干扰编号，n(t)表示t时刻接收数据中的噪声成分；

（1.2）对步骤（1.1）中的阵列接收数据x(t)估计其M×M维协方差矩阵

对

做特征分解，得到

的特征值与特征值对应的特征向量：

$[U,\mathrm{\Λ}]=\mathrm{eig}\left({\hat{R}}_x\right),$

其中

它表示协方差矩阵的估计，N表示采样快拍数，H表示取共轭转置，eig(·)表示特征分解，Λ＝diag([λ₁ λ₂...λ_M])为特征值对角矩阵，diag(·)表示将矢量对角化，不失一般性特征值由大到小排列有λ₁≥λ₂≥...≥λ_M，U＝[u₁,u₂,…,u_M]为特征矢量矩阵，u_i为第i个特征矢量，i＝1,2,…,M，M为阵元数；

（1.3）由M.Wax和T.Kailath在文献Detection of Signal by Information TheoretiedCriteria，IEEE，Trans.ASSP-33No.2，1985中提出的信息论AIC信源数估计准则，得到干扰个数Ψ的估计值：

$\mathrm{\Ψ}=\arg \underset{k}{\min }\left\{\mathrm{AIC}\left(k\right)\right\}=\arg \underset{k}{\min }\{\mathrm{LLF}\left(k\right)+P_{\mathrm{AIC}}(N,M,k)\}$

其中 $\mathrm{LLF}\left(k\right)=N(M-k)\log (\frac{1}{M-k}\underset{i=k+1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\λ}}_i/{\left(\underset{i=k+1}{\overset{M}{\mathrm{\Π}}}{\mathrm{\λ}}_i\right)}^{\frac{1}{M-k}}),$ 表示k的对数似然函数，P_AIC(N,M,k)＝k(2M-k)，表示惩罚函数，用于对LLF(k)进行约束，λ_i为协方差矩阵

的特征值，M为阵元数，N为采样快拍数，k＝1,2,…，M-1，表示对k寻优，AIC(k)最小时对应的k为干扰个数Ψ的估计值；

（1.4）根据步骤（1.3）估计的干扰个数Ψ，得到干扰正交补空间B：

$B=I-U_J{\left(U_J^H{U}_J\right)}^{-1}{U}_J^H,$

其中U_J＝[u₁,u₂，...u_Ψ]，它是干扰特征矢量组成的矩阵，Ψ为干扰个数的估计值，H表示取共轭转置，I为M维单位阵，M是阵元数；

（1.5）将步骤（1.1）中的阵列接收数据x(t)向步骤（1.4）得到的干扰正交补空间B投影，投影后数据记为x₁(t)：

x₁(t)＝Bx(t)。

步骤2．利用导航卫星信号的自相关特性，寻找数据x₁(t)中卫星的最优自相关矢量，利用该自相关矢量对数据x₁(t)进行波束形成，获得卫星信号空间增益。

（2.1）设置外循环变量p=1，p表示卫星序号，本发明以全球定位系统GPS导航卫星为例，共对32颗卫星进行处理；

（2.2）产生本地第p颗卫星2ms的GPS卫星L1载波对应的民用伪随机码，即C/A码序列作为本地参考信号

并设多普勒频率为：

f_d(q)＝[-5+(q-1)×0.5]KHz，

其中

上标ref表示该信号为本地参考信号，q为内循环变量，表示第p颗卫星信号多普勒频率序号，初始化内循环变量q=1，多普勒频率f_d从-5KHz变化到5KHz，以0.5KHz为步长，共21个频率点，循环变量q≤21；

（2.3）将数据x₁(t)做频移后与步骤（2.2）产生的本地参考信号

做相关运算，记录相关结果中最大峰值矢量和次大峰值矢量至最大峰值矩阵的第q列和次大峰值矩阵

的第q列：

$Z_1^p(:,q)=\mathrm{peak}\left\{\right[x_1\left(t\right)e^{j2\mathrm{\π}(f_d\left(q\right)+f_0)t}]*\mathrm{conj}[s_p^{\mathrm{ref}}\left(t\right)],1\},$

$Z_2^p(:,q)=\mathrm{peak}\left\{\right[x_1\left(t\right)e^{j2\mathrm{\π}(f_d\left(q\right)+f_0)t}]*\mathrm{conj}[s_p^{\mathrm{ref}}\left(t\right)],2\},$

其中

表示将x₁(t)各行数据乘以f_d(q)加上中频f₀＝7.5MHz的复单音信号实现频移，*表示做相关运算，conj表示取共轭，peak{·,k}表示取第k高峰值矢量，

和

分别表示最大峰值矩阵

的第q列和次大峰值矩阵的第q列，该相关运算是指求两个随机变量之间的协方差，用来表征两个信号之间的相似程度；

（2.4）更新内循环变量q＝q+1，当q≤21时，返回步骤（2.3），否则执行步骤（2.5）；

（2.5）在步骤（2.3）中得到的最大峰值矩阵

的第1行中寻找最大值的列标q₀，计算卫星p在多普勒频率点q₀处的相关峰主副比矢量r_p：

$r_p=Z_1^p(:,q_0)./Z_2^p(:,q_0),$

其中./表示对应元素相除，

表示取矩阵

的第1行中最大值的列标，表示最大峰值矩阵

的第1行；

（2.6）根据步骤（2.5）得到的相关峰主副比矢量r_p判断是否需要对第p颗卫星进行信号增强，如果r_p中所有元素都大于设定门限T_h说明视野中存在卫星p，需要对其增强，执行步骤（2.7），否则跳至步骤（2.8）；

（2.7）由步骤（2.5）中的列标q₀得到第p颗卫星的归一化最佳相关峰矢量b_p，用b_p对步骤（1）中得到的干扰抑制后数据x₁(t)进行波束形成，获得第p颗卫星的空间增益：

$y_p\left(t\right)=b_p^H{x}_1\left(t\right),$

其中：

它表示第p颗卫星的归一化最佳相关峰矢量，

表示通道1在q₀频率点处的相关峰值，上标H表示对矢量取共轭转置；

（2.8）更新循环变量p＝p+1，如果p≤32，返回步骤（2.2），否则结束。

本发明的效果可以通过以下实测数据处理结果进一步说明。

1．实验环境

本发明采集数据的天线阵列排布图如图2所示，图中4个阵元在xoy平面呈均匀圆阵排列，阵元间距0.095m，该间距是GPS卫星L1载频1.57542GHz对应的半波长；GPS卫星信号经阵列天线到低噪放和下变频器后，获得55dB增益；图2所示采集数据的天线阵列，它的下变频中心频率为7.5MHz，数据采集部分频率为10MHz。使用图2天线阵列在外场录取数据，在该阵列接收数据中分别加点频干扰和宽带干扰对本发明方法进行验证，接收数据如图3和图6所示。

2．实验内容与结果

实验1，在阵列接收数据中加1个频率为1.57542GHz的-50dBm的点频干扰。使用常规方法与本发明方法，对图3所示接收数据进行抗干扰处理，抗干扰效果如图4所示。其中图4（a）是用常规LMS方法进行抗干扰后的效果图，图4（b）是用常规SMI方法抗干扰效果图，图4（c）是本发明方法抗干扰效果图。使用本发明对抗干扰后数据进行信号增强，与常规方法相比，输出信号相关峰值情况如图5所示。

由图4（a）可以看出，常规LSM方法抑制图3中的干扰时在干扰频率点处形成凹陷，卫星信号有损失。

由图4（b）可以看出，用常规SMI方法可有效抑制图3中的点频干扰。

由图4（c）可以看出，本发明方法可以有效抑制图3中的点频干扰。

对图5所示信号相关峰进行统计，统计结果如表1所示。由表1可以看出，与常规方法相比，本发明提高了卫星信号相关峰，获得了信号的空间增益。

表1 点频干扰下各方法处理后卫星信号相关峰统计结果（单位dB）

实验2，阵列接收数据中加1个带宽为2MHz，中心频率是1.57542GHz的-40dBm宽带干扰。使用常规方法与本发明方法，对图6所示接收数据进行抗干扰处理，抗干扰效果如图7所示。其中图7（a）是用常规SMI方法抗干扰效果图，图7（b）是本发明方法抗干扰效果图。使用本发明对抗干扰后数据进行信号增强，与常规方法相比，输出信号相关峰值情况如图8所示。

由图7（a）可以看出，常规SMI方法可以有效抑制图6中的宽带干扰。

由图7（b）可以看出，本发明也可有效抑制图6中的宽带干扰。

对图8所示信号相关峰进行统计，统计结果如表2所示。由表2可以看出，与常规方法相比，本发明提高了卫星信号相关峰，获得了信号的空间增益。

表2 宽带干扰下各方法处理后卫星信号相关峰统计结果（单位dB）

综上，由实验1和实验2结果可知，本发明方法能够有效的抑制针对导航信号的点频干扰和宽带干扰，并有效提高导航信号强度，获得导航信号的空间增益。

Claims (2)

1.一种导航信号子空间干扰抑制与信号增强方法，包括如下步骤：

(1)对阵列接收数据x(t)进行子空间投影，得到抑制干扰后的数据x₁(t)；

(2)利用导航信号自相关特性，寻找x₁(t)中视野内卫星的最优自相关矢量，利用该相关矢量对数据x₁(t)进行波束形成，获得该卫星信号空间增益：

(2.1)设置外循环变量p＝1，p表示卫星序号；

(2.2)产生本地第p颗卫星2ms的C/A码序列作为本地参考信号

并设多普勒频率为：

f_d(q)＝[-5+(q-1)×0.5]KHz，

其中

上标ref表示该信号为本地参考信号，q为内循环变量，表示卫星信号多普勒频率序号，初始化q＝1，多普勒频率f_d从-5KHz变化到5KHz，以0.5KHz为步长；

(2.3)将步骤(1)中产生的数据x₁(t)做频移后与步骤(2.2)产生的

做相关运算，记录相关结果中最大峰值矢量和次大峰值矢量至矩阵

的第q列和矩阵的第q列：

$Z_1^p(:,q)=\mathrm{peak}\left\{\right[x_1\left(t\right)e^{j2\mathrm{\π}(f_d\left(q\right)+f_0)t}]*\mathrm{conj}[s_p^{\mathrm{ref}}\left(t\right)],1\},$

$Z_2^p(:,q)=\mathrm{peak}\left\{\right[x_1\left(t\right)e^{j2\mathrm{\π}(f_d\left(q\right)+f_0)t}]*\mathrm{conj}[s_p^{\mathrm{ref}}\left(t\right)],2\},$

其中：

表示将x₁(t)各行数据乘以f_d(q)加上中频f₀＝7.5MHz的复单音信号实现频移，*表示相关运算，conj表示取共轭，peak{·，k}表示取第k高峰值矢量，

和分别表示最大峰值矩阵

的第q列和次大峰值矩阵

的第q列；

(2.4)更新内循环变量q＝q+1，当q≤21时，返回步骤(2.3)，否则执行步骤(2.5)；

(2.5)在步骤(2.3)中得到的矩阵

的第1行中寻找最大值的列标q₀，计算卫星p在多普勒频率点q₀处的相关峰主副比矢量r_p：

$r_p=Z_1^p(:,q_0)./Z_2^p(:,q_0),$

其中./表示对应元素相除，

表示取矩阵

的第1行中最大值的列标；

(2.6)根据步骤(2.5)得到的相关峰主副比矢量r_p判断是否需要对第p颗卫星进行信号增强，如果r_p中所有元素都大于设定门限T_h说明视野中存在卫星p，需要对其增强，执行步骤(2.7)，否则跳至步骤(2.8)；

(2.7)由步骤(2.5)中的列标q₀得到第p颗卫星的归一化最佳相关峰矢量b_p，用b_p对步骤(1)中得到的干扰抑制后数据x₁(t)进行波束形成，获得第p颗卫星的空间增益：

$y_p\left(t\right)=b_p^H{x}_1\left(t\right),$

其中：

它表示第p颗卫星的归一化最佳相关峰矢量，

表示通道1在q₀频率点处的相关峰值，上标H表示对矢量取共轭转置；

(2.8)更新循环变量p＝p+1，如果p≤32，返回步骤(2.2)，否则结束。

2.根据权利要求1所述的方法，其特征在于步骤(1)所述的对阵列接收数据x(t)进行子空间投影，按如下步骤进行：

(1.1)利用M个阵元组成的阵列天线接收导航数据，得到t时刻一批包含L颗卫星信号和J个干扰的阵列接收数据x(t)：

其中

表示第l颗卫星的导向矢量，

分别表示第l颗卫星的方位角和俯仰角，s_l(t)表示t时刻第l颗卫星信号，此处l为卫星编号，

表示第k个干扰的导向矢量，

分别表示第k个干扰的方位角和俯仰角，j_k(t)表示t时刻的第k个干扰，此处k为干扰编号，n(t)表示t时刻接收数据中的噪声成分；

(1.2)对步骤(1.1)中的阵列接收数据x(t)估计其M×M维协方差矩阵对

做特征分解，得到其特征值与相应的特征向量：

$[U,\mathrm{\Λ}]=\mathrm{eig}\left({\hat{R}}_x\right),$

其中

它表示协方差矩阵的估计，N表示采样快拍数，H表示取共轭转置，eig(·)表示特征分解，A＝diag([λ₁ λ₂ ...λ_M])为特征值对角矩阵，diag(·)表示将矢量对角化，不失一般性特征值由大到小排列有λ₁≥λ₂≥...≥λ_M，U＝[u₁，u₂，…，u_M]为特征矢量矩阵，u_i为第i个特征矢量，i＝1，2，…，M，M为阵元数；

(1.3)由信息论AIC信源数估计准则，得到干扰个数Ψ的估计值：

$\mathrm{\Ψ}=\arg \underset{k}{\min }\left\{\mathrm{AIC}\left(k\right)\right\}=\arg \underset{k}{\min }\{\mathrm{LLF}\left(k\right)+P_{\mathrm{AIC}}(N,M,k)\}$

其中 $\mathrm{LLF}\left(k\right)=N(M-k)\log (\frac{1}{M-k}\underset{i=k+1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\λ}}_i/{\left(\underset{i=k+1}{\overset{M}{\mathrm{\Π}}}{\mathrm{\λ}}_i\right)}^{\frac{1}{M-k}}),$ 表示对数似然函数，P_AIC(N，M，k)＝k(2M-k)，表示惩罚函数，λ_i为协方差矩阵

的特征值，M为阵元数，N为采样快拍数，k＝1，2，...，M-1，

表示对k寻优，AIC(k)最小时对应的k为Ψ的估计值；

(1.4)根据步骤(1.3)估计的干扰个数Ψ，得到干扰正交补空间B，阵列接收数据x(t)向干扰正交补空间B投影，投影后数据记为x₁(t)：

x₁(t)＝Bx(t)，

其中

它表示干扰正交补空间，U_J＝[u₁，u₂，...u_Ψ]为干扰特征矢量组成的矩阵，H表示取共轭转置，I为M维单位阵。

Images