UAV PPK Post-Processing Method and Device Incorporating Fusion Filtering and Ambiguity Processing

Abstract

The present disclosure provides a method and device for post-processing PPK of an unmanned aerial vehicle that integrates filtering and ambiguity processing. The positioning data is filtered multiple times, and the results of multiple filterings are integrated in subsequent processing steps. In the integration process, the present disclosure not only integrates the result types, but also integrates the numerical values. Compared with the prior art that only integrates the result types, the accuracy and anti-interference ability of the integration result are effectively improved. At the same time, the present disclosure performs composite fixation in the ambiguity fixation step, and the fixation rate is significantly better than the existing solutions.

Description

Technical Field

[0001] The present disclosure relates to the field of data processing, and particularly to the field of UAV positioning. A UAV PPK post-processing method and device that integrate filtering and ambiguity processing are disclosed. Background Art

[0002] The UAV PPK positioning technology is an advanced positioning method that allows the UAV to obtain high-precision position information during flight. The PPK technology, namely Post-Processing Kinematic, synchronously receives GPS satellite signals through a reference station and at least one rover (UAV), and then processes the data in a computer to determine the precise three-dimensional coordinates of the rover.

[0003] However, the current solution only fuses the results of two-way filtering and cannot truly exert the advantages of post-processing. At the same time, the ambiguity fixing in the current solution is only a simple partial ambiguity fixing, without a composite fixing method, and the fixing rate is low. In addition, the result fusion technology in the current solution is limited to the comparison of result types and does not examine more detailed parameters. Summary of the Invention

[0004] The present disclosure provides at least a UAV PPK post-processing method and device that integrate filtering and ambiguity processing to solve at least one of the foregoing technical problems.

[0005] According to one aspect of the present disclosure, a UAV PPK post-processing method that integrates filtering and ambiguity processing is provided, including:

[0006] Obtain the first positioning data received by the base station, and convert the first positioning data and the second positioning data into data in a predetermined data format; wherein, the second positioning data is the positioning data received by the UAV;

[0007] Perform N times of filtering processing on the first positioning data and the second positioning data converted according to the predetermined data format; wherein, N is a positive integer greater than or equal to 4;

[0008] Perform cycle slip detection processing, floating-point solution processing, integer ambiguity fixing processing, baseline solution, coordinate transformation, and adjustment processing on the first positioning data and the second positioning data after each filtering respectively, to obtain the UAV positioning result corresponding to each filtering;

[0009] For the UAV positioning results corresponding to each filtering, if the types of the N UAV positioning results are the same, then use the type of the N UAV positioning results as the initial positioning result type; if the types of two UAV positioning results are the same, then use the types of the two UAV positioning results as the initial positioning result type; if the types of the N UAV positioning results are all different, then sort the UAV positioning results according to the filtering order, and use the types of the middle two UAV positioning results as the initial positioning result type; fuse the values of the N UAV positioning results to obtain the target value of the UAV positioning result; perform up - grading and down - grading processing on the types of the N UAV positioning results to obtain the target positioning result type of the UAV positioning result.

[0010] Among them, performing up - grading and down - grading processing on the types of the N UAV positioning results includes:

[0011] Take the difference between the values of every two UAV positioning results among the types of the N UAV positioning results to obtain multiple differences. If all the obtained differences exceed the initial threshold, then downgrade the initial positioning result type;

[0012] Or, take the difference between the values of two UAV positioning results. If both differences exceed the initial threshold, then downgrade the initial positioning result type;

[0013] Among them, fusing the values of the N UAV positioning results includes:

[0014] According to the values of the N UAV positioning results, and using the acceptance residual sigma as the inverse proportional weighting, obtain the target value of the UAV positioning result;

[0015] Or, according to two UAV positioning results of the same type, and using the a posteriori residual sigma as the inverse proportional weighting, obtain the target value of the UAV positioning result.

[0016] In a possible implementation manner, the performing N - time filtering processing on the first positioning data and the second positioning data converted according to the predetermined data format includes:

[0017] Perform 4 - time filtering processing on the first positioning data and the second positioning data converted according to the predetermined data format.

[0018] In a possible implementation manner, it further includes the following steps:

[0019] The first filtering is forward Kalman filtering, and the first epoch is the initial state, and the initial state matrix and the initial variance matrix are the parameters of the initial state; where the epoch is a frame of data.

[0020] In a possible implementation manner, it further includes the following steps:

[0021] The second filtering is Kalman filtering, and the last epoch is the initial state, and the initial state matrix and the initial variance matrix are the parameters of the initial state.

[0022] In a possible implementation manner, the following steps are further included:

[0023] The third filtering is Kalman filtering. Among them, the first epoch of the third filtering is the last epoch, and the initial state matrix and the initial variance matrix are the results of the last epoch of the first filtering. The third filtering is carried out for half a course, and when it reaches the middle time point epoch, the third filtering ends.

[0024] In a possible implementation manner, the following steps are further included:

[0025] The fourth filtering is Kalman filtering. Among them, the first epoch of the fourth filtering is the first epoch, and the initial state matrix and the initial variance matrix are the results of the last epoch of the second filtering. The fourth filtering is only carried out for half a course, and when it reaches the middle time point epoch, the fourth filtering ends.

[0026] In a possible implementation manner, the integer ambiguity fixing process includes:

[0027] Screen the satellite frequencies of all full-frequency navigation satellites, select the navigation satellites of the first frequency, and perform ambiguity search and fixation on the floating-point solutions of the corresponding ambiguities obtained by floating-point solution processing;

[0028] Arrange the parameters of all full-frequency points, assign scores to the parameters according to weights respectively, arrange all the scores, perform star-kicking operations according to multiple preset ratios respectively, and perform ambiguity search and fixation on the floating-point solutions of the ambiguities obtained by floating-point solution processing corresponding to the remaining navigation satellites;

[0029] For the remaining navigation satellites, perform single-satellite elimination and cyclic single-satellite elimination until the number of eliminated navigation satellites accounts for the first preset percentage of the total number of all navigation satellites. Then perform ambiguity search and fixation on the floating-point solutions of the ambiguities obtained by floating-point solution processing corresponding to the remaining navigation satellites; for the navigation satellites that cannot be fixed among the remaining navigation satellites, perform two-satellite elimination, and so on, until the number of eliminated navigation satellites accounts for the second preset percentage of the total number of all navigation satellites. Then perform ambiguity search and fixation on the floating-point solutions of the ambiguities obtained by floating-point solution processing corresponding to the remaining navigation satellites; among them, the parameters include the elevation angle, signal-to-noise ratio, azimuth angle, and continuous tracking epoch number of the navigation satellite; the ambiguity search and fixation are used to determine the ambiguity combination; the ambiguity combination includes the integer solutions of the floating-point solutions of each navigation satellite obtained by fixation.

[0030] Based on probability theory, through statistical tests, a target integer ambiguity combination is searched from the ambiguity combinations.

[0031] According to another aspect of the present disclosure, a UAV PPK post-processing device integrating filtering and ambiguity processing is provided, including:

[0032] A data receiving and converting module, configured to obtain first positioning data received by a base station, and convert the first positioning data and second positioning data into data in a predetermined data format; wherein, the second positioning data is positioning data received by the UAV;

[0033] A filtering processing module, configured to perform N filtering processes on the first positioning data and the second positioning data converted according to the predetermined data format; wherein, N is a positive integer greater than or equal to 4;

[0034] A solution and fixed positioning module, configured to perform cycle slip detection processing, floating point solution processing, integer ambiguity fixing processing, baseline solution, coordinate transformation, and adjustment processing on the first positioning data and the second positioning data after each filtering respectively, to obtain a UAV positioning result corresponding to each filtering;

[0035] A filtering fusion module, for the UAV positioning results corresponding to each filtering, if the types of the N UAV positioning results are all the same, then use the type of the N UAV positioning results as the initial positioning result type; if the types of two UAV positioning results are the same, then use the types of the two UAV positioning results as the initial positioning result type; if the types of the N UAV positioning results are all different, then sort the UAV positioning results according to the filtering order, and use the types of the middle two UAV positioning results as the initial positioning result type; fuse the values of the N UAV positioning results to obtain a target value of the UAV positioning result; perform up and down grading processing on the types of the N UAV positioning results to obtain a target positioning result type of the UAV positioning result;

[0036] Wherein, performing up and down grading processing on the types of the N UAV positioning results includes:

[0037] Taking the difference between the values of two-by-two UAV positioning results among the types of the N UAV positioning results to obtain a plurality of differences, if all the obtained differences exceed an initial threshold, then downgrade the initial positioning result type;

[0038] Or, take the difference between the values of two UAV positioning results, if both differences exceed the initial threshold, then downgrade the initial positioning result type;

[0039] Wherein, fusing the values of the N UAV positioning results includes:

[0040] Based on the numerical values of the positioning results of N drones, and using the acceptance residual sigma as inverse proportional weighting, the target numerical value of the drone positioning results is obtained;

[0041] Alternatively, according to the positioning results of two drones of the same type, and using the a posteriori residual sigma as inverse proportional weighting, the target numerical value of the drone positioning results is obtained.

[0042] According to another aspect of the present disclosure, there is provided an electronic device, including a memory, a processor, and a computer program stored on the memory, and the processor implements the method described in any one of the above when executing the computer program.

[0043] According to another aspect of the present disclosure, there is provided a computer-readable storage medium, in which a computer program is stored, and the computer program implements the method described in any one of the above when executed by a processor.

[0044] The method and device for post-processing drone PPK with fusion filtering and ambiguity processing according to the present disclosure perform multiple filtering on the positioning data, and fuse the results of multiple filtering in subsequent processing steps. In the fusion process of the present disclosure, not only the result types are fused, but also the numerical values are fused. Compared with the prior art that only fuses the result types, the accuracy and anti-interference ability of the fusion result are effectively improved. At the same time, the present disclosure performs composite fixation in the ambiguity fixation step, and the fixation rate is significantly better than the existing solutions.

[0045] It should be understood that the content described in this part is not intended to identify the key or important features of the embodiments of the present disclosure, nor is it used to limit the scope of the present disclosure. Other features of the present disclosure will become easily understood through the following description. Description of the Drawings

[0046] The drawings are used to better understand the present solution and do not constitute a limitation to the present disclosure. Among them:

[0047] Figure 1 is a flowchart of the method for post-processing drone PPK with fusion filtering and ambiguity processing according to the present disclosure;

[0048] Figure 2A is a flowchart of the operation of cycle slip detection in the embodiment of the present disclosure;

[0049] Figure 2B is a flowchart of single-station data preprocessing in the embodiment of the present disclosure;

[0050] Figure 2C is a flowchart of the operation of the polynomial fitting algorithm in the embodiment of the present disclosure;

[0051] Figure 3 It is a schematic structural diagram of a UAV PPK post - processing device that integrates filtering and ambiguity processing according to the present disclosure;

[0052] Figure 4 It is a schematic structural diagram of an electronic device according to the present disclosure. Specific Embodiments

[0053] The following describes exemplary embodiments of the present disclosure with reference to the accompanying drawings. Various details of the embodiments of the present disclosure are included to facilitate understanding, and they should be considered merely exemplary. Therefore, those of ordinary skill in the art should recognize that various changes and modifications can be made to the embodiments described herein without departing from the scope and spirit of the present disclosure. Similarly, descriptions of well - known functions and structures are omitted for clarity and conciseness.

[0054] In the current UAV PPK post - processing solution, there is only the result fusion of two - way filtering, which cannot truly exert the advantages of post - processing. Also, in the current solution, the ambiguity fixing is only a simple partial ambiguity fixing, without a composite fixing method, and the fixing rate is low. Moreover, the result fusion technology in the current solution is limited to the comparison and fusion of result types, without examining more detailed parameters. The present disclosure provides a UAV PPK post - processing method and device that integrate filtering and ambiguity processing. The present disclosure performs multiple filterings on the positioning data of the UAV and fuses the results of multiple filterings in subsequent processing steps. In the fusion process, the present disclosure not only fuses the result types but also fuses the numerical values, effectively improving the accuracy and anti - interference ability of the fusion result compared with only fusing the result types in the prior art. At the same time, the present disclosure performs composite fixing in the ambiguity fixing step, and the fixing rate is significantly better than the existing solutions.

[0055] The technical solutions of the present disclosure will be described below through specific embodiments.

[0056] As Figure 1 shown, it is a flowchart of the UAV PPK post - processing method that integrates filtering and ambiguity processing in this embodiment. The execution subject of this embodiment is a computing device or component with data - processing capabilities. Specifically, the method of this embodiment may include the following steps:

[0057] S110. Obtain the first positioning data received by the base station and convert the first positioning data and the second positioning data into data in a predetermined data format; wherein, the second positioning data is the positioning data received by the UAV.

[0058] The UAV can be regarded as a rover station and receives the second positioning data of the UAV sent by the navigation satellite.

[0059] The base station here is a station with a fixed position, which is used to receive the first positioning data of the UAV sent by the navigation satellite.

[0060] The predetermined format here can be the RINEX format.

[0061] S120. Perform N filtering processes on the first positioning data and the second positioning data after being converted according to the predetermined data format; where N is a positive integer greater than or equal to 4.

[0062] Specifically, 4 filtering processes can be performed on the first positioning data and the second positioning data after being converted according to the predetermined data format. The first filtering is a forward Kalman filter, and the first epoch is the initial state, and the initial state matrix and the initial variance matrix are the parameters of the initial state. The second filtering is a Kalman filter, and the last epoch is the initial state, and the initial state matrix and the initial variance matrix are the parameters of the initial state. The third filtering is a Kalman filter. Among them, the first epoch of this filtering is the last epoch of the positioning data, and the initial state matrix and the initial variance matrix are the results of the last epoch of the first filtering. The third filtering is carried out for half of the process, and when it reaches the middle time point epoch, the third filtering ends. The fourth filtering is a Kalman filter. Among them, the first epoch of this filtering is the first epoch of the positioning data, and the initial state matrix and the initial variance matrix are the results of the last epoch of the second filtering. The fourth filtering only performs half of the process, and when it reaches the middle time point epoch, the fourth filtering ends.

[0063] The parameters will be different for the overall filtering accuracy and convergence time. Therefore, special processing is required for the number of filterings and the settings of the initial state matrix and the initial variance matrix. The setting method of the present disclosure can effectively improve the filtering effect and is beneficial to improving the subsequent positioning accuracy of the UAV.

[0064] S130. Perform cycle slip detection processing, floating point solution processing, integer ambiguity fixing processing, baseline solution, coordinate transformation, and adjustment processing on the first positioning data and the second positioning data after each filtering respectively, to obtain the UAV positioning results corresponding to each filtering.

[0065] The cycle slip detection processing, floating point solution processing, baseline solution, coordinate transformation, and adjustment processing are all the same as the processing steps in the prior art.

[0066] The difference lies in the above integer ambiguity fixing processing, which specifically may include the following steps:

[0067] Step 1. Screen the satellite frequency points of all full-frequency navigation satellites, select the navigation satellites of the first frequency point, and perform ambiguity search and fixing on the floating point solutions of the corresponding ambiguities obtained by the floating point solution processing.

[0068] Step 2: Arrange the parameters of all full frequency points, assign scores to the parameters according to the weights respectively, arrange all the scores, perform star kicking operations according to multiple preset ratios respectively, and perform ambiguity search and fixing on the floating point solutions of the ambiguities obtained by corresponding floating point solution processing of the remaining navigation satellites.

[0069] Step 3: For the remaining navigation satellites, perform single satellite elimination and cyclic single satellite elimination until the number of eliminated navigation satellites accounts for the first preset percentage of the total number of all navigation satellites. Then, perform ambiguity search and fixing on the floating point solutions of the ambiguities obtained by corresponding floating point solution processing of the remaining navigation satellites; for the navigation satellites that cannot be fixed among the remaining navigation satellites, perform two satellite elimination, and so on, until the number of eliminated navigation satellites accounts for the second preset percentage of the total number of all navigation satellites. Then, perform ambiguity search and fixing on the floating point solutions of the ambiguities obtained by corresponding floating point solution processing of the remaining navigation satellites; wherein, the parameters include the elevation angle of the navigation satellite, signal to noise ratio, azimuth angle, and number of consecutive tracking epochs; the ambiguity search and fixing is used to determine the ambiguity combination; the ambiguity combination includes the integer solutions of the floating point solutions of each navigation satellite obtained by fixing.

[0070] The above second preset percentage can be set according to the actual situation, for example, set to 25%. The first preset percentage is higher than the second preset percentage.

[0071] Step 4: Based on probability theory, through statistical tests, search for the target integer ambiguity combination from the ambiguity combinations.

[0072] S140: Perform fusion processing on the UAV positioning results corresponding to each filtering: For the UAV positioning results corresponding to each filtering, if the types of the N UAV positioning results are all the same, then use the type of the N UAV positioning results as the initial positioning result type; if the types of two UAV positioning results are the same, then use the types of the two UAV positioning results as the initial positioning result type; if the types of the N UAV positioning results are all different, then sort the UAV positioning results according to the filtering order and use the type of the middle UAV positioning result as the initial positioning result type; fuse the values of the N UAV positioning results to obtain the target value of the UAV positioning result; perform up - grading and down - grading processing on the types of the N UAV positioning results to obtain the target positioning result type of the UAV positioning result;

[0073] Among them, performing up - grading and down - grading processing on the types of the N UAV positioning results includes:

[0074] The numerical values of pairwise UAV positioning results among the types of N UAV positioning results are mutually subtracted to obtain a plurality of differences. If all the obtained differences exceed the initial threshold, the initial positioning result type is downgraded.

[0075] Alternatively, the numerical values of two UAV positioning results are mutually subtracted. If both differences exceed the initial threshold, the initial positioning result type is downgraded.

[0076] The result types from high to low are fixed solution, float solution, and single point solution, respectively.

[0077] Among them, the fusion of the numerical values of N UAV positioning results includes:

[0078] Based on the numerical values of N UAV positioning results and using the acceptance residual sigma as the inverse proportional weight, the target numerical value of the UAV positioning result is obtained.

[0079] Alternatively, according to the UAV positioning results of two identical types and using the a posteriori residual sigma as the inverse proportional weight, the target numerical value of the UAV positioning result is obtained.

[0080] After the receiver captures the navigation satellite signal, as long as the tracking does not interrupt (loss of lock), the receiver will automatically give the change in the integer number of carrier phase cycles during the tracking. However, in the actual process, due to the temporary blockage of the navigation satellite signal or the influence of external interference factors, the satellite signal tracking is often temporarily interrupted, resulting in the cycle slip phenomenon. When the cycle slip occurs, it will seriously reduce the carrier phase ranging accuracy, thus causing the RTK algorithm to lose the centimeter-level measurement accuracy. This step performs processing such as cycle slip identification and repair, improving the data accuracy and being beneficial to improving the subsequent UAV positioning accuracy.

[0081] Among them, the ambiguity search and fixation of the navigation satellite can be realized by the following steps:

[0082] The WL combination is used to perform ambiguity search and fixation on the navigation satellite; among them, the WL combination is the wide lane combination, with a wavelength of 0.86 m and the ambiguity being an integer.

[0083] The above steps can be implemented using specific modules. Specifically, the conversion of the data format in step S110 can be implemented using a data conversion module, which is mainly used to convert the raw data of the GNSS board into data in RINEX format for subsequent algorithms to read and process the data. The filtering in step S120 can be completed using a filtering processing module. The cycle slip detection in step S130 can be completed using a data preprocessing module. In the floating-point solution processing in step S130, it is necessary to determine the common-view satellites, and this step can be completed using a common-view satellite processing module. In the floating-point solution processing in step S130, it is necessary to determine the pseudorange observation equation, and this step can be completed using a double-difference pseudorange processing module. In the floating-point solution processing in step S130, it is necessary to determine the single-difference pseudorange observation equation, and this step can be completed using a double-difference carrier phase processing module. In the floating-point solution processing in step S130, it is necessary to determine the double-difference pseudorange equation, and combine the double-difference pseudorange equation and the double-difference equation to perform floating-point solution calculation, and this step can be completed using a floating-point solution calculation processing module. The integer ambiguity fixing processing in step S130 can be completed by an ambiguity search and fixing module. The verification of the fixed solution in step S130 can be completed by an ambiguity fixed solution verification module. Step S140 can be completed by a filtering result fusion module.

[0084] The data preprocessing module can be used to complete the following processing:

[0085] For the cycle slip problem, the present disclosure adopts a method combining single-station undifferenced data preprocessing (TurboEdit algorithm), differential data processing (polynomial fitting algorithm), and a posteriori residual analysis. The cycle slips are preliminarily detected separately in the undifferenced and differential observation data processing stages, and the detected cycle slips are marked. After the differential carrier phase solution, a posteriori residual analysis is used for judgment to determine whether there are undetected cycle slips. When the ambiguity has been fixed and the residual exceeds the limit, the satellites with possible cycle slips are excluded.

[0086] According to the signal frequencies involved in the solution, the carrier phase differential algorithm can be divided into two categories. One is single-frequency differential positioning, and the other is dual-frequency differential positioning. The main difference lies in the methods of carrier cycle slip detection and the strategies for searching for integer ambiguities. The original observables of the dual-frequency differential positioning algorithm are diverse, including dual-frequency carrier phases and dual-frequency pseudorange values, which can form various combinations of dual-frequency phases and pseudoranges, facilitating cycle slip detection and repair and accelerating the speed of ambiguity search.

[0087] The following gives several commonly used dual-frequency carrier linear combinations, and their characteristics are described taking the GPS system as an example. Assume that the phase observations of the L1 and L2 carriers of the GPS system are φ1 and φ2 respectively, and the general form of the linear combination is φ n,m = n·φ1 + m·φ2. m and n are the coefficients of the linear combination. The corresponding frequency f of the linear combination observation valuen,m , wavelength λ n,m , integer ambiguity N n,m and measurement noise The relationship with the L1 and L2 response values is:

[0088]

[0089] f1 and f2 respectively represent the frequencies of the two carriers, N1 and N2 respectively represent the integer ambiguities of the two carriers, c represents the speed of light, σ φ1 , σ φ2 are respectively the measurement noises of the two carriers.

[0090] Using the above method, various linear combinations can be performed on the dual-frequency observables. Different combined observables have different application characteristics. The linear combinations of the commonly used observables are shown in the following table:

[0091] GPS Dual-Frequency Carrier Phase Observation Value Linear Combinations

[0092]

[0093] The LI combination is called the ionosphere-free combination. It eliminates the first-order ionospheric effect and can significantly improve the accuracy of medium- and long-baseline solutions. For medium- and long-baselines (baseline distance greater than 10 km), the residual ionospheric delay error after single-differencing is still relatively large and is the main error source affecting the positioning accuracy. In this disclosure, this combination is used for baseline calculation. For short baselines (baseline distance within 10 km), the ionospheric delay error can be basically eliminated after single-differencing, and the observation noise is the main error. However, the observation noise of the LI combination is relatively large. Therefore, in this disclosure, single-frequency observation values are used for baseline calculation.

[0094] The LG combination is called the ionospheric residual combination. This combination is independent of the geometric distance from the receiver to the navigation satellite. At the same time, it eliminates the geometric distance from the receiver to the navigation satellite, orbit error, receiver clock error, satellite clock error, and tropospheric delay error, and is only related to the ionospheric delay, combined ambiguity, and observation value noise. Therefore, this observable should be very smooth in theory. In this disclosure, the LG combination is used for the detection and repair of dual-frequency cycle slips.

[0095] The WL combination is called the wide-lane combination, with a wavelength of 0.86 m, and the ambiguity is an integer, and the ionospheric delay is not large. In this disclosure, the LW combination is used to accelerate the fixation of the ambiguity.

[0096] The LN combination is called the narrow-lane combination, with a wavelength of 0.11 m, the relative observation noise is close to the noise of L1, and the ambiguity is also an integer. The ionospheric influence in LN is equal in magnitude and opposite in sign to that in LW. Although the relative noise of LN is close to that of L1, the ratio of the ionospheric error to the wavelength in LN is too large, and it is only used for high-precision positioning of short baselines. This combination is not used in this disclosure.

[0097] In addition, by combining dual-frequency carriers and dual-frequency pseudoranges, various combination forms can be formed, such as the MW combination, the GF combination, etc., which can be used for the detection of carrier cycle slips.

[0098] In the stage of processing the observation data, a polynomial fitting algorithm is adopted for single-frequency data to detect cycle slips in the single-frequency carrier phase data of double differences; for dual-frequency data, since there are relatively more observables and various dual-frequency combinations can be adopted to improve the detection efficiency, the TurboEdit algorithm is adopted in the design. Combining the dual-frequency pseudorange data, the cycle slips of the undifferenced carrier phase data are detected, and then the polynomial fitting method is used to detect the cycle slips of the L4 and L5 combinations of single differences.

[0099] In the carrier phase positioning algorithm, after detecting the carrier cycle slips, the method of cycle slip repair can be further adopted for processing. However, there are often risks in carrier cycle slip repair. Therefore, another solution adopted in this disclosure is: once a cycle slip is detected for a certain satellite, the ambiguity of this satellite is set as a new unknown. In this case, it is not necessary to repair the magnitude of the cycle slip value, avoiding the influence of incorrect determination of the cycle slip value on subsequent calculations. At the same time, there are also some drawbacks to this method. When the data quality deteriorates and there are too many satellites with cycle slips, the number of new ambiguity unknowns also increases, which will affect the positioning result. However, considering the errors caused by incorrect cycle slip repair, they are often more unacceptable. Therefore, here, cycle slip detection and repair (marking) are performed on single-station data, while in the processing of differential data, only detection is performed, but no repair is carried out. This not only enables the algorithm to have a certain ability to repair cycle slips but also reduces the risk brought by incorrect cycle slip repair. As Figure 2A shown, it is the workflow diagram of cycle slip detection before solution.

[0100] TurboEdit algorithm

[0101] Currently, TurboEdit is generally used for cycle slip detection and repair in dual-frequency RTK positioning technology. General measurement-type dual-frequency GPS receivers usually include the following observables: CA code pseudorange C1; L1 precise code pseudorange P1; L2 precise code pseudorange P2; L1 carrier phase L2 carrier phase Doppler observations D1, D2. The workflow of cycle slip detection under different combinations will be introduced below in combination with the above observables.

[0102] (1) Melbourne-Wübbena combination

[0103] Properties and characteristics of the Melbourne-Wübbena combined observations: ① It eliminates the influence of errors such as ionosphere, troposphere, receiver and satellite clock errors, and station geometry, and is only affected by observation noise and multipath effects; ② It has a relatively long wavelength (about 86 cm) and small measurement noise; ③ The calculation result of this combined observation only contains the wide-lane ambiguity parameter. Therefore, the MW combination can well complete the detection and repair of wide-lane cycle slips, rejection of gross errors, etc. The MW combination value and its variance are:

[0104]

[0105] where L7(i) is the MW combined observation value at epoch i, L1(i) and L2(i) are the dual-frequency carrier phase observation values respectively, f1 and f2 are the dual-frequency carrier frequencies, λ7 is the MW combined carrier wavelength, and N7(i) is the MW combined carrier integer ambiguity; are the measurement variances of the MW combined observation value, L1 observation value, L2 observation value, P1 observation value, and P2 observation value at epoch i respectively.

[0106] Therefore, the integer ambiguity of the combined observation value and its variance are respectively expressed as:

[0107] N7(i) = L7(i) / λ7 (4)

[0108]

[0109] In actual calculation, a recursive method can be used to calculate the predicted value <N7> of the integer ambiguity for each epoch i and its variance value:

[0110]

[0111] For the observation data of the i-th epoch, if |N7(i) - <N7> i-1 | < 4σ i-1 , it is considered that no cycle slip occurs at epoch i, and then the mean and variance of epoch i+1 are continued to be recursively calculated.

[0112] If |N7(i) - <N7> i-1 | ≥ 4σ i-1 , then a cycle slip may occur at epoch i or it may be an outlier point. If a cycle slip occurs, by calculating the difference ΔN7 of the integer ambiguities of multiple epochs before and after the cycle slip, ΔN7 is the cycle slip value between epochs, and there is the following relationship between the cycle slip ΔN7 and the cycle slips ΔN1 and ΔN2 of the L1 and L2 observation values:

[0113] ΔN7 = ΔN1 - ΔN2 (8)

[0114] It should be noted that if the cycle slips on L1 and L2 are of equal magnitude, i.e., ΔN1 = ΔN2, then ΔN7 = 0 at this time, and the cycle slips in the observed values cannot be detected using the MW combination. Therefore, other dual-frequency combination methods need to be combined for cycle slip detection.

[0115] (2) Geometry-free GF combination

[0116] By re-combining the ionospheric residual combination L4 of dual-frequency phases and the ionospheric residual combination P4 of dual-frequency code measurements, the geometry-free combined observed value LP4 of dual-frequency phases and pseudoranges is formed, and its formula is as follows:

[0117] LP4 = L4 + P4 = (L1 - L2) + (P1 - P2) = λ2N2 - λ1N1 (9)

[0118] This combination eliminates the influences of the ionosphere, troposphere, receiver and satellite clocks, and the satellite geometry at the station, etc. The influence result is the difference in ambiguities between L1 and L2. Since the integer ambiguities remain unchanged in the case of no cycle slips. Therefore, this combined observed value is suitable for data processing work such as gross error rejection, cycle slip detection and repair, etc., and when the cycle slips on L1 and L2 are of equal magnitude, i.e., ΔN1 = ΔN2, the cycle slips can also be detected using this combination.

[0119] After rejecting the gross error observed values and based on the multi-epoch observed data before and after the cycle slip, the cycle slip value ΔN4 of L4 at the cycle slip can be calculated. The relationship between ΔN4 and the cycle slips ΔN1 and ΔN2 of the L1 and L2 observed values is as follows:

[0120] ΔN4 = λ2ΔN2 - λ1ΔN1 (10)

[0121] Based on ΔN7 and ΔN4, ΔN1 and ΔN2 can be obtained as follows:

[0122]

[0123] (3) Ionosphere-free IF combination

[0124] The ionosphere-free combination L3 of dual-frequency phases and the ionosphere-free combination P3 of dual-frequency code pseudoranges are re-combined to form the observed value LP3:

[0125]

[0126] This combined observed value eliminates the influences of geometric distance, troposphere, ionosphere, etc., and only contains the influence of noise. Its disadvantage is that the noise is significantly amplified (3 times the noise of the P1 code), but the gross errors caused by the system errors of the receiver itself can be detected using LP3, and the relatively poor observed values that were not deleted in the MW combination can be rejected.

[0127] (4) Based on the above three traditional algorithms for combining carrier and pseudorange, a simplified TurboEdit algorithm is adopted in the present disclosure, that is, the MW combination and the simplified ionospheric residual combination LG combination are used for cycle slip detection respectively. The LG combination method is given in Section 3.2.1, that is:

[0128] L4 = L1 - L2 = λ2N2 - λ1N1 + ΔION L4 (13)

[0129] where ΔION L4 is the dual-frequency ionospheric residual, and its value is the difference between the dual-frequency ionospheric errors ΔION L1 and ΔION L2 , that is:

[0130]

[0131] The LG combination does not utilize the observables of P1 and P2, and its result no longer depends on the measurement accuracy of the P code, improving the measurement accuracy from the decimeter level to the millimeter level. However, at the same time, the ionospheric residual value is introduced, but the dual-frequency ionospheric residual value changes very little between adjacent epochs. The following is the estimation of the change rate of the ionospheric residual value:

[0132] The change rate of the ionospheric residual between adjacent epochs Δε I and the change rate of the single-frequency L1 ionospheric error Δε IL1 are related as:

[0133]

[0134] The period of the MEO satellite is 11h56min. During this period, the satellite rotates 360° around the earth. Then, the change amount of the observed elevation angle of the GPS satellite within 1 second is about 0.5′ (the corresponding radian value is ). Taking GPS L1 as an example, when the satellite moves Δφ, the corresponding change rate of the ionospheric delay Δε IL1 is approximately:

[0135]

[0136] Since Δφ is a very small quantity, it can be obtained that:

[0137]

[0138] The value of the single-frequency L1 ionospheric error ΔION L1 is affected by factors such as the satellite elevation angle and the electron concentration, and the error value is approximately 5m - 100m; α is the satellite elevation angle, and the elevation angle is greater than 15°. For the extreme case, the ionospheric error is 100m and the elevation angle is 15°, it can be obtained that:

[0139]

[0140] Not I ≈0.7Δε IL1 = 0.0384 m (19)

[0141] It can be seen that the change rate of the ionospheric residual of the ionospheric residual of the LG combined observation value within 1 second does not exceed 4 cm. Therefore, the error of the LG combined observation value is also at the centimeter level. Compared with the GF combination, it is easier to detect cycle slips. In this disclosure, the threshold for cycle slip detection of the LG combination is set to 5 cm. The LG combination value can detect most cycle slips, except for some cycle slips that occur on L1 and L2 and satisfy the special relationship

[0142] |λ1ΔN1 - λ2ΔN2| < 5 cm

[0143] When encountering cycle slips that cannot be detected by the GF combination, through the cooperation of the MW combination, effective detection of single-site undifferenced data can be basically achieved

[0144] (5) Through the above combined observations, cycle slips on the L1 and L2 carriers and their magnitudes ΔN1 and ΔN2 can be detected to a great extent. After cycle slip detection, the processing of carrier observations at cycle slip points mainly includes two methods: the first is to correct the carrier phase observations after the cycle slip occurs based on the cycle slips ΔN1 and ΔN2; the second is to reset the ambiguity parameters at each cycle slip point (currently, many GNSS post-processing software adopt this method, such as Bernese 5.0 software).

[0145] In this project, when the receiver observations have dual-frequency P1, P2, L1, and L2, a preprocessing of single-site observation data is performed based on the simplified TurboEdit algorithm. When a cycle slip is detected, the cycle slip values on L1 and L2 can be solved according to the corresponding formula, the carrier phase observations are corrected, and the corresponding satellites are marked

[0146] After preprocessing the data through the simplified TurboEdit, the observation data entering the main RTK algorithm is relatively "clean" and carries the corresponding "cycle slip repair mark". At this time, the main RTK algorithm can process these "marked" satellites according to the corresponding specific requirements, and adopt a polynomial fitting method based on single-difference observations to detect, repair, and mark cycle slips, thereby reconfirming the cycle slip repair of the TurboEdit algorithm. For single-frequency receivers, since there are no corresponding P2, L2, etc. observations, the single-site preprocessing step will be skipped

[0147] Taking GPS processing as an example, the single-site data preprocessing process based on the TurboEdit algorithm is asFigure 2B as shown

[0148] Polynomial fitting algorithm

[0149] Based on the observation data of a single station with dual frequencies (including P1 and P2), simplified TurboEdit data preprocessing can be carried out, and relatively "clean" carrier observation data can be basically obtained. However, due to the limitations of P-code observation noise and the simplified TurboEdit algorithm, the effect of preprocessing cannot be fully guaranteed. Therefore, after inter-station data differencing, a cycle slip detection and marking module based on the polynomial fitting algorithm is designed and implemented, and ambiguity parameters are newly added at the cycle slips.

[0150] The principle and method of polynomial fitting are as follows:

[0151] Substitute m cycle-slip-free differential observation values into the following formula:

[0152]

[0153] In the formula, n is the polynomial fitting order; t0 is the time reference; t i is the time variation; is the fitted carrier phase differential observation value corresponding to the moment of t i where i = 1, 2, 3,..., m (m > n + 1).

[0154] The first step: Use the least squares adjustment principle

[0155]

[0156] to obtain the polynomial coefficients a0, a1, a2, a3 in the formula, and calculate the mean square error according to the residuals v i after fitting

[0157]

[0158] The second step: Use the obtained polynomial coefficients to extrapolate the double-difference carrier phase value of the next epoch and compare it with the actual observation value When the difference between the two is less than 3δ, it is considered that there is no cycle slip in this double-difference observation value. Remove the earliest observation value, add this actual double-difference observation value, and then return to the first step to continue polynomial fitting and extrapolate the next epoch.

[0159] When the difference between the extrapolated value and the actual double-difference observation value is greater than 3δ, it is considered that the actual double-difference observation value contains a cycle slip. At this time, the integer week of the extrapolated value should be used to replace the integer week number of the observation value with cycle slip error. Then remove the earliest observation value and add the corrected Return to the first step to continue calculating the polynomial coefficients and extrapolate the difference values for the next epoch. Generally, the change in the fitting value above the third order is very small and has limited contribution to the fitting accuracy. Therefore, the current algorithm uses a second-order polynomial fitting based on the observations of five epochs for cycle slip detection. As Figure 2C shown, it is the flowchart of the polynomial fitting algorithm.

[0160] When the observed data is single-frequency, the receiver clock error cannot be eliminated by the single difference between stations. Therefore, the double difference method is adopted to further eliminate this error and improve the accuracy of cycle slip detection. When the observed data is dual-frequency, through the selected dual-frequency combination The single difference can eliminate various main errors including the receiver clock error. Therefore, cycle slip detection by polynomial fitting can be performed at the single difference stage. At this stage, a strategy of only detecting without repairing is adopted to "mark" the differential observation data that may have cycle slips, adding new ambiguity parameters. Subsequently, for the double difference ambiguities corresponding to the "marked" satellites, they will be reset to the parameters to be estimated and the ambiguity parameters will be solved again.

[0161] The double difference pseudorange processing module is used to complete the following processing:

[0162] Pseudorange difference uses the pseudorange observation as the basic input. Assume that the reference station and the rover station synchronously observe a group of navigation satellites, and obtain the pseudorange observations of n common-view satellites synchronously observed by the reference station (Station A) and the rover station (Station B) and i = 1, 2, 3......n. Select the navigation satellite r with the highest elevation angle among these n navigation satellites as the reference star to form the pseudorange observation equation for any navigation satellite j (j = 1, 2, 3......n and j ≠ r) and the reference star r at any t i moment as follows:

[0163]

[0164] The meanings of the parameters in the above formula are as follows:

[0165] c: speed of light (m / s);

[0166] λ: carrier wavelength of the satellite navigation signal (m);

[0167] f: carrier frequency of the satellite navigation signal (Hz);

[0168] t i The pseudorange observation value of the navigation satellite r observed by the reference station at the moment (m);

[0169] t i The pseudorange observation value of the navigation satellite j observed by the reference station at the moment (m);

[0170] t i The geometric distance (m) between the time reference station and navigation satellite r;

[0171] t i The geometric distance (m) between the time reference station and navigation satellite j;

[0172] δt A (t i ): t i The receiver clock error of the time reference station (s);

[0173] δt r (t i ): t i The satellite clock error of navigation satellite r at time t (s);

[0174] δt j (t i ): t i The satellite clock error of navigation satellite j at time t (s);

[0175] t i The earth rotation error (m) between the time reference station and navigation satellite r at time t;

[0176] t i The earth rotation error (m) between the time reference station and navigation satellite j at time t;

[0177] t i The ionospheric delay error (m) between the time reference station and navigation satellite r at time t;

[0178] t i The ionospheric delay error (m) between the time reference station and navigation satellite j at time t;

[0179] t i The tropospheric delay error (m) between the time reference station and navigation satellite r at time t;

[0180] t i The tropospheric delay error (m) between the time reference station and navigation satellite j at time t;

[0181] t i The pseudorange measurement thermal noise (m) of the time reference station receiver with respect to navigation satellite r at time t;

[0182] t iPseudorange measurement thermal noise (m) of the reference station receiver with respect to navigation satellite j;

[0183] t i Pseudorange observation value (m) of navigation satellite r observed by the rover station at time t;

[0184] t i Pseudorange observation value (m) of navigation satellite j observed by the rover station at time t;

[0185] t i Geometric distance (m) between the rover station and navigation satellite r at time t;

[0186] t i Geometric distance (m) between the rover station and navigation satellite j at time t;

[0187] δt B (t i ) : t i Rover station receiver clock error (s);

[0188] t i Earth rotation error (m) between the rover station and navigation satellite r at time t;

[0189] t i Earth rotation error (m) between the rover station and navigation satellite j at time t;

[0190] t i Ionospheric delay error (m) between the rover station and navigation satellite r at time t;

[0191] t i Ionospheric delay error (m) between the rover station and navigation satellite j at time t;

[0192] t i Tropospheric delay error (m) between the rover station and navigation satellite r at time t;

[0193] t i Tropospheric delay error (m) between the rover station and navigation satellite j at time t;

[0194] t i Pseudorange measurement thermal noise (m) of the rover station receiver with respect to navigation satellite r;

[0195] ti Pseudorange measurement thermal noise (m) of the rover receiver at a certain moment with respect to navigation satellite j;

[0196] Linearization:

[0197] In the above formula,

[0198]

[0199] In the formula, (X r , Y r , Z r ) is the satellite position, and (X A , Y A , Z A ) is the receiver position

[0200] Let (X A0 , Y A0 , Z A0 ), (δX A (t), δY A (t), δZ A (t)) be the approximate value and correction of the observation station coordinates respectively. After Taylor expansion of equation (a) at (X A0 , Y A0 , Z A0 ), the linearized observation equation can be obtained as follows:

[0201]

[0202] In the formula,

[0203]

[0204] Then the formula in (23) can be transformed into:

[0205]

[0206] The double-difference carrier phase processing module is used to complete the following processing:

[0207] Due to the existence of various errors such as the user receiver clock error, navigation satellite clock error, earth rotation error, ionospheric delay, and tropospheric delay, the pseudorange observation equation and the carrier phase observation equation have many unknown parameters, complex solutions, and poor solution accuracy. Since the common errors such as the navigation satellite clock error, earth rotation error, ionospheric delay, and tropospheric delay between the reference station and the rover are highly correlated, the influence of the common errors can be eliminated by taking the difference between the observation equations to improve the solution accuracy. By taking the difference between the two pseudorange observation equations of the reference station and the rover corresponding to the same navigation satellite, the corresponding pseudorange single-difference observation equation can be obtained as follows:

[0208]

[0209] The meanings of the parameters in the above formula are as follows:

[0210] t i The single-difference pseudorange between the reference station and the rover at time t with respect to navigation satellite r (m);

[0211] t i The single-difference pseudorange between the reference station and the rover at time t with respect to navigation satellite j (m);

[0212] δt AB (t i ): t i The difference between the receiver clock error of the reference station and the receiver clock error of the rover at time t (s);

[0213] t i The single-difference pseudorange noise between the reference station and the rover at time t with respect to navigation satellite r (m); t i The single-difference pseudorange noise between the reference station and the rover at time t with respect to navigation satellite j (m);

[0214] The single-difference observation equation eliminates the influence of common errors such as navigation satellite clock error, ionospheric delay error, and tropospheric delay error. However, due to the existence of the user receiver clock error in the single-difference observation equation, the complexity of the solution is increased, affecting the solution accuracy. By taking the difference of the single-difference observation equation again, the influence of the user receiver clock error can be eliminated, further improving the relative position solution accuracy.

[0215] Taking the difference between the single-difference pseudorange observation equations of navigation satellite j observed by the reference station and the rover and the single-difference pseudorange observation equation of the reference satellite r respectively, n - 1 pseudorange double-difference equations are obtained as follows:

[0216]

[0217] The meanings of the parameters in the above formula are as follows:

[0218] t i The pseudorange measurement residual between the reference station and the rover at time t (m);

[0219] t i The double-difference measurement coefficient matrix of the rover with respect to navigation satellites r and j at time t;

[0220] t i The position correction vector between the rovers at time t (m);

[0221] The constant term (m) of the pseudorange measurements of the base station and the rover station with respect to the navigation satellites r and j;

[0222] The double-difference pseudorange observation equation can be simplified as:

[0223] L ρ =AX+V ρ (30)

[0224] The meanings of the parameters in the above formula are as follows:

[0225] X: position correction vector of the rover;

[0226] L ρ : Pseudorange double difference measurement constant vector between the base station and the rover;

[0227] A: relative position vector coefficient matrix between the base station and the mobile station;

[0228] V ρ : Pseudorange double difference measurement residual between the base station and the rover;

[0229] The floating point solution processing module is used to complete the following processing:

[0230] For n GPS solution satellites, the carrier phase double difference observation equation contains 3 unknown relative position parameters and n-1 unknown integer ambiguity parameters. The number of unknown parameters is greater than the number of equations. In order to solve the above carrier phase double difference observation equation, it is necessary to comprehensively utilize the two observation quantities of pseudorange and carrier. The pseudorange carrier phase joint solution process is given below.

[0231] Combining the double-difference pseudorange equation and the carrier phase double-difference equation, we can get the following set of equations:

[0232]

[0233] The meanings of the parameters in the above formula are as follows:

[0234] a: relative position correction vector between the base station and the rover;

[0235] b: integer ambiguity vector between the base station and the rover;

[0236] L ρ : Pseudorange double difference measurement constant vector between the base station and the rover;

[0237] The carrier phase double difference measurement constant vector between the base station and the rover;

[0238] A: relative position vector coefficient matrix between the base station and the mobile station;

[0239] B: integer ambiguity vector coefficient matrix between the base station and the mobile station;

[0240] V ρ : The pseudorange double-difference measurement residual between the reference station and the rover station;

[0241] The carrier-phase double-difference measurement residual between the reference station and the rover station;

[0242] P ρ : The pseudorange double-difference measurement weight coefficient matrix between the reference station and the rover station;

[0243] The carrier-phase double-difference measurement weight coefficient matrix between the reference station and the rover station;

[0244] P: The basic weight coefficient matrix, n is the number of satellites for solution, lin and col represent the rows and columns of matrix P, is the measurement noise variance;

[0245] q: The pseudorange double-difference measurement weight coefficient factor between the reference station and the rover station, whose magnitude depends on the pseudorange measurement accuracy and error level, and the typical value is 10-4 to 10-6;

[0246] The above equations can be expressed in the following matrix form:

[0247]

[0248] Adding the coefficient matrix and using the least squares method, we can obtain:

[0249]

[0250] The above equation can be simplified to:

[0251]

[0252] Let N aa = A T PA, N ab = A T PB, N ba = B T PA, N bb = B T PB, Then the above equation can be briefly written as:

[0253]

[0254] Using the least squares method, the relative position and the floating-point solution result of the integer ambiguity can be obtained:

[0255]

[0256]

[0257] The ambiguity search and fixation module can perform the following data processing:

[0258] The three-step method is adopted to fix the ambiguity: 1. Perform ambiguity search only for the first frequency point of all ambiguities; 2. Use information such as satellite elevation angle, signal-to-noise ratio, and azimuth angle for control, screen satellites, and perform partial ambiguity fixation; 3. Finally, use the method of partially excluding satellites to perform cyclic satellite kicking, and then perform ambiguity fixation.

[0259] The first step is to screen the satellite frequency points of all satellites in the whole system and all frequency points, and only select the satellites of the first frequency point for ambiguity search, which can not only reduce the time of ambiguity search and fixation, but also take into account the positioning accuracy of the fixed solution.

[0260] The second step is to screen the satellite frequency points using information such as satellite elevation angle, signal-to-noise ratio, and number of consecutive tracking epochs. First, arrange the information such as satellite elevation angle, signal-to-noise ratio, and number of consecutive tracking epochs of all satellites in the whole system and all frequency points, and assign scores to these parameters as a whole according to weights respectively. Arrange the overall scores, and perform overall satellite kicking according to three ratios of 90%, 85%, and 80% respectively, and search and fix the remaining satellites.

[0261] The third step is to use the method of partially excluding satellites to perform ambiguity search and fixation. In the first step, single satellites are excluded, and single satellites are cyclically excluded, and then ambiguity search and fixation are performed. If this method cannot fix, then 2 satellites are excluded, and so on, until the number of excluded satellites reaches 25% of all satellites, and no more satellite exclusion processing is performed.

[0262] The ambiguity fixed solution verification module is used to complete the following processing:

[0263] Fixed solution verification method

[0264] Ambiguity confirmation is the last link in the integer ambiguity solution. Its main function is to select the optimal one from the set of possible ambiguity combinations obtained during the ambiguity search process and judge its correctness. Generally speaking, all ambiguity solution methods select the optimal one from the set of ambiguity combinations through a certain test. Therefore, the choice of the test method becomes a crucial issue. If the conditions are too loose, it will lead to an extended time for the whole process, which is not conducive to the rapid solution of the ambiguity; if the conditions are too strict, it is very likely to exclude the correct ambiguity.

[0265] The ambiguity confirmation algorithm adopted here is based on probability theory. Through three statistical tests, it is finally confirmed that the integer ambiguity combination searched out is the correct integer ambiguity combination to be sought. The solution obtained after substituting this set of integer ambiguities into the normal equation is the correct fixed solution.

[0266] (1) Consistency test of the baseline vectors obtained from the integer solution and the initial solution.

[0267] Let the baseline vector obtained from the integer solution be X, and the baseline vector obtained from the initial solution be The corresponding cofactor matrix is If the following equation holds, then and X are consistent and compatible from the perspective of statistical tests:

[0268]

[0269] where β = ξ F (u, f, 1 - α) is the one-tailed quantile value of the Fisher distribution with a confidence level of 1 - α, degrees of freedom of f and u; u is the number of unknown parameters; f is the degrees of freedom in parameter estimation; σ0 is the mean square error of unit weight in the initial solution.

[0270] Since the ambiguity parameters are already within the corresponding confidence intervals, as long as the baseline vectors are also consistent, it means that both the integer solution and the initial solution vectors are consistent.

[0271] (2) Consistency test of the mean square error of unit weight between the integer solution and the initial solution.

[0272] Let the mean square error of unit weight of the integer solution be σ A , and the mean square error of unit weight of the initial solution be σ0. If the following equation holds, it means they are consistent from the perspective of statistical tests:

[0273]

[0274] The above test is also called the χ 2 test of the variance factor. The meanings of the symbols in the above equation are the same as above.

[0275] (3) Ratio test: Significance test between the minimum mean square error of unit weight σ and the second minimum mean square error of unit weight σ 次 in the integer solution.

[0276] The Ratio test is one of the most commonly used test methods. Its basic idea is to compare the minimum sum of squared residuals and the second minimum sum of squared residuals. Since the double-difference residuals calculated using the correct ambiguity group are significantly smaller than those of other incorrect ones, a threshold can be set according to various error factors such as the measurement error and multipath error level. If it satisfies:

[0277]

[0278] That is, if it is greater than the preset threshold value, it is considered to pass the inspection and the fixed solution calculation is correct; otherwise, it is considered that there is still a deviation in the solution. In this design, the ratio value is set to ≥3.

[0279] If any one of the above three inspections fails, it means that the search fails, and rework and retesting are required or other measures need to be taken (such as expanding the range of the confidence interval, etc.).

[0280] The filtering result fusion module is used to complete the following processing:

[0281] Since there are four filtering results to be fused in this disclosure, the arrangement of the result location types is carried out first. There will be four positioning results for each epoch, and the preliminary determination of the result types can be carried out according to the following three schemes: First, define the result types first. If the four result types are the same, it is tentatively set as this positioning result type; Second, if two are the same, temporarily use these two same types as the positioning result type; If the four result types are all different, take the middle result in order as the standard.

[0282] Fusion of the result values: First, according to the existing four results, use the acceptance residual sigma as the inverse proportional weight, and finally obtain the positioning result; Second, according to the positioning results of two same types, use the a posteriori residual sigma as the inverse proportional weight, and finally obtain the positioning result.

[0283] Upgrading or downgrading of the result types: First, calculate the mutual difference between every two of the four positioning results to obtain 12 differences. If all 12 differences exceed the initial threshold, downgrade the result; Second, calculate the mutual difference between two positioning results. If both differences exceed the initial threshold, downgrade the result; Third, do not process the result types of the data.

[0284] So far, the target positioning result type and the target value of the UAV are obtained, that is, the positioning result of the UAV is obtained.

[0285] In the present disclosure, in the PPK positioning of an unmanned aerial vehicle (UAV), after the rover, i.e., the UAV, obtains the relevant observation data of the reference station through post-processing and combines it with the local original observation data, a double-difference equation is formed between the two to eliminate most of the correlation errors. In this way, the carrier-phase double-difference ambiguity can be fixed relatively quickly to obtain a high-precision positioning result. To ensure positioning accuracy and reliability, according to the actual application requirements of RTK positioning of GNSS survey receivers. The algorithm mainly includes parsing and obtaining dual-station original observation data (this data includes GPS: L1L1L5 / BDS: B1I B2I B3I B1C B2A / QZSS: L1 L2 L5 / GAL: E1 E5aE5b / GLN: LG1 G2), and the original data includes pseudorange, carrier phase, Doppler, signal-to-noise ratio information, etc., single-station data preprocessing based on simplified TurboEdit, single-point positioning and single-point velocity measurement of the base station and the rover, pseudorange differential positioning, cycle slip detection, repair and marking based on polynomial fitting algorithm, establishment and adjustment of carrier-phase double-difference observation, solution of floating-point ambiguity and covariance matrix, LAMBDA ambiguity search and fixation, baseline solution based on fixed ambiguity, a posteriori residual estimation and analysis, and result output, etc. steps. In the PPK solution process of the present disclosure, the filtering results of four-way filtering are used and the results are fused. Different from the existing two-way filtering method, in the PPK solution of the present disclosure, the four-way result fusion in the PPK solution comprehensively considers parameters such as baseline length, ratio, nb, rms, etc. At the same time, in the PPK solution of the present disclosure, the ambiguity fixation method adds a composite fixation method to the existing traditional fixation methods, and the fixation rate is significantly better than the existing methods.

[0286] Based on the same inventive concept, the present disclosure provides a UAV PPK post-processing device for fusion filtering and ambiguity processing. The steps executed by the components of this device are the same as or similar to those of the above method, so the similar parts will not be elaborated here. As Figure 3 shown, the UAV PPK post-processing device for fusion filtering and ambiguity processing of this embodiment includes:

[0287] A data receiving and converting module 310, configured to obtain the first positioning data received by the base station and convert the first positioning data and the second positioning data into data in a predetermined data format; wherein, the second positioning data is the positioning data received by the UAV.

[0288] A filtering processing module 320, configured to perform N times of filtering processing on the first positioning data and the second positioning data converted according to the predetermined data format; wherein, N is a positive integer greater than or equal to 4.

[0289] The solution fixed positioning module 330 is used to perform cycle slip detection processing, floating point solution processing, integer ambiguity fixing processing, baseline solution, coordinate transformation, and adjustment processing on the first positioning data and the second positioning data after each filtering respectively, so as to obtain the UAV positioning result corresponding to each filtering;

[0290] The filtering and fusion module 340 is used for the UAV positioning results corresponding to each filtering. If the types of the N UAV positioning results are all the same, then use the type of the N UAV positioning results as the initial positioning result type; if the types of two UAV positioning results are the same, then use the types of the two UAV positioning results as the initial positioning result type; if the types of the N UAV positioning results are all different, then sort the UAV positioning results according to the filtering order, and use the types of the middle two UAV positioning results as the initial positioning result type; fuse the values of the N UAV positioning results to obtain the target value of the UAV positioning result; perform up - grading and down - grading processing on the types of the N UAV positioning results to obtain the target positioning result type of the UAV positioning result;

[0291] Among them, performing up - grading and down - grading processing on the types of N UAV positioning results includes:

[0292] Take the difference between the values of every two UAV positioning results among the types of the N UAV positioning results to obtain multiple differences. If all the obtained differences exceed the initial threshold, then downgrade the initial positioning result type;

[0293] Or, take the difference between the values of two UAV positioning results. If both differences exceed the initial threshold, then downgrade the initial positioning result type;

[0294] Among them, fusing the values of N UAV positioning results includes:

[0295] According to the values of the N UAV positioning results, and using the acceptance residual sigma as the inverse proportional weighting, obtain the target value of the UAV positioning result;

[0296] Or, according to two UAV positioning results of the same type, and using the a posteriori residual sigma as the inverse proportional weighting, obtain the target value of the UAV positioning result.

[0297] In some embodiments, when the filtering processing module 320 performs N - time filtering processing on the first positioning data and the second positioning data converted according to the predetermined data format, it is specifically used for:

[0298] Perform 4 - time filtering processing on the first positioning data and the second positioning data converted according to the predetermined data format.

[0299] In some embodiments, the first filtering is forward Kalman filtering, and the first epoch is the initial state, and the initial state matrix and the initial variance matrix are parameters of the initial state; wherein the epoch is a frame of data.

[0300] In some embodiments, the second filtering is Kalman filtering, and the last epoch is the initial state, and the initial state matrix and the initial variance matrix are parameters of the initial state.

[0301] In some embodiments, the third filtering is Kalman filtering, wherein the first epoch of the third filtering is the last epoch, and the initial state matrix and the initial variance matrix are the results of the last epoch of the first filtering. The third filtering is performed for half the process and ends at the epoch of the middle time point.

[0302] In some embodiments, the fourth filtering is Kalman filtering, wherein the first epoch of the fourth filtering is the first epoch, and the initial state matrix and the initial variance matrix are the results of the last epoch of the second filtering. The fourth filtering is only performed for half the process and ends at the epoch of the middle time point.

[0303] In some embodiments, when the solution fixed positioning module 330 performs integer ambiguity fixing processing, it is specifically used for:

[0304] Screen the satellite frequency points of all full-frequency navigation satellites, select the navigation satellites of the first frequency point, and perform ambiguity search and fixing on the floating-point solutions of the corresponding ambiguities obtained by floating-point solution processing;

[0305] Arrange the parameters of all full-frequency points, perform score assignment processing on the parameters according to the weights respectively, arrange all the scores, perform satellite kicking operations according to multiple preset ratios respectively, and perform ambiguity search and fixing on the floating-point solutions of the ambiguities obtained by the corresponding floating-point solution processing of the remaining navigation satellites;

[0306] For the remaining navigation satellites, single satellite rejection and cyclic single satellite rejection are performed until the number of rejected navigation satellites accounts for a first preset percentage of the total number of all navigation satellites. Then, ambiguity search and fixing are performed on the floating point solutions of the ambiguities obtained by corresponding floating point solution processing of the remaining navigation satellites. For the navigation satellites that cannot be fixed among the remaining navigation satellites, two satellites are rejected, and so on, until the number of rejected navigation satellites accounts for a second preset percentage of the total number of all navigation satellites. Then, ambiguity search and fixing are performed on the floating point solutions of the ambiguities obtained by corresponding floating point solution processing of the remaining navigation satellites. Among them, the parameters include the elevation angle, signal-to-noise ratio, azimuth angle, and continuous tracking epoch number of the navigation satellites. The ambiguity search and fixing are used to determine the ambiguity combination. The ambiguity combination includes the integer solutions of the floating point solutions of each navigation satellite obtained by fixing.

[0307] Based on probability theory, through statistical tests, the target integer ambiguity combination is searched from the ambiguity combinations.

[0308] According to an embodiment of the present disclosure, the present disclosure also provides an electronic device and a computer-readable storage medium.

[0309] Figure 4 FIG. shows a schematic block diagram of an exemplary electronic device 400 that can be used to implement the embodiments of the present disclosure. The electronic device is intended to represent various forms of digital computers, such as, laptop computers, desktop computers, workstations, personal digital assistants, servers, blade servers, mainframe computers, and other suitable computers. The electronic device can also represent various forms of mobile devices, such as, personal digital processing, cellular phones, smart phones, wearable devices, and other similar computing devices. The components shown herein, their connections and relationships, and their functions are merely examples and are not intended to limit the implementation of the present disclosure described and / or claimed herein.

[0310] As Figure 4 shown, the device 400 includes a computing unit 410, which can perform various appropriate actions and processes according to the computer program stored in the read-only memory (ROM) 420 or the computer program loaded from the storage unit 480 into the random access memory (RAM) 430. In the RAM 430, various programs and data required for the operation of the device 400 can also be stored. The computing unit 410, the ROM 420, and the RAM 430 are connected to each other through a bus 440. The input / output (I / O) interface 450 is also connected to the bus 440.

[0311] Multiple components in device 400 are connected to I / O interface 450, including: input unit 460, such as a keyboard, mouse, etc.; output unit 470, such as various types of displays, speakers, etc.; storage unit 480, such as a disk, optical disc, etc.; and communication unit 490, such as a network card, modem, wireless communication transceiver, etc. Communication unit 490 allows device 400 to exchange information / data with other devices via a computer network such as the Internet and / or various telecommunication networks.

[0312] Computing unit 410 can be various general-purpose and / or special-purpose processing components with processing and computing capabilities. Some examples of computing unit 410 include, but are not limited to, a central processing unit (CPU), a graphics processing unit (GPU), various dedicated artificial intelligence (AI) computing chips, various computing units running machine learning model algorithms, a digital signal processor (DSP), and any suitable processor, controller, microcontroller, etc. Computing unit 410 executes the various methods and processes described above. For example, in some embodiments, any of the above methods can be implemented as a computer software program tangibly embodied in a machine-readable medium, such as storage unit 480. In some embodiments, part or all of the computer program can be loaded and / or installed onto device 400 via ROM 420 and / or communication unit 490. When the computer program is loaded into RAM 430 and executed by computing unit 410, one or more steps of any of the above-described methods can be performed. Alternatively, in other embodiments, computing unit 410 can be configured to execute any of the above-described methods in any other suitable manner (e.g., by means of firmware).

[0313] The various embodiments of the systems and techniques described above can be implemented in digital electronic circuitry, integrated circuit systems, field programmable gate arrays (FPGAs), application specific integrated circuits (ASICs), application specific standard products (ASSPs), systems on a chip (SOCs), complex programmable logic devices (CPLDs), computer hardware, firmware, software, and / or combinations thereof. These various embodiments can include: implemented in one or more computer programs that can be executed and / or interpreted on a programmable system including at least one programmable processor, which can be a special or general-purpose programmable processor that receives data and instructions from a storage system, at least one input device, and at least one output device, and transmits the data and instructions to the storage system, the at least one input device, and the at least one output device.

[0314] The program code for implementing the methods of the present disclosure can be written in any combination of one or more programming languages. These program codes can be provided to a processor or controller of a general-purpose computer, a special-purpose computer, or other programmable data processing devices, such that when the program codes are executed by the processor or controller, the functions / operations specified in the flowchart and / or block diagram are implemented. The program codes can be executed entirely on the machine, partially on the machine, executed partially on the machine and partially on a remote machine as an independent software package, or executed entirely on a remote machine or server.

[0315] In the context of the present disclosure, a machine-readable medium can be a tangible medium that can contain or store a program for use by or in connection with an instruction execution system, apparatus, or device. A machine-readable medium can be a machine-readable signal medium or a machine-readable storage medium. A machine-readable medium can include, but is not limited to, electronic, magnetic, optical, electromagnetic, infrared, or semiconductor systems, apparatus, or devices, or any suitable combination of the foregoing. More specific examples of a machine-readable storage medium would include an electrical connection based on one or more wires, a portable computer disk, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), an optical fiber, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing.

[0316] In order to provide interaction with a user, the systems and techniques described herein can be implemented on a computer having: a display device (e.g., a CRT (cathode ray tube) or LCD (liquid crystal display) monitor) for displaying information to the user; and a keyboard and a pointing device (e.g., a mouse or a trackball) through which the user can provide input to the computer. Other kinds of devices can also be used to provide interaction with the user; for example, the feedback provided to the user can be any form of sensory feedback (e.g., visual feedback, auditory feedback, or tactile feedback); and the input from the user can be received in any form (including acoustic input, voice input, or tactile input).

[0317] The systems and techniques described herein can be implemented in a computing system including backend components (e.g., as a data server), or a computing system including middleware components (e.g., an application server), or a computing system including frontend components (e.g., a user computer having a graphical user interface or a web browser through which a user can interact with an implementation of the systems and techniques described herein), or a computing system including any combination of such backend components, middleware components, or frontend components. The components of the system can be interconnected to each other by digital data communication in any form or medium (e.g., a communication network). Examples of communication networks include: local area network (LAN), wide area network (WAN), and the Internet.

[0318] A computer system can include a client and a server. The client and the server are generally far from each other and typically interact through a communication network. The client-server relationship is created by computer programs running on the respective computers and having a client-server relationship with each other. The server can be a cloud server, a server of a distributed system, or a server incorporating blockchain.

[0319] It should be understood that various forms of the processes shown above can be used, steps can be reordered, added, or deleted. For example, the steps recited in this disclosure can be executed in parallel, sequentially, or in a different order, as long as the desired results of the technical solutions disclosed in this disclosure can be achieved, and this is not limited herein.

[0320] The above specific embodiments do not constitute a limitation on the protection scope of this disclosure. Those skilled in the art should understand that various modifications, combinations, sub-combinations, and substitutions can be made according to design requirements and other factors. Any modifications, equivalent substitutions, and improvements made within the spirit and principle of this disclosure shall be included within the protection scope of this disclosure.

Claims

1. A UAV PPK post-processing method integrating filtering and fuzzy processing, characterized in that: include: Acquire first positioning data received by the base station, and convert the first positioning data and second positioning data into data in a predetermined data format; wherein the second positioning data is positioning data received by the drone; Performing N filtering processes on the first positioning data and the second positioning data converted according to the predetermined data format; wherein N is a positive integer greater than or equal to 4; Perform cycle slip detection processing, floating point solution processing, integer ambiguity fixation processing, baseline solution, coordinate conversion and adjustment processing on the first positioning data and the second positioning data after each filtering, and obtain the UAV positioning result corresponding to each filtering; For the drone positioning results corresponding to each filtering, if the types of N drone positioning results are the same, the types of the N drone positioning results are used as the initial positioning result type; if the types of two drone positioning results are the same, the types of the two drone positioning results are used as the initial positioning result type; if the types of N drone positioning results are different, the drone positioning results are sorted according to the filtering order, and the types of the two middle drone positioning results are used as the initial positioning result type; the values ​​of the N drone positioning results are merged to obtain the target value of the drone positioning result; the types of the N drone positioning results are upgraded and downgraded to obtain the target positioning result type of the drone positioning result; Among them, the types of N drone positioning results are upgraded and downgraded, including: Subtract the values ​​of the positioning results of two drones in the types of N drone positioning results to obtain multiple difference values. If all the obtained difference values ​​exceed the initial threshold, the type of the initial positioning result is downgraded; Alternatively, the values ​​of the positioning results of the two drones are subtracted from each other, and if both differences exceed the initial threshold, the type of the initial positioning result is downgraded; Among them, the numerical values ​​of N UAV positioning results are fused, including: According to the values ​​of the positioning results of N drones, and using the acceptance residual sigma as the inverse proportional weighting, the target value of the drone positioning result is obtained; Alternatively, the target value of the UAV positioning result is obtained by taking the positioning results of two UAVs of the same type and using the posterior residual sigma as an inverse proportional weight.

2. The method according to claim 1, characterized in that The performing N times of filtering on the first positioning data and the second positioning data converted according to the predetermined data format comprises: The first positioning data and the second positioning data converted according to the predetermined data format are subjected to four filtering processes.

3. The method according to claim 2, characterized in that Also includes: The first filtering is a forward Kalman filter, and the first epoch is the initial state, and the initial state matrix and the initial variance matrix are parameters of the initial state; wherein the epoch is a frame of data.

4. The method according to claim 3, characterized in that Also includes: The second filtering is Kalman filtering, and the last epoch is the initial state, and the initial state matrix and initial variance matrix are the parameters of the initial state.

5. The method according to claim 4, characterized in that Also includes: The third filtering is Kalman filtering, in which the first epoch of the third filtering is the last epoch, and the initial state matrix and the initial variance matrix are the results of the last epoch of the first filtering. The third filtering is performed halfway, and the third filtering ends at the middle time point epoch.

6. The method according to claim 5, characterized in that Also includes: The fourth filtering is Kalman filtering, in which the first epoch of the fourth filtering is the first epoch, and the initial state matrix and the initial variance matrix are the results of the last epoch of the second filtering. The fourth filtering is only performed halfway, and the fourth filtering ends at the middle time point epoch.

7. The method according to claim 1, characterized in that The integer ambiguity fixing process comprises: Screening satellite frequencies of all full-frequency navigation satellites, selecting a navigation satellite of the first frequency, and performing ambiguity search and fixing on the corresponding ambiguity floating-point solutions obtained by floating-point solution processing; Arrange the parameters of all full-frequency points, assign scores to the parameters according to weights, arrange all scores, perform star-kicking operations according to multiple preset ratios, and perform ambiguity search and fixation on the floating-point solutions of the ambiguities obtained by the floating-point solution processing corresponding to the remaining navigation satellites; For the remaining navigation satellites, single satellites are eliminated, and single satellites are eliminated cyclically until the number of eliminated navigation satellites accounts for a first preset percentage of the number of all navigation satellites, and then ambiguity search and fixation are performed on the floating-point solutions of ambiguities obtained by floating-point solution processing corresponding to the remaining navigation satellites; for the navigation satellites that cannot be fixed among the remaining navigation satellites, two satellites are eliminated, and so on, until the number of eliminated navigation satellites accounts for a second preset percentage of the number of all navigation satellites, and then ambiguity search and fixation are performed on the floating-point solutions of ambiguities obtained by floating-point solution processing corresponding to the remaining navigation satellites; wherein the parameters include navigation satellite altitude angle, signal-to-noise ratio, azimuth angle, and number of continuous tracking epochs; the ambiguity search and fixation are used to determine an ambiguity combination; the ambiguity combination includes integer solutions of the floating-point solutions of each navigation satellite obtained by fixing; Based on probability theory and through statistical testing, a target integer ambiguity combination is searched out from the ambiguity combinations.

8. A UAV PPK post-processing device integrating filtering and fuzzy processing, characterized in that: include: A data receiving and converting module, used for acquiring first positioning data received by the base station, and converting the first positioning data and the second positioning data into data in a predetermined data format; wherein the second positioning data is positioning data received by the drone; A filtering processing module, used to perform N filtering processes on the first positioning data and the second positioning data converted according to the predetermined data format; wherein N is a positive integer greater than or equal to 4; The fixed positioning solution module is used to perform cycle slip detection processing, floating point solution processing, integer ambiguity fixation processing, baseline solution, coordinate conversion and adjustment processing on the first positioning data and the second positioning data after each filtering, so as to obtain the UAV positioning result corresponding to each filtering; The filtering fusion module is used for the UAV positioning results corresponding to each filtering. If the types of N UAV positioning results are the same, the types of the N UAV positioning results are used as the initial positioning result type; if the types of two UAV positioning results are the same, the types of the two UAV positioning results are used as the initial positioning result type; if the types of N UAV positioning results are different, the UAV positioning results are sorted according to the filtering order, and the types of the two middle UAV positioning results are used as the initial positioning result type; the values ​​of the N UAV positioning results are fused to obtain the target value of the UAV positioning result; the types of the N UAV positioning results are upgraded and upgraded to obtain the target positioning result type of the UAV positioning result; Among them, the types of N drone positioning results are upgraded and downgraded, including: Subtract the values ​​of the positioning results of two drones in the types of N drone positioning results to obtain multiple difference values. If all the obtained difference values ​​exceed the initial threshold, the type of the initial positioning result is downgraded; Alternatively, the values ​​of the positioning results of the two drones are subtracted from each other, and if both differences exceed the initial threshold, the type of the initial positioning result is downgraded; Among them, the numerical values ​​of N UAV positioning results are fused, including: According to the values ​​of the positioning results of N drones, and using the acceptance residual sigma as the inverse proportional weighting, the target value of the drone positioning result is obtained; Alternatively, the target value of the UAV positioning result is obtained by taking the positioning results of two UAVs of the same type and using the posterior residual sigma as an inverse proportional weight.

Images