GPS / BDS dual-mode satellite positioning method based on LS and IBAS-AUKF combined calculation

Abstract

The invention discloses a GPS / BDS dual-mode satellite positioning method based on LS and IBAS-AUKF combined solution, which comprises the following steps of: respectively establishing a GPS satellite positioning system model and a BDS satellite positioning system model, realizing pseudo-range measurement through intermediate frequency signals of a GPS satellite and a BDS satellite received by a receiver, obtaining a combined positioning model based on the pseudo-range measurement, and solving by utilizing a least square method model to obtain a positioning result; obtaining a first-order least square estimation value of the position solution; initializing a Kalman filter by using the primary least square estimation value, and selecting an improved unscented Kalman filter as a dual-mode navigation filter to carry out state estimation to obtain a state estimation value; an improved beetle antennae search algorithm is introduced to update a process noise and observation noise covariance matrix; and obtaining relative position information of the receiver according to the state estimation value. The method provided by the invention can improve the positioning precision of the system in a dynamic environment, improve the convergence speed, solve the problem of uncertainty of noise covariance, and have good anti-interference capability.

Description

Technical Field

[0001] This invention relates to the field of satellite positioning technology, and in particular to a GPS / BDS dual-mode satellite positioning method based on joint calculation of LS and IBAS-AUKF. Background Technology

[0002] With the rapid development of navigation and positioning systems, multi-system GNSS satellite joint positioning has gradually become one of the main research trends. Compared with a single satellite navigation and positioning system, multi-constellation joint positioning will significantly increase the number of visible satellites, greatly improve the satellite ground space geometry, thereby weakening the ill-conditioned nature of the coefficient matrix of the normal equation of a single system, and effectively improving the accuracy, continuity and reliability of satellite positioning.

[0003] Currently, the BDS / GPS combined operating mode is the commonly used application mode. However, existing dual-mode software receiver positioning algorithms have low positioning accuracy and large fluctuations in dynamic environments. Therefore, there is an urgent need for a positioning method with higher accuracy in dual-mode systems. Summary of the Invention

[0004] This invention provides a GPS / BDS dual-mode satellite positioning method based on joint calculation of LS and IBAS-AUKF, which aims to effectively solve the problem of low positioning accuracy and large fluctuation in the existing dual-mode software receiver positioning algorithm in dynamic environment, and overcome the shortcomings of using the least squares method and Kalman filter method alone in single-point positioning solution.

[0005] This invention provides a GPS / BDS dual-mode satellite positioning method based on joint calculation of LS and IBAS-AUKF, the method comprising:

[0006] The digital intermediate frequency (IF) signals of BDS and GPS satellite signals are acquired, and the digital IF signals are mediated to obtain navigation information and pseudorange measurement values.

[0007] GPS satellite positioning system model and BDS satellite positioning system model are established respectively. Based on the navigation information and pseudorange measurement values, the pseudorange equation of GPS satellite and pseudorange equation of BDS satellite are solved by pseudorange equation solving method.

[0008] A combined positioning model is obtained based on the pseudorange equations of the GPS satellites and the BDS satellites.

[0009] The combined positioning model is solved using the least squares method to obtain the first least squares estimate of the position solution;

[0010] The Kalman filter is initialized using the first least squares estimate, and the improved unscented Kalman filter is selected as the dual-mode navigation filter for state estimation to obtain the state estimate.

[0011] The relative position information of the receiver is obtained based on the state estimate.

[0012] Furthermore, the combined positioning model is as follows:

[0013]

[0014] Where (AX=b) is the least squares form of the combined localization model, [X i0 ,Y i0 Z i0 [This represents the approximate location of the receiver.] and Let ΔX represent the geometric distances from the GPS satellite and the BDS satellite to the receiver, respectively. i ,ΔY i ,ΔZ i It is the difference between the approximate and actual coordinates of the receiver, where c is the speed of light, and dt is the distance between the approximate and actual positions of the receiver. Gi and dt Ci These are the receiver clock errors for GPS satellites and BDS satellites, respectively, where G is the GPS satellite, C is the BDS satellite, and b = P. i -ρ i +c·dt G +c·dt C -T i -I i -ε i I i T i ε i These represent ionospheric delay, tropospheric delay, and other delay error terms, dt G It is the star clock error, P i ρ is the pseudorange measurement from the satellite to receiver i. i It is the geometric distance from the receiver to the satellite, [x Gm ,y Gm ,z Gm ] and [x Cn ,y Cn ,z Cn These are the GPS satellite positions and BDS satellite positions, respectively.

[0015] Further, the step of solving the combined positioning model using the least squares method to obtain a first least squares estimate of the location solution includes:

[0016] The GPS / BDS pseudorange linear equation is calculated according to the preset least squares calculation formula. The equation is iterated repeatedly to minimize b, and the position correction value that minimizes b is used as the first least squares estimate of the position solution.

[0017] dx0=(A T A) -1 A T b;

[0018] x1=x k-1 +dx0;

[0019] In the formula, dx0 is the correction amount. The initial value is corrected using the correction amount dx0 to obtain the corrected value x1. k-1 Let A represent the state estimate at time k-1. A is a known matrix representing the coefficients in the localization model.

[0020] Furthermore, the step of selecting an improved unscented Kalman filter as a dual-mode navigation filter for state estimation to obtain state estimates includes:

[0021] The initial position value calculated by least squares is assigned to the Kalman filter initialization, and the state estimate is obtained according to the following state equation and measurement equation:

[0022] Take error state Establish discrete state equations:

[0023]

[0024] Where x is the actual state quantity at the location, Let F be the initial position value after initialization, and w be the state transition matrix. k The system excitation noise sequence has statistical properties derived from the process noise covariance matrix Q. k describe;

[0025] The measurement equation is a linearized GPS / BDS pseudorange equation:

[0026] δρ=H·δx(k)+v(k);

[0027] Here, the row vector of H corresponds to the partial derivative of the observation equation of each satellite at the initial position, and v(k) is the measurement noise sequence.

[0028] Furthermore, the step of selecting an improved unscented Kalman filter as a dual-mode navigation filter for state estimation to obtain state estimates further includes:

[0029] The system state at time k is updated using an improved unscented Kalman filter algorithm, and the state estimate is obtained based on the state update result.

[0030] Furthermore, the method of updating the system state at time k using the improved unscented Kalman filter algorithm includes:

[0031] S1 is obtained by calculating the Sigma point and its weight using the following formula:

[0032]

[0033]

[0034] In the formula, Let λ represent the i-th column of the square root of the matrix, where λ = α 2 (n+κ)-n, where λ is the scaling parameter, the value of α depends on the distribution of the sampling points, κ is a parameter, n is the dimension of the state vector, and X0 represents the initial state value. X represents the mean of a Gaussian distribution at time k-1. i (k-1) represents the weighting coefficients corresponding to the Sigma point, P k Represents the covariance matrix. β is 2, where, This represents the weights required to calculate the covariance. This represents the weights needed to calculate the mean. This represents the weights required to calculate the mean at the initial time step;

[0035] S2, Calculate the one-step predicted state vector X (i) (k+1|k) and one-step predicted state covariance

[0036] X (i) (k+1|k)=f(X (i) (k|k),u k-1 );

[0037]

[0038]

[0039] Among them, X (i) (k+1|k) represents the sample points for one-step prediction, X (i) (k|k) represents the i-th sampling point at time k. It is a one-step prediction of the estimated state at time k, u k-1 For deterministic control terms, Q k It is the process noise covariance. This represents the state prediction covariance at time k;

[0040] The measurement update formula is as follows:

[0041]

[0042]

[0043]

[0044]

[0045]

[0046]

[0047] Among them, y (i) (k+1|k) represents the measurement prediction sample points, P yy,k and P xy,k These are the prediction measurement covariance matrix and the cross covariance matrix, respectively. To estimate the state, P(k+1|k+1) is the estimated covariance matrix, and K... k+1 R is the gain matrix; k To observe the noise covariance, Let y(k+1) be the one-step predicted mean of the measured variable, y(k+1) be the measured value at time k+1, and P(k+1|k) be the estimated covariance matrix.

[0048] S3, repeat S1-S2 until the preset iteration termination condition is reached, and obtain the LS-UKF joint estimated receiver position as the receiver relative position information.

[0049] Furthermore, the method of using the improved unscented Kalman filter algorithm to update the system state at time k also includes:

[0050] The process noise covariance and the observation noise covariance are updated using the beetle whisker search algorithm to obtain the updated process noise covariance and observation noise covariance.

[0051] Furthermore, the update of the process noise covariance and the observation noise covariance using the beetle whisker search algorithm includes:

[0052] Vector normalization was performed on the head orientation of the longhorn beetle:

[0053]

[0054] Where t represents the number of unknowns in the problem to be solved, and r is the result of vector normalization;

[0055] Obtain the spatial coordinates of the longhorn beetle's left and right tentacles:

[0056]

[0057] In the formula, x t It is the location of the centroid of the longhorn beetle at the t-th iteration. These are the coordinates of the left tentacle's position during the t-th iteration. d is the position coordinate of the right tentacle during the t-th iteration. t It is the distance between the two tentacles of the longhorn beetle;

[0058] Construct a fitness function f(·), using the absolute cumulative error value of the innovation as the fitness function f(·) of the BAS algorithm, to obtain the optimal process noise covariance Q. k and observation noise covariance R k :

[0059]

[0060] In the formula, the row vector of H corresponds to the partial derivative of the observation equation of each satellite at the initial position;

[0061] The fitness values ​​of the locations of the left and right antennae of the longhorn beetle are calculated based on the fitness function f(·). and And update the location:

[0062]

[0063] In the formula, δ t Let x represent the search step size in the t-th iteration. t Let x be the location of the longhorn beetle at time t. t+1 Let t+1 be the location of the longhorn beetle;

[0064] The distance d between the two tentacles of the longhorn beetle t and the search step size δ of the t-th iteration t Update:

[0065]

[0066] In the formula, d e It is the attenuation factor of the distance between the two whiskers, and α represents the step size reduction factor.

[0067] Furthermore, the fitness values ​​of the locations of the left and right antennae of the longhorn beetle are calculated based on the fitness function f(·). and After updating the location, the method also includes:

[0068] The simulated annealing algorithm is used to determine whether the updated position meets the preset loop termination condition, including:

[0069] The probability P corresponding to the fitness value of the location of the longhorn beetle is calculated according to the preset probability formula:

[0070]

[0071] In the formula, f(x) t+1 f(x) t ) are the fitness functions at different times, where T is the current temperature, and when f(x) t+1 )<f(x t When x is in the specified position, it indicates that x is in the specified position. t+1 This is the current optimal solution; we can directly accept x. t+1 For the new solution, when f(x) t+1 )≥f(x t When x is given, accept the inferior solution x with probability P. t+1 This is a new solution, which is a new process noise covariance and observation noise covariance.

[0072] Furthermore, the step size reduction factor α is updated in the following manner to obtain the updated step size reduction factor β:

[0073]

[0074] In the formula, f t T is the fitness function in the BAS algorithm at the current iteration number t. max This represents the total number of iterations in the BAS algorithm, α is the original step size reduction factor, and f min It is the minimum fitness function value.

[0075] Through one or more embodiments of the above-described embodiments of the present invention, at least the following technical effects can be achieved: Combining LS (Least Squares) with UKF overcomes the problems of LS obtaining user pseudorange information, where the information at each moment is independent and cannot be jointly utilized, and the positioning error is relatively large in the presence of noise interference. Meanwhile, UKF's positioning results are relatively smooth, and initializing UKF with LS avoids positioning divergence caused by unreasonable initial value settings. Considering the more complex structure of the dual-mode navigation system and the nonlinear nature of the pseudorange observation equation, the unscented Kalman filter method is used as the solution model to improve the positioning accuracy of the system in dynamic environments.

[0076] In addition, based on the LS-UKF joint localization algorithm, the Improved Beetle Antennae Search Algorithm (IBAS) is introduced to update the noise and observation noise covariance. Simulated annealing algorithm and adaptive factor are added to the BAS algorithm to avoid the BAS algorithm getting trapped in local optima, while improving the convergence speed, solving the noise covariance uncertainty problem, and having good anti-interference ability. Attached Figure Description

[0077] The technical solution and other beneficial effects of the present invention will become apparent from the following detailed description of specific embodiments of the invention, in conjunction with the accompanying drawings.

[0078] Figure 1 This is one of the flowcharts for the GPS / BDS dual-mode satellite positioning method based on joint calculation of LS and IBAS-AUKF provided in the embodiments of the present invention;

[0079] Figure 2 This is the second flowchart of the GPS / BDS dual-mode satellite positioning method based on the joint calculation of LS and IBAS-AUKF provided in the embodiments of the present invention. Detailed Implementation

[0080] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0081] In the description of this invention, it should be noted that, unless otherwise explicitly specified and limited, the term "and / or" in this document is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, and B existing alone. Furthermore, the character " / " in this document, unless otherwise specified, generally indicates that the preceding and following related objects have an "or" relationship.

[0082] See Figure 1 and Figure 2 This application provides a GPS / BDS dual-mode satellite positioning method based on joint calculation of LS and IBAS-AUKF, which includes the following steps:

[0083] S101, acquire the digital intermediate frequency signals of BDS and GPS satellite signals, perform signal modulation on the digital intermediate frequency signals, and obtain navigation information and pseudorange measurement values.

[0084] In this step, the BDS / GPS dual-mode receiver RF front-end processing module receives all visible satellite signals through the antenna, down-converts them to intermediate frequency (IF) signals via mixing, and finally converts the IF signals into discrete-time BDS and GPS digital IF signals via an analog-to-digital (A / D) converter. IF (n), and respectively the digital intermediate frequency signals S of BDS and GPS IF (n) Perform signal conditioning to obtain their respective navigation information and pseudorange measurements. The navigation information includes clock errors, satellite position and velocity, ionospheric and tropospheric delays, etc.

[0085] S102, establish GPS satellite positioning system model and BDS satellite positioning system model respectively, and use pseudorange equation solving method to obtain GPS satellite pseudorange equation and BDS satellite pseudorange equation based on the navigation information and pseudorange measurement value.

[0086] Specifically, the pseudorange equations P for GPS satellites are obtained according to the following pseudorange formula. G Pseudorange equation P of BDS satellite C :

[0087] P=ρ+c(dt u,k -dt k (s) )+I k +T k +ε ρ,k ;

[0088]

[0089] Where ρ is the observation pseudorange, the receiver coordinates are (x, y, z), and the satellite coordinates are (x, y, z). n ,y n ,z n ), rk is the geometric distance between the satellite and the receiver, δt u,k and These are the receiver clock bias and the satellite clock bias, respectively. k and T k These are ionospheric delay and tropospheric delay, respectively. c For the speed of light, ε ρ, k represents the measurement error of the pseudorange.

[0090] GPS satellite pseudorange equation P G =ρ G +c(dt G -dt G )+I G +T G +ε ρ,G The pseudorange equation P of the BDS satellite C =ρ C +c(dt C -dt C )+I C +T C +ε ρ,C In the formula, P G The pseudorange observation equation from GPS satellite to receiver, ρ G dt is the geometric distance from the receiver to the GPS satellite, c is the speed of light, and dt is the distance between the receiver and the GPS satellite. G It is the receiver clock error of the GPS satellites, dt GIt is the clock error of GPS satellites, I G T G ε ρ,G It is the ionospheric delay, tropospheric delay, and other delay error terms of GPS satellites, P C It is the pseudorange observation equation from the BDS satellite to the receiver, ρ C It is the geometric distance dt from the receiver to the BDS satellite. C This is the receiver clock error of the BDS satellite, dt C It is the clock error of the BDS satellite, I C T C ε ρ,C It includes the ionospheric delay, tropospheric delay, and other delay error terms of the BDS satellite.

[0091] S103, A combined positioning model is obtained based on the pseudorange equations of the GPS satellites and the BDS satellites.

[0092] It should be noted that since the receiver clock has different deviations relative to BeiDou (i.e., BDS satellites) and GPS, the receiver clock error must be calculated separately during the combined positioning process. Therefore, at least 5 satellites are required to achieve combined positioning in this step.

[0093] S104, Solve the combined positioning model using the least squares method to obtain a first least squares estimate of the position solution. Specifically, use the least squares method to calculate the GPS / BDS pseudorange linear equation, iterate repeatedly to minimize b, calculate the relative coordinate correction value, and then obtain a first least squares estimate of the position solution.

[0094] S105, the Kalman filter is initialized using the first least squares estimate, and an improved unscented Kalman filter is selected as a dual-mode navigation filter for state estimation to obtain the state estimate.

[0095] S106, Obtain the relative position information of the receiver based on the state estimate.

[0096] In some embodiments of this application, the combined positioning model is:

[0097]

[0098] Where (AX=b) is the least squares form of the combined localization model, [X i0 ,Y i0 Z i0 [This represents the approximate location of the receiver.] and Let ΔX represent the geometric distances from the GPS satellite and the BDS satellite to the receiver, respectively. i ,ΔY i,ΔZ i It is the difference between the approximate and actual coordinates of the receiver, where c is the speed of light, and dt is the distance between the approximate and actual positions of the receiver. Gi and dt Ci These are the receiver clock errors for GPS satellites and BDS satellites, respectively, where G is the GPS satellite, C is the BDS satellite, and b = P. i -ρ i +c·dt G +c·dt C -T i -I i -ε i I i T i ε i These represent ionospheric delay, tropospheric delay, and other delay error terms, dt G It is the star clock error, P i ρ is the pseudorange measurement from the satellite to receiver i. i It is the geometric distance from the receiver to the satellite, [x Gm ,y Gm ,z Gm ] and [x Cn ,y Cn ,z Cn [x] represents the GPS satellite position and the BDS satellite position, respectively. Gm ,y Gm ,z Gm ] and [x Cn ,y Cn ,z Cn The pseudorange observation ρ can be calculated from the ephemeris data at the time of signal transmission. Gm ρ Bn As known quantities, they can be obtained from navigation messages.

[0099] In some embodiments of this application, the step of solving the combined positioning model using the least squares method to obtain a first least squares estimate of the location solution includes:

[0100] Based on the preset least squares calculation of the GPS / BDS pseudorange linear equation, the process iterates repeatedly to minimize b, and uses the position correction value that minimizes b as the first least squares estimate of the position solution:

[0101] dx0=(A T A) -1 A T b;

[0102] x1=x k-1 +dx0;

[0103] In the formula, dx0 is the correction amount. The initial value is corrected using the correction amount dx0 to obtain the corrected value x1.k-1 Let represent the state estimate at time k-1. Finally, the position correction value that minimizes b is taken as the first least squares estimate of the position solution. A is a known matrix representing the coefficients in the positioning model.

[0104] In some embodiments of this application, the step of selecting an improved unscented Kalman filter as a dual-mode navigation filter for state estimation to obtain state estimates includes:

[0105] The initial position value calculated by least squares is assigned to the Kalman filter initialization, and the state estimate is obtained according to the following state equation and measurement equation:

[0106] Take error state Establish discrete state equations:

[0107]

[0108] Where x is the actual state quantity at the location, Let F be the initial position value after initialization, and w be the state transition matrix. k The system excitation noise sequence has statistical properties derived from the process noise covariance matrix Q. k describe;

[0109] The measurement equation is a linearized GPS / BDS pseudorange equation:

[0110] δρ=H·δx(k)+v(k);

[0111] Here, the row vector of H corresponds to the partial derivative of the observation equation of each satellite at the initial position, and v(k) is the measurement noise sequence.

[0112] In some embodiments of this application, the step of selecting an improved unscented Kalman filter as a dual-mode navigation filter for state estimation to obtain a state estimate further includes:

[0113] The system state at time k is updated using an improved unscented Kalman filter algorithm, and the state estimate is obtained based on the state update result.

[0114] In some embodiments of this application, the method of using an improved unscented Kalman filter algorithm to update the system state at time k includes:

[0115] S1 is obtained by calculating the Sigma point and its weight using the following formula:

[0116]

[0117]

[0118] In the formula, Let λ represent the i-th column of the square root of the matrix, where λ = α 2 (n+k)-n, where λ is the scaling parameter, the value of α is related to the distribution of the sampling points, κ is often 0 or 3-n, and is selected according to the actual situation, n is the dimension of the state vector, and X0 represents the initial state value. X represents the mean of a Gaussian distribution at time k-1. i (k-1) represents the weighting coefficients corresponding to the Sigma point, P k Let represent the covariance matrix. β is 2, where, This represents the weights required to calculate the covariance. This represents the weights needed to calculate the mean. This represents the weights required to calculate the mean at the initial time.

[0119] S2, Calculate the one-step predicted state vector X (i) (k+1|k) and one-step predicted state covariance

[0120] X (i) (k+1|k)=f(X (i) (k|k),u k-1 );

[0121]

[0122]

[0123] Among them, X (i) (k+1|k) represents the sample points for one-step prediction, X (i) (k|k) represents the i-th sampling point at time k. It is a one-step prediction of the estimated state at time k, u k-1 For deterministic control terms, Q k It is the process noise covariance. This represents the state prediction covariance at time k.

[0124] The measurement update formula is as follows:

[0125]

[0126]

[0127]

[0128]

[0129]

[0130]

[0131] Among them, y (i) (k+1|k) represents the measurement prediction sample points, P yy,k and P xy,k These are the prediction measurement covariance matrix and the cross covariance matrix, respectively. To estimate the state, P(k+1|k+1) is the estimated covariance matrix, and K... k+1 R is the gain matrix; k To observe the noise covariance, Let y(k+1) be the one-step predicted mean of the measured variable, y(k+1) be the measured value at time k+1, and P(k+1|k) be the estimated covariance matrix.

[0132] S3, repeat S1-S2 until the preset iteration termination condition is reached, and obtain the LS-UKF joint estimated receiver position, that is, the receiver relative position information.

[0133] In some embodiments of this application, considering that process noise and observation noise covariance are two important parameters in the UKF algorithm positioning solution, and are generally given a fixed value based on experience, this leads to a certain error in position estimation. From the perspective of intelligent search algorithms, the Beetle Antennae Search Algorithm (BAS) is introduced to optimize the two noise covariances, thereby updating the noise covariance. The Beetle Antennae Search Algorithm is a mathematical model abstracted from the living habits of longhorn beetles and applied to various problems requiring optimal solutions.

[0134] In this example, the process noise covariance and the observation noise covariance are updated using the beetle whisker search algorithm to obtain the updated process noise covariance and observation noise covariance.

[0135] In some embodiments of this application, updating the process noise covariance and the observation noise covariance using the beetle whisker search algorithm includes:

[0136] Vector normalization was performed on the head orientation of the longhorn beetle:

[0137]

[0138] Where t represents the number of unknowns in the problem to be solved, and r is the result of vector normalization.

[0139] Specifically, since the direction of the longhorn beetle's head is uncertain before each move, it can be randomly generated using the random function rand(·), with the positive direction of the vector defined as from the right whisker to the left whisker.

[0140] Obtain the spatial coordinates of the longhorn beetle's left and right tentacles:

[0141]

[0142] In the formula, x t It is the location of the centroid of the longhorn beetle at the t-th iteration. These are the coordinates of the left tentacle's position during the t-th iteration. d is the position coordinate of the right tentacle during the t-th iteration. t It is the distance between the two tentacles of the longhorn beetle.

[0143] Construct a fitness function f(·), using the absolute cumulative error value of the innovation as the fitness function f(·) of the BAS algorithm, to obtain the optimal process noise covariance Q. k and observation noise covariance R k :

[0144]

[0145] In the formula, the row vector of H corresponds to the partial derivative of the observation equation of each satellite at the initial position.

[0146] The fitness values ​​of the locations of the left and right antennae of the longhorn beetle are calculated based on the fitness function f(·). and And update the location:

[0147]

[0148] In the formula, δ t Let x represent the search step size in the t-th iteration. t Let x be the location of the longhorn beetle at time t. t+1 Let t+1 be the location of the longhorn beetle.

[0149] The fitness value mentioned above represents the odor intensity in the search area. In this step, the position is updated by comparing the odor intensity at the left and right tentacle positions to determine whether to move to the left or right.

[0150] The distance d between the two tentacles of the longhorn beetle t and the search step size δ of the t-th iteration t Update:

[0151]

[0152] In the formula, d e It is the attenuation factor of the distance between the two whiskers, and α represents the step size reduction factor, which is generally set to 0.95.

[0153] It should be noted that this application uses the IBAS algorithm to update the process noise and observation noise covariance. The initial value of the noise covariance is input into the IBAS algorithm, and after automatic optimization, the optimal noise covariance at the current time is obtained. Then, the optimal noise covariance is passed to UKF to update the noise covariance, thus completing the estimation of the position at the current time.

[0154] In some embodiments of this application, since the BAS algorithm (i.e., the beetle whisker search algorithm) uses only one beetle, it is prone to getting trapped in local optima. Therefore, simulated annealing is introduced. A core principle of simulated annealing is the Metropolis criterion, which causes the algorithm to accept a poor solution with a certain probability. In the early stages of the algorithm, this probability is relatively high, so the algorithm is very likely to accept a poor solution and will not get trapped in local optima. As the temperature decreases, the algorithm gradually converges to the optimal solution, and this probability gradually approaches zero, at which point the algorithm accepts a poor solution again.

[0155] Specifically, the fitness values ​​of the locations of the left and right antennae of the longhorn beetle are calculated based on the fitness function f(·). and After updating the location, the method also includes:

[0156] The simulated annealing algorithm is used to determine whether the updated position meets the preset loop termination condition.

[0157] In some embodiments of this application, determining whether the updated position satisfies a preset loop termination condition using the simulated annealing algorithm includes:

[0158] The probability P corresponding to the fitness value of the location of the longhorn beetle is calculated according to the preset probability formula:

[0159]

[0160] In the formula, f(x) t+1 f(x) t ) are the fitness functions at different times, where T is the current temperature, and when f(x) t+1 )<f(x t When x is in the specified position, it indicates that x is in the specified position. t+1 This is the current optimal solution; we can directly accept x. t+1 For the new solution, when f(x) t+1 )≥f(x t When x is given, accept the inferior solution x with probability P. t+1 This is a new solution, which is a new process noise covariance and observation noise covariance.

[0161] Furthermore, in the BAS algorithm, the step size decreases linearly, which prolongs the convergence time. An adaptive factor β is introduced to further adjust the step size reduction factor α, thereby shortening the convergence time. The expression is as follows:

[0162]

[0163] In the formula, f t This is the fitness function in the BAS algorithm at the current iteration number t, where t is the current iteration number in the BAS algorithm, and T is the fitness function. max f is the total number of iterations of the BAS algorithm. min It is the minimum fitness function value.

[0164] Based on the above adaptive factor β, probability P, and the updated distance d between the two tentacles of the longhorn beetle. t+1 and the updated search step size δ for the (t+1)th iteration t+1 Updated position x t+1 Obtain the two noise covariances Q after adaptive update k and R k And adaptively update the two noise covariances Q k and R k The information is imported into a UKF filter to calculate the relative position information of the receiver.

[0165] The GPS / BDS dual-mode satellite positioning method based on joint calculation of LS and IBAS-AUKF provided in this application combines LS (least squares) with UKF. This overcomes the problems of LS, which calculates current information based on user pseudorange independently at each moment, making joint utilization impossible and resulting in relatively large positioning errors in the presence of noise interference. Meanwhile, UKF provides smoother positioning results, and LS is used to initialize UKF, avoiding positioning divergence caused by unreasonable initial value settings. Considering the more complex structure of dual-mode navigation systems and the nonlinear nature of pseudorange observation equations, an unscented Kalman filter is used as the solution model to improve the positioning accuracy of the system in dynamic environments.

[0166] In addition, based on the LS-UKF joint localization algorithm, the Improved Beetle Antennae Search Algorithm (IBAS) is introduced to update the noise and observation noise covariance. Simulated annealing algorithm and adaptive factor are added to the BAS algorithm to avoid the BAS algorithm getting trapped in local optima, while improving the convergence speed, solving the noise covariance uncertainty problem, and having good anti-interference ability.

[0167] Based on any of the above embodiments, another embodiment of the present invention provides an electronic device, which may include: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor may invoke logical instructions in the memory to execute the above method.

[0168] Furthermore, when the logical instructions in the aforementioned memory can be implemented as software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, essentially, or the part that contributes to the prior art, or a part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0169] On the other hand, embodiments of the present invention also provide a storage medium storing a plurality of instructions adapted for loading by a processor to execute the GPS / BDS dual-mode satellite positioning method based on joint LS and IBAS-AUKF calculations as provided in the above embodiments.

[0170] On the other hand, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0171] The device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Those skilled in the art can understand and implement this without any creative effort.

[0172] Through the above description of the embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus necessary general-purpose hardware platforms, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solutions, in essence or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the methods described in the various embodiments or some parts of the embodiments.

[0173] In summary, although the present invention has been disclosed above with reference to preferred embodiments, the above preferred embodiments are not intended to limit the present invention. Those skilled in the art can make various modifications and refinements without departing from the spirit and scope of the present invention. Therefore, the scope of protection of the present invention shall be determined by the scope defined in the claims.

Claims

1. A GPS / BDS dual-mode satellite positioning method based on joint calculation of LS and IBAS-AUKF, characterized in that, include: The digital intermediate frequency (IF) signals of BDS and GPS satellite signals are acquired, and the digital IF signals are mediated to obtain navigation information and pseudorange measurement values. GPS satellite positioning system model and BDS satellite positioning system model are established respectively. Based on the navigation information and pseudorange measurement values, the pseudorange equation of GPS satellite and pseudorange equation of BDS satellite are solved by pseudorange equation solving method. A combined positioning model is obtained based on the pseudorange equations of the GPS satellites and the BDS satellites. The combined positioning model is solved using the least squares method to obtain the first least squares estimate of the position solution; The Kalman filter is initialized using the first least squares estimate, and the improved unscented Kalman filter is selected as the dual-mode navigation filter for state estimation to obtain the state estimate. The relative position information of the receiver is obtained based on the state estimate.

2. The method as described in claim 1, characterized in that, The combined positioning model is as follows: Where (AX=b) is the least squares form of the combined localization model, [X i0 ,Y i0 Z i0 [This represents the approximate location of the receiver.] and Let ΔX represent the geometric distances from the GPS satellite and the BDS satellite to the receiver, respectively. i ,ΔY i ,ΔZ i It is the difference between the approximate and actual coordinates of the receiver, where c is the speed of light, and dt is the distance between the approximate and actual positions of the receiver. Gi and dt Ci These are the receiver clock errors for GPS satellites and BDS satellites, respectively, where G is the GPS satellite, C is the BDS satellite, and b = P. i -ρ i +c·dt G +c·dt C -T i -I i -ε i I i T i ε i These represent ionospheric delay, tropospheric delay, and other delay error terms, dt G It is the star clock error, P i ρ is the pseudorange measurement from the satellite to receiver i. i It is the geometric distance from the receiver to the satellite, [x Gm ,y Gm ,z Gm ] and [x Cn ,y Cn ,z Cn These are the GPS satellite positions and BDS satellite positions, respectively.

3. The method as described in claim 2, characterized in that, The step of solving the combined positioning model using the least squares method to obtain a first-order least squares estimate of the location solution includes: The GPS / BDS pseudorange linear equation is calculated according to the preset least squares calculation formula. The equation is iterated repeatedly to minimize b, and the position correction value that minimizes b is used as the first least squares estimate of the position solution. dx0=(A T A) -1 A T b; x1=x k-1 +dx0; In the formula, dx0 is the correction amount. The initial value is corrected using the correction amount dx0 to obtain the corrected value x1. k-1 Let A represent the state estimate at time k-1. A is a known matrix representing the coefficients in the localization model.

4. The method as described in claim 2, characterized in that, The selection of an improved unscented Kalman filter as a dual-mode navigation filter for state estimation to obtain state estimates includes: The initial position value calculated by least squares is assigned to the Kalman filter initialization, and the state estimate is obtained according to the following state equation and measurement equation: Take error state Establish discrete state equations: Where x is the actual state quantity at the location, Let F be the initial position value after initialization, and w be the state transition matrix. k The system excitation noise sequence has statistical properties derived from the process noise covariance matrix Q. k describe; The measurement equation is a linearized GPS / BDS pseudorange equation: δρ=H·δx(k)+v(k); Here, the row vector of H corresponds to the partial derivative of the observation equation of each satellite at the initial position, and v(k) is the measurement noise sequence.

5. The method as described in claim 4, characterized in that, The step of selecting an improved unscented Kalman filter as a dual-mode navigation filter for state estimation to obtain state estimates further includes: The system state at time k is updated using an improved unscented Kalman filter algorithm, and the state estimate is obtained based on the state update result.

6. The method as described in claim 5, characterized in that, The method of updating the system state at time k using an improved unscented Kalman filter algorithm includes: S1 is obtained by calculating the Sigma point and its weight using the following formula: In the formula, Let λ represent the i-th column of the square root of the matrix, where λ = α 2 (n+κ)-n, where λ is the scaling parameter, the value of α depends on the distribution of the sampling points, κ is a parameter, n is the dimension of the state vector, and X0 represents the initial state value. X represents the mean of a Gaussian distribution at time k-1. i (k-1) represents the weighting coefficients corresponding to the Sigma point, P k Represents the covariance matrix. β is 2, where, This represents the weights required to calculate the covariance. This represents the weights needed to calculate the mean. This represents the weights required to calculate the mean at the initial time step; S2, Calculate the one-step predicted state vector X (i) (k+1|k) and one-step predicted state covariance Among them, X (i) (k+1|k) represents the sample points for one-step prediction, X (i) (k|k) represents the i-th sampling point at time k. It is a one-step prediction of the estimated state at time k, u k-1 For deterministic control terms, Q k It is the process noise covariance. This represents the state prediction covariance at time k; The measurement update formula is as follows: Among them, y (i) (k+1|k) represents the measurement prediction sample points, P yy,k and P xy,k These are the prediction measurement covariance matrix and the cross covariance matrix, respectively. To estimate the state, P(k+1|k+1) is the estimated covariance matrix, and K... k+1 R is the gain matrix; k To observe the noise covariance, Let y(k+1) be the one-step predicted mean of the measured variable, y(k+1) be the measured value at time k+1, and P(k+1|k) be the estimated covariance matrix. S3, repeat S1-S2 until the preset iteration termination condition is reached, and obtain the LS-UKF joint estimated receiver position as the receiver relative position information.

7. The method as described in claim 6, characterized in that, The method of using an improved unscented Kalman filter algorithm to update the system state at time k also includes: The process noise covariance and the observation noise covariance are updated using the beetle whisker search algorithm to obtain the updated process noise covariance and observation noise covariance.

8. The method as described in claim 7, characterized in that, The update of the process noise covariance and the observation noise covariance using the beetle whisker search algorithm includes: Vector normalization was performed on the head orientation of the longhorn beetle: Where t represents the number of unknowns in the problem to be solved, and r is the result of vector normalization; Obtain the spatial coordinates of the longhorn beetle's left and right tentacles: In the formula, x t It is the location of the centroid of the longhorn beetle at the t-th iteration. These are the coordinates of the left tentacle's position during the t-th iteration. d is the position coordinate of the right tentacle during the t-th iteration. t It is the distance between the two tentacles of the longhorn beetle; Construct a fitness function f(·), using the absolute cumulative error value of the innovation as the fitness function f(·) of the BAS algorithm, to obtain the optimal process noise covariance Q. k and observation noise covariance R k : In the formula, the row vector of H corresponds to the partial derivative of the observation equation of each satellite at the initial position; The fitness values ​​of the locations of the left and right antennae of the longhorn beetle are calculated based on the fitness function f(·). and And update the location: In the formula, δ t Let x represent the search step size in the t-th iteration. t Let x be the location of the longhorn beetle at time t. t+1 Let t+1 be the location of the longhorn beetle; The distance d between the two tentacles of the longhorn beetle t and the search step size δ of the t-th iteration t Update: In the formula, d e It is the attenuation factor of the distance between the two whiskers, and α represents the step size reduction factor.

9. The method as described in claim 8, characterized in that, The fitness values ​​of the locations of the left and right antennae of the longhorn beetle are calculated based on the fitness function f(·). and After updating the location, the method also includes: The simulated annealing algorithm is used to determine whether the updated position meets the preset loop termination condition, including: The probability P corresponding to the fitness value of the location of the longhorn beetle is calculated according to the preset probability formula: In the formula, f(x) t+1 f(x) t ) are the fitness functions at different times, where T is the current temperature, and when f(x) t+1 )<f(x t When xt+1 is the current optimal solution, we can directly accept x. t+1 For the new solution, when f(x) t+1 )≥f(x t When x is given, accept the inferior solution x with probability P. t+ 1 represents the new solution, which is a new process noise covariance and observation noise covariance.

10. The method as described in claim 8, characterized in that, The step size reduction factor α is updated in the following manner to obtain the updated step size reduction factor β: In the formula, f t T is the fitness function in the BAS algorithm at the current iteration number t. max This represents the total number of iterations in the BAS algorithm, α is the original step size reduction factor, and f min It is the minimum fitness function value.

Images