Correction of floating point coordinate solution method for differential equation of RTK mobile station

Abstract

The present invention relates to a method for solving floating-point solution coordinates of a modified RTK mobile station differential equation, characterized in that the method comprises the following steps: (1) obtaining the moving speed of a receiver by using Doppler frequency shift velocity measurement; (2) calculating the average moving speed of the receiver between the last solution time and the current time by using the fixed solution coordinates obtained in the previous solution and the floating-point solution coordinates of the current RTK mobile station differential equation; and (3) correcting the floating-point solution coordinates of the current differential equation by using the difference between the moving speed of the receiver obtained in step (1) and the average moving speed calculated in step (2). The present invention uses Doppler frequency shift velocity measurement to obtain the precise speed of the receiver, and uses the difference between the two speeds to correct the floating-point solution coordinates of the current differential equation, so that the floating-point solution coordinates are closer to the fixed solution coordinates, thereby improving the efficiency of solving integer ambiguity.

Description

Technical Field

[0001] The present invention relates to a precise positioning technology of an RTK differential positioning system, in particular to a method for calculating floating-point solution coordinates of a modified RTK mobile station differential equation. Background Art

[0002] RTK positioning technology offers the advantages of high operational efficiency, high positioning accuracy, secure and reliable data, and minimal error accumulation. Furthermore, this positioning technology requires minimal operating conditions, is simple to operate, is highly integrated, and possesses strong data processing capabilities. However, RTK technology has some drawbacks. Limited by satellite status, satellite signals can be blocked or blocked for extended periods in densely populated urban areas or in mountainous canyon environments, causing satellite lock loss and impacting the initialization of integer ambiguity, resulting in low positioning accuracy or even failure. Furthermore, due to the influence of the sky environment, severe ionospheric interference can reduce the number of shared satellites between the rover and base stations, making positioning impossible. Furthermore, due to interference and limitations on the data transmission link, such as mountains and tall buildings, or interference from high-frequency signal sources within the data transmission link, the signal between the base and rover stations can be severely attenuated during transmission, impacting positioning accuracy.

[0003] The Doppler velocity measurement method is based on a strict Doppler velocity measurement mathematical model. The Doppler observation value is numerically expressed as the rate of change of the instantaneous carrier phase and is independent of the observation value of the carrier phase. It does not change due to cycle slips in the carrier phase. Therefore, the Doppler observation value is used to obtain the receiver movement speed.

[0004] The accuracy of Doppler velocity measurements is affected not only by the rate of change of the receiver clock error, the rate of change of the satellite clock error, and the rate of change of ionospheric and tropospheric delays, but also by observation noise. When the carrier is in motion, obstruction by surrounding objects such as buildings and trees typically results in greater noise in the observation data than when the carrier is stationary. Observation noise is a major source of error affecting velocity measurements. To achieve high-precision velocity, it is necessary to reduce the impact of observation noise on Doppler velocity measurements.

[0005] In actual RTK surveying, the RTK floating-point solution mode is primarily caused by the failure to resolve integer ambiguities when the GPS signal is poor or the star conditions are poor. The core of the RTK positioning algorithm is the carrier phase integer ambiguity fixation technique, but this technique is susceptible to factors such as navigation satellite failures and changes in the operating environment, which can lead to ambiguity fixation failures. Traditional RTK positioning algorithms select all carrier phase integer ambiguities fixed in the observation epoch. However, in actual applications, due to navigation satellite failures or deterioration of the satellite navigation system operating environment, the measurement accuracy of the satellite navigation observations decreases, making it impossible for the LAMBDA algorithm to completely and correctly fix all carrier phase integer ambiguities.

[0006] Variance, as a measure of observation accuracy, also reflects the closeness of the ambiguity float solution to the fixed solution. A smaller variance indicates a closer ambiguity float solution to the fixed solution, and a greater likelihood that the integer ambiguity solution will be correctly fixed. Summary of the Invention

[0007] In order to solve the above problems, the present invention proposes a method for solving the floating-point solution coordinates of the RTK mobile station differential equation to make the floating-point solution coordinates closer to the fixed solution coordinates, making the receiver's precise positioning more stable.

[0008] In order to achieve the above-mentioned object, the present invention proposes a method for solving the floating-point solution coordinates of the modified RTK mobile station differential equation, which is characterized by comprising the following steps: (1) using Doppler frequency shift velocity measurement to obtain the moving speed of the receiver; (2) using the fixed solution coordinates obtained in the last solution and the floating-point solution coordinates of the current RTK mobile station differential equation to calculate the average moving speed of the receiver between the last solution time and the current time; (3) using the difference between the moving speed of the receiver obtained in step (1) and the average moving speed calculated in step (2) to correct the floating-point solution coordinates of the current differential equation.

[0009] The present invention adopts Doppler frequency shift velocity measurement to obtain the precise velocity of the receiver, and uses the difference between the two velocities to correct the floating-point solution coordinates of the current differential equation, so that the floating-point solution coordinates are closer to the fixed solution coordinates, thereby improving the efficiency of resolving integer ambiguity. In this way, in low-dynamic conditions, the RTK differential system does not add other auxiliary positioning technologies, fully utilizes the data information obtained from satellite signals, reduces the impact of complex environments on positioning accuracy, and realizes accurate real-time positioning of the carrier. DETAILED DESCRIPTION

[0010] The method for calculating the floating-point solution coordinates of the differential equation of the modified RTK mobile station includes the following steps.

[0011] (1) Use Doppler frequency shift to measure the speed of the receiver.

[0012] A stationary signal tower broadcasts a signal with a frequency of f. If the carrier mobile receiving station moves at a speed of v, then the signal frequency received by the carrier mobile receiving station is f. r It is not the transmission frequency f of the signal, but f+f d , this signal receiving frequency f r The phenomenon that changes with the relative motion between the signal transmitting source and the carrier mobile receiving station is called the Doppler effect. d It is called the Doppler frequency shift value. d Equal to the frequency f of the satellite signal received by the receiver r The difference between the frequency f and the satellite's signal transmission frequency is:

[0013] f d =f r -f (1)

[0014] Based on the basic theory of electromagnetic propagation, the Doppler frequency shift value f d The calculation formula is:

[0015]

[0016] Among them, v s is the operating speed of the satellite; λ is the wavelength of the carrier signal transmitted by the satellite; v is the speed of the receiver; is the rate of change of the distance from the receiver to the satellite.

[0017] The pseudorange observation equation is:

[0018]

[0019] Among them, ρ n is the pseudorange observation value from the nth satellite to the receiver station; r n is the true geometric distance from the nth satellite to the receiver; δt r is the receiver clock error; δt n is the clock error of the nth satellite; I n is the ionospheric delay on the propagation path from the nth satellite to the receiver; T n is the tropospheric delay on the propagation path from the nth satellite signal to the receiver; is the sum of errors, including satellite orbit error, multipath error and receiver internal noise.

[0020] Taking the derivative of the pseudorange observation equation with respect to time, we get:

[0021]

[0022] in, is the rate of change of the pseudorange observation value from the nth satellite to the receiver station; is the rate of change of the geometric distance from the nth satellite to the receiver; is the unknown receiver clock frequency drift; is the clock frequency drift of the nth satellite; is the rate of change of ionospheric delay; is the rate of change of tropospheric delay.

[0023] The rate of change of the geometric distance from the nth satellite to the receiver The relationship between the speed of the receiver relative to the satellite is:

[0024]

[0025] v n is the speed of the nth satellite; l n is the direction cosine from the receiver to the nth satellite; v = (v x ,v y ,v z ) is the moving speed of the receiver to be solved.

[0026] The receiver receives the satellite signal, captures the pseudo code phase and Kepler frequency shift, and calculates the satellite velocity, direction cosine from the receiver to the satellite, and satellite clock velocity offset data obtained from the navigation message to solve the three-dimensional coordinates of the Earth-centered Earth-fixed coordinate system of the nth satellite as (x n ,y n ,z n ).

[0027] The receiver's three-dimensional floating-point coordinate position (x, y, z) in the Earth-centered Earth-fixed coordinate system is obtained by differentially solving the observation data transmitted from the base station. Calculate the length of the observation vector from the nth satellite to the receiver:

[0028]

[0029] The direction cosines from the receiver to the nth satellite are:

[0030]

[0031] The receiver accurately measures the Doppler shift of the nth satellite Can reflect the rate of change of pseudorange observations The size relationship between the two is:

[0032] Substituting equations (4) and (5) into equation (8), the ionospheric delay variation rate and the tropospheric delay variation rate can be ignored. The constant speed equation is obtained:

[0033]

[0034] The unknown quantities in the above formula are v, and The other parameters are known.

[0035] When the receiver receives N (N≥4) satellite signals and obtains N original Doppler measurement values, the N equations form a constant speed equation matrix equation:

[0036]

[0037] In the formula

[0038]

[0039] Among them, G is obtained by the coordinates of the satellite and the receiver in the geodetic coordinate system, and the satellite speed and the direction cosine l from the receiver to the satellite are calculated using the navigation message. n , the satellite clock difference velocity offset is calculated

[0040] When the number of satellites exceeds four, the equation for the receiver velocity is obtained by Kalman filtering and least squares method:

[0041]

[0042] Use Kalman filtering to solve the speed and reduce the noise impact. Kalman filtering is a recursive linear unbiased minimum variance estimation method that uses the state transition equation of the system to recursively estimate the new state estimate based on the state estimate at the previous moment and the observation value at the current moment, thereby reducing the total error. impact.

[0043] Assume the state quantity at time k

[0044]

[0045] In the formula, (v x v y v z ) is the speed at time k, (a x a y a z ) is the acceleration at time k.

[0046] Initialize the state x = [0,0,0,0,0,0,0], the posterior mean square error matrix P k Initialized as a 7-order square matrix with a main diagonal value of 100.

[0047] The satellite elevation angle weighted model is used to set the measurement noise n k The covariance matrix R is:

[0048]

[0049] Where, coefficient a is the carrier phase noise coefficient. According to the empirical value, coefficient a=3mm,θ i is the elevation angle of the ith satellite, then the measurement noise is only a function of the satellite angle. Since the time difference will amplify the noise variance to twice the original value, the coefficient needs to be multiplied by 2.

[0050] The process noise of each satellite is selected to obey a normal distribution that is independent of each other and has a mean of 0 and a variance of 0.005, that is,

[0051]

[0052] Establish the state transfer equation:

[0053]

[0054] Where x k|k-1 State prediction value, x k-1|k-1 is the best estimate of the Kalman filter at time t-1, w k-1 The variance is The process noise, The state transition matrix is ​​7×7:

[0055]

[0056] Where I represents the identity matrix and Δt represents the sampling interval.

[0057] Establish the observation equation:

[0058] z k =Hx k|k +n k (14)

[0059] Where H is the measurement matrix, x k|k is the Kalman filter error estimate at time t-1, n k The variance is The measurement noise, is the measurement residual.

[0060] The calculation formula of the measurement matrix H is:

[0061]

[0062] Assuming that the noise of each satellite observation is the same, calculate the last error covariance P k|k-1 The calculation formula is:

[0063]

[0064] Calculate the Kalman gain K k :

[0065] K k =P k|k-1 H T (HP k|k-1 H T +R) -1 (16)

[0066] Among them, P k|k-1 is the state quantity estimation covariance matrix, which is calculated from the measurement value at time k-1, and R is the measurement noise nk Covariance matrix.

[0067] The state update equation is corrected and updated by Kalman filtering to update the state quantity x at time k. k|k :

[0068] x k|k =x k|k-1 +K k ·[z k -Hx k|k-1 ] (17)

[0069] Where K k Kalman filter gain, z k is the actual error measurement value, x k|k-1 Predict the state value for one step.

[0070] Linearize the state equation to get the covariance matrix P k|k , update the covariance matrix P k|k :

[0071] P k|k =(lK k H)P k|k-l (18)

[0072] Where I is the 7*7 identity matrix.

[0073] By iterating the above Kalman filter prediction and correction recursive process, the receiver's velocity state v can be solved as [v x v y v z ].

[0074] (2) Using the fixed solution coordinates of the previous solution and the floating-point solution coordinates of the current RTK mobile station differential equation, calculate the average movement speed of the receiver between the previous solution time and the current time.

[0075] The receiver saves the last fixed solution coordinate r1 = [X1 Y1 z1] T The floating-point solution coordinates of the RTK mobile station differential equation at the corresponding time t1 and t2 Calculate the average velocity of the receiver between time t1 and time t2 in the WGS-84 coordinate system:

[0076]

[0077] (3) Using the difference between the moving speed of the receiver calculated in step (1) and the average moving speed calculated in step (2), the floating-point solution coordinates of the current differential equation are corrected.

[0078] When the interval between t1 and t2 is small and the dynamic situation of the carrier is not very large, the moving speed v of the receiver obtained by Doppler frequency shift velocity measurement in step (1) is [v x v y v z ] can be used as the receiver’s precise instantaneous velocity. The receiver’s moving velocity v = [v x v y v z ] and the average velocity between time t1 and time t2 in the WGS-84 coordinate system calculated in step (2) The speed error between them is:

[0079]

[0080] The relationship between speed error and position error correction is:

[0081]

[0082] Rewriting the above formula into the component form in the WSG-84 coordinate system is:

[0083]

[0084] Thus, the corrected floating point solution coordinates r2 = [X2 Y2 Z2] are obtained T , use the corrected floating-point solution coordinates to search for integer ambiguities.

[0085] The present invention adopts Doppler frequency shift to measure speed and reduces the influence of noise through Kalman filtering to obtain the precise speed of the receiver. At the same time, the receiver uses the fixed solution coordinates calculated last time and the floating-point solution coordinates of the current differential equation to calculate the average movement speed of the receiver. The difference between the two speeds is used to correct the floating-point solution coordinates of the current differential equation, so that the floating-point solution coordinates are closer to the fixed solution coordinates, thereby improving the efficiency of solving the integer ambiguity. In this way, in low-dynamic conditions, the RTK differential system does not add other auxiliary positioning technologies, fully utilizes the data information obtained from the satellite signal, reduces the influence of the complex environment on the positioning accuracy, and realizes the precise real-time positioning of the carrier.

Claims

1. A method for calculating floating-point solution coordinates of a modified RTK mobile station differential equation, characterized in that: The following steps are involved: (1) Using Doppler frequency shift to measure the speed of the receiver; (2) Using the fixed solution coordinates calculated at the last solution time and the floating-point solution coordinates of the RTK mobile station differential equation at the current time, calculate the average motion speed of the receiver between the last solution time and the current time; (3) Using the difference between the moving speed of the receiver obtained in step (1) and the average moving speed calculated in step (2), the floating-point solution coordinates of the current differential equation are corrected.

2. The method for solving the floating-point solution coordinates of the modified RTK mobile station differential equation according to claim 1, characterized in that: In step (2), the fixed solution coordinates r1 = [X1Y1Z1] calculated at the last solution time t1 T , the floating point solution coordinates of the difference equation at the current time t2 The average speed of the receiver between t1 and t2 Calculate as follows: The moving speed v of the receiver obtained in step (1) and the average speed The speed error between them is: Correct position error with velocity error: Rewriting the above formula into the component form in the WSG-84 coordinate system is:

3. The method for solving the floating-point solution coordinates of the modified RTK mobile station differential equation according to claim 2, characterized in that: In step (1), a stationary signal transmission tower broadcasts a signal with a frequency of f. If the carrier mobile receiving station moves at a speed of v, the signal frequency f received by the carrier mobile receiving station is r It is not the transmission frequency f of the signal, but f+f d , this signal receiving frequency f r The phenomenon that changes with the relative motion between the signal transmitting source and the carrier mobile receiving station is called the Doppler effect. d It is called the Doppler frequency shift value, the Doppler frequency shift value f d Equal to the frequency f of the satellite signal received by the receiver r The difference between the frequency f and the satellite's signal transmission frequency is: f d =f r -f (1) Based on the basic theory of electromagnetic propagation, the Doppler frequency shift value f d The calculation formula is: Among them, v s is the operating speed of the satellite; λ is the wavelength of the carrier signal transmitted by the satellite; v is the speed of the receiver; is the rate of change of the distance from the receiver to the satellite; The pseudorange observation equation is: Among them, ρ n is the pseudorange observation value from the nth satellite to the receiver station; r n is the true geometric distance from the nth satellite to the receiver; δt r is the receiver clock error; δt n is the clock error of the nth satellite; I n is the ionospheric delay on the propagation path from the nth satellite to the receiver; T n is the tropospheric delay on the propagation path from the nth satellite signal to the receiver; is the sum of errors, including satellite orbit error, multipath error and receiver internal noise; Taking the derivative of the pseudorange observation equation with respect to time, we get: in, is the rate of change of the pseudorange observation value from the nth satellite to the receiver station; is the rate of change of the geometric distance from the nth satellite to the receiver; is the unknown receiver clock frequency drift; is the clock frequency drift of the nth satellite; is the rate of change of ionospheric delay; is the rate of change of tropospheric delay; The rate of change of the geometric distance from the nth satellite to the receiver The relationship between the speed of the receiver relative to the satellite is: r n =(v n -v)·l n (5) v n is the speed of the nth satellite; l n is the direction cosine from the receiver to the nth satellite; v = (v x ,v y ,v x ) is the moving speed of the receiver to be solved; The receiver receives the satellite signal, captures the pseudo code phase and Kepler frequency shift, and calculates the satellite velocity, direction cosine from the receiver to the satellite, and satellite clock velocity offset data obtained from the navigation message to solve the three-dimensional coordinates of the Earth-centered Earth-fixed coordinate system of the nth satellite as (x n ,y n ,z n ); The receiver's three-dimensional floating-point coordinate position (x, y, z) in the Earth-centered Earth-fixed coordinate system is obtained by differentially solving the observation data transmitted from the base station to obtain the three-dimensional floating-point coordinate position (x, y, z) in the Earth-centered Earth-fixed coordinate system. The length of the observation vector from the nth satellite to the receiver is calculated as follows: The direction cosines from the receiver to the nth satellite are: The receiver accurately measures the Doppler shift of the nth satellite Can reflect the rate of change of pseudorange observations The size relationship between the two is: Substituting equations (4) and (5) into (8), the ionospheric delay variation rate and the tropospheric delay variation rate can be ignored, and the constant speed equation is obtained: The unknown quantities in the above formula are v, and All other parameters are known; When the receiver receives N (N≥4) satellite signals and obtains N original Doppler measurement values, the N equations form a constant speed equation matrix equation: In the formula Among them, G is obtained by the coordinates of the satellite and the receiver in the geodetic coordinate system, and the satellite speed and the direction cosine l from the receiver to the satellite are calculated using the navigation message. n , the satellite clock difference velocity offset is calculated When the number of satellites exceeds four, the equation for the receiver velocity is obtained by Kalman filtering and least squares method: Use Kalman filtering to solve the speed and reduce the influence of noise. Kalman filtering is a recursive linear unbiased minimum variance estimation method. Through the state transition equation of the system, the new state estimate is recursively estimated based on the state estimate at the previous moment and the observation value at the current moment, thereby reducing the total error. the impact of; Assume the state quantity at time k In the formula, (v x v y v z ) is the speed at time k, (a x a y a z ) is the acceleration at time k; Initialize the state x = [0,0,0,0,0,0,0], the posterior mean square error matrix P k Initialize it to a 7-order square matrix with a main diagonal value of 100; The satellite elevation angle weighted model is used to set the measurement noise n k The covariance matrix R is: Where, coefficient a is the carrier phase noise coefficient. According to the empirical value, coefficient a=3mm,θ i is the elevation angle of the ith satellite, then the measurement noise is only a function of the satellite angle. Since time difference will amplify the noise variance to twice the original value, the coefficient needs to be multiplied by 2; The process noise of each satellite is selected to obey a normal distribution that is independent of each other and has a mean of 0 and a variance of 0.005, that is, Establish the state transfer equation: Where x k|k-1 State prediction value, x k-1|k-1 is the best estimate of the Kalman filter at time t-1, w k-1 The variance is The process noise, The state transition matrix is ​​7×7: Where I represents the identity matrix, Δt represents the sampling interval; Establish the observation equation: z k =Hx k|k +n k (14) Where H is the measurement matrix, x k|k is the Kalman filter error estimate at time t-1, n k The variance is The measurement noise, is the measurement residual; The calculation formula of the measurement matrix H is: Assuming that the noise of each satellite observation is the same, calculate the last error covariance P k|k-1 The calculation formula is: Calculate the Kalman gain K k : K k =P k|k-1 H T (HP k|k-1 H T +T) -1 (16) Among them, P k|k-1 is the state quantity estimation covariance matrix, which is calculated from the measurement value at time k-1, and R is the measurement noise n k covariance matrix; The state update equation is corrected and updated by Kalman filtering to update the state quantity x at time k. k|k : x k|k =x k|k-1 +K k ·[z k -Hx k|k-1 ] (17) Where K k Kalman filter gain, z k is the actual error measurement value, x k|k-1 Predict the state value for one step; Linearize the state equation to obtain the covariance matrix p k|k , update the covariance matrix P k|k : P k|k =(I-K k H)P k|k-1 (18) Where I is the 7*7 identity matrix; By iterating the above Kalman filter prediction and correction recursive process, the receiver's velocity state v can be solved as [v x v y v z ].

Images