Disclosure of Invention
The invention provides a method for monitoring the near-real-time millimeter-scale deformation by using a long base line, which improves the processing precision of the near-real-time long base line from the aspects of posterior residual editing iteration, sliding window processing, multi-base station combined processing and ambiguity fixing strategy processing; in order to improve the high-efficiency processing of massive users, in addition to distributed processing, optimization is also performed in the aspects of non-differential mode parallel processing, common reference station processing and normal equation superposition.
The invention solves the following technical problems:
1) the deployment quantity of the domestic reference stations can be effectively reduced, and the resources of the existing reference stations are fully utilized.
2) And a high-efficiency data processing strategy is provided, which comprises a short arc equation combination and non-differential distributed processing mode, so that the computing resources can be effectively reduced, near-real-time deformation monitoring service across the country is provided, and the service cost is reduced.
The technical scheme adopted by the invention is as follows:
a Multi-GNSS long baseline near real-time deformation monitoring method based on multiple base stations is characterized by comprising the following steps:
according to the coordinates of the monitoring points, a reference station analysis is combined, a reference station selection strategy is utilized, and a reference station list participating in resolving is given;
dividing subtasks according to monitoring point data participating in resolving, and executing cloud platform distributed processing;
the cloud platform distributed processing utilizes an ultra-fast Multi-GNSS satellite orbit product to reduce the influence of orbit errors on near-real-time long baseline solution of the Multi-GNSS.
Further, the reference station selection strategy specifically includes the following steps:
step 1, screening a reference station within a distance of 100km according to the coordinates of the monitoring points;
step 2, screening the reference station with the absolute coordinate point repetition rate less than 3mm horizontally and 6mm vertically according to the coordinate stability information of the reference station;
step 3, screening the reference station meeting the admission standard according to the data quality information of the reference station;
step 4, dividing the candidate reference stations into six areas according to the azimuth angles, and selecting the reference station closest to the monitoring point in each area;
and 5, giving a reference station list participating in resolving.
Further, the candidate reference stations are divided into six zones according to the following method:
zone 1, 0-60 °;
zone 2, 60-120 °;
zone 3, 120-180 °;
zone 4, 180-240 °;
zone 5, 240-300 °;
zone 6, 300-360 deg.
Further, the Multi-GNSS near real-time long baseline solution specifically includes the following steps:
performing non-differential PPP (Point-to-Point protocol) processing on a reference station and a monitoring point by using an ultra-fast Multi-GNSS satellite orbit product, performing posterior residual editing iteration, removing abnormal data, and generating normal equation information of a current arc section, wherein m-1 normal equations are arranged in front of the normal equation of the current arc section;
combining the normal equation of the current arc section and the previous m-1 normal equations into a long arc normal equation by utilizing sliding window batch processing;
adopting an ambiguity fixing strategy, adopting independent base line pairs for a plurality of base stations, adopting independent satellite pairs for satellite pairs on the base line pairs, and screening the ambiguity fixing probability of the wide lane/narrow lane of the satellite pairs on the base line pairs;
and (3) carrying out ambiguity fixing by adopting a wide lane and narrow lane combined mode, adding ambiguity fixing information into a long arc equation to update the long arc equation, and repeatedly carrying out ambiguity screening and ambiguity fixing treatment on the wide lane/narrow lane.
Further, the parameters to be estimated for non-differential PPP processing of the reference station and the monitoring point include: satellite clock error, receiver clock error, coordinates, zenith troposphere delay, inter-satellite system delay, and inter-satellite frequency delay.
Furthermore, the ambiguity fixing strategy adopts double-difference ambiguity constraint and adds a long arc equation to realize the joint processing of ambiguity fixing and multiple base stations.
Further, the double-difference ambiguity constraint is equivalent to adding constraint conditions to four non-difference ambiguities to define a virtual observation value V of the double-difference ambiguityB：
Wherein,double-difference ambiguity for fixed ionospheric-free combinations, PBThe weight value of the virtual observed value of the double-difference ambiguity is shown, D is a projection matrix of the non-difference ambiguity, and X is a parameter vector to be estimated;
wherein,andrepresenting two satellites, the two satellites are provided with a single satellite,which represents two of the receivers, is,is composed ofStation andthe degree of non-differential ambiguity between satellites,is composed ofStation andthe degree of non-differential ambiguity between satellites,is composed ofStation andthe degree of non-differential ambiguity between satellites,is composed ofStation andnon-differential ambiguity between satellites;
based on the constraints of equation (11) and equation (12), the form of the least squares solution is rewritten as equation (13), implementing a transition from non-differential to double-differential baseline processing:
where A is a coefficient matrix, P is a weight matrix, L is an observed value minus a model calculated value, NX ═ W, N is a normal equation, and X is a weight matrixnewFor X, N with double interpolation ambiguity constraintnewFor adding double-interpolation ambiguity constraint method N, WnewW with double-interpolated ambiguity constraints added.
Further, if double-difference ambiguity constraints of multiple baselines are applied, the multi-reference station combined fusion processing can be realized.
Further, the degree of blur needs to reach 95% in the degree of blur fixing process.
Further, the cloud platform distributed processing specifically includes the following steps:
setting parameter and non-difference PPP processing subtask division of n monitoring points, setting parameter and non-difference PPP processing subtask division of q reference stations, and setting parameter and non-difference PPP processing subtask division of p cloud servers;
processing tasks of the monitoring points and the reference station are divided to different cloud processing servers for parallel processing through a non-differential PPP processing subtask dividing module;
recovering and combining short arc equations of the monitoring points and the reference station, and dividing the equations to the cloud server for parallel processing;
and carrying out ambiguity fixing processing.
The invention has the following beneficial effects:
1) a near real-time long baseline resolving algorithm is provided, the problem of low near real-time processing precision of a long baseline is solved, and the method has the core innovation points that: the method has the advantages that the influence of orbit errors on baseline resolution is reduced by using an ultra-fast Multi-GNSS (Multi-Global navigation satellite system) satellite orbit product, baseline processing precision and stability can be improved by using multiple residual iterations, near-real-time high-precision baseline resolution can be realized by using a sliding window batch processing mode, baseline resolution precision and stability can be improved by using a Multi-reference station combined processing mode, and long baseline ambiguity fixing rate can be improved by using an innovative ambiguity processing strategy, so that the precision is further improved.
2) Based on a non-differential processing mode, double-differential ambiguity constraint is applied in a virtual observation value mode, a multi-reference station processing technology can be flexibly realized, distributed subtask division can be performed in a finer granularity mode, the distributed computing advantages can be fully utilized, the subtask division is flexible, and the data processing efficiency is improved; the method for combining the short arc method equation into the long arc method equation is provided, so that redundant processing in batch processing of the sliding window can be avoided, and the resolving efficiency is further improved.
Detailed Description
In order to improve the stability and the precision of long baseline solution, the invention adopts a multi-reference station joint processing mode. Due to data interruption of a single reference station, periodic change of satellite observation geometry, different stability of the reference station, different data quality of the reference station, different geometric distance between the reference station and a monitoring point and different geometric distribution of the reference station relative to the monitoring point, the influence of selecting an improper reference station on the monitoring precision is large. The invention provides a selection strategy of multiple reference stations around a monitoring point, and the reference stations participating in resolving are automatically selected according to a certain threshold value, so that the stability and the precision of long-baseline resolving are preliminarily ensured.
In order to improve the long baseline resolving precision, the method utilizes an ultra-fast Multi-GNSS satellite orbit product to reduce the influence of orbit errors on baseline resolving; editing and iterating by utilizing a plurality of posterior residual errors, and putting forward the influence of data abnormity on baseline calculation as much as possible; by utilizing a sliding window batch processing mode, near-real-time high-precision baseline resolution can be realized; by utilizing a multi-reference station combined processing mode, the accuracy and stability of baseline calculation can be improved; and an innovative ambiguity processing strategy is adopted, so that the baseline resolving precision is further improved.
In order to improve the data processing efficiency, the method and the system utilize a method equation superposition technology, fully utilize the prior method equation information, and avoid the resolving repeatability; processing by utilizing a common reference station, and avoiding repeated processing of the same reference station; by utilizing the non-differential processing mode, the distributed computing advantages can be fully utilized, and the data processing efficiency is improved.
The invention is further illustrated below with reference to the figures and examples.
Fig. 1 is a general flow chart of Multi-base-station Multi-GNSS near real-time long baseline deformation monitoring of the present invention, which includes the following steps: according to the coordinates of the monitoring points, a reference station analysis is combined, a reference station selection strategy is utilized, and a reference station list participating in resolving is given; and performing subtask division according to the monitoring point data participating in resolving, executing cloud platform distributed processing, and improving processing efficiency.
Fig. 2 is a flow chart of reference station selection of the present invention, comprising the steps of: firstly, screening a reference station within a distance of 100km according to the coordinates of the monitoring points; secondly, screening the reference station with the absolute coordinate point repetition rate less than 3mm horizontally and 6mm vertically according to the coordinate stability information of the reference station; thirdly, screening the reference station meeting the admission standard according to the data quality information of the reference station; fourthly, dividing the candidate reference stations into six zones according to the azimuth angle, and selecting the reference station closest to the monitoring point in each zone as shown in figure 3 (zone 1, 0-60 degrees, zone 2, 60-120 degrees, zone 3, 120-180 degrees, zone 4, 180-240 degrees, zone 5, 240-300 degrees, zone 6, 300-360 degrees); and finally, giving a reference station list participating in the resolving.
The preferred embodiment of the invention: the monitoring precision is 2mm horizontally and 4mm vertically at the baseline of 70 kilometers, and the calculation updating frequency is 1 hour.
FIG. 4 is a flow chart of the Multi-GNSS near real-time long baseline solution of the present invention. Firstly, carrying out non-differential PPP processing on a reference station and a monitoring point by using an ultra-fast Multi-GNSS satellite orbit product, wherein the parameters are estimated as follows: satellite clock error, receiver clock error, coordinates, zenith troposphere delay, satellite inter-system delay (ISB) and satellite inter-frequency delay (IFB), and carrying out posterior residual editing iteration to remove abnormal data and generate normal equation information (formula, detailed description of matrix and detailed description of each important parameter) of the current arc segment:
wherein,is a GPS ionosphere-free carrier-phase observation,is an ionospheric-free carrier-phase observation for GLONASS (russian GLONASS satellite navigation system),for BDS (China Beidou satellite navigation System) ionosphere-free carrier phase observed value, G represents GPS, R represents GLONASS, C represents BDS, i represents a receiver, j represents a satellite number, rho is the geometric distance from the satellite to the receiver, and is the satellite position [ x ]jyjzj]TWith receiver coordinate [ x ]iyizi]TM is a troposphere model mapAs a function, T is the zenith tropospheric delay, c is the speed of light, λ is the wavelength,is the receiver GPS clock time-out, dtjIs the clock error of the satellite or the like,is the signal propagation time, ISB, between satellite j and receiver iR-GAnd ISBC-GIs the intersystem deviation,. ljIs the frequency label of the GLONASS satellite, IFBiIs the carrier phase inter-frequency offset of the receiver, biAnd bjRespectively, the fractional deviation of the phase of the receiver and the satellite, and B is an ambiguity parameter.
The invention adopts ultra-fast Multi-GNSS satellite orbit product pair [ x ] when non-differential PPP processing is carried outjyjzj]TFixing; fixing the ambiguity by using a double-difference constraint modeiAnd bjUnifying into the non-differential ambiguity parameter, the later mentioned ambiguity contains the hardware phase fractional deviation term. Parameters that need to be estimated for non-differential PPP processing include, the receiver coordinate [ x ]iyizi]TClock error of receiverSatellite clock difference dtkInter-system bias ISBR-GAnd ISBC-GInter-frequency offset IFBiZenith tropospheric delay T, and ambiguity parameter B.
Equation (2) is an error equation, a is a coefficient matrix, X is a parameter vector to be estimated, V is a residual vector,is a virtual observed value of the prior value, L is the observed value minus the calculated value of the model, P is the weight matrix,X0andis the prior value and value array of the parameter vector, the corresponding normal equation is, NX ═ W,n is the normal equation. In order to improve the precision and stability of the solution, the posterior residual error needs to be iteratively edited to remove abnormal data, the basic principle is as follows, the variation of the posterior residual error of the kth epoch relative to the residual error of the k-1 epoch can be expressed as,
Δvk＝v(k)-v(k-1) (3)
the residual mean and mean square error of each epoch within successive observation arcs can be expressed as,
after the maximum residual variation max (delta v) is removed, the mean value and mean square deviation of the residual variation of each epoch are calculated again,is the number of the observed values,
if the following conditions are satisfied,
＞k0new
max(Δv)＞J0(6)
wherein k is0、J0If the current calendar is the set threshold value, cycle slip and gross error of epoch k can occur, and the (k + 1) th calendar is continuously checkedAnd (5) Yuan. If the k +1 th epoch also exceeds the limit and the residual variance also exceeds the limit, the k epoch is considered as gross error. Otherwise, the cycle slip is considered to have occurred at the kth epoch. And taking the kth epoch as the starting epoch of the new residual arc segment, and judging again.
The method comprises the steps of obtaining a clean normal equation N through multiple times of non-differential PPP (precision Point Positioning) processing and posterior residual editing, performing long arc section solution in order to enhance the stability of the solution, avoiding redundancy processing in order to improve the solution efficiency, performing batch processing by using a sliding window, merging the previous m-1 normal equations with the current normal equation, and combining the boundary constraint conditions of the two adjacent methods as shown in a formula (7), wherein V isk，k+1Is a virtual observation of a boundary constraint, XkAnd Xk+1Is the estimated parameter correction value of two arc segments, Pk，k+1Is the weight value of the virtual observation, Lk，k+1Is the observed value minus the model calculated value. The adjacent method equation combination is actually the adjustment problem with additional virtual observation value constraint, and the minimum variance estimation is adopted, so that the performance function is minimum, and the equation (8) is satisfied.
Vk，k+1＝Φ(tk+1，tk)Xk-Xk+1，Pk，k+1-Lk，k+1(7)
Wherein, tkIs the time of epoch k, tk+1Is the time of epoch k +1, Φ (t)k+1，tk) Is a state transition matrix, VkAnd PkIs k epoch residual vector and weight value, Vk+1And Pk+1Is the k +1 epoch residual vector and weight value, min represents the minimum value expressed by this formula.
According toIt is possible to obtain,
wherein N iskIs the normal equation for k epochs, Nk+1Is the normal equation for the k +1 epoch. The k epoch method equation, the k +1 epoch method equation and constraint information among epochs are considered, the merged method equation is shown in a formula (10), the merging method can fully utilize the stability of long arc resolving and can avoid redundant non-differential PPP processing to improve the processing efficiency.
In order to avoid the fault of a single reference station, the influence of the satellite geometry on the result is weakened, the stability of the result is enhanced, a multi-reference station mode is adopted, and the joint processing effect of ambiguity fixing and the multi-reference station is realized in a method equation added by double-difference ambiguity constraint. The double-difference ambiguity problem is equivalent to defining a virtual observation value V of double-difference ambiguity by adding proper constraint conditions to four non-difference ambiguitiesB，
Wherein,double-difference ambiguity for fixed ionospheric-free combinations, PBIs the weight of the virtual observations, D is the projection matrix of the non-differential ambiguities,
wherein,andrepresenting two satellites, the two satellites are provided with a single satellite,which represents two of the receivers, is,is composed ofStation andnon-differential ambiguity between satellites. In consideration of the above constraint, the form of the least square solution can be rewritten as formula (13), the conversion from non-difference to double-difference baseline processing is realized, if double-difference ambiguity constraint of a plurality of baselines is applied, the fusion processing mode of the multi-reference station can be realized, and the reliability of the processing result is enhanced.
Wherein, Xnew、Nnew、WnewX, N, W with double-interpolated ambiguity constraints.
Compared with other double-difference ambiguity fixing strategies, the fixed ambiguity can be eliminated from a normal equation by using the algorithm of the invention, and the calculation efficiency can be greatly improved.
The inventive ambiguity fixing strategy is adopted, the independent base line pair is adopted by the multi-base station, the independent satellite pair is also adopted by the satellite pair on the base line pair, the independent satellite pair is selected by adopting the similar strategy of selecting the independent base line, and the ambiguity fixing probability of the Wide Lane (WL)/Narrow Lane (NL) of the satellite pair on the base line pair is screened.
Carrying out ambiguity fixing by adopting a wide lane and narrow lane combined mode, and firstly fixingWidth-fixed lane ambiguity Bw，B1Is the ambiguity at the L1 frequency, B2Ambiguity at L2 frequency, wide-lane double-difference ambiguityIs composed of four non-differential ambiguities, as shown in formula (15), and wide lane double-differential ambiguity variance valueThe error propagation equation (16) yields. The fixed probability Prob of double-difference wide-lane ambiguity can be obtained according to the formula (17)wWhere B is the ambiguity parameter, σ is the estimated variance value,is the nearest integer value of B.
Bw＝B1-B2(14)
By using the formulas (19) and (20), the estimated value of the narrow-lane double-difference ambiguity can be obtainedAnd estimating the varianceWhereinAndthe fixed probability Prob of the narrow lane double-difference ambiguity can be obtained by a formula (17) through the ionosphere-free combined double-difference ambiguity value and the estimated variance valuen，
Probw＞Prob0，Probn＞Prob0(21)
Ambiguity fixing for satellite pairs passing the screening formula (21), where Prob0And if the ambiguity fixed probability threshold is the ambiguity fixed probability threshold, adding the ambiguity fixed information into the normal equation to update the normal equation, and repeatedly performing ambiguity screening and fixing until 95% of ambiguities are fixed.
FIG. 5 is a schematic diagram of the distributed solution principle of the present invention. As shown in FIG. 5, monitor points Mon participating in a certain process are set1，Mon2，...，MonnA total of n reference stations Ref participating in the process1，Ref2，...，RefqQ total cloud server Ser participating in processing1，Ser2，...，SerPThe invention adopts a non-differential PPP processing mode, so that the processing tasks of the monitoring points and the reference station can be divided to different cloud processing servers for parallel processing through a non-differential PPP processing subtask dividing module, short arc equations of the monitoring points and the reference station are combined and can be divided to the cloud processing servers for parallel processing, and double differences can be applied according to the equations (11) to (13) at the stage of ambiguity fixingAmbiguity constraint realizes multi-reference station joint coordinate calculation of a single monitoring point, and by the above strategies, distributed calculation advantages can be fully utilized, and the method is more suitable for scenes of mass data processing.
The core processing code of the invention utilizes C/C + +, and the automated control utilizes Python language.
Although the present invention has been described with reference to the preferred embodiments, it is not intended to limit the present invention, and those skilled in the art can make variations and modifications of the present invention without departing from the spirit and scope of the present invention by using the methods and technical contents disclosed above.