Positioning solution and trajectory correction method of geological disaster monitoring terminal

Abstract

The present application relates to the field of radio navigation positioning technology, and discloses a positioning solution and trajectory correction method for a geological disaster monitoring terminal, comprising: collecting original radio observation data and generating relative motion trajectory according to Doppler frequency shift cumulative integration, subtracting Doppler frequency shift from geometric prediction value to calculate Doppler transient residual error, estimating receiving position space deflection vector through smoothing filtering and shifting search axis center to construct a deflection-removed search cylinder, determining double-difference carrier phase integer ambiguity inside the cylinder to correct spatial coordinates, and resetting integration starting point and widening search radius by using double-difference pseudo-range when signal interruption occurs, wherein the present application eliminates phase lag deviation caused by forest refraction, prevents integer ambiguity solution divergence caused by search space deflection, and ensures smooth output of monitoring terminal positioning trajectory under complex shielding.

Description

Technical Field

[0001] This invention belongs to the field of radio navigation and positioning technology, and in particular relates to a method for positioning calculation and trajectory correction of a geological disaster monitoring terminal. Background Technology

[0002] High-precision positioning technology using Global Navigation Satellite Systems (GNSS) has been widely applied to terrain deformation monitoring. Positioning terminals typically calculate their position based on carrier phase observations and double-difference pseudorange observations, generating a three-dimensional displacement trajectory. To narrow the search range for carrier phase integer ambiguity, conventional methods also extract Doppler frequency shift observations, eliminate tropospheric and ionospheric delays through epoch-level difference analysis, and then integrate along the time axis to obtain the relative motion trajectory of the receiver phase center. This relative motion trajectory can serve as a spatial geometric constraint to limit the search range of candidate solutions for integer ambiguity. When the monitoring terminal is deployed in a densely forested canyon environment... The swaying vegetation canopy can block radio signals, causing the receiver tracking loop to lose lock and generating frequent cycle slips. Especially when the leaves are shaking and wet, the carrier signal passing through the leaves will exhibit phase lag, resulting in a phase lag bias in the observations that has a non-zero mean and drifts over time. This bias cannot be treated as ordinary Gaussian white noise. In this case, the center of the relative motion trajectory obtained by Doppler integration will deviate from the true carrier coordinate center, and the original symmetric geometric constraints will also be biased. The true integer ambiguity solution may fall outside the search boundary, leading to grid search failure.

[0003] Existing correction methods typically reduce signal interference by adjusting the observation noise covariance of the Kalman filter or by weighting the observations according to the satellite elevation angle. However, when there is a non-zero mean bias in the observations, these methods are insufficient to offset the resulting systematic spatial skew. If the bias is only covered by expanding the search boundary, the number of integer ambiguity candidate combinations will increase exponentially, which not only increases search time but also easily leads to incorrect fixation and trajectory jumps, making it difficult to meet the continuous observation requirements of high-precision deformation monitoring. Furthermore, in addition to the inherent limitations of hardware configuration and site selection, existing technologies also face bottlenecks at the monitoring algorithm level. Existing processing schemes mostly rely on preset geometric models, which have limited adaptability to sudden interference under complex operating conditions. For example, Chinese invention patent CN115235330B discloses a deformation monitoring data processing method and a deformation monitoring device. It achieves physical shock resistance through integrated devices, but its data calculation logic heavily relies on a specific double-difference observation model and LAMBDA algorithm. In complex field environments, once non-Gaussian noise interference or instantaneous interruption of satellite signals is encountered, this method is difficult to achieve adaptive adjustment of dynamic model parameters, often leading to failure of integer ambiguity resolution or a sharp drop in convergence speed, which in turn causes intermittent breakpoints in monitoring data. It is difficult to meet the high reliability early warning requirements of high-risk landslides at critical moments. This invention aims to improve the dynamic adaptive capability of the data processing algorithm and improve the accuracy and stability of deformation calculation in dynamic environments.

[0004] Therefore, the technical problem to be solved by this invention is how to correct the spatial skew of the Doppler relative trajectory center in response to the non-zero mean phase lag bias caused by vegetation refraction, fix the integer ambiguity of the double-difference carrier phase within a narrow search boundary, and maintain trajectory continuity when Doppler observation is interrupted. Summary of the Invention

[0005] To address the problems in the background art, the technical solution of the present invention is as follows: A method for positioning calculation and trajectory correction of a geological disaster monitoring terminal, comprising: Step S1: Obtain the raw radio observation data of the positioning terminal at the current epoch. The raw radio observation data includes dual-difference pseudorange, carrier phase, Doppler frequency shift, and carrier-to-noise ratio. Step S2: Extract a set of high-quality satellites with a carrier-to-noise ratio greater than a preset threshold. Use the Doppler frequency shift of the high-quality satellite set between epochs to determine the three-dimensional velocity vector and accumulate the integral to generate the relative motion trajectory. Step S3: Subtract the actual Doppler frequency shift from the geometric prediction value calculated based on the relative motion trajectory to determine the Doppler transient residual of each high-quality satellite at the epoch; Step S4: Project each Doppler transient residual onto the receiver coordinate system and smooth and filter it within a pre-set sliding time-domain window to estimate the spatial skew vector of the receiver position caused by the refraction of the tree canopy. Step S5: Calculate the vector sum of the relative motion trajectory and the spatial skew vector to determine the three-dimensional search center. Using this as the central axis and the standard deviation of the error as the search radius, construct a deskew search cylinder. Step S6: Determine the integer ambiguity of the double-difference carrier phase within the deskew search cylinder to solve the absolute position coordinates and output the deformation trajectory. At the same time, feed back the absolute position coordinates to clear the cumulative drift of the relative motion trajectory. If the Doppler frequency shift of the high-quality satellite set is continuously interrupted for 3 epochs, use the double-difference pseudorange to solve the absolute spatial position as the integration starting point, and widen the search radius to 0.35m to redetermine the integer ambiguity.

[0006] Preferably, in step S4, estimating the spatial skew vector of the receiving position caused by the refraction of the tree canopy includes the following sub-steps: Step S41, extracting the three-dimensional projection vector of the Doppler transient residual of the historical epoch within a pre-set sliding time-domain window in the receiver coordinate system; Step S42, using a time-decreasing discrete exponential weighting factor to weight and sum each three-dimensional projection vector to obtain a weighted mean vector; Step S43, using a pre-set Doppler phase deviation mapping scaling factor to scale the weighted mean vector to generate the spatial skew vector at the current moment.

[0007] Preferably, step S2 extracts a set of high-quality satellites, including the following sub-steps: Step S21, extracts the high-order temporal fluctuation characteristics of the carrier-to-noise ratio of a single satellite within the sliding time window in real time, and calculates the standard deviation of the carrier-to-noise ratio within the sliding time window as an objective measure of the degree of multipath interference; Step S22, weights the inverse of the standard deviation to convert the fluctuation characteristics into variance weighting factors of the corresponding satellite carrier phase observation data, reducing the weight ratio of interfered satellites with drastic fluctuations in carrier-to-noise ratio in differential positioning solutions; Step S23, introduces variance weighting factors to weight the observations when constructing the geometric observation design matrix of the current epoch satellite constellation, so that high-quality satellite signals dominate the solution of the three-dimensional velocity vector.

[0008] Preferably, step S6, determining the integer ambiguity of the double-difference carrier phase within the deskew search cylinder, includes the following sub-steps: Step S61, using the spatial grid points covered by the deskew search cylinder as candidate integer ambiguity combinations and performing ambiguity downcorrelation transformation to compress the unknown space to be searched; Step S62, performing a spatial tree-like optimization search within the transformed unknown space to obtain the first ambiguity candidate vector with the smallest residual sum of squares and the second ambiguity candidate vector with the second smallest residual sum of squares; Step S63, calculating the ratio of the residual sum of squares of the second ambiguity candidate vector to the residual sum of squares of the first ambiguity candidate vector. When the ratio is greater than a preset threshold of 3.0, the first ambiguity candidate vector is determined to have passed the rationality check and the integer ambiguity is locked.

[0009] Preferably, in step S6, when the Doppler frequency shift of the high-quality satellite set is continuously interrupted for 3 epochs, the absolute spatial position is calculated using the double-difference pseudorange as the starting point for the integration of the relative motion trajectory, including the following sub-steps: Step S64, extract the pseudorange observation values ​​of the multi-frequency global navigation satellite system under the interrupted epoch and smooth them in the time domain to filter out high-frequency abrupt noise; Step S65, introduce the pseudorange observation data of the remote reference station, construct the double-difference pseudorange positioning equation between the positioning terminal and the reference station, and calculate the temporary absolute position coordinates with decimeter-level accuracy by least squares estimation; Step S66, use the temporary absolute position coordinates as the new spatial origin of the Doppler integration, and update the state covariance matrix at the same time.

[0010] Preferably, the generation of the relative motion trajectory in step S2 further includes the following sub-steps: Step S24, introducing a priori constraint rule for the movement velocity of the geological hazard body to determine the maximum allowable step size upper limit of the displacement state vector between adjacent epochs; Step S25, comparing the three-dimensional displacement increment calculated in the current epoch with the step size upper limit during the state prediction stage, and determining that there is a carrier phase abrupt cycle slip caused by external interference in the corresponding epoch when the three-dimensional displacement increment exceeds the step size upper limit; Step S26, removing abnormal displacement increments that exceed the step size upper limit to reduce the influence of abnormal radio measurement jumps on the relative motion trajectory.

[0011] Preferably, in step S1, the raw radio observation data of the positioning terminal at the current epoch is obtained, wherein the satellites include a combination of multiple constellations of the BeiDou Navigation Satellite System, the Global Positioning System, and the GLONASS Navigation Satellite System, and the raw radio observation data covers multi-frequency radio signals.

[0012] Preferably, the sampling interval for acquiring the raw radio observation data of the positioning terminal at the current epoch is preset to 0.1s to 1.0s, and the length of the preset sliding time domain window in step S4 is preset to 10 epochs to 20 epochs.

[0013] Preferably, after outputting the deformation trajectory in step S6, the following sub-steps are also included: Step S67, the deformation trajectory output in consecutive epochs is introduced into the time series analysis model to calculate the cumulative deformation, deformation rate, and deformation acceleration at the current moment; Step S68, the deformation rate and deformation acceleration are compared with a preset judgment threshold; Step S69, when the deformation rate exceeds the preset warning threshold for three consecutive epochs and the deformation acceleration is greater than zero, the radio signal transmission frequency of the positioning terminal is adjusted to change the reporting cycle of trajectory positioning.

[0014] Compared with existing technologies, the positioning calculation and trajectory correction method of the geological disaster monitoring terminal of this invention has the following advantages: 1. In the positioning calculation of the geological disaster monitoring terminal, the spatial skew vector of the receiving position is estimated based on the Doppler transient residual, and the deviation between the relative motion trajectory and the true carrier phase coordinate is corrected accordingly. This allows the search center of the integer ambiguity to be adjusted according to the phase lag caused by the refraction of the vegetation canopy. This reduces the impact of asymmetric signal interference on ambiguity search, avoids the true solution from falling outside the search range due to the offset of the search center, and improves the reliability of integer ambiguity fixation and the stability of the positioning trajectory in complex forest cover environments.

[0015] 2. Constructing a deskipation search cylinder using the corrected 3D search center can limit the search range of integer ambiguity to a region closer to the true position, reducing invalid candidate combinations, shortening search time, and lowering the possibility of incorrect fixation. When Doppler observations are interrupted due to external obstruction, a temporary absolute position can be obtained using the double-difference pseudorange and the integration starting point can be re-established. At the same time, appropriately widening the search range can enable faster recovery of the positioning solution, reduce position jumps caused by re-initialization, and maintain the continuity and smoothness of the receiver coordinate sequence.

[0016] 3. Combining the relative motion trajectory formed by Doppler frequency shift with dual-difference carrier phase positioning can provide continuous spatial constraints for integer ambiguity search, narrowing the search range. At the same time, the absolute position coordinates obtained after fixing the integer ambiguity can be used to correct the drift accumulated over time in the relative motion trajectory. The two observation results correct each other, which is beneficial to improving the positioning stability under leaf jitter, cycle slip and multipath interference conditions, enabling the monitoring terminal to continuously output a more accurate deformation trajectory. Attached Figure Description

[0017] Figure 1 This is an overall flowchart of the monitoring terminal positioning calculation and trajectory correction of the present invention; Figure 2 This is a logical association diagram of the monitoring terminal positioning calculation and trajectory correction of the present invention. Detailed Implementation

[0018] The technical solutions of the embodiments of this application will be clearly described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of this application are within the scope of protection of this application.

[0019] A method for positioning calculation and trajectory correction of a geological disaster monitoring terminal, comprising: Step S1: Obtain the raw radio observation data of the positioning terminal at the current epoch. The raw radio observation data includes dual-difference pseudorange, carrier phase, Doppler frequency shift, and carrier-to-noise ratio. Step S2: Extract a set of high-quality satellites with a carrier-to-noise ratio greater than a preset threshold. Use the Doppler frequency shift of the high-quality satellite set between epochs to determine the three-dimensional velocity vector and accumulate the integral to generate the relative motion trajectory. Step S3: Subtract the actual Doppler frequency shift from the geometric prediction value calculated based on the relative motion trajectory to determine the Doppler transient residual of each high-quality satellite at the epoch; Step S4: Project each Doppler transient residual onto the receiver coordinate system and smooth and filter it within a pre-set sliding time-domain window to estimate the spatial skew vector of the receiver position caused by the refraction of the tree canopy. Step S5: Calculate the vector sum of the relative motion trajectory and the spatial skew vector to determine the three-dimensional search center. Using this as the central axis and the standard deviation of the error as the search radius, construct a deskew search cylinder. Step S6: Determine the integer ambiguity of the double-difference carrier phase within the deskew search cylinder to solve the absolute position coordinates and output the deformation trajectory. At the same time, feed back the absolute position coordinates to clear the cumulative drift of the relative motion trajectory. If the Doppler frequency shift of the high-quality satellite set is continuously interrupted for 3 epochs, use the double-difference pseudorange to solve the absolute spatial position as the integration starting point, and widen the search radius to 0.35m to redetermine the integer ambiguity.

[0020] Preferably, in step S4, estimating the spatial skew vector of the receiving position caused by the refraction of the tree canopy includes the following sub-steps: Step S41, extracting the three-dimensional projection vector of the Doppler transient residual of the historical epoch within a pre-set sliding time-domain window in the receiver coordinate system; Step S42, using a time-decreasing discrete exponential weighting factor to weight and sum each three-dimensional projection vector to obtain a weighted mean vector; Step S43, using a pre-set Doppler phase deviation mapping scaling factor to scale the weighted mean vector to generate the spatial skew vector at the current moment.

[0021] Preferably, step S2 extracts a set of high-quality satellites, including the following sub-steps: Step S21, extracts the high-order temporal fluctuation characteristics of the carrier-to-noise ratio of a single satellite within the sliding time window in real time, and calculates the standard deviation of the carrier-to-noise ratio within the sliding time window as an objective measure of the degree of multipath interference; Step S22, weights the inverse of the standard deviation to convert the fluctuation characteristics into variance weighting factors of the corresponding satellite carrier phase observation data, reducing the weight ratio of interfered satellites with drastic fluctuations in carrier-to-noise ratio in differential positioning solutions; Step S23, introduces variance weighting factors to weight the observations when constructing the geometric observation design matrix of the current epoch satellite constellation, so that high-quality satellite signals dominate the solution of the three-dimensional velocity vector.

[0022] Preferably, step S6, determining the integer ambiguity of the double-difference carrier phase within the deskew search cylinder, includes the following sub-steps: Step S61, using the spatial grid points covered by the deskew search cylinder as candidate integer ambiguity combinations and performing ambiguity downcorrelation transformation to compress the unknown space to be searched; Step S62, performing a spatial tree-like optimization search within the transformed unknown space to obtain the first ambiguity candidate vector with the smallest residual sum of squares and the second ambiguity candidate vector with the second smallest residual sum of squares; Step S63, calculating the ratio of the residual sum of squares of the second ambiguity candidate vector to the residual sum of squares of the first ambiguity candidate vector. When the ratio is greater than a preset threshold of 3.0, the first ambiguity candidate vector is determined to have passed the rationality check and the integer ambiguity is locked.

[0023] Preferably, in step S6, when the Doppler frequency shift of the high-quality satellite set is continuously interrupted for 3 epochs, the absolute spatial position is calculated using the double-difference pseudorange as the starting point for the integration of the relative motion trajectory, including the following sub-steps: Step S64, extract the pseudorange observation values ​​of the multi-frequency global navigation satellite system under the interrupted epoch and smooth them in the time domain to filter out high-frequency abrupt noise; Step S65, introduce the pseudorange observation data of the remote reference station, construct the double-difference pseudorange positioning equation between the positioning terminal and the reference station, and calculate the temporary absolute position coordinates with decimeter-level accuracy by least squares estimation; Step S66, use the temporary absolute position coordinates as the new spatial origin of the Doppler integration, and update the state covariance matrix at the same time.

[0024] Preferably, the generation of the relative motion trajectory in step S2 further includes the following sub-steps: Step S24, introducing a priori constraint rule for the movement velocity of the geological hazard body to determine the maximum allowable step size upper limit of the displacement state vector between adjacent epochs; Step S25, comparing the three-dimensional displacement increment calculated in the current epoch with the step size upper limit during the state prediction stage, and determining that there is a carrier phase abrupt cycle slip caused by external interference in the corresponding epoch when the three-dimensional displacement increment exceeds the step size upper limit; Step S26, removing abnormal displacement increments that exceed the step size upper limit to reduce the influence of abnormal radio measurement jumps on the relative motion trajectory.

[0025] Preferably, in step S1, the raw radio observation data of the positioning terminal at the current epoch is obtained, wherein the satellites include a combination of multiple constellations of the BeiDou Navigation Satellite System, the Global Positioning System, and the GLONASS Navigation Satellite System, and the raw radio observation data covers multi-frequency radio signals.

[0026] Preferably, the sampling interval for acquiring the raw radio observation data of the positioning terminal at the current epoch is preset to 0.1s to 1.0s, and the length of the preset sliding time domain window in step S4 is preset to 10 epochs to 20 epochs.

[0027] Preferably, after outputting the deformation trajectory in step S6, the following sub-steps are also included: Step S67, the deformation trajectory output in consecutive epochs is introduced into the time series analysis model to calculate the cumulative deformation, deformation rate, and deformation acceleration at the current moment; Step S68, the deformation rate and deformation acceleration are compared with a preset judgment threshold; Step S69, when the deformation rate exceeds the preset warning threshold for three consecutive epochs and the deformation acceleration is greater than zero, the radio signal transmission frequency of the positioning terminal is adjusted to change the reporting cycle of trajectory positioning.

[0028] Example 1: Currently, long-term deformation monitoring points for geological hazards such as landslides, collapses, and debris flows are mostly located in high mountains and valleys. When the positioning terminal antenna is blocked by dense forests, the wind-induced shaking of branches and leaves and surface rainwater will aggravate the multipath interference of radio signals. The carrier signal passing through the dynamic tree leaves will generate asymmetric phase lag deviation, which can easily cause tracking loop lock-up and high-frequency cycle slip. Conventional double-difference carrier phase positioning needs to go through a refinement and convergence process when re-determining the integer ambiguity. Simply increasing the observation noise covariance of the Kalman filter will reduce the positioning accuracy. In severe cases, it will cause trajectory divergence, resulting in measurement pseudo-abrupt changes, false alarms, and breakpoints in the output deformation trajectory.

[0029] The positioning terminal's navigation processor collects dual-difference pseudorange, carrier phase, Doppler frequency shift, and carrier-to-noise ratio at sampling intervals to form raw radio observation data. For each visible satellite within the current epoch, the navigation processor compares the real-time carrier-to-noise ratio with a pre-set threshold, extracts satellites with a carrier-to-noise ratio greater than the threshold, and forms a set of high-quality satellites. This threshold is pre-set based on the minimum carrier-to-noise ratio that the receiver tracking loop can maintain a stable lock and the noise level during static observation at the monitoring point.

[0030] For each satellite in the high-quality satellite set, the navigation processor extracts the temporal high-order fluctuation characteristics of its carrier-to-noise ratio (CNR) within the CNR sliding time window. Specifically, the deviation of each CNR sample value relative to the window mean is taken as the fluctuation sample, and the standard deviation of the CNR sequence is used to characterize the temporal high-order fluctuation characteristics. The inverse of the standard deviation of each satellite is normalized and used as the variance weighting factor. When constructing the geometric observation design matrix of the current epoch satellite constellation, the variance weighting factor is used to weight the observations of the corresponding satellites. Satellites with smaller CNR fluctuations receive larger weights, while satellites with larger CNR fluctuations receive smaller weights, thus enabling satellites with higher signal quality to dominate the solution of the three-dimensional velocity vector.

[0031] The navigation processor utilizes a high-quality satellite array in the current epoch. Compared with the previous era The Doppler frequency shift between them determines the three-dimensional velocity vector, and is based on the sampling interval. By accumulating and integrating, a relative motion trajectory is generated. The calculation formula is shown in equation (1): ,in, For the current epoch The relative motion trajectory, For the previous epoch The relative motion trajectory, For the current epoch The corresponding three-dimensional velocity vector, The sampling interval is denoted as .

[0032] To reduce abnormal radio measurement jumps caused by physical collisions, equipment vibrations, or momentary obstruction from tree branches and leaves, the navigation processor determines the maximum permissible step size upper limit of the displacement state vector between adjacent epochs based on the prior constraint rule of the movement velocity of the geological hazard body. This upper limit is based on the product of the pre-determined maximum permissible movement velocity of the monitored object and the sampling interval, and takes into account the uncertainty margin of the positioning solution. The navigation processor compares the magnitude of the three-dimensional displacement increment of the current epoch with the maximum permissible step size upper limit. When the three-dimensional displacement increment exceeds this upper limit, it is determined that there is a carrier phase abrupt cycle slip caused by external interference in the corresponding epoch. During the determination process, the carrier phase epoch difference and the continuous change of the double-difference pseudorange of the corresponding epoch are checked at the same time. When the carrier phase abruptly changes, and the change of the double-difference pseudorange cannot support the same amplitude of the real displacement, the corresponding three-dimensional displacement increment is recorded as an abnormal displacement increment and discarded, and is not used to update the relative motion trajectory.

[0033] When the double-difference carrier phase integer ambiguity of the current epoch is not yet fixed, the navigation processor uses the absolute position coordinates obtained from the previous epoch as the reference space origin, and superimposes the relative motion trajectory displacement of the current epoch onto this reference space origin to obtain the predicted absolute position coordinates of the receiver phase center at the current epoch. Subsequently, based on the satellite orbit prediction information, the three-dimensional position coordinates and velocity vectors of each satellite at the current epoch are calculated, and the predicted geometric distance vector from the receiver to each satellite is determined by the predicted absolute position coordinates and the three-dimensional position coordinates of the satellites. The relative velocity between the satellite and the receiver is projected onto the direction of the predicted geometric distance vector to obtain the predicted distance change rate after eliminating the influence of receiver clock drift. Then, the predicted distance change rate is multiplied by the carrier center frequency and divided by the speed of light to obtain the geometric Doppler prediction value of each high-quality satellite at the current epoch.

[0034] The relative motion trajectory formed by the Doppler frequency shift integral due to the refraction effect of the tree canopy A dynamic skew will occur between the actual Doppler frequency shift and the true carrier phase coordinate center. The navigation processor subtracts the geometric Doppler prediction value from the actual Doppler frequency shift of each high-quality satellite to obtain the Doppler transient residual of the corresponding satellite at the current epoch. For each Doppler transient residual, according to the line-of-sight unit vector of the corresponding satellite in the receiver coordinate system, the residual is projected onto the three coordinate axes of the receiver coordinate system to form a three-dimensional projection vector.

[0035] Within a Doppler transient residual sliding time-domain window consisting of 10 to 20 epochs, the navigation processor reads the three-dimensional projection vectors corresponding to each historical epoch. It then assigns time-decreasing discrete exponential weighting factors based on the distance between the epoch and the current time, normalizing the sum of all discrete exponential weighting factors to 1. A weighted mean vector is obtained by weighted summation of the three-dimensional projection vectors. This weighted mean vector is then scaled using a pre-calibrated Doppler-to-phase-bias mapping coefficient to generate the spatial skew vector of the receiving position at the current time. .

[0036] relative motion trajectory Spatial skew vector of the receiving location Adding them together, we obtain the corrected 3D search center. The calculation formula is shown in the figure: ,in, The corrected 3D search center, The relative trajectory at the current moment, This is the spatial skew vector of the receiving position at the current moment.

[0037] Navigation processor with 3D search center Using the central axis as the search radius, the standard deviation of the error for the current cumulative time period is calculated based on the standard deviation of the Doppler integral single-epoch solution and the cumulative number of time epochs since the last fixation of the double-difference carrier phase integer ambiguity. Construct a deflection search cylinder, and calculate the search radius using the formula: ,in, To determine the search radius of the cylinder for deskewing, Calculate the standard deviation for the single-epoch solution of the Doppler integral. The cumulative number of time epochs since the last successful fixation of the double-difference carrier phase integer ambiguity is given. The constant 3 is the standard deviation expansion factor. By translating the three-dimensional search center, the central axis of the de-skew search cylinder is synchronously adjusted with the spatial skew vector of the receiving position so that it covers the carrier phase coordinate area affected by the refraction of the tree canopy.

[0038] The navigation processor traverses the spatial grid points covered by the deskew search cylinder. For each spatial grid point, it calculates the corresponding double-difference geometric distance based on the geometric relationship between the positioning terminal, the remote reference station, the reference satellite, and other high-quality satellites. It then substitutes the double-difference geometric distance into the double-difference carrier phase observation equation to inversely calculate the integer number of each satellite relative to the reference satellite. Based on integer constraints, it generates candidate integer ambiguity combinations corresponding to the spatial grid point, and retains only one set of duplicate combinations.

[0039] A fuzziness downcorrelation transformation is performed on the candidate integer fuzziness combinations to compress the unknown space to be searched. Then, a spatial tree-like optimization search is performed in the transformed unknown space to obtain the first fuzziness candidate vector with the smallest residual sum of squares and the second fuzziness candidate vector with the second smallest residual sum of squares. The ratio of the residual sum of squares of the second fuzziness candidate vector to the residual sum of squares of the first fuzziness candidate vector is calculated. When the ratio is greater than 3.0, the first fuzziness candidate vector is determined to have passed the rationality test, and the double-difference carrier phase integer fuzziness is locked accordingly.

[0040] After the ambiguity is locked for the entire cycle, the navigation processor calculates the absolute position coordinates of the current epoch and outputs the deformation trajectory. At the same time, it feeds back the absolute position coordinates to the integral iteration process of the relative motion trajectory. The absolute position coordinates are used to correct the current position state of the relative motion trajectory, clear the accumulated drift formed by the Doppler integral, and enable the Doppler frequency shift relative integral and the double-difference carrier phase absolute positioning to form a two-way correction.

[0041] When a sudden change in the external environment or a sudden strong interference causes a continuous interruption in the Doppler frequency shift of the high-quality satellite array for 3 epochs, the navigation processor stops the original accumulation and integration, performs time-domain smoothing on the pseudorange observations of the multi-frequency global navigation satellite system up to the interruption epoch, filters out high-frequency sudden noise, introduces pseudorange observation data from the remote reference station, constructs a double-difference pseudorange positioning equation between the positioning terminal and the reference station, and solves the temporary absolute position coordinates with decimeter-level accuracy through least squares estimation.

[0042] The navigation processor uses the temporary absolute position coordinates as the new spatial origin for the Doppler integral. Simultaneously, it updates the state covariance matrix describing the uncertainties of the position and three-dimensional velocity states. During the update, the position covariance estimated by double-difference pseudorange least squares is written into the submatrix corresponding to the position state in the state covariance matrix, and the cross-covariance between the position state and other states before the reset is cleared, ensuring the state covariance matrix is ​​consistent with the new integration starting point. Subsequently, the search radius of the deskew search cylinder is... The width was temporarily increased to 0.35m.

[0043] After the positioning terminal reacquires the satellite and restores Doppler frequency shift tracking and carrier phase lock, the navigation processor performs a new grid search within the widened deskipation search cylinder, relocks integer ambiguities within 2 to 3 epochs, and after regaining the absolute position coordinates, restarts the calculation. The search radius is gradually reduced according to the standard deviation of the restored Doppler integral single epoch solution, so that it returns to the millimeter range, thereby reducing the position step caused by the integral reset.

[0044] Through the above processing, the relative motion trajectory generated by the Doppler frequency shift integral continuously provides spatial constraints for integer ambiguity search. The spatial skew vector of the received position is used to correct the search center, and the absolute position coordinates obtained by the double-difference carrier phase solution are used to eliminate the cumulative drift of the relative motion trajectory. Even when tree obstruction causes multipath interference and frequent cycle slips, the positioning terminal can still continuously output smooth millimeter-level deformation trajectories. After the Doppler frequency shift is continuously interrupted, the double-difference pseudorange can also be used to re-establish the integration starting point and restore the integer ambiguity solution.

[0045] Example 2: This example uses a high-precision radio signal simulation source and a Global Navigation Satellite System (GNSS) positioning receiver test bench to build a physical experimental platform. The high-precision radio signal simulation source has an RF signal output range of 1.2 GHz to 1.6 GHz, a signal power control accuracy of 0.1 dB, and a sampling frequency of 20 MHz. It is used to simulate the geometric distribution of visible satellite constellations and RF signal transmission in mountainous and canyon areas. The GNSS positioning receiver test bench has sub-millimeter level three-dimensional absolute coordinate reference positioning accuracy. It is used to carry the positioning terminal to be tested and measure its real-time relative motion displacement.

[0046] When the processor acquires raw radio observation data, the sampling interval is... Based on the upper limit of the landslide creep velocity and the time-domain changes of the tracking loop lockout transient, the creep condition in the early stage of a landslide requires a high trajectory sampling frequency. To reduce the distortion of the relative motion trajectory caused by undersampling, this embodiment adjusts the sampling interval. Set to 0.1s, the navigation processor determines the three-dimensional velocity vector based on the Doppler frequency shift of the high-quality satellite set, multiplies the three-dimensional velocity vector by the sampling interval to obtain the displacement increment of the current epoch, and then adds it to the relative motion trajectory of the previous epoch to generate the relative motion trajectory of the current epoch.

[0047] The length of the sliding time-domain window used to estimate the Doppler transient residual of the received location spatial skew vector. The window length is determined based on the wind-induced shaking period of the tree leaf layer and the transient fluctuation frequency of the carrier-to-noise ratio. If the window is too short, the number of residual samples in the historical epochs is insufficient, and the smoothing result is easily affected by transient clutter. If the window is too long, the smoothing result will be delayed in response to the dynamic changes of the tree leaf layer. To balance residual smoothing and transient response, this embodiment sets the window length to... It is set between 10 and 20 epochs.

[0048] During the experiment, a high-precision radio signal simulation source injected multipath clutter interference with an envelope following a Rayleigh distribution into the RF output signal. The dynamic range of the time delay was set to 10ns to 50ns, and a random load-to-noise ratio drop of 3dBHz to 8dBHz was superimposed as an operating condition disturbance. The impact of the noise level was determined by whether the spatial skew vector of the receiver position was estimated and the window length. Whether within the specified range, sample group 1, sample group 2, sample group 3, control group 1, control group 2, control group 3, out-of-range control group 1, out-of-range control group 2, and out-of-range control group 3 of this invention are set up respectively. Each group operates under the same disturbance conditions. The navigation processor calculates the received position spatial skew vector from the Doppler transient residual. The relative motion trajectory and the spatial skew vector of the receiving position are vectored together to determine the three-dimensional search center, and the search radius is determined based on the standard deviation of the Doppler integral single epoch solution, the cumulative number of time epochs, and the 3-fold expansion coefficient. Subsequently, the double-difference carrier phase integer ambiguity was determined within the deskew search cylinder.

[0049] When the carrier-to-noise ratio perturbation amplitude is 3.1 dBHz, the window length of sample group 1 of this invention is... The spatial skew vector of the receiver position is calculated using a time interval of 10 epochs. The modulus is 0.06m, and the search radius is... With a radius of 0.008m, after the spatial skew vector of the receiving position is shifted, the success rate of fixing the integer ambiguity of the double-difference carrier phase is 99.1%, and the three-dimensional absolute positioning error is 1.3mm. Control group 1 uses the same 10 epoch windows and a search radius of 0.008m, but does not calculate the spatial skew vector of the receiving position, and its magnitude is kept at 0.00m. The three-dimensional search center is not shifted accordingly, and the success rate of fixing the integer ambiguity of the double-difference carrier phase is 72.3%, and the three-dimensional absolute positioning error is 34.2mm.

[0050] When the carrier-to-noise ratio perturbation amplitude increases to 5.2 dBHz, sample group 2 of the present invention still maintains the window length. The spatial skew vector of the receiver position is calculated using a time interval of 10 epochs. The modulus is 0.12m, and the search radius is 0.012m. After three-dimensional search center translation and integer ambiguity search, the success rate of fixing the integer ambiguity of the double-difference carrier phase was 98.2%, and the three-dimensional absolute positioning error was 2.4 mm. The window length of control group 2 was also 10 epochs, and the search radius was also 0.012 m, but the spatial skew vector of the receiving position was different. With a modulus of 0.00m, the success rate of fixing the integer ambiguity of the double-difference carrier phase drops to 41.5%, and the three-dimensional absolute positioning error increases to 112.5mm.

[0051] Under the same 5.2dBHz carrier-to-noise ratio perturbation amplitude and 0.012m search radius, the out-of-range control group 1 will have a window length of... Shortened to 5 epochs, the calculated receiver location spatial skew vector is obtained due to the reduction in the number of historical Doppler transient residual samples involved in smoothing. With a modulus of 0.08m, the success rate of fixing double-difference carrier phase integer ambiguity was 61.4%, and the three-dimensional absolute positioning error was 54.2mm. The window length exceeded the range of control group 2. Extended to 25 epochs, receiving the position space skew vector The modulus is 0.19m; the vector is affected by the window response delay and cannot reflect the transient skew changes caused by the tree leaf layer in a timely manner. The success rate of fixing the double-difference carrier phase integer ambiguity is 55.8%, and the three-dimensional absolute positioning error is 72.6mm.

[0052] When the carrier-to-noise ratio perturbation amplitude further increases to 7.8 dBHz, the window length of sample group 3 of this invention will be... The spatial skew vector of the receiver position is calculated using a time interval of 15 epochs. The modulus is 0.21m, and the search radius is... With a radius of 0.018m, after the navigation processor translates the 3D search center based on this vector, the success rate of fixing the integer ambiguity of the double-difference carrier phase is 95.4%, and the 3D absolute positioning error is 4.6mm. Control group 3 uses the same 15 epoch windows and a search radius of 0.018m, but the received position spatial skew vector... With a modulus of 0.00m, the success rate of fixing the double-difference carrier phase integer ambiguity is 12.1%, and the three-dimensional absolute positioning error is 358.1mm.

[0053] The out-of-range control group 3, under a carrier-to-noise ratio perturbation amplitude of 7.8 dBHz, adjusted the window length. Shortened to 4 epochs and the search radius The spatial skew vector of the receiving position is calculated using a short window with a setting of 0.018m. With a modulus of 0.09m, it is greatly affected by instantaneous multipath clutter. The success rate of fixing the phase integer ambiguity of the double-difference carrier is 48.7%, and the three-dimensional absolute positioning error is 93.4mm.

[0054] In the above tests, the sample group of the present invention exponentially weights the three-dimensional projection vector of the Doppler transient residual according to the epochal order, and uses the weighting result to calculate the spatial skew vector of the receiving position. The system then translates the 3D search center accordingly. The control group, operating at the same search radius, did not perform this translation. The out-of-range control group used a window length outside the defined range. Test results show that the received position spatial skew vector can correct the search center offset caused by tree refraction. A window length within the defined range can simultaneously maintain residual smoothing and response speed to transient skew changes. As the carrier-to-noise ratio disturbance amplitude increases, the navigation processor updates the received position spatial skew vector. And adjust the search radius This ensures a high success rate in fixing double-difference carrier phase integer ambiguity and controls the three-dimensional absolute positioning error to the millimeter level.

[0055] In the above experiment, the positioning terminal used the original hardware structure. It completed trajectory correction through Doppler frequency shift differential, residual smoothing and double-difference carrier phase calculation. The relative motion trajectory provides continuous spatial constraints for integer ambiguity search. The received position spatial skew vector is used to correct the central axis of the deskew search cylinder. The absolute position coordinates obtained after fixing the integer ambiguity are fed back to the integration process of the relative motion trajectory to clear the accumulated drift. Under the conditions of high frequency cycle slip and dynamic multipath interference, the positioning terminal can continuously output smooth millimeter-level deformation trajectory.

[0056] Example 3: In this example, the geological disaster monitoring terminal is installed on steep slopes and in areas with dense shrub cover. When the leaves of the trees sway in the wind, the radio signal received by the positioning terminal will generate non-uniform multipath diffraction. As a result, the numerically controlled oscillator in the carrier phase tracking loop will experience frequency jitter and phase tracking lag. When the cycle slip detector fails to identify the sudden change of high dynamic signal in time, the carrier phase observation data stored in the receiver's dynamic random access memory will contain nonlinear phase deviation. This deviation is reflected in both the carrier phase and Doppler frequency shift observation results.

[0057] Navigation processor specifications are as follows OK, The column's floating-point temporary buffer, where the sliding time-domain window length Set to 10, the total number of visible satellites in the current high-quality satellite set. The floating-point buffer is set to 8. It stores the difference between the actual Doppler frequency shift and the geometric Doppler prediction value of each high-quality satellite in epoch order to form a Doppler transient residual sequence. This sequence comes from the baseband digital signal output by the antenna low-noise amplifier and after down-conversion processing. The embedded processor implementing this process operates at a frequency of no less than 200MHz and can complete the residual projection, weighted summation and subsequent matrix dimensionality reduction operations within a single epoch period of 0.1s.

[0058] After the floating-point buffer completes the data writing for the current window, the navigation processor reads the Doppler transient residuals of each historical epoch. For each high-quality satellite, the corresponding Doppler transient residual is first projected onto the three orthogonal axes of the receiver coordinate system based on the three-dimensional line-of-sight unit vector of the satellite relative to the receiver, thus obtaining the three-dimensional projection vector at that epoch. Then, in order from the most recent to the oldest time, discrete exponential weighting factors that decay with time are assigned to each three-dimensional projection vector within the sliding time window. The sum of each weighting factor is 1. After weighted summation of the three-dimensional projection vectors, the weighted mean vector of the current epoch is obtained.

[0059] The navigation processor utilizes a pre-calibrated Doppler phase deviation mapping scaling factor. Scale the weighted mean vector to generate a spatial skew vector of the receiving location. Doppler phase deviation mapping scaling factor The dimension of the value is m / Hz. In this embodiment, a precision servo slide rail is used to control the positioning terminal to move at a constant speed of 1m / s in a windless environment. The known Doppler residual and carrier phase lag are collected synchronously, and the coefficient is determined based on the statistical fitting results of the two. The value of the coefficient is 0.19.

[0060] Obtain the spatial skew vector of the receiving position Then, the navigation processor adds it to the current relative motion trajectory vector to obtain the three-dimensional search center after skew correction. This three-dimensional search center serves as the central axis of the deskew search cylinder. The search radius is determined based on the standard deviation of the Doppler integral single epoch solution and the cumulative number of time epochs. Thus, the non-zero mean deviation caused by the tree canopy is converted into the translation of the search center. The deskew search cylinder moves synchronously with the spatial skew vector of the receiving position, covering the carrier phase coordinate area affected by the refraction of the tree canopy.

[0061] Within the deskew search cylinder, the navigation processor converts spatial grid points into candidate integer ambiguity combinations and sequentially performs ambiguity downcorrelation transformation and spatial tree-based optimization search to obtain the first ambiguity candidate vector with the smallest residual sum of squares and the second ambiguity candidate vector with the second smallest residual sum of squares. When the ratio of the residual sum of squares of the second ambiguity candidate vector to that of the first ambiguity candidate vector is greater than 3.0, the double-difference carrier phase integer ambiguity corresponding to the first ambiguity candidate vector is locked, the absolute position coordinates of the current epoch are calculated, and the deformation trajectory is output.

[0062] After the absolute position coordinates are fed back to the integration process of the relative motion trajectory, they are used to clear the cumulative drift formed by the Doppler frequency shift accumulation. Under dynamic multipath interference, the navigation processor can still adjust the received position spatial skew vector according to the continuously updated Doppler transient residual, so that the central axis of the deskew search cylinder moves with the actual signal skew direction, reducing the situation where the true solution of integer ambiguity falls outside the search range. The positioning terminal thus continuously outputs the three-dimensional deformation trajectory, reducing the possibility that the phase lag caused by tree refraction is misjudged as a sudden displacement of geological disaster body.

[0063] Example 4: This example combines Figures 1 to 2 The method for positioning calculation and trajectory correction of geological disaster monitoring terminals is explained, such as... Figure 1As shown, the process involves acquiring the raw radio observation data of the positioning terminal at the current epoch. This raw radio observation data includes double-difference pseudorange, carrier phase, Doppler frequency shift, and carrier-to-noise ratio (CNR). Step S1 extracts a set of high-quality satellites with a CNR greater than a pre-set threshold. The three-dimensional velocity vector is determined using the Doppler frequency shift of the high-quality satellite set between the current epoch and the previous epoch, and then integrated to generate a relative motion trajectory. Step S2 subtracts the actual Doppler frequency shift from the geometric prediction value calculated based on the relative motion trajectory to determine the Doppler transient residual of each high-quality satellite at the current epoch. Step S3 projects each Doppler transient residual onto the receiver coordinate system and smooths it within a pre-set sliding time window to estimate the relative motion trajectory caused by the forest. The spatial skew vector of the receiving position caused by canopy refraction is calculated in step S4. The vector sum of the relative motion trajectory and the spatial skew vector is used to determine the three-dimensional search center. Using this as the central axis and the standard deviation of the error as the search radius, a deskew search cylinder is constructed. In step S5, the double-difference carrier phase integer ambiguity is determined within the deskew search cylinder to solve the absolute position coordinates and output the deformation trajectory. At the same time, the absolute position coordinates are fed back to clear the cumulative drift of the relative motion trajectory. If the Doppler frequency shift of the high-quality satellite set is continuously interrupted for 3 epochs, the absolute spatial position is solved using the double-difference pseudorange as the integration starting point, and the search radius is widened to 0.35m to redetermine the integer ambiguity. This process corresponds to step S6.

[0064] like Figure 2 As shown, when generating a relative motion trajectory, if the three-dimensional displacement increment exceeds the upper limit of the step size, an operation to remove abnormal displacement increments is performed. After successfully generating the relative motion trajectory, the absolute position coordinates are calculated by locking the integer ambiguity of the dual-difference carrier phase. The absolute position coordinates obtained by calculating the absolute position coordinates are fed back through a feedback path to clear the accumulated drift, thus affecting the step of generating the relative motion trajectory. If the Doppler frequency shift is interrupted continuously for three epochs during system operation, the operation of calculating temporary absolute position coordinates is performed. After restoring Doppler frequency shift tracking and carrier phase locking to relock the integer ambiguity, the step of calculating absolute position coordinates is switched back. Finally, after calculating the absolute position coordinates and outputting the deformation trajectory, if the deformation rate exceeds the warning threshold for three consecutive epochs and the deformation acceleration is greater than zero, an adjustment operation to change the reporting cycle of trajectory positioning is triggered.

[0065] Example 5: Before collecting raw radio observation data, the navigation processor first performs zero-creep field initial adjustment. After the positioning terminal is fixed at the static reference point, it continuously latches the Doppler frequency shift observation values ​​of available satellites, with an observation length of 2000 epochs. Subsequently, it reads the double-difference pseudorange and carrier phase observation values ​​from the temporary buffer, calculates the initial baseline vector of the three-dimensional coordinates of the positioning terminal antenna phase center relative to the remote reference station, and uses this baseline vector to determine the absolute spatial starting point of the cumulative integration of the relative motion trajectory.

[0066] The navigation processor calculates the first-order time difference of the Doppler frequency shift and obtains the mean square error within the sliding time-domain window of the Doppler transient residual, thereby calibrating the standard deviation of the Doppler integral single-epoch solution. The initial value, in this embodiment, is the standard deviation of the Doppler integral single-epoch solution measured under static operating conditions. The value is 1.8 mm, which is used to determine the initial search range for the deflection search cylinder.

[0067] After the positioning terminal was installed on-site, the antenna receiving module continuously tracked satellite signals for 60 seconds in a stationary state. The navigation processor extracted the time-domain variation sequence of the carrier-to-noise ratio (CNR) for each visible satellite, calculated the standard deviation of each signal channel within the CNR sliding time window, and then calculated the average of the standard deviations. When the average value exceeded 2.5 dBHz, it indicated that the CNR fluctuation caused by the shading of the tree canopy had reached a level requiring a wider search range. The navigation processor then accumulated the time epochs. Zero out and set the search radius. The initial value was set to 0.035m to cover centimeter-level positional deviations caused by multipath interference.

[0068] After the double-difference carrier phase integer ambiguity is locked within a single epoch, the navigation processor feeds back the calculated absolute position coordinates to the integration process of the relative motion trajectory, clearing the accumulated drift caused by the Doppler frequency shift integration. After the feedback is completed, the absolute position coordinate feedback loop is closed, and the search radius is gradually reduced through the parameter down-adjustment channel. This allows the dynamic convergence to within 3mm, enabling the positioning terminal to continuously output a smooth deformation trajectory in high mountain and canyon forest environments.

[0069] Search radius under different levels of signal degradation The adjustment is tiered based on whether the Doppler frequency shift is continuously tracked. When the Doppler frequency shift of the high-quality satellite set is continuously interrupted for 3 epochs, the cumulative integration of the relative motion trajectory loses continuity. The navigation processor uses dual-difference pseudorange to calculate temporary absolute position coordinates with decimeter-level accuracy and uses them as the new integration starting point. At this time, in order to cover the decimeter-level position uncertainty caused by the dual-difference pseudorange reset, the search radius is adjusted. The cylinder was temporarily widened to 0.35m, and the double-difference carrier phase integer ambiguity was redefined within the widened deskew search cylinder.

[0070] When the Doppler frequency shift does not experience continuous interruption, but the average standard deviation of the carrier-to-noise ratio fluctuation exceeds 2.5 dBHz, the relative motion trajectory remains continuous, affected only by centimeter-level deviations caused by multipath interference from trees. Under this condition, the search radius is... Setting the search radius to 0.035m covers the deviation range of the relative motion trajectory. The two search radii of 0.35m and 0.035m correspond to the decimeter-level uncertainty when resetting the integration starting point with double-difference pseudorange, and the centimeter-level uncertainty when continuously tracking Doppler frequency shift, respectively. The navigation processor selects the appropriate search range accordingly, controlling the number of spatial grid points while covering the true solution of integer ambiguity.

[0071] Example 6: Before deploying the current geological disaster monitoring terminal in the field, the on-site benchmark calibration and radio frequency channel calibration under the condition of multipath interference of forest trees are completed to determine the quantitative mapping relationship between Doppler transient residual and carrier phase lag deviation. During calibration, the positioning terminal is fixed on a static benchmark pier with known geographical coordinates and stable structure, so that the geometric distribution of the visible satellite constellation in the antenna field of view meets the positioning solution requirements. The navigation processor observes at the monitoring point for 300s, continuously collects the Doppler frequency shift observation value and carrier phase observation value of each satellite, calculates the Doppler frequency shift differential mean of continuous epochs, and forms a carrier phase double difference residual sequence.

[0072] The navigation processor uses the known three-dimensional coordinates of the remote reference station and the positioning terminal antenna, combined with the satellite orbit information of the current epoch, to calculate the transient geometric distances from the reference station to each visible satellite and from the positioning terminal to each visible satellite. It then performs inter-station and inter-satellite differences on these transient geometric distances to obtain the double-difference geometric distance for the current epoch. Finally, it subtracts the double-difference geometric distance of the previous epoch from the current epoch's double-difference geometric distance to obtain the double-difference geometric distance increment. The difference between the double-difference carrier phase observations of the current epoch and the previous epoch is used to obtain the double-difference carrier phase increment. The double-difference Doppler shift increment is obtained by subtracting the double-difference Doppler shift observations from those of the current epoch and the previous epoch. .

[0073] Doppler phase deviation mapping scaling factor The baseline value is calculated according to the following formula: ,in, The unit is m / Hz. The unit is m. The unit is m. The unit is Hz. During calibration, the positioning terminal remains stationary, and the change in geometric distance caused by receiver motion is zero or within the static observation error range. Based on this, the navigation processor separates the phase lag corresponding to the Doppler frequency shift from the carrier phase double-difference residual and determines the phase lag based on the calculation results of continuous epochs. The baseline value.

[0074] After calibrating the Doppler-to-phase deviation mapping scaling factor, the navigation processor reads the raw radio observation data stored in the buffer. It subtracts the double-difference pseudorange observation value from the previous epoch from the current epoch's double-difference pseudorange observation value to obtain the double-difference pseudorange change. This change is then subtracted from the corresponding satellite's double-difference carrier phase change to obtain the deviation residual. Finally, the deviation residual is used in conjunction with... The product of the two is used to perform frequency preset compensation on the receiver's digitally controlled oscillator, and the control word of the digitally controlled oscillator is adjusted to restore phase lock of the carrier phase tracking loop of the baseband electrical channel.

[0075] The navigation processor simultaneously reads the carrier-to-noise ratio (CNR) sampling sequence within the CNR sliding time window, calculates the time-domain fluctuation standard deviation of the CNR for each signal channel, and when the time-domain fluctuation standard deviation of a signal channel exceeds 2.5 dBHz, it adjusts the search radius in the status recording register. The control parameter was set to 0.035m to cover the positional deviation caused by multipath interference under dense vegetation conditions.

[0076] After calibration, the navigation processor will map the Doppler to the phase deviation scaling factor. and search radius The basic values ​​are written into non-volatile memory to form a reference mapping table between the terminal physical parameters and the algorithm control parameters. After the positioning terminal enters the working state, the navigation processor adjusts the bandwidth coefficient of the carrier phase tracking loop according to the collected values ​​of external power supply voltage and ambient temperature, so that the receiver can maintain continuous satellite observation and positioning calculation when the power supply voltage or ambient temperature changes.

[0077] The embodiments of this application have been described above with reference to the accompanying drawings. Unless otherwise specified, the embodiments and features in the embodiments of this application can be combined with each other. This application is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms under the guidance of this application without departing from the spirit of this application and the scope of protection of this invention, and all of these forms are within the protection scope of this application.

Claims

1. A method for positioning calculation and trajectory correction of a geological disaster monitoring terminal, characterized in that, include: Step S1: Obtain the raw radio observation data of the positioning terminal at the current epoch. The raw radio observation data includes dual-difference pseudorange, carrier phase, Doppler frequency shift, and carrier-to-noise ratio. Step S2: Extract a set of high-quality satellites with a carrier-to-noise ratio greater than a preset threshold. Use the Doppler frequency shift of the high-quality satellite set between epochs to determine the three-dimensional velocity vector and accumulate the integral to generate the relative motion trajectory. Step S3: Subtract the actual Doppler frequency shift from the geometric prediction value calculated based on the relative motion trajectory to determine the Doppler transient residual of each high-quality satellite at the epoch; Step S4: Project each Doppler transient residual onto the receiver coordinate system and smooth and filter it within a pre-set sliding time-domain window to estimate the spatial skew vector of the receiver position caused by the refraction of the tree canopy. Step S5: Calculate the vector sum of the relative motion trajectory and the spatial skew vector to determine the three-dimensional search center. Using this as the central axis and the standard deviation of the error as the search radius, construct a deskew search cylinder. Step S6: Determine the integer ambiguity of the double-difference carrier phase within the deskew search cylinder to solve the absolute position coordinates and output the deformation trajectory. At the same time, feed back the absolute position coordinates to clear the cumulative drift of the relative motion trajectory. If the Doppler frequency shift of the high-quality satellite set is continuously interrupted for 3 epochs, use the double-difference pseudorange to solve the absolute spatial position as the integration starting point, and widen the search radius to 0.35m to redetermine the integer ambiguity.

2. The method for positioning calculation and trajectory correction of a geological disaster monitoring terminal according to claim 1, characterized in that, Step S4 estimates the spatial skew vector of the receiving position caused by the refraction of the tree canopy, including the following sub-steps: Step S41, extract the three-dimensional projection vector of the Doppler transient residual of the historical epoch within the pre-set sliding time-domain window in the receiver coordinate system; Step S42, use a time-decreasing discrete exponential weighting factor to weight and sum the three-dimensional projection vectors to obtain a weighted mean vector; Step S43, use a pre-set Doppler phase deviation mapping scaling factor to scale the weighted mean vector to generate the spatial skew vector at the current moment.

3. The method for positioning calculation and trajectory correction of a geological disaster monitoring terminal according to claim 1, characterized in that, Step S2 extracts a set of high-quality satellites, including the following sub-steps: Step S21, extracts the high-order temporal fluctuation characteristics of the carrier-to-noise ratio of a single satellite within the sliding time window in real time, and calculates the standard deviation of the carrier-to-noise ratio within the sliding time window as an objective measure of the degree of multipath interference. Step S22: The inverse of the standard deviation is weighted to convert the fluctuation characteristics into the variance weighting factor of the corresponding satellite carrier phase observation data, thereby reducing the weight ratio of the interfered satellite with drastic fluctuations in carrier-to-noise ratio in the differential positioning solution. Step S23: When constructing the geometric observation design matrix of the current epoch satellite constellation, a variance weighting factor is introduced to weight the observations so that high-quality satellite signals dominate the solution of the three-dimensional velocity vector.

4. The method for positioning calculation and trajectory correction of a geological disaster monitoring terminal according to claim 1, characterized in that, Step S6 involves determining the integer ambiguity of the double-difference carrier phase within the deskew search cylinder, including the following sub-steps: Step S61: Using the spatial grid points covered by the deskew search cylinder as candidate integer ambiguity combinations and performing ambiguity downcorrelation transformation to compress the unknown space to be searched; Step S62: Performing a spatial tree-like optimization search within the transformed unknown space to obtain the first ambiguity candidate vector with the smallest residual sum of squares and the second ambiguity candidate vector with the second smallest residual sum of squares; Step S63: Calculating the ratio of the residual sum of squares of the second ambiguity candidate vector to the residual sum of squares of the first ambiguity candidate vector. When the ratio is greater than a preset threshold of 3.0, the first ambiguity candidate vector is determined to have passed the rationality check and the integer ambiguity is locked.

5. The method for positioning calculation and trajectory correction of a geological disaster monitoring terminal according to claim 1, characterized in that, In step S6, when the Doppler frequency shift of the high-quality satellite set is continuously interrupted for 3 epochs, the absolute spatial position is calculated using double-difference pseudorange as the starting point for the integration of the relative motion trajectory. This includes the following sub-steps: Step S64: Extract the pseudorange observation values ​​of the multi-frequency global navigation satellite system under the interrupted epoch and smooth them in the time domain to filter out high-frequency abrupt noise; Step S65: Introduce the pseudorange observation data of the remote reference station, construct the double-difference pseudorange positioning equation between the positioning terminal and the reference station, and calculate the temporary absolute position coordinates with decimeter-level accuracy through least squares estimation; Step S66: Use the temporary absolute position coordinates as the new spatial origin of the Doppler integration, and update the state covariance matrix at the same time.

6. The method for positioning calculation and trajectory correction of a geological disaster monitoring terminal according to claim 1, characterized in that, The generation of relative motion trajectory in step S2 also includes the following sub-steps: Step S24, introducing the prior constraint rule of the movement velocity of the geological hazard body, and determining the upper limit of the maximum allowable step size of the displacement state vector between adjacent epochs; Step S25: In the state prediction stage, compare the three-dimensional displacement increment calculated in the current epoch with the upper limit of the step size. When the three-dimensional displacement increment exceeds the upper limit of the step size, it is determined that there is a carrier phase abrupt cycle slip caused by external interference in the corresponding epoch. Step S26: Remove abnormal displacement increments that exceed the upper limit of the step size to reduce the impact of abnormal radio measurement jumps on the relative motion trajectory.

7. The method for positioning calculation and trajectory correction of a geological disaster monitoring terminal according to claim 1, characterized in that, In step S1, the raw radio observation data of the positioning terminal at the current epoch is obtained. The satellites include a combination of multiple constellations in the BeiDou Navigation Satellite System, the Global Positioning System, and the GLONASS Navigation Satellite System. The raw radio observation data covers multiple frequency radio signals.

8. The method for positioning calculation and trajectory correction of a geological disaster monitoring terminal according to claim 1, characterized in that, The sampling interval for acquiring the raw radio observation data of the positioning terminal at the current epoch is preset to 0.1s to 1.0s, and the length of the preset sliding time domain window in step S4 is preset to 10 epochs to 20 epochs.

9. The method for positioning calculation and trajectory correction of a geological disaster monitoring terminal according to claim 1, characterized in that, After outputting the deformation trajectory in step S6, the following sub-steps are also included: Step S67, the deformation trajectory output from the continuous epochs is introduced into the time series analysis model to calculate the cumulative deformation, deformation rate and deformation acceleration at the current moment. Step S68: Compare the deformation rate and deformation acceleration with the preset judgment threshold; Step S69: When the deformation rate exceeds the preset warning threshold for 3 consecutive epochs and the deformation acceleration is greater than zero, adjust the radio signal transmission frequency of the positioning terminal to change the reporting cycle of trajectory positioning.

Images