Efficient dynamic and static information fusion method for airborne all-source satellite-inertial tight integration navigation

Abstract

The present application relates to an efficient dynamic and static information fusion method for airborne all-source satellite-inertial tightly combined navigation. The method combines the improved information filtering method with the tightly combined navigation method, and reduces the high-dimensional Kalman filter measurement update to two low-dimensional fusion processes, namely, the dynamic measurement update of the carrier phase observation subset and the static measurement update of the pseudorange and pseudorange rate observation subset, and uses the dynamic and static information filtering form to ensure high computational efficiency under the condition of global optimal estimation; the information fusion contribution rate of GNSS observation anomalies is properly suppressed through the dynamic and static filtering new information adaptive weighting strategy; the method effectively solves the problems of low computational efficiency and long computational time of tightly combined Kalman measurement updates caused by high data rate and high-dimensional GNSS observations, and overcomes the adverse effects of GNSS gross errors on navigation accuracy, and improves the real-time performance and accuracy of airborne combined navigation calculations under the condition of limited airborne computing resources.

Description

Technical Field

[0001] The present application relates to the technical field of inertial / satellite combined navigation, and in particular to an efficient dynamic and static information fusion method for airborne all-source satellite-inertial tight combined navigation. Background Art

[0002] Due to limited load capacity and energy supply, existing small unmanned aerial vehicles cannot be equipped with traditional high-precision autonomous navigation systems such as high-precision inertial navigation systems, precision approach radars, and astronomical navigation star sensors. The resulting lack of autonomous navigation capabilities has largely caused difficulties in autonomous flight technology. Existing engineering solutions tend to use cost-effective airborne satellite navigation (Global Navigation Satellite System, GNSS) modules and visual sensors to cooperate with remote control links to complete manual navigation and control of unmanned aerial vehicles. On the one hand, this solution makes it impossible to carry out flight missions under severe weather conditions with poor visibility such as dense fog, heavy rain and snow, greatly reducing the operating efficiency of unmanned aerial vehicles; on the other hand, this also places high demands on the quality of radio links and the operating skills of ground personnel, increasing the workload of ground stations. In recent years, the emergence of high-performance Micro-Electro-Mechanical System (MEMS) strapdown inertial navigation system (SINS) has brought a turning point for enhancing the autonomous navigation capability of small unmanned aerial vehicles. MEMS-SINS is a miniature, low-power, low-cost sensor integration system based on silicon chips. It can be directly installed on small unmanned aerial vehicles without significantly increasing the weight, size and energy consumption of the airborne navigation system. It forms a miniaturized, high-performance satellite / inertial combined navigation system with the existing GNSS module, providing accurate and reliable navigation information such as position, speed, attitude, angular velocity and angular velocity in real time, enhancing the autonomous navigation and autonomous flight capabilities of small unmanned aerial vehicles.

[0003] Satellite / inertial integrated navigation systems can generally be divided into three architectures: loosely coupled, tightly coupled, and ultra-tightly coupled (UTC). Although ultra-tight coupling has shown superior performance under high carrier dynamics and weak GNSS signal conditions, its design and development involves GNSS receiver baseband signal processing and is usually not open to commercial modules. Therefore, loose coupling and tight coupling are widely used in airborne navigation. The difference between the two lies mainly in the form of GNSS measurement data fusion. Loose coupling uses the position and velocity information provided by the GNSS navigation module to fuse with inertial navigation, while tight coupling directly uses GNSS raw measurement data (such as pseudorange and Doppler shift measurements) for deeper information fusion. Tight coupling is not affected by satellite geometry, partial GNSS signal interruption (i.e., the number of visible satellites is less than 4), and cascade filter correlation, and is usually superior to loose coupling in terms of navigation accuracy.

[0004] In engineering, the satellite / inertial combined navigation system uses the extended Kalman filter method to achieve data fusion. With the full construction and operation of the Beidou navigation satellite system, the existing multi-GNSS system has brought about a surge in the number of available observations, and formed high-dimensional GNSS observation information based on multi-system, multi-frequency, and multi-type GNSS observations (referred to as all-source GNSS measurements), which greatly improved the accuracy performance that can be achieved by traditional satellite / inertial combined navigation. However, the price that comes with it is the huge amount of calculations brought about by the inversion of the high-dimensional observation matrix for Kalman filter measurement updates. At the same time, high-speed and high-dynamic carriers (such as small UAVs) usually require high data update rates for navigation calculations and outputs, resulting in a sharp increase in the amount of high-dimensional matrix multiplication and inversion operations in a short period of time, further exacerbating the burden of onboard computing. Unfortunately, the scarce onboard computing resources of small unmanned aerial vehicles lead to lower navigation computing efficiency, causing significant computing delays under high-dimensional filtering conditions, threatening the real-time performance of onboard navigation and flight safety. Therefore, in order not to lose the potential real-time navigation accuracy, it is necessary to improve the efficiency of the combined navigation algorithm under high-dimensional GNSS observation conditions, improve the real-time performance of the calculation, and provide airborne computing power margin for computing tasks such as GNSS precise relative positioning, integrity risk estimation, and high sampling rate strapdown inertial solution. Summary of the invention

[0005] Based on this, it is necessary to provide an efficient dynamic and static information fusion method for airborne all-source satellite-inertial tight integrated navigation to address the above technical problems.

[0006] An airborne all-source satellite-inertial tight integrated navigation efficient dynamic and static information fusion method, the method comprising:

[0007] Build an airborne all-source satellite-inertial integrated navigation system. The integrated navigation system acquires the aircraft's inertial navigation measurement data and multi-system, multi-frequency, and multi-type GNSS navigation measurement data in real time during the carrier's motion.

[0008] Perform inertial navigation calculation based on inertial navigation measurement data, establish system error state equation, and then perform Kalman time update;

[0009] Search the synchronous ephemeris for the GNSS navigation measurement data according to the GNSS observation data time, check the validity of the observation data and synchronous ephemeris, and determine whether the GNSS navigation measurement data is received normally;

[0010] If it is judged to be abnormal, the strapdown inertial navigation result data is directly output;

[0011] If it is judged to be normal, a high-dimensional GNSS observation vector is constructed according to the GNSS navigation measurement data, and the high-dimensional GNSS observation vector is divided into a carrier phase observation subset and a pseudorange and pseudorange rate observation subset; a robust dynamic and static filtering algorithm based on successive fusion of independent observation information is constructed to perform information fusion of satellite / inertial tightly integrated navigation, and calculate the navigation error, use the navigation error estimate to correct the navigation result, and output the final navigation result data; the robust dynamic and static filtering algorithm based on successive fusion of independent observation information means that in any epoch, robust dynamic Kalman filtering is performed based on the dynamic model and the carrier phase observation subset, and then the pseudorange and pseudorange rate observation subsets are sequentially added to perform robust static Kalman filtering to obtain the global optimal fusion solution.

[0012] The above-mentioned airborne all-source satellite-inertial tightly integrated navigation efficient dynamic and static information fusion method combines the improved information filtering method and the tightly integrated navigation method, reduces the high-dimensional Kalman filter measurement update to two low-dimensional fusion processes, namely, the dynamic measurement update of the carrier phase observation subset and the static measurement update of the pseudorange and pseudorange rate observation subset, and uses the dynamic and static information filtering form to ensure high computational efficiency under the condition of global optimal estimation; the information fusion contribution rate of GNSS observation anomalies is properly suppressed through the dynamic and static filtering new information adaptive weighting strategy; the method effectively solves the problems of low computational efficiency and long computational time of tightly combined Kalman measurement updates brought about by high data rate and high-dimensional GNSS observations, while overcoming the adverse effects of GNSS gross errors on navigation accuracy, and improving the real-time and accuracy of airborne integrated navigation calculations under the condition of limited airborne computing resources. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Figure 1 It is a flowchart of an efficient dynamic and static information fusion method of airborne all-source satellite-inertial tight integrated navigation in one embodiment;

[0014] Figure 2 It is a flow chart of an efficient information fusion method for airborne all-source satellite / inertial tight integrated navigation in one embodiment;

[0015] Figure 3 A three-dimensional motion trajectory diagram of an unmanned aerial vehicle in one embodiment;

[0016] Figure 4 A comparison diagram of the time taken for dynamic and static information filtering, information filtering and sequential Kalman filtering measurement update in one embodiment;

[0017] Figure 5 A schematic diagram of the relationship between the airborne test integrated navigation measurement update calculation time and the number of GNSS observation values ​​in one embodiment;

[0018] Figure 6 It is a curve diagram of the three-dimensional positioning error of the satellite / inertial combined navigation of the unmanned aerial vehicle in one embodiment. DETAILED DESCRIPTION

[0019] In order to make the purpose, technical solution and advantages of the present application more clearly understood, the present application is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

[0020] In one embodiment, Figure 1 As shown, a method for efficient dynamic and static information fusion of airborne all-source satellite-inertial tight integrated navigation is provided, and the method comprises the following steps:

[0021] Step 100: Build an airborne all-source satellite-inertial integrated navigation system. The integrated navigation system acquires the aircraft's inertial navigation measurement data and multi-system, multi-frequency, and multi-type GNSS navigation measurement data in real time during the carrier's motion.

[0022] Specifically, inertial navigation measurement data and GNSS navigation measurement data are the core information sources of the satellite / inertial combined navigation system. The all-source GNSS navigation measurement information includes satellite broadcast ephemeris data corresponding to multiple GNSS systems (such as GPS, BDS, GLONASS, GALILEO, etc.) and multi-frequency and multi-type satellite navigation observation data (such as Beidou B1 / B2 / B3 triple-frequency pseudorange, pseudorange change rate, carrier phase observation values, GPS L1 / L2 / L5 triple-frequency pseudorange, pseudorange change rate, carrier phase observation values, etc.).

[0023] Inertial navigation measurement data include GNSS time scale, apparent acceleration (or velocity increment) and angular velocity (or angle increment); GNSS navigation measurement data include satellite broadcast ephemeris data corresponding to multiple GNSS systems and multi-frequency and multi-type satellite navigation observation data.

[0024] Step 102: Perform inertial navigation solution according to the inertial navigation measurement data, establish the system error state equation, and then perform Kalman time update;

[0025] Specifically, the information matrix, information error state vector and two auxiliary calculation matrices of the integrated navigation system are defined. According to the discrete system error state equation, the information related to the system error state at the previous moment is used to predict the information error state vector and information matrix at the current moment. This process is called the time update equation of information filtering.

[0026] The system error state equation is established based on the navigation error of the carrier where the tightly integrated navigation system is located, and the inertial navigation solution provides real-time prediction values ​​of the current motion parameters of the carrier.

[0027] Step 104: searching for synchronous ephemeris for the GNSS navigation measurement data according to the GNSS observation data time, performing validity check on the observation data and synchronous ephemeris, and determining whether the GNSS navigation measurement data is received normally.

[0028] Specifically, the received GNSS navigation measurement data is checked for data validity, including ephemeris data synchronization check, observation data satellite number check, observation data frequency point completeness check, etc. Under the condition of ensuring correct reception of observation data and synchronous ephemeris data, multi-frequency multi-system GNSS pseudorange, pseudorange change rate and carrier phase observation equations are established, and then Kalman measurement is performed to update and calculate navigation errors to achieve integrated navigation information fusion.

[0029] Step 106: If it is judged to be abnormal, the strapdown inertial navigation result data is directly output.

[0030] Step 108: If it is judged to be normal, a high-dimensional GNSS observation vector is constructed according to the GNSS navigation measurement data, and the high-dimensional GNSS observation vector is divided into a carrier phase observation sub-set and a pseudorange and pseudorange rate variation observation sub-set; a robust dynamic and static filtering algorithm based on successive fusion of independent observation information is constructed to perform information fusion of satellite / inertial tight integrated navigation, and calculate the navigation error, and use the navigation error estimate to correct the navigation result, and output the final navigation result data; the robust dynamic and static filtering algorithm based on the successive fusion of independent observation information means that at any epoch, robust dynamic Kalman filtering is performed based on the dynamic model and the carrier phase observation sub-set, and then the pseudorange and pseudorange rate variation observation sub-set is sequentially added to perform robust static Kalman filtering to obtain the global optimal fusion solution.

[0031] Specifically, multi-system, multi-frequency and multi-type GNSS navigation measurement data mainly include Beidou B1 / B2 / B3 three-frequency pseudorange, pseudorange rate, carrier phase observations, GPS L1 / L2 / L5 three-frequency pseudorange, pseudorange rate, carrier phase observations, etc., which are constructed from GNSS observations of different systems, frequencies and types to form a high-dimensional GNSS observation vector. Pseudorange and carrier phase observations of different systems and frequencies are independently measured by different hardware loops, and it can be considered that there is no correlation between the observations. Therefore, the high-dimensional GNSS observation vector can be divided into two subsets according to statistical correlation, namely "carrier phase observation quantum set" and "pseudorange and pseudorange rate observation quantum set". In this way, the high-dimensional GNSS observation vector can be reduced in dimension to improve the efficiency of subsequent information processing.

[0032] The tightly combined measurement equation is based on the GNSS navigation measurement equation of pseudorange and carrier phase. In order to achieve single-station independent navigation without relying on the reference station, multi-system multi-frequency epoch inter-satellite differential carrier phase observations and their observation equations are often used in engineering.

[0033] Based on the calculation results of inertial navigation and time update, the inter-epoch inter-satellite differential carrier phase measurement information is introduced to measure and update the information error state vector and information matrix, and the process of obtaining the information error state vector and information matrix after measurement update is the measurement update process of information filtering. In the measurement update process, the inertial navigation dynamics prediction information and redundant carrier phase observation information are fully utilized to construct the new information test statistic to effectively detect GNSS observation anomalies, and the weight adaptive adjustment method is used to optimize the Kalman filter estimation effect to obtain the optimal state estimation value of the navigation error.

[0034] In order to eliminate the systematic error of pseudorange measurement, dual-frequency ionospheric-free combined pseudorange observations and their observation equations are often used in engineering. Based on the inertial navigation prediction and GNSS pseudorange and pseudorange rate observations at the current moment, an adaptive static information filter with robustness is used to perform information fusion of satellite / inertial tight integrated navigation and calculate the navigation error. Since the inertial navigation dynamics prediction information has been used in the dynamic information filtering stage, in order to avoid repeated use, the remaining GNSS pseudorange and pseudorange rate observation information needs to be statically filtered. At this time, the state prediction information provided by the error state equation is no longer used, but the information error state vector and information matrix obtained by the previous step of dynamic information filtering are used as the state prediction information, which is the measurement update process of static information filtering. Specifically, in the measurement update process, the static equation information and redundant pseudorange and pseudorange rate observation information are fully utilized to construct the relevant new information test statistics to effectively detect GNSS observation anomalies, and the weight adaptive adjustment method is used to optimize the Kalman filter estimation effect, and finally the optimal state estimate of the navigation error is obtained. If the inter-epoch inter-satellite differential carrier phase measurement information is not introduced for measurement update in the previous step, then at this time, based on the inertial navigation and time update calculation results, the multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate measurement information is directly introduced to measure and update the information error state vector and information matrix to obtain the information error state vector and information matrix after measurement update, which is the measurement update process degenerated into dynamic information filtering.

[0035] The navigation error estimate is used to correct the navigation result and output the final navigation result data. It can be understood as follows: the optimal estimate of the system error state vector at the current moment is used to correct the navigation result (i.e., carrier position, speed, attitude and other parameters) at the current moment.

[0036] In the above-mentioned airborne all-source satellite-inertial tightly integrated navigation efficient dynamic and static information fusion method, the method combines the improved information filtering method and the tightly integrated navigation method, and reduces the high-dimensional Kalman filter measurement update to two low-dimensional fusion processes, namely, the dynamic measurement update of the carrier phase observation subset and the static measurement update of the pseudorange and pseudorange rate observation subset, and uses the dynamic and static information filtering form to ensure high computational efficiency under the condition of global optimal estimation; the information fusion contribution rate of GNSS observation anomalies is properly suppressed through the dynamic and static filtering new information adaptive weighting strategy; the method effectively solves the problems of low computational efficiency and long computational time of tightly combined Kalman measurement updates brought about by high data rate and high-dimensional GNSS observations, and overcomes the adverse effects of GNSS gross errors on navigation accuracy, and improves the real-time and accuracy of airborne integrated navigation calculations under the condition of limited airborne computing resources.

[0037] In one embodiment, step 102 includes: constructing an error state equation based on inertial navigation measurement data with position error, velocity error, misalignment angle, gyro bias and table bias as state quantities:

[0038] (1)

[0039] in, represents the position error, represents the speed error, represents the misalignment angle, Indicates the gyro bias, Indicates the addition of zero bias. and They represent the first-order derivatives of position error, velocity error, misalignment angle, gyro bias and table bias respectively. Represents the angular velocity of the Earth's rotation projected on the "North-East-Earth" navigation coordinate system; express The angular velocity of the system relative to the earth coordinate system is projected in the "north-east-earth" navigation coordinate system; ; is the cross product matrix; Represents the ground velocity vector projected on the "North-East-Earth" navigation coordinate system; Represents the transformation matrix from the system to the navigation system; , and Represents three intermediate error vectors, which can be expressed as , and Function of represents the apparent acceleration vector in the body coordinate system; and are the random noises measured by the gyro and the meter, respectively, and are described by a first-order Gauss-Markov process, where is unit Gaussian white noise, and represents the first-order Gauss-Markov model coefficients obtained by inertial navigation calibration, subscript and They represent the addition table and the top respectively;

[0040] The error state vector is , the error state equation is transformed into a discrete error state equation:

[0041] (2)

[0042] in, express The system error state vector at time ; express kSystem error state vector at time -1; Indicates from Time to k The state transfer matrix at time; represents the system noise transfer matrix; express White noise vector at time t;

[0043] Defining the Information Matrix is the error state covariance matrix The inverse matrix of , the information error state vector is To achieve the information filtering time update, the first auxiliary calculation matrix is ​​defined and the second auxiliary calculation matrix , the first auxiliary calculation matrix and the second auxiliary calculation matrix They are:

[0044] (3)

[0045] (4)

[0046] in, represents the optimal estimate of the information matrix obtained at the last moment, that is, , represents the system process noise covariance matrix, ;

[0047] According to the information error state vector, information matrix estimation value, first auxiliary calculation matrix and second auxiliary calculation matrix at the previous moment, the information error state vector and information matrix at the current moment are predicted to perform Kalman information filter time update, namely:

[0048] (5)

[0049] (6)

[0050] in, Represents the information matrix forecast value obtained at the current moment, that is, ; represents the identity matrix, represents the information matrix, , represents the information error state vector at the previous moment, , Represents the information error state vector at the current moment.

[0051] In one embodiment, a robust dynamic and static filtering algorithm based on successive fusion of independent observation information is constructed to perform information fusion of satellite / inertial tight integrated navigation, calculate navigation error, use navigation error estimation value to correct navigation result, and output final navigation result data, including the following steps:

[0052] Step 200: Determine whether the carrier phase observation quantum set is empty;

[0053] Step 202: if the carrier phase observation quantum set is a non-empty set, then a multi-system multi-frequency epoch inter-satellite differential carrier phase observation equation is established according to the carrier phase observation quantum set; according to the inertial navigation prediction quantity and the GNSS carrier phase observation quantity at the current moment, an anti-error adaptive dynamic information filter is used to perform information fusion of satellite / inertial tight integrated navigation, and the navigation error is calculated; it is determined whether the pseudorange and pseudorange rate observation quantum set is an empty set; if the pseudorange and pseudorange rate observation quantum set is an empty set, the navigation result is corrected by using the navigation error estimation value, and the final navigation result data is output; if the pseudorange and pseudorange rate observation quantum set is a non-empty set, then a multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate observation equation is established; based on the inertial navigation prediction quantity and the GNSS pseudorange and pseudorange rate observation quantity at the current moment, an anti-error adaptive static information filter is used to perform information fusion of satellite / inertial tight integrated navigation, and the navigation error is calculated, the navigation result is corrected by using the navigation error estimation value, and the final navigation result data is output.

[0054] Step 204: if the carrier phase observation quantum set is an empty set, determine whether the pseudorange and pseudorange rate observation quantum set is an empty set; if the pseudorange and pseudorange rate observation quantum set is an empty set, directly output the navigation result; if the pseudorange and pseudorange rate observation quantum set is a non-empty set, establish a multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate observation equation; based on the inertial navigation prediction quantity and the GNSS pseudorange and pseudorange rate observation quantity at the current moment, use an anti-error adaptive static information filter to perform satellite / inertial tight integrated navigation information fusion, and calculate the navigation error; use the navigation error estimate to correct the navigation result, and output the final navigation result data.

[0055] In one embodiment, a high-dimensional GNSS observation vector is constructed based on GNSS navigation measurement data, and the high-dimensional GNSS observation vector is divided into a carrier phase observation sub-set and a pseudorange and pseudorange rate observation sub-set, including:

[0056] The first The corresponding pseudorange, pseudorange change rate and carrier phase measurements of GNSS satellite No. , and ; Subscript symbol " " indicates the GNSS system and frequency to which the observed value belongs. Taking the Beidou B1 frequency measurement value as an example, Inter-satellite differential pseudo-range observations corresponding to the B1 frequency point of Beidou satellite , inter-satellite differential pseudo-range rate observation and inter-epoch intersatellite differential carrier phase observations for:

[0057] (7)

[0058] in, Indicates The corresponding pseudo-range observation value of the B1 frequency point corresponding to the GNSS satellite, Indicates The pseudo-range rate observation of the B1 frequency point corresponding to the GNSS satellite, Indicates r The pseudo-range observation value of the B1 frequency point corresponding to the GNSS satellite No. Indicates r The pseudo-range rate observation of the B1 frequency point corresponding to the GNSS satellite No. Indicates the reference satellite number. Usually the GNSS satellite with the largest altitude angle is selected as the reference satellite. represents the time difference operator, , Respectively represent n The B1 frequency point corresponding to the GNSS satellite is and The GNSS carrier phase measurement value at time, and Indicates the reception time of two consecutive GNSS carrier phase measurement signals. , Indicates the number of available satellites; Respectively represent n No. and r The GNSS carrier phase time difference measurement value corresponding to the B1 frequency point of the GNSS satellite.

[0059] The first The calculation formulas for the inter-satellite differential pseudorange, inter-satellite differential pseudorange change rate and inter-epoch inter-satellite differential carrier phase observations corresponding to the B1 frequency point of the Beidou satellite No. 1 are as follows:

[0060] (8)

[0061] in, , and They are respectively The inter-satellite differential pseudo-range, inter-satellite differential pseudo-range variation rate and inter-epoch inter-satellite differential carrier phase observations corresponding to the B1 frequency point of the BeiDou satellite No. and Indicates The position and velocity vector of the GNSS satellite, and Indicates The position and velocity vector of the reference satellite can be obtained by broadcast ephemeris calculation; and is the position and velocity vector of the phase center of the receiver GNSS antenna predicted by the inertial navigation; Indicates that from The unit direction vector from the satellite position to the receiver position, Indicates r The unit vector from the satellite position to the receiver position, and , calculated from the broadcast ephemeris and the inertial navigation predicted position; Indicates the inertial navigation speed; Inertial navigation forecast Time from To The direction cosine matrix of the system; The arm vector from the GNSS antenna center to the inertial navigation center is System projection; the construction process of pseudorange, pseudorange change rate and carrier phase observations of a certain frequency point of other GNSS systems is similar to the construction of the three types of observations of the Beidou B1 frequency point mentioned above, and will not be repeated here.

[0062] The Kalman filter high-dimensional GNSS observation vector is The Kalman filter high-dimensional GNSS observation vector constructed with the observations of the five frequencies of the Beidou and GPS dual systems (i.e., B1, B2, B3, L1, and L2 frequencies) is:

[0063] (9)

[0064] in, represents the Kalman filtered high-dimensional GNSS observation vector; represents the inter-satellite differential pseudo-range observation column vector, subscript Indicates the type of GNSS system used, The value of is BDS or GPS, which is used to indicate the type of GNSS system used; and denote the inter-satellite differential pseudorange variation rate of each frequency point and the inter-epoch inter-satellite differential carrier phase observation column vector, respectively. Indicates the type of GNSS system frequency used. The value of is B1, B2, B3, L1 or L2, which is used to indicate the type of GNSS system frequency used.

[0065] Taking Beidou B1 frequency as an example, the column vector , and The calculation formulas are:

[0066] (10)

[0067] (11)

[0068] (12)

[0069] in, and They represent the dual-frequency ionospheric-free combined pseudorange observations of Beidou B1 / B3 and GPS L1 / L2 respectively. It should be noted that the pseudorange observation is a linear combination of the dual-frequency observations of Beidou and GPS to eliminate the large ionospheric delay error in the pseudorange. Therefore, the multi-frequency pseudorange measurement values ​​of each system are actually combined into a pseudorange observation value without ionospheric delay error. and ; Respectively expressed as r Satellite No. 1 is the reference satellite, and satellite No. 1, No. 2, and No. m use Beidou's dual-frequency ionospheric elimination combined pseudo-range measurement values. Respectively expressed as r Satellite No. is the reference satellite corresponding to No. 1, No. 2, m The satellite uses BeiDou's dual-frequency ionospheric-free combined pseudo-range observation. , Respectively expressed as r Satellite No. is the reference satellite corresponding to No. 1, No. 2, m The inter-satellite differential pseudo-range rate measurement value of satellite No. 1 at the B1 frequency point, Respectively expressed as r Satellite No. is the reference satellite corresponding to No. 1, No. 2, m The inter-satellite differential pseudo-range rate observation of satellite No. 1 at the B1 frequency point, Respectively expressed as r Satellite No. is the reference satellite corresponding to No. 1, No. 2, m The inter-epoch inter-satellite differential carrier phase measurement value of satellite No. at the B1 frequency point, Respectively expressed as r Satellite No. is the reference satellite corresponding to No. 1, No. 2, m The inter-epoch inter-satellite differential carrier phase observations of satellite No. 1 at the B1 frequency point.

[0070] In one embodiment, a multi-system multi-frequency epoch inter-satellite differential carrier phase observation equation is established based on a carrier phase observation quantum set, including:

[0071] Taking Beidou B1 frequency as an example, the inter-satellite differential carrier phase observation equation between single-system single-frequency epochs is constructed. The construction process of the inter-satellite differential carrier phase observation equation between single-system single-frequency epochs includes:

[0072] In obtaining the The inter-satellite differential carrier phase observations corresponding to the epoch Based on this observation and the combined navigation error state The relationship between them is:

[0073] (13)

[0074] in, represents the inter-satellite differential carrier phase corresponding to the epoch of the nth satellite, Indicates the relative position of GNSS satellite n r The inter-satellite difference vector of the satellite-to-ground unit direction vector of the GNSS satellite No. n The unit vector pointing from the satellite position to the receiver position With r The unit vector pointing from the satellite position to the receiver position difference, is the inter-epoch inter-satellite differential carrier phase measurement noise; and They are:

[0075] (14)

[0076] in, Indicates the inertial navigation sampling time interval; Indicates time interval Number of internal inertial navigation measurement data; and They represent the constant matrices used to extract velocity error and attitude error from the error state vector respectively; represents the system dynamics model matrix, Indicates t k The SINS attitude matrix corresponding to the moment, Indicates t k-1 The SINS attitude matrix corresponding to the moment, It represents the projection of the arm vector from the user SINS center to the GNSS satellite antenna center in the carrier coordinate system. Represents the identity matrix.

[0077] The inter-satellite differential carrier phase observation vector between epochs of Beidou B1 frequency point can be obtained. Combined navigation error status The relationship between them is:

[0078] (15)

[0079] (16)

[0080] (17)

[0081] in, represents the inter-satellite differential carrier phase observation vector between epochs of BeiDou B1 frequency point, Indicates the integrated navigation error status; represents the noise vector of the inter-satellite differential carrier phase observation between epochs of the BeiDou B1 frequency point, , , They represent the inter-satellite epoch differential carrier phase measurement noise corresponding to the 1st, 2nd, and mth GNSS satellites with the rth GNSS satellite as the reference satellite, respectively. They represent the satellite-to-ground unit direction vectors and the inter-satellite difference vectors of the 1st, 2nd, and mth GNSS satellites respectively, with the rth GNSS satellite as the reference satellite. Represents the design matrix of inter-satellite differential carrier phase observation between epochs at the B1 frequency point.

[0082] Based on the construction process of the inter-satellite differential carrier phase observation equation between single-system single-frequency epochs, the inter-satellite differential carrier phase observation equation between multi-system multi-frequency epochs can be constructed; specifically, the inter-epoch inter-satellite differential carrier phase observation equation between five frequencies (i.e., B1, B2, B3, L1, and L2) of the Beidou and GPS dual systems is taken as an example. Formula (14) is used to establish the inter-epoch inter-satellite differential carrier phase observation matrices for the B1, B2, and B3 frequencies of the Beidou system respectively. , and ; Then use formula (14) to establish the inter-epoch inter-satellite differential carrier phase observation matrix of the GPS system L1 and L2 frequency points respectively: and The inter-satellite differential carrier phase observation equation for multi-system multi-frequency epochs is constructed as:

[0083] (18)

[0084] (19)

[0085] in, represents the multi-system multi-frequency epoch inter-satellite differential carrier phase observation vector, represents the noise vector of the inter-satellite differential carrier phase observation between multi-system and multi-frequency epochs, They represent the inter-satellite differential carrier phase observation vectors of the BeiDou B2 and B3 frequency points respectively. They represent the inter-satellite differential carrier phase observation vectors between epochs of GPS L1 and L2 frequency points, They represent the inter-epoch inter-satellite differential carrier phase observation design matrices for the B2 and B3 frequency points, Denote the inter-epoch inter-satellite differential carrier phase observation design matrices of L1 and L2 frequency points, They represent the inter-satellite differential carrier phase observation noise vectors of BeiDou B2 and B3 frequency points respectively. They represent the inter-satellite differential carrier phase observation noise vectors between epochs of GPS L1 and L2 frequency points, Design matrices for multi-system multi-frequency inter-epoch intersatellite differential carrier phase observations.

[0086] In one embodiment, the process of robust adaptive dynamic information filtering includes: using the current moment inertial navigation prediction and GNSS carrier phase observations , and then detect GNSS observation anomalies channel by channel, and then use the variance matching method to perform weighted information fusion on the GNSS observation anomalies; , observation matrix and the observation noise variance matrix By dividing into blocks, we can get:

[0087] (20)

[0088] (twenty one)

[0089] in, Indicates the current moment GNSS channel observations; Indicates the corresponding The row vector of the observation design matrix corresponding to the channels; Indicates The channel observation noise variance is , represents the number of available observations at the current moment, Represents the design matrix of multi-system multi-frequency inter-epoch inter-satellite differential carrier phase observation.

[0090] Current time The anomaly test statistic of GNSS observations corresponding to the channels for:

[0091] (twenty two)

[0092] in, Indicates GNSS channel updates, , Represents the row vector of the observation design matrix corresponding to the i-th channel at the current moment, Indicates the current SINS forecast error state value; Indicates GNSS channel innovation variance, , represents the SINS forecast error state covariance matrix at the current moment, For the corresponding i The column vector of the observation design matrix corresponding to the channels.

[0093] like , then judge the The GNSS channels are observed normally, and the settings ;in Represents the standard normal distribution Quantile;

[0094] like , then judge the GNSS channel observation anomalies, set the diagonal elements of the weight matrix for:

[0095] (twenty three)

[0096] in, are two preset constants, ;

[0097] Traverse all After observing channels, construct the weight matrix ; Use the fusion formula to perform information filtering and information fusion to calculate the navigation error; the fusion formula is:

[0098] (twenty four)

[0099] (25)

[0100] (26)

[0101] (27)

[0102] in, , , and Respectively represent the current time after the dynamic information filter measurement is updated The information error state vector, information matrix, error state covariance matrix and error state vector, Indicates the current time k The information error state vector state prediction information, Represents the information matrix forecast value obtained at the current moment.

[0103] In one embodiment, the step of establishing multi-system multi-frequency inter-satellite differential pseudorange and pseudorange change rate observation equations includes:

[0104] GNSS pseudorange observations usually use dual-frequency ionosphere-eliminating combined pseudorange observations. Taking the Beidou single system B1 and B3 frequency point ionosphere-eliminating combined pseudorange observations as an example, after obtaining the first Inter-satellite differential B1 / B3 dual-frequency ionospheric-free combined pseudo-range observations corresponding to satellite Based on the pseudo-range observation Combined navigation error status The relationship between them is:

[0105] (28)

[0106] in, represents the inter-satellite differential pseudo-range measurement noise, It represents the Beidou B1 / B3 dual-frequency ionospheric-free combined pseudorange observation. It represents the inter-satellite difference vector of the satellite-to-ground unit direction vector of BDS satellite n relative to BDS satellite r, represents the direction cosine matrix from the local navigation coordinate system to the Earth-centered Earth-fixed system (can be provided by SINS), represents the attitude matrix output by SINS, It represents the projection of the arm vector from the GNSS antenna center to the inertial navigation center in the b system, represents the BDS dual-frequency linear combination pseudorange measurement noise, represents the SINS attitude misalignment angle vector, Represents the SINS navigation position error vector.

[0107] Taking the BeiDou B1 frequency point as an example, a single-system single-frequency inter-satellite differential pseudo-range rate observation equation is constructed. A single-system single-frequency inter-satellite differential pseudo-range rate observation equation is constructed at the BeiDou B1 frequency point; Inter-satellite differential pseudo-range rate observations corresponding to satellite Combined navigation error status The relationship between them is:

[0108] (29)

[0109] in, is the inter-satellite differential pseudorange rate measurement noise; It represents the inter-satellite differential pseudo-range rate measurement value of the B1 frequency point corresponding to the nth Beidou satellite, with the rth Beidou satellite as the reference satellite. Represents the projection vector of the angular velocity vector of the local navigation coordinate system relative to the inertial coordinate system in the local navigation coordinate system, Represents the projection vector of the angular velocity vector of the carrier coordinate system measured by the gyro relative to the inertial coordinate system in the body coordinate system, represents the SINS velocity error vector in the local navigation coordinate system, represents the SINS gyro bias error vector, Represents the noise vector of the inter-satellite differential pseudorange rate measurement at the B1 frequency point.

[0110] The BeiDou B1 / B3 dual-frequency intersatellite differential pseudo-range observation vector can be obtained. , B1 frequency point inter-satellite differential pseudo-range rate observation vector Combined navigation error status The relationships between them are:

[0111] (30)

[0112] (31)

[0113] (32)

[0114] (33)

[0115] (34)

[0116] in, represents the BDS dual-frequency linear combination pseudo-range measurement value vector, represents the BDS dual-frequency linear combination pseudorange measurement noise vector, Respectively represent No. 1, No. 2, m The inter-satellite difference vector of the satellite-to-ground unit direction vector of BDS satellite No. 1 relative to BDS satellite No. 2, represents the direction cosine matrix from the local navigation coordinate system to the Earth-centered Earth-fixed system, represents the attitude matrix output by SINS, It represents the projection of the arm vector from the GNSS antenna center to the inertial navigation center in the b system, Respectively represent No. 1, No. 2, m BDS satellite relative to r The BDS dual-frequency linear combination pseudo-range measurement noise of BDS satellite No. represents the inter-satellite differential pseudorange rate observation matrix of the B1 frequency point, represents the inter-satellite differential pseudo-range rate observation noise at the B1 frequency point, Respectively represent No. 1, No. 2, m BDS satellite relative to rThe noise vector of the inter-satellite differential pseudo-range rate observation of BDS satellite No. 1 at the B1 frequency point;

[0117] Based on the steps of constructing the single-system single-frequency inter-satellite differential pseudorange rate observation equation at the BeiDou B1 frequency point, the multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate observation equations are established. Specifically, taking the BeiDou and GPS dual-system five-frequency inter-satellite differential pseudorange and pseudorange rate observation equations (i.e., B1, B2, B3, L1 and L2 frequency points) as an example, the BeiDou B1 / B3 ionospheric elimination combined pseudorange observation matrix is ​​established using formulas (29) to (31) respectively , and the inter-satellite differential pseudorange rate observation matrix of B1, B2 and B3 frequency points , and ; Then use formula (29) to formula (31) to establish the GPS L1 / L2 ionospheric elimination combined pseudorange observation matrix , and the inter-satellite differential pseudo-range rate observation matrix of L1 and L2 frequency points and The observation equations for constructing multi-system multi-frequency inter-satellite differential pseudorange and pseudorange variation rate are:

[0118] (35)

[0119] (36)

[0120] in, represents the multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate observation vector, represents the multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate observation design matrix, represents the multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate observation error vector, represents the BeiDou B1 / B3 ionospheric-free combined pseudorange observation matrix, , and They represent the inter-satellite differential pseudorange rate observation matrices of B1, B2 and B3 frequency points respectively, represents the GPS L1 / L2 ionospheric-free combined pseudorange observation matrix, and Respectively represent the inter-satellite differential pseudorange rate observation matrices of L1 and L2 frequency points; represents the BeiDou B1 / B3 ionospheric-free combined pseudorange observation error, and represent the inter-satellite differential pseudorange rate observation errors of B1, B2 and B3 frequency points respectively, represents the GPS L1 / L2 ionospheric-free combined pseudorange observation error, and They represent the inter-satellite differential pseudo-range rate observation errors of L1 and L2 frequency points, represents the BeiDou B1 / B3 dual-frequency inter-satellite differential pseudo-range observation vector, represents the GPS dual-frequency inter-satellite differential pseudo-range observation vector, , They represent the inter-satellite differential pseudorange rate observation vectors of B1, B2, B3, L1 and L2 frequency points respectively.

[0121] In one embodiment, the process of robust adaptive static information filtering includes: if the "carrier phase observation quantum set" is not empty and the robust adaptive dynamic information filtering calculation step has been completed, then the inertial navigation prediction quantity at the current moment is used and GNSS pseudorange and pseudorange rate observations , perform robust adaptive static information filtering.

[0122] First, GNSS observation anomalies are detected channel by channel, and then the variance matching method is used to reduce the weight of GNSS observation anomalies. The state vector obtained by the first step of dynamic filtering is used as the state prediction information for static filtering; that is: , .

[0123] After using the observations, observation matrix and noise matrix related to GNSS pseudorange and pseudorange rate observations, we set , construct the weight matrix ;in, represents the observation vector, represents the observation matrix, represents the multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate observation vector, Represents the design matrix of multi-system multi-frequency inter-satellite differential pseudorange and pseudorange variability observations.

[0124] A fusion formula is used to fuse information and calculate the navigation error; the fusion formula is shown in formula (24) to formula (27).

[0125] In one embodiment, Figure 2 As shown, a method for efficient dynamic and static information fusion of airborne all-source satellite-inertial tight integrated navigation is provided, comprising the following steps:

[0126] Step 1: Build an airborne all-source satellite / inertial integrated navigation system. The integrated navigation system acquires the aircraft's inertial navigation measurement data and GNSS navigation measurement data in real time during the carrier's motion. The inertial navigation measurement data includes GNSS time scale, apparent acceleration (or velocity increment) and angular velocity (or angle increment); the GNSS navigation measurement data includes satellite broadcast ephemeris data corresponding to multiple GNSS systems and multi-frequency and multi-type satellite navigation observation data.

[0127] Step 2: Perform inertial navigation solution based on airborne inertial navigation measurement data, establish the system error state equation, implement Kalman filter time update, and go to step 3.

[0128] Step 3: Based on the received GNSS navigation measurement data, search for the synchronous ephemeris according to the GNSS observation data time, check the validity of the observation data and the synchronous ephemeris, and determine whether the GNSS navigation measurement data reception is normal. If it is normal, go to step 4; if it is abnormal, determine that the GNSS observation set is an empty set, and go to step 9.

[0129] Step 4: Construct a high-dimensional GNSS observation vector for the multi-system, multi-frequency, and multi-type GNSS navigation measurement data obtained in step 3, and divide the high-dimensional observations into "carrier phase observation sub-set" and "pseudorange and pseudorange rate observation sub-set" according to statistical correlation. Determine whether the "carrier phase observation sub-set" is an empty set. If it is not empty, go to step 5; if it is an empty set, go to step 7.

[0130] Step 5: For the "carrier phase observation quantum set" obtained in step 4, establish the multi-system multi-frequency epoch inter-satellite differential carrier phase observation equation, and go to step 6.

[0131] Step 6: Based on the inertial navigation prediction and GNSS carrier phase observation at the current moment, an adaptive dynamic information filter is used to perform satellite / inertial tight integrated navigation information fusion and calculate the navigation error.

[0132] Step 7: Determine whether the "pseudorange and pseudorange rate observation quantum set" is an empty set. If it is an empty set, go to step 9; if it is not empty, establish the multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate observation equations, and go to step 8.

[0133] Step 8: Based on the current inertial navigation prediction and GNSS pseudorange and pseudorange rate observations, a robust adaptive static information filter is used to perform satellite / inertial tight integrated navigation information fusion and calculate the navigation error.

[0134] Step 9: If both subsets are empty sets, the navigation result is directly output; otherwise, the navigation result is corrected using the navigation error estimate and the final navigation result data is output.

[0135] It should be understood that although Figure 1 The steps in the flowchart are shown in sequence as indicated by the arrows, but these steps are not necessarily executed in the order indicated by the arrows. Unless otherwise specified in this document, there is no strict order restriction for the execution of these steps, and these steps can be executed in other orders. Moreover, Figure 1At least part of the steps may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily executed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential, but can be executed in turn or alternately with other steps or at least part of the sub-steps or stages of other steps.

[0136] In another specific embodiment, the airborne all-source satellite / inertial tight integrated navigation efficient dynamic and static information fusion method is further described:

[0137] Use an RF cable to connect the airborne integrated navigation receiver to the GNSS antenna, turn on the power, set the integrated navigation system GNSS information utilization mode to "dual system five frequencies (BDS B1\B2\B3+GPS L1\L2)", and set the arm vector from the airborne inertial navigation center to the GNSS antenna center.

[0138] Based on the inertial navigation solution and the system error state time update operation, the position, velocity, attitude and system error state transfer matrix of the carrier are calculated. , system noise transfer matrix , system error state covariance matrix and the system noise variance matrix .

[0139] Taking a small unmanned aerial vehicle airborne test as an example, select The inertial navigation and GNSS measurement data near the time are discretized at an interval of 8ms. First, the state transfer matrix in the system error state equation is calculated. and the system noise transfer matrix :

[0140] (37)

[0141] (38)

[0142] The information matrix and information error state vector after time update are given by the formula and After the time is updated, it is obtained from the inertial navigation forecast The carrier position at the moment is , the carrier speed is , the carrier posture is .

[0143] When receiving the time-synchronized satellite navigation data, the satellite navigation observation data and synchronous ephemeris validity check are first performed, including satellite number check (the number of available satellites for a single system is greater than 5), observation data frequency point completeness check, etc. At this time, the system enters the measurement update. At this time, the satellite navigation observation data time is At this moment, the visible satellites include 10 Beidou satellites and 9 GPS satellites. Beidou satellites have 3 frequency points (B1, B2 and B3), and GPS satellites have 2 frequency points (L1 and L2). The satellite signal cutoff elevation angle is set to 15°, and the number of available observation dimensions that meet the requirements of complete frequency points and elevation angle is 67, including 8 BDS B1 / B3 ionospheric-free combined pseudorange observations, 5 GPS L1 / L2 ionospheric-free combined pseudorange observations, 8 BDS B1 Doppler shift observations, 7 BDS B2 Doppler shift observations, 8 BDS B3 Doppler shift observations, 5 GPSL1 Doppler shift observations, 5 GPS L2 Doppler shift observations, 4 BDS B1 inter-epoch inter-satellite carrier phase difference observations, 3 BDS B2 inter-epoch inter-satellite carrier phase difference observations, 4 BDS B3 inter-epoch inter-satellite carrier phase difference observations, 5 GPS L1 inter-epoch inter-satellite carrier phase difference observations, and 5 GPS L2 inter-epoch inter-satellite carrier phase difference observations. All observation quantities constituting the 68-dimensional GNSS observation vector are divided into a 21-dimensional "carrier phase observation quantum set" and a 46-dimensional "pseudorange and pseudorange rate observation quantum set".

[0144] First, using the 21-dimensional “carrier phase observation quantum set”, the inter-epoch inter-satellite differential carrier phase observation equation for a single frequency point is established according to formulas (13) to (14), forming the observation vector described in formulas (15) to (14): and the observation design matrix Then, according to formula (20) to formula (23), the GNSS observation anomaly is detected and the GNSS observation anomaly value is downgraded. Finally, according to formula (24) to formula (25), the information matrix and the information error state vector Perform dynamic filtering measurement update, and the information error state vector after measurement update is:

[0145] (39)

[0146] Then, using the 46-dimensional “pseudorange and pseudorange rate observation quantum set”, the inter-satellite differential pseudorange and pseudorange rate observation equations for a single frequency point are established according to formulas (27) and (28), forming the observation vector shown in formula (35): and the observation design matrix .in accordance with and Complete the static equation forecast, detect GNSS observation anomalies channel by channel according to formula (21) to formula (23), and downgrade the GNSS observation anomalies. Finally, according to formula (24) to formula (25), the information matrix and the information error state vector Perform static filtering measurement update, and the information error state vector after measurement update is:

[0147] (40)

[0148] According to formula (26) to formula (27), the information error state vector Calculate the system error state for:

[0149] (41)

[0150] The system error state vector updated using dynamic and static measurements The carrier position, speed and attitude are corrected. After calculation, the carrier position after dynamic and static filtering measurement is updated: The carrier speed is , the carrier posture is On this basis, the inertial navigation measurement data can continue to be read and the time update process can be carried out forward.

[0151] In this embodiment, based on the performance analysis of the tightly integrated navigation system after the improved dynamic and static information filtering algorithm, by processing a set of trajectories such as Figure 3 , 80 minutes of airborne test data, to compare and verify the computational efficiency of the classic batch Kalman filter algorithm, the classic sequential Kalman filter algorithm, the information filter and the dynamic and static information filter algorithm used in the present invention. The four filtering algorithms are all developed in ANSI C and run on a laptop computer with Windows 11 system and Intel Core(TM) i7-1185G7 @ 3.00 GHz configuration; in the Debug mode of Visual Studio 2013 environment, the WIN32 functions QueryPerformanceFrequency() and QueryPerformanceCounter() are used to measure the computational time. By Figure 4 and Figure 5 It can be seen that the information filtering algorithm used in this method takes an average of about 3.5% of the time of the classical Kalman filtering algorithm, about 6% of the time of the classical sequential Kalman filtering algorithm, and 76% of the time of the classical information filtering algorithm during a single measurement update, which greatly reduces the computational burden of the embedded computing chip. The specific time statistics are shown in Table 1.

[0152] Table 1 Average calculation time for different combined navigation filter measurement updates in airborne tests

[0153]

[0154] In this embodiment, the inertial navigation data update rate is 125Hz, and the corresponding time interval is 8ms. In order to ensure the real-time performance of aircraft navigation, the Kalman filter time update and measurement update calculation must be completed within 8ms, and the measurement update calculation usually takes a long time. As can be seen from Table 1, among the four algorithms, the calculation time of the classic Kalman filter batch processing and the sequential Kalman filter algorithm is too long, while the dynamic and static information filtering algorithm of the present invention takes the shortest time, will not affect the normal Kalman time update step, and has the best real-time performance.

[0155] In order to verify the accuracy of this method, the real-time navigation results of an airborne test were analyzed, and the carrier position trajectory calculated by RTKLIB software was used as a benchmark to calculate the absolute position error. Based on the satellite / inertial combined navigation data post-processing software, three algorithms, namely classic batch Kalman filtering, sequential Kalman filtering and information filtering, were used to post-process the same set of airborne data, and the absolute position errors corresponding to different methods were calculated. Finally, the position error curves of the four methods changing with time were drawn. Figure 6 In the figure, the horizontal axis represents the epoch time in seconds, and the vertical axis represents the three-dimensional position error. The real-time position error output by the dynamic and static information filtering of the combined navigation system is represented by a solid line, and the real-time position error output by the classic batch Kalman filter, sequential Kalman filter and information filter is represented by a dotted line, a dashed line and a cross line, respectively. It can be seen that the four positioning error curves are basically consistent. In addition, due to misalignment at the initial moment, the positioning error is large, and the three-dimensional position error is basically within 2m after the late filtering converges. At the same time, it can be seen from Table 2 that the root mean square error and the maximum error of the positioning error of the four methods are basically the same. This shows that the dynamic and static information filtering algorithm of the present invention does not lose any positioning accuracy while reducing the amount of calculation, and can be applied to airborne real-time navigation scenarios.

[0156] Table 2 Comparison of positioning RMSE and maximum error of airborne experimental tight integration navigation

[0157]

[0158] The technical features of the above embodiments may be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0159] The above-described embodiments only express several implementation methods of the present application, and the descriptions thereof are relatively specific and detailed, but they cannot be understood as limiting the scope of the present application. It should be pointed out that, for a person of ordinary skill in the art, several variations and improvements can be made without departing from the concept of the present application, and these all belong to the protection scope of the present application. Therefore, the protection scope of the present application shall be subject to the attached claims.

Claims

1. An airborne all-source satellite-inertial tight integrated navigation efficient dynamic and static information fusion method, characterized in that: The method comprises: Build an airborne all-source satellite-inertial integrated navigation system. The integrated navigation system acquires the aircraft's inertial navigation measurement data and multi-system, multi-frequency, and multi-type GNSS navigation measurement data in real time during the carrier's motion. Performing inertial navigation solution according to the inertial navigation measurement data, establishing a system error state equation, and then performing Kalman time update; Searching for synchronized ephemeris for the GNSS navigation measurement data according to the GNSS observation data time, performing validity check of the observation data and synchronized ephemeris, and determining whether the GNSS navigation measurement data is received normally; If it is judged to be abnormal, the strapdown inertial navigation result data is directly output; If it is judged to be normal, a high-dimensional GNSS observation vector is constructed according to the GNSS navigation measurement data, and the high-dimensional GNSS observation vector is divided into a carrier phase observation quantum set and a pseudorange and pseudorange rate observation quantum set; a robust dynamic and static filtering algorithm based on successive fusion of independent observation information is constructed to perform information fusion of satellite / inertial tight integrated navigation, and calculate the navigation error, and use the navigation error estimate to correct the navigation result, and output the final navigation result data; the robust dynamic and static filtering algorithm based on successive fusion of independent observation information refers to performing robust dynamic Kalman filtering based on the dynamic model and the carrier phase observation quantum set at any epoch, and then sequentially adding the pseudorange and pseudorange rate observation quantum set to perform robust static Kalman filtering to obtain the global optimal fusion solution; Among them, a robust dynamic and static filtering algorithm based on the successive fusion of independent observation information is constructed to perform information fusion of satellite / inertial tight integrated navigation, calculate the navigation error, use the navigation error estimation value to correct the navigation result, and output the final navigation result data, including: Determining whether the carrier phase observation quantum set is empty; If the carrier phase observation quantum set is a non-empty set, a multi-system multi-frequency inter-epoch satellite differential carrier phase observation equation is established according to the carrier phase observation quantum set; according to the inertial navigation prediction amount and the GNSS carrier phase observation amount at the current moment, an anti-error adaptive dynamic information filter is used to perform satellite / inertial tightly integrated navigation information fusion, and the navigation error is calculated; it is determined whether the pseudorange and pseudorange rate observation quantum set is an empty set; if the pseudorange and pseudorange rate observation quantum set is an empty set, the navigation result is corrected by using the navigation error estimation value, and the final navigation result data is output; if the pseudorange and pseudorange rate observation quantum set is a non-empty set, a multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate observation equation is established; based on the inertial navigation prediction amount and the GNSS pseudorange and pseudorange rate observation amount at the current moment, an anti-error adaptive static information filter is used to perform satellite / inertial tightly integrated navigation information fusion, and the navigation error is calculated, the navigation result is corrected by using the navigation error estimation value, and the final navigation result data is output; If the carrier phase observation quantum set is an empty set, it is determined whether the pseudorange and pseudorange rate observation quantum set is an empty set; if the pseudorange and pseudorange rate observation quantum set is an empty set, the navigation result is directly output; if the pseudorange and pseudorange rate observation quantum set is a non-empty set, a multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate observation equation is established; based on the inertial navigation prediction quantity at the current moment and the GNSS pseudorange and pseudorange rate observation quantity, an anti-error adaptive static information filter is used to perform information fusion of satellite / inertial tight integrated navigation, and the navigation error is calculated; the navigation result is corrected using the navigation error estimate, and the final navigation result data is output.

2. The method according to claim 1, characterized in that According to the inertial navigation measurement data, an inertial navigation solution is performed, and a system error state equation is established, and then a Kalman time update is performed, including: According to the inertial navigation measurement data, the error state equation is constructed with position error, velocity error, misalignment angle, gyro bias and table bias as state quantities: ; in, represents the position error, represents the speed error, represents the misalignment angle, Indicates the gyro bias, Indicates the addition table zero bias, , , , and They represent the first-order derivatives of position error, velocity error, misalignment angle, gyro bias and table bias respectively. Represents the angular velocity of the Earth's rotation projected in the "North-East-Earth" navigation coordinate system; It represents the projection of the angular velocity of the n system relative to the earth coordinate system in the "north-east-earth" navigation coordinate system; ; represents the cross product matrix; Represents the ground velocity vector projected on the "North-East-Earth" navigation coordinate system; Represents the transformation matrix from the system to the navigation system; , and represents three intermediate error vectors; represents the apparent acceleration vector in the body coordinate system; and Denote the random noise of gyro and meter measurement respectively, which are described by a first-order Gauss-Markov process, where is unit Gaussian white noise, and represents the first-order Gauss-Markov model coefficients obtained by inertial navigation calibration, and the subscripts a and g represent the addition table and gyro, respectively; The error state vector is , the error state equation is converted into a discrete error state equation: ; in, express k The system error state vector at time ; express k System error state vector at time -1; Indicates from k -1 time to k The state transfer matrix at time; represents the system noise transfer matrix; express k -1 time instant white noise vector; Define the first auxiliary calculation matrix and the second auxiliary calculation matrix They are: ; ; in, represents the optimal estimate of the information matrix obtained at the last moment, represents the system process noise covariance matrix, ; According to the information error state vector, information matrix estimation value, first auxiliary calculation matrix and second auxiliary calculation matrix at the previous moment, the information error state vector and information matrix at the current moment are predicted to perform Kalman information filter time update, namely: ; ; in, Represents the information matrix forecast value obtained at the current moment, that is, ; represents the identity matrix, represents the information matrix, , represents the information error state vector at the previous moment, , Represents the information error state vector at the current moment.

3. The method according to claim 1, characterized in that Constructing a high-dimensional GNSS observation vector according to the GNSS navigation measurement data, and dividing the high-dimensional GNSS observation vector into a carrier phase observation sub-set and a pseudorange and pseudorange rate observation sub-set, including: Inter-satellite differential pseudo-range observations of BeiDou satellite n corresponding to B1 frequency point , intersatellite differential pseudorange rate observation and inter-epoch intersatellite differential carrier phase observations for: ; in, Indicates n The corresponding pseudo-range observation value of the B1 frequency point corresponding to the GNSS satellite, It represents the pseudo-range rate observation value corresponding to the B1 frequency point of the nth GNSS satellite. Indicates r The pseudo-range observation value of the B1 frequency point corresponding to the GNSS satellite No. Indicates r The pseudo-range rate observation of the B1 frequency point corresponding to the GNSS satellite, represents the time difference operator, , Respectively represent n The B1 frequency point corresponding to the GNSS satellite is and The GNSS carrier phase measurement value at time, and Indicates the reception time of two consecutive GNSS carrier phase measurement signals. , m Indicates the number of available satellites; , Respectively represent n No. and r The GNSS carrier phase time difference measurement value corresponding to the B1 frequency point of the GNSS satellite; The calculation formulas for the inter-satellite differential pseudorange, inter-satellite differential pseudorange change rate and inter-epoch inter-satellite differential carrier phase observations corresponding to the B1 frequency point of the nth Beidou satellite predicted by the inertial navigation are as follows: ; in, , and They represent the inter-satellite differential pseudorange, inter-satellite differential pseudorange variation rate and inter-epoch inter-satellite differential carrier phase observations corresponding to the B1 frequency point of the nth Beidou satellite predicted by the inertial navigation, and Respectively represent n The position and velocity vector of the GNSS satellite, and Indicates r The position and velocity vector of the reference satellite can be obtained by broadcast ephemeris calculation; and They represent the position and velocity vector at the phase center of the receiver GNSS antenna predicted by the inertial navigation respectively; Indicates that from n The unit direction vector from the satellite position to the receiver position, Indicates r The unit vector from the satellite position to the receiver position, and , calculated from the broadcast ephemeris and the inertial navigation predicted position; Indicates the inertial navigation speed; Inertial navigation forecast The direction cosine matrix from the b system to the n system at any moment; It represents the projection of the arm vector from the GNSS antenna center to the inertial navigation center in the b system; The Kalman filter high-dimensional GNSS observation vector constructed with the five-frequency observations of the Beidou and GPS dual systems is: in, represents the Kalman filtered high-dimensional GNSS observation vector; represents the inter-satellite differential pseudo-range observation column vector, subscript s Indicates the type of GNSS system used, s The value is BDS or GPS, which indicates the type of GNSS system used; and denote the inter-satellite differential pseudorange variation rate of each frequency point and the inter-epoch inter-satellite differential carrier phase observation column vector, respectively. c Indicates the type of GNSS system frequency used. c The value of is B1, B2, B3, L1 or L2.

4. The method according to claim 1, characterized in that A multi-system multi-frequency epoch inter-satellite differential carrier phase observation equation is established according to the carrier phase observation quantum set, including: The construction process of the inter-satellite differential carrier phase observation equation for a single-system single-frequency epoch includes: After obtaining the inter-satellite differential carrier phase observation corresponding to the nth satellite Based on this observation and the combined navigation error state The relationship between ; in, represents the inter-satellite differential carrier phase corresponding to the epoch of the nth satellite, Indicates the relative position of GNSS satellite n r The inter-satellite difference vector of the satellite-to-ground unit direction vector of the GNSS satellite No. n The unit vector pointing from the satellite position to the receiver position With r The unit vector pointing from the satellite position to the receiver position difference, , is the inter-epoch inter-satellite differential carrier phase measurement noise; and They are: ; in, δt Indicates the inertial navigation sampling time interval; Indicates time interval Number of internal inertial navigation measurement data; and denote the constant matrices used to extract velocity error and attitude error from the error state vector, respectively. represents the system dynamics model matrix, Indicates t k The SINS attitude matrix corresponding to the moment, Indicates t k-1 The SINS attitude matrix corresponding to the moment, Represents the projection of the arm vector from the user SINS center to the GNSS satellite antenna center in the carrier coordinate system; Represents the identity matrix Inter-satellite differential carrier phase observation vector between epochs of BeiDou B1 frequency point Combined navigation error status The relationship between them is: ; ; in, represents the inter-satellite differential carrier phase observation vector between epochs of BeiDou B1 frequency point, Indicates the integrated navigation error state, represents the noise vector of the inter-satellite differential carrier phase observation between epochs of the BeiDou B1 frequency point, , , Respectively expressed as r GNSS satellite No. is the reference satellite corresponding to No. 1, No. 2, m The inter-epoch differential carrier phase measurement noise of the GNSS satellites, , , Respectively expressed as r GNSS satellite No. is the reference satellite corresponding to No. 1, No. 2, m The satellite-to-ground unit direction vector of the GNSS satellite is the inter-satellite difference vector, represents the inter-satellite differential carrier phase observation design matrix between epochs of the B1 frequency point; Based on the construction process of the inter-satellite differential carrier phase observation equation between single-system single-frequency epochs, the inter-satellite differential carrier phase observation equation between multi-system multi-frequency epochs is constructed as follows: ; in, represents the multi-system multi-frequency epoch inter-satellite differential carrier phase observation vector, represents the noise vector of the inter-satellite differential carrier phase observation between multi-system and multi-frequency epochs, , They represent the inter-satellite differential carrier phase observation vectors of the BeiDou B2 and B3 frequency points respectively. , They represent the inter-satellite differential carrier phase observation vectors between epochs of GPS L1 and L2 frequency points, , They represent the inter-epoch inter-satellite differential carrier phase observation design matrices for the B2 and B3 frequency points, , Denote the inter-epoch inter-satellite differential carrier phase observation design matrices of L1 and L2 frequency points, , They represent the inter-satellite differential carrier phase observation noise vectors of BeiDou B2 and B3 frequency points respectively. , They represent the inter-satellite differential carrier phase observation noise vectors between epochs of GPS L1 and L2 frequency points, Design matrices for multi-system multi-frequency inter-epoch intersatellite differential carrier phase observations.

5. The method according to claim 1, characterized in that The process of robust adaptive dynamic information filtering includes: Use the current moment inertial navigation to predict the amount and GNSS carrier phase observations , and then detect GNSS observation anomalies channel by channel, and then use the variance matching method to perform weighted information fusion on the GNSS observation anomalies; The observation vector , observation matrix and the observation noise variance matrix By dividing into blocks, we can get: ; in, Indicates the current moment i GNSS channel observations; For the corresponding i The row vector of the observation design matrix corresponding to the channels; Indicates i The channel observation noise variance is i =1,2,……, m , m is the number of available observations at the current moment; represents the design matrix of multi-system multi-frequency epoch inter-satellite differential carrier phase observation; Current time i The GNSS observation anomaly test statistic corresponding to each channel is: ; in, Indicates i GNSS channel updates, , Represents the row vector of the observation design matrix corresponding to the i-th channel at the current moment, Indicates the current SINS forecast error state value; Indicates i GNSS channel innovation variance, , represents the SINS forecast error state covariance matrix at the current moment, For the corresponding i The column vector of the observation design matrix corresponding to the channels; like , then judge the i The GNSS channels are observed normally, and the settings ;in Represents the standard normal distribution Quantile; like , then judge the i GNSS channel observation anomalies, set the diagonal elements of the weight matrix for: ; in, , are two preset constants, i =1,2,……, m ; Traverse all m After observing channels, construct the weight matrix ; A fusion formula is used to perform information filtering and information fusion to calculate the navigation error; the fusion formula is: ; ; ; ; in, , , and Respectively represent the current time after the dynamic information filter measurement is updated k The information error state vector, information matrix, error state covariance matrix and error state vector, Indicates the current time k The information error state vector state prediction information, Represents the information matrix forecast value obtained at the current moment.

6. The method according to claim 1, characterized in that The steps of establishing the multi-system multi-frequency inter-satellite differential pseudorange and pseudorange variation rate observation equations include: After obtaining the inter-satellite differential B1 / B3 dual-frequency ionospheric elimination combined pseudo-range observation corresponding to the nth satellite Based on the pseudo-range observation Combined navigation error status The relationship between them is: ; in, is the inter-satellite differential pseudo-range measurement noise, It represents the Beidou B1 / B3 dual-frequency ionospheric-free combined pseudorange observation. It represents the inter-satellite difference vector of the satellite-to-ground unit direction vector of BDS satellite n relative to BDS satellite r, represents the direction cosine matrix from the local navigation coordinate system to the Earth-centered Earth-fixed system, represents the attitude matrix output by SINS, It represents the projection of the arm vector from the GNSS antenna center to the inertial navigation center in the b system, represents the BDS dual-frequency linear combination pseudorange measurement noise, represents the SINS attitude misalignment angle vector, represents the SINS navigation position error vector; The single-system single-frequency inter-satellite differential pseudo-range rate observation equation is constructed at the BeiDou B1 frequency point; n The inter-satellite differential pseudo-range rate observation corresponding to the satellite Combined navigation error status The relationship between them is: ; in, is the inter-satellite differential pseudorange rate measurement noise; It represents the measured value of the inter-satellite differential pseudo-range rate of the B1 frequency point corresponding to the nth Beidou satellite with the rth Beidou satellite as the reference satellite. Represents the projection vector of the angular velocity vector of the local navigation coordinate system relative to the inertial coordinate system in the local navigation coordinate system, Represents the projection vector of the angular velocity vector of the carrier coordinate system measured by the gyro relative to the inertial coordinate system in the body coordinate system, represents the SINS velocity error vector in the local navigation coordinate system, represents the SINS gyro bias error vector, represents the inter-satellite differential pseudo-range rate measurement noise vector of the B1 frequency point; BeiDou B1 / B3 dual-frequency intersatellite differential pseudo-range observation vector , B1 frequency point inter-satellite differential pseudo-range rate observation vector Combined navigation error status The relationships between them are: ; ; in, represents the BDS dual-frequency linear combination pseudo-range measurement value vector, represents the BDS dual-frequency linear combination pseudorange measurement noise vector, , , Respectively represent No. 1, No. 2, m The inter-satellite difference vector of the satellite-to-ground unit direction vector of BDS satellite No. 1 relative to BDS satellite No. 2, represents the direction cosine matrix from the local navigation coordinate system to the Earth-centered Earth-fixed system, represents the attitude matrix output by SINS, It represents the projection of the arm vector from the GNSS antenna center to the inertial navigation center in the b system, , , Respectively represent No. 1, No. 2, m BDS satellite relative to r The BDS dual-frequency linear combination pseudo-range measurement noise of BDS satellite No. represents the inter-satellite differential pseudorange rate observation matrix of the B1 frequency point, represents the inter-satellite differential pseudo-range rate observation noise at the B1 frequency point, , , Respectively represent No. 1, No. 2, m BDS satellite relative to r The noise vector of the inter-satellite differential pseudo-range rate observation of BDS satellite No. 1 at the B1 frequency point; Based on the steps of constructing the single-system single-frequency inter-satellite differential pseudorange rate observation equation at the Beidou B1 frequency point, the multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate observation equations are established as follows: ; in, represents the multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate observation vector, represents the multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate observation design matrix, represents the multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate observation error vector, represents the BeiDou B1 / B3 ionospheric-free combined pseudorange observation matrix, , and They represent the inter-satellite differential pseudorange rate observation matrices of B1, B2 and B3 frequency points respectively, represents the GPS L1 / L2 ionospheric-free combined pseudorange observation matrix, and They represent the inter-satellite differential pseudorange rate observation matrices of L1 and L2 frequency points, represents the BeiDou B1 / B3 ionospheric-free combined pseudorange observation error, , and represent the inter-satellite differential pseudorange rate observation errors of B1, B2 and B3 frequency points respectively, represents the GPS L1 / L2 ionospheric-free combined pseudorange observation error, and They represent the inter-satellite differential pseudo-range rate observation errors of L1 and L2 frequency points, represents the BeiDou B1 / B3 dual-frequency inter-satellite differential pseudo-range observation vector, represents the GPS dual-frequency inter-satellite differential pseudo-range observation vector, , , , , They represent the inter-satellite differential pseudorange rate observation vectors of B1, B2, B3, L1 and L2 frequency points respectively.

7. The method according to claim 1, characterized in that The process of robust adaptive static information filtering includes: Detect GNSS observation anomalies channel by channel, and then use the variance matching method to reduce the weight of GNSS observation anomalies. Use the state vector obtained by the first step of dynamic filtering as the state prediction information for static filtering. Using the observations, observation matrix and noise matrix related to GNSS pseudorange and pseudorange rate observations, we set , construct the weight matrix ;in, represents the observation vector, represents the observation matrix, represents the multi-system multi-frequency inter-satellite differential pseudorange and pseudorange rate observation vector, represents the design matrix of multi-system multi-frequency inter-satellite differential pseudorange and pseudorange variation rate observation; The static filtering calculation uses a fusion formula to fuse information and calculate the navigation error; the fusion formula is: ; ; ; ; in, , , and Respectively represent the current time after the dynamic information filter measurement is updated k The information error state vector, information matrix, error state covariance matrix and error state vector, represents the observation noise variance matrix, Indicates the current time k The information error state vector state prediction information, Represents the information matrix forecast value obtained at the current moment.

Images