Water conservancy multi-parameter monitoring method based on GNSS

Abstract

The invention provides a water conservancy multi-parameter monitoring method based on a GNSS. The water conservancy multi-parameter monitoring method based on the GNSS comprises the following steps: 1) receiving a GNSS satellite signal based on a GNSS receiving antenna and processing the GNSS satellite signal: extracting a direct signal; extracting a reflected signal; 2) monitoring displacement deformation of the water conservancy project based on the extracted direct signal; 3) monitoring the water level based on the extracted reflected signal; and 4) performing real-time inversion on the atmospheric precipitable water amount based on the GNSS troposphere delay information. Different signal processing is carried out by adopting the same receiving device, rapid acquisition and high-precision monitoring of deformation displacement, water vapor and water level parameter information of the water conservancy project are realized, the water conservancy monitoring perception capability is comprehensively improved, and scientific data support can be provided for projects such as flood control and flood drainage.

Description

Technical Field

[0001] The present invention relates to the technical field of GNSS water conservancy engineering and hydrological monitoring, and in particular to a GNSS-based water conservancy multi-parameter monitoring method. Background Art

[0002] In the field of water conservancy projects and hydrological monitoring, traditional global satellite navigation system (GNSS) water conservancy monitoring technology can only process the displacement information of the point based on the pseudo-range and carrier phase observations of the GNSS direct signal. Therefore, traditional water conservancy monitoring methods usually rely on a single displacement parameter and are difficult to provide comprehensive monitoring data, resulting in delayed or inaccurate early warnings.

[0003] Moreover, the complex environment of water conservancy projects often interferes with GNSS signals, leading to problems such as multipath, non-line-of-sight signals, continuous gross errors and multiple gross errors, and frequent cycle slips in the observations. In severe cases, the ambiguity is difficult to fix correctly, which directly reduces the accuracy and reliability of deformation monitoring.

[0004] Especially in the discovery of safety risks and hidden dangers in water conservancy projects and disaster monitoring and early warning, traditional GNSS water conservancy monitoring has obvious problems of low displacement monitoring accuracy and single monitoring information in complex environments such as mountains and trees. It is necessary to propose a GNSS-based multi-parameter monitoring technology to achieve rapid acquisition and high-precision monitoring of parameter information such as deformation displacement, water vapor, and water level of water conservancy projects.

[0005] As a global satellite navigation system independently developed by my country, the Beidou satellite system has the advantages of high precision, wide coverage, and multi-source data integration. It is of great significance to develop multi-parameter Beidou water conservancy measurement technology to quickly and synchronously monitor multiple parameters and improve the measurement accuracy. Summary of the invention

[0006] The present invention aims at the problems of low GNSS displacement monitoring accuracy and single Beidou monitoring information in the discovery of safety risk hazards and disaster monitoring and early warning of water conservancy projects under complex environments such as mountains and trees, and proposes a GNSS-based water conservancy multi-parameter monitoring method, which uses the same receiving device to perform different signal processing to achieve rapid acquisition and high-precision monitoring of parameter information such as deformation displacement of water conservancy projects, atmospheric precipitation, and water level, and comprehensively improves the water conservancy monitoring perception ability. The method of the present invention is applicable to the Beidou satellite system.

[0007] The objective of the present invention is achieved through the following technical solutions:

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

[0009] Step 1: Receive GNSS satellite signals based on a GNSS receiving antenna and process them: extract direct signals; extract signal-to-noise ratio data, remove the direct signal part in the signal-to-noise ratio data, and only retain the reflected signal;

[0010] Step 2: Monitor the displacement and deformation of water conservancy projects 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 judge the affected type of the signal, and select the observation variance calculation model accordingly; bring the calculated observation variance into RTK (Real-Time Kinematic, i.e. real-time dynamic differential positioning) processing to obtain coordinates, and adjust the weight of the quality indicator in the observation variance calculation model according to the posterior variance, so as to construct a comprehensive random model suitable for GNSS positioning in complex environments to calculate the coordinates of water conservancy projects and realize displacement and deformation monitoring;

[0011] Step 3: Monitor the water level based on the reflection signal extracted in step 1: use the reverse modeling method to invert the water level, and use the B-spline curve to fit the water level change;

[0012] Step 4: Based on the direct signal extracted in step 1, the GNSS tropospheric delay is estimated using precise point positioning technology, and then the tropospheric delay information is inverted into atmospheric precipitable water in real time.

[0013] Further optimization, the GNSS receiving antenna described in step 1 is a geodetic measurement antenna.

[0014] Furthermore, step 2 includes the following operations:

[0015] S21 collects GNSS raw observation data of 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;

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

[0017] S23 uses the difference between the measured value of the quality indicator and the template value to judge the affected type of the signal and determine the segment weighting strategy to establish a segment weighting model, substitutes the standardized difference data between the measured value of the quality indicator and the template value into the segment weighting model, selects the observation variance calculation model based on the segment, and calculates the variance;

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

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

[0020] In a further step S24, after obtaining the variance of the observation value, it is also necessary to propagate the variance of the non-difference (original) observation value to the double-difference observation value via the error propagation law, and then combine the function model to complete the RTK positioning.

[0021] Further, in step three: performing water level monitoring including post-inversion or real-time inversion;

[0022] The post-inversion: 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 B-spline curves to model the water level height that changes with time as a smooth and continuous function;

[0023] The real-time inversion is as follows: first, the real-time RTCM data stream is received from the receiver, and data format conversion, SNR data extraction and detrending items are performed in the software. The reflection height data is calculated using the satellite arc segment. Then, the water level curve in each window is fitted using the inverse modeling method according to the window size and node interval set by the user, and the result is output. The last node water level in each window is the real-time inversion 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.

[0024] Furthermore, step 4 includes the following steps:

[0025] S41 uses GPS / BDS precise point positioning, and uses the Kalman filter method to calculate the coordinates, carrier phase ambiguity, tropospheric zenith wet delay, and clock error as unknowns to estimate 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.

[0026] S42 determines the altitude, pressure and temperature data of four grid points near a GNSS station:

[0027]

[0028] Where: and represents the adjusted pressure and temperature of the ith grid point; and Indicates the pressure values ​​of the two levels closest to the station; and Indicates the closest two-level temperature value; and Represents the potential height of the two closest levels; Hs is the station elevation.

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

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

[0031] Beneficial effects of the present invention:

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

[0033] (1) High precision and high reliability: The displacement monitoring method of the present invention adopts advanced signal processing and error correction methods to solve the common problems of multipath, non-line-of-sight signals, gross errors, multiple gross errors, frequent cycle slips, etc. in complex engineering environments, ensuring 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 adopting evenly spaced node settings, the influence of insufficient observation values ​​or gross errors is effectively avoided, ensuring the accuracy of the fitting results, and making the Beidou continuous static monitoring accuracy stable at a high standard of 3 mm (RMS).

[0034] (2) High timeliness and real-time performance: This method focuses on the timeliness of monitoring results on the basis of high precision. When dealing with near-real-time water level inversion, the innovative 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 to meet the needs of emergency disaster response.

[0035] (3) Precise meteorological data processing and water vapor inversion: To address the problem of insufficient accuracy in BeiDou PWV (atmospheric precipitable water volume) inversion, a quality control method for zenith tropospheric delay (ZTD) calculation with precise single-point positioning was proposed, and a wet delay (ZWD) and PWV conversion coefficient model was established to make the water vapor inversion process more accurate. This solved the key problem of missing ground meteorological data, allowing effective water vapor monitoring through the BeiDou system even when the data is incomplete. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] The drawings constituting a part of this application are used to provide a further understanding of the present invention. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation of the present invention. In the drawings:

[0037] Figure 1 This is a schematic diagram of the five-segment weight determination 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 This is the post-inversion water level result of the demonstration point in Example 1;

[0040] Figure 4 This is the near real-time inversion water level result of the demonstration point in Example 1;

[0041] Figure 5 The PWV time series diagram inverted by four GNSS stations in the Hetao area of ​​Bayannur, Inner Mongolia in Example 1;

[0042] Figure 6 The ZTD comparison result based on Beidou B2b and post-precision ephemeris estimation in Example 1 (Bayannur Hetao JZ01);

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

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

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

[0046] Fig.10 This is the comparison between the PWV estimated by NCEP-GFS and the empirical GPT model and the sounding data in Example 1 (Bayannur JC04 station). DETAILED DESCRIPTION

[0047] Example 1

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

[0049] Step 1: Receive and process GNSS satellite signals based on a GNSS receiving antenna: extract direct signals; extract signal-to-noise ratio data, remove the direct signal part in the signal-to-noise ratio data, and only retain the reflected signal; extract GNSS tropospheric delay information; the GNSS receiving antenna is a geodetic measurement antenna.

[0050] Step 2: Monitor the displacement and deformation of water conservancy projects 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 judge the affected type of the signal, and select the observation variance calculation model accordingly; bring the calculated observation variance into the RTK processing to obtain coordinates, and adjust the weight of the quality indicator in the observation variance calculation model according to the posterior variance, so as to construct a comprehensive random model suitable for GNSS positioning in complex environments to calculate the coordinates of water conservancy projects and realize displacement and deformation monitoring; the specific operations of this embodiment are as follows:

[0051] S21 collects GNSS raw observation data of different complex scenarios, fits the functional relationship between observation variance and quality index, and selects key quality indexes. In this embodiment, the selected key quality indexes include altitude angle, carrier-to-noise ratio and PLD.

[0052] First, the original observation data of different complex scenes are collected, and the functional relationship between the observation variance and the quality index is fitted.

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

[0054]

[0055] Among them, a and b are empirical coefficients, obtained by fitting method, or empirically configured as a = 3 mm, b = 3 mm, θ represents the satellite altitude angle, i represents the satellite identifier, σ {·} is the standard deviation of the observations.

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

[0057]

[0058] Where C / N0 is the carrier-to-noise ratio in dB-Hz, and C is a constant obtained by fitting or empirical setting. For example, it can be set to 0.00224m for GPSL1 observations. 2 Hz.

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

[0060]

[0061] in is the phase observation variance, and the subscripts G and C represent GPS and BDS respectively. The corresponding pseudorange variance can be obtained by taking the ratio of the variance to the phase variance:

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

[0063]

[0064] S22 establishes template functions 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 altitude angle as a variable is:

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

[0067] +2.06×10 -5 ×θ 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×θ-5.51×10 -3 ×θ 2

[0071] +2.46×10 -5 ×θ 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×θ+6.96×10 -4 ×θ 2

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

[0076]

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

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

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

[0080]

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

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

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

[0084]

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

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

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

[0088]

[0089] S23 uses the difference between the measured value of the quality indicator and the template value to judge the affected type of the signal and determine the segmented weighting strategy to establish a segmented weighting model. The standardized difference data between the measured value of the quality indicator and the template value is substituted into the segmented weighting model. Based on the segment, the observation variance calculation model is selected and the variance is calculated.

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

[0091] Calculate the normalized carrier-to-noise ratio and PLD:

[0092]

[0093] The subscript mea represents the observed value, and the subscript nom represents the template value.

[0094] Substitute it into the following five-stage weighting model, and the strategy diagram is as follows Figure 1 As shown, where n 0 ,n 1 ,m 0 ,m 1 is an empirical constant, representing the quantile of the normal distribution, reflecting the abnormality of the indicator. 0 and m 0 Set to 1~2, and set n 1 and m 1 For settings 2 to 3, each segment model is obtained by fitting the measured data:

[0095] 1) Weighted model section 1: It is assumed that the observation value is not affected by multipath, diffraction, etc., and the altitude angle model is used:

[0096]

[0097] 2) Weighted model segment 2: When the standardized index falls into this range, it is considered that the signal is interfered by multipath and diffraction, and the elevation angle model is inaccurate, and the following combined model is used:

[0098]

[0099] 3) Fixed-weight model segment 3: When the standardized index falls into this interval, it is considered that the signal is mainly affected by cycle slips and unstable phase tracking, and the following comprehensive model is used:

[0100]

[0101] 4) Weighted model segment 4: When the standardized 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 point, it is considered that the signal is affected by strong multipath, cycle slips and tracking loop instability, and the following model is used:

[0104]

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

[0106] After obtaining the variance of the observation value, it is also necessary to propagate the variance of the non-differenced (original) observation value to the double-differenced observation value through the error propagation law, and then combine it with the function model to complete the RTK positioning. The specific process is:

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

[0108]

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

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

[0111]

[0112] In the formula, O is the undifference observation vector and E is the unit matrix. represents the variance of a monitoring station corresponding to a certain satellite, Represents the variance of the reference station corresponding to the same satellite. Based on the inter-station 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 vector and its variance covariance matrix can be expressed as:

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

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

[0115]

[0116] Where I is a vector whose elements are 1.

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

[0118]

[0119] In the formula, is the double difference operator, ρ is the satellite-to-ground distance, ε is the observation noise, and N is the integer ambiguity. The above are the double difference function model and double difference random model required for RTK positioning. With these, the coordinates of the mobile station can be obtained by following the steps of Kalman filtering, ambiguity fixation and verification.

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

[0121] In order to achieve better robustness, the Huber variance inflation function is used to adjust the observation weights based on the five-segment model according to the posterior residuals:

[0122]

[0123] Among them, γ ii is the variance inflation factor, c 0 is a predetermined empirical constant, such as 2.5 to 3.0, is the standardized residual.

[0124] Experimental verification:

[0125] In order to verify the performance of the proposed random model for RTK positioning, 1 hour of GPS+Beidou L1 / B1 data was collected in a tree-shaded environment. Static data collection was performed using a special multi-frequency multi-GNSS receiver that can output PLD measurements. 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., the five-segment fixed weight model + Huber) and the comprehensive random model (Compre, i.e., the five-segment fixed weight model) have higher fuzzy success rate and positioning accuracy than other random models.

[0126] In terms of fuzzy success rate, the comprehensive random model improved by 50.04%, 17.62% and 8.46% compared with the altitude angle, carrier-to-noise ratio and PLD random models, respectively. The fuzzy success rate of the robust comprehensive random model increased by 4.65% compared with the comprehensive random model. In particular, the fuzzy success rate of the altitude angle random model was 0%, which may be due to the presence of a large amount of occlusion and cycle slips in the shaded environment observation data. In this case, it is inappropriate to use only the altitude angle to determine the random model. In terms of positioning accuracy, when the ambiguity is successfully fixed, the results of all random models (except the altitude angle random model) obtain centimeter-level positioning accuracy in three directions, such as Figure 2 and Table 1. The model combines three data quality indicators, namely, altitude angle, carrier-to-noise ratio and PLD. The weighting method is adjusted by the difference between the measured value of each indicator and the template value. At the same time, the residual information is used to adjust the pre-verification weight, which can effectively solve the problem of failure of weighting of a certain indicator in a complex environment and ensure that a reasonable weight can be given to the observation value during the entire solution period. By comparing and analyzing the results of RTK positioning with different random models, the model can effectively improve the success rate of ambiguity fixation while ensuring positioning accuracy.

[0127] Table 1 Statistical comparison of RTK positioning effects of different random models in shaded environments

[0128]

[0129]

[0130] Step 3: Monitor the water level based on the reflection signal extracted in step 1: Use the inverse modeling method to invert the water level, and use the B-spline curve to fit the water level change.

[0131] This embodiment reads the pre-processed GNSS signal-to-noise ratio (SNR) data and selects post-inversion or near-real-time inversion methods according to user needs. The common method for Beidou-R water level inversion is to analyze the spectral characteristics of SNR, the results of which are greatly affected by noise, and the time distribution of the output results is uneven. The present invention adopts an inverse modeling method for water level inversion, which uses B-spline curves to fit water level changes, which can significantly improve the accuracy and stability of the inversion results.

[0132] ① Post-inversion: The LSP spectrum analysis method reads the SNR sequence of a single satellite, performs spectrum analysis to calculate the peak frequency, and then performs data quality control. The inverse modeling method does not require satellite-by-satellite analysis, but uses B-spline curves to model the water level that changes over time as a smooth and continuous function.

[0133] The specific operations in this embodiment are as follows:

[0134] After receiving the GNSS satellite signal data, the signal-to-noise ratio data contained therein is extracted, and the direct signal part in the signal-to-noise ratio data is removed using a low-order polynomial, leaving only the reflected signal part, and establishing a mathematical model between it and the reflection height. First, the signal-to-noise ratio data after the trend term is removed is modeled as:

[0135]

[0136] In the formula, δSNR is the signal-to-noise ratio data after removing the direct signal; λ is the wavelength of the satellite signal; h is the reflection height; θ is the satellite altitude angle; k is the wave number; s is the roughness parameter of the reflection surface; In the formula, C 1 With C 2 In-phase and out-of-phase components are used to replace the amplitude A and phase The conversion relationship is:

[0137]

[0138] Before applying the inverse modeling method, the δSNR data should be preliminarily analyzed using the LSP spectrum analysis method to determine the C 1 , C 2 ,h,s 2 Initial values ​​of each parameter.

[0139] Before using the B-spline curve to fit the water level change curve, it is necessary to select n curve nodes at equal time intervals, and 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 parameter, where x is the unknown parameter, that is, C 1 , C 2 , P i , roughness parameter s 2 , node P i , the function can be expressed as:

[0140]

[0141] In the formula, i represents any integer between 1 and n, indicating the ordinal number of the node; δSNR i is the observed data, that is, the signal-to-noise ratio data after removing the direct signal. The initial reflection height is calculated by 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 results of LSP, and the undetermined parameters are estimated using the nonlinear least squares method. Through repeated iterations, y i The optimal parameter solution is obtained when the residual sum of squares of (x) is minimized, namely:

[0142]

[0143] Where n is the number of δSNR data. After the iteration is completed, the B-spline interpolation can be performed using the estimated value h at the node to obtain a uniform water surface height change sequence. In the present invention, the function of the B-spline curve of the water (tide) level change can be expressed as:

[0144]

[0145] Where h(t) is the B-spline curve of the water surface changing with time t; is a p-order B-spline basis function, which can be expressed as:

[0146]

[0147] Where u is the parameter in the node space, that is, u∈[t 0 ,t n ]; t i represents the time window [t 0 ,t n ]; p represents the order of the spline curve. Since the water surface height changes continuously with time, the order of the B-spline basis function is selected as the second order in the present invention.

[0148] The present invention intends to utilize the above method to improve the inversion accuracy of BeiDou-R post-inversion and make breakthroughs in high-precision BeiDou-R water level post-inversion technology.

[0149] ② Real-time inversion: First, receive the real-time RTCM data stream from the receiver, perform data format conversion, SNR data extraction and detrending items in the software, and use satellite arc segments to calculate reflection height and other data preprocessing operations. 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 results are output. The last node water level in each window is the real-time inversion 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 nodes.

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

[0151] The JC01 station is located at 40.66°N and 107.26°E. The observation environment of JC01 is good, with an open monitoring range. There are no tall buildings and trees around it, so it can capture abundant water surface reflection signals, ensuring the adequacy of the data source. Considering the actual situation, it is necessary to artificially limit the range of satellite elevation angle and satellite azimuth angle. According to the Fresnel reflection zone of the JC01 station under the map perspective and the distribution of colored strips, it can be determined that the reflection signal from the water area is received in the azimuth range of 135°~220° and the elevation range of 7°~25°.

[0152] 1. Post-event inversion demonstration effect

[0153] The post-inversion results are as follows Figure 3 As shown in the figure, the horizontal axis represents time, specifically expressed in the form of annual accumulated days, and the vertical axis represents the change of the water surface. The blue curve is the measured water level value with a sampling rate of 5 minutes at the tide gauge station, and the red curve is the water level value inverted using the post-inversion algorithm. It can be found that the post-inversion experimental results and the measured data of the tide gauge station are in good consistency. The water level inversion error of this station is 2.3cm, and the correlation coefficient is 0.99, which meets the project requirements and can effectively reflect the water level change trend of this water area.

[0154] (II) Near real-time inversion demonstration effect

[0155] Figure 4 This is the inversion result of the near-real-time algorithm. The accuracy of the inversion result of the near-real-time algorithm is 4.3 cm, and the correlation coefficient is 0.97. Due to the short length of data used by the near-real-time algorithm and the small area of ​​water in this area, there is a certain deviation in the inversion result, but it can still achieve a relatively accurate match with the measured data of the tide gauge station.

[0156] Step 4: Based on the direct signal extracted in step 1, the GNSS tropospheric delay is estimated using precise point positioning technology, and then the tropospheric delay information is inverted into atmospheric precipitable water in real time.

[0157] First, based on the precise single-point positioning technology, the coordinates, carrier phase ambiguity, tropospheric zenith wet delay, clock error, etc. are used as unknowns to calculate the tropospheric delay ZTD using the Kalman filter method. 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] Where, 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, m gN / gE is the corresponding projection function. In precise point positioning, the tropospheric delay correction model is used to correct the tropospheric dry delay, and the random walk process is used to estimate the zenith tropospheric wet delay. The zenith tropospheric delay can be decomposed into two parts: the zenith static delay ZHD and the zenith wet delay ZWD. Among them, ZHD can be accurately obtained using the empirical model supplemented by ground pressure data. The specific formula is:

[0160]

[0161] Where P s : surface air pressure, Station latitude, Hs: station elevation.

[0162] ZTD-ZHD=ZWD

[0163] Then, ZWD can be directly converted to PWV using the conversion factor:

[0164] PWV=Π×ZWD

[0165] Where Π is the dimensionless conversion coefficient, which is related to the weighted average temperature T m related.

[0166]

[0167] Among them, the physical constant k 3 and k′ 2 They are 3.7546×10 5 K 2 / hPa and 22.9721K / hPa. m It is defined as the integral of humidity and 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. Vertical humidity and temperature profiles can be extracted from sounding and reanalysis data. However, in most cases, atmospheric vertical profile meteorological data are difficult to obtain in real time. m Estimation models are widely used in the field of meteorology.

[0170] As mentioned above, the surface pressure P is required to convert the tropospheric delay ZTD obtained by GNSS inversion into atmospheric precipitable water PWV. s and weighted temperature T m Two meteorological data. When the GNSS station lacks meteorological sensors, ZTD-PWV inversion can usually be performed with the help of empirical meteorological parameter models, but the accuracy is relatively low. The present invention proposes a PWV real-time inversion technology based on numerical forecast model products.

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

[0172]

[0173] Where: and represents the adjusted pressure and temperature of the ith grid point; and Indicates the pressure values ​​of the two levels closest to the station; and Indicates the closest two-level temperature value; and It represents the potential height of the two closest levels. Once the pressure and temperature data of the four nearby grid points are obtained, bilinear interpolation is performed to determine the ground pressure and weighted temperature within the station, and then the ZTD data observed by the GNSS station can be converted into PWV.

[0174] Four demonstration sites (JZ01, JC01, JCC02, and JC04) have been set up in the Hetao area of ​​Bayannur, Inner Mongolia, and four Beidou / GNSS receivers have been installed. Figure 5 The PWV time series graphs inverted from four GNSS stations are shown. The region is dry and rainy with high altitude, and the atmospheric water vapor content is relatively low.

[0175] Based on 10 days of data collected by the GNSS base station in the Hetao area of ​​Bayannur, Inner Mongolia (October 29 to November 7, 2024), the ZTD estimated by the post-precision ephemeris PPP was compared with the ZTD product estimated by the B2b output of the receiver. Figure 6) shows that the RMS value of the ZTD product estimated based on B2b is 8.1mm, which proves that the ZTD data accuracy of the real-time PPP estimated based on Beidou B2b product can meet the needs of this project.

[0176] In addition, the atmospheric precipitable water PWV data inverted by the Bayannur Linhe sounding station was used to evaluate the PWV accuracy inverted by four GNSS demonstration stations (JZ01, JC01, JC02, and JC04). The sounding station is up to 12 kilometers away from the GNSS station, which is very suitable for evaluating the water vapor accuracy of GNSS inversion. Since the four GNSS demonstration stations do not have co-located meteorological observation data, the project uses the meteorological parameter forecast products of the NCEP-GFS numerical forecast model to achieve ZTD-PWV conversion. As follows Figure 7 As shown, at the JZ01 station, the RMS value of the PWV product converted based on the numerical forecast 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, with the deviation exceeding 3 mm at some moments. Figure 8-10 It is shown that at the JC01, JC02 and JC04 stations, the RMS values ​​of the PWV products converted based on the numerical forecast model NCEP-GFS reached 1.08mm, 1.16mm and 1.50mm respectively, and the corresponding RMS values ​​of the PWV products converted based on the empirical meteorological model GPT2w were 2.22mm, 2.37mm and 2.81mm respectively. It can be seen that the real-time PWV estimation accuracy of the method of the present invention based on the Beidou B2b ephemeris product and the numerical forecast model NCEP-GFS fully meets the 3mm accuracy requirement. Although the overall accuracy of the PWV converted based on the empirical meteorological model GPT2w can also meet the 3mm accuracy requirement, the deviation exceeds 3mm at some moments.

[0177] Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

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 based on a GNSS receiving antenna and process them: extract direct signals; extract signal-to-noise ratio data, remove the direct signal part in the signal-to-noise ratio data, and only retain the reflected signal; Step 2: Monitor the displacement and deformation of water conservancy projects 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 affected type of the signal, and select the observation variance calculation model accordingly; bring the calculated observation variance into the RTK processing to obtain the coordinates, and adjust the weight of the quality indicator in the observation variance calculation model according to the posterior variance, so as to construct a comprehensive random model suitable for GNSS positioning in complex environments to calculate the coordinates of water conservancy projects and realize displacement and deformation monitoring; Step 3: Monitor the water level based on the reflection signal extracted in step 1: use the reverse modeling method to invert the water level, 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 is characterized in that: Step 2 includes the following operations: S21 collects GNSS raw observation data of 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 judge the affected type of the signal and determine the segment weighting strategy to establish a segment weighting model, substitutes the standardized difference data between the measured value of the quality indicator and the template value into the segment weighting model, selects the observation variance calculation model based on the segment, and calculates the variance; S24 substitutes the variance calculated in S23 into the positioning equation to obtain the coordinates to complete the RTK positioning; S25 uses the posterior variance to adjust the weights and back-substitutes step S24 to update the coordinates until the coordinates converge.

4. The GNSS-based water conservancy multi-parameter monitoring method according to claim 3 is 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 the function model to complete the RTK positioning.

5. The GNSS-based water conservancy multi-parameter monitoring method according to claim 1 is characterized in that: In step 3: water level monitoring including post-inversion or real-time inversion; The post-inversion: 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 B-spline curves to model the water level height that changes with time as a smooth and continuous function; The real-time inversion is as follows: first, real-time RTCM data stream is received from the receiver, data format conversion, SNR data extraction and detrending items are performed in the software, and reflection height data is calculated using satellite arc segments. 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 results are output. At the same time, the node records in the window are saved, and finally a long water level time series is constructed through all nodes.

6. The GNSS-based water conservancy multi-parameter monitoring method according to claim 1, characterized in that: Step 4 includes the following steps: S41 uses GPS / BDS precise point positioning, and uses the Kalman filter method to calculate the coordinates, carrier phase ambiguity, tropospheric zenith wet delay, and clock error as unknowns to estimate 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 ith grid point; and Indicates the pressure values ​​of the two levels closest to the station; and Indicates the closest two-level temperature value; and It represents the potential height of the two closest levels; Hs is the height of the measuring station; After obtaining the pressure and temperature data of four grid points near the GNSS station, bilinear interpolation is performed to determine the ground pressure and weighted temperature within the station; S43 uses the ground 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