Positioning processing method suitable for observation data under complex conditions of GNSS (Global Navigation Satellite System)

Abstract

The invention discloses a positioning processing method suitable for observation data under a GNSS complex condition. The positioning processing method comprises the following steps: step 1, setting residual check thresholds of a pseudo-range positioning observation value and a Doppler velocity observation value; step 2, on the basis of a residual check threshold of a pseudo-range positioning observation value, calculating an optimal estimated value of pseudo-range positioning observation and a corresponding variance-covariance matrix according to residual satellites for satellites with abnormal pseudo-range positioning observation values; step 3, based on a residual error test threshold value of the Doppler velocity observation value and a detection result in the step 2, eliminating satellites with abnormal Doppler velocity observation values, and calculating an optimal estimated value of Doppler velocity observation and a corresponding variance-covariance matrix according to the residual satellites; and step 4, updating the Kalman filtering state vector of the GNSS based on the results of the step 2 and the step 3. The method can improve the estimation precision of the initial values of the position and the speed, and further improves the positioning precision of the GNSS.

Description

Technical Field

[0001] The present invention relates to a GNSS observation data processing method, and in particular to a positioning processing method for GNSS observation data under complex conditions. Background Art

[0002] GNSS (Global Navigation Satellite System) plays a vital role in modern navigation and positioning, especially in real-time kinematic (RTK) and precise point positioning (PPP) applications. However, with the development of cities, the demand for high-precision positioning using low-cost GNSS systems is increasing. In complex urban environments, satellite signals are frequently lost due to obstruction by tall buildings, tree shades and other objects, resulting in problems such as gross observation errors and non-line-of-sight signals, which in turn affects the performance of GNSS precision positioning.

[0003] Under complex observation conditions, data preprocessing and quality control are crucial for high-precision positioning of RTK and PPP. In recent years, research on data preprocessing and quality control has continued to emerge, including multipath effect processing, signal shielding processing, ionospheric and tropospheric delay correction, observation noise filtering and data fusion, and data quality assessment and control. Multipath effect processing is mainly achieved by detecting and eliminating multipath signals, designing anti-multipath antennas, and modeling and correcting multipath errors. Signal shielding processing relies on shielding detection and compensation technology and multi-constellation combined positioning methods to improve positioning stability. The correction of ionospheric and tropospheric delays uses dual-frequency observation data, ionospheric and tropospheric models, and real-time meteorological data to improve correction accuracy. In addition, filtering algorithms such as Kalman filtering, particle filtering, and multi-source data fusion technology play an important role in reducing observation noise. Data quality assessment and control can ensure the accuracy and reliability of positioning results through real-time monitoring of observation data quality, residual analysis, and multipath effect detection and elimination.

[0004] However, the above research mainly focuses on the optimization of hardware and GNSS precision positioning models. The accuracy and convergence performance of precision positioning depend on the quality of observations on the one hand and the accuracy of initial positioning values ​​on the other. Affected by the observation environment, pseudorange positioning and Doppler velocity observations under complex observation conditions will have anomalies, resulting in poor accuracy of the initial position and velocity values ​​of the pseudorange single-point positioning solution, which in turn leads to abnormal Kalman filter state updates for PPP and RTK precision positioning, and ultimately leads to significant deviations in receiver state changes between epochs. Therefore, eliminating abnormal observations and improving positioning accuracy through the data preprocessing stage is of great significance to promoting the advancement and popularization of high-precision positioning technology. Summary of the invention

[0005] The present invention provides a positioning processing method suitable for GNSS observation data under complex conditions, so as to solve the problem that the observation data obtained by the prior art GNSS has abnormalities due to the complex observation conditions, resulting in poor accuracy of initial position and velocity solution, and further leading to abnormal Kalman filter state update.

[0006] In order to achieve the above object, the technical solution adopted by the present invention is:

[0007] The positioning processing method for observation data under complex GNSS conditions includes the following steps:

[0008] Step 1: setting a residual check threshold of a pseudorange positioning observation value and a residual check threshold of a Doppler velocity observation value for a satellite receiver;

[0009] Step 2: Based on the residual test threshold of the pseudorange positioning observation value obtained in step 1, after finding and eliminating the satellites with abnormal pseudorange positioning observation values ​​among the satellites corresponding to the receiver, the optimal estimate of the pseudorange positioning observation and the corresponding variance covariance matrix are calculated according to the pseudorange positioning observation values ​​of the remaining satellites;

[0010] Step 3: Based on the residual test threshold of the Doppler velocity observation value obtained in step 1 and the detection result of the satellite with abnormal pseudorange positioning observation value in step 2, find and eliminate the satellite with abnormal Doppler velocity observation value among the satellites corresponding to the receiver, and then calculate the optimal estimate of the Doppler velocity observation and the corresponding variance covariance matrix according to the Doppler velocity observation values ​​of the remaining satellites;

[0011] Step 4: Based on the optimal estimate of the pseudorange positioning observation obtained in step 2, and using the optimal estimate of the Doppler velocity observation obtained in step 3 as a constraint, a state matrix is ​​obtained, thereby updating the Kalman filter state vector of the GNSS.

[0012] Furthermore, in the residual check threshold of the pseudorange positioning observation value and the residual check threshold of the Doppler velocity observation value described in step 1, each residual check threshold is composed of an adaptive factor and an empirical threshold constant, wherein the adaptive factor is determined according to the number of satellites in two adjacent epochs of the receiver, and the empirical threshold constant is set based on experience according to the type of receiver or the observation quality.

[0013] Furthermore, in step 2, first, a mathematical model for pseudo-range single-point positioning of multiple satellites is established, and the mathematical model of pseudo-range single-point positioning is corrected for spherical delay, ionospheric delay, and satellite clock error to obtain a normal equation matrix; then, the normal equation matrix is ​​solved to obtain the posterior residual of the pseudo-range positioning observation value of each satellite; then, based on the posterior residual of the pseudo-range positioning observation value of each satellite and the corresponding covariance information, the standardized residual of the pseudo-range positioning observation value of each satellite is obtained; finally, according to the standardized residual of the pseudo-range positioning observation value of each satellite and the residual check threshold of the pseudo-range positioning observation value obtained in step 1, the satellites with abnormal pseudo-range positioning observation values ​​among the satellites corresponding to the receiver are eliminated.

[0014] Furthermore, in step 3, firstly, a Doppler observation method equation is established for the receiver to perform Doppler velocity observations on multiple satellites; then, the Doppler observation method equation is solved to obtain the posterior residual of the Doppler velocity observation value of each satellite; then, based on the posterior residual of the Doppler velocity observation value of each satellite, the residual test threshold of the Doppler velocity observation value obtained in step 1 is used, combined with the detection result of the satellite with abnormal pseudorange positioning observation value in step 2, to eliminate the satellites with abnormal Doppler velocity observation values ​​among the satellites corresponding to the receiver.

[0015] Furthermore, in step 4, a consistency check is also performed based on the Euclidean geometric distance between the pseudorange positioning observation result of step 2 and the Doppler velocity observation result of step 3.

[0016] The present invention proposes a preprocessing method for GNSS observation data under complex observation conditions such as tree cover and urban canyons. The method sets a residual check threshold of a pseudorange positioning observation value and a residual check threshold of a Doppler velocity observation value for a receiver, thereby eliminating corresponding abnormal values ​​through a multidimensional gross error detection method in the pseudorange positioning and Doppler velocity measurement stages, thereby improving the estimation accuracy of initial position and velocity values. Therefore, when the Kalman filter state is updated, the obtained optimal estimate of the Doppler velocity observation can be used to perform an inter-epoch state update of precise positioning, thereby improving the positioning accuracy of the GNSS. DETAILED DESCRIPTION

[0017] The present invention is further described below in conjunction with embodiments.

[0018] This embodiment discloses a positioning processing method for observation data under complex GNSS conditions, including the following steps:

[0019] Step 1: Set the residual error test threshold τ of the pseudo-range positioning observation value for the satellite receiver rc , the residual test threshold τ of Doppler velocity observation rd .

[0020] Since different types of receivers, such as survey-grade receivers, low-cost professional GNSS receivers, and mass-market devices such as smartphones, exhibit different observation qualities and noise levels, using a fixed threshold for outlier detection may not yield stable and reliable positioning results. To address this issue and ensure effective identification and removal of measurement outliers, this embodiment proposes an adaptive and dynamically adjustable residual test threshold τ r , the residual test threshold τ r As shown in formula (1):

[0021] τ r =α·c0 (1)

[0022] In formula (1), α is the adaptive factor; c0 is an empirical threshold constant related to the receiver noise level of pseudorange, Doppler and carrier phase residuals, which is set according to the receiver type or observation quality. The adaptive factor α is shown in formula (2):

[0023]

[0024] In formula (2): Nsat k-1 and Nsat k Refers to the number of satellites in two adjacent epochs of the receiver. In challenging environments where the number of satellites tracked by the receiver is limited or the observation quality is poor, using a larger residual threshold can prevent more satellites from being eliminated to maintain the continuity of the positioning result. Therefore, this embodiment constructs an adaptive and dynamically adjustable residual test threshold τ r When , the adaptive factor α is determined according to the number of satellites in two adjacent epochs of the receiver.

[0025] The empirical threshold constant c0 is shown in formula (3):

[0026]

[0027] In formula (3), Rtype represents the receiver type, where Rtype = 0 represents a measurement receiver, which mainly refers to a receiver that can track satellite observation data of multiple systems and multiple frequency bands and has good observation value quality; Rtype = 1 represents a popular receiver, which mainly refers to tracking single-frequency or L1 / L5 dual-frequency observation data, and is generally used in vehicle-mounted environments; Rtype = 2 represents the observation value output by the receiver chip inside the smartphone.

[0028] Based on the above process, in this embodiment, the user needs to pre-set the receiver type Rtype before processing the data. In the subsequent processing, the adaptive factor α is calculated according to the number of satellites tracked in the current epoch and the previous epoch, and the empirical threshold constant c0 is determined according to the receiver type Rtype. Then, the residual test threshold τ of the pseudorange positioning observation value can be calculated by formula (1). rc and the residual test threshold τ of the Doppler velocity observation rd .

[0029] Step 2: Based on the residual error test threshold τ of the pseudo-range positioning observation value obtained in step 1 rc After finding and removing the satellites with abnormal pseudo-range positioning observation values ​​among the satellites corresponding to the receiver, the optimal estimate of the pseudo-range positioning observation and the corresponding variance covariance matrix are calculated based on the pseudo-range positioning observation values ​​of the remaining satellites. The specific process is as follows:

[0030] Pseudo-range single point positioning (SPP) can provide initial values ​​for real-time differential positioning (RTK) and precise point positioning (PPP) solutions. First, the mathematical model of pseudo-range single point positioning is established as shown in formula (4):

[0031] P=ρ+c·(δt r -δt s )+I i +T i +ε i (4)

[0032] In formula (4), P is the pseudorange measurement in meters; ρ is the geometric distance between the phase center of the satellite antenna when the signal is transmitted and the phase center of the receiver antenna when the signal is received, in meters; δt r and δt s I represents the receiver clock error and satellite clock error respectively; i and T i denote the ionospheric delay and tropospheric delay respectively; ε i represents observation noise; c is the speed of light.

[0033] After correcting the spherical delay, ionospheric delay and satellite clock error for the mathematical model of pseudo-range single point positioning, the normal equation matrix can be obtained as shown below:

[0034]

[0035] The above normal equation matrix formula can be further simplified as shown in formula (5):

[0036] l p =A·X p +ε p ,Q p (5)

[0037] In formula (5): p represents the satellite observation prior residual vector, l p =(lp1, lp2, ..., lpn) T , lp1, lp2, ..., lpn are the prior residual vectors of the observation values ​​of the 1st to nth satellites;

[0038] represents the design matrix, where (x1,y1,z1)…(x n ,y n ,z n ) represent the spatial rectangular coordinates of the 1st to nth satellites observed by the observation station; (x r ,y r ,z r ) represents the approximate value of the spatial rectangular coordinates of the user's receiver position; ρ1, ρ2…, ρ n The geometric distances calculated for the approximate coordinates of the receiver and the coordinates of the 1st to nth satellites respectively;

[0039] X p represents the correction value of the user position and receiver clock error to be estimated, X p =(x r ,y r ,z r ,ct r ), where x r ,y r ,z r Indicates the user coordinates to be estimated; ct r Indicates the correction value of the receiver clock error.

[0040] ε p represents the pseudorange observation noise, ε p =(∈ p1 ,∈ p2 ,...,∈ pn ),∈ p1 ,∈ p2 ,...,∈ pn represents the observation noise of the pseudorange observation value of the nth satellite;

[0041] Q p Represents the variance-covariance information of the pseudorange observations.

[0042] Then, based on the adjustment algorithm, the simplified form of the normal equation shown in formula (5) is solved to obtain the estimated parameter X p The optimal valuation vector The optimal estimate vector Substituting (5) into the posterior residual matrix v of the pseudorange positioning observation value can be calculated: p, the posterior residual matrix v of the pseudorange positioning observation p As shown in formula (6):

[0043]

[0044] In formula (6): v p Represents the posterior residual matrix of the pseudorange positioning observations.

[0045] Assume that there are ns satellites in each epoch, namely s1,s2,...,s n , then according to formula (6) we can get the i-th satellite s i The posterior residual v of the pseudorange positioning observation pi , as shown in the following formula (7):

[0046]

[0047] Where: l pi represents the OMC observation value of the i-th satellite; A i represents the i-th row of the design matrix coefficient A; Represents the parameter to be estimated X p The optimal estimate vector of .

[0048] Then, based on the posterior residual v of each satellite pi and its corresponding covariance information σ pi , and then obtain the standardized residual and where σ pi Determined based on the pseudorange observation noise and its corresponding altitude angle.

[0049] Finally, according to the standardized residual of each satellite pseudorange positioning observation And the residual test threshold τ of the pseudorange positioning observation value obtained in step 1 rc , find and remove the satellites with abnormal pseudo-range positioning observation values ​​among the satellites corresponding to the receiver. The specific process is as follows:

[0050] (2.1) Create a list of satellites with abnormal pseudorange observation values, sat_outlier, and a list of pseudorange satellites, sat_outlier_code, that have been detected and marked as gross errors.

[0051] Initialize the satellite list sat_outlier with abnormal pseudorange observation values ​​and the satellite list sat_outlier_code with gross errors of pseudorange to be empty;

[0052] Let the satellite with the largest residual in each epoch be sat_max, and the corresponding maximum residual be v max, initialize the satellite number sat_max with the largest residual in each epoch as sat_max = 0; the corresponding maximum residual is v max = 0.

[0053] (2.2) Based on the residual test threshold τ obtained in step 1 for the pseudorange positioning observation values related to the receiver rc , perform a residual test on each satellite in turn. The process is as follows:

[0054] For the i-th satellite, if then it is considered that the pseudorange residual of the i-th satellite s i is abnormal, and the satellite number of the i-th satellite s i is added to the satellite list sat_outlier of abnormal pseudorange observations;

[0055] If the pseudorange residual of the i-th satellite s i simultaneously satisfies then update and record its number sat_max = i;

[0056] For all satellites i = 1, 2, …… n, repeat steps (2.1)-(2.2) to test the residuals of each satellite. After the residual test, the satellite list sat_outlier of abnormal pseudorange observations contains all the satellites marked as having abnormal pseudorange residuals; find the maximum residual and the corresponding satellite number among them, and store them in v max and sat_max respectively;

[0057] (2.3) If |v max | > 10 * τ rc , then it is considered that the satellite with the satellite number sat_max corresponding to v max has a significant abnormality; at this time, delete the remaining satellites from the satellite list sat_outlier of abnormal pseudorange observations, and only represent the satellite with the number sat_max as abnormal and add it to the gross error satellite list sat_outlier_code of pseudorange.

[0058] If it satisfies |v max | < 10 * τ rc , and the total number of satellites in the satellite list sat_outlier of abnormal pseudorange observations < ns / 3, where ns is the total number of satellites in this epoch, then all the satellites in the satellite list sat_outlier of abnormal pseudorange observations are added to the gross error satellite list sat_outlier_code of pseudorange.

[0059] If it satisfies |v max | < 10 * τ rc, and if the satellites in the satellite list sat_outlier with abnormal pseudorange observation values ​​are >ns / 3, then determine the navigation system to which each satellite in the satellite list sat_outlier with abnormal pseudorange observation values ​​belongs, and the determination result is one of GPS, GAL (Galileo), GLO (GLONASS), and BDS, and add one to the gross error count GPS_noutlier, GAL_noutlier, BDS_noutlier, and GLO_noutlier of the corresponding navigation system. After the cyclic determination of all satellites in the satellite list sat_outlier with abnormal pseudorange observation values ​​is completed, compare the maximum values ​​of GPS_noutlier, GAL_noutlier, BDS_noutlier, and GLO_noutlier, and add the satellite number corresponding to the maximum value to the gross error satellite list sat_outlier_code of pseudorange;

[0060] (2.4) Remove the satellites in the sat_outlier_code list of the gross error satellites added to the pseudorange in step (2.3), and then use the least squares method to solve the residual matrix of the remaining satellites according to formula (6); and repeat steps (2.1)-(2.4) for iteration.

[0061] After the above iterative steps (2.1)-(2.4), the estimated parameter X is obtained by formula (7) after the iteration. p The optimal valuation vector That is, the optimal estimate of the pseudorange positioning measurement value can be obtained and its corresponding variance-covariance information Q p .

[0062] Step 3: Based on the residual test threshold of the Doppler velocity observation value obtained in step 1 and the detection result of the satellite with abnormal pseudo-range positioning observation value in step 2, find and remove the satellites with abnormal Doppler velocity observation value among the satellites corresponding to the receiver, and then calculate the optimal estimate of the Doppler velocity observation and the corresponding variance covariance matrix based on the Doppler velocity observation value of the remaining satellites. The specific process is as follows:

[0063] The mathematical model of the receiver's Doppler velocity observation of multiple satellites is shown in formula (8):

[0064]

[0065] In formula (8), subscripts r and s are receiver r and satellite s respectively; f d is the Doppler observation value corresponding to the phase observation value; e is the receiver-satellite line-of-sight unit vector; and Represent the speed of satellite and receiver respectively; c represents the speed of light; represents the receiver clock drift; ε represents the multipath and noise of Doppler observation. Among them, the satellite speed and satellite clock error The parameters can be calculated using the broadcast ephemeris parameters broadcast by the satellite. Therefore, the receiver velocity parameters for a single epoch can be solved based on the Doppler velocity observations.

[0066] In one observation epoch, the receiver can synchronously observe multiple satellites. Therefore, this embodiment first establishes the Doppler observation method equation for the receiver to observe the Doppler velocity of multiple satellites, as shown in formula (9):

[0067] l d =A·X v +ε d ,Q d (9)

[0068] In formula (9): d is the OMC (Observed-Minus-Computed) observation vector of the Doppler observation;

[0069] A is the design matrix corresponding to the receiver state parameters to be estimated;

[0070] X v represents the parameter vector to be estimated, X v =(x v ,y v ,z v ,ct r ) T , x v ,y v ,z v Respectively represent the movement speed of the receiver; ct r Indicates the receiver clock error.

[0071] ε d represents the Doppler observation noise;

[0072] Q d represents the Doppler variance and covariance matrices.

[0073] Then, based on the adjustment algorithm, the Doppler observation method equation shown in formula (9) is solved to obtain the estimated parameter X v The optimal valuation vector As shown in formula (10):

[0074]

[0075] The parameter to be estimated X v The optimal valuation vector After substituting into formula (9), the posterior residual matrix v of the Doppler velocity observation value can be calculated: d As shown in formula (11):

[0076]

[0077] In formula (11), v d Represents the posterior residual matrix of the Doppler velocity observations.

[0078] Assume there are ns satellites in total, namely s1, s2, ..., s n , then the i-th satellite s i The posterior residual v of the Doppler velocity observation di It can be expressed as shown in formula (12):

[0079]

[0080] In formula (12): v di represents the i-th satellite s i The posterior residual of the Doppler velocity observation; di A represents the Doppler OMC value of the i-th satellite; i represents the i-th row of the design matrix A; represents the optimal estimate vector of the velocity parameter;

[0081] Next, the posterior residual v based on the Doppler velocity observations of each satellite di and the corresponding variance σ di , and obtain the standardized residual and Doppler velocity observation value of any i-th satellite The variance of the velocity of the i-th satellite is di According to the variance-covariance matrix Q d The diagonal elements are obtained.

[0082] Finally, the posterior residual v based on the Doppler velocity observations of each satellite di , the residual test threshold τ of the Doppler velocity observation obtained in step 1 rd , combined with the detection results of the satellites with abnormal pseudo-range positioning observation values ​​in step 2, a residual test is performed to eliminate the satellites with abnormal Doppler velocity observation values ​​among the satellites corresponding to the receiver. The process is as follows:

[0083] (3.1) Initialize the satellite list dop_outMark with Doppler observation value anomalies and the satellite list sat_outlier_doppler with Doppler observation value gross errors, both of which are empty.

[0084] Suppose the satellite with the largest residual in the epoch is numbered sat_max, and its corresponding maximum residual is vmax , initialize sat_max = 0; its corresponding residual v max =0.

[0085] According to the receiver type and the corresponding number of satellites pre-configured by the user, the Doppler residual error test adaptive threshold τ is obtained based on step 1 rd .

[0086] (3.2) Perform residual test on any i-th satellite in turn, the process is as follows:

[0087] Determine the standardized residual of the Doppler velocity observation of the i-th satellite if Then add its satellite number to the dop_outMark exception list.

[0088] Further determine the i-th satellite s i Does it exist in the pseudorange gross error satellite list sat_outlier_code? If it does, the satellite is directly marked as abnormal and stored in the Doppler observation gross error satellite list sat_outlier_doppler; if it does not exist, its residual is equal to the maximum residual record value v max Contrast, if greater than v max , then update sat_max with its satellite number and update v with its residual value max .

[0089] Repeat steps (3.1)-(3.2) for all satellites in each epoch, and then perform the subsequent steps;

[0090] (3.3) If the maximum value of the Doppler velocity estimation residual |v max |>10*τ rd , it is considered that the satellite observation value has a significant anomaly. At this time, the satellite numbered sat_max is marked as a gross error, and the satellite numbered sat_max is added to the Doppler observation gross error satellite list sat_outlier_doppler.

[0091] If |v max |<10*τ rd , all satellites in the Doppler observation value abnormal satellite list dop_outMark are directly added to the Doppler observation value gross error satellite list sat_outlier_doppler.

[0092] (3.4) Eliminate the satellites in the satellite list sat_outlier_doppler that has been added to the Doppler observation gross error in steps (3.2) and (3.3), and then perform the least squares solution according to formula (9). Repeat steps (3.1)-(3.4) until all standardized residuals and satellites are eliminated, then stop the iteration and retain the remaining satellites.

[0093] After eliminating all residuals, the optimal estimate of the Doppler velocity vector can be obtained using formula (11):

[0094] Therefore, this embodiment sets the residual error test threshold τ of the pseudorange positioning observation value rc , the residual test threshold τ of Doppler velocity observation rd , satellites with abnormal pseudorange positioning observations and abnormal Doppler velocity observations are eliminated, thereby calculating the optimal estimate of the pseudorange positioning observation and the corresponding variance covariance matrix, and then obtaining the optimal estimate of the Doppler velocity observation and the corresponding variance covariance matrix.

[0095] Step 4: Update the Kalman filter state vector of the GNSS.

[0096] In this embodiment, the calculation results of steps 1-3 can be used to update the Kalman filter state vector in the GNSS.

[0097] In high-precision differential positioning (RTK) or precise point positioning (PPP) solutions, the extended Kalman filter (EKF) is generally used for parameter estimation. k-1 The state vector of the parameters to be estimated under and its covariance matrix P k-1 , the state update matrix can be expressed as shown in formulas (13) and (14):

[0098]

[0099] in: and P k-1 are the posterior state vector and covariance matrix at epoch k respectively. If k=0, they are the initial state and covariance matrix.

[0100] and is the prior state vector and covariance matrix at time k. The state of this epoch is obtained from the forecast of the previous epoch, and the observation update has been performed before obtaining new observations.

[0101] and are the state transition matrix and the system noise covariance matrix respectively.

[0102] is the state transfer matrix The transpose of .

[0103] When new observations are available, the state quantity is updated based on the new observations, which can be expressed as shown in formulas (15), (16), and (17):

[0104]

[0105] Where: K k is the gain matrix.

[0106] and is the posterior state vector and covariance matrix updated using the new information observation at time k.

[0107] y k is the residual vector of the observed value before the test, H k is the measurement matrix, R k is the measurement noise matrix.

[0108] I is the identity matrix.

[0109] At any epoch k, the correction value of the user position and receiver clock error obtained based on the pseudorange positioning solution in step 2 is The corresponding variance-covariance-covariance information is Q p When performing Kalman filtering, the Kalman state matrix is ​​initialized at epoch k = 0, as shown in formula (18):

[0110]

[0111] Where: X p,0 and Q p,0 They represent the position and variance-covariance matrix obtained by the initial epoch pseudorange positioning solution respectively; and P0 represent the initial state vector of the Kalman filter and its corresponding covariance matrix respectively.

[0112] Based on the pseudorange and carrier phase observations of this epoch, the state of the Kalman filter can be updated according to formulas (15)-(17) to calculate the state vector of the estimated parameter of this epoch:

[0113] For dynamic positioning, when epoch k>=1, the position parameters in the Kalman filter state vector and their corresponding variance-covariance information are generally initialized using the pseudorange positioning results, as shown in formula (19):

[0114]

[0115] Where: X p,k and Q p,k They represent the position and variance-covariance matrix obtained by the pseudorange positioning solution of the kth epoch in step 2 respectively; and P k,0 They represent the initial values ​​of the Kalman filter state vector and covariance matrix of the kth epoch respectively.

[0116] The above is the Kalman filter state vector update process based on pseudo-range positioning, but the pseudo-range positioning provides The position accuracy is generally in the meter level, and the error under complex observation conditions may reach tens of meters or even hundreds of meters. In order to obtain higher positioning accuracy, this embodiment proposes to obtain a state matrix based on the above pseudo-range positioning and using the optimal estimate of the Doppler velocity observation as a constraint, thereby updating the Kalman filter state vector of the GNSS. The specific process is as follows:

[0117] (4.1) First, the system transfer matrix for state update at epochs k-1 and k is formed using the velocity constraint and the system noise covariance matrix As shown in formula (20):

[0118]

[0119] Where: I 3×3 is a unit vector matrix of size 3; τ r is the sampling interval between epochs k and k+1; I (2m-2)×(2m-2) is the state transition unit matrix of the double difference ambiguity of size 2m-2; diag(…) represents the diagonal matrix identifier; Q v is the covariance matrix of velocity; (σ ve ,σ vn ,σ vu ) is the variance-covariance matrix information of the velocity vector estimated according to step (3); E r represents the identity matrix; τ r =t k+1 -t k represents the epoch time t k and t k+1 The sampling interval between .

[0120] If epoch k>=1, the state vector calculated after Kalman updating of equations (15)-(17) at the last epoch k-1 is The corresponding variance-covariance matrix is Then, the state vector calculated according to the velocity constraint at the kth epoch can be expressed as shown in formulas (21) and (22):

[0121]

[0122] in: represents the state vector of the receiver after filtering at k-1 epochs; represents the Doppler observation value vector calculated according to step 3; N represents other ambiguity parameter information to be estimated in precise positioning;

[0123] and They represent the receiver state vector and its variance covariance vector at epoch k obtained after updating based on the velocity constraint respectively;

[0124] Different from formula (19) that performs Kalman state initialization based on pseudorange positioning, this embodiment performs Kalman state initialization based on the following formula (23):

[0125]

[0126] in: and They represent the position and variance-covariance matrix obtained by updating the velocity constraint at the kth epoch, respectively; and P k,0 They represent the initial values ​​of the Kalman filter state vector and covariance matrix of the kth epoch respectively.

[0127] In this embodiment, the position accuracy predicted based on Doppler velocity constraints Better than the accuracy of direct solution of pseudo-range positioning p,k , which can improve the filtering performance of precision positioning. However, if the Kalman filter is abnormal due to the influence of undetected abnormalities during the high-precision positioning solution, the filtering solution will be There are large gross errors. If they cannot be detected in time, the filtering update process based on formula (23) will diverge. Based on this, this embodiment further proposes to perform a consistency check, as shown in the following formula:

[0128]

[0129] Where: d refers to the three-dimensional error between the two; |·| means calculating the Euclidean geometric distance between two vectors; represents the position of the Doppler velocity constraint prediction calculated according to formula (21); X p,k represents the position directly solved by pseudorange positioning according to formula (5).

[0130] Compare d with the set error 100. If d exceeds 100, it is considered that the predicted position and the single point positioning result have poor consistency and there is a gross error. Reinitialize the initial state vector of the Kalman filter. And the corresponding variance covariance information P0.

[0131] The preferred embodiments of the present invention are described in detail above. The embodiments described in the present invention are merely descriptions of the preferred embodiments of the present invention, and are not intended to limit the concept and scope of the present invention. The various specific technical features described in the above specific embodiments can be combined in any suitable manner without contradiction, and such combinations, as long as they do not violate the concept of the present invention, should also be regarded as the contents disclosed in the present disclosure. In order to avoid unnecessary repetition, the present invention will not further describe various possible combinations.

[0132] The present invention is not limited to the specific details of the above-mentioned embodiments. Within the scope of the technical concept of the present invention and without departing from the design concept of the present invention, various modifications and improvements made to the technical solution of the present invention by technical personnel in this field should fall within the protection scope of the present invention. The technical contents for which protection is sought in the present invention have been fully recorded in the claims.

Claims

1. A positioning processing method for observation data under complex GNSS conditions, characterized in that: The following steps are involved: Step 1: setting a residual check threshold of a pseudorange positioning observation value and a residual check threshold of a Doppler velocity observation value for a satellite receiver; Step 2: Based on the residual test threshold of the pseudorange positioning observation value obtained in step 1, after finding and eliminating the satellites with abnormal pseudorange positioning observation values ​​among the satellites corresponding to the receiver, the optimal estimate of the pseudorange positioning observation and the corresponding variance covariance matrix are calculated according to the pseudorange positioning observation values ​​of the remaining satellites; Step 3: Based on the residual test threshold of the Doppler velocity observation value obtained in step 1 and the detection result of the satellite with abnormal pseudorange positioning observation value in step 2, find and eliminate the satellite with abnormal Doppler velocity observation value among the satellites corresponding to the receiver, and then calculate the optimal estimate of the Doppler velocity observation and the corresponding variance covariance matrix according to the Doppler velocity observation values ​​of the remaining satellites; Step 4: Based on the optimal estimate of the pseudorange positioning observation obtained in step 2, and using the optimal estimate of the Doppler velocity observation obtained in step 3 as a constraint, a state matrix is ​​obtained, thereby updating the Kalman filter state vector of the GNSS.

2. The positioning processing method for observation data under complex GNSS conditions according to claim 1, characterized in that: In the residual check threshold of the pseudorange positioning observation value and the residual check threshold of the Doppler velocity observation value described in step 1, each residual check threshold is composed of an adaptive factor and an empirical threshold constant, wherein the adaptive factor is determined according to the number of satellites in two adjacent epochs of the receiver, and the empirical threshold constant is set based on experience according to the type of receiver or the observation quality.

3. The positioning processing method for observation data under complex GNSS conditions according to claim 1 is characterized in that: In step 2, firstly, a mathematical model for pseudo-range single-point positioning of multiple satellites is established, and the mathematical model of pseudo-range single-point positioning is corrected for spherical delay, ionospheric delay and satellite clock error to obtain a normal equation matrix; then, the normal equation matrix is ​​solved to obtain the posterior residual of the pseudo-range positioning observation value of each satellite; then, based on the posterior residual of the pseudo-range positioning observation value of each satellite and the corresponding covariance information, the standardized residual of the pseudo-range positioning observation value of each satellite is obtained; finally, according to the standardized residual of the pseudo-range positioning observation value of each satellite and the residual check threshold of the pseudo-range positioning observation value obtained in step 1, the satellites with abnormal pseudo-range positioning observation values ​​among the satellites corresponding to the receiver are eliminated.

4. The positioning processing method for observation data under complex GNSS conditions according to claim 1, characterized in that: In step 3, firstly, a Doppler observation method equation is established for the receiver to perform Doppler velocity observations on multiple satellites; then, the Doppler observation method equation is solved to obtain the post-test residual of the Doppler velocity observation value of each satellite; then, based on the post-test residual of the Doppler velocity observation value of each satellite, the residual test threshold of the Doppler velocity observation value obtained in step 1 is used, combined with the detection result of the satellite with abnormal pseudorange positioning observation value in step 2, to eliminate the satellites with abnormal Doppler velocity observation values ​​among the satellites corresponding to the receiver.

5. The positioning processing method for observation data under complex GNSS conditions according to claim 1, characterized in that: In step 4, a consistency check is also performed based on the Euclidean geometric distance between the pseudorange positioning observation result of step 2 and the Doppler velocity observation result of step 3.

Images