CN103792550B - 一种基于阵列天线和gps/sins的联合抗干扰方法 - Google Patents

Abstract

本发明提供的是一种基于阵列天线和GPS/SINS的联合抗干扰方法。初始化载体的位置和姿态之后，建立GPS/SINS组合导航状态方程和量测方程；GPS/SINS组合导航实时提供载体的位置和姿态，根据卫星星历信息计算出当前卫星的位置，获得卫星到载体之间的导向矢量；所述导向矢量作为多约束最小方差空时自适应处理算法的先验信息，在空域、时域同时抑制宽带干扰和窄带干扰。本发明能够同时在多颗可视卫星方向形成波束，在干扰方向形成零陷，从而增强卫星信号的同时抑制干扰信号。本发明采用圆形结构天线阵列，GPS/SINS组合导航为阵列天线的波束形成提供载体的位置和姿态，采用卫星星历提供卫星的位置，从而为波束形成提供先验信息。

Description

一种基于阵列天线和GPS/SINS的联合抗干扰方法

技术领域

本发明涉及的是一种组合导航技术，具体地说是一种阵列天线和GPS（全球定位***）/SINS（捷联惯性导航***）组合导航技术。

背景技术

GPS导航***能够全天候的为其用户提供覆盖全球的精确位置、速度与时间信息，其应用价值越来越高。然而，随着人为干扰技术的提高，卫星只依靠其扩频体制进行抗干扰已经不能满足用户需求。根据ICD-200，商用GPS接收机的抗干扰容限不超过24dB（取决于噪声电平）。即如果干信比大于24dB，商用GPSC/A码接收机就无法保持对信号的跟踪。试验表明，功率为1W的干扰机可以使85公里以内的C/A码接收机无法工作。另外，现有GPS信号是在众所周知的频率上发射的，其调制特征广为人知，信噪比又比较低，因而很容易进行干扰或欺骗。因此，研究GPS抗干扰技术已成为热点和重点。

目前，基于信号处理的抗干扰方法是研究最为活跃的领域，常见的技术包括时域滤波技术和空域滤波技术。自适应时域滤波技术是一种窄带干扰抑制方法，它将接收到的有用信号、干扰和噪声通过自适应算法使所需的代价函数最小化来移除干扰信号。时域滤波在用于有界干扰源时是十分稳定的，这是因为它能同时提供复杂带阻滤波器准则，并且它被看作是GPS接收机前后处理的嵌入部分或是GPS接收机前的独立嵌入部分。在对尺寸影响很小的情况下，这种技术对窄带干扰的抑制大于30dB。对于消除窄带干扰，时域滤波技术可用于复杂窄带和连续波干扰源，但这又会受到残余计算带宽的限制，这种残余计算带宽会妨碍有效的GPS信号处理。空域滤波技术是通过自适应天线阵来实现，它可以有效抑制相干干扰和宽带干扰。自适应天线与时域自适应滤波有点类似，也要使得某个代价函数最小化，它有个很大的缺陷就是如果有用信号和干扰信号的入射方向彼此靠近时，自适应天线在抑制干扰信号的同时对有用信号也会造成影响。它的实现需要消耗比较大的功率、更多的花费，而且需要一个比较大的工作平台。空域滤波主要有两种类型:零陷和波束形成。零陷技术对有用信号信息量的需求是最小的，它比波束形成技术更容易实现。零陷技术最小化信号的输出功率，它可以分为窄带和宽带两种情况。波束形成技术要求比较多的有关有用信号的信息，而且实现起来也比较复杂。波束形成技术最大化信号的输出信噪比，波束形成技术也分为窄带和宽带两种情况。

单纯的时域滤波和单纯的空域各自具有优缺点，但这两种处理方法的优劣恰好可以互补，可以将二者结合应用，进而形成一种联合抗干扰技术，即时空二维联合处理（STAP）抗干扰技术。空时二维联合处理的空域处理能力比单纯的空域滤波要更强。当需要对宽带信号进行零陷的时候，由于信号带宽不可忽略，会在空间中造成所谓“发散”现象。一般的空域滤波手段都需要对这种宽带引发的空间“发散”予以特别的考虑，这也是单纯的空域方法处理宽带信号难以获得满意效果的原因。空时二维联合处理将频率信息和空域信息结合使用，自然地会在空时平面上将发散的各信号分量完整的表现出来，从而可以最大限度地对干扰进行消除和抑制。并且，STAP技术还具有内在的波束形成、内在的信号均衡和内在的抗多径干扰能力等诸多优点，因而可以在增强信号的同时实现抗干扰处理。

另外，GPS/SINS组合导航也能够提高GPS接收机的抗干扰能力。接收机的环路带宽需要在抗干扰能力和动态跟踪性能之间折中设计，即接收机为了提高自身抑制干扰和噪声的能力，需要将环路带宽减小，而为了跟踪载体的动态性能又需要将环路带宽增大。而引入惯性信息后，SINS可以精确估计出载体的速度，计算出载体相对于卫星的多普勒频移，从而减小接收机的环路带宽，增加接收机的抗干扰能力。

然而，以上传统的抗干扰方法都存在一定的缺陷。单纯的时域处理方法不能够抑制宽带干扰，空域滤波中的零陷方法不能在卫星信号方向提供增益，卫星信号也有可能被抑制。空域滤波中的波束形成方法和空时自适应处理往往需要卫星到达接收机的角度。GPS/SINS组合抗干扰方法需要GPS接收机能够正确输出伪距、伪距率或位置、速度信息。

发明内容

本发明的目的在于提供一种抑制干扰能力强的基于阵列天线和GPS/SINS的联合抗干扰方法。

本发明的目的是这样实现的：

初始化载体的位置和姿态之后，建立GPS/SINS组合导航状态方程和量测方程；GPS/SINS组合导航实时提供载体的位置和姿态，同时根据卫星星历信息计算出当前卫星的位置，从而获得卫星到载体之间的导向矢量，即卫星到达载体的方向角和俯仰角；然后，所述导向矢量作为多约束最小方差空时自适应处理(MCMV-STAP)算法的先验信息，在空域、时域同时抑制宽带干扰和窄带干扰。具体包括如下步骤：

步骤1：初始化载***置、速度和姿态信息

在大地坐标系中，设定载体初始时刻的坐标：纬度L、经度λ和高度h；初始化载体在东北天坐标系中的速度：东向速度V_E、北向速度V_N和天向速度V_U；初始化载体的姿态角，包括俯仰角θ、横滚角γ和方位角ψ；然后设定载体的飞行路径，可以设定为静止、直线运动或圆周运动等轨迹。从而获得理想陀螺仪和加速度计的输出，f_E、f_N和f_U分别表示为加速度计在东向、北向和天向的输出。

步骤2：初始化导航滤波器的参数信息

在GPS/SINS组合导航滤波器中，采用反馈校正方式。导航滤波器的状态量为姿态角误差、速度误差、位置误差、陀螺仪误差和加速度计误差。φ_E、φ_N和φ_U分别表示为载体俯仰角误差、横滚角误差和方位角误差，δV_E、δV_N和δV_U分别表示为载体的东向速度误差、北向速度误差和天向速度误差，δL、δλ和δh分别表示为载体的纬度误差、经度误差和高度误差。ε_E、ε_N和ε_U分别表示为陀螺仪在东向、北向和天向的漂移。和分别表示为加速度计在东向、北向和天向的输出误差。

步骤3：计算陀螺仪和加速度计的误差变化率

步骤4：根据步骤1-3中的参数，计算载体姿态角误差变化率、速度误差变化率和位置误差变化率

步骤5：引入GPS伪距量测信息，采用反馈校正方式对SINS输出信息进行校正，获得当前准确的位置和姿态。

根据GPS接收机跟踪环路所测量的码相位误差计算出卫星到载体的之间的伪距，伪距作为量测信息来更新导航滤波器的状态量，从而预测出当前SINS所量测状态量的误差。然后采用反馈校正方式对SINS输出信息进行校正，获得当前准确的位置和姿态。

步骤6：计算载体在地心地固坐标系中的坐标

由步骤5计算出的载***置位于大地坐标系中，即(L,λ,h)，将其转化为地心地固坐标系中，其位置坐标可以表示为

$X_p^e={[(R+h)\cos L\cos \mathrm{\λ},(R+h)\cos L\sin \mathrm{\λ},(R+h)\sin L]}^T---\left(1\right)$

式(1)中，R为地球半径。

步骤7：计算载体坐标系(b系)到导航坐标系(n系)的转换矩阵、导航坐标系到地心地固坐标系(e系)的转换矩阵

由步骤5计算出的载体姿态角，可以计算出载体坐标系到导航坐标系的转换矩阵为

$C_b^n=\left[\begin{array}{ccc}\cos \mathrm{\γ}\cos \mathrm{\ψ}-\sin \mathrm{\ψ}\sin \mathrm{\θ}\sin \mathrm{\γ}& -\sin \mathrm{\ψ}\cos \mathrm{\θ}& \sin \mathrm{\γ}\cos \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}\\ \cos \mathrm{\γ}\sin \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}& \cos \mathrm{\ψ}\cos \mathrm{\θ}& \sin \mathrm{\γ}\sin \mathrm{\ψ}-\cos \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}\\ -\sin \mathrm{\γ}\cos \mathrm{\θ}& \sin \mathrm{\θ}& \cos \mathrm{\θ\γ}\cos \end{array}\right]---\left(2\right)$

由步骤5计算出的载***置，可以计算出导航坐标系到地心地固坐标系的转换矩阵为

$C_n^e=\left[\begin{array}{ccc}-\sin \mathrm{\λ}& -\sin L\cos \mathrm{\λ}& \cos L\cos \mathrm{\λ}\\ \cos \mathrm{\λ}& -\sin L\sin \mathrm{\λ}& \cos L\sin \mathrm{\λ}\\ 0& \cos L& \sin L\end{array}\right]---\left(3\right)$

步骤8：计算载体到卫星之间的导向矢量

由步骤7中的(2)、(3)可以求得载体坐标系到地心地固坐标系的转化矩阵

$C_b^e=C_n^e{C}_b^n---\left(4\right)$

利用卫星星历解算出卫星在地心地固坐标系中的位置并结合公式(1)(4),可以求得载体到卫星之间的导向矢量

${\stackrel{\→}{r}}^b={\left(C_b^e\right)}^T(X_s^e-X_p^e)---\left(5\right)$

步骤9：在载体坐标系中，计算卫星到达天线的方位角和俯仰角

若将在载体坐标系中的坐标定义为则卫星到达天线的方位角α和俯仰角β分别可以表示为

$\mathrm{\α}=\arctan (\hat{x},\hat{y})---\left(6\right)$

$\mathrm{\β}=\arctan (\hat{z},\sqrt{{\hat{x}}^2+{\hat{y}}^2})---\left(7\right)$

步骤10：设计圆形阵列天线结构

为了能够在方位角和俯仰角方向同时控制波束指向可视卫星方向，本发明采用圆形阵列天线结构。6个阵元均匀分布于圆周上的圆阵，令圆半径为r，则圆周上相邻阵元间的间隔也为r。阵元间隔的选取，与时域采样间隔满足奈奎斯特定理一样，空域采样间隔d应小于卫星载波波长λ的1/2。由卫星信号频率f=1575.42×10⁶MHz，所以阵元间距为

$d\≤\frac{\mathrm{\λ}}{2}=\frac{1}{2}\·\frac{c}{f_{l1}}=\frac{1}{2}\·\frac{3\×{10}^8}{1575.42\×{10}^6}=\frac{1}{2}\·0.19=0.095m=9.5\mathrm{cm}---\left(8\right)$

式(8)中，c为光速。

为了使主波束宽度越窄，副瓣越低，分辨率高，需要使阵元间隔尽量的大，所以取圆半径r=d=9.5cm，整个阵列天线的直径约为19cm。

步骤11：计算卫星信号到达天线各阵元之间的时间延迟

根据步骤9中的(6)(7)，可以将α和β表示单位矢量

e(α,β)=(sinαcosβ,cosαcosβ,sinβ)^T(9)

因此，卫星信号到达第i个天线阵元与到第一个参考阵元之间的时间差τ_i可以表示为

τ_i=e^T·(x_i-x₁)/ci=1,2,…M-1(10)

式(10)中，M代表天线阵元数目。

步骤12：建立阵列天线接收信号模型

用户一般可以接收4颗以上的卫星信号，因此，需要形成多个波束指向对应的卫星。假设阵列天线接收到了P个卫星信号，Q个干扰信号，则天线接收到的信号模型可以表示为

$U\left(t\right)=\underset{k=1}{\overset{P}{\mathrm{\Σ}}}a({\mathrm{\α}}_k,{\mathrm{\β}}_k,T_s)s_k\left(t\right)+\underset{k=P+1}{\overset{P+Q}{\mathrm{\Σ}}}a({\mathrm{\α}}_k,{\mathrm{\β}}_k,T_s)j_{k-P}\left(t\right)+n\left(t\right)---\left(11\right)$

在公式(11)中，s(t)和j(t)分别表示接收到的卫星信号和干扰信号，a(α_k,β_k,T_s)是第k个目标信号(卫星信号或干扰信号)的导向矢量。T_s是时域延迟线间隔，α_k和β_k分别表示为第k个目标信号到达阵列天线的方位角和俯仰角。n(t)表示高斯白噪声，它的功率谱密度表示为N₀/2。

在公式(11)中，a(α_k,β_k,T_s)表示空时二维目标矢量，即时间矢量a_s(T_s)和空间方向矢量a_s(α_k,β_k)的克罗奈克积（KroneckerProduct），并可以表示为：

$a({\mathrm{\α}}_k,{\mathrm{\β}}_k,T_s)=a({\mathrm{\α}}_k,{\mathrm{\β}}_k)\⊗a_s\left(T_s\right)---\left(12\right)$

令每个射频通道的时间延迟单元数目为N，则

$a_s({\mathrm{\α}}_k,{\mathrm{\β}}_k)=\left[\begin{array}{cccc}1& e^{-j2\mathrm{\πf}{\mathrm{\τ}}_1}& \·\·\·& e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{M-1}}\end{array}\right]---\left(13\right)$

$a_t=\left[\begin{array}{cccc}1& e^{-j2\mathrm{\πfT}}& \·\·\·& e^{-j2\mathrm{\πf}(N-1)T_s}\end{array}\right]---\left(14\right)$

式(13)中，f表示卫星信号频率，τ_i(i=1,2,…M-1)由式(10)给出。式(14)中，T_s表示延迟单元的时间间隔，其值应该小于信号带宽。

步骤13：计算天线阵列权值矢量

采用多约束最小方差空时自适应处理(MCMV-STAP)算法对可视卫星信号进行约束，然后使得阵列输出功率最小，从而保护卫星信号的同时抑制干扰信号。该算法需要求得阵列的最优权值，表达如下：

$\left\{\begin{array}{c}\mathrm{\ω}=\arg \min {\mathrm{\ω}}^H{R}_U\mathrm{\ω}\\ \mathrm{subjectto}:{\mathrm{\ω}}^{H}A=F^T\end{array}\right.---\left(15\right)$

式(15)中，ω表示阵列的权值矢量，H表示矩阵的共轭转置，R_U是输入信号的阵列协方差矩阵，可以表示为

R_U=E{UU^H}(16)

式(15)中，A表示卫星信号的约束矩阵，F表示与A相对应的约束矢量，A和F分别表示为

A=[a(α₁,β₁,T_s),a(α₂,β₂,T_s)…a(α_P,β_P,T_s)](17)

f^T=[1,1…1]_1×p(18)采用拉格朗日乘子法，阵列权值矢量ω可以表示为

$\mathrm{\ω}=R_U^{-1}A{\left(A^H{R}_U^{-1}A\right)}^{-1}f---\left(19\right)$

在获得阵列权值之后，可以得到阵列的输出表达是为

y(t)=ω^HU(t)。

传统的空时自适应处理方法，通常采用盲自适应波束形成算法，即不需要知道卫星到达阵列天线的方位角和俯仰角，只是简单的将某一天线阵元作为参考阵元进行权值约束。然而，由于未在卫星信号方向进行约束，不能在卫星信号方向形成波束。因此，该类盲波束形成算法在抑制干扰信号的同时也将削弱期望的卫星信号，卫星信号功率将被减弱，不能够得到信号的最大信干噪比。

本发明主要利用GPS/SINS组合导航***提供的载***置和姿态，GPS星历提供的卫星位置，从而获得卫星到载体之间的导向矢量，即卫星到达载体的方向角和俯仰角。然后，该导向矢量作为多约束最小方差空时自适应处理(MCMV-STAP)算法的先验信息，在空域、时域同时抑制宽带干扰和窄带干扰。

本发明所述的计算卫星到载体之间的导向矢量方法，其特征在于GPS/SINS所解算出载体的位置位于大地坐标系中，姿态位于东北天坐标系中，利用卫星星历解算的卫星位置在地心地固坐标系中，而波束形成所需的导向矢量是应该在载体坐标系下提供的。本发明将所求得的先验信息(卫星的位置、载体的位置和姿态)经过坐标变化后，在载体坐标系上获得卫星到达载体的方向角和俯仰角。

本发明采用阵列天线与GPS/SINS相结合的方式来提高导航接收机的抗干扰能力。在阵列天线中，采用多约束最小方差准则的空时自适应波束形成算法(MCMV-STAP)，该算法能够在多颗可视卫星方向形成波束，在干扰方向形成零陷，从而增强卫星信号的同时抑制干扰信号。为了获得卫星信号的最大功率，该算法需要知道卫星信号到达载体的方位角和俯仰角。因此，本发明采用圆形阵列天线结构，引入GPS/SINS组合导航为阵列天线的波束形成提供载体的位置和姿态，采用卫星星历提供卫星的位置，从而为波束形成提供先验信息，即可视卫星到达载体的方位角和俯仰角。

附图说明

图1是基于阵列天线和GPS/SINS的联合抗干扰方法的实施流程图；

图2是阵列天线、东北天坐标系和载体坐标系的示意图；

图3是多约束最小方差空时二维自适应处理(MCMV-STAP)的结构图。

具体实施方式

下面结合附图对本发明做进一步的阐述。

步骤1：初始化载***置、速度和姿态信息

如附图1所示，需要先初始化载体的位置与姿态。在大地坐标系中，设定载体初始时刻的坐标：纬度L=45°、经度λ=126°和高度h=200m。初始化载体的姿态角，包括俯仰角θ=0°、横滚角γ=0°和方位角ψ=45°。然后设定载体的飞行路径，可以设定为静止、直线运动或圆周运动等轨迹。根据要求的载体飞行路径设定载体的东向速度V_E、北向速度V_N和天向速度V_U。从而获得理想陀螺仪和加速度计的输出，f_E、f_N和f_U分别表示为加速度计在东向、北向和天向的输出。

步骤2：初始化导航滤波器的参数信息

在GPS/SINS组合导航滤波器中，采用反馈校正方式。导航滤波器的状态量为姿态角误差、速度误差、位置误差、陀螺仪误差和加速度计误差。在初始时刻，载体俯仰角误差φ_E=0.1°、横滚角误差φ_N=0.1°和方位角误差φ_U=1°。载体的东向速度误差δV_E=0.2m/s、北向速度误差δV_N=0.2m/s和天向速度误差δV_U=0.5m/s。载体的纬度误差δL=0°、经度误差δλ=0°和高度误差δh=20m。陀螺仪常值漂移ε_b=0.1°/h,白噪声ε_g=0.05°/h，加速度计常值误差白噪声w_a=5×10^-4g，g为重力加速度。陀螺仪误差ε在东向、北向和天向的漂移分别表示为ε_E、ε_N和ε_U。加速度计在东向、北向和天向的输出误差分别可以表示为 和

步骤3：计算陀螺仪和加速度计的误差变化率

陀螺仪的常值漂移可以以用随机常数描述为

${\stackrel{\·}{\mathrm{\ϵ}}}_b=0---\left(1\right)$

陀螺仪的噪声ε_g来可以用狄拉克函数描述。因此，陀螺误差变化率可以表示为

$\stackrel{\·}{\mathrm{\ϵ}}={\stackrel{\·}{\mathrm{\ϵ}}}_b+{\mathrm{\ϵ}}_g---\left(2\right)$

对于加速度计误差变化率可以将其考虑为一阶马尔柯夫过程，误差模型取为

$\stackrel{\·}{\▿}=-\frac{1}{T_a}\▿+w_a---\left(3\right)$

其中，T_a表示相关时间，w_a为白噪声过程。

步骤4：根据步骤1-3中的参数，计算载体姿态角误差变化率、速度误差变化率和位置误差变化率

1、计算载体姿态角误差变化率

和分别表示为载体俯仰角误差、横滚角误差和方位角误差的变化率，载体姿态角误差变化率可以表示为

$\begin{array}{c}{\stackrel{\·}{\mathrm{\φ}}}_E=-\frac{{\mathrm{\δV}}_N}{R+h}+({\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_E}{R+h}\tan L){\mathrm{\φ}}_N-({\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R+h}\tan L){\mathrm{\φ}}_U+{\mathrm{\ϵ}}_E\\ {\stackrel{\·}{\mathrm{\φ}}}_N=\frac{{\mathrm{\δV}}_E}{R+h}-{\mathrm{\ω}}_{\mathrm{ie}}\sin \mathrm{L\δL}-({\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_E}{R+h}\tan L){\mathrm{\φ}}_E-\frac{V_E}{R+h}{\mathrm{\φ}}_U+{\mathrm{\ϵ}}_N\\ {\stackrel{\·}{\mathrm{\φ}}}_U=\frac{{\mathrm{\δV}}_E}{R+h}\tan L+({\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R+h}{\sec }^{2}L)\mathrm{\δL}+({\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R+h}){\mathrm{\φ}}_E+\frac{V_N}{R+h}{\mathrm{\φ}}_N+{\mathrm{\ϵ}}_U\end{array}---\left(4\right)$

式(4)中，R表示地球半径，ω_ie表示地球旋转角速率。

2、计算载体速度误差的变化率

和分别表示为载体的东向速度误差、北向速度误差和天向速度误差的变化率，载体速度误差变化率可以表示为

$\begin{array}{c}{\stackrel{\·}{\mathrm{\δV}}}_E=(\frac{V_N}{R+h}\tan L-\frac{V_U}{R+h}){\mathrm{\δV}}_E+({2\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_E}{R+h}\tan L){\mathrm{\δV}}_N\\ -(-{\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R+h}){\mathrm{\δV}}_U+({2\mathrm{\ω}}_{\mathrm{ie}}{V}_N\cos L+\frac{V_E{V}_N}{R+h}{\sec }^{2}L+2{\mathrm{\ω}}_{\mathrm{ie}}{V}_U\sin L)\mathrm{\δL}\\ +{\▿}_E+f_N{\mathrm{\φ}}_U-f_U{\mathrm{\φ}}_N\\ {\stackrel{\·}{\mathrm{\δV}}}_N=-2({\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_E}{R+h}\tan L){\mathrm{\δV}}_E-\frac{V_U}{R+h}{\mathrm{\δV}}_N-\frac{V_N}{R+h}{\mathrm{\δV}}_U\\ -(2{\mathrm{\ω}}_{\mathrm{ie}}{V}_E\cos L+\frac{V_E^2{\sec }^{2}L}{R+h})\mathrm{\δL}+{\▿}_N-f_E{\mathrm{\φ}}_U+f_U{\mathrm{\φ}}_E\\ {\stackrel{\·}{\mathrm{\δV}}}_U=(2{\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R+h}){\mathrm{\δV}}_E+\frac{{2V}_N}{R+h}{\mathrm{\δV}}_N-2{\mathrm{\ω}}_{\mathrm{ie}}{V}_E\sin \mathrm{L\δL}+{\▿}_U+f_E{\mathrm{\φ}}_N-f_N{\mathrm{\φ}}_E\end{array}---\left(5\right)$

3、计算载***置误差的变化率

和分别表示为载体的纬度误差、经度误差和高度误差的变化率。载***置误差变化率可以表示为

$\begin{array}{c}\mathrm{\δ}\stackrel{\·}{L}=\frac{{\mathrm{\δV}}_N}{R+h}\\ \mathrm{\δ}\stackrel{\·}{\mathrm{\λ}}=\frac{{\mathrm{\δV}}_E}{R+h}\\ \mathrm{\δ}\stackrel{\·}{h}=\mathrm{\δ}{V}_U\end{array},\sec L+\frac{V_E}{R+h}\sec \mathrm{LtgL\δL}---\left(6\right)$

步骤5：引入GPS伪距量测信息，采用反馈校正方式对SINS输出信息进行校正，获得当前

准确的位置和姿态。

根据GPS接收机跟踪环路所测量的码相位误差计算出卫星到载体的之间的伪距，伪距作为量测信息来更新导航滤波器的状态量，从而预测出当前SINS所量测的位置误差(δL,δλ,δh)和姿态误差(δθ,δγ,δψ)。然后采用反馈校正方式对SINS输出的位置(L,λ,h)和姿态(θ,γ,ψ)进行校正，获得当前准确的位置和姿态。

L=L+δL

λ=λ+δλ

h=h+δh

θ=θ+δθ(7)

γ=γ+δγ

ψ=ψ+δψ

步骤6：计算载体在地心地固坐标系中的坐标

由步骤5计算出的载***置位于大地坐标系中，即(L,λ,h)，将其转化为地心地固坐标系中，其位置坐标可以表示为

$X_p^e={[(R+h)\cos L\cos \mathrm{\λ},(R+h)\cos L\sin \mathrm{\λ},(R+h)\sin L]}^T---\left(8\right)$

步骤7：计算载体坐标系(b系)到导航坐标系(n系)的转换矩阵、导航坐标系到地心地固坐标系(e系)的转换矩阵

由步骤5计算出的载体姿态可以计算出载体坐标系到导航坐标系的转换矩阵为

$C_b^n=\left[\begin{array}{ccc}\cos \mathrm{\γ}\cos \mathrm{\ψ}-\sin \mathrm{\ψ}\sin \mathrm{\θ}\sin \mathrm{\γ}& -\sin \mathrm{\ψ}\cos \mathrm{\θ}& \sin \mathrm{\γ}\cos \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}\\ \cos \mathrm{\γ}\sin \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}& \cos \mathrm{\ψ}\cos \mathrm{\θ}& \sin \mathrm{\γ}\sin \mathrm{\ψ}-\cos \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}\\ -\sin \mathrm{\γ}\cos \mathrm{\θ}& \sin \mathrm{\θ}& \cos \mathrm{\θ\γ}\cos \end{array}\right]---\left(9\right)$

由步骤5计算出的载***置，可以计算出导航坐标系到地心地固坐标系的转换矩阵为

$C_n^e=\left[\begin{array}{ccc}-\sin \mathrm{\λ}& -\sin L\cos \mathrm{\λ}& \cos L\cos \mathrm{\λ}\\ \cos \mathrm{\λ}& -\sin L\sin \mathrm{\λ}& \cos L\sin \mathrm{\λ}\\ 0& \cos L& \sin L\end{array}\right]---\left(10\right)$

步骤8：计算载体到卫星之间的导向矢量

由步骤7中的(9)、(10)可以求得载体坐标系到地心地固坐标系的转化矩阵

$C_b^e=C_n^e{C}_b^n---\left(11\right)$

利用卫星星历解算出卫星在地心地固坐标系中的位置并结合公式(8)(11),可以求得载体到卫星之间的导向矢量

${\stackrel{\→}{r}}^b={\left(C_b^e\right)}^T(X_s^e-X_p^e)---\left(12\right)$

步骤9：在载体坐标系中，计算卫星到达天线的方位角和俯仰角

若将在载体坐标系中的坐标定义为则卫星到达天线的方位角α和俯仰角β分别可以表示为

$\mathrm{\α}=\arctan (\hat{x},\hat{y})---\left(13\right)$

$\mathrm{\β}=\arctan (\hat{z},\sqrt{{\hat{x}}^2+{\hat{y}}^2})---\left(14\right)$

步骤10：设计圆形阵列天线结构

为了能够在方位角和俯仰角方向同时控制波束指向可视卫星方向，本发明采用圆形阵列天线结构，如附图2所示。6个阵元均匀分布于圆周上的圆阵，令圆半径为r，则圆周上相邻阵元间的间隔也为r。阵元间隔的选取，与时域采样间隔满足奈奎斯特定理一样，空域采样间隔d应小于卫星载波波长λ的1/2。由卫星信号频率f=1575.42×10⁶MHz，所以阵元间距为

$d\≤\frac{\mathrm{\λ}}{2}=\frac{1}{2}\·\frac{c}{f_{l1}}=\frac{1}{2}\·\frac{3\×{10}^8}{1575.42\×{10}^6}=\frac{1}{2}\·0.19=0.095m=9.5\mathrm{cm}---\left(15\right)$

式(15)中，c为光速。

为了使主波束宽度越窄，副瓣越低，分辨率高，需要使阵元间隔尽量的大，所以取圆半径r=d=9.5cm，整个阵列天线的直径约为19cm。

以天线阵所在平面为xoy平面，以坐标原点为参考，x轴阵元为1号阵元，则6阵元的极坐标分别为：(r,0)、(r,π/3)、(r,2π/3)、(r,π)、(r,4π/3)、(r,5π/3)。

步骤11：计算卫星信号到达天线各阵元之间的时间延迟

根据步骤9中的(13)(14)，可以将α和β表示单位矢量

e(α,β)=(sinαcosβ,cosαcosβ,sinβ)^T(16)

因此，卫星信号到达第i个天线阵元与到第一个参考阵元之间的时间差τi可以表示为

τ_i=e^T·(x_i-x₁)/ci=1,2,…M-1(17)式(17)中，M代表天线阵元数目。

步骤12：建立阵列天线接收信号模型

用户一般可以接收4颗以上的卫星信号，因此，需要形成多个波束指向对应的卫星。假设阵列天线接收到了P个卫星信号，Q个干扰信号，则天线接收到的信号模型可以表示为

$U\left(t\right)=\underset{k=1}{\overset{P}{\mathrm{\Σ}}}a({\mathrm{\α}}_k,{\mathrm{\β}}_k,T_s)s_k\left(t\right)+\underset{k=P+1}{\overset{P+Q}{\mathrm{\Σ}}}a({\mathrm{\α}}_k,{\mathrm{\β}}_k,T_s)j_{k-P}\left(t\right)+n\left(t\right)---\left(18\right)$

在公式(18)中，s(t)和j(t)分别表示接收到的卫星信号和干扰信号，a(α_k,β_k,T_s)是第k个目标信号(卫星信号或干扰信号)的导向矢量。T_s是时域延迟线间隔，α_k和β_k分别表示为第k个目标信号到达阵列天线的方位角和俯仰角。n(t)表示高斯白噪声，它的功率谱密度表示为N₀/2。

在公式(18)中，a(α_k,β_k,T_s)表示空时二维目标矢量，即时间矢量a_s(T_s)和空间方向矢量a_s(α_k,β_k)的克罗奈克积（KroneckerProduct），并可以表示为：

$a({\mathrm{\α}}_k,{\mathrm{\β}}_k,T_s)=a_s({\mathrm{\α}}_k,{\mathrm{\β}}_k)\⊗a_s\left(T_s\right)---\left(19\right)$

令每个射频通道的时间延迟单元数目为N，则

$a_s({\mathrm{\α}}_k,{\mathrm{\β}}_k)=\left[\begin{array}{cccc}1& e^{-j2\mathrm{\πf}{\mathrm{\τ}}_1}& \·\·\·& e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{M-1}}\end{array}\right]---\left(20\right)$

$a_t=\left[\begin{array}{cccc}1& e^{-j2\mathrm{\πfT}}& \·\·\·& e^{-j2\mathrm{\πf}(N-1)T_s}\end{array}\right]---\left(21\right)$

式(20)中，f表示卫星信号频率，τ_i(i=1,2,…M-1)由式(17)给出。式(21)中，T_s表示延迟单元的时间间隔，其值应该小于信号带宽。

步骤13：计算天线阵列权值矢量

附图3为多约束最小方差空时二维自适应处理(MCMV-STAP)的结构图。从每个阵元的通道来看，各级延时构成了FIR滤波器，可在时域去除干扰；从相同的时间延迟节点看，不同的阵元构成了空域的自适应滤波，可以分辨空间干扰源进而形成空域零陷从空域上抑制干扰。而空域的处理也可以进一步利用时域处理后的反馈信息，空时处理也因此具有在空时二维域同时剔除干扰的能力。

采用MCMV-STAP算法对可视卫星信号进行约束，然后使得阵列输出功率最小，从而保护卫星信号的同时抑制干扰信号。该算法需要求得阵列的最优权值，表达如下：

$\left\{\begin{array}{c}\mathrm{\ω}=\arg \min {\mathrm{\ω}}^H{R}_U\mathrm{\ω}\\ \mathrm{subjectto}:{\mathrm{\ω}}^{H}A=F^T\end{array}\right.---\left(22\right)$

式(22)中，ω表示阵列的权值矢量，H表示矩阵的共轭转置，R_U是输入信号的阵列协方差矩阵，可以表示为

R_U=E{UU^H}(23)

式(23)中，A表示卫星信号的约束矩阵，F表示与A相对应的约束矢量，A和F分别表示为

A=[a(α₁,β₁,T_s),a(α₂,β₂,T_s)…a(α_P,β_P,T_s)](24)

f^T=[1,1…1]_1×p(25)采用拉格朗日乘子法，阵列权值矢量ω可以表示为

$\mathrm{\ω}=R_U^{-1}A{\left(A^H{R}_U^{-1}A\right)}^{-1}f---\left(26\right)$

在获得阵列权值之后，可以得到阵列的输出表达是为

y(t)=ω^HU(t)。(27)。

Claims (1)

1.一种基于阵列天线和GPS/SINS的联合抗干扰方法，初始化载体的位置和姿态之后，建立GPS/SINS组合导航状态方程和量测方程；GPS/SINS组合导航实时提供载体的位置和姿态，同时根据卫星星历信息计算出当前卫星的位置，从而获得卫星到载体之间的导向矢量，即卫星到达载体的方向角和俯仰角；然后，所述导向矢量作为多约束最小方差空时自适应处理算法的先验信息，在空域、时域同时抑制宽带干扰和窄带干扰；其特征是具体包括如下步骤：

步骤1：初始化载***置、速度和姿态信息；

在大地坐标系中，设定载体初始时刻的坐标：纬度L、经度λ和高度h；初始化载体在东北天坐标系中的速度：东向速度V_E、北向速度V_N和天向速度V_U；初始化载体的姿态角，包括俯仰角θ、横滚角γ和方位角ψ；然后设定载体的飞行路径为静止、直线运动或圆周运动轨迹中的一种；从而获得陀螺仪和加速度计的输出，f_E、f_N和f_U分别表示为加速度计在东向、北向和天向的输出；

步骤2：初始化导航滤波器的参数信息；

在GPS/SINS组合导航滤波器中，采用反馈校正方式；导航滤波器的状态量为姿态角误差、速度误差、位置误差、陀螺仪误差和加速度计误差，φ_E、φ_N和φ_U分别表示为载体俯仰角误差、横滚角误差和方位角误差，δV_E、δV_N和δV_U分别表示为载体的东向速度误差、北向速度误差和天向速度误差，δL、δλ和δh分别表示为载体的纬度误差、经度误差和高度误差，ε_E、ε_N和ε_U分别表示为陀螺仪在东向、北向和天向的漂移，▽_E、▽_N和▽_U分别表示为加速度计在东向、北向和天向的输出误差；

步骤3：计算陀螺仪和加速度计的误差变化率；

步骤4：根据步骤1-3中的参数，计算载体姿态角误差变化率、速度误差变化率和位置误差变化率；

步骤5：引入GPS伪距量测信息，采用反馈校正方式对SINS输出信息进行校正，获得当前准确的位置和姿态；

根据GPS接收机跟踪环路所测量的码相位误差计算出卫星到载体的之间的伪距，伪距作为量测信息来更新导航滤波器的状态量，从而预测出当前SINS所量测状态量的误差，然后采用反馈校正方式对SINS输出信息进行校正，获得当前准确的位置和姿态；

步骤6：计算载体在地心地固坐标系中的坐标；

由步骤5计算出的载***置位于大地坐标系中，即(L,λ,h)，将其转化为地心地固坐标系中，其位置坐标表示为

$X_p^e={\[(R+h)\cos Lcos\mathrm{\λ},(R+h)\cos Lsin\mathrm{\λ},(R+h)\sin L\]}^T$

其中，R为地球半径；

步骤7：计算载体坐标系即b系到导航坐标系即n系的转换矩阵、导航坐标系到地心地固坐标系即e系的转换矩阵；

由步骤5计算出的载体姿态角，计算出载体坐标系到导航坐标系的转换矩阵为

$C_b^n=\left[\begin{array}{ccc}\cos \mathrm{\γ}\cos \mathrm{\ψ}-\sin \mathrm{\ψ}\sin \mathrm{\θ}\sin \mathrm{\γ}& -\sin \mathrm{\ψ}\cos \mathrm{\θ}& \sin \mathrm{\γ}\cos \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}\\ \cos \mathrm{\γ}\sin \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}& \cos \mathrm{\ψ}\cos \mathrm{\θ}& \sin \mathrm{\γ}\sin \mathrm{\ψ}-\cos \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}\\ -\sin \mathrm{\γ}\cos \mathrm{\θ}& \sin \mathrm{\θ}& \cos \mathrm{\θ}\cos \mathrm{\γ}\end{array}\right]$

由步骤5计算出的载***置，计算出导航坐标系到地心地固坐标系的转换矩阵为

$C_n^e=\left[\begin{array}{ccc}-\sin \mathrm{\λ}& -\sin L\cos \mathrm{\λ}& \cos L\cos \mathrm{\λ}\\ \cos \mathrm{\λ}& -\sin L\cos \mathrm{\λ}& \cos L\sin \mathrm{\λ}\\ 0& \cos L& \sin L\end{array}\right];$

步骤8：计算载体到卫星之间的导向矢量；

由步骤7中的转换矩阵转换矩阵求得载体坐标系到地心地固坐标系的转化矩阵

$C_b^e=C_n^e{C}_b^n$

利用卫星星历解算出卫星在地心地固坐标系中的位置并结合公式 $X_p^e={\[(R+h)\cos Lcos\mathrm{\λ},(R+h)\cos L\sin \mathrm{\λ},(R+h)\sin L\]}^T$ 和 $C_b^e=C_n^e{C}_b^n,$ 求得载体到卫星之间的导向矢量

${\stackrel{\→}{r}}^b={\left(C_b^e\right)}^T(X_s^e-X_p^e)$

步骤9：在载体坐标系中，计算卫星到达天线的方位角和俯仰角；

若将在载体坐标系中的坐标定义为则卫星到达天线的方位角α和俯仰角β分别可以表示为

$\mathrm{\α}=arctan(\hat{x},\hat{y})$

$\mathrm{\β}=arctan(\hat{z},\sqrt{{\hat{x}}^2+{\hat{y}}^2})$

步骤10：设计圆形阵列天线结构；

6个阵元均匀分布于圆周上的圆阵，令圆半径为r，则圆周上相邻阵元间的间隔也为r，阵元间隔的选取，与时域采样间隔满足奈奎斯特定理一样，空域采样间隔d应小于卫星载波波长λ的1/2，由卫星信号频率f＝1575.42×10⁶MHz，所以阵元间距为

$d\≤\frac{\mathrm{\λ}}{2}=\frac{1}{2}\·\frac{c}{f_{l1}}=\frac{1}{2}\·\frac{3\×{10}^8}{1575.42\×{10}^6}=\frac{1}{2}\·0.19=0.095m=9.5cm$

其中，c为光速；

步骤11：计算卫星信号到达天线各阵元之间的时间延迟；

根据步骤9中的 $\mathrm{\α}=arctan(\hat{x},\hat{y}),$ $\mathrm{\β}=arctan(\hat{z},\sqrt{{\hat{x}}^2+{\hat{y}}^2}),$ 将α和β表示单位矢量

e(α,β)＝(sinαcosβ,cosαcosβ,sinβ)^T

因此，卫星信号到达第i个天线阵元与到第一个参考阵元之间的时间差τ_i可以表示为

τ_i＝e^T·(x_i-x₁)/ci＝1,2,…M-1

其中，M代表天线阵元数目；

步骤12：建立阵列天线接收信号模型

设阵列天线接收到了P个卫星信号，Q个干扰信号，则天线接收到的信号模型表示为

$U\left(t\right)=\underset{k=1}{\overset{P}{\Σ}}a({\mathrm{\α}}_k,{\mathrm{\β}}_k,T_s)s_k\left(t\right)+\underset{k=P+1}{\overset{P+Q}{\Σ}}a({\mathrm{\α}}_k,{\mathrm{\β}}_k,T_s)j_{k-P}\left(t\right)+n\left(t\right)$

其中，s(t)和j(t)分别表示接收到的卫星信号和干扰信号，a(α_k,β_k,T_s)是第k个目标信号的导向矢量，T_s是时域延迟线间隔，α_k和β_k分别表示为第k个目标信号到达阵列天线的方位角和俯仰角，n(t)表示高斯白噪声、它的功率谱密度表示为N₀/2；

a(α_k,β_k,T_s)表示空时二维目标矢量，即时间矢量a_s(T_s)和空间方向矢量a_s(α_k,β_k)的克罗奈克积，并表示为：

$a({\mathrm{\α}}_k,{\mathrm{\β}}_k,T_s)=a_s({\mathrm{\α}}_k,{\mathrm{\β}}_k)\⊗a_s\left(T_s\right)$

令每个射频通道的时间延迟单元数目为N，则

$a_s({\mathrm{\α}}_k,{\mathrm{\β}}_k)=\left[\begin{array}{cccc}1& e^{-j2{\mathrm{\πf\τ}}_1}& \mathrm{...}& e^{-j2{\mathrm{\πf\τ}}_{M-1}}\end{array}\right]$

a_t＝[1e^-j2πfT…e^{-j2πf(N-1)Ts}]

f表示卫星信号频率，T_s表示延迟单元的时间间隔、其值应该小于信号带宽；

步骤13：计算天线阵列权值矢量；

采用多约束最小方差空时自适应处理算法对可视卫星信号进行约束，然后使得阵列输出功率最小，阵列的最优权值的表达如下：

$\left\{\begin{array}{c}\mathrm{\ω}=\arg {\mathrm{min\ω}}^H{R}_U\mathrm{\ω}\\ subjectto:{\mathrm{\ω}}^{H}A=F^T\end{array}\right.$

其中，ω表示阵列的权值矢量，H表示矩阵的共轭转置，R_U是输入信号的阵列协方差矩阵、表示为

R_U＝E{UU^H}

A表示卫星信号的约束矩阵，F表示与A相对应的约束矢量，A和F分别表示为

A＝[a(α₁,β₁,T_s),a(α₂,β₂,T_s)…a(α_P,β_P,T_s)]

f^T＝[1,1…1]_1×p

采用拉格朗日乘子法，阵列权值矢量ω表示为

$\mathrm{\ω}=R_U^{-1}A{\left(A^H{R}_U^{-1}A\right)}^{-1}f$

在获得阵列权值之后，得到阵列的输出表达是为

y(t)＝ω^HU(t)。