A Data Processing Method for Deformation Monitoring of Beidou / GNSS Sluices in Complex Occlusion Environments

Abstract

To solve the problems of few available satellites, poor data quality, and insufficient monitoring accuracy in the complex occlusion environment during the deformation monitoring of the sluice, the present invention provides a Beidou / GNSS sluice deformation monitoring data processing method for the complex occlusion environment. According to the characteristics that the sluice deformation monitoring in the complex occlusion environment is seriously affected by the multipath effect and slow displacements are generated due to the influence of the upstream and downstream water levels and seasonal temperatures, a random model constrained by the carrier-to-noise ratio, robust estimation of the slow displacement constrained observations, and a half-sky model multipath correction method are used to obtain the solution results of the dynamic period coordinate sequence, and the refined processing is carried out on the dynamic time series results by extracting the period displacement based on the smooth prior analysis to obtain the period solution for real-time, continuous, and long-term monitoring of the deformation and displacement of the sluice. The sluice deformation monitoring method in the complex environment of the present invention can ultimately achieve high-precision sluice displacement monitoring with an accuracy of about 1 mm level.

Description

Technical Field

[0001] The present invention relates to the field of monitoring based on satellite positioning technology, and in particular to a Beidou / GNSS sluice deformation monitoring data processing method for complex shielding environments. Background Art

[0002] Sluice gates are a common hydraulic structure and an important facility for ensuring flood control and drainage safety. Their operation is affected by seasonal temperature changes and upstream and downstream water level changes, which will cause a certain amount of slow deformation. The use of Beidou / GNSS (Global Navigation Satellite System) technology can achieve real-time, time-based, and all-weather monitoring, and can provide timely warnings when sluice gates are deformed, providing a basis for decision-making.

[0003] The monitoring points of sluice gates are generally installed on load-bearing piers, with roads paved above or structures built for the management of sluice gates. This will lead to serious obstruction of the surrounding environment of the monitoring points, resulting in a reduction in the number of visible satellites, a deterioration in the satellite spatial distribution structure, a reduction in the quality of satellite data, and serious multipath errors. At the same time, the deformation process of the sluice gate is relatively slow, which places high requirements on the accuracy and precision of Beidou / GNSS positioning.

[0004] Therefore, a new sluice deformation data processing method is proposed to solve the above problems. Summary of the invention

[0005] The embodiment of the present application provides a Beidou / GNSS sluice deformation monitoring data processing method for complex obstruction environments, which is used to solve the problems of few available satellites, poor data quality, and insufficient monitoring accuracy in complex obstruction environments, and realize the monitoring function of the sluice with 1 mm time period solution.

[0006] The embodiment of the present application provides a Beidou / GNSS sluice deformation monitoring data processing method for complex occlusion environments, including:

[0007] S1, collect observation data of monitoring equipment for no less than 24 hours in an open environment, and analyze the functional relationship between the carrier-to-noise ratio and altitude angle of each satellite in the open environment;

[0008] S2, monitoring work is carried out in a sluice environment with complex actual obstruction, and the relationship between the download-to-noise ratio and the altitude angle in an open environment is used to determine in real time whether the satellite is available in the obstruction environment;

[0009] S3, screen healthy satellites of monitoring stations and reference stations, select common-view satellites, use the altitude sine function to express the random model, perform anti-error processing on the observation values ​​according to the slowly changing characteristics of the sluice gate, and form a double-difference observation equation;

[0010] S4, after performing Kalman filtering, LAMBDA search, and ambiguity fixing operations on the double difference observation equation, the multipath double difference residual is obtained;

[0011] S5, according to the single-double difference conversion formula, the multipath double difference residual is converted into a single difference residual, and the altitude angle, azimuth angle, and multipath single difference residual values ​​of different satellites at different times are output;

[0012] S6. Divide the visible range into 80×360 grid areas according to the satellite altitude angle and azimuth angle, and the size of each grid is 1°×1°; extract the trend of all data in each grid point, and average the trend items to obtain the multipath single difference correction value of the current grid point, which is used to correct the double difference equation of the next satellite cycle;

[0013] S7, receiving the observation data of the same position of the monitoring station in the next cycle in real time, judging whether the satellite is available and preprocessing the data according to the relationship between the download-to-noise ratio and the altitude angle in the open environment, and determining the healthy and available satellite;

[0014] S8, read the multipath correction value of the semi-spherical model, and search the model correction value according to the satellite altitude angle and azimuth angle;

[0015] S9, based on the basic solution model of Beidou / GNSS dynamic deformation monitoring, multipath model correction is performed to form a double difference equation, and Kalman filtering, LAMBDA search, and ambiguity fixing algorithm are performed to output the solution result of the current epoch;

[0016] S10, using the methods from S1 to S9 above, taking the actual observation data of the sluice as an example for analysis, three schemes are carried out to verify the effectiveness of the above-mentioned methods;

[0017] S11, according to the time series solution results, the smoothing prior analysis method is used to extract the trend;

[0018] S12, using the method in S11, analyzes the data of the sluice for one month, and conducts two schemes to verify the effectiveness of the above method.

[0019] Among them, in step S1, a polynomial fitting is used to express the relationship between the carrier-to-noise ratio and the altitude angle in an open environment. In the monitoring solution, if the difference between the carrier-to-noise ratio of a satellite with a fixed altitude angle and the carrier-to-noise ratio calculated by fitting in an open environment is greater than a threshold, the current satellite is eliminated. The functional relationship between the carrier-to-noise ratio and the altitude angle can be expressed as follows:

[0020] SNR(E)=α 0 +α 1 E+α 2 E 2

[0021] Where SNR(E) is the function relationship between carrier-to-noise ratio and altitude angle; α i is the fitting parameter; E is the satellite altitude angle.

[0022] Among them, in step S2, the observation data is received in the actual environment, the carrier-to-noise ratio SNR is analyzed and the altitude angle E is calculated, and the theoretical carrier-to-noise ratio value SNR corresponding to the current altitude angle is calculated according to the functional relationship between the carrier-to-noise ratio and the altitude angle. m and compare it with the actual carrier-to-noise ratio value SNR; if the actual carrier-to-noise ratio value SNR is lower than the minimum carrier-to-noise ratio threshold Th SNR If the actual carrier-to-noise ratio SNR is equal to the theoretical carrier-to-noise ratio SNR calculated by the model, it will be directly eliminated. Otherwise, it will proceed to the next step of judgment. m The difference is greater than the set threshold Th m If yes, remove it directly, otherwise go to the next step.

[0023] Among them, in step S3, Beidou / GNSS dynamic data processing adopts inter-station and inter-satellite double difference to solve and process, and the double difference observation equation is expressed as:

[0024]

[0025] In the formula, is the double difference factor; superscript i, j is the satellite number; subscript r, q is the station number, and p is the pseudorange observation value; is the carrier observation value; ρ is the geometric distance from the station to the satellite; λ is the wavelength; N is the integer ambiguity; T is the tropospheric delay error; I is the ionospheric delay error; ε p is the random error of pseudorange observation; is the random error of the carrier observation value; for the above pseudorange and carrier observation equations, as well as the satellite after the carrier-to-noise ratio constraint is eliminated, the altitude angle sine function is used to express its random model:

[0026] σ 2 =a 2 +b 2 / sin 2 (E)

[0027] Where a and b are empirical constants determined based on the receiver’s pseudorange and carrier observations; E is the satellite elevation angle; σ 2 is the variance; according to the slow transformation characteristics of the sluice itself, the previous epoch coordinates are used to verify the pseudorange and carrier error observation values ​​of the current epoch. The specific equation can be expressed as:

[0028] V k =B k X k-1 -L k

[0029] Where k is the epoch mark; V is the prior residual corresponding to the pseudorange and carrier; B is the design matrix of the current k epoch; X is the parameter estimation matrix of the previous epoch; L is the observation matrix of the current k epoch; after determining the observation model and the random model, the double difference observation values ​​of the healthy satellite are constructed, and finally the double difference observation equation is formed.

[0030] Among them, in step S4, after the double difference observation equation is formed, Kalman filtering is used for parameter estimation; the coordinate parameters are transitioned using process noise whose values ​​meet the preset conditions, and the ambiguity is kept as a time-invariant parameter; the LAMBDA algorithm is used for ambiguity search, and after the double difference satellite ambiguity is determined, the satellite double difference residual is calculated and can be expressed as:

[0031] V=BX-L

[0032] Where V is the double difference residual matrix; B is the design matrix; X is the parameter estimation matrix; and L is the observation matrix.

[0033] Among them, in step S5, according to the single-double difference conversion formula, the multipath double difference residual is converted into a single difference residual, and the altitude angle, azimuth angle, and multipath single difference residual values ​​of different satellites at different times are output; the multipath single difference residual sequence in a complex occlusion environment can be calculated according to the following method:

[0034] ZS=V

[0035] Where Z is the single-difference double-difference conversion matrix, and S is the single-difference residual matrix. For N satellites observed, only N-1 double-difference observations can be formed. Therefore, when converting double-difference observations to single-difference observations, it is impossible to perform an inversion operation through the single-double-difference conversion rectangle. Constraints need to be added to make the equation reversible. The specific matrix form is shown in the following formula:

[0036]

[0037] In the formula, p i is the weight factor corresponding to different satellites; is the single difference residual of different satellites; is the double difference residual with the reference satellite as serial number 1; after inverting the conversion matrix, the double difference residual V can be converted into the single difference residual S, and the azimuth, altitude angle and single difference residual of all satellites in all epochs are output and saved.

[0038] Among them, in step S6, the saved single difference residual sequence is read satellite by satellite, and all the data are distributed to different grids; the grid division basis of the semi-spherical model is based on the monitoring station as the center, the altitude angle of 90 degrees is directly above, the horizontal direction is 0 degrees, and it is divided into 90 parts according to the altitude angle; the azimuth angle is 0 degrees in the due north direction, and it increases in the clockwise direction until 360 degrees, and it is divided into 360 parts according to the azimuth angle to form a 1 degree × 1 degree grid point; the azimuth angle and altitude angle corresponding to the single difference residual value of the current satellite are rounded off and distributed to the corresponding grid; the trend of the sequence in each grid is extracted, and the trend items are averaged as the model correction value of the current grid point of the semi-spherical model,

[0039]

[0040] Where, E is the elevation angle of satellite i; A is the azimuth angle of satellite i; MHM i (E, A) is the multipath single-difference model correction value of satellite i at the grid point with altitude angle E and azimuth angle A; N is the number of epochs that satellite i stays on the (E, A) grid; V i (E,A) k is the kth single-difference residual of satellite i on the (E,A) grid.

[0041] Among them, in step S7, the principle of establishing the semi-celestial model for multipath error correction is to extract the multipath residual of one period for correcting the multipath error of subsequent periods based on the characteristics that the multipath is only related to the surrounding environment of the measuring station and the spatial position of the satellite and has temporal and spatial repeatability; obtain the correction value of the semi-celestial multipath model of each grid point established, perform multipath model correction and baseline solution on the basis of the Beidou / GNSS dynamic deformation monitoring basic solution model; obtain the base station data, monitoring station data, broadcast ephemeris and semi-celestial model multipath single difference correction value; perform data preprocessing and determine the available satellites.

[0042] Among them, in step S8, the altitude angle and azimuth angle of each satellite are calculated, and after rounding, they are searched and matched in the multipath correction value list to obtain the multipath single difference correction value of the current satellite at the current moment; the specific acquisition can be expressed by the following formula:

[0043] MHM=f(A,E)

[0044] In the formula, MHM i is the multipath single difference correction value of satellite i; fi is the search function of the current satellite.

[0045] Among them, in step S9, multipath model correction is performed on the basis of the basic solution model of Beidou / GNSS dynamic deformation monitoring to obtain the double difference equation, and then the double difference observation equation is solved by Kalman filtering, LAMBDA search, and ambiguity fixing algorithm, and then the current epoch solution is output, all epochs are traversed and solved, and finally all epoch results of the current period are output. The specific formula of the double difference observation equation can be expressed as:

[0046]

[0047] The Beidou / GNSS sluice deformation monitoring data processing method for complex shielding environments in the embodiment of the present application has the following beneficial effects:

[0048] This application solves the problems of few available satellites, poor data quality and insufficient monitoring accuracy in complex obstruction environments during sluice deformation monitoring, and realizes the monitoring function of sluices with 1 mm time period resolution. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 This is a flowchart of a Beidou / GNSS sluice deformation monitoring data processing method for a complex shielding environment according to an embodiment of the present application;

[0050] Figure 2 The carrier-to-noise ratio spatial map of the satellite visible to the sluice base station;

[0051] Figure 3 The spatial spectrum of the carrier-to-noise ratio of the satellite visible to the monitoring station upstream of the sluice gate;

[0052] Figure 4 The spatial spectrum of the carrier-to-noise ratio of the satellite visible to the monitoring station downstream of the sluice;

[0053] Figure 5 This is the spatial spectrum of single-difference residual multipath correction of the semi-celestial model at the monitoring point upstream of the sluice gate;

[0054] Figure 6 This is the spatial spectrum of single-difference residual multipath correction of the semi-celestial model at the monitoring point downstream of the sluice gate;

[0055] Figure 7 It is the error sequence diagram of three groups of experiments at the upstream monitoring points of the sluice gate;

[0056] Figure 8 It is the error sequence diagram of three groups of experiments at the monitoring points downstream of the sluice gate;

[0057] Fig. 9 Solve the error series for two groups of experimental periods at the monitoring point upstream of the sluice gate;

[0058] Fig.10 Solve the error sequence for two experimental periods at the monitoring point downstream of the sluice gate. DETAILED DESCRIPTION

[0059] The present application is further described below in conjunction with the accompanying drawings and embodiments.

[0060] In the following introduction, the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance. The following introduction provides multiple embodiments of the present invention, and different embodiments can be replaced or combined, so this application can also be considered to include all possible combinations of the same and / or different embodiments recorded. Therefore, if one embodiment includes features A, B, and C, and another embodiment includes features B and D, then this application should also be considered to include embodiments containing one or more of all other possible combinations of features A, B, C, and D, even though the embodiment may not be clearly recorded in the following text.

[0061] According to the characteristics that the deformation monitoring of sluice gates in complex shielding environments is seriously affected by the multipath effect, and the upstream and downstream water levels and seasonal temperatures produce slow displacements, the present invention adopts a random model with carrier-to-noise ratio constraints, anti-error of slow displacement constraint observation values, and a multipath correction method of a semi-celestial model to obtain the dynamic time period coordinate sequence solution results, and refines the time period displacement extraction based on smooth prior analysis of the dynamic time series results to obtain a time period solution for real-time, continuous and long-term monitoring of the deformation and displacement of the sluice gate. The experimental results using the measured data of the GNSS monitoring of the Wudingmen Regulating Gate as the solution case show that the GNSS deformation monitoring method of the sluice gate in a complex environment of the present application can ultimately achieve high-precision sluice gate displacement monitoring with an accuracy of about 1 mm.

[0062] like Figure 1 As shown, the Beidou / GNSS sluice deformation monitoring data processing method for complex occlusion environment of the present application includes the following steps:

[0063] S1. Collect observation data of monitoring equipment for no less than 24 hours in an open environment, and analyze the functional relationship between the carrier-to-noise ratio and altitude angle of each satellite in the open environment.

[0064] S2. Monitoring is carried out in a sluice environment with complex actual obstructions. The relationship between the download-to-noise ratio and the altitude angle in an open environment is used to determine in real time whether the satellite is available in the obstructed environment.

[0065] S3. Screen healthy satellites of monitoring stations and reference stations, select common-view satellites, use the sine function of the altitude angle to express the random model, perform anti-error processing on the observation values ​​according to the slowly changing characteristics of the sluice, and form a double-difference observation equation.

[0066] S4. After performing Kalman filtering, LAMBDA search, and ambiguity fixing operations on the double difference observation equation, the multipath double difference residual is obtained by solving.

[0067] S5. According to the single-double difference conversion formula, the multipath double difference residual is converted into a single difference residual, and the altitude angle, azimuth angle, and multipath single difference residual values ​​of different satellites at different times are output.

[0068] S6. Divide the visible range into 80×360 grid areas according to the satellite altitude angle (10-89°) and azimuth angle (0-359°), and the size of each grid is 1°×1°. Extract the trend of all data in each grid point, and average the trend items to obtain the multipath single difference correction value of the current grid point, which is used to correct the double difference equation of the next satellite cycle.

[0069] S7. Receive the observation data of the same position of the monitoring station in the next cycle in real time, judge whether the satellite is available according to the relationship between the download-to-noise ratio and the altitude angle in an open environment, perform data preprocessing such as gross error detection and elimination, cycle slip detection and elimination, and determine healthy and available satellites.

[0070] S8. Read the multipath correction value of the semi-celestial model and search for the model correction value according to the satellite altitude angle and azimuth angle.

[0071] S9. Based on the basic solution model of Beidou / GNSS dynamic deformation monitoring, multipath model correction is performed to form a double difference equation, and the solution result of the current epoch is output after performing Kalman filtering, LAMBDA search, ambiguity fixation and other algorithms.

[0072] S10. Using the methods from S1 to S9 above, the actual observation data of the sluice gate is analyzed as an example, and three schemes are carried out to verify the effectiveness of the above-mentioned methods.

[0073] S11. According to the time series solution results, the smoothing prior analysis (SPA) method is used to extract the trend.

[0074] S12. Using the method in S11, analyze the data of the sluice for one month, and conduct two schemes to verify the effectiveness of the above method.

[0075] A detailed introduction is given below.

[0076] S1. Collect observation data of monitoring equipment for no less than 24 hours in an open environment, and analyze the functional relationship between the carrier-to-noise ratio and the altitude angle of each satellite in the open environment. Use polynomial fitting to express the relationship between the carrier-to-noise ratio and the altitude angle in the open environment. If the difference between the carrier-to-noise ratio of a satellite with a fixed altitude angle and the carrier-to-noise ratio calculated by fitting in the open environment is greater than the threshold value in the monitoring solution, the current satellite will be eliminated. The functional relationship between the carrier-to-noise ratio and the altitude angle can be expressed as follows:

[0077] SNR(E)=α 0 +α 1 E+α 2 E2 (1)

[0078] Where SNR(E) is the function relationship between carrier-to-noise ratio and altitude angle; α i is the fitting parameter; E is the satellite altitude angle.

[0079] S2, Figure 2 is the reference station carrier-to-noise ratio spatial spectrum, Figure 3 It is the spatial spectrum of the carrier-to-noise ratio of the upstream monitoring point under the occlusion environment. Figure 4 is the spatial spectrum of the carrier-to-noise ratio of the downstream monitoring point under the occlusion environment, such as Figure 2-4 As shown in the figure, monitoring work is carried out in a sluice environment with complex obstructions, and the relationship between the carrier-to-noise ratio and the altitude angle in an open environment is used to determine in real time whether the satellite is available in the obstructed environment. Observation data is received in the actual environment, the carrier-to-noise ratio is analyzed and the altitude angle is calculated. According to the functional relationship between the carrier-to-noise ratio and the altitude angle, the theoretical carrier-to-noise ratio value corresponding to the current altitude angle is calculated and compared with the actual carrier-to-noise ratio value. If the actual carrier-to-noise ratio value is lower than the minimum carrier-to-noise ratio threshold, it is directly eliminated, otherwise it proceeds to the next step of judgment; if the difference between the actual carrier-to-noise ratio value and the theoretical carrier-to-noise ratio value calculated by the model is greater than the set threshold, it is directly eliminated, otherwise it proceeds to the next step.

[0080] S3, Beidou / GNSS dynamic data processing uses inter-station and inter-satellite double difference to solve and process. The double difference observation equation weakens the satellite orbit error, ionosphere and troposphere atmospheric delay error, eliminates the receiver and satellite clock error, reduces the parameters to be estimated, and improves the solution accuracy. Its double difference observation equation is expressed as:

[0081]

[0082] In the formula, is the double difference factor; superscript i, j is the satellite number; subscript r, q is the station number; p is the pseudorange observation value, in meters; is the carrier observation value, in cycles; ρ is the geometric distance from the station to the satellite, in meters; λ is the wavelength, in meters; N is the integer ambiguity, in cycles; T is the tropospheric delay error, in meters; I is the ionospheric delay error, in meters; ε p is the random error of pseudorange observation, in meters; is the random error of the carrier observation value, in meters. For the above pseudorange and carrier observation equations, as well as the satellite after the carrier-to-noise ratio constraint is eliminated, the altitude angle sine function is used to express its random model:

[0083] σ 2 =a 2 +b 2 / sin 2 (E) (3)

[0084] Where a and b are empirical constants determined based on the receiver’s pseudorange and carrier observations; E is the satellite elevation angle in radians; σ 2 According to the slow transformation characteristics of the sluice gate itself, the previous epoch coordinates are used to verify the pseudorange and carrier error observation values ​​of the current epoch. The specific equation can be expressed as:

[0085] V k =B k X k-1 -L k (4)

[0086] Where k is the epoch mark; V is the prior residual corresponding to the pseudorange and carrier; B is the design matrix of the current k epoch; X is the parameter estimation matrix of the previous epoch; and L is the observation matrix of the current k epoch. After determining the observation model and the random model, double-difference observations are constructed for the healthy satellite, and finally the double-difference observation equation is formed.

[0087] S4. Multipath error is caused by the superposition of the reflection and diffraction signals of the Beidou satellite signal caused by the surrounding environment of the measuring station, resulting in phase delay, which leads to deviation between the received observation value and the true value, resulting in multipath error. Generally, the higher the Beidou satellite altitude angle, the stronger the Beidou satellite signal, and the smaller the multipath effect. The multipath error is strongly related to the surrounding environment and the spatial position of the satellite, and cannot be weakened or eliminated by inter-satellite differential or inter-station differential. The impact of multipath error on the observation value can reach up to one-quarter of the wavelength. When performing RTK relative differential positioning, the maximum value of the multipath error can reach half the wavelength, which seriously affects the positioning accuracy. Therefore, for the high-precision positioning of the Wudingmen control gate in the complex obstruction environment, the satellite space domain sidereal day filter is used for correction.

[0088] After forming the double-difference observation equation, Kalman filtering is used for parameter estimation. Considering the small displacement characteristics between epochs, the coordinate parameters are transitioned using process noise with a small value, and the ambiguity is kept as a time-invariant parameter. The LAMBDA algorithm is used for ambiguity search. After determining the double-difference satellite ambiguity, the satellite double-difference residual can be expressed as:

[0089] V=BX-L (5)

[0090] Where V is the double difference residual matrix; B is the design matrix; X is the parameter estimation matrix; and L is the observation matrix.

[0091] S5. According to the single-double difference conversion formula, the multipath double difference residual is converted into a single difference residual, and the altitude angle, azimuth angle, and multipath single difference residual values ​​of different satellites at different times are output. The multipath single difference residual sequence in a complex occlusion environment can be calculated as follows:

[0092] ZS=V (6)

[0093] In the formula, Z is the single-difference double-difference conversion matrix, and S is the single-difference residual matrix. For N satellites observed, only N-1 double-difference observations can be formed. Therefore, when converting double-difference observations to single-difference observations, it is impossible to perform an inversion operation through the single-double-difference conversion rectangle, and constraints need to be added to make the equation reversible. The specific matrix form is shown in the following formula:

[0094]

[0095] In the formula, p i is the weight factor corresponding to different satellites; is the single difference residual of different satellites; is the double difference residual with the reference satellite as serial number 1. There are enough satellites observed in the current epoch, and the multipath utility error in the observation value can be regarded as random noise. The weighted average of the single difference residuals of the common-view satellites should be zero in theory, that is, Therefore, after inverting the conversion matrix, the double difference residual V can be converted into a single difference residual S, and the azimuth, altitude angle and single difference residual of all satellites in all epochs can be output and saved.

[0096] S6, Figure 5 It is the spatial spectrum of single-difference residual multipath correction of the semi-celestial model at the upstream monitoring point. Figure 6 It is the spatial spectrum of single-difference residual multipath correction of the semi-celestial model at the downstream monitoring point, such as Figure 5-6 As shown, read the saved single difference residual sequence satellite by satellite, and distribute all the data to different grids. The grid division basis of the semi-spherical model is based on the monitoring station as the center, with the altitude angle of 90 degrees directly above and 0 degrees horizontally, and divided into 90 parts according to the altitude angle; the azimuth angle is 0 degrees in the north direction, and increases in the clockwise direction until 360 degrees, and divided into 360 parts according to the azimuth angle, forming a 1 degree × 1 degree grid point. The azimuth and altitude angles corresponding to the single difference residual value of the current satellite are rounded off and distributed to the corresponding grid. The trend of the sequence in each grid is extracted, and the trend items are averaged as the model correction value of the current grid point of the semi-spherical model.

[0097]

[0098] Where, E is the elevation angle of satellite i; A is the azimuth angle of satellite i; MHM i (E, A) is the multipath single-difference model correction value of satellite i at the grid point with altitude angle E and azimuth angle A; N is the number of epochs that satellite i stays on the (E, A) grid; V i (E,A) k is the kth single-difference residual of satellite i on the (E,A) grid.

[0099] S7. The principle of establishing the semi-celestial model for multipath error correction is mainly based on the fact that multipath is only related to the surrounding environment of the measuring station and the spatial position of the satellite, and has the characteristics of temporal and spatial repeatability. The multipath residual of one cycle is extracted to correct the multipath error of the subsequent cycle and improve the dynamic positioning accuracy. The correction values ​​of the semi-celestial multipath model of each grid point are obtained, and the multipath model correction and baseline solution are performed based on the basic solution model of Beidou / GNSS dynamic deformation monitoring. The base station data, monitoring station data, broadcast ephemeris and semi-celestial model multipath single difference correction values ​​are obtained; data preprocessing is performed according to the steps of carrier-to-noise ratio verification, gross error detection and elimination, cycle slip detection and elimination, etc., to determine the available satellites.

[0100] S8. Calculate the altitude and azimuth of each satellite, round them off to the nearest integer, and search and match them in the multipath correction value list to obtain the multipath single difference correction value of the current satellite at the current moment. The specific acquisition can be expressed by the following formula:

[0101] MHM=f(A,E) (9)

[0102] In the formula, MHM i is the multipath single difference correction value of satellite i; f i Search function for the current satellite.

[0103] S9. Based on the basic solution model of Beidou / GNSS dynamic deformation monitoring, multipath model correction is performed to obtain the double difference equation. Then, the double difference observation equation is solved by Kalman filtering, LAMBDA search, ambiguity fixation and other algorithms, and then the current epoch solution is output. All epochs are traversed and solved, and finally all epoch results of the current period are output. The specific formula of the double difference observation equation can be expressed as:

[0104]

[0105] S10, Figure 7 is the error sequence diagram of three groups of experiments at the upstream monitoring points, Figure 8 is the error sequence diagram of three groups of experiments at downstream monitoring points, such as Figure 7-8 As shown in the figure, according to the established relationship between the carrier-to-noise ratio and the altitude angle in the open environment and the correction value of the semi-celestial multipath model in the spatial domain, the observation data of the upstream monitoring point and the downstream monitoring point of a sluice from June 1 to 2, 2023 are calculated, and the accuracy before and after the correction is compared and analyzed. Three schemes are used for experimental solution. Scheme 1 is not to perform pre-test judgment based on the slowly changing characteristics of the sluice, and no semi-celestial model correction is performed; Scheme 2 is to perform pre-test verification judgment based on the slowly changing characteristics of the sluice, and no semi-celestial model correction is performed; Scheme 3 is to perform carrier-to-noise ratio verification judgment based on the slowly changing characteristics of the sluice, and perform semi-celestial model correction. The specific statistical values ​​can be expressed as Table 1:

[0106] Table 1 Comparison of RTK solution results of three schemes

[0107]

[0108] It can be seen from the table that the statistical values ​​of ENU in the upstream and downstream directions of Scheme 1 are 21.5mm, 7.4mm, 33.5mm and 6.8mm, 4.4mm, 12.9mm respectively; the accuracy of the results of Scheme 2 at the upstream and downstream monitoring points is improved by more than 33% and 13% in each direction compared with Scheme 1; the accuracy of the results of Scheme 3 at the upstream and downstream monitoring points is improved by more than 30% and 40% in each direction compared with Scheme 2. Therefore, anti-error processing based on the slow deformation characteristics of the sluice and multi-path correction of the semi-celestial model are important parts of the dynamic solution of deformation monitoring under complex occlusion environments.

[0109] S11. Beidou / GNSS deformation monitoring time series can be regarded as a trend signal polluted by noise and gross errors. In order to accurately extract deformation signals and objectively reflect deformation trends, the smoothing prior analysis method is used to extract time period displacement of time series data, and judge the overall displacement trend of the time period, and finally obtain the deformation value of the current time period. In the Beidou / GNSS deformation monitoring time series, the noise and gross error parts show high-frequency characteristics, while the trend change part (including deformation trend and multipath effect) shows low-frequency long-period changes. Separating deformation trend and positioning error by differences in deformation monitoring data in frequency domain and other aspects is more resistant to the influence of noise and gross errors than traditional time domain analysis. At the same time, considering the multipath effect and the slow displacement characteristics of the sluice, the smoothing prior analysis method (SPA) is used to extract time period displacement. As an effective method for separating nonlinear trend items of signals, the smoothing prior analysis method can separate trend items of different smoothness degrees by adjusting only one parameter. The refined processing of the time period solution based on the dynamic coordinate sequence is based on days. Outlier detection and elimination and SPA smoothing prior analysis are implemented on the dynamic coordinate time series to obtain a daily smoothed sequence, and then the comprehensive solution for each day is analyzed and statistically analyzed. The daily solution results within one month are then calculated.

[0110] After removing the outliers, the time series is smoothed by SPA prior analysis and the trend item is extracted. The displacement change of each increasing or decreasing interval in the trend item is analyzed according to the extreme points of the trend item. If the change is too large, it will be treated as a gross error and removed. Otherwise, the average is taken as the change in the current interval, and the average of all interval changes is taken as the current time period solution.

[0111] S12, such as Figure 9-10As shown, the deformation monitoring system of the Wudingmen sluice gate has been running for a long time. The time period when the water level is relatively stable and no overall deformation occurs is selected for accuracy statistics. Finally, 30 groups of benchmark station observation data, upstream monitoring observation data, and downstream monitoring observation data are collected from June 1 to 30, 2023. The multi-path error correction method in a complex occlusion environment is used to obtain the time series results of the upstream monitoring points and the downstream monitoring points. Two schemes are used to experimentally analyze the accuracy of the results before and after the refinement of the time period solution based on the dynamic coordinate sequence. Scheme 1 is to directly average the time result sequence within one day as the comprehensive solution, and scheme 2 is to average the time period displacement extraction after SPA smoothing prior analysis as the comprehensive solution. The accuracy of the 30 groups of results from the two groups of experiments is analyzed. The specific data processing results are shown as follows:

[0112] Table 2 Comparison of monitoring point results of two schemes

[0113]

[0114] The effectiveness of the proposed method for fine-tuning the time period solution based on the dynamic coordinate sequence was verified through two sets of experiments. The accuracy of the upstream and downstream planes of Scheme 1 is within 1 mm, and the elevation accuracy is 1.1 mm and 1.0 mm respectively. The accuracy of the upstream and downstream planes of Scheme 2 is improved by 11% or more compared with Scheme 1, and the accuracy is within 0.8 mm, and the elevation accuracy is improved by 18% or more, and the accuracy is 0.9 mm and 0.7 mm respectively.

[0115] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A Beidou / GNSS sluice deformation monitoring data processing method for complex occlusion environment, characterized in that: include: S1, collect observation data of monitoring equipment for no less than 24 hours in an open environment, and analyze the functional relationship between the carrier-to-noise ratio and altitude angle of each satellite in the open environment; S2, monitoring work is carried out in a sluice environment with complex actual obstruction, and the relationship between the download-to-noise ratio and the altitude angle in an open environment is used to determine in real time whether the satellite is available in the obstruction environment; S3, screen healthy satellites of monitoring stations and reference stations, select common-view satellites, use the altitude sine function to express the random model, perform anti-error processing on the observation values ​​according to the slowly changing characteristics of the sluice gate, and form a double-difference observation equation; S4, after performing Kalman filtering, LAMBDA search, and ambiguity fixing operations on the double difference observation equation, the multipath double difference residual is obtained; S5, according to the single-double difference conversion formula, the multipath double difference residual is converted into a single difference residual, and the altitude angle, azimuth angle, and multipath single difference residual values ​​of different satellites at different times are output; S6. Divide the visible range into 80×360 grid areas according to the satellite altitude angle and azimuth angle, and the size of each grid is 1°×1°; extract the trend of all data in each grid point, and average the trend items to obtain the multipath single difference correction value of the current grid point, which is used to correct the double difference equation of the next satellite cycle; S7, receiving the observation data of the same position of the monitoring station in the next cycle in real time, judging whether the satellite is available and preprocessing the data according to the relationship between the download-to-noise ratio and the altitude angle in the open environment, and determining the healthy and available satellite; S8, read the multipath correction value of the semi-spherical model, and search the model correction value according to the satellite altitude angle and azimuth angle; S9, based on the basic solution model of Beidou / GNSS dynamic deformation monitoring, multipath model correction is performed to form a double difference equation, and Kalman filtering, LAMBDA search, and ambiguity fixing algorithm are performed to output the solution result of the current epoch; S10, using the methods of S1 to S9 above, taking the actual observation data of the sluice as an example to analyze and verify the effectiveness of the above methods; S11, according to the time series solution results, the smoothing prior analysis method is used to extract the trend; S12, using the method in S11, analyze the data of the sluice for one month and verify the effectiveness of the above method.

2. The Beidou / GNSS sluice deformation monitoring data processing method for complex shielding environments according to claim 1 is characterized in that: In step S1, a polynomial fitting is used to express the relationship between the carrier-to-noise ratio and the altitude angle in an open environment. In the monitoring solution, if the difference between the carrier-to-noise ratio of a satellite with a fixed altitude angle and the carrier-to-noise ratio calculated by fitting in an open environment is greater than a threshold, the current satellite is eliminated. The functional relationship between the carrier-to-noise ratio and the altitude angle can be expressed as follows: SNR(E)=α0+α1E+α2E 2 Where SNR(E) is the function relationship between carrier-to-noise ratio and altitude angle; α i is the fitting parameter; E is the satellite altitude angle.

3. The Beidou / GNSS sluice deformation monitoring data processing method for complex shielding environments according to claim 1 or 2 is characterized in that: In step S2, the observation data is received in the actual environment, the carrier-to-noise ratio SNR is analyzed and the altitude angle E is calculated. According to the functional relationship between the carrier-to-noise ratio and the altitude angle, the theoretical carrier-to-noise ratio value SNR corresponding to the current altitude angle is calculated. m and compare it with the actual carrier-to-noise ratio value SNR; if the actual carrier-to-noise ratio value SNR is lower than the minimum carrier-to-noise ratio threshold Th SNR If the actual carrier-to-noise ratio SNR is equal to the theoretical carrier-to-noise ratio SNR calculated by the model, it will be directly eliminated. Otherwise, it will proceed to the next step of judgment. m The difference is greater than the set threshold Th m If yes, remove it directly, otherwise go to the next step.

4. The Beidou / GNSS sluice deformation monitoring data processing method for complex shielding environments according to claim 1 or 2 is characterized in that: In step S3, Beidou / GNSS dynamic data processing uses inter-station and inter-satellite double differences for solution processing, and the double difference observation equation is expressed as: Where ▽Δ is the double difference factor; superscript i, j is the satellite number; subscript r, q is the station number, and p is the pseudorange observation value; is the carrier observation value; ρ is the geometric distance from the station to the satellite; λ is the wavelength; N is the integer ambiguity; T is the tropospheric delay error; I is the ionospheric delay error; ε p is the random error of pseudorange observation; is the random error of the carrier observation value; for the above pseudorange and carrier observation equations, as well as the satellite after the carrier-to-noise ratio constraint is eliminated, the altitude angle sine function is used to express its random model: s 2 =a 2 +b 2 / sin 2 (E) Where a and b are empirical constants determined based on the receiver’s pseudorange and carrier observations; E is the satellite elevation angle; σ 2 is the variance; according to the slow transformation characteristics of the sluice itself, the previous epoch coordinates are used to verify the pseudorange and carrier error observation values ​​of the current epoch. The specific equation can be expressed as: V k =B k X k-1 -L k Where k is the epoch mark; V is the prior residual corresponding to the pseudorange and carrier; B is the design matrix of the current k epoch; X is the parameter estimation matrix of the previous epoch; L is the observation matrix of the current k epoch; After determining the observation model and the random model, double-difference observations are constructed for the healthy satellite, and finally the double-difference observation equation is formed.

5. The Beidou / GNSS sluice deformation monitoring data processing method for complex shielding environments according to claim 1 or 2 is characterized in that: In step S4, after forming the double difference observation equation, Kalman filtering is used for parameter estimation; the coordinate parameters are transitioned using process noise whose values ​​meet the preset conditions, and the ambiguity is kept as a time-invariant parameter; the LAMBDA algorithm is used for ambiguity search, and after determining the double difference satellite ambiguity, the satellite double difference residual is calculated and can be expressed as: V=BX-L Where V is the double difference residual matrix; B is the design matrix; X is the parameter estimation matrix; and L is the observation matrix.

6. The Beidou / GNSS sluice deformation monitoring data processing method for complex shielding environments according to claim 1 or 2 is characterized in that: In step S5, the multipath double difference residual is converted into a single difference residual according to the single-double difference conversion formula, and the altitude angle, azimuth angle, and multipath single difference residual values ​​of different satellites at different times are output; the multipath single difference residual sequence in a complex occlusion environment can be calculated as follows: ZS=V Where Z is the single-difference double-difference conversion matrix, and S is the single-difference residual matrix. For N satellites observed, only N-1 double-difference observations can be formed. Therefore, when converting double-difference observations to single-difference observations, it is impossible to perform an inversion operation through the single-double-difference conversion rectangle. Constraints need to be added to make the equation reversible. The specific matrix form is shown in the following formula: In the formula, p i is the weight factor corresponding to different satellites; is the single difference residual of different satellites; is the double difference residual with the reference satellite as serial number 1; after inverting the conversion matrix, the double difference residual V can be converted into the single difference residual S, and the azimuth, altitude angle and single difference residual of all satellites in all epochs are output and saved.

7. The Beidou / GNSS sluice deformation monitoring data processing method for complex shielding environments according to claim 1 or 2 is characterized in that: In step S6, the stored single difference residual sequence is read satellite by satellite, and all data are allocated to different grids; the grid division basis of the semi-celestial model is based on the monitoring station as the center, the altitude angle of 90 degrees is directly above, the horizontal direction is 0 degrees, and the altitude angle is divided into 90 parts; the azimuth angle is 0 degrees in the due north direction, and increases in the clockwise direction until 360 degrees, and the azimuth angle is divided into 360 parts to form a 1 degree × 1 degree grid point; the azimuth angle and altitude angle corresponding to the single difference residual value of the current satellite are rounded to integers and allocated to the corresponding grid; The trend of the sequence in each grid is extracted, and the average of the trend items is taken as the model correction value of the current grid point of the semi-celestial model. Where, E is the elevation angle of satellite i; A is the azimuth angle of satellite i; MHM i (E,A) is the multipath single-difference model correction value of satellite i at the grid point with altitude angle E and azimuth angle A; N is the number of epochs that satellite i stays on the (E,A) grid; Vi(E,A) k is the kth single-difference residual of satellite i on the (E,A) grid.

8. The Beidou / GNSS sluice deformation monitoring data processing method for complex shielding environments according to claim 1 or 2 is characterized in that: In step S7, the principle of establishing a semi-celestial sphere model for multipath error correction is to extract the multipath residual of one cycle to correct the multipath error of subsequent cycles based on the fact that the multipath is only related to the surrounding environment of the measuring station and the spatial position of the satellite and has the characteristics of temporal and spatial repeatability; obtain the correction value of the semi-celestial sphere multipath model of each grid point established, perform multipath model correction and baseline solution on the basis of the Beidou / GNSS dynamic deformation monitoring basic solution model; Obtain reference station data, monitoring station data, broadcast ephemeris, and multipath single-difference correction values ​​of the semi-celestial model; Perform data preprocessing and determine available satellites.

9. The Beidou / GNSS sluice deformation monitoring data processing method for complex shielding environments according to claim 1 or 2, characterized in that: In step S8, the altitude angle and azimuth angle of each satellite are calculated, and after rounding, a search and match is performed in the multipath correction value list to obtain the multipath single difference correction value of the current satellite at the current moment; the specific acquisition can be expressed by the following formula: MHM=f(A,E) In the formula, MHM i is the multipath single difference correction value of satellite i; f i Search function for the current satellite.

10. The Beidou / GNSS sluice deformation monitoring data processing method for complex shielding environments according to claim 1 or 2, characterized in that: In step S9, Based on the basic solution model of Beidou / GNSS dynamic deformation monitoring, multipath model correction is performed to obtain the double difference equation. Then, the double difference observation equation is solved by Kalman filtering, LAMBDA search, and ambiguity fixing algorithm, and then the current epoch solution is output. All epochs are traversed and solved, and finally all epoch results of the current period are output. The specific formula of the double difference observation equation can be expressed as: In the formula, ▽Δ is the double difference factor; the superscript ij is the satellite number; the subscript rq is the station number; is the carrier observation value; ρ is the geometric distance from the station to the satellite; λ is the wavelength; N is the integer ambiguity; T is the tropospheric delay error; I is the ionospheric delay error; MHM is the multipath single difference correction value; is the random error of the carrier observation value.

Images