A water conservancy multi-parameter monitoring method based on GNSS

Abstract

The application provides a water conservancy multi-parameter monitoring method based on GNSS, comprising the following steps: 1) receiving GNSS satellite signals based on a GNSS receiving antenna and processing: extracting direct signals; extracting reflected signals; 2) monitoring water conservancy engineering displacement deformation based on the extracted direct signals; 3) monitoring water level based on the extracted reflected signals; and 4) inversing atmospheric precipitable water in real time based on GNSS troposphere delay information. The same receiving device is used for different signal processing, water conservancy engineering deformation displacement, water vapor, water level parameter information are quickly acquired and high-precision monitoring is realized, and the water conservancy monitoring sensing capability is comprehensively improved, so that scientific data support can be provided for flood control, drainage and other projects.

Description

TECHNICAL FIELD

[0001] The application relates to the field of GNSS water conservancy and hydrological monitoring technology, and particularly relates to a water conservancy multi-parameter monitoring method based on GNSS. BACKGROUND

[0002] In the field of water conservancy and hydrological monitoring, the traditional global satellite navigation system (GNSS) water conservancy monitoring technology can only obtain the displacement information of a point based on the pseudorange and carrier phase observation values of the GNSS direct signal. Therefore, the traditional water conservancy monitoring method usually relies on a single displacement parameter and is difficult to provide comprehensive monitoring data, resulting in delayed or inaccurate early warning.

[0003] Moreover, the complex environment of water conservancy projects often interferes with GNSS signals, resulting in problems such as multipath, non-line-of-sight signals, continuous gross errors and multiple gross errors, frequent cycle slips and the like in the observation values, which makes it difficult to correctly fix the ambiguity and directly reduces the precision and reliability of deformation monitoring.

[0004] Especially in the discovery of safety risk hidden dangers and disaster monitoring and early warning of water conservancy projects, the traditional GNSS water conservancy monitoring has the problems of low displacement monitoring precision and single monitoring information in complex environments such as mountain and tree shielding, and a multi-parameter GNSS monitoring technology needs to be proposed to realize rapid acquisition and high-precision monitoring of parameter information such as water conservancy project deformation displacement, water vapor and water level.

[0005] As a global satellite navigation system independently developed by China, the Beidou satellite system has the advantages of high precision, wide coverage and multi-source data integration, and it is of great significance to develop multi-parameter Beidou water conservancy measurement technology to realize rapid multi-parameter synchronous monitoring and improve the precision of measurement. SUMMARY

[0006] The application is aimed at the problems of low GNSS displacement monitoring precision and single Beidou monitoring information in complex environments such as mountain and tree shielding in the discovery of safety risk hidden dangers and disaster monitoring and early warning of water conservancy projects, and proposes a water conservancy multi-parameter monitoring method based on GNSS, which uses the same receiving device for different signal processing to realize rapid acquisition and high-precision monitoring of parameter information such as water conservancy project deformation displacement, atmospheric precipitation and water level, and comprehensively improves the water conservancy monitoring sensing capability. The method of the application is suitable for the Beidou satellite system.

[0007] The object of the application is realized by the following technical scheme:

[0008] A water conservancy multi-parameter monitoring method based on GNSS, comprising the following steps:

[0009] Step one, receiving GNSS satellite signals based on GNSS receiving antenna and processing: extracting direct signals; extracting signal-to-noise ratio data, removing direct signal part in signal-to-noise ratio data, and only keeping reflected signals;

[0010] Step two, water conservancy engineering displacement deformation monitoring based on direct signals extracted in step one: collecting GNSS observation data and screening key quality indicators; establishing template functions between different key quality indicators; judging the affected type of signals by the difference between the measured value and the template value of the quality indicators, and selecting the observation variance calculation model accordingly; bringing the calculated observation variance into the RTK (Real-Time Kinematic, i.e. real-time dynamic differential positioning) processing to obtain coordinates, and adjusting the weight of the quality indicators in the observation variance calculation model according to the post-adjustment variance, so as to build a comprehensive random model suitable for GNSS positioning in complex environments to calculate the coordinates of water conservancy projects and realize displacement deformation monitoring;

[0011] Step three, water level monitoring based on reflected signals extracted in step one: using inverse modeling method for water level inversion, and using B-spline curve to fit water level changes;

[0012] Step four, based on the direct signals extracted in step one, using precise point positioning technology to estimate GNSS tropospheric delay, and then inversely calculating the tropospheric delay information into atmospheric precipitable water content.

[0013] Further optimization, the GNSS receiving antenna in step one is a geodetic survey antenna.

[0014] Further, step two includes the following operations:

[0015] S21 collects GNSS raw observation data in different complex scenes, fits the functional relationship between observation variance and quality indicators, screens key quality indicators, and the quality indicators include elevation angle, carrier-to-noise ratio and PLD;

[0016] S22 establishes template functions between different key quality indicators based on different GNSS;

[0017] S23 judges the affected type of signals by the difference between the measured value and the template value of the quality indicators and determines the segmented weighting strategy to establish a segmented weighting model, substitutes the standardized difference data between the measured value and the template value of the quality indicators into the segmented weighting model, selects the observation variance calculation model based on the segment it is in, and calculates the variance;

[0018] S24 substitutes the variance calculated in S23 into the positioning equation to obtain coordinates to complete RTK positioning;

[0019] S25 adjusts the weight using the post-adjustment variance and updates the coordinates back to step S24 until the coordinates converge.

[0020] Further, in step S24, after obtaining the variance of the observation value, the variance of the non-difference (original) observation value also needs to be propagated to the double-difference observation value via the error propagation law, and combined with the function model to complete the RTK positioning.

[0021] Further, in step three, the water level monitoring includes post-inversion or real-time inversion.

[0022] The post-inversion: the LSP spectrum analysis method reads the SNR sequence of a single satellite, performs spectrum analysis to calculate the peak frequency, performs data quality control, and uses a B-spline curve to model the water level height changing with time as a smooth and continuous function.

[0023] The real-time inversion: first, accept real-time RTCM data stream from the receiver, perform data format conversion, SNR data extraction and detrending in software, calculate the reflected height data using satellite arcs, then according to the user set window size and node interval, use the inverse modeling method to fit the water level curve in each window and output the result, the last node water level in each window is the real-time inversion result; at the same time, record the nodes in the window, and finally construct a long water level time series through all the nodes.

[0024] Further, in step four, the following steps are included:

[0025] S41 According to the GPS / BDS precise point positioning, the coordinates, carrier phase ambiguity, troposphere zenith wet delay, and clock error are used as unknowns together with the Kalman filtering method to calculate and estimate the zenith troposphere delay ZTD; the zenith troposphere delay is decomposed into zenith hydrostatic delay ZHD and zenith wet delay ZWD two parts;

[0026] S42 Determine the height pressure and temperature data of four grid points near a GNSS station:

[0027]

[0028] In the formula: And Pi and Ti represent the adjusted pressure and temperature of the i-th grid point; And Pi-1 and Pi+1 represent the pressure values of the two closest horizontal points to the station; And Ti-1 and Ti+1 represent the two closest temperature values; And Pi-1 and Pi+1 represent the two closest horizontal potential heights; Hs is the elevation of the station.

[0029] After obtaining the pressure and temperature data of the four grid points near the GNSS station, the bilinear interpolation is performed to determine the ground pressure and weighted temperature in the station.

[0030] S43 Convert the ZTD data observed by the GNSS station into PWV using the ground pressure and weighted temperature obtained in the previous step.

[0031] Advantages of the present application:

[0032] The present application provides a GNSS-based water conservancy multi-parameter monitoring technology, which uses the same receiving device for different signal processing. Based on GNSS direct signal observation values, the displacement parameters of the point can be obtained by processing. In addition, the water level parameters of the nearby water surface can be obtained by GNSS reflection signal inversion. Furthermore, the atmospheric precipitable water content information can be extracted from the GNSS tropospheric delay information. At the same time, the Beidou deformation monitoring, Beidou water level inversion, and atmospheric water vapor inversion are introduced to monitor water conservancy together, which improves the accuracy and reliability of water conservancy GNSS (especially Beidou) monitoring in complex environmental conditions. It can provide scientific data support for flood control, drainage and other projects, and more efficient and low-cost solutions. Specifically, there are the following points:

[0033] (1) High precision and high reliability: The displacement monitoring method of the present application is aimed at the problems of multipath, non-line-of-sight signal, gross error, multiple gross error, and frequent cycle slip commonly found in complex engineering environments. Advanced signal processing and error correction methods are used to ensure that GNSS (especially Beidou system) can still provide high-precision and high-reliability monitoring results in complex environments. By optimizing the use of observation data and using uniformly spaced nodes, the lack of observation values or gross error is effectively avoided, ensuring the accuracy of the fitting results, so that the Beidou continuous static monitoring accuracy is stable at a high standard of 3 millimeters (RMS).

[0034] (2) High timeliness and real-time: Based on high precision, the method focuses on the timeliness of the monitoring results. When processing the near-real-time inversion of water level, the innovation of water level inversion technology enables the inversion results to be calculated and output in a very short time, greatly improving the inversion accuracy and timeliness. Compared with traditional monitoring methods, it can provide water level change data more quickly and accurately, meeting the needs of emergency disaster response.

[0035] (3) Precise weather data processing and water vapor inversion: To solve the problem of insufficient Beidou PWV (atmospheric precipitable water content) inversion accuracy, a zenith tropospheric delay (ZTD) calculation quality control method of precise point positioning is proposed, and a wet delay (ZWD) and PWV conversion coefficient model is established, making the water vapor inversion process more accurate. Solving the key problem of missing ground weather data, even in the case of incomplete data, Beidou system can still effectively monitor water vapor. BRIEF DESCRIPTION OF DRAWINGS

[0036] The accompanying drawings, which form a part of this application, are included to provide a further understanding of the application, illustrate the preferred embodiments of the application and assist in

[0037] Figure 1 The five-stage decision-making strategy in Example 1;

[0038] Figure 2 Comparison of RTK results based on different random models in the shade environment in Example 1;

[0039] Figure 3 Post-inversion water level results for the demonstration point in Example 1;

[0040] Figure 4 Near real-time inversion water level results for the demonstration point in Example 1;

[0041] Figure 5 PWV time series graph of four GNSS stations in the Hetao area of Bayannur, Inner Mongolia, in Example 1;

[0042] Figure 6 Comparison of ZTD estimated based on Beidou B2b and post-precise ephemeris (Bayannur Hetao JZ01) in Example 1;

[0043] Figure 7 Comparison of PWV estimated based on NCEP-GFS and empirical GPT model with sounding data (Bayannur JZ01 station) in Example 1;

[0044] Figure 8 Comparison of PWV estimated based on NCEP-GFS and empirical GPT model with sounding data (Bayannur JC01 station) in Example 1;

[0045] Figure 9 Comparison of PWV estimated based on NCEP-GFS and empirical GPT model with sounding data (Bayannur JC02 station) in Example 1;

[0046] Figure 10 Comparison of PWV estimated based on NCEP-GFS and empirical GPT model with sounding data (Bayannur JC04 station) in Example 1. DETAILED DESCRIPTION

[0047] Example 1

[0048] A GNSS-based water conservancy multi-parameter monitoring method, comprising the following steps:

[0049] Step one, receiving GNSS satellite signals based on GNSS receiving antenna and processing: extracting direct signal; extracting signal-to-noise ratio data, removing direct signal part in signal-to-noise ratio data, and only keeping reflected signal; extracting GNSS troposphere delay information; the GNSS receiving antenna is a geodetic survey antenna.

[0050] Step two, water conservancy engineering displacement deformation monitoring based on the direct signal extracted in step one: collecting GNSS observation data and screening key quality indicators; establishing template functions between different key quality indicators; using the difference between the measured value and the template value of the quality indicator to judge the affected type of the signal, and selecting the observation variance calculation model accordingly; bringing the calculated observation variance into the RTK processing to obtain the coordinates, and adjusting the weight of the quality indicator in the observation variance calculation model according to the post-check variance, so as to build a comprehensive random model suitable for GNSS positioning in complex environments to calculate the coordinates of water conservancy projects and realize displacement deformation monitoring; the specific operation of this embodiment is as follows:

[0051] S21 collects GNSS original observation data in different complex scenes, fits the functional relationship between observation variance and quality indicators, and screens key quality indicators; the key quality indicators screened in this embodiment include elevation angle, carrier-to-noise ratio and PLD.

[0052] First, collect original observation data in different complex scenes, and fit the functional relationship between observation variance and quality indicators.

[0053] ① The observation variance function based on the elevation angle is:

[0054]

[0055] Wherein, a and b are empirical coefficients, which are obtained by fitting method or configured as a=3mm, b=3mm, θ represents the satellite elevation angle, i represents the satellite identifier, and σ {·} is the standard deviation of the observation value.

[0056] ② The observation variance function based on the carrier-to-noise ratio is:

[0057]

[0058] Wherein, C / N0 is the carrier-to-noise ratio, the unit is dB-Hz, C is a constant, which is fitted or empirically set, for example, for GPS L1 observation, it can be set as 0.00224m 2 Hz.

[0059] ③ The observation variance function based on PLD is:

[0060]

[0061] Wherein For phase observation variance, subscript G and C represent GPS and BDS respectively. The corresponding pseudo-range variance can be obtained by using the ratio of phase variance:

[0062] ratio G = (-12.57 x PLD + 4236.48) / (PLD + 22.59 x PLD)

[0063]

[0064] S22 establishes a template function between different key quality indicators based on different GNSS;

[0065] In this embodiment, the GPS MEO carrier-to-noise ratio and its standard deviation template with elevation angle as the variable are:

[0066] C / N0 nom (θ) = 37.86 + 0.43 x θ - 5.02 x 10 -3 x θ 2

[0067] + 2.06 x 10 -5 x θ 3

[0068]

[0069] The BDS MEO satellite carrier-to-noise ratio and its standard deviation template are:

[0070] C / N0 nom (θ) = 39.45 + 0.44 x θ - 5.51 x 10 -3 x θ 2

[0071] + 2.46 x 10 -5 x θ 3

[0072]

[0073] The BDS IGSO satellite carrier-to-noise ratio and its standard deviation template are:

[0074] C / N0 nom (θ) = 39.14 + 0.17 x θ + 6.96 x 10 -4 x θ 2

[0075] - 1.37 x 10 -5 x θ 3

[0076]

[0077] In this embodiment, the GPS MEO PLD and its standard deviation template with carrier-to-noise ratio as a variable are as follows:

[0078] PLD nom (C / N0) = 4276.36 - 338.18 x C / N0 + 10.05 x C / N0 2

[0079] -0.13 x C / N0 3 +6.57 x 10 -4 x C / N0 4

[0080]

[0081] The BDS MEO satellite PLD and its standard deviation template are as follows:

[0082] PLD nom (C / N0) = 4160.74 - 330.82 x C / N0 + 9.89 x C / N0 2

[0083] -0.13 x C / N0 3 +6.57 x 10 -4 x C / N0 4

[0084]

[0085] The BDS IGSO satellite PLD and its standard deviation template are as follows:

[0086] PLD nom (C / N0) = 4630.39 - 371.82 x C / N0 + 15.06 x C / N0 2

[0087] -0.15 x C / N0 3 +7.58 x 10 -4 x C / N0 4

[0088]

[0089] S23 judges the affected type of signal and determines the segmented weighting strategy by using the difference between the measured value and the template value of the quality index, establishes a segmented weighting model, substitutes the standardized difference data between the measured value and the template value of the quality index into the segmented weighting model, selects an observation variance calculation model based on the segment where it is located, and calculates the variance.

[0090] Based on the correlation between the quality characteristics such as cycle slip ratio, data integrity rate, pseudorange multipath and the geometric quality index of elevation angle, and the baseband signal processing quality indexes such as carrier-to-noise ratio and PLD, it is concluded that the carrier-to-noise ratio can better reflect the phase tracking situation, the PLD can better reflect the cycle slip situation, and the elevation angle can reflect the trend that the observation error increases with the decrease of the elevation angle. Based on the above assumptions, a five-section modeling strategy is adopted: ① In open environment, the observation value is not disturbed by multipath, diffraction and other errors, and the function relationship between the observation value variance and the elevation angle is established; ② When the signal is disturbed by multipath and diffraction, the elevation angle model becomes inaccurate, and the function relationship between the observation value variance and the PLD and the carrier-to-noise ratio is established; ③ When the signal is mainly affected by cycle slip and unstable phase tracking, the function relationship between the observation value variance and the PLD, the elevation angle and the carrier-to-noise ratio is established; ④ When the signal is mainly affected by multipath, cycle slip and unstable phase tracking, the function relationship between the observation value variance and the PLD and the carrier-to-noise ratio is established; ⑤ When the signal is affected by strong multipath, cycle slip and unstable tracking loop, a larger observation value variance should be used to reduce the influence of errors on positioning.

[0091] The normalized carrier-to-noise ratio and PLD are calculated as follows:

[0092]

[0093] Wherein, the subscript mea represents the observation value, and the subscript nom represents the template value.

[0094] Substitute it into the following five-section weighting model, and the strategy diagram is shown in Figure 1 , wherein n0, n1, m0, m1 are empirical constants representing normal distribution quantile values, reflecting the abnormal degree of the index, n0 and m0 can be set to 1-2, n1 and m1 can be set to 2-3, and each section model is obtained by fitting the measured data:

[0095] 1) The first section of the weighting model: it is considered that the observation value is not disturbed by multipath, diffraction and other errors, and the elevation angle model is used:

[0096]

[0097] 2) The second section of the weighting model: when the normalized index falls into this interval, it is considered that the signal is disturbed by multipath and diffraction, and the elevation angle model is not accurate, and the following combined model is used:

[0098]

[0099] 3) The third section of the weighting model: when the normalized index falls into this interval, it is considered that the signal is mainly affected by cycle slip and unstable phase tracking, and the following comprehensive model is used:

[0100]

[0101] 4) Weighted model section 4: when the normalized index falls into this interval, it is considered that the signal is mainly affected by multipath, cycle slip and unstable phase tracking, and the following comprehensive model is used:

[0102]

[0103] 5) Weighted model section 5: at this time, it is considered that the signal is affected by strong multipath, cycle slip and unstable tracking loop, and the following model is used:

[0104]

[0105] S24 substitutes the variance calculated by S23 into the positioning equation to obtain coordinates to complete PKT positioning;

[0106] After obtaining the variance of the observation value, the variance of the non-difference (original) observation value needs to be propagated to the double-difference observation value through the error propagation law, and then combined with the function model to complete the RTK positioning. The specific process is as follows:

[0107] In RTK data processing, double-difference data is often used for data solving. Therefore, after obtaining the weight of a single observation data, the error propagation law is used to construct the weight matrix of the double-difference observation. The variance and covariance matrix of the non-difference observation value can be represented as:

[0108]

[0109] In the formula, n is the number of observed satellites. At this time, the single-difference observation vector SD between stations and its variance and covariance matrix D SD can be represented as:

[0110] SD=C·O, where C=(-E E)

[0111]

[0112] In the formula, O is the non-difference observation vector, and E is the unit matrix. represents the variance of the corresponding satellite of the monitoring station, represents the variance of the corresponding satellite of the reference station. Based on the single-difference observation equation, a reference satellite is selected for inter-satellite difference. Taking the first satellite as the reference satellite as an example, the double-difference observation value vector and its variance and covariance matrix can be represented as:

[0113] DD=C d ·SD,C d =(-I E)

[0114] D DD =C d ·D SD ·C d T

[0115]

[0116] where I is a vector with elements of 1.

[0117] Neglecting multipath and residual atmospheric errors, the double-difference pseudorange and carrier phase observation equations between stations a, b and satellites k and l can be expressed as:

[0118]

[0119] where, is the double-difference operator, ρ is the geodetic distance, ε is the observation noise, and N is the integer ambiguity. The above are the double-difference function model and double-difference stochastic model required for RTK positioning. With these, subsequent steps such as Kalman filtering, ambiguity fixing, and checking can be followed to obtain the coordinates of the rover station.

[0120] S25 uses the posterior variance to adjust the weight and back-substitutes step S24 to update the coordinates until the coordinates converge.

[0121] In order to achieve better robustness, Huber variance inflation function is used, and on the basis of the five-section model, the observation weight is adjusted according to the posterior residual:

[0122]

[0123] where, γ ii is the variance inflation factor, c0 is a predetermined empirical constant, for example, 2.5 to 3.0, is the standardized residual.

[0124] Experimental verification:

[0125] In order to verify the performance of the proposed stochastic model for RTK positioning, one hour of GPS+Beidou L1 / B1 data was collected in a tree-shaded environment. A special multi-frequency multi-GNSS receiver that can output PLD measurements was used for static data collection. The three panels from left to right in the figure represent the errors in the east, north, and up directions, respectively. It can be observed that the robust comprehensive model (R-compre, i.e., five-section weighting model + Huber) and the comprehensive stochastic model (Compre, i.e., five-section weighting model) have higher ambiguity success rate and positioning accuracy than other stochastic models.

[0126] Compared with the elevation angle, carrier-to-noise ratio and PLD random models, the fuzzy success rate of the comprehensive random model is increased by 50.04%, 17.62% and 8.46%, respectively. Compared with the comprehensive random model, the fuzzy success rate of the robust comprehensive random model is increased by 4.65%. In particular, the fuzzy success rate of the elevation angle random model is 0%, which may be due to the existence of a large number of occlusions and cycle slips in the observation data of the tree shade environment. In this case, it is not appropriate to use only the elevation angle to determine the random model. In terms of positioning accuracy, when the ambiguity is successfully fixed, all random models (except the elevation angle random model) achieve centimeter-level positioning accuracy in three directions, such as Figure 2 and Table 1. The model integrates the three data quality indicators of elevation angle, carrier-to-noise ratio and PLD, adjusts the weight determination method according to the difference between the measured value and the template value of each indicator, and adjusts the priori weight using residual information, which can effectively solve the problem of invalid weight determination of a certain indicator in a complex environment, ensure that the observation values are given reasonable weights during the entire calculation period, and through comparative analysis of the RTK positioning results of different random models, the model can effectively improve the ambiguity fixing success rate while ensuring the positioning accuracy.

[0127] Table 1 RTK positioning effect statistics comparison of different random models in tree shade environment

[0128]

[0129]

[0130] Step three, water level monitoring based on the reflected signal extracted in step one: using the inverse modeling method to perform water level inversion, and using B-spline curve to fit the water level change.

[0131] In this embodiment, the preprocessed GNSS signal-to-noise ratio (SNR) data is read, and the post-processing or near real-time inversion method is selected according to user requirements. The common method for Beidou-R water level inversion is to analyze the spectral characteristics of SNR, and the results are greatly affected by noise, and the time distribution of the output results is uneven. The present application uses the inverse modeling method to perform water level inversion, which uses B-spline curve to fit the water level change, which can significantly improve the accuracy and stability of the inversion results.

[0132] ① Post-processing: the LSP spectral analysis method reads the SNR sequence of a single satellite, calculates the peak frequency by spectral analysis, and then performs data quality control. The inverse modeling method does not need to analyze each satellite, but uses B-spline curve to model the water level height changing with time as a smooth and continuous function.

[0133] In this embodiment, the specific operation is as follows:

[0134] After receiving the GNSS satellite signal data, the signal-to-noise ratio data contained therein is extracted. A low-order polynomial is used to remove the direct signal portion of the signal-to-noise ratio data, retaining only the reflected signal portion. A mathematical model is then established between the reflected signal portion and the reflection height. First, the signal-to-noise ratio data after detrending is modeled as follows:

[0135]

[0136] Where δSNR is the signal-to-noise ratio after removing the direct signal; λ is the satellite signal wavelength; h is the reflection height; θ is the satellite elevation angle; k is the wave number; s is the roughness parameter of the reflecting surface; C1 and C2 are the in-phase and out-of-phase components, which are used to replace the amplitude A and phase The conversion relationship is:

[0137]

[0138] Before applying the inverse modeling method, it is necessary to use the LSP spectrum analysis method to perform a preliminary analysis of the δSNR data to determine the C1, C2, h, s in the mathematical model. 2 Initial values ​​of each parameter.

[0139] Before using B-spline curve to fit the water level change curve, it is necessary to select n curve nodes at equal time intervals. The nodes are represented by P i Using the observed data and the previously established δSNR mathematical model, a function y is fitted. i (x) Estimate the unknown parameters, where x is the unknown parameter, namely C1, C2, P i , roughness parameter s 2 , node P i , the function can be expressed as:

[0140]

[0141] Where i represents any integer between 1 and n, indicating the ordinal number of the node; δSNR i is the observed data, i.e. the signal-to-noise ratio data after removing the direct signal. The initial reflection height is calculated using step LSP spectrum analysis, and then the node spacing of the B-spline curve is determined. The parameters at the nodes are initialized according to the LSP results, and the undetermined parameters are estimated using the nonlinear least squares method. Through repeated iterations, y i The optimal parameter solution is obtained after the residual sum of squares of (x) is minimized, that is:

[0142]

[0143] In the formula, n is the number of SNR data, and after iteration is completed, a uniform water level change sequence can be obtained by B-spline interpolation using the estimated value h at the node. In the application, the function of the B-spline curve of the water (tide) level change can be expressed as:

[0144]

[0145] In the formula, h(t) is the B-spline curve of the water level change with time t; is the B-spline base function of p order, which can be expressed as:

[0146]

[0147] In the formula, u is a parameter in the node space, that is, u [t0, t n ] ; t i represents the time in the time window [t0, t n ]; p represents the order of the spline curve, and since the water level is continuously changed with time, the order of the B-spline base function is selected as two in the application.

[0148] The application intends to improve the inversion accuracy of the Beidou-R post-processing inversion by using the above method, and break through the high-precision Beidou-R water level post-processing inversion technology.

[0149] ②Real-time inversion: first, accept real-time RTCM data stream from the receiver, perform data format conversion, SNR data extraction and detrending in the software, and perform data preprocessing operations such as satellite arc segment calculation and reflection height calculation. Then, according to the window size and node interval set by the user, the water level curve in each window is fitted using the inverse modeling method and the result is output, and the last node water level of each window is the real-time inversion result. At the same time, the nodes in the window are recorded and saved, and finally a long water level time sequence is constructed by all the nodes.

[0150] The data of the demonstration point JC01 constructed in the second canal of the Hetao Irrigation District main canal are used to verify the post-processing algorithm and near real-time algorithm. The data of four frequency bands B1, B2a, B2b and B3 of the Beidou system are used for water level inversion examples.

[0151] The JC01 station is located at north latitude 40.66° and east longitude 107.26°, and the observation environment of JC01 is good, has an open monitoring range, and there is no high building and tree shelter around, which can capture rich water surface reflection signals and ensure the sufficiency of data sources. Considering the actual situation, it is necessary to artificially limit the range of satellite elevation angle and satellite azimuth angle. According to the distribution of color bands in the Fresnel reflection zone of JC01 station under the map view, it can be determined that the reflection signals from the water area are received in the azimuth angle range of 135°-220° and the elevation angle range of 7°-25°.

[0152] (One)Demonstration effect of post inversion

[0153] The post inversion results are shown in Figure 3 , the horizontal coordinate represents time, specifically in the form of accumulated days, and the vertical coordinate represents the change of water level, the blue curve is the measured water level value with a sampling rate of 5 minutes at the tide station, and the red curve is the water level value inverted by the post inversion algorithm. It can be found that the post inversion experimental results and the measured data of the tide station have good consistency. The water level inversion error of this site is 2.3 cm, and the correlation coefficient reaches 0.99, which meets the project requirement standard and can effectively reflect the water level change trend of the water area.

[0154] (Two)Demonstration effect of near real-time inversion

[0155] Figure 4 The near real-time algorithm inversion result is shown in the figure. The inversion result accuracy of the near real-time algorithm is 4.3 cm, and the correlation coefficient is 0.97. Due to the shorter data length used by the near real-time algorithm and the smaller area of the water area in this region, the inversion result has a certain deviation, but it can still achieve accurate matching with the measured data of the tide station.

[0156] Step four, based on the direct signal extracted in step one, the GNSS tropospheric delay is estimated by using precise point positioning technology, and then the tropospheric delay information is inverted into atmospheric precipitable water content in real time.

[0157] First, for GNSS / Beidou multi-frequency carrier phase observation data, based on precise point positioning technology, the coordinates, carrier phase ambiguity, tropospheric zenith wet delay, clock error, etc. are calculated together using Kalman filtering method to estimate the tropospheric delay ZTD. The GNSS tropospheric delay correction can be expressed as:

[0158] ZTD = d dry ·m dry +d wet ·m wet +d gN ·m gN +d gE ·m gE

[0159] In the formula, d dry / wet is the zenith tropospheric dry / wet component delay, m dry / wet is the dry / wet component projection function, d gN / gE is the horizontal gradient delay, and m gN / gE is the dry / wet component projection function.for the corresponding projection function. In precise point positioning, the tropospheric dry delay is corrected by the tropospheric delay correction model, and the zenith tropospheric wet delay is estimated by the random walk process. The zenith tropospheric delay can be decomposed into two parts, the zenith hydrostatic delay ZHD and the zenith wet delay ZWD. Among them, ZHD can be accurately obtained by using the empirical model supplemented by ground pressure data. The specific formula is:

[0160]

[0161] where P s : ground pressure, latitude of the station, Hs: elevation of the station.

[0162] ZTD-ZHD=ZWD

[0163] Further, the conversion coefficient can be used to directly convert ZWD into PWV:

[0164] PWV=Π×ZWD

[0165] where Π is a dimensionless conversion coefficient, which is related to the weighted mean temperature T m of the atmosphere.

[0166]

[0167] where the physical constants k3 and k'2 are 3.7546×10 5 K 2 / hPa and 22.9721K / hPa, respectively. T m is defined as the integral of the humidity and the corresponding temperature at different heights of the vertical atmospheric column above the station, and the calculation formula is:

[0168]

[0169] where P v is the water vapor pressure, T is the absolute temperature, and z is the vertical height. The vertical humidity and temperature profile can be extracted from sounding and reanalysis data. However, in most cases, it is difficult to obtain real-time atmospheric vertical profile data, and the T m estimation model based on ground meteorological parameters is widely used in the meteorological field.

[0170] As described above, in the process of converting the tropospheric delay ZTD obtained by GNSS inversion into atmospheric precipitable water vapor PWV, two meteorological data, ground pressure P s and weighted temperature T m , are required. When the GNSS station lacks meteorological sensors, the ZTD-PWV inversion can usually be performed with the aid of empirical meteorological parameter models, but the accuracy is relatively low. The present application proposes a real-time PWV inversion technology based on numerical prediction model products.

[0171] The NCEP-GFS forecast product data is distributed in the form of a grid on the global latitude and longitude grid, and the PWV inverted from GNSS data reflects the atmospheric precipitable water above the site. Therefore, the NCEP GFS forecast PWV value for a GNSS station is obtained by bilinear interpolation of the GFS data of the four adjacent grid points, and the surface pressure and temperature are determined by vertical adjustment and horizontal interpolation. First, the height pressure and temperature data of the four grid points near the station are determined as follows:

[0172]

[0173] In the formula: and represent the adjusted pressure and temperature of the i-th grid point; and represent the two closest horizontal pressure values from the station; and represent the two closest temperature values; and represent the two closest horizontal potential heights. Once the pressure and temperature data of the four adjacent grid points are obtained, bilinear interpolation is performed to determine the surface pressure and weighted temperature at the station, and then the ZTD data observed by the GNSS station can be converted to PWV.

[0174] Four demonstration points (JZ01, JC01, JCC02, and JC04) have been deployed in the Hetao area of Bayannur, Inner Mongolia, and four Beidou / GNSS receivers have been installed. Figure 5 The PWV time series inverted from the data of four GNSS stations is shown. The region is dry and rainy, and the altitude is high, so the atmospheric water vapor content is small.

[0175] Based on the 10-day data (from October 29 to November 7, 2024) collected by the GNSS base station in the Hetao area of Bayannur, Inner Mongolia, the ZTD estimated by PPP using precise ephemeris is compared with the ZTD product estimated by B2b receiver output. The results (see below Figure 6 ) show that the RMS value of the ZTD product estimated based on B2b is 8.1 mm, which confirms that the accuracy of the real-time PPP estimated ZTD data based on Beidou B2b product can meet the needs of this project.

[0176] In addition, the atmospheric precipitable water vapor (PWV) data retrieved from the Bayannur sounding station is used to evaluate the accuracy of the PWV retrieved from the four GNSS demonstration stations (JZ01, JC01, JC02, and JC04). The sounding station is 12 kilometers away from the GNSS stations, which is very suitable for evaluating the accuracy of GNSS retrieval of water vapor. Since the four GNSS demonstration stations do not have co-located meteorological observation data, the project uses the numerical weather prediction model NCEP-GFS meteorological parameter prediction product to realize the conversion of ZTD-PWV. As shown below, at the JZ01 station, the RMS value of the PWV product converted based on the numerical prediction model NCEP-GFS is 0.88 mm, while the RMS value of the PWV product converted based on the empirical meteorological model GPT2w is 2.13 mm, and the deviation at some time exceeds 3 mm. Figure 7 As shown below, at the JZ01 station, the RMS value of the PWV product converted based on the numerical prediction model NCEP-GFS is 0.88 mm, while the RMS value of the PWV product converted based on the empirical meteorological model GPT2w is 2.13 mm, and the deviation at some time exceeds 3 mm. Figures 8-10 As shown below, at the JZ01 station, the RMS value of the PWV product converted based on the numerical prediction model NCEP-GFS is 0.88 mm, while the RMS value of the PWV product converted based on the empirical meteorological model GPT2w is 2.13 mm, and the deviation at some time exceeds 3 mm.

[0177] Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present application shall be included in the protection scope of the present application.

Claims

1. A GNSS-based water conservancy multi-parameter monitoring method, characterized in that: The following steps are included Step 1: Receive GNSS satellite signals using a GNSS receiving antenna and process them: extract direct signals; extract signal-to-noise ratio data, remove the direct signal portion from the signal-to-noise ratio data, and retain only the reflected signals; Step 2: Monitor the displacement and deformation of the water conservancy project based on the direct signal extracted in step 1: Collect GNSS observation data and screen key quality indicators; establish template functions between different key quality indicators; use the difference between the measured value of the quality indicator and the template value to determine the type of signal affected, and select the observation variance calculation model accordingly; bring the calculated observation variance into RTK processing to obtain coordinates, and adjust the weight of the quality indicator in the observation variance calculation model according to the posterior variance. In this way, a comprehensive random model suitable for GNSS positioning in complex environments is constructed to calculate the coordinates of the water conservancy project and realize displacement and deformation monitoring; The following operations are included: S21 collects GNSS raw observation data from different complex scenarios, fits the functional relationship between observation variance and quality indicators, and screens key quality indicators, including altitude angle, carrier-to-noise ratio, and PLD; S22 establishes template functions between different key quality indicators based on different GNSS; S23 uses the difference between the measured value of the quality indicator and the template value to determine the affected type of the signal and determine the segment weighting strategy to establish a segment weighting model. Substitute the standardized difference data between the measured value of the quality indicator and the template value into the segment weighting model, select an observation variance calculation model based on the segment, and calculate the variance; S24 substitutes the variance calculated in S23 into the positioning equation to obtain the coordinates and complete the RTK positioning; S25 uses the posterior variance to adjust the weights and back-substitute step S24 to update the coordinates until the coordinates converge; Step 3: Monitor the water level based on the reflection signal extracted in step 1: use the inverse modeling method to perform water level inversion, and use the B-spline curve to fit the water level change; Step 4: Based on the direct signal extracted in step 1, estimate the GNSS tropospheric delay, and then invert the tropospheric delay information into atmospheric precipitable water in real time.

2. The GNSS-based water conservancy multi-parameter monitoring method according to claim 1, characterized in that: The GNSS receiving antenna described in step 1 is a geodetic measurement antenna.

3. The GNSS-based water conservancy multi-parameter monitoring method according to claim 1, characterized in that: In step S24, after obtaining the variance of the observation value, it is also necessary to propagate the variance of the non-difference observation value to the double-difference observation value through the error propagation law, and then combine it with the function model to complete the RTK positioning.

4. The GNSS-based water conservancy multi-parameter monitoring method according to claim 1, characterized in that: In step 3: water level monitoring including post-inversion or real-time inversion; The post-inversion method is as follows: the LSP spectrum analysis method reads the SNR sequence of a single satellite, performs spectrum analysis to calculate the peak frequency, performs data quality control, and uses a B-spline curve to model the water level height that varies with time as a smooth and continuous function; The real-time inversion method first receives the real-time RTCM data stream from the receiver, performs data format conversion, SNR data extraction and detrending in the software, calculates the reflection height data using the satellite arc segment, and then fits the water level curve in each window using the inverse modeling method according to the window size and node interval set by the user and outputs the result. At the same time, the node records in the window are saved, and finally a long water level time series is constructed through all the nodes.

5. The GNSS-based water conservancy multi-parameter monitoring method according to claim 1, characterized in that: Step 4 includes the following steps: Based on GPS / BDS precise point positioning, the S41 uses the coordinates, carrier phase ambiguity, zenith wet tropospheric delay, and clock error as unknowns and uses the Kalman filter method to calculate the zenith tropospheric delay (ZTD). The zenith tropospheric delay is decomposed into two parts: the zenith static delay (ZHD) and the zenith wet delay (ZWD). S42 determines the altitude, pressure and temperature data of four grid points near a GNSS station: Where: and represents the adjusted pressure and temperature of the i-th grid point; and Indicates the pressure values ​​of the two levels closest to the station; and Indicates the two closest temperature values; and Indicates the potential height of the two closest levels; Hs is the elevation of the measuring station; Once the pressure and temperature data of four grid points near the GNSS station are obtained, bilinear interpolation is performed to determine the ground pressure and weighted temperature within the station; S43 uses the surface pressure and weighted temperature calculated in the previous step to convert the ZTD data observed by the GNSS station into atmospheric precipitable water PWV.

Images