CN103792550B

CN103792550B - A kind of associating anti-interference method based on array antenna and GPS/SINS

Info

Publication number: CN103792550B
Application number: CN201410047886.5A
Authority: CN
Inventors: 王伟; 李强; 徐定杰; 沈锋; 王咸鹏; 刘明凯; 范岳; 刘海峰; 宋金阳; 桑静
Original assignee: Harbin Engineering University
Current assignee: Harbin Engineering University
Priority date: 2014-02-11
Filing date: 2014-02-11
Publication date: 2015-12-02
Anticipated expiration: 2034-02-11
Also published as: CN103792550A

Abstract

The invention provides a joint anti-jamming method based on array antenna and GPS/SINS. After initializing the position and attitude of the carrier, establish the GPS/SINS integrated navigation state equation and measurement equation; GPS/SINS integrated navigation provides the position and attitude of the carrier in real time, calculates the current satellite position according to the satellite ephemeris information, and obtains the satellite to carrier The steering vector between them; the steering vector is used as the prior information of the multi-constrained minimum variance space-time adaptive processing algorithm, and simultaneously suppresses wideband interference and narrowband interference in space domain and time domain. The invention can simultaneously form beams in the direction of multiple visible satellites and form nulls in the direction of interference, thereby enhancing satellite signals and suppressing interference signals. The invention adopts a circular structure antenna array, GPS/SINS integrated navigation provides the position and attitude of the carrier for the beam forming of the array antenna, and uses the satellite ephemeris to provide the position of the satellite, thereby providing prior information for the beam forming.

Description

A joint anti-jamming method based on array antenna and GPS/SINS

技术领域 technical field

本发明涉及的是一种组合导航技术，具体地说是一种阵列天线和GPS（全球定位系统）/SINS（捷联惯性导航系统）组合导航技术。 The present invention relates to an integrated navigation technology, specifically an array antenna and GPS (Global Positioning System)/SINS (Strapdown Inertial Navigation System) integrated navigation technology.

背景技术 Background technique

GPS导航系统能够全天候的为其用户提供覆盖全球的精确位置、速度与时间信息，其应用价值越来越高。然而，随着人为干扰技术的提高，卫星只依靠其扩频体制进行抗干扰已经不能满足用户需求。根据ICD-200，商用GPS接收机的抗干扰容限不超过24dB（取决于噪声电平）。即如果干信比大于24dB，商用GPSC/A码接收机就无法保持对信号的跟踪。试验表明，功率为1W的干扰机可以使85公里以内的C/A码接收机无法工作。另外，现有GPS信号是在众所周知的频率上发射的，其调制特征广为人知，信噪比又比较低，因而很容易进行干扰或欺骗。因此，研究GPS抗干扰技术已成为热点和重点。 The GPS navigation system can provide its users with accurate position, speed and time information covering the whole world around the clock, and its application value is getting higher and higher. However, with the improvement of man-made jamming technology, the satellite can no longer meet the needs of users only relying on its spread spectrum system for anti-jamming. According to ICD-200, commercial GPS receivers have an anti-jamming tolerance of no more than 24dB (depending on the noise level). That is, if the interference-to-signal ratio is greater than 24dB, the commercial GPSC/A code receiver cannot keep track of the signal. Tests have shown that a jammer with a power of 1W can disable the C/A code receiver within 85 kilometers. In addition, existing GPS signals are transmitted on well-known frequencies with well-known modulation characteristics and relatively low signal-to-noise ratios, making them easy to jam or spoof. Therefore, research on GPS anti-jamming technology has become a hot spot and focus.

目前，基于信号处理的抗干扰方法是研究最为活跃的领域，常见的技术包括时域滤波技术和空域滤波技术。自适应时域滤波技术是一种窄带干扰抑制方法，它将接收到的有用信号、干扰和噪声通过自适应算法使所需的代价函数最小化来移除干扰信号。时域滤波在用于有界干扰源时是十分稳定的，这是因为它能同时提供复杂带阻滤波器准则，并且它被看作是GPS接收机前后处理的嵌入部分或是GPS接收机前的独立嵌入部分。在对尺寸影响很小的情况下，这种技术对窄带干扰的抑制大于30dB。对于消除窄带干扰，时域滤波技术可用于复杂窄带和连续波干扰源，但这又会受到残余计算带宽的限制，这种残余计算带宽会妨碍有效的GPS信号处理。空域滤波技术是通过自适应天线阵来实现，它可以有效抑制相干干扰和宽带干扰。自适应天线与时域自适应滤波有点类似，也要使得某个代价函数最小化，它有个很大的缺陷就是如果有用信号和干扰信号的入射方向彼此靠近时，自适应天线在抑制干扰信号的同时对有用信号也会造成影响。它的实现需要消耗比较大的功率、更多的花费，而且需要一个比较大的工作平台。空域滤波主要有两种类型:零陷和波束形成。零陷技术对有用信号信息量的需求是最小的，它比波束形成技术更容易实现。零陷技术最小化信号的输出功率，它可以分为窄带和宽带两种情况。波束形成技术要求比较多的有关有用信号的信息，而且实现起来也比较复杂。波束形成技术最大化信号的输出信噪比，波束形成技术也分为窄带和宽带两种情况。 At present, the anti-jamming method based on signal processing is the most active field of research, and common techniques include time domain filtering technology and spatial domain filtering technology. Adaptive time-domain filtering technology is a narrow-band interference suppression method, which uses an adaptive algorithm to minimize the required cost function of the received useful signal, interference and noise to remove the interference signal. Time-domain filtering is very stable when used with bounded interferers, since it can simultaneously provide complex band-stop filter criteria, and it is considered as an embedded part of GPS receiver pre- and post-processing or as a GPS receiver pre- and post-processing The independent embedded part of . This technique suppresses narrowband interference by more than 30dB with little impact on size. For narrowband interference removal, time-domain filtering techniques can be used for complex narrowband and continuous wave interferers, but this is again limited by the residual computational bandwidth that prevents effective GPS signal processing. The spatial filtering technology is realized through an adaptive antenna array, which can effectively suppress coherent interference and broadband interference. The adaptive antenna is somewhat similar to the time-domain adaptive filtering, and it also needs to minimize a certain cost function. It has a big defect that if the incident directions of the useful signal and the interference signal are close to each other, the adaptive antenna suppresses the interference signal At the same time, it will also affect the useful signal. Its implementation requires relatively large power consumption, more costs, and a relatively large working platform. There are two main types of spatial filtering: nulling and beamforming. Null notching requires minimal amount of useful signal information, and it is easier to implement than beamforming. The null trap technology minimizes the output power of the signal, which can be divided into narrowband and wideband situations. Beamforming technology requires more information about useful signals, and it is more complicated to implement. Beamforming technology maximizes the output signal-to-noise ratio of the signal, and beamforming technology is also divided into narrowband and broadband situations.

单纯的时域滤波和单纯的空域各自具有优缺点，但这两种处理方法的优劣恰好可以互补，可以将二者结合应用，进而形成一种联合抗干扰技术，即时空二维联合处理（STAP）抗干扰技术。空时二维联合处理的空域处理能力比单纯的空域滤波要更强。当需要对宽带信号进行零陷的时候，由于信号带宽不可忽略，会在空间中造成所谓“发散”现象。一般的空域滤波手段都需要对这种宽带引发的空间“发散”予以特别的考虑，这也是单纯的空域方法处理宽带信号难以获得满意效果的原因。空时二维联合处理将频率信息和空域信息结合使用，自然地会在空时平面上将发散的各信号分量完整的表现出来，从而可以最大限度地对干扰进行消除和抑制。并且，STAP技术还具有内在的波束形成、内在的信号均衡和内在的抗多径干扰能力等诸多优点，因而可以在增强信号的同时实现抗干扰处理。 Simple time-domain filtering and pure air-domain filtering have their own advantages and disadvantages, but the advantages and disadvantages of these two processing methods can be complementary, and the two can be combined to form a joint anti-jamming technology, that is, spatial-temporal two-dimensional joint processing ( STAP) anti-jamming technology. The spatial processing ability of space-time two-dimensional joint processing is stronger than that of pure spatial filtering. When it is necessary to zero-notch a wideband signal, the so-called "divergence" phenomenon will be caused in space because the signal bandwidth cannot be ignored. The general spatial domain filtering methods need to give special consideration to the spatial "divergence" caused by this broadband, which is also the reason why it is difficult to obtain satisfactory results when processing broadband signals with pure spatial domain methods. Space-time two-dimensional joint processing uses frequency information and space domain information together, and will naturally completely display the divergent signal components on the space-time plane, so that interference can be eliminated and suppressed to the greatest extent. Moreover, STAP technology also has many advantages such as internal beamforming, internal signal equalization and internal anti-multipath interference capability, so it can realize anti-interference processing while enhancing the signal.

另外，GPS/SINS组合导航也能够提高GPS接收机的抗干扰能力。接收机的环路带宽需要在抗干扰能力和动态跟踪性能之间折中设计，即接收机为了提高自身抑制干扰和噪声的能力，需要将环路带宽减小，而为了跟踪载体的动态性能又需要将环路带宽增大。而引入惯性信息后，SINS可以精确估计出载体的速度，计算出载体相对于卫星的多普勒频移，从而减小接收机的环路带宽，增加接收机的抗干扰能力。 In addition, GPS/SINS integrated navigation can also improve the anti-jamming capability of GPS receivers. The loop bandwidth of the receiver needs to be designed as a compromise between anti-interference ability and dynamic tracking performance, that is, in order to improve the receiver’s ability to suppress interference and noise, the loop bandwidth needs to be reduced, and in order to track the dynamic performance of the carrier The loop bandwidth needs to be increased. After introducing inertial information, SINS can accurately estimate the velocity of the carrier and calculate the Doppler frequency shift of the carrier relative to the satellite, thereby reducing the loop bandwidth of the receiver and increasing the anti-interference ability of the receiver.

然而，以上传统的抗干扰方法都存在一定的缺陷。单纯的时域处理方法不能够抑制宽带干扰，空域滤波中的零陷方法不能在卫星信号方向提供增益，卫星信号也有可能被抑制。空域滤波中的波束形成方法和空时自适应处理往往需要卫星到达接收机的角度。GPS/SINS组合抗干扰方法需要GPS接收机能够正确输出伪距、伪距率或位置、速度信息。 However, the above traditional anti-jamming methods all have certain defects. The pure time-domain processing method cannot suppress broadband interference, and the nulling method in air-domain filtering cannot provide gain in the direction of satellite signals, and satellite signals may also be suppressed. Beamforming methods and space-time adaptive processing in spatial filtering often require the angle at which the satellite arrives at the receiver. The GPS/SINS combined anti-jamming method requires that the GPS receiver can correctly output pseudo-range, pseudo-range rate or position and velocity information.

发明内容 Contents of the invention

本发明的目的在于提供一种抑制干扰能力强的基于阵列天线和GPS/SINS的联合抗干扰方法。 The purpose of the present invention is to provide a joint anti-interference method based on array antenna and GPS/SINS with strong interference suppression capability.

本发明的目的是这样实现的： The purpose of the present invention is achieved like this:

初始化载体的位置和姿态之后，建立GPS/SINS组合导航状态方程和量测方程；GPS/SINS组合导航实时提供载体的位置和姿态，同时根据卫星星历信息计算出当前卫星的位置，从而获得卫星到载体之间的导向矢量，即卫星到达载体的方向角和俯仰角；然后，所述导向矢量作为多约束最小方差空时自适应处理(MCMV-STAP)算法的先验信息，在空域、时域同时抑制宽带干扰和窄带干扰。具体包括如下步骤： After initializing the position and attitude of the carrier, establish the GPS/SINS integrated navigation state equation and measurement equation; the GPS/SINS integrated navigation provides the position and attitude of the carrier in real time, and calculates the current satellite position according to the satellite ephemeris information, thereby obtaining the satellite The steering vector between the carrier and the carrier, that is, the azimuth and pitch angle of the satellite arriving at the carrier; then, the steering vector is used as the prior information of the multi-constrained minimum variance space-time adaptive processing (MCMV-STAP) algorithm, in the airspace, time The domain suppresses both wideband and narrowband interference. Specifically include the following steps:

步骤1：初始化载体位置、速度和姿态信息 Step 1: Initialize carrier position, velocity and attitude information

在大地坐标系中，设定载体初始时刻的坐标：纬度L、经度λ和高度h；初始化载体在东北天坐标系中的速度：东向速度V_E、北向速度V_N和天向速度V_U；初始化载体的姿态角，包括俯仰角θ、横滚角γ和方位角ψ；然后设定载体的飞行路径，可以设定为静止、直线运动或圆周运动等轨迹。从而获得理想陀螺仪和加速度计的输出，f_E、f_N和f_U分别表示为加速度计在东向、北向和天向的输出。 In the geodetic coordinate system, set the coordinates of the carrier at the initial moment: latitude L, longitude λ, and height h; initialize the velocity of the carrier in the northeast sky coordinate system: eastward velocity V _E , northward velocity V _N and skyward velocity V _U ; Initialize the attitude angle of the carrier, including pitch angle θ, roll angle γ and azimuth ψ; then set the flight path of the carrier, which can be set as static, linear motion or circular motion. In order to obtain the output of the ideal gyroscope and accelerometer, f _E , f _N and f _U are respectively expressed as the output of the accelerometer in the east, north and sky directions.

步骤2：初始化导航滤波器的参数信息 Step 2: Initialize the parameter information of the navigation filter

在GPS/SINS组合导航滤波器中，采用反馈校正方式。导航滤波器的状态量为姿态角误差、速度误差、位置误差、陀螺仪误差和加速度计误差。φ_E、φ_N和φ_U分别表示为载体俯仰角误差、横滚角误差和方位角误差，δV_E、δV_N和δV_U分别表示为载体的东向速度误差、北向速度误差和天向速度误差，δL、δλ和δh分别表示为载体的纬度误差、经度误差和高度误差。ε_E、ε_N和ε_U分别表示为陀螺仪在东向、北向和天向的漂移。和分别表示为加速度计在东向、北向和天向的输出误差。 In the GPS/SINS integrated navigation filter, the feedback correction method is adopted. The state quantities of the navigation filter are attitude angle error, velocity error, position error, gyroscope error and accelerometer error. φ _E , φ _N and φ _U are respectively the carrier pitch angle error, roll angle error and azimuth error, and δV _E , δV _N and δV _U are respectively the carrier’s eastward velocity error, northward velocity error and skyward velocity Errors, δL, δλ, and δh are respectively the latitude error, longitude error, and height error of the carrier. ε _E , ε _N and ε _U represent the drift of the gyroscope in the east, north and sky directions, respectively. and Expressed as the output errors of the accelerometer in the east, north and celestial directions, respectively.

步骤3：计算陀螺仪和加速度计的误差变化率 Step 3: Calculate the Rate of Change of Error for the Gyroscope and Accelerometer

步骤4：根据步骤1-3中的参数，计算载体姿态角误差变化率、速度误差变化率和位置误差变化率 Step 4: According to the parameters in steps 1-3, calculate the rate of change of carrier attitude angle error, rate of change of velocity error and rate of change of position error

步骤5：引入GPS伪距量测信息，采用反馈校正方式对SINS输出信息进行校正，获得当前准确的位置和姿态。 Step 5: Introduce the GPS pseudo-range measurement information, and use the feedback correction method to correct the SINS output information to obtain the current accurate position and attitude.

根据GPS接收机跟踪环路所测量的码相位误差计算出卫星到载体的之间的伪距，伪距作为量测信息来更新导航滤波器的状态量，从而预测出当前SINS所量测状态量的误差。然后采用反馈校正方式对SINS输出信息进行校正，获得当前准确的位置和姿态。 Calculate the pseudo-range between the satellite and the carrier according to the code phase error measured by the tracking loop of the GPS receiver. The pseudo-range is used as measurement information to update the state quantity of the navigation filter, thereby predicting the state quantity measured by the current SINS error. Then use the feedback correction method to correct the SINS output information to obtain the current accurate position and attitude.

步骤6：计算载体在地心地固坐标系中的坐标 Step 6: Calculate the coordinates of the carrier in the earth-centered earth-fixed coordinate system

由步骤5计算出的载体位置位于大地坐标系中，即(L,λ,h)，将其转化为地心地固坐标系中，其位置坐标可以表示为 The carrier position calculated by step 5 is located in the geodetic coordinate system, namely (L, λ, h), which is converted into the earth-centered and ground-fixed coordinate system, and its position coordinates It can be expressed as

${X x}_{p p}^{e e} = = {[[((R R + + h h)) cos cos L L cos cos λ λ,, ((R R + + h h)) cos cos L L sin sin λ λ,, ((R R + + h h)) sin sin L L]]}^{T T} - - - - - - ((11))$

式(1)中，R为地球半径。 In formula (1), R is the radius of the earth.

步骤7：计算载体坐标系(b系)到导航坐标系(n系)的转换矩阵、导航坐标系到地心地固坐标系(e系)的转换矩阵 Step 7: Calculate the conversion matrix from the carrier coordinate system (b system) to the navigation coordinate system (n system), and the conversion matrix from the navigation coordinate system to the earth-centered ground-fixed coordinate system (e system)

由步骤5计算出的载体姿态角，可以计算出载体坐标系到导航坐标系的转换矩阵为 From the carrier attitude angle calculated in step 5, the transformation matrix from the carrier coordinate system to the navigation coordinate system can be calculated for

${C C}_{b b}^{n no} = = [\begin{matrix} cos cos γ γ cos cos ψ ψ - - sin sin ψ ψ sin sin θ θ sin sin γ γ & - - sin sin ψ ψ cos cos θ θ & sin sin γ γ cos cos ψ ψ + + cos cos γ γ sin sin θ θ sin sin ψ ψ \\ cos cos γ γ sin sin ψ ψ + + sin sin γ γ sin sin θ θ cos cos ψ ψ & cos cos ψ ψ cos cos θ θ & sin sin γ γ sin sin ψ ψ - - cos cos γ γ sin sin θ θ cos cos ψ ψ \\ - - sin sin γ γ cos cos θ θ & sin sin θ θ & cos cos θγ θγ cos cos \end{matrix}] - - - - - - ((22))$

由步骤5计算出的载体位置，可以计算出导航坐标系到地心地固坐标系的转换矩阵为 From the position of the carrier calculated in step 5, the transformation matrix from the navigation coordinate system to the earth-centered earth-fixed coordinate system can be calculated for

${C C}_{n no}^{e e} = = [\begin{matrix} - - sin sin λ λ & - - sin sin L L cos cos λ λ & cos cos L L cos cos λ λ \\ cos cos λ λ & - - sin sin L L sin sin λ λ & cos cos L L sin sin λ λ \\ 00 & cos cos L L & sin sin L L \end{matrix}] - - - - - - ((33))$

步骤8：计算载体到卫星之间的导向矢量 Step 8: Calculate the steering vector between the vehicle and the satellite

由步骤7中的(2)、(3)可以求得载体坐标系到地心地固坐标系的转化矩阵 From (2) and (3) in step 7, the transformation matrix from the carrier coordinate system to the earth-centered and ground-fixed coordinate system can be obtained

${C C}_{b b}^{e e} = = {C C}_{n no}^{e e} {C C}_{b b}^{n no} - - - - - - ((44))$

利用卫星星历解算出卫星在地心地固坐标系中的位置并结合公式(1)(4),可以求得载体到卫星之间的导向矢量 Calculating the Satellite's Position in the Earth-centered Earth-Fixed Coordinate System Using Satellite Ephemeris Combined with formulas (1) (4), the steering vector between the carrier and the satellite can be obtained

${\overset{&RightArrow; &Right Arrow;}{r r}}^{b b} = = {(({C C}_{b b}^{e e}))}^{T T} (({X x}_{s the s}^{e e} - - {X x}_{p p}^{e e})) - - - - - - ((55))$

步骤9：在载体坐标系中，计算卫星到达天线的方位角和俯仰角 Step 9: In the carrier coordinate system, calculate the azimuth and elevation angles at which the satellite arrives at the antenna

若将在载体坐标系中的坐标定义为则卫星到达天线的方位角α和俯仰角β分别可以表示为 If will The coordinates in the carrier coordinate system are defined as Then the azimuth α and elevation angle β of the satellite arriving at the antenna can be expressed as

$α α = = arctan arctan ((\overset{^^}{x x},, \overset{^^}{y the y})) - - - - - - ((66))$

$β β = = arctan arctan ((\overset{^^}{z z},, \sqrt{{\overset{^^}{x x}}^{22} + + {\overset{^^}{y the y}}^{22}})) - - - - - - ((77))$

步骤10：设计圆形阵列天线结构 Step 10: Design the Circular Array Antenna Structure

为了能够在方位角和俯仰角方向同时控制波束指向可视卫星方向，本发明采用圆形阵列天线结构。6个阵元均匀分布于圆周上的圆阵，令圆半径为r，则圆周上相邻阵元间的间隔也为r。阵元间隔的选取，与时域采样间隔满足奈奎斯特定理一样，空域采样间隔d应小于卫星载波波长λ的1/2。由卫星信号频率f=1575.42×10⁶MHz，所以阵元间距为 In order to be able to simultaneously control the beam to point to the visible satellite direction in the direction of azimuth angle and elevation angle, the present invention adopts a circular array antenna structure. A circular array with 6 array elements evenly distributed on the circumference, if the radius of the circle is r, then the interval between adjacent array elements on the circumference is also r. The selection of the array element interval is the same as the sampling interval in the time domain satisfies the Nyquist theorem, and the sampling interval d in the air domain should be less than 1/2 of the satellite carrier wavelength λ. The satellite signal frequency f=1575.42×10 ⁶ MHz, so the array element spacing is

$d d \leq \leq \frac{λ λ}{22} = = \frac{11}{22} \cdot &Center Dot; \frac{c c}{{f f}_{l l 11}} = = \frac{11}{22} \cdot &Center Dot; \frac{33 \times \times 1010^{88}}{1575.42 1575.42 \times \times 1010^{66}} = = \frac{11}{22} \cdot &Center Dot; 0.19 0.19 = = 0.095 0.095 m m = = 9.5 9.5 cm cm - - - - - - ((88))$

式(8)中，c为光速。 In formula (8), c is the speed of light.

为了使主波束宽度越窄，副瓣越低，分辨率高，需要使阵元间隔尽量的大，所以取圆半径r=d=9.5cm，整个阵列天线的直径约为19cm。 In order to make the main beam width narrower, the sidelobe lower, and the resolution higher, it is necessary to make the array element spacing as large as possible, so take the circle radius r=d=9.5cm, and the diameter of the entire array antenna is about 19cm.

步骤11：计算卫星信号到达天线各阵元之间的时间延迟 Step 11: Calculate the time delay between the arrival of the satellite signal at each element of the antenna

根据步骤9中的(6)(7)，可以将α和β表示单位矢量 According to (6)(7) in step 9, α and β can be expressed as unit vectors

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

因此，卫星信号到达第i个天线阵元与到第一个参考阵元之间的时间差τ_i可以表示为 Therefore, the time difference τ _i between the arrival of the satellite signal at the i-th antenna element and the first reference element can be expressed as

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

式(10)中，M代表天线阵元数目。 In formula (10), M represents the number of antenna elements.

步骤12：建立阵列天线接收信号模型 Step 12: Build an array antenna receiving signal model

用户一般可以接收4颗以上的卫星信号，因此，需要形成多个波束指向对应的卫星。假设阵列天线接收到了P个卫星信号，Q个干扰信号，则天线接收到的信号模型可以表示为 Generally, users can receive signals from more than 4 satellites. Therefore, multiple beams need to be formed to point to corresponding satellites. Assuming that the array antenna receives P satellite signals and Q interference signals, the signal model received by the antenna can be expressed as

$U u ((t t)) = = {Σ Σ}_{k k = = 11}^{P P} a a (({α α}_{k k},, {β β}_{k k},, {T T}_{s the s})) {s the s}_{k k} ((t t)) + + {Σ Σ}_{k k = = P P + + 11}^{P P + + Q Q} a a (({α α}_{k k},, {β β}_{k k},, {T T}_{s the s})) {j j}_{k k - - P P} ((t t)) + + n no ((t t)) - - - - - - ((1111))$

在公式(11)中，s(t)和j(t)分别表示接收到的卫星信号和干扰信号，a(α_k,β_k,T_s)是第k个目标信号(卫星信号或干扰信号)的导向矢量。T_s是时域延迟线间隔，α_k和β_k分别表示为第k个目标信号到达阵列天线的方位角和俯仰角。n(t)表示高斯白噪声，它的功率谱密度表示为N₀/2。 In formula (11), s(t) and j(t) represent the received satellite signal and interference signal respectively, and a(α _k , β _k , T _s ) is the kth target signal (satellite signal or interference signal ) of the steering vector. T _s is the delay line interval in the time domain, and α _k and β _k represent the azimuth and elevation angles of the kth target signal arriving at the array antenna, respectively. n(t) represents Gaussian white noise, and its power spectral density is expressed as N ₀ /2.

在公式(11)中，a(α_k,β_k,T_s)表示空时二维目标矢量，即时间矢量a_s(T_s)和空间方向矢量a_s(α_k,β_k)的克罗奈克积（KroneckerProduct），并可以表示为： In formula (11), a(α _k ,β _k ,T _s ) represents the space-time two-dimensional target vector, that is, the gram of the time vector a _s (T _s ) and the space direction vector a _s (α _k ,β _k ) Ronecker Product (KroneckerProduct), and can be expressed as:

$a a (({α α}_{k k},, {β β}_{k k},, {T T}_{s the s})) = = a a (({α α}_{k k},, {β β}_{k k})) &CircleTimes; &CircleTimes; {a a}_{s the s} (({T T}_{s the s})) - - - - - - ((1212))$

令每个射频通道的时间延迟单元数目为N，则 Let the number of time delay units of each RF channel be N, then

${a a}_{s the s} (({α α}_{k k},, {β β}_{k k})) = = [\begin{matrix} 11 & {e e}^{- - j j 22 πf πf {τ τ}_{11}} & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & {e e}^{- - j j 22 πf πf {τ τ}_{M m - - 11}} \end{matrix}] - - - - - - ((1313))$

${a a}_{t t} = = [\begin{matrix} 11 & {e e}^{- - j j 22 πfT πfT} & \cdot \cdot \cdot &Center Dot; \cdot &Center Dot; & {e e}^{- - j j 22 πf πf ((N N - - 11)) {T T}_{s the s}} \end{matrix}] - - - - - - ((1414))$

式(13)中，f表示卫星信号频率，τ_i(i=1,2,…M-1)由式(10)给出。式(14)中，T_s表示延迟单元的时间间隔，其值应该小于信号带宽。 In Equation (13), f represents the satellite signal frequency, and τ _i (i=1,2,...M-1) is given by Equation (10). In formula (14), T _s represents the time interval of the delay unit, and its value should be smaller than the signal bandwidth.

步骤13：计算天线阵列权值矢量 Step 13: Compute Antenna Array Weight Vector

采用多约束最小方差空时自适应处理(MCMV-STAP)算法对可视卫星信号进行约束，然后使得阵列输出功率最小，从而保护卫星信号的同时抑制干扰信号。该算法需要求得阵列的最优权值，表达如下： The multi-constrained minimum variance space-time adaptive processing (MCMV-STAP) algorithm is used to constrain the visible satellite signals, and then minimize the output power of the array, so as to protect the satellite signals and suppress the interference signals at the same time. The algorithm needs to find the optimal weight of the array, expressed as follows:

$\{\begin{matrix} ω ω = = arg arg min min {ω ω}^{H h} {R R}_{U u} ω ω \\ subjectto subject to : : {ω ω}^{H h} A A = = {F f}^{T T} \end{matrix} - - - - - - ((1515))$

式(15)中，ω表示阵列的权值矢量，H表示矩阵的共轭转置，R_U是输入信号的阵列协方差矩阵，可以表示为 In formula (15), ω represents the weight vector of the array, H represents the conjugate transpose of the matrix, and R _U is the array covariance matrix of the input signal, which can be expressed as

R_U=E{UU^H}(16) R _U =E{UU ^H }(16)

式(15)中，A表示卫星信号的约束矩阵，F表示与A相对应的约束矢量，A和F分别表示为 In formula (15), A represents the constraint matrix of the satellite signal, F represents the constraint vector corresponding to A, and A and F are expressed as

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

f^T=[1,1…1]_1×p(18)采用拉格朗日乘子法，阵列权值矢量ω可以表示为 f ^T =[1,1…1] _1×p (18) Using the Lagrange multiplier method, the array weight vector ω can be expressed as

$ω ω = = {R R}_{U u}^{- - 11} A A {(({A A}^{H h} {R R}_{U u}^{- - 11} A A))}^{- - 11} f f - - - - - - ((1919))$

在获得阵列权值之后，可以得到阵列的输出表达是为 After obtaining the array weights, the output expression of the array can be obtained as

y(t)=ω^HU(t)。 y(t)=ω ^H U(t).

传统的空时自适应处理方法，通常采用盲自适应波束形成算法，即不需要知道卫星到达阵列天线的方位角和俯仰角，只是简单的将某一天线阵元作为参考阵元进行权值约束。然而，由于未在卫星信号方向进行约束，不能在卫星信号方向形成波束。因此，该类盲波束形成算法在抑制干扰信号的同时也将削弱期望的卫星信号，卫星信号功率将被减弱，不能够得到信号的最大信干噪比。 The traditional space-time adaptive processing method usually adopts the blind adaptive beamforming algorithm, that is, it does not need to know the azimuth and elevation angles of the satellite arrival array antenna, but simply uses a certain antenna array element as a reference array element for weight constraints . However, since there is no constraint in the direction of the satellite signal, the beam cannot be formed in the direction of the satellite signal. Therefore, this type of blind beamforming algorithm will also weaken the desired satellite signal while suppressing the interference signal, the power of the satellite signal will be weakened, and the maximum signal-to-interference-noise ratio of the signal cannot be obtained.

本发明主要利用GPS/SINS组合导航系统提供的载体位置和姿态，GPS星历提供的卫星位置，从而获得卫星到载体之间的导向矢量，即卫星到达载体的方向角和俯仰角。然后，该导向矢量作为多约束最小方差空时自适应处理(MCMV-STAP)算法的先验信息，在空域、时域同时抑制宽带干扰和窄带干扰。 The present invention mainly utilizes the position and attitude of the carrier provided by the GPS/SINS integrated navigation system, and the satellite position provided by the GPS ephemeris, thereby obtaining the guiding vector between the satellite and the carrier, that is, the direction angle and pitch angle at which the satellite arrives at the carrier. Then, the steering vector is used as the prior information of the Multi-Constrained Minimum Variance Space-Time Adaptive Processing (MCMV-STAP) algorithm to suppress wideband interference and narrowband interference in both space and time domains.

本发明所述的计算卫星到载体之间的导向矢量方法，其特征在于GPS/SINS所解算出载体的位置位于大地坐标系中，姿态位于东北天坐标系中，利用卫星星历解算的卫星位置在地心地固坐标系中，而波束形成所需的导向矢量是应该在载体坐标系下提供的。本发明将所求得的先验信息(卫星的位置、载体的位置和姿态)经过坐标变化后，在载体坐标系上获得卫星到达载体的方向角和俯仰角。 The method for calculating the steering vector between the satellite and the carrier of the present invention is characterized in that the position of the carrier calculated by GPS/SINS is located in the geodetic coordinate system, the attitude is located in the northeast sky coordinate system, and the satellite ephemeris is used to solve the problem. The position is in the earth-centered ground-fixed coordinate system, and the steering vector required for beamforming should be provided in the vehicle coordinate system. The present invention obtains the azimuth angle and pitch angle at which the satellite arrives at the carrier on the carrier coordinate system after the obtained prior information (position of the satellite, position and attitude of the carrier) is changed through coordinates.

本发明采用阵列天线与GPS/SINS相结合的方式来提高导航接收机的抗干扰能力。在阵列天线中，采用多约束最小方差准则的空时自适应波束形成算法(MCMV-STAP)，该算法能够在多颗可视卫星方向形成波束，在干扰方向形成零陷，从而增强卫星信号的同时抑制干扰信号。为了获得卫星信号的最大功率，该算法需要知道卫星信号到达载体的方位角和俯仰角。因此，本发明采用圆形阵列天线结构，引入GPS/SINS组合导航为阵列天线的波束形成提供载体的位置和姿态，采用卫星星历提供卫星的位置，从而为波束形成提供先验信息，即可视卫星到达载体的方位角和俯仰角。 The invention adopts the method of combining the array antenna and GPS/SINS to improve the anti-jamming capability of the navigation receiver. In the array antenna, the multi-constrained minimum variance criterion space-time adaptive beamforming algorithm (MCMV-STAP) is used. This algorithm can form beams in the direction of multiple visible satellites and form nulls in the interference direction, thereby enhancing the satellite signal. Simultaneously suppresses interfering signals. In order to obtain the maximum power of the satellite signal, the algorithm needs to know the azimuth and elevation angle at which the satellite signal arrives at the carrier. Therefore, the present invention adopts a circular array antenna structure, introduces GPS/SINS integrated navigation to provide the position and attitude of the carrier for the beam forming of the array antenna, and uses the satellite ephemeris to provide the position of the satellite, thereby providing prior information for the beam forming, that is, The azimuth and elevation angles at which the satellite arrives at the carrier.

附图说明 Description of drawings

图1是基于阵列天线和GPS/SINS的联合抗干扰方法的实施流程图； Fig. 1 is the implementation flowchart of the joint anti-jamming method based on array antenna and GPS/SINS;

图2是阵列天线、东北天坐标系和载体坐标系的示意图； Fig. 2 is a schematic diagram of the array antenna, the northeast sky coordinate system and the carrier coordinate system;

图3是多约束最小方差空时二维自适应处理(MCMV-STAP)的结构图。 Fig. 3 is a structural diagram of the multi-constrained minimum variance space-time two-dimensional adaptive processing (MCMV-STAP).

具体实施方式 Detailed ways

下面结合附图对本发明做进一步的阐述。 The present invention will be further elaborated below in conjunction with the accompanying drawings.

如附图1所示，需要先初始化载体的位置与姿态。在大地坐标系中，设定载体初始时刻的坐标：纬度L=45°、经度λ=126°和高度h=200m。初始化载体的姿态角，包括俯仰角θ=0°、横滚角γ=0°和方位角ψ=45°。然后设定载体的飞行路径，可以设定为静止、直线运动或圆周运动等轨迹。根据要求的载体飞行路径设定载体的东向速度V_E、北向速度V_N和天向速度V_U。从而获得理想陀螺仪和加速度计的输出，f_E、f_N和f_U分别表示为加速度计在东向、北向和天向的输出。 As shown in Figure 1, the position and attitude of the carrier need to be initialized first. In the geodetic coordinate system, set the coordinates of the carrier at the initial moment: latitude L=45°, longitude λ=126° and height h=200m. Initialize the attitude angle of the carrier, including pitch angle θ=0°, roll angle γ=0° and azimuth angle ψ=45°. Then set the flight path of the carrier, which can be set as static, linear motion or circular motion and other trajectories. Set the eastward velocity V _E , northward velocity V _N and skyward velocity V _U of the carrier according to the required flight path of the carrier. In order to obtain the output of the ideal gyroscope and accelerometer, f _E , f _N and f _U are respectively expressed as the output of the accelerometer in the east, north and sky directions.

在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。加速度计在东向、北向和天向的输出误差分别可以表示为和 In the GPS/SINS integrated navigation filter, the feedback correction method is adopted. The state quantities of the navigation filter are attitude angle error, velocity error, position error, gyroscope error and accelerometer error. At the initial moment, the carrier pitch angle error φ _E = 0.1°, roll angle error φ _N = 0.1° and azimuth angle error φ _U = 1°. The carrier's eastward velocity error δV _E =0.2m/s, northward velocity error δV _N =0.2m/s and skyward velocity error δV _U =0.5m/s. The carrier's latitude error δL=0°, longitude error δλ=0° and height error δh=20m. Gyroscope constant drift ε _b =0.1°/h, white noise ε _g =0.05°/h, accelerometer constant error White noise w _a =5×10 ^-4 g, g is the gravitational acceleration. The drifts of the gyroscope error ε in the east, north and sky directions are denoted as ε _E , ε _N and ε _U , respectively. The output errors of the accelerometer in the east direction, north direction and sky direction can be expressed as and

陀螺仪的常值漂移可以以用随机常数描述为 The constant drift of a gyroscope can be described by a random constant as

${\overset{\cdot \cdot}{ϵ ϵ}}_{b b} = = 00 - - - - - - ((11))$

陀螺仪的噪声ε_g来可以用狄拉克函数描述。因此，陀螺误差变化率可以表示为 The noise ε _g of the gyroscope can be described by the Dirac function. Therefore, the rate of change of gyro error It can be expressed as

$\overset{\cdot \cdot}{ϵ ϵ} = = {\overset{\cdot \cdot}{ϵ ϵ}}_{b b} + + {ϵ ϵ}_{g g} - - - - - - ((22))$

对于加速度计误差变化率可以将其考虑为一阶马尔柯夫过程，误差模型取为 For the accelerometer error rate of change It can be considered as a first-order Markov process, and the error model is taken as

$\overset{\cdot \cdot}{&dtri; &dtri;} = = - - \frac{11}{{T T}_{a a}} &dtri; &dtri; + + {w w}_{a a} - - - - - - ((33))$

其中，T_a表示相关时间，w_a为白噪声过程。 Among them, T _a represents the correlation time, and w _a is the white noise process.

1、计算载体姿态角误差变化率 1. Calculate the rate of change of carrier attitude angle error

和分别表示为载体俯仰角误差、横滚角误差和方位角误差的变化率，载体姿态角误差变化率可以表示为 and Expressed as the change rate of carrier pitch angle error, roll angle error and azimuth error respectively, the change rate of carrier attitude angle error can be expressed as

$\begin{matrix} {\overset{\cdot \cdot}{φ φ}}_{E E.} = = - - \frac{{δV δV}_{N N}}{R R + + h h} + + (({ω ω}_{ie ie} sin sin L L + + \frac{{V V}_{E E.}}{R R + + h h} tan the tan L L)) {φ φ}_{N N} - - (({ω ω}_{ie ie} cos cos L L + + \frac{{V V}_{E E.}}{R R + + h h} tan the tan L L)) {φ φ}_{U u} + + {ϵ ϵ}_{E E.} \\ {\overset{\cdot \cdot}{φ φ}}_{N N} = = \frac{{δV δV}_{E E.}}{R R + + h h} - - {ω ω}_{ie ie} sin sin LδL LδL - - (({ω ω}_{ie ie} sin sin L L + + \frac{{V V}_{E E.}}{R R + + h h} tan the tan L L)) {φ φ}_{E E.} - - \frac{{V V}_{E E.}}{R R + + h h} {φ φ}_{U u} + + {ϵ ϵ}_{N N} \\ {\overset{\cdot \cdot}{φ φ}}_{U u} = = \frac{{δV δV}_{E E.}}{R R + + h h} tan the tan L L + + (({ω ω}_{ie ie} cos cos L L + + \frac{{V V}_{E E.}}{R R + + h h} {sec sec}^{22} L L)) δL δ L + + (({ω ω}_{ie ie} cos cos L L + + \frac{{V V}_{E E.}}{R R + + h h})) {φ φ}_{E E.} + + \frac{{V V}_{N N}}{R R + + h h} {φ φ}_{N N} + + {ϵ ϵ}_{U u} \end{matrix} - - - - - - ((44))$

式(4)中，R表示地球半径，ω_ie表示地球旋转角速率。 In formula (4), R represents the radius of the earth, and ω _ie represents the angular rate of the earth's rotation.

2、计算载体速度误差的变化率 2. Calculate the rate of change of carrier velocity error

和分别表示为载体的东向速度误差、北向速度误差和天向速度误差的变化率，载体速度误差变化率可以表示为 and Expressed as the rate of change of the carrier’s eastward velocity error, northward velocity error, and skyward velocity error, respectively, the carrier velocity error rate of change can be expressed as

$\begin{matrix} {\overset{\cdot &Center Dot;}{δV δV}}_{E E.} = = ((\frac{{V V}_{N N}}{R R + + h h} tan the tan L L - - \frac{{V V}_{U u}}{R R + + h h})) {δV δV}_{E E.} + + (({22 ω ω}_{ie ie} sin sin L L + + \frac{{V V}_{E E.}}{R R + + h h} tan the tan L L)) {δV δV}_{N N} \\ - - ((- - {ω ω}_{ie ie} cos cos L L + + \frac{{V V}_{E E.}}{R R + + h h})) {δV δV}_{U u} + + (({22 ω ω}_{ie ie} {V V}_{N N} cos cos L L + + \frac{{V V}_{E E.} {V V}_{N N}}{R R + + h h} {sec sec}^{22} L L + + 22 {ω ω}_{ie ie} {V V}_{U u} sin sin L L)) δL δ L \\ + + {&dtri; &dtri;}_{E E.} + + {f f}_{N N} {φ φ}_{U u} - - {f f}_{U u} {φ φ}_{N N} \\ {\overset{\cdot &Center Dot;}{δV δV}}_{N N} = = - - 22 (({ω ω}_{ie ie} sin sin L L + + \frac{{V V}_{E E.}}{R R + + h h} tan the tan L L)) {δV δV}_{E E.} - - \frac{{V V}_{U u}}{R R + + h h} {δV δV}_{N N} - - \frac{{V V}_{N N}}{R R + + h h} {δV δV}_{U u} \\ - - ((22 {ω ω}_{ie ie} {V V}_{E E.} cos cos L L + + \frac{{V V}_{E E.}^{22} {sec sec}^{22} L L}{R R + + h h})) δL δ L + + {&dtri; &dtri;}_{N N} - - {f f}_{E E.} {φ φ}_{U u} + + {f f}_{U u} {φ φ}_{E E.} \\ {\overset{\cdot \cdot}{δV δV}}_{U u} = = ((22 {ω ω}_{ie ie} cos cos L L + + \frac{{V V}_{E E.}}{R R + + h h})) {δV δV}_{E E.} + + \frac{{22 V V}_{N N}}{R R + + h h} {δV δV}_{N N} - - 22 {ω ω}_{ie ie} {V V}_{E E.} sin sin LδL LδL + + {&dtri; &dtri;}_{U u} + + {f f}_{E E.} {φ φ}_{N N} - - {f f}_{N N} {φ φ}_{E E.} \end{matrix} - - - - - - ((55))$

3、计算载体位置误差的变化率 3. Calculate the rate of change of the carrier position error

和分别表示为载体的纬度误差、经度误差和高度误差的变化率。载体位置误差变化率可以表示为 and Expressed as the rate of change of the carrier's latitude error, longitude error, and altitude error, respectively. The rate of change of carrier position error can be expressed as

$\begin{matrix} δ δ \overset{\cdot \cdot}{L L} = = \frac{{δV δV}_{N N}}{R R + + h h} \\ δ δ \overset{\cdot \cdot}{λ λ} = = \frac{{δV δV}_{E E.}}{R R + + h h} \\ δ δ \overset{\cdot \cdot}{h h} = = δ δ {V V}_{U u} \end{matrix}, sec sec L L + + \frac{{V V}_{E E.}}{R R + + h h} sec sec LtgLδL LtgLδL - - - - - - ((66))$

步骤5：引入GPS伪距量测信息，采用反馈校正方式对SINS输出信息进行校正，获得当前 Step 5: Introduce the GPS pseudo-range measurement information, and use the feedback correction method to correct the SINS output information to obtain the current

准确的位置和姿态。 Accurate position and attitude.

根据GPS接收机跟踪环路所测量的码相位误差计算出卫星到载体的之间的伪距，伪距作为量测信息来更新导航滤波器的状态量，从而预测出当前SINS所量测的位置误差(δL,δλ,δh)和姿态误差(δθ,δγ,δψ)。然后采用反馈校正方式对SINS输出的位置(L,λ,h)和姿态(θ,γ,ψ)进行校正，获得当前准确的位置和姿态。 The pseudo-range between the satellite and the carrier is calculated according to the code phase error measured by the tracking loop of the GPS receiver, and the pseudo-range is used as measurement information to update the state quantity of the navigation filter, thereby predicting the position measured by the current SINS Errors (δL, δλ, δh) and attitude errors (δθ, δγ, δψ). Then the position (L, λ, h) and attitude (θ, γ, ψ) output by SINS are corrected by feedback correction method to obtain the current accurate position and attitude.

L=L+δL L=L+δL

λ=λ+δλ λ=λ+δλ

h=h+δh h=h+δh

θ=θ+δθ(7) θ=θ+δθ(7)

γ=γ+δγ γ=γ+δγ

ψ=ψ+δψ ψ=ψ+δψ

${X x}_{p p}^{e e} = = {[[((R R + + h h)) cos cos L L cos cos λ λ,, ((R R + + h h)) cos cos L L sin sin λ λ,, ((R R + + h h)) sin sin L L]]}^{T T} - - - - - - ((88))$

由步骤5计算出的载体姿态可以计算出载体坐标系到导航坐标系的转换矩阵为 The transformation matrix from the carrier coordinate system to the navigation coordinate system can be calculated from the carrier attitude calculated in step 5 for

${C C}_{b b}^{n no} = = [\begin{matrix} cos cos γ γ cos cos ψ ψ - - sin sin ψ ψ sin sin θ θ sin sin γ γ & - - sin sin ψ ψ cos cos θ θ & sin sin γ γ cos cos ψ ψ + + cos cos γ γ sin sin θ θ sin sin ψ ψ \\ cos cos γ γ sin sin ψ ψ + + sin sin γ γ sin sin θ θ cos cos ψ ψ & cos cos ψ ψ cos cos θ θ & sin sin γ γ sin sin ψ ψ - - cos cos γ γ sin sin θ θ cos cos ψ ψ \\ - - sin sin γ γ cos cos θ θ & sin sin θ θ & cos cos θγ θγ cos cos \end{matrix}] - - - - - - ((99))$

${C C}_{n no}^{e e} = = [\begin{matrix} - - sin sin λ λ & - - sin sin L L cos cos λ λ & cos cos L L cos cos λ λ \\ cos cos λ λ & - - sin sin L L sin sin λ λ & cos cos L L sin sin λ λ \\ 00 & cos cos L L & sin sin L L \end{matrix}] - - - - - - ((1010))$

由步骤7中的(9)、(10)可以求得载体坐标系到地心地固坐标系的转化矩阵 From (9) and (10) in step 7, the transformation matrix from the carrier coordinate system to the earth-centered earth-fixed coordinate system can be obtained

${C C}_{b b}^{e e} = = {C C}_{n no}^{e e} {C C}_{b b}^{n no} - - - - - - ((1111))$

利用卫星星历解算出卫星在地心地固坐标系中的位置并结合公式(8)(11),可以求得载体到卫星之间的导向矢量 Calculating the Satellite's Position in the Earth-centered Earth-Fixed Coordinate System Using Satellite Ephemeris Combined with formulas (8) (11), the steering vector between the carrier and the satellite can be obtained

${\overset{&RightArrow; &Right Arrow;}{r r}}^{b b} = = {(({C C}_{b b}^{e e}))}^{T T} (({X x}_{s the s}^{e e} - - {X x}_{p p}^{e e})) - - - - - - ((1212))$

$α α = = arctan arctan ((\overset{^^}{x x},, \overset{^^}{y the y})) - - - - - - ((1313))$

$β β = = arctan arctan ((\overset{^^}{z z},, \sqrt{{\overset{^^}{x x}}^{22} + + {\overset{^^}{y the y}}^{22}})) - - - - - - ((1414))$

为了能够在方位角和俯仰角方向同时控制波束指向可视卫星方向，本发明采用圆形阵列天线结构，如附图2所示。6个阵元均匀分布于圆周上的圆阵，令圆半径为r，则圆周上相邻阵元间的间隔也为r。阵元间隔的选取，与时域采样间隔满足奈奎斯特定理一样，空域采样间隔d应小于卫星载波波长λ的1/2。由卫星信号频率f=1575.42×10⁶MHz，所以阵元间距为 In order to be able to simultaneously control the beam to point to the visible satellite in both the azimuth and elevation directions, the present invention adopts a circular array antenna structure, as shown in Figure 2. A circular array with 6 array elements evenly distributed on the circumference, if the radius of the circle is r, then the interval between adjacent array elements on the circumference is also r. The selection of the array element interval is the same as the sampling interval in the time domain satisfies the Nyquist theorem, and the sampling interval d in the air domain should be less than 1/2 of the satellite carrier wavelength λ. The satellite signal frequency f=1575.42×10 ⁶ MHz, so the array element spacing is

$d d \leq \leq \frac{λ λ}{22} = = \frac{11}{22} \cdot &Center Dot; \frac{c c}{{f f}_{l l 11}} = = \frac{11}{22} \cdot \cdot \frac{33 \times \times 1010^{88}}{1575.42 1575.42 \times \times 1010^{66}} = = \frac{11}{22} \cdot &Center Dot; 0.19 0.19 = = 0.095 0.095 m m = = 9.5 9.5 cm cm - - - - - - ((1515))$

式(15)中，c为光速。 In formula (15), c is the speed of light.

以天线阵所在平面为xoy平面，以坐标原点为参考，x轴阵元为1号阵元，则6阵元的极坐标分别为：(r,0)、(r,π/3)、(r,2π/3)、(r,π)、(r,4π/3)、(r,5π/3)。 Taking the plane where the antenna array is located as the xoy plane, taking the origin of the coordinates as a reference, and the x-axis array element as the No. 1 array element, the polar coordinates of the 6 array elements are: (r, 0), (r, π/3), ( r,2π/3), (r,π), (r,4π/3), (r,5π/3).

根据步骤9中的(13)(14)，可以将α和β表示单位矢量 According to (13)(14) in step 9, α and β can be expressed as unit vectors

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

因此，卫星信号到达第i个天线阵元与到第一个参考阵元之间的时间差τi可以表示为 Therefore, the time difference τi between the arrival of the satellite signal at the i-th antenna element and the first reference element can be expressed as

τ_i=e^T·(x_i-x₁)/ci=1,2,…M-1(17)式(17)中，M代表天线阵元数目。 τ _i =e ^T ·( _xi -x ₁ )/ci=1,2,...M-1 (17) In formula (17), M represents the number of antenna elements.

$U u ((t t)) = = {Σ Σ}_{k k = = 11}^{P P} a a (({α α}_{k k},, {β β}_{k k},, {T T}_{s the s})) {s the s}_{k k} ((t t)) + + {Σ Σ}_{k k = = P P + + 11}^{P P + + Q Q} a a (({α α}_{k k},, {β β}_{k k},, {T T}_{s the s})) {j j}_{k k - - P P} ((t t)) + + n no ((t t)) - - - - - - ((1818))$

在公式(18)中，s(t)和j(t)分别表示接收到的卫星信号和干扰信号，a(α_k,β_k,T_s)是第k个目标信号(卫星信号或干扰信号)的导向矢量。T_s是时域延迟线间隔，α_k和β_k分别表示为第k个目标信号到达阵列天线的方位角和俯仰角。n(t)表示高斯白噪声，它的功率谱密度表示为N₀/2。 In formula (18), s(t) and j(t) represent the received satellite signal and interference signal respectively, and a(α _k , β _k , T _s ) is the kth target signal (satellite signal or interference signal ) of the steering vector. T _s is the delay line interval in the time domain, and α _k and β _k represent the azimuth and elevation angles of the kth target signal arriving at the array antenna, respectively. n(t) represents Gaussian white noise, and its power spectral density is expressed as N ₀ /2.

在公式(18)中，a(α_k,β_k,T_s)表示空时二维目标矢量，即时间矢量a_s(T_s)和空间方向矢量a_s(α_k,β_k)的克罗奈克积（KroneckerProduct），并可以表示为： In formula (18), a(α _k ,β _k ,T _s ) represents the space-time two-dimensional target vector, that is, the gram of the time vector a _s (T _s ) and the space direction vector a _s (α _k ,β _k ) Ronecker Product (KroneckerProduct), and can be expressed as:

$a a (({α α}_{k k},, {β β}_{k k},, {T T}_{s the s})) = = {a a}_{s the s} (({α α}_{k k},, {β β}_{k k})) &CircleTimes; &CircleTimes; {a a}_{s the s} (({T T}_{s the s})) - - - - - - ((1919))$

${a a}_{s the s} (({α α}_{k k},, {β β}_{k k})) = = [\begin{matrix} 11 & {e e}^{- - j j 22 πf πf {τ τ}_{11}} & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & {e e}^{- - j j 22 πf πf {τ τ}_{M m - - 11}} \end{matrix}] - - - - - - ((2020))$

${a a}_{t t} = = [\begin{matrix} 11 & {e e}^{- - j j 22 πfT πfT} & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & {e e}^{- - j j 22 πf πf ((N N - - 11)) {T T}_{s the s}} \end{matrix}] - - - - - - ((21 twenty one))$

式(20)中，f表示卫星信号频率，τ_i(i=1,2,…M-1)由式(17)给出。式(21)中，T_s表示延迟单元的时间间隔，其值应该小于信号带宽。 In Equation (20), f represents the satellite signal frequency, and τ _i (i=1,2,...M-1) is given by Equation (17). In formula (21), T _s represents the time interval of the delay unit, and its value should be smaller than the signal bandwidth.

附图3为多约束最小方差空时二维自适应处理(MCMV-STAP)的结构图。从每个阵元的通道来看，各级延时构成了FIR滤波器，可在时域去除干扰；从相同的时间延迟节点看，不同的阵元构成了空域的自适应滤波，可以分辨空间干扰源进而形成空域零陷从空域上抑制干扰。而空域的处理也可以进一步利用时域处理后的反馈信息，空时处理也因此具有在空时二维域同时剔除干扰的能力。 Accompanying drawing 3 is the structural diagram of multi-constrained minimum variance space-time two-dimensional adaptive processing (MCMV-STAP). From the perspective of the channel of each array element, the delays at all levels constitute an FIR filter, which can remove interference in the time domain; from the perspective of the same time delay node, different array elements constitute an adaptive filter in the spatial domain, which can distinguish spatial The interference source then forms an airspace null to suppress interference from the airspace. The processing in the space domain can also further utilize the feedback information after the time domain processing, so the space-time processing has the ability to simultaneously eliminate interference in the space-time two-dimensional domain.

采用MCMV-STAP算法对可视卫星信号进行约束，然后使得阵列输出功率最小，从而保护卫星信号的同时抑制干扰信号。该算法需要求得阵列的最优权值，表达如下： The MCMV-STAP algorithm is used to constrain the visible satellite signals, and then the output power of the array is minimized, so as to protect the satellite signals and suppress the interference signals at the same time. The algorithm needs to find the optimal weight of the array, expressed as follows:

$\{\begin{matrix} ω ω = = arg arg min min {ω ω}^{H h} {R R}_{U u} ω ω \\ subjectto subject to : : {ω ω}^{H h} A A = = {F f}^{T T} \end{matrix} - - - - - - ((22 twenty two))$

式(22)中，ω表示阵列的权值矢量，H表示矩阵的共轭转置，R_U是输入信号的阵列协方差矩阵，可以表示为 In formula (22), ω represents the weight vector of the array, H represents the conjugate transpose of the matrix, and R _U is the array covariance matrix of the input signal, which can be expressed as

R_U=E{UU^H}(23) R _U =E{UU ^H }(23)

式(23)中，A表示卫星信号的约束矩阵，F表示与A相对应的约束矢量，A和F分别表示为 In formula (23), A represents the constraint matrix of the satellite signal, F represents the constraint vector corresponding to A, and A and F are expressed as

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

f^T=[1,1…1]_1×p(25)采用拉格朗日乘子法，阵列权值矢量ω可以表示为 f ^T =[1,1…1] _1×p (25) Using the Lagrange multiplier method, the array weight vector ω can be expressed as

$ω ω = = {R R}_{U u}^{- - 11} A A {(({A A}^{H h} {R R}_{U u}^{- - 11} A A))}^{- - 11} f f - - - - - - ((2626))$

y(t)=ω^HU(t)。(27)。 y(t)=ω ^H U(t). (27).

Claims

1., based on an associating anti-interference method of array antenna and GPS/SINS, after the position of initialization carrier and attitude, set up GPS/SINS integrated navigation state equation and measurement equation; GPS/SINS integrated navigation provides position and the attitude of carrier in real time, calculates the position of present satellites according to satellite ephemeris information simultaneously, thus obtains the steering vector between satellite to carrier, and namely satellite arrives deflection and the angle of pitch of carrier; Then, described steering vector as the prior imformation of multiple constraint minimum variance space-time adaptive Processing Algorithm, in spatial domain, time domain suppress simultaneously broadband interference and arrowband interference; It is characterized in that specifically comprising the steps:

Step 1: initialization carrier positions, speed and attitude information;

In earth coordinates, the coordinate of setting carrier initial time: latitude L, longitude λ and height h; The speed of initialization carrier northeastward in sky coordinate system: east orientation speed V _e, north orientation speed V _nwith sky to speed V _u; The attitude angle of initialization carrier, comprises pitching angle theta, roll angle γ and position angle ψ; Then the flight path setting carrier is the one in static, rectilinear motion or circular motion track; Thus obtain the output of gyroscope and accelerometer, f _e, f _nand f _ube expressed as accelerometer in east orientation, north orientation and sky to output;

Step 2: the parameter information of initialization Navigation Filter;

In GPS/SINS integrated navigation wave filter, adopt feedback compensation mode; The quantity of state of Navigation Filter is attitude error, velocity error, site error, gyro error and accelerometer error, φ _e, φ _nand φ _ube expressed as carrier angle of pitch error, roll angle error and azimuth angle error, δ V _e, δ V _nwith δ V _ube expressed as the east orientation velocity error of carrier, north orientation velocity error and sky to velocity error, δ L, δ λ and δ h are expressed as the latitude error of carrier, longitude error and height error, ε _e, ε _nand ε _ube expressed as gyroscope in east orientation, north orientation and sky to drift, ▽ _e, ▽ _nand ▽ _ube expressed as accelerometer in east orientation, north orientation and sky to output error;

Step 3: the error rate of computing gyroscope and accelerometer;

Step 4: according to the parameter in step 1-3, calculates attitude of carrier angle error rate of change, velocity error rate of change and site error rate of change;

Step 5: introduce GPS pseudo range measurement information, adopts feedback compensation mode to correct SINS output information, obtains current position accurately and attitude;

Code phase error measured by GPS track loop calculates the pseudorange between satellite to carrier, pseudorange upgrades the quantity of state of Navigation Filter as measurement information, thus dope the error of current SINS institute state quantity measurement amount, then adopt feedback compensation mode to correct SINS output information, obtain current position accurately and attitude;

Step 6: calculate the coordinate of carrier in ECEF coordinate system;

The carrier positions calculated by step 5 is arranged in earth coordinates, i.e. (L, λ, h), is translated in ECEF coordinate system, its position coordinates be expressed as

X_{p}^{e} = {[(R + h) \cos L c o s λ, (R + h) \cos L s i n λ, (R + h) \sin L]}^{T}

Wherein, R is earth radius;

Step 7: calculate the transition matrix that carrier coordinate system and b are tied to the transition matrix of navigational coordinate system and n system, navigation coordinate is tied to ECEF coordinate system and e system;

The attitude of carrier angle calculated by step 5, calculates the transition matrix of carrier coordinate system to navigational coordinate system for

C_{b}^{n} = [\begin{matrix} \cos γ \cos ψ - \sin ψ \sin θ \sin γ & - \sin ψ \cos θ & \sin γ \cos ψ + \cos γ \sin θ \sin ψ \\ \cos γ \sin ψ + \sin γ \sin θ \cos ψ & \cos ψ \cos θ & \sin γ \sin ψ - \cos γ \sin θ \cos ψ \\ - \sin γ \cos θ & \sin θ & \cos θ \cos γ \end{matrix}]

The carrier positions calculated by step 5, calculates the transition matrix that navigation coordinate is tied to ECEF coordinate system for

C_{n}^{e} = [\begin{matrix} - \sin λ & - \sin L \cos λ & \cos L \cos λ \\ \cos λ & - \sin L \cos λ & \cos L \sin λ \\ 0 & \cos L & \sin L \end{matrix}];

Step 8: calculate the steering vector between carrier to satellite;

By the transition matrix in step 7 transition matrix try to achieve the transformed matrix of carrier coordinate system to ECEF coordinate system

C_{b}^{e} = C_{n}^{e} C_{b}^{n}

Satellite ephemeris is utilized to calculate the position of satellite in ECEF coordinate system and in conjunction with formula

X_{p}^{e} = {[(R + h) \cos L c o s λ, (R + h) \cos L \sin λ, (R + h) \sin L]}^{T}

With

C_{b}^{e} = C_{n}^{e} C_{b}^{n},

Try to achieve the steering vector between carrier to satellite

{\overset{&RightArrow;}{r}}^{b} = {(C_{b}^{e})}^{T} (X_{s}^{e} - X_{p}^{e})

Step 9: in carrier coordinate system, calculates position angle and the angle of pitch that satellite arrives antenna;

If will coordinate definition in carrier coordinate system is then satellite arrives the azimuth angle alpha of antenna and angle of pitch β can be expressed as respectively

α = a r c t a n (\hat{x}, \hat{y})

β = a r c t a n (\hat{z}, \sqrt{{\hat{x}}^{2} + {\hat{y}}^{2}})

Step 10: design circular array antenna structure;

6 array element is uniformly distributed in round battle array circumferentially, radius of circle is made to be r, interval then circumferentially between adjacent array element is also r, choosing of array element interval, and it is the same that time-domain sampling interval meets Nyquist's theorem, Space domain sampling interval d should be less than 1/2 of satellite carrier wavelength X, by satellite frequency f=1575.42 × 10 ⁶mHz, so array element distance is

d \leq \frac{λ}{2} = \frac{1}{2} \cdot \frac{c}{f_{l 1}} = \frac{1}{2} \cdot \frac{3 \times 10^{8}}{1575.42 \times 10^{6}} = \frac{1}{2} \cdot 0.19 = 0.095 m = 9.5 c m

Wherein, c is the light velocity;

Step 11: calculate the time delay between each array element of satellite-signal arrival antenna;

According in step 9

α = a r c t a n (\hat{x}, \hat{y}),

β = a r c t a n (\hat{z}, \sqrt{{\hat{x}}^{2} + {\hat{y}}^{2}}),

By α and β representation unit vector

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

Therefore, satellite-signal arrives i-th bay and to the mistiming τ between first reference array element _ican be expressed as

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

Wherein, M representative antennas array element number;

Step 12: set up array antenna received signals model

If array antenna received has arrived P satellite-signal, Q undesired signal, then the signal model that antenna receives has been expressed as

U (t) = Σ_{k = 1}^{P} a (α_{k}, β_{k}, T_{s}) s_{k} (t) + Σ_{k = P + 1}^{P + Q} a (α_{k}, β_{k}, T_{s}) j_{k - P} (t) + n (t)

Wherein, s (t) and j (t) represents the satellite-signal and undesired signal that receive, a (α respectively _k, β _k, T _s) be the steering vector of a kth echo signal, T _stime domain delay line interval, α _kand β _kbe expressed as position angle and the angle of pitch that kth echo signal arrives array antenna, n (t) represents white Gaussian noise, its power spectrum density is expressed as N ₀/ 2;

A (α _k, β _k, T _s) represent space-time two-dimensional target vector, i.e. time vector a _s(T _s) and direction in space vector a _s(α _k, β _k) Crow Neck amass, and to be expressed as:

a (α_{k}, β_{k}, T_{s}) = a_{s} (α_{k}, β_{k}) &CircleTimes; a_{s} (T_{s})

The time delay unit number of each radio-frequency channel is made to be N, then

a_{s} (α_{k}, β_{k}) = [\begin{matrix} 1 & e^{- j 2 {πfτ}_{1}} & ... & e^{- j 2 {πfτ}_{M - 1}} \end{matrix}]

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

F represents satellite frequency, T _srepresent that the time interval, its value of delay cell should be less than signal bandwidth;

Step 13: calculate Antenna array weight vector;

Adopt multiple constraint minimum variance space-time adaptive Processing Algorithm to retrain satellites in view signal, then make array output power minimum, being expressed as follows of the best initial weights of array:

\{\begin{matrix} ω = \arg {minω}^{H} R_{U} ω \\ s u b j e c t t o : ω^{H} A = F^{T} \end{matrix}

Wherein, ω represents the weighted vector of array, the conjugate transpose of H representing matrix, R _ube the array covariance matrix of input signal, be expressed as

R _U＝E{UU ^H}

A represents the constraint matrix of satellite-signal, and F represents the constraint vector corresponding with A, A and F is expressed as

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

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

Adopt method of Lagrange multipliers, array weight vector ω is expressed as

ω = R_{U}^{- 1} A {(A^{H} R_{U}^{- 1} A)}^{- 1} f

Acquisition array weight after, obtain array output express be for

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