[go: up one dir, main page]

CN115096303A - GNSS multi-antenna and INS tightly-combined positioning and attitude determination method and equipment - Google Patents

GNSS multi-antenna and INS tightly-combined positioning and attitude determination method and equipment Download PDF

Info

Publication number
CN115096303A
CN115096303A CN202211024566.9A CN202211024566A CN115096303A CN 115096303 A CN115096303 A CN 115096303A CN 202211024566 A CN202211024566 A CN 202211024566A CN 115096303 A CN115096303 A CN 115096303A
Authority
CN
China
Prior art keywords
vector
matrix
antenna
double
difference
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202211024566.9A
Other languages
Chinese (zh)
Other versions
CN115096303B (en
Inventor
余文坤
李鑫
戴吾蛟
潘林
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202211024566.9A priority Critical patent/CN115096303B/en
Publication of CN115096303A publication Critical patent/CN115096303A/en
Application granted granted Critical
Publication of CN115096303B publication Critical patent/CN115096303B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/10—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
    • G01C21/165—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation combined with non-inertial navigation instruments
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/20—Instruments for performing navigational calculations
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • G01S19/48—Determining position by combining or switching between position solutions derived from the satellite radio beacon positioning system and position solutions derived from a further system
    • G01S19/49—Determining position by combining or switching between position solutions derived from the satellite radio beacon positioning system and position solutions derived from a further system whereby the further system is an inertial position system, e.g. loosely-coupled

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

The invention discloses a GNSS multi-antenna and INS tightly combined positioning and attitude determining method and equipment, wherein the method comprises the following steps: constructing a double-difference carrier and double-difference pseudorange rate measurement equation represented by a position vector and a velocity vector at an inertial navigation center through a lever arm vector and a rotation matrix; considering the correlation when the multi-antenna GNSS observation values form double-difference observation values, and constructing a covariance matrix of the observation values of the GNSS multi-antenna and INS tight combination system by using an error propagation law; according to the double-difference carrier wave, the double-difference pseudo range rate equation and the covariance matrix, a measurement equation based on variable correction of a tight combination system is constructed; and taking errors such as position, speed, attitude and the like as state vectors of the system, establishing a system state equation by adopting a first-order Gaussian Markov process, and carrying out Kalman filtering solution to obtain the optimal estimation of the position, the speed and the attitude of the inertial navigation center at each moment. The invention has high stability and precision of positioning and attitude determination.

Description

一种GNSS多天线与INS紧组合定位定姿方法和设备A kind of GNSS multi-antenna and INS tight combination positioning and attitude determination method and equipment

技术领域technical field

本发明属于导航定位技术领域，具体涉及一种GNSS多天线与INS紧组合定位定姿方法和设备。The invention belongs to the technical field of navigation and positioning, and in particular relates to a method and equipment for positioning and attitude determination by tight combination of GNSS multi-antenna and INS.

背景技术Background technique

准确可靠的位置和姿态等信息在载体导航、制导和控制中发挥着关键作用，GNSS（Global Navigation Satellite System，全球卫星导航系统）和INS（InertialNavigation System，惯性导航系统）是获取位姿信息的两大主流手段。GNSS系统能够在全球范围为用户提供高精度的定位导航服务，但是其在动态环境中容易受到干扰造成卫星信号的失锁，导致定位结果不可用。INS系统具有较强的自主性和抗干扰能力，但是其每一时刻的结果都由前一时刻递推而来，长时间工作时误差会随时间累积，导致导航结果不可靠。目前，常将这两个系统组合起来进行定位定姿，以提高导航系统的稳定性和可靠性。Accurate and reliable information such as position and attitude plays a key role in carrier navigation, guidance and control. GNSS (Global Navigation Satellite System, global satellite navigation system) and INS (Inertial Navigation System, inertial navigation system) are two ways to obtain position and attitude information. Mainstream means. The GNSS system can provide users with high-precision positioning and navigation services on a global scale, but it is susceptible to interference in a dynamic environment, resulting in loss of satellite signal lock, resulting in unusable positioning results. The INS system has strong autonomy and anti-interference ability, but the results of each moment are recursively derived from the previous moment, and errors will accumulate over time when working for a long time, resulting in unreliable navigation results. At present, these two systems are often combined for positioning and attitude determination to improve the stability and reliability of the navigation system.

现有的GNSS与INS定位定姿方法都存在其技术问题：The existing GNSS and INS positioning and attitude determination methods have their technical problems:

公开号为CN103245963A的专利方案提供一种双天线GNSS/INS深组合导航方法及装置，其使用双天线GNSS差分解算的姿态对INS系统姿态进行约束。但是由双天线GNSS差分解算的基线组合对姿态进行约束属于松组合，未充分利用数据，抗干扰的能力差，而且首先需要固定基线向量的模糊度，如果模糊度不能固定或者错误固定，会将在误差带入计算的姿态中，影响最终的导航结果。The patent solution with the publication number CN103245963A provides a dual-antenna GNSS/INS deep integrated navigation method and device, which uses the dual-antenna GNSS differentially calculated attitude to constrain the attitude of the INS system. However, the attitude constraint of the baseline combination calculated by the dual-antenna GNSS differential decomposition is a loose combination, the data is not fully utilized, and the anti-interference ability is poor, and the ambiguity of the baseline vector needs to be fixed first. If the ambiguity cannot be fixed or wrongly fixed, it will Bring the error into the calculated attitude and affect the final navigation result.

公开号为CN114252077A的专利方案提供一种基于联邦滤波器的双天线GPS/SINS的组合导航方法及系统，其采用联邦滤波器将两个GPS系统的观测分别与INS组合进行滤波，以提高系统的健壮性，但是该方案对GPS与INS的组合也是采用的松组合模式，抗干扰的能力和解算的精度都不及紧组合模式，同时，一旦某个子系统出现故障，其观测信息均会被删除，没有实现对多个天线的冗余观测的有效利用。The patent scheme with publication number CN114252077A provides a combined navigation method and system of dual-antenna GPS/SINS based on a federated filter, which uses a federated filter to filter the observations of the two GPS systems and the INS respectively, so as to improve the system's performance. Robustness, but this scheme also adopts the loose combination mode for the combination of GPS and INS. The anti-interference ability and the accuracy of the solution are not as good as the tight combination mode. At the same time, once a subsystem fails, its observation information will be deleted. Effective utilization of redundant observations of multiple antennas is not achieved.

柴艳菊等人提供了一种多天线GNSS/INS组合导航算法及结果分析（柴艳菊,胡付帅,钟世明. 多天线GNSS/INS组合导航算法及结果分析[C],2021:146-150.），采用选权自适应滤波实现多天线GNSS观测信息对INS误差快速校正，但是其对于GNSS和INS观测信息的组合仍然采用的是松组合的模式，在GNSS系统受到干扰时，该方法的导航精度会下降甚至不可用。Chai Yanju et al. provided a multi-antenna GNSS/INS integrated navigation algorithm and result analysis (Chai Yanju, Hu Fushuai, Zhong Shiming. Multi-antenna GNSS/INS integrated navigation algorithm and result analysis [C], 2021:146-150.), using the selected The weighted adaptive filtering can quickly correct the INS error by the multi-antenna GNSS observation information, but it still adopts a loose combination mode for the combination of GNSS and INS observation information. When the GNSS system is interfered, the navigation accuracy of this method will decrease or even unavailable.

在GNSS与INS的组合定位定姿问题中，如何设计合理的函数模型、充分有效地利用多个天线的冗余观测信息是一个关键的问题。传统的解决方案是首先求解主天线和从天线之间的基线向量，利用基线向量与载体姿态之间的关系建立多GNSS天线与INS组合的量测方程，但这种方式通常属于松组合，并且严重依赖于整周模糊度的正确固定。In the problem of combined positioning and attitude determination of GNSS and INS, how to design a reasonable function model and fully and effectively utilize the redundant observation information of multiple antennas is a key issue. The traditional solution is to first solve the baseline vector between the master antenna and the slave antenna, and use the relationship between the baseline vector and the carrier attitude to establish the measurement equation of the combination of multiple GNSS antennas and INS, but this method is usually a loose combination, and Heavy reliance on correct fixation of the ambiguity of the whole week.

发明内容SUMMARY OF THE INVENTION

针对上述GNSS与INS的组合定位定姿的问题，本发明提供一种GNSS多天线与INS紧组合定位定姿方法和设备，将多个天线的相位中心统一归算到惯导中心位置，构建了一种待估坐标参数少、模型几何强度大的观测模型，从而实现多GNSS天线与INS的紧组合。Aiming at the above-mentioned problem of combined positioning and attitude determination of GNSS and INS, the present invention provides a method and device for combined positioning and attitude determination of GNSS multi-antenna and INS, which uniformly reduces the phase centers of multiple antennas to the position of the inertial navigation center, and constructs a An observation model with few coordinate parameters to be estimated and high geometric strength of the model, so as to realize the tight combination of multiple GNSS antennas and INS.

为实现上述技术目的，本发明采用如下技术方案：For realizing the above-mentioned technical purpose, the present invention adopts following technical scheme:

一种GNSS多天线与INS紧组合定位定姿方法，包括：A GNSS multi-antenna and INS tight combination positioning and attitude determination method, comprising:

步骤1，通过杠杆臂向量和旋转矩阵，将多个GNSS天线的坐标参数归化到惯导中心坐标参数，构建以惯导中心处的位置向量和速度向量表示的双差载波和双差伪距率量测方程；Step 1: Through the lever arm vector and rotation matrix, the coordinate parameters of multiple GNSS antennas are normalized to the coordinate parameters of the inertial navigation center, and the double-difference carrier and double-difference pseudorange represented by the position vector and velocity vector at the inertial navigation center are constructed. rate measurement equation;

步骤2，顾及多天线GNSS观测值组成双差观测值时的相关性，利用误差传播定律构建GNSS多天线与INS紧组合系统观测值的协方差矩阵；Step 2: Considering the correlation when the multi-antenna GNSS observations form double-difference observations, the error propagation law is used to construct the covariance matrix of the observations of the GNSS multi-antenna and INS compact combination system;

步骤3，根据步骤1得到的双差载波和双差伪距率量测方程，以及步骤2得到的多天线观测协方差矩阵，构建得到GNSS多天线与INS紧组合系统的基于变量改正数的量测方程；Step 3, according to the double-difference carrier and double-difference pseudorange rate measurement equations obtained in step 1, and the multi-antenna observation covariance matrix obtained in step 2, construct and obtain the variable correction number based on the GNSS multi-antenna and INS compact combination system. measuring equation;

步骤4，以惯导中心处的位置、速度、姿态、传感器零偏及比例因子的误差作为系统的状态向量，采用一阶高斯马尔科夫过程建立系统状态方程；Step 4, taking the position, velocity, attitude, sensor bias and scale factor error at the center of inertial navigation as the state vector of the system, and using the first-order Gaussian Markov process to establish the system state equation;

步骤5，根据上述量测方程和系统状态方程进行卡尔曼滤波解算，得到惯导中心在每个时刻的位置、速度、姿态的最优估计。In step 5, Kalman filtering is performed according to the above measurement equation and the system state equation to obtain the optimal estimation of the position, velocity and attitude of the inertial navigation center at each moment.

进一步地，将GNSS天线的坐标参数归化到惯导中心坐标参数的方法为：Further, the method for normalizing the coordinate parameters of the GNSS antenna to the coordinate parameters of the inertial navigation center is:

Figure DEST_PATH_IMAGE001
Figure DEST_PATH_IMAGE001

其中，

Figure DEST_PATH_IMAGE002
、
Figure DEST_PATH_IMAGE003
分别表示天线
Figure DEST_PATH_IMAGE004
的位置向量和速度向量，
Figure DEST_PATH_IMAGE005
、v分别表示惯导中心 的位置向量和速度向量；
Figure DEST_PATH_IMAGE006
是方向余弦矩阵，即旋转矩阵；
Figure DEST_PATH_IMAGE007
是天线
Figure 472470DEST_PATH_IMAGE004
的杠杆臂向量；
Figure DEST_PATH_IMAGE008
是IMU输出的三维角速度
Figure DEST_PATH_IMAGE009
的反对称矩阵，即
Figure DEST_PATH_IMAGE010
，
Figure DEST_PATH_IMAGE011
为反对称阵算子；
Figure DEST_PATH_IMAGE012
是地 球自转角速度的反对称矩阵。 in,
Figure DEST_PATH_IMAGE002
,
Figure DEST_PATH_IMAGE003
Respectively represent the antenna
Figure DEST_PATH_IMAGE004
The position vector and velocity vector of ,
Figure DEST_PATH_IMAGE005
, v represent the position vector and velocity vector of the inertial navigation center, respectively;
Figure DEST_PATH_IMAGE006
is the direction cosine matrix, that is, the rotation matrix;
Figure DEST_PATH_IMAGE007
is the antenna
Figure 472470DEST_PATH_IMAGE004
the lever arm vector of ;
Figure DEST_PATH_IMAGE008
is the three-dimensional angular velocity output by the IMU
Figure DEST_PATH_IMAGE009
The antisymmetric matrix of , that is
Figure DEST_PATH_IMAGE010
,
Figure DEST_PATH_IMAGE011
is an antisymmetric matrix operator;
Figure DEST_PATH_IMAGE012
is the antisymmetric matrix of the angular velocity of the Earth's rotation.

进一步地，构建以惯导中心处的位置向量和速度向量表示的双差载波和双差伪距率方程，表示为：Further, construct the double-difference carrier and double-difference pseudorange rate equations represented by the position vector and velocity vector at the center of the inertial navigation, which are expressed as:

Figure DEST_PATH_IMAGE013
Figure DEST_PATH_IMAGE013

其中，

Figure DEST_PATH_IMAGE014
为天线
Figure 717506DEST_PATH_IMAGE004
处星间单差视线向量矩阵，它是一个
Figure DEST_PATH_IMAGE015
维的矩阵，
Figure DEST_PATH_IMAGE016
为观测卫星总数减去所采用的卫星系统数；
Figure DEST_PATH_IMAGE017
为中间量，
Figure DEST_PATH_IMAGE018
为由相应载波波长构成的非零对角阵，
Figure DEST_PATH_IMAGE019
表示天线
Figure DEST_PATH_IMAGE020
之间双差观测噪声；
Figure DEST_PATH_IMAGE021
表示天线
Figure 575872DEST_PATH_IMAGE020
之间双差载波观测值，
Figure DEST_PATH_IMAGE022
表示天线
Figure 502239DEST_PATH_IMAGE020
之间的双差伪距率观测值。 in,
Figure DEST_PATH_IMAGE014
for the antenna
Figure 717506DEST_PATH_IMAGE004
single-difference line-of-sight vector matrix between the stars, which is a
Figure DEST_PATH_IMAGE015
dimensional matrix,
Figure DEST_PATH_IMAGE016
Subtract the number of satellite systems used from the total number of observation satellites;
Figure DEST_PATH_IMAGE017
is an intermediate quantity,
Figure DEST_PATH_IMAGE018
is a non-zero diagonal matrix composed of the corresponding carrier wavelengths,
Figure DEST_PATH_IMAGE019
Indicates the antenna
Figure DEST_PATH_IMAGE020
between double-difference observation noise;
Figure DEST_PATH_IMAGE021
Indicates the antenna
Figure 575872DEST_PATH_IMAGE020
between double-difference carrier observations,
Figure DEST_PATH_IMAGE022
Indicates the antenna
Figure 502239DEST_PATH_IMAGE020
The double-difference pseudorange rate between observations.

进一步地，步骤3得到的量测方程表示为：Further, the measurement equation obtained in step 3 is expressed as:

Figure DEST_PATH_IMAGE023
Figure DEST_PATH_IMAGE023

其中，

Figure DEST_PATH_IMAGE024
表示变量的改正数；
Figure DEST_PATH_IMAGE025
表示GNSS天线组合系统的观测向量，包括GNSS 系统观测得到的载波和伪距率，
Figure DEST_PATH_IMAGE026
表示观测向量
Figure 365547DEST_PATH_IMAGE025
的改正数，由GNSS的载波和伪距率与 惯导推算的载波和伪距率做差得到；
Figure DEST_PATH_IMAGE027
表示状态向量
Figure DEST_PATH_IMAGE028
的改正；
Figure 100002_DEST_PATH_IMAGE029
表示观测向量改正数
Figure 984747DEST_PATH_IMAGE026
与状态向量改正数
Figure 517359DEST_PATH_IMAGE027
之间的联系矩阵；
Figure DEST_PATH_IMAGE030
表示白噪声；
Figure DEST_PATH_IMAGE031
表示观测值的协方差矩阵； in,
Figure DEST_PATH_IMAGE024
Represents the correction number of the variable;
Figure DEST_PATH_IMAGE025
represents the observation vector of the GNSS antenna combination system, including the carrier and pseudorange rates observed by the GNSS system,
Figure DEST_PATH_IMAGE026
represents the observation vector
Figure 365547DEST_PATH_IMAGE025
The correction number is obtained by the difference between the carrier and pseudorange rates of GNSS and the carrier and pseudorange rates estimated by inertial navigation;
Figure DEST_PATH_IMAGE027
Represents a state vector
Figure DEST_PATH_IMAGE028
correction;
Figure 100002_DEST_PATH_IMAGE029
Indicates the observation vector correction number
Figure 984747DEST_PATH_IMAGE026
with the state vector correction number
Figure 517359DEST_PATH_IMAGE027
The relationship matrix between;
Figure DEST_PATH_IMAGE030
represents white noise;
Figure DEST_PATH_IMAGE031
represents the covariance matrix of observations;

状态向量

Figure 247418DEST_PATH_IMAGE028
为： state vector
Figure 247418DEST_PATH_IMAGE028
for:

Figure DEST_PATH_IMAGE032
Figure DEST_PATH_IMAGE032

其中，

Figure 181876DEST_PATH_IMAGE005
为载体三维位置向量，定为惯导参考中心；
Figure DEST_PATH_IMAGE033
代表速度向量；
Figure DEST_PATH_IMAGE034
代表姿态 向量；
Figure DEST_PATH_IMAGE035
代表加速度计零偏；
Figure DEST_PATH_IMAGE036
代表陀螺零偏；
Figure DEST_PATH_IMAGE037
代表加速度计比例因子；
Figure DEST_PATH_IMAGE038
代表陀 螺比例因子；
Figure DEST_PATH_IMAGE039
代表双差模糊度向量。 in,
Figure 181876DEST_PATH_IMAGE005
is the three-dimensional position vector of the carrier, and is set as the inertial navigation reference center;
Figure DEST_PATH_IMAGE033
represents the velocity vector;
Figure DEST_PATH_IMAGE034
represents the pose vector;
Figure DEST_PATH_IMAGE035
represents the zero offset of the accelerometer;
Figure DEST_PATH_IMAGE036
Represents the zero bias of the gyro;
Figure DEST_PATH_IMAGE037
represents the accelerometer scale factor;
Figure DEST_PATH_IMAGE038
represents the gyro scale factor;
Figure DEST_PATH_IMAGE039
represents the double-difference ambiguity vector.

进一步地，观测向量改正数

Figure 253868DEST_PATH_IMAGE026
和联系矩阵
Figure 273777DEST_PATH_IMAGE029
分别为： Further, the observation vector correction number
Figure 253868DEST_PATH_IMAGE026
and contact matrix
Figure 273777DEST_PATH_IMAGE029
They are:

Figure DEST_PATH_IMAGE040
Figure DEST_PATH_IMAGE040

其中，

Figure DEST_PATH_IMAGE041
表示主天线，
Figure DEST_PATH_IMAGE042
表示
Figure DEST_PATH_IMAGE043
个辅天线，
Figure DEST_PATH_IMAGE044
表示
Figure DEST_PATH_IMAGE045
个基准站；
Figure 100002_DEST_PATH_IMAGE046
表示下标
Figure DEST_PATH_IMAGE047
之间的双差载波观测值，
Figure DEST_PATH_IMAGE048
表示下标
Figure 853532DEST_PATH_IMAGE047
之间的双差伪距率观测 值；
Figure DEST_PATH_IMAGE049
表示下标
Figure DEST_PATH_IMAGE050
相关的双差向量改正数与状态向量改正数之间的联系矩阵，且有： in,
Figure DEST_PATH_IMAGE041
represents the main antenna,
Figure DEST_PATH_IMAGE042
express
Figure DEST_PATH_IMAGE043
an auxiliary antenna,
Figure DEST_PATH_IMAGE044
express
Figure DEST_PATH_IMAGE045
base stations;
Figure 100002_DEST_PATH_IMAGE046
Indicates subscript
Figure DEST_PATH_IMAGE047
The double-difference carrier observations between,
Figure DEST_PATH_IMAGE048
Indicates subscript
Figure 853532DEST_PATH_IMAGE047
The double-difference pseudorange rate between observations;
Figure DEST_PATH_IMAGE049
Indicates subscript
Figure DEST_PATH_IMAGE050
The correlation matrix between the related double-difference vector corrections and the state vector corrections, and has:

Figure DEST_PATH_IMAGE051
Figure DEST_PATH_IMAGE051

其中，

Figure DEST_PATH_IMAGE052
表示
Figure 100002_DEST_PATH_IMAGE053
维的零矩阵，
Figure DEST_PATH_IMAGE054
为由相应载波波长构成的非零对 角阵，各中间变量
Figure DEST_PATH_IMAGE055
的具体表达如下： in,
Figure DEST_PATH_IMAGE052
express
Figure 100002_DEST_PATH_IMAGE053
dimensional zero matrix,
Figure DEST_PATH_IMAGE054
is a non-zero diagonal matrix formed by the corresponding carrier wavelength, each intermediate variable
Figure DEST_PATH_IMAGE055
The specific expression is as follows:

Figure 100002_DEST_PATH_IMAGE056
。
Figure 100002_DEST_PATH_IMAGE056
.

进一步地，利用误差传播定律构建GNSS多天线与INS紧组合系统观测值的协方差矩阵R的方法为：Further, the method of constructing the covariance matrix R of the observations of the GNSS multi-antenna and INS compact combination system using the law of error propagation is:

首先，将天线

Figure 252283DEST_PATH_IMAGE020
和卫星
Figure DEST_PATH_IMAGE057
,
Figure 213286DEST_PATH_IMAGE045
之间的双差观测值表示为
Figure DEST_PATH_IMAGE058
，将任意天线
Figure 720491DEST_PATH_IMAGE004
对任 意卫星
Figure 995614DEST_PATH_IMAGE057
的原始观测值表示为
Figure DEST_PATH_IMAGE059
，则由误差传播定律可得原始非差观测值
Figure 416581DEST_PATH_IMAGE059
与双差观测 值
Figure 17327DEST_PATH_IMAGE058
之间的关系为： First, place the antenna
Figure 252283DEST_PATH_IMAGE020
and satellite
Figure DEST_PATH_IMAGE057
,
Figure 213286DEST_PATH_IMAGE045
The double-differenced observations between are expressed as
Figure DEST_PATH_IMAGE058
, placing any antenna
Figure 720491DEST_PATH_IMAGE004
to any satellite
Figure 995614DEST_PATH_IMAGE057
The raw observations of are expressed as
Figure DEST_PATH_IMAGE059
, then the original non-difference observations can be obtained from the law of error propagation
Figure 416581DEST_PATH_IMAGE059
with double difference observations
Figure 17327DEST_PATH_IMAGE058
The relationship between is:

Figure DEST_PATH_IMAGE060
Figure DEST_PATH_IMAGE060

然后，根据上述关系式计算其中的转化矩阵

Figure DEST_PATH_IMAGE061
为： Then, according to the above relationship, calculate the transformation matrix in it
Figure DEST_PATH_IMAGE061
for:

Figure DEST_PATH_IMAGE062
Figure DEST_PATH_IMAGE062

其中，

Figure DEST_PATH_IMAGE063
表示
Figure DEST_PATH_IMAGE064
维的零矩阵，
Figure DEST_PATH_IMAGE065
的具体表达为： in,
Figure DEST_PATH_IMAGE063
express
Figure DEST_PATH_IMAGE064
dimensional zero matrix,
Figure DEST_PATH_IMAGE065
The specific expression is:

Figure DEST_PATH_IMAGE066
Figure DEST_PATH_IMAGE066

最后，根据转化矩阵

Figure 74145DEST_PATH_IMAGE061
以及原始非差观测值的协方差矩阵
Figure DEST_PATH_IMAGE067
，确定组合系统 观测值的协方差矩阵为： Finally, according to the transformation matrix
Figure 74145DEST_PATH_IMAGE061
and the covariance matrix of the original non-differenced observations
Figure DEST_PATH_IMAGE067
, determine the covariance matrix of the combined system observations as:

Figure DEST_PATH_IMAGE068
。
Figure DEST_PATH_IMAGE068
.

进一步地，观测向量改正数

Figure 434850DEST_PATH_IMAGE026
的计算方法为：根据惯导中心更新的位置、速度和 姿态信息以及卫星星历，推算双差载波观测值和双差伪距率观测值；将其与GNSS系统观测 得到的载波和伪距率做差，作为量测方程中观测向量的改正数。 Further, the observation vector correction number
Figure 434850DEST_PATH_IMAGE026
The calculation method is: according to the updated position, velocity and attitude information and satellite ephemeris of the inertial navigation center, calculate the double-difference carrier observation value and the double-difference pseudorange rate observation value; compare it with the carrier and pseudorange rate observed by the GNSS system. The difference is used as the correction number for the observation vector in the measurement equation.

进一步地，根据对多天线GNSS系统观测值构建双差观测方程，通过解算双差观测方程得到组合系统的位置、速度和姿态信息，将其用于惯导中心的初始化。Furthermore, a double-difference observation equation is constructed according to the observations of the multi-antenna GNSS system, and the position, velocity and attitude information of the combined system are obtained by solving the double-difference observation equation, which is used for the initialization of the inertial navigation center.

进一步地，采用一阶高斯马尔科夫过程建立系统的状态方程为：Further, the state equation of the system established by the first-order Gaussian Markov process is:

Figure DEST_PATH_IMAGE069
Figure DEST_PATH_IMAGE069

其中

Figure DEST_PATH_IMAGE070
是状态向量
Figure 463986DEST_PATH_IMAGE028
在时刻
Figure DEST_PATH_IMAGE071
的估值，上标“—”表示某变量还未被最新 的观测所更新，一旦更新，相应的上标将变为“+”；
Figure DEST_PATH_IMAGE072
为状态转移矩阵，
Figure DEST_PATH_IMAGE073
为过程噪声向 量，假定服从零均值正态分布，
Figure DEST_PATH_IMAGE074
为过程噪声的转移矩阵，
Figure DEST_PATH_IMAGE075
为过程噪声的协方差矩阵；
Figure 297950DEST_PATH_IMAGE072
的具体表达为： in
Figure DEST_PATH_IMAGE070
is the state vector
Figure 463986DEST_PATH_IMAGE028
at the moment
Figure DEST_PATH_IMAGE071
The superscript "—" indicates that a variable has not been updated by the latest observation. Once updated, the corresponding superscript will become "+";
Figure DEST_PATH_IMAGE072
is the state transition matrix,
Figure DEST_PATH_IMAGE073
is the process noise vector, assuming a zero-mean normal distribution,
Figure DEST_PATH_IMAGE074
is the transition matrix of process noise,
Figure DEST_PATH_IMAGE075
is the covariance matrix of the process noise;
Figure 297950DEST_PATH_IMAGE072
The specific expression is:

Figure DEST_PATH_IMAGE076
Figure DEST_PATH_IMAGE076

Figure DEST_PATH_IMAGE077
Figure DEST_PATH_IMAGE077

Figure DEST_PATH_IMAGE078
Figure DEST_PATH_IMAGE078

其中，状态转移矩阵

Figure 29014DEST_PATH_IMAGE072
中各元素下标中的3表示3×3的矩阵，n表示n×n的矩阵，3 ×n表示3行n列的矩阵，n×3表示n行3列的矩阵；
Figure 177099DEST_PATH_IMAGE005
为惯导中心的三维位置向量；
Figure DEST_PATH_IMAGE079
为地 表与地心之间的距离；
Figure DEST_PATH_IMAGE080
为重力向量；
Figure DEST_PATH_IMAGE081
是加速度计输出的三维比力向量；
Figure 60741DEST_PATH_IMAGE009
为IMU输出的 三维角速度；
Figure DEST_PATH_IMAGE082
为以
Figure DEST_PATH_IMAGE083
为对角线元素的对角矩阵；
Figure DEST_PATH_IMAGE084
均为中间变量；
Figure DEST_PATH_IMAGE085
为加速度计 零偏的相关时间，
Figure DEST_PATH_IMAGE086
为陀螺零偏的相关时间，
Figure DEST_PATH_IMAGE087
为加速度计比例因子的相关时 间，
Figure DEST_PATH_IMAGE088
为陀螺比例因子的相关时间；
Figure DEST_PATH_IMAGE089
为加速度计和陀螺测量时间间隔；
Figure DEST_PATH_IMAGE090
是地球自 转角速度的反对称矩阵； Among them, the state transition matrix
Figure 29014DEST_PATH_IMAGE072
The 3 in the subscript of each element in the subscript represents a 3×3 matrix, n represents an n×n matrix, 3×n represents a matrix with 3 rows and n columns, and n×3 represents a matrix with n rows and 3 columns;
Figure 177099DEST_PATH_IMAGE005
is the three-dimensional position vector of the inertial navigation center;
Figure DEST_PATH_IMAGE079
is the distance between the surface and the center of the earth;
Figure DEST_PATH_IMAGE080
is the gravity vector;
Figure DEST_PATH_IMAGE081
is the three-dimensional specific force vector output by the accelerometer;
Figure 60741DEST_PATH_IMAGE009
is the three-dimensional angular velocity output by the IMU;
Figure DEST_PATH_IMAGE082
for
Figure DEST_PATH_IMAGE083
is a diagonal matrix of diagonal elements;
Figure DEST_PATH_IMAGE084
are intermediate variables;
Figure DEST_PATH_IMAGE085
is the relative time of the zero offset of the accelerometer,
Figure DEST_PATH_IMAGE086
is the relative time of the zero bias of the gyro,
Figure DEST_PATH_IMAGE087
is the correlation time of the accelerometer scale factor,
Figure DEST_PATH_IMAGE088
is the relative time of the gyro scale factor;
Figure DEST_PATH_IMAGE089
Measure time intervals for accelerometers and gyroscopes;
Figure DEST_PATH_IMAGE090
is the antisymmetric matrix of the angular velocity of the Earth's rotation;

根据量测方程和系统状态方程进行卡尔曼滤波解算的具体递推计算公式为：The specific recursive calculation formula for Kalman filter solution based on the measurement equation and the system state equation is:

Figure DEST_PATH_IMAGE091
Figure DEST_PATH_IMAGE091

其中，

Figure DEST_PATH_IMAGE092
是状态向量的方差协方差阵；
Figure DEST_PATH_IMAGE093
为Kalman增益矩阵；
Figure 144235DEST_PATH_IMAGE031
为观测向量的方 差协方差阵；
Figure DEST_PATH_IMAGE094
为单位阵。 in,
Figure DEST_PATH_IMAGE092
is the variance-covariance matrix of the state vector;
Figure DEST_PATH_IMAGE093
is the Kalman gain matrix;
Figure 144235DEST_PATH_IMAGE031
is the variance covariance matrix of the observation vector;
Figure DEST_PATH_IMAGE094
is a unit matrix.

一种电子设备，包括存储器及处理器，所述存储器中存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器实现上述任一项技术方案中所述的GNSS多天线与INS紧组合定位定姿方法。An electronic device, comprising a memory and a processor, wherein a computer program is stored in the memory, and when the computer program is executed by the processor, the processor is made to implement the GNSS multiplexer described in any of the above technical solutions. The antenna and INS are closely combined positioning and attitude determination method.

有益效果beneficial effect

传统GNSS/INS定位定姿方法首先求解主天线和从天线之间的基线向量，利用基线向量求解载体姿态，据此建立多GNSS天线与INS组合的量测方程，这种方式依赖于先正确固定整周模糊度后的基线值，属于松组合。本发明中的方法避免了这一问题，以观测值层面的紧组合有效利用了多天线的几何构型信息和冗余观测信息，可以极大程度提高集成系统定位定姿结果的精度与可靠性。The traditional GNSS/INS positioning and attitude determination method firstly solves the baseline vector between the main antenna and the slave antenna, and uses the baseline vector to solve the carrier attitude. Based on this, the measurement equation of the combination of multiple GNSS antennas and INS is established. This method relies on the correct fixation first. The baseline value after the ambiguity of the whole week belongs to the loose combination. The method in the present invention avoids this problem, and effectively utilizes the geometric configuration information and redundant observation information of multiple antennas with the tight combination of the observation value level, which can greatly improve the accuracy and reliability of the positioning and attitude determination results of the integrated system .

附图说明Description of drawings

图1是本申请实施例应用GNSS多天线与INS紧组合的定位定姿方法总体框架图。FIG. 1 is an overall framework diagram of a positioning and attitude determination method using a tight combination of GNSS multi-antenna and INS according to an embodiment of the present application.

具体实施方式Detailed ways

下面对本发明的实施例作详细说明，本实施例以本发明的技术方案为依据开展，给出了详细的实施方式和具体的操作过程，对本发明的技术方案作进一步解释说明。The embodiments of the present invention are described in detail below. This embodiment is carried out on the basis of the technical solutions of the present invention, and provides a detailed implementation manner and a specific operation process, and further explains the technical solutions of the present invention.

本实施例提供一种GNSS多天线与INS紧组合定位定姿方法，包括：This embodiment provides a GNSS multi-antenna and INS tight combination positioning and attitude determination method, including:

步骤1，通过杠杆臂向量和旋转矩阵，将多个GNSS天线的坐标参数归化到惯导中心坐标参数，构建以惯导中心处的位置向量和速度向量表示的双差载波和双差伪距率量测方程。Step 1: Through the lever arm vector and rotation matrix, the coordinate parameters of multiple GNSS antennas are normalized to the coordinate parameters of the inertial navigation center, and the double-difference carrier and double-difference pseudorange represented by the position vector and velocity vector at the inertial navigation center are constructed. rate measurement equation.

（1）对姿态矩阵进行初始化。考虑到短基线和模糊度固定的精度限制，采用GNSS天线的位置和速度对INS进行初始化，包括惯导中心的位置、速度和姿态，姿态包括了俯仰角、横滚角和航向角；无GNSS固定解时姿态的初始化采用INS自对准的方式，利用加速度调平计算俯仰角和横滚角，陀螺罗盘计算航向角。(1) Initialize the attitude matrix. Considering the short baseline and the accuracy limitations of fixed ambiguity, the INS is initialized with the position and velocity of the GNSS antenna, including the position, velocity and attitude of the inertial navigation center, and the attitude includes pitch, roll and heading angles; no GNSS The attitude initialization of the fixed solution adopts the INS self-alignment method, and uses the acceleration leveling to calculate the pitch angle and roll angle, and the gyro compass calculates the heading angle.

（2）通过杠杆臂向量和姿态矩阵分别建立每个天线的位置向量和速度向量与惯导中心处的位置向量和速度向量的关系：(2) The relationship between the position vector and velocity vector of each antenna and the position vector and velocity vector at the center of inertial navigation is established by lever arm vector and attitude matrix respectively:

Figure DEST_PATH_IMAGE095
Figure DEST_PATH_IMAGE095

其中，

Figure 113328DEST_PATH_IMAGE002
、
Figure 533945DEST_PATH_IMAGE003
分别表示天线
Figure 272094DEST_PATH_IMAGE004
的位置向量和速度向量，
Figure 903320DEST_PATH_IMAGE005
、v分别表示惯导中心 的位置向量和速度向量；
Figure 828550DEST_PATH_IMAGE006
是方向余弦矩阵，即旋转矩阵；
Figure 318438DEST_PATH_IMAGE007
是天线
Figure 645514DEST_PATH_IMAGE004
的杠杆臂向量；
Figure 929865DEST_PATH_IMAGE008
是IMU输出的三维角速度
Figure 607971DEST_PATH_IMAGE009
的反对称矩阵，即
Figure 167128DEST_PATH_IMAGE010
，
Figure 614290DEST_PATH_IMAGE011
为反对称阵算子；
Figure 69542DEST_PATH_IMAGE012
是地 球自转角速度的反对称矩阵。 in,
Figure 113328DEST_PATH_IMAGE002
,
Figure 533945DEST_PATH_IMAGE003
Respectively represent the antenna
Figure 272094DEST_PATH_IMAGE004
The position vector and velocity vector of ,
Figure 903320DEST_PATH_IMAGE005
, v represent the position vector and velocity vector of the inertial navigation center, respectively;
Figure 828550DEST_PATH_IMAGE006
is the direction cosine matrix, that is, the rotation matrix;
Figure 318438DEST_PATH_IMAGE007
is the antenna
Figure 645514DEST_PATH_IMAGE004
the lever arm vector of ;
Figure 929865DEST_PATH_IMAGE008
is the three-dimensional angular velocity output by the IMU
Figure 607971DEST_PATH_IMAGE009
The antisymmetric matrix of , that is
Figure 167128DEST_PATH_IMAGE010
,
Figure 614290DEST_PATH_IMAGE011
is an antisymmetric matrix operator;
Figure 69542DEST_PATH_IMAGE012
is the antisymmetric matrix of the angular velocity of the Earth's rotation.

（3）对于天线

Figure 969365DEST_PATH_IMAGE004
和天线
Figure DEST_PATH_IMAGE096
构成的短基线，两天线之间的电离层延迟和对流层延迟具 有较强的空间相关性，其双差载波方程和伪距率方程近似为： (3) For the antenna
Figure 969365DEST_PATH_IMAGE004
and antenna
Figure DEST_PATH_IMAGE096
The short baseline formed by the two antennas has a strong spatial correlation between the ionospheric delay and the tropospheric delay. The double-difference carrier equation and pseudorange rate equation are approximated as:

Figure DEST_PATH_IMAGE097
Figure DEST_PATH_IMAGE097

其中，

Figure 82945DEST_PATH_IMAGE014
为天线
Figure 384614DEST_PATH_IMAGE004
处星间单差视线向量矩阵，它是一个
Figure 10767DEST_PATH_IMAGE015
维的矩阵，
Figure 663465DEST_PATH_IMAGE016
为观测卫星总数减去所采用的卫星系统数；
Figure DEST_PATH_IMAGE098
为由相应载波波长构成的非零对角阵,
Figure 298846DEST_PATH_IMAGE019
表示天线
Figure 455021DEST_PATH_IMAGE020
之间双差观测噪声。另外，下标中的
Figure 517655DEST_PATH_IMAGE020
除了区别构成主天线和辅天线的 不同天线，还用于区分辅天线和基准站。即：当i表示主天线时，j表示基准站；当i表示辅天 线时，j表示主天线；
Figure 126491DEST_PATH_IMAGE046
表示下标
Figure 814830DEST_PATH_IMAGE047
之间的双差载波观测值，
Figure 91090DEST_PATH_IMAGE048
表示下标
Figure 59046DEST_PATH_IMAGE047
之间的双差伪 距率观测值； in,
Figure 82945DEST_PATH_IMAGE014
for the antenna
Figure 384614DEST_PATH_IMAGE004
single-difference line-of-sight vector matrix between the stars, which is a
Figure 10767DEST_PATH_IMAGE015
dimensional matrix,
Figure 663465DEST_PATH_IMAGE016
Subtract the number of satellite systems used from the total number of observation satellites;
Figure DEST_PATH_IMAGE098
is a non-zero diagonal matrix composed of corresponding carrier wavelengths,
Figure 298846DEST_PATH_IMAGE019
Indicates the antenna
Figure 455021DEST_PATH_IMAGE020
between double-difference observations. In addition, in the subscript
Figure 517655DEST_PATH_IMAGE020
In addition to distinguishing the different antennas that constitute the main antenna and the auxiliary antenna, it is also used to distinguish the auxiliary antenna from the reference station. That is: when i represents the main antenna, j represents the reference station; when i represents the secondary antenna, j represents the main antenna;
Figure 126491DEST_PATH_IMAGE046
Indicates subscript
Figure 814830DEST_PATH_IMAGE047
The double-difference carrier observations between,
Figure 91090DEST_PATH_IMAGE048
Indicates subscript
Figure 59046DEST_PATH_IMAGE047
The double-difference pseudorange rate between observations;

（4）将上述双差观测值中每个天线的位置向量和速度向量用惯导中心处的位置向量和速度向量代替，得到以惯导中心处的位置向量和速度向量表示的双差载波和双差伪距率方程：(4) Replace the position vector and velocity vector of each antenna in the above double-difference observations with the position vector and velocity vector at the center of inertial navigation, and obtain the double-difference carrier sum represented by the position vector and velocity vector at the center of inertial navigation. The double-difference pseudorange rate equation:

Figure DEST_PATH_IMAGE099
Figure DEST_PATH_IMAGE099

其中，

Figure DEST_PATH_IMAGE100
。 in,
Figure DEST_PATH_IMAGE100
.

步骤2，顾及多天线GNSS观测值组成双差观测值时的相关性，利用误差传播定律构建多天线观测协方差矩阵。Step 2: Considering the correlation when the multi-antenna GNSS observations form double-difference observations, the multi-antenna observation covariance matrix is constructed by using the law of error propagation.

考虑到多天线GNSS观测值组成双差观测值时的相关性，其协方差矩阵可以由误差传播定律得到，以n个天线和m颗卫星形成的站星间双差观测值为例：Considering the correlation of multi-antenna GNSS observations to form double-differenced observations, the covariance matrix can be obtained from the error propagation law. Taking the inter-station-to-satellite double-difference observations formed by n antennas and m satellites as an example:

Figure DEST_PATH_IMAGE101
Figure DEST_PATH_IMAGE101

其中，

Figure 951916DEST_PATH_IMAGE058
表示天线
Figure 929099DEST_PATH_IMAGE020
和卫星
Figure 59866DEST_PATH_IMAGE057
,
Figure 198724DEST_PATH_IMAGE045
之间的双差观测值，
Figure 47731DEST_PATH_IMAGE059
表示天线
Figure 313758DEST_PATH_IMAGE004
对卫星
Figure 564611DEST_PATH_IMAGE057
的原始观测值，
Figure 874370DEST_PATH_IMAGE061
为原始非差观测值与双差观测值之间的转化矩阵，
Figure 210673DEST_PATH_IMAGE061
的具体表达为： in,
Figure 951916DEST_PATH_IMAGE058
Indicates the antenna
Figure 929099DEST_PATH_IMAGE020
and satellite
Figure 59866DEST_PATH_IMAGE057
,
Figure 198724DEST_PATH_IMAGE045
The double-difference observations between ,
Figure 47731DEST_PATH_IMAGE059
Indicates the antenna
Figure 313758DEST_PATH_IMAGE004
to satellite
Figure 564611DEST_PATH_IMAGE057
the raw observations of ,
Figure 874370DEST_PATH_IMAGE061
is the transformation matrix between the original non-differenced observations and the double-differenced observations,
Figure 210673DEST_PATH_IMAGE061
The specific expression is:

Figure DEST_PATH_IMAGE102
Figure DEST_PATH_IMAGE102

其中，

Figure 529659DEST_PATH_IMAGE063
表示
Figure 635018DEST_PATH_IMAGE064
维的零矩阵，
Figure DEST_PATH_IMAGE103
的具体表达为： in,
Figure 529659DEST_PATH_IMAGE063
express
Figure 635018DEST_PATH_IMAGE064
dimensional zero matrix,
Figure DEST_PATH_IMAGE103
The specific expression is:

Figure DEST_PATH_IMAGE104
Figure DEST_PATH_IMAGE104

则GNSS多天线与INS紧组合系统观测值的协方差矩阵为：Then the covariance matrix of the observations of the GNSS multi-antenna and INS compact system is:

Figure DEST_PATH_IMAGE105
Figure DEST_PATH_IMAGE105

其中

Figure 912416DEST_PATH_IMAGE067
为原始非差观测值的协方差矩阵，它是一个对角阵，对角线元素值可根 据载噪比或者卫星高度角等确定。 in
Figure 912416DEST_PATH_IMAGE067
is the covariance matrix of the original non-differenced observations, which is a diagonal matrix, and the value of the diagonal elements can be determined according to the carrier-to-noise ratio or the altitude angle of the satellite.

步骤3，根据步骤1得到的双差载波和双差伪距率量测方程，以及步骤2得到的多天线观测协方差矩阵，构建得到GNSS多天线与INS紧组合系统的基于变量改正数的量测方程。Step 3, according to the double-difference carrier and double-difference pseudorange rate measurement equations obtained in step 1, and the multi-antenna observation covariance matrix obtained in step 2, construct and obtain the variable correction number based on the GNSS multi-antenna and INS compact combination system. measuring equation.

令GNSS多天线与INS紧组合系统的量测方程为：Let the measurement equation of the GNSS multi-antenna and INS compact system be:

Figure 204857DEST_PATH_IMAGE023
Figure 204857DEST_PATH_IMAGE023

其中，

Figure 573872DEST_PATH_IMAGE024
表示变量的改正数；
Figure 533737DEST_PATH_IMAGE025
表示组合系统的观测向量，包括GNSS系统观测 得到的载波和伪距率，
Figure 450878DEST_PATH_IMAGE026
表示GNSS观测向量
Figure 496194DEST_PATH_IMAGE025
的改正数，由GNSS的载波和伪距率与惯导 推算的载波和伪距率做差得到；
Figure 156983DEST_PATH_IMAGE027
表示状态向量
Figure 971355DEST_PATH_IMAGE028
的改正；
Figure 324976DEST_PATH_IMAGE029
表示两者之间的联系矩 阵；
Figure 592009DEST_PATH_IMAGE030
表示白噪声；
Figure 56489DEST_PATH_IMAGE031
表示GNSS多天线与INS紧组合系统观测值的协方差矩阵。 in,
Figure 573872DEST_PATH_IMAGE024
Represents the correction number of the variable;
Figure 533737DEST_PATH_IMAGE025
represents the observation vector of the combined system, including the carrier and pseudorange rates observed by the GNSS system,
Figure 450878DEST_PATH_IMAGE026
Represents a GNSS observation vector
Figure 496194DEST_PATH_IMAGE025
The correction number is obtained by the difference between the carrier and pseudorange rates of GNSS and the carrier and pseudorange rates estimated by inertial navigation;
Figure 156983DEST_PATH_IMAGE027
Represents a state vector
Figure 971355DEST_PATH_IMAGE028
correction;
Figure 324976DEST_PATH_IMAGE029
Represents the connection matrix between the two;
Figure 592009DEST_PATH_IMAGE030
represents white noise;
Figure 56489DEST_PATH_IMAGE031
Represents the covariance matrix of observations from a GNSS multi-antenna and INS compact system.

状态向量

Figure 725367DEST_PATH_IMAGE028
为：state vector
Figure 725367DEST_PATH_IMAGE028
for:

Figure DEST_PATH_IMAGE106
Figure DEST_PATH_IMAGE106

其中,

Figure 622DEST_PATH_IMAGE005
为载体三维位置向量，定为惯导参考中心；
Figure 20531DEST_PATH_IMAGE033
代表速度向量；
Figure 757542DEST_PATH_IMAGE034
代表姿态 向量；
Figure 546507DEST_PATH_IMAGE035
代表加速度计零偏；
Figure 241930DEST_PATH_IMAGE036
代表陀螺零偏；
Figure 483556DEST_PATH_IMAGE037
代表加速度计比例因子；
Figure 24259DEST_PATH_IMAGE038
代表陀 螺比例因子；
Figure 933309DEST_PATH_IMAGE039
代表双差模糊度向量。 in,
Figure 622DEST_PATH_IMAGE005
is the three-dimensional position vector of the carrier, and is set as the inertial navigation reference center;
Figure 20531DEST_PATH_IMAGE033
represents the velocity vector;
Figure 757542DEST_PATH_IMAGE034
represents the pose vector;
Figure 546507DEST_PATH_IMAGE035
represents the zero offset of the accelerometer;
Figure 241930DEST_PATH_IMAGE036
Represents the zero bias of the gyro;
Figure 483556DEST_PATH_IMAGE037
represents the accelerometer scale factor;
Figure 24259DEST_PATH_IMAGE038
represents the gyro scale factor;
Figure 933309DEST_PATH_IMAGE039
represents the double-difference ambiguity vector.

由于because

Figure DEST_PATH_IMAGE107
Figure DEST_PATH_IMAGE107

其中，

Figure 799634DEST_PATH_IMAGE024
表示变量的改正数，
Figure 528555DEST_PATH_IMAGE034
代表姿态向量，将上述步骤1得到的双差载波 和双差伪距率方程写成量测方程的形式，可得观测向量改正数与状态向量改正数之间的联 系矩阵为： in,
Figure 799634DEST_PATH_IMAGE024
represents the correction number of the variable,
Figure 528555DEST_PATH_IMAGE034
Representing the attitude vector, the double-difference carrier and double-difference pseudorange rate equations obtained in the above step 1 are written in the form of measurement equations, and the relationship matrix between the correction number of the observation vector and the correction number of the state vector can be obtained as:

Figure DEST_PATH_IMAGE108
Figure DEST_PATH_IMAGE108

其中，

Figure DEST_PATH_IMAGE109
表示
Figure 387796DEST_PATH_IMAGE053
维的零矩阵，
Figure DEST_PATH_IMAGE110
为由相应载波波长构成的非零对角 阵，
Figure DEST_PATH_IMAGE111
的具体表达如下： in,
Figure DEST_PATH_IMAGE109
express
Figure 387796DEST_PATH_IMAGE053
dimensional zero matrix,
Figure DEST_PATH_IMAGE110
is a non-zero diagonal matrix composed of the corresponding carrier wavelengths,
Figure DEST_PATH_IMAGE111
The specific expression is as follows:

Figure 416932DEST_PATH_IMAGE056
Figure 416932DEST_PATH_IMAGE056

在组成整个GNSS多天线与INS紧组合系统的量测方程时，通常选取一个移动天线 作为主天线，其余移动天线作为辅天线，主天线与基准站组成双差观测方程，为整个组合系 统提供位置基准，同时在主天线与辅天线之间组成双差观测方程，以达到对多天线冗余观 测值的充分利用。以p个移动天线（主天线下标为

Figure DEST_PATH_IMAGE112
，辅天线下标分别为
Figure DEST_PATH_IMAGE113
） 和q个基准站(下标分别为
Figure DEST_PATH_IMAGE114
)为例，其观测向量改正
Figure 719737DEST_PATH_IMAGE026
为： When forming the measurement equation of the entire GNSS multi-antenna and INS tight combination system, one mobile antenna is usually selected as the main antenna, and the other mobile antennas are used as auxiliary antennas. The main antenna and the base station form a double-difference observation equation to provide the position for the entire combined system. At the same time, a double-difference observation equation is formed between the main antenna and the auxiliary antenna, so as to make full use of the redundant observation values of multiple antennas. Take p mobile antennas (the main antenna is subscripted as
Figure DEST_PATH_IMAGE112
, the subscripts of the auxiliary antennas are
Figure DEST_PATH_IMAGE113
) and q base stations (the subscripts are
Figure DEST_PATH_IMAGE114
) as an example, its observation vector is corrected
Figure 719737DEST_PATH_IMAGE026
for:

Figure DEST_PATH_IMAGE115
Figure DEST_PATH_IMAGE115

其中，各元素表示下标对应天线之间的双差观测向量的改正数，如

Figure DEST_PATH_IMAGE116
表示 主天线
Figure 217845DEST_PATH_IMAGE112
与基准站
Figure DEST_PATH_IMAGE117
之间的双差载波观测向量改正数。Among them, each element represents the correction number of the double-difference observation vector between the subscripts corresponding to the antennas, such as
Figure DEST_PATH_IMAGE116
Indicates the main antenna
Figure 217845DEST_PATH_IMAGE112
with the base station
Figure DEST_PATH_IMAGE117
The double-difference carrier observation vector correction number between.

其联系矩阵

Figure 365930DEST_PATH_IMAGE029
为： its contact matrix
Figure 365930DEST_PATH_IMAGE029
for:

Figure DEST_PATH_IMAGE118
Figure DEST_PATH_IMAGE118

其中，各元素表示下标对应天线之间双差观测向量改正数与状态向量改正数之间的联系矩阵。Among them, each element represents the relationship matrix between the double-difference observation vector correction number and the state vector correction number between the corresponding antennas.

步骤4，以惯导中心处的位置、速度、姿态、传感器零偏及比例因子的误差作为系统的状态向量，采用一阶高斯马尔科夫过程建立系统状态方程。In step 4, the position, velocity, attitude, sensor bias and scale factor errors at the center of the inertial navigation are used as the state vector of the system, and a first-order Gaussian Markov process is used to establish the system state equation.

采用一阶高斯马尔科夫过程建立系统的状态方程为：The state equation of the system established by the first-order Gaussian Markov process is:

Figure 249572DEST_PATH_IMAGE069
Figure 249572DEST_PATH_IMAGE069

其中

Figure 192121DEST_PATH_IMAGE070
是状态向量
Figure 895634DEST_PATH_IMAGE028
在时刻
Figure 581831DEST_PATH_IMAGE071
的估值，上标“—”表示某变量还未被最新 的观测所更新，一旦更新，相应的上标将变为“+”；
Figure 54400DEST_PATH_IMAGE072
为状态转移矩阵，
Figure 167850DEST_PATH_IMAGE073
为过程噪声向 量，假定服从零均值正态分布，
Figure 358660DEST_PATH_IMAGE074
为过程噪声的转移矩阵，
Figure 582968DEST_PATH_IMAGE075
为过程噪声的协方差矩阵；
Figure 427820DEST_PATH_IMAGE072
的具体表达为： in
Figure 192121DEST_PATH_IMAGE070
is the state vector
Figure 895634DEST_PATH_IMAGE028
at the moment
Figure 581831DEST_PATH_IMAGE071
The superscript "—" indicates that a variable has not been updated by the latest observation. Once updated, the corresponding superscript will become "+";
Figure 54400DEST_PATH_IMAGE072
is the state transition matrix,
Figure 167850DEST_PATH_IMAGE073
is the process noise vector, assuming a zero-mean normal distribution,
Figure 358660DEST_PATH_IMAGE074
is the transition matrix of process noise,
Figure 582968DEST_PATH_IMAGE075
is the covariance matrix of the process noise;
Figure 427820DEST_PATH_IMAGE072
The specific expression is:

Figure DEST_PATH_IMAGE119
Figure DEST_PATH_IMAGE119

Figure 243330DEST_PATH_IMAGE077
Figure 243330DEST_PATH_IMAGE077

Figure 655856DEST_PATH_IMAGE078
Figure 655856DEST_PATH_IMAGE078

其中，

Figure 683855DEST_PATH_IMAGE005
为惯导中心的三维位置向量；
Figure 131017DEST_PATH_IMAGE079
为地表与地心之间的距离；
Figure 586269DEST_PATH_IMAGE080
为重力向 量；
Figure 751671DEST_PATH_IMAGE081
是加速度计输出的三维比力向量；
Figure 583361DEST_PATH_IMAGE009
为IMU输出的三维角速度；
Figure DEST_PATH_IMAGE120
为以
Figure 885029DEST_PATH_IMAGE009
为对角线元素的对角矩阵；
Figure DEST_PATH_IMAGE121
为相关时间，可以通过Allan方差分析得到；
Figure 527495DEST_PATH_IMAGE089
为加速度 计和陀螺测量时间间隔。 in,
Figure 683855DEST_PATH_IMAGE005
is the three-dimensional position vector of the inertial navigation center;
Figure 131017DEST_PATH_IMAGE079
is the distance between the surface and the center of the earth;
Figure 586269DEST_PATH_IMAGE080
is the gravity vector;
Figure 751671DEST_PATH_IMAGE081
is the three-dimensional specific force vector output by the accelerometer;
Figure 583361DEST_PATH_IMAGE009
is the three-dimensional angular velocity output by the IMU;
Figure DEST_PATH_IMAGE120
for
Figure 885029DEST_PATH_IMAGE009
is a diagonal matrix of diagonal elements;
Figure DEST_PATH_IMAGE121
is the correlation time, which can be obtained by Allan variance analysis;
Figure 527495DEST_PATH_IMAGE089
Measure time intervals for accelerometers and gyroscopes.

步骤5，对上述量测方程和系统状态方程进行卡尔曼滤波解算，得到惯导中心在每个时刻的位置、速度、姿态的最优估计。Step 5: Perform Kalman filter calculation on the above-mentioned measurement equation and system state equation to obtain the optimal estimation of the position, velocity and attitude of the inertial navigation center at each moment.

卡尔曼滤波具体递推计算公式为：The specific recursive calculation formula of Kalman filter is:

Figure 180193DEST_PATH_IMAGE091
Figure 180193DEST_PATH_IMAGE091

其中，

Figure 549994DEST_PATH_IMAGE070
是状态向量
Figure 971748DEST_PATH_IMAGE028
在时刻
Figure 34382DEST_PATH_IMAGE071
的估值，上标“—”表示某变量还未被最 新的观测所更新，一旦更新，相应的上标将变为“+”；
Figure 643218DEST_PATH_IMAGE072
为状态转移矩阵，
Figure 82290DEST_PATH_IMAGE073
为过程噪声向 量，假定服从零均值正态分布，
Figure 358550DEST_PATH_IMAGE074
为过程噪声的转移矩阵，
Figure DEST_PATH_IMAGE122
为过程噪声的协方差矩阵；
Figure 592085DEST_PATH_IMAGE092
是状态向量的方差协方差阵；
Figure 203064DEST_PATH_IMAGE093
为Kalman增益矩阵；
Figure 649089DEST_PATH_IMAGE031
为观测向量的方差协方差阵；
Figure 45435DEST_PATH_IMAGE094
为 单位阵。 in,
Figure 549994DEST_PATH_IMAGE070
is the state vector
Figure 971748DEST_PATH_IMAGE028
at the moment
Figure 34382DEST_PATH_IMAGE071
The superscript "—" indicates that a variable has not been updated by the latest observation. Once updated, the corresponding superscript will become "+";
Figure 643218DEST_PATH_IMAGE072
is the state transition matrix,
Figure 82290DEST_PATH_IMAGE073
is the process noise vector, assuming a zero-mean normal distribution,
Figure 358550DEST_PATH_IMAGE074
is the transition matrix of process noise,
Figure DEST_PATH_IMAGE122
is the covariance matrix of the process noise;
Figure 592085DEST_PATH_IMAGE092
is the variance-covariance matrix of the state vector;
Figure 203064DEST_PATH_IMAGE093
is the Kalman gain matrix;
Figure 649089DEST_PATH_IMAGE031
is the variance covariance matrix of the observation vector;
Figure 45435DEST_PATH_IMAGE094
is a unit matrix.

应用本发明GNSS多天线与INS紧组合的定姿定位方法，参考图1所示，完整流程为：Applying the method for determining the attitude and positioning of the GNSS multi-antenna and INS tight combination of the present invention, with reference to Fig. 1, the complete process is as follows:

（1）对多天线原始观测值进行处理，选取一个移动天线作为主天线，其余移动天线作为辅天线，主天线与基准站组成双差观测方程，同时在主天线与辅天线之间组成双差观测方程，进行相对定位解算，得到组合系统的位置、速度和姿态信息。(1) Process the original observations of multiple antennas, select one mobile antenna as the main antenna, and the other mobile antennas as auxiliary antennas. The main antenna and the base station form a double-difference observation equation, and at the same time form a double-difference between the main antenna and the auxiliary antenna. Observing the equation, performing relative positioning calculation, and obtaining the position, velocity and attitude information of the combined system.

（2）利用上述信息对IMU进行初始对准，当无GNSS固定解时，姿态的初始化采用INS自对准的方式，利用加速度调平计算俯仰角和横滚角，陀螺罗盘计算航向角。(2) Use the above information to initially align the IMU. When there is no GNSS fixed solution, the attitude initialization adopts the INS self-alignment method, and the acceleration leveling is used to calculate the pitch angle and roll angle, and the gyro compass calculates the heading angle.

（3）根据IMU的输出对组合系统的位置、速度和姿态信息进行更新，首先利用陀螺输出的角速率完成姿态更新，在此基础上对加速度计输出的比力进行坐标转换，最后更新速度和位置信息。(3) Update the position, velocity and attitude information of the combined system according to the output of the IMU. First, use the angular rate output by the gyro to complete the attitude update, and on this basis, perform coordinate transformation on the specific force output by the accelerometer. location information.

（4）组成卡尔曼滤波的状态方程和量测方程，根据上述更新后的位置、速度和姿态信息和卫星星历推算载波和伪距率观测值，并与GNSS系统观测得到的载波和伪距率做差作为量测方程中观测向量的改正数，同时通过杠杆臂向量和姿态矩阵分别建立每个天线的位置向量和速度向量与惯导中心处的位置向量和速度向量的关系，将量测方程中每个天线的位置向量和速度向量用惯导中心处的位置向量和速度向量代替，以惯导中心处的位置、速度、姿态、传感器零偏等误差作为系统的状态向量，同时采用一阶高斯马尔科夫过程建立系统的状态方程。(4) The state equation and measurement equation of the Kalman filter are composed, and the carrier and pseudorange rate observations are calculated according to the above updated position, velocity and attitude information and satellite ephemeris, and are compared with the carrier and pseudorange observed by the GNSS system. The rate difference is used as the correction number of the observation vector in the measurement equation. At the same time, the relationship between the position vector and velocity vector of each antenna and the position vector and velocity vector at the center of the inertial navigation is established through the lever arm vector and the attitude matrix. The position vector and velocity vector of each antenna in the equation are replaced by the position vector and velocity vector at the inertial navigation center, and the position, velocity, attitude, sensor bias and other errors at the inertial navigation center are used as the state vector of the system, and a The order Gaussian Markov process establishes the equation of state of the system.

（5）进行卡尔曼滤波解算，得到每个时刻惯导中心的位置、速度、姿态的最优估计，同时传感器误差进行反馈校正。(5) Perform Kalman filter calculation to obtain the optimal estimation of the position, velocity and attitude of the inertial navigation center at each moment, and at the same time, the sensor error is corrected by feedback.

综上所述，本发明的GNSS多天线与INS紧组合定位定姿新方法利用了多天线的几何构型信息和冗余观测信息，缓解了单天线GNSS在恶劣环境下，卫星可见性及信号抗干扰能力差的问题，提高了系统的容错性，相对于传统基于模糊度固定基线的GNSS/INS定位定姿方法，充分地利用了观测数据，有利于提高集成系统的稳定性和精度。To sum up, the new method for positioning and attitude determination of GNSS multi-antenna and INS tight combination of the present invention utilizes the geometric configuration information and redundant observation information of the multi-antenna, and relieves the single-antenna GNSS in harsh environments. The problem of poor anti-interference ability improves the fault tolerance of the system. Compared with the traditional GNSS/INS positioning and attitude determination method based on the fixed baseline of ambiguity, the observation data is fully utilized, which is beneficial to improve the stability and accuracy of the integrated system.

以上实施例为本申请的优选实施例，本领域的普通技术人员还可以在此基础上进行各种变换或改进，在不脱离本申请总的构思的前提下，这些变换或改进都应当属于本申请要求保护的范围之内。The above embodiments are the preferred embodiments of the application, and those of ordinary skill in the art can also carry out various transformations or improvements on this basis. Without departing from the general concept of the application, these transformations or improvements should belong to the present application. within the scope of the application for protection.

Claims (10)

1. A GNSS multi-antenna and INS tightly combined positioning and attitude determination method is characterized by comprising the following steps:
step 1, normalizing coordinate parameters of a plurality of GNSS antennas to inertial navigation center coordinate parameters through a lever arm vector and a rotation matrix, and constructing a double-difference carrier and double-difference pseudo-range rate measurement equation represented by a position vector and a speed vector at an inertial navigation center;
step 2, considering the correlation when the multi-antenna GNSS observation values form a double-difference observation value, and constructing a covariance matrix of the observation values of the GNSS multi-antenna and INS tight combination system by using an error propagation law;
step 3, according to the double difference carrier wave and double difference pseudo range rate measurement equation obtained in the step 1 and the multi-antenna observation covariance matrix obtained in the step 2, a variable correction-based measurement equation of the GNSS multi-antenna and INS tight combination system is constructed;
step 4, taking the position, the speed, the attitude, the zero offset of the sensor and the error of the scale factor at the inertial navigation center as the state vector of the system, and establishing a system state equation by adopting a first-order Gaussian Markov process;
and 5, performing Kalman filtering solution according to the measurement equation and the system state equation to obtain the optimal estimation of the position, the speed and the attitude of the inertial navigation center at each moment.
2. The method for positioning and determining the attitude according to claim 1, wherein the method for normalizing the coordinate parameters of the GNSS antenna to the inertial navigation center coordinate parameters comprises:
Figure 654149DEST_PATH_IMAGE001
wherein,
Figure 894637DEST_PATH_IMAGE002
、
Figure 845275DEST_PATH_IMAGE003
respectively represent an antenna
Figure 960737DEST_PATH_IMAGE004
The position vector and the velocity vector of (c),
Figure 47641DEST_PATH_IMAGE005
、
Figure 724610DEST_PATH_IMAGE006
respectively representing a position vector and a velocity vector of an inertial navigation center;
Figure 896966DEST_PATH_IMAGE007
is a direction cosine matrix, i.e. a rotation matrix;
Figure 317583DEST_PATH_IMAGE008
is an antenna
Figure 524573DEST_PATH_IMAGE004
The lever arm vector of (a);
Figure 372444DEST_PATH_IMAGE009
is the three-dimensional angular velocity of the IMU output
Figure 32095DEST_PATH_IMAGE010
Of antisymmetric matrices, i.e.
Figure 256403DEST_PATH_IMAGE011
，
Figure 317900DEST_PATH_IMAGE012
Is an antisymmetric array operator;
Figure 336671DEST_PATH_IMAGE013
is an antisymmetric matrix of the rotational angular velocity of the earth.
3. A method for positioning and attitude determination according to claim 2, characterized by constructing double-differenced carrier and double-differenced pseudorange rate equations expressed in a position vector and a velocity vector at the inertial navigation center, expressed as:
Figure 483619DEST_PATH_IMAGE014
wherein,
Figure 246039DEST_PATH_IMAGE015
is an antenna
Figure 162042DEST_PATH_IMAGE004
A single difference sight line vector matrix between the satellites, which is
Figure 617294DEST_PATH_IMAGE016
A matrix of dimensions is formed by a matrix of dimensions,
Figure 251538DEST_PATH_IMAGE017
subtracting the number of adopted satellite systems from the total number of the observation satellites;
Figure 552069DEST_PATH_IMAGE018
the amount of the carbon dioxide is the intermediate amount,
Figure 322579DEST_PATH_IMAGE019
is a non-zero diagonal matrix of corresponding carrier wavelengths,
Figure 948732DEST_PATH_IMAGE020
antenna for representation
Figure 804693DEST_PATH_IMAGE021
Double difference observation noise between the two;
Figure 174495DEST_PATH_IMAGE022
antenna for representation
Figure 65090DEST_PATH_IMAGE021
The two-difference carrier observations in between,
Figure 596566DEST_PATH_IMAGE023
antenna for representation
Figure 205401DEST_PATH_IMAGE021
Double-differenced pseudorange rate observations between.
4. The method according to claim 1, wherein the measurement equation obtained in step 3 is expressed as:
Figure 346271DEST_PATH_IMAGE024
wherein,
Figure 91373DEST_PATH_IMAGE025
a correction number representing a variable;
Figure 59329DEST_PATH_IMAGE026
the observation vector representing the GNSS antenna combined system comprises the carrier wave and the pseudo range rate observed by the GNSS system,
Figure 889881DEST_PATH_IMAGE027
representing observation vectors
Figure 335906DEST_PATH_IMAGE026
The correction number of the inertial navigation system is obtained by subtracting the carrier wave and the pseudo-range rate of the GNSS and the carrier wave and the pseudo-range rate calculated by the inertial navigation system;
Figure 466673DEST_PATH_IMAGE028
representing state vectors
Figure DEST_PATH_IMAGE029
Correcting;
Figure 74372DEST_PATH_IMAGE030
indicating the number of corrections of the observation vector
Figure 392221DEST_PATH_IMAGE027
And state vector correction number
Figure 641937DEST_PATH_IMAGE031
A matrix of relationships between;
Figure 361631DEST_PATH_IMAGE032
representing white noise;
Figure 671390DEST_PATH_IMAGE033
a covariance matrix representing the observations;
state vector
Figure 210955DEST_PATH_IMAGE029
Comprises the following steps:
Figure 998783DEST_PATH_IMAGE034
wherein,
Figure 838563DEST_PATH_IMAGE005
determining a carrier three-dimensional position vector as an inertial navigation reference center;
Figure 53644DEST_PATH_IMAGE035
representing a velocity vector;
Figure 346085DEST_PATH_IMAGE036
representing a pose vector;
Figure 937603DEST_PATH_IMAGE037
represents the accelerometer zero offset;
Figure 631890DEST_PATH_IMAGE038
represents a gyro zero bias;
Figure 283451DEST_PATH_IMAGE039
representing an accelerometer scale factor;
Figure 797609DEST_PATH_IMAGE040
represents a gyro scale factor;
Figure 927239DEST_PATH_IMAGE041
representing a double-difference ambiguity vector.
5. The method according to claim 4, wherein the number of correction of the observation vector is set to be equal to or less than the number of correction of the observation vector
Figure 741611DEST_PATH_IMAGE042
And a matrix of associations
Figure 298494DEST_PATH_IMAGE043
Respectively as follows:
Figure 821921DEST_PATH_IMAGE044
wherein,
Figure 755242DEST_PATH_IMAGE045
a main antenna is shown, which is,
Figure DEST_PATH_IMAGE046
to represent
Figure 158541DEST_PATH_IMAGE047
A plurality of auxiliary antennas are arranged on the base plate,
Figure 151905DEST_PATH_IMAGE048
represent
Figure 640655DEST_PATH_IMAGE049
A reference station;
Figure 377667DEST_PATH_IMAGE050
denotes a subscript
Figure 901052DEST_PATH_IMAGE051
A double-difference carrier observation between them,
Figure 799738DEST_PATH_IMAGE052
denotes a subscript
Figure 775785DEST_PATH_IMAGE051
Double-difference pseudorange rate observations between;
Figure DEST_PATH_IMAGE053
represent a subscript
Figure 50908DEST_PATH_IMAGE051
A matrix of associations between the associated double difference vector corrections and the state vector corrections, and having:
Figure 163221DEST_PATH_IMAGE054
wherein,
Figure 498387DEST_PATH_IMAGE055
to represent
Figure DEST_PATH_IMAGE056
The zero matrix of the dimension(s) is,
Figure 696150DEST_PATH_IMAGE057
for non-zero diagonal arrays of corresponding carrier wavelengths, intermediate variables
Figure 40544DEST_PATH_IMAGE058
The specific expression of (A) is as follows:
Figure 7363DEST_PATH_IMAGE059
。
6. the method according to claim 4, wherein the covariance matrix of observations of the GNSS multi-antenna and INS tightly-combined system is constructed using the law of error propagation
Figure 513430DEST_PATH_IMAGE060
The method comprises the following steps:
first, the antenna is mounted
Figure 198490DEST_PATH_IMAGE021
And satellite
Figure 80995DEST_PATH_IMAGE061
,
Figure 167900DEST_PATH_IMAGE049
The double difference between the observed values is expressed as
Figure 844869DEST_PATH_IMAGE062
Will be arbitrary antenna
Figure 282803DEST_PATH_IMAGE004
For any satellite
Figure 201956DEST_PATH_IMAGE061
Is expressed as
Figure 408946DEST_PATH_IMAGE063
Then the original non-differential observed value can be obtained from the law of error propagation
Figure 256816DEST_PATH_IMAGE063
And double difference observed value
Figure 916468DEST_PATH_IMAGE062
The relationship between them is:
Figure 140776DEST_PATH_IMAGE064
then, the conversion matrix is calculated according to the relational expression
Figure 202273DEST_PATH_IMAGE065
Comprises the following steps:
Figure 955465DEST_PATH_IMAGE066
wherein,
Figure 367992DEST_PATH_IMAGE067
represent
Figure 130411DEST_PATH_IMAGE068
A zero matrix of the dimensions is formed,
Figure 46415DEST_PATH_IMAGE069
the specific expression of (A) is as follows:
Figure 501667DEST_PATH_IMAGE070
finally, according to the transformation matrix
Figure 870331DEST_PATH_IMAGE065
And covariance matrix of original non-differential observations
Figure 170863DEST_PATH_IMAGE071
Determining the covariance matrix of the combined system observations as:
Figure 472531DEST_PATH_IMAGE072
。
7. the method according to claim 4, wherein the correction of the observation vector is performed by using a correction value of the observation vector
Figure 833105DEST_PATH_IMAGE042
The calculation method comprises the following steps: calculating a double-difference carrier observation value and a double-difference pseudo range rate observation value according to position, speed and attitude information updated by an inertial navigation center and a satellite ephemeris; and (4) subtracting the carrier wave and the pseudo-range rate obtained by the observation of the GNSS system, and using the difference as the correction number of the observation vector in the measurement equation.
8. The method according to claim 7, wherein a double-difference observation equation is constructed from the observed values of the multi-antenna GNSS system, and the position, velocity and attitude information of the combined system is obtained by solving the double-difference observation equation and is used for initializing the inertial navigation center.
9. The method according to claim 4, wherein the state equation of the system is established by a first-order Gaussian Markov process as follows:
Figure 689066DEST_PATH_IMAGE073
wherein,
Figure 58867DEST_PATH_IMAGE074
is a state vector
Figure 949463DEST_PATH_IMAGE029
At the moment of time
Figure 480938DEST_PATH_IMAGE075
The superscript "-" indicates that a variable has not been updated by the most recent observation, and once updated, the corresponding superscript will become "+";
Figure 89774DEST_PATH_IMAGE076
in order to be a state transition matrix,
Figure 732108DEST_PATH_IMAGE077
for the process noise vector, assuming a normal distribution following zero mean,
Figure 477210DEST_PATH_IMAGE078
is a transfer matrix for the process noise,
Figure 445166DEST_PATH_IMAGE079
a covariance matrix that is the process noise;
Figure 275719DEST_PATH_IMAGE076
the specific expression of (A) is as follows:
Figure 987323DEST_PATH_IMAGE080
Figure 85467DEST_PATH_IMAGE081
Figure 958745DEST_PATH_IMAGE082
wherein the state transition matrix
Figure 276594DEST_PATH_IMAGE076
The 3 in each element subscript in (a) represents a 3 x 3 matrix,nto representn×n3 is generatednRepresents 3 linesnA matrix of the columns is formed,nx 3 representsnA matrix of rows and 3 columns;
Figure 791889DEST_PATH_IMAGE005
the three-dimensional position vector of the inertial navigation center is obtained;
Figure 246004DEST_PATH_IMAGE083
the distance between the earth surface and the earth center;
Figure 290183DEST_PATH_IMAGE084
is a gravity vector;
Figure 360907DEST_PATH_IMAGE085
is the three-dimensional specific force vector output by the accelerometer;
Figure 148735DEST_PATH_IMAGE010
a three-dimensional angular velocity output for the IMU;
Figure 722936DEST_PATH_IMAGE086
to be composed of
Figure 203596DEST_PATH_IMAGE010
A diagonal matrix being diagonal elements;
Figure 496037DEST_PATH_IMAGE087
are all intermediate variables;
Figure 87555DEST_PATH_IMAGE088
for the correlation time of the accelerometer zero-offset,
Figure 516262DEST_PATH_IMAGE089
is the relative time of the gyro zero offset,
Figure 167823DEST_PATH_IMAGE090
is the correlation time of the accelerometer scale factor,
Figure 947561DEST_PATH_IMAGE091
the correlation time of the gyro scale factor;
Figure 77191DEST_PATH_IMAGE092
measuring time intervals for the accelerometer and gyroscope;
Figure 625984DEST_PATH_IMAGE093
is an antisymmetric matrix of the rotational angular velocity of the earth;
the specific recursion calculation formula for performing Kalman filtering solution according to the measurement equation and the system state equation is as follows:
Figure 182867DEST_PATH_IMAGE094
wherein,
Figure 449900DEST_PATH_IMAGE095
is a state directionVariance covariance matrix of the quantities;
Figure 383221DEST_PATH_IMAGE096
is a Kalman gain matrix;
Figure 52100DEST_PATH_IMAGE060
a variance covariance matrix for the observation vector;
Figure 779884DEST_PATH_IMAGE097
is a unit array.
10. An electronic device comprising a memory and a processor, the memory having stored therein a computer program, wherein the computer program, when executed by the processor, causes the processor to implement the method of any of claims 1-9.
CN202211024566.9A 2022-08-25 2022-08-25 A GNSS multi-antenna and INS tightly combined positioning and attitude determination method and device Active CN115096303B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211024566.9A CN115096303B (en) 2022-08-25 2022-08-25 A GNSS multi-antenna and INS tightly combined positioning and attitude determination method and device

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211024566.9A CN115096303B (en) 2022-08-25 2022-08-25 A GNSS multi-antenna and INS tightly combined positioning and attitude determination method and device

Publications (2)

Publication Number Publication Date
CN115096303A true CN115096303A (en) 2022-09-23
CN115096303B CN115096303B (en) 2022-11-22

Family

ID=83300374

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211024566.9A Active CN115096303B (en) 2022-08-25 2022-08-25 A GNSS multi-antenna and INS tightly combined positioning and attitude determination method and device

Country Status (1)

Country Link
CN (1) CN115096303B (en)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115900527A (en) * 2023-01-06 2023-04-04 中南大学 Deformation Monitoring Method Based on GNSS System Error Recursive Semiparametric Modeling
CN117289322A (en) * 2023-11-24 2023-12-26 江苏领创星通卫星通信科技有限公司 High-precision positioning algorithm based on IMU, GPS and UWB
CN117388900A (en) * 2023-12-13 2024-01-12 深圳大学 A method for constructing a GNSS/INS joint ocean dynamic base station
CN119618207A (en) * 2025-02-11 2025-03-14 华航导控(天津)科技有限公司 Self-adaptive combined navigation method suitable for multiple scenes
CN120008643A (en) * 2025-04-17 2025-05-16 中南大学 A navigation and positioning method based on Kalman filtering and related equipment
CN120721128A (en) * 2025-08-14 2025-09-30 毫厘智能科技(江苏)有限公司 Antenna arm calibration method, device, and storage medium
US12474484B2 (en) 2022-11-11 2025-11-18 Guilin University Of Technology POS high-accuracy positioning method under occlusion
CN120991915A (en) * 2025-10-24 2025-11-21 中交第二公路勘察设计研究院有限公司 A method and system for estimating error parameters of a sensor
CN120991915B (en) * 2025-10-24 2026-02-06 中交第二公路勘察设计研究院有限公司 Error parameter estimation method and system of sensor

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6496778B1 (en) * 2000-09-14 2002-12-17 American Gnc Corporation Real-time integrated vehicle positioning method and system with differential GPS
US20050234644A1 (en) * 2004-04-17 2005-10-20 Ching-Fang Lin Positioning and navigation method and system thereof
GB0723810D0 (en) * 2006-12-05 2008-01-16 Boeing Co Ultra-tightly coupled global navigation satellite system space borne receiver system
US20140152493A1 (en) * 2012-11-30 2014-06-05 Applanix Inc. Quasi tightly coupled gnss-ins integration process
US20170350988A1 (en) * 2014-12-26 2017-12-07 Furuno Electric Co., Ltd. State calculating device, method of calculating state, and state calculating program
CN111288990A (en) * 2020-03-19 2020-06-16 云南电网有限责任公司电力科学研究院 Combined attitude measurement method for overhead maintenance robot

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6496778B1 (en) * 2000-09-14 2002-12-17 American Gnc Corporation Real-time integrated vehicle positioning method and system with differential GPS
US20050234644A1 (en) * 2004-04-17 2005-10-20 Ching-Fang Lin Positioning and navigation method and system thereof
GB0723810D0 (en) * 2006-12-05 2008-01-16 Boeing Co Ultra-tightly coupled global navigation satellite system space borne receiver system
US20140152493A1 (en) * 2012-11-30 2014-06-05 Applanix Inc. Quasi tightly coupled gnss-ins integration process
US20170350988A1 (en) * 2014-12-26 2017-12-07 Furuno Electric Co., Ltd. State calculating device, method of calculating state, and state calculating program
CN111288990A (en) * 2020-03-19 2020-06-16 云南电网有限责任公司电力科学研究院 Combined attitude measurement method for overhead maintenance robot

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
YONG WANG 等: "The Influence of Attitude Dilution of Precision on the Observable Degree and Observability Analysis With Different Numbers of Visible Satellites in a Multi-Antenna GNSS/INS Attitude Determination System", 《IEEE ACCESS》 *
任钊 等: "复杂环境下GNSS/INS紧组合垂直振动监测精度评估", 《导航定位学报》 *
刘猛奎: "多天线GNSS/INS组合导航算法研究", 《中国优秀硕士学位论文全文数据库基础科学辑》 *
田源: "GNSS/INS定位测姿模型构建与算法研究", 《中国优秀硕士学位论文全文数据库基础科学辑》 *

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US12474484B2 (en) 2022-11-11 2025-11-18 Guilin University Of Technology POS high-accuracy positioning method under occlusion
CN115900527A (en) * 2023-01-06 2023-04-04 中南大学 Deformation Monitoring Method Based on GNSS System Error Recursive Semiparametric Modeling
CN117289322A (en) * 2023-11-24 2023-12-26 江苏领创星通卫星通信科技有限公司 High-precision positioning algorithm based on IMU, GPS and UWB
CN117388900A (en) * 2023-12-13 2024-01-12 深圳大学 A method for constructing a GNSS/INS joint ocean dynamic base station
CN117388900B (en) * 2023-12-13 2024-03-08 深圳大学 GNSS/INS combined ocean dynamic reference station construction method
CN119618207A (en) * 2025-02-11 2025-03-14 华航导控(天津)科技有限公司 Self-adaptive combined navigation method suitable for multiple scenes
CN120008643A (en) * 2025-04-17 2025-05-16 中南大学 A navigation and positioning method based on Kalman filtering and related equipment
CN120008643B (en) * 2025-04-17 2025-07-01 中南大学 Navigation positioning method based on Kalman filtering and related equipment
CN120721128A (en) * 2025-08-14 2025-09-30 毫厘智能科技(江苏)有限公司 Antenna arm calibration method, device, and storage medium
CN120991915A (en) * 2025-10-24 2025-11-21 中交第二公路勘察设计研究院有限公司 A method and system for estimating error parameters of a sensor
CN120991915B (en) * 2025-10-24 2026-02-06 中交第二公路勘察设计研究院有限公司 Error parameter estimation method and system of sensor

Also Published As

Publication number Publication date
CN115096303B (en) 2022-11-22

Similar Documents

Publication Publication Date Title
CN115096303A (en) GNSS multi-antenna and INS tightly-combined positioning and attitude determination method and equipment
CN113203418B (en) GNSSINS visual fusion positioning method and system based on sequential Kalman filtering
US6496778B1 (en) Real-time integrated vehicle positioning method and system with differential GPS
US5543804A (en) Navagation apparatus with improved attitude determination
EP0856747A1 (en) Method and apparatus for attitude determination utilizing an inertial measurement unit and a plurality of satellite transmitters
CA2273158C (en) Attitude determination method and system
CN114167472B (en) INS-assisted GNSS PPP precise dynamic navigation positioning method and system
US6697736B2 (en) Positioning and navigation method and system thereof
CN109613585A (en) A real-time direction finding method for ultra-short baseline GNSS dual-antenna base station antennas
CN110133692B (en) Inertial navigation technology-assisted high-precision GNSS dynamic inclination measurement system and method
CN110779521A (en) Multi-source fusion high-precision positioning method and device
CN107037469A (en) Based on the self-alignment double antenna combined inertial nevigation apparatus of installation parameter
CN111965685A (en) A low-orbit satellite/inertial integrated navigation and positioning method based on Doppler information
CN108873034A (en) A kind of implementation method of inertial navigation subcarrier ambiguity resolution
De Agostino et al. Performances comparison of different MEMS-based IMUs
CN117804443A (en) A Beidou satellite multi-mode fusion positioning monitoring method and system
CA2716597A1 (en) Apparatus and method for driftless attitude determination and reliable localization of vehicles
CN118443024A (en) Self-adaptive robust filtering method and system for RTK/INS (real-time kinematic/inertial navigation system) tightly integrated navigation
CN118688839A (en) A multi-level joint enhanced positioning method and system assisted by a low-orbit satellite constellation
CN112697154A (en) Self-adaptive multi-source fusion navigation method based on vector distribution
CN114895340B (en) Positioning method and device of dual-antenna GNSS/INS integrated navigation system
CN116879927B (en) Ship satellite compass heading determination method based on three-antenna collinear common clock architecture
CN114994730B (en) A High-Precision Positioning Method and System Based on Carrier Phase-Based Satellite-Inertial Integrated Navigation
CN115657095A (en) Beidou high-precision positioning method of self-adaptive weight matrix
CN115507848B (en) A positioning method based on wheel speed inertial navigation joint pre-integration and tight integration with RTK

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant