Real-time co-seismic displacement velocity solving method and system based on beidou PPP-B2b service

Abstract

The present application belongs to the technical field of information technology service, and particularly relates to a real-time co-seismic displacement velocity solving method and system based on Beidou PPP-B2b service, which uses the orbit, clock error and code bias correction number broadcast by the PPP-B2b service to solve the co-seismic displacement and velocity in real time, obtains the PPP-B2b text data broadcast by the Beidou GEO satellite, combines the broadcast ephemeris data to restore the precise clock error and orbit based on the PPP-B2b correction number, obtains the GNSS observation data, applies the PPP-B2b precise clock error and orbit to the ionosphere-free combined PPP observation equation, weights according to the satellite elevation angle, and constructs a random model. The present application provides a stable and reliable, low-cost real-time displacement, velocity and acceleration calculation method for the GNSS earthquake monitoring station in the strong earthquake near field without network communication based on the Beidou PPP-B2b service and the constant acceleration kinetic model, and provides observation data for the earthquake monitoring and early warning research, and can be widely applied in the satellite measurement field.

Description

Technical Field

[0001] This invention belongs to the field of information technology service technology, and in particular relates to a method and system for real-time coseismic displacement velocity calculation based on BeiDou PPP-B2b service. Background Technology

[0002] Currently, with the rapid development of GNSS (Global Navigation Satellite System), GNSS high-precision positioning technology has been widely applied in the field of earthquake monitoring. The high-precision GNSS surface coseismic displacement sequence acquired is a prerequisite for inverting the source rupture process and can provide technical support for earthquake emergency response.

[0003] There are two main methods for high-precision coseismic displacement extraction using GNSS: relative positioning and absolute positioning. Real-Time Kinematic (RTK) positioning requires one or more reference stations to remove common errors and recover double-difference ambiguities, achieving high-precision centimeter-level positioning. However, when the reference station is located in the seismic zone, the coseismic phenomena will interfere with the rover's displacement sequence. Furthermore, as the length of the GNSS baseline increases, the correlation between satellite ephemeris and atmospheric errors gradually weakens, affecting positioning accuracy.

[0004] Precise Point Positioning (PPP) can acquire high-precision absolute coseismic displacement sequences using only a single receiver. However, due to the fixed ambiguity, PPP and PPP-AR methods have long initialization times. If the earthquake happens to occur during the ambiguity convergence period, the accuracy of the acquired coseismic displacement sequence will be significantly affected. The Variomertic Approach Displacement Analysis Stand-alone Engine (VADASE) uses a single receiver to estimate velocity based on carrier phase observations differentiated between epochs. Its epoch-differential method effectively reduces the influence of ambiguity and atmospheric residual errors and requires no initialization time. However, the VADASE method requires integration to obtain displacement values, resulting in significant drift, which needs to be eliminated using additional linear or spatial filtering. Temporal Point Positioning (TPP) methods no longer limit the difference between epochs to adjacent epochs, avoiding the integration process and effectively mitigating the cumulative effect of residual errors in epoch-level differences. Furthermore, it requires high-precision receiver positioning, and when data is discontinuous, cycle slip detection and corrected values ​​need to be superimposed onto the GNSS observations obtained through temporal point difference. Compared to TPP, Precise Point Positioning (PPP) Velocity Estimation (PPP-VE) and Acceleration Estimation (PPP-AE) methods can directly estimate the station's velocity and acceleration vectors, and then obtain the coseismic displacement sequence through integration, effectively overcoming the shortcomings of the aforementioned GNSS coseismic displacement extraction methods. The International GNSS Service (IGS) has been providing Real-Time Service (RTS) via the internet since 2013. Using RTS products, users can implement real-time PPP applications. However, this service relies on a fixed internet connection. In actual earthquake monitoring, communication base stations may be damaged due to strong ground movement, making it difficult for IGS RTS precision ephemeris products to be transmitted over the network, thus making it difficult to obtain high-precision coseismic displacement sequences.

[0005] Besides services provided via internet communication, satellite broadcasting services have been advertised by some commercial companies at high costs. However, equipping large-scale GNSS seismic monitoring stations with satellite communication equipment is prohibitively expensive and difficult to implement in real time. In contrast, the PPP-B2b service of the BeiDou-3 system (BDS-3) is a free and open-source service, a significant innovation of the BeiDou system. This service can directly broadcast precise orbits and clock corrections via geostationary orbit (GEO) satellites, unaffected by network disruptions. It can also calculate high-precision displacement and velocity sequences in real time during strong earthquakes when the network is interrupted, making it applicable to earthquake monitoring and early warning systems.

[0006] Based on the above analysis, existing technologies for obtaining coseismic displacement velocity include commercial technologies such as RTK, RTS PPP, and satellite-based differential methods. The problems and defects are as follows: In the RTK method, the reference station is affected by the earthquake, which leads to an increase in the error of the station's solution and makes it impossible to obtain objective displacement velocity information; the RTS PPP method relies on the network, and there is a risk of network interruption during strong earthquakes, which will prevent coseismic displacement calculation; commercial technologies such as satellite-based differential methods are expensive and cannot be applied in real time in earthquake monitoring. Summary of the Invention

[0007] To address the problems existing in the prior art, this invention provides a method and system for real-time coseismic displacement velocity calculation based on BeiDou PPP-B2b service.

[0008] This invention is implemented as follows: a real-time coseismic displacement velocity calculation method based on BeiDou PPP-B2b service, the implementation steps of which include:

[0009] Step 1: Combine broadcast ephemeris to complete the precise ephemeris calculation method based on BeiDou-3 PPP-B2b service, mainly including the calculation of precise orbit and clock bias;

[0010] Step 2: Construct a real-time coseismic displacement calculation method, which mainly includes the preprocessing and error correction of GNSS observation data, the establishment of PPP observation and stochastic models, the construction of constant acceleration dynamic models, and the estimation methods of three-dimensional displacement, velocity and acceleration of monitoring stations.

[0011] Furthermore, in step one, in order to ensure the correlation between the information content broadcast by different information types, as well as the correlation between the correction data and the broadcast ephemeris, the information is identified by four version numbers: Spatial State Expression Data Period Number (IOD SSR), Pseudo-Random Noise Code Data Period Number (IODP), Navigation Data Period Number (IODN), and Correction Data Identification and Data Period Number (IOD Corr), to facilitate matching and use;

[0012] The orbit correction data in the PPP-B2b message has an epoch time interval of 48 seconds, and the broadcast orbit correction information includes the orbit correction vector. The satellite position correction vector can be calculated using the orbit correction values. After the IODN of the orbit correction information successfully matches the IODE in the wide-interpolation ephemeris navigation message, the precise orbit correction data in the satellite orbit coordinate system can be expressed as the precise orbit correction data in the geocentric-geocentric coordinate system.

[0013]

[0014] Among them, [δO r δO a δO c ] T These are the orbital corrections in the radial, tangential, and normal star-fixed (ACR) coordinate systems; [δO x δO y δO z ] T This refers to the orbital corrections for the satellite within the Earth-Centered Earth Fixed (ECEF) reference frame. r e a e c The unit vectors corresponding to the radial, tangential, and normal directions, respectively, can be expressed as:

[0015]

[0016] In the formula, r and To calculate satellite position and satellite velocity using broadcast ephemeris;

[0017] Calculations are performed on precision orbits based on PPP-B2b services:

[0018]

[0019] Among them, [XYZ] T These are the coordinates of the satellite within the Earth Center Earth Fixed (ECEF) reference frame.

[0020] Furthermore, in step one, the epoch time interval of the clock correction number in the PPP-B2b clock correction information is 6 seconds; after the clock correction information is successfully matched with the broadcast ephemeris, the clock correction parameter C0 contained in the message can be used to calculate the clock error parameter obtained from the broadcast ephemeris. Corrections are made to calculate the precision clock bias based on PPP-B2b services:

[0021]

[0022] in, The precise clock error recovered from the PPP-B2b signal, C0 is the clock error correction in the line-of-sight direction calculated using a third-order polynomial; when the track correction and clock error correction in the PPP-B2b message are the same as the IOD Corr in the clock error correction information, and their IOD can match the IODE in the navigation message, the track and clock error correction can be used to perform error correction.

[0023] Furthermore, in step two, the commonly used GNSS primitive pseudorange and phase observations in the PPP positioning model can be expressed as:

[0024]

[0025] In the formula, P i and L i The original code / phase observations; ρ represents the geometric distance from the receiver to the satellite; c represents the speed of light; δt r and δt s These represent the clock difference between the receiver and the satellite, respectively; T is the tropospheric delay; f i λ is the frequency value; I1 is the ionospheric delay of L1; i For wavelength, N i b represents integer ambiguity; r,i and These represent the receiver and satellite phase hardware offsets, respectively; B r,i and These represent the hardware deviations of the receiver and the satellite pseudorange, respectively. and These represent the unmodeled errors in pseudorange and carrier phase, including thermal noise and multipath error. Furthermore, errors such as unspecified relativistic effects, Sagnac effects, Shapiro delay, phase entanglement, and tidal effects have been corrected using existing models.

[0026] In the PPP-VE positioning model, a dual-frequency ionospheric de-de-spheric combination is used to eliminate the influence of the first-order ionospheric delay. The established observation equation for the dual-frequency ionospheric de-de-spheric combination can be expressed as:

[0027]

[0028] In the formula,

[0029] The weighting is performed using a satellite elevation angle stochastic model, expressed as follows:

[0030]

[0031] In the formula, w is the calculated observation weight matrix, p jFor the observation weighting factor of each satellite, el j Let be the satellite elevation angle of satellite j.

[0032] Furthermore, in step two, the dynamic model of PPP-VE during filtering calculations is generally based on Newton's second law, typically using dynamic position vectors, velocities, and accelerations to describe the vehicle's trajectory. In actual calculations, finite-order differentials are often used for approximation, and the constant acceleration (CA) dynamic model can be expressed as:

[0033]

[0034] In the formula, k and k-1 represent the current time and the previous time, respectively; τ is the sampling interval; w represents the state noise matrix; q a q is the power density of acceleration. a Set to 0.01 m·s -5 / 2 I is a three-dimensional identity matrix.

[0035] Because of the use of precise orbits and precise satellite clock biases in BeiDou PPP-B2b, satellite orbital and clock errors can be considered eliminated, while receiver and satellite pseudorange hardware deviations are absorbed by receiver and satellite clock biases. Tropospheric delay can be divided into dry delay and wet delay. The zenith tropospheric dry delay is generally corrected using the Saastamonien model, while the zenith tropospheric wet delay component needs to be estimated. The linearized GNSS pseudorange and carrier phase for each epoch are expressed in the observation equations under the CA model as follows:

[0036]

[0037] In the formula,

[0038] and To absorb the ambiguity caused by hardware deviations in receiver and satellite pseudorange; x and e represent the station position increment and its corresponding coefficient matrix, respectively; The receiver clock bias is the result of absorbing hardware deviations in both the receiver and satellite pseudorange; v and a represent the receiver's three-dimensional velocity and acceleration, respectively; zwd is the zenith wet delay; M w This is the projection function of the tropospheric wet component;

[0039] The above equation can be simplified to the observation equation in Kalman filtering:

[0040] L k =A k X k +ν k ,ν k~N(0,R k );

[0041] In the formula, k represents the current time, L is the observation vector, A is the coefficient matrix of the observation equation, ν is the observation noise matrix, and R is the observation noise variance matrix, which can be expressed in simplified matrix form as follows:

[0042]

[0043] In the formula, and Let represent the observation noise of the carrier phase and pseudorange, respectively. In the above observation function model, the estimated state parameters of PPPVE can be expressed as:

[0044]

[0045] Furthermore, in step two, the PPP-VE model uses a state transition matrix to link the system's velocity and acceleration with other parameters, and its state transition equation is as follows:

[0046] X k =F k,k-1 X k-1 +w k ,w k ~N(0,Q) w );

[0047] In the formula, k-1 represents the previous time step, w is the state noise matrix; F is the state transition matrix; Q w Let F and Q be the variance matrix of the state noise. w The specific expression is as follows:

[0048]

[0049]

[0050] In the formula, m t and m z These represent the power densities of the receiver clock and the zenith wet delay, respectively. Set to 100m 2 ·s -1 , Set to 10 -9 m 2 ·s -1 ;

[0051] Finally, the observation epochs are calculated using the Kalman filtering method, and the coseismic displacement of the station can be obtained by integration from the following equation:

[0052]

[0053] In the formula, s is the coseismic displacement value obtained by integration, and k0 is the initial epoch.

[0054] Another objective of this invention is to provide a real-time coseismic displacement velocity calculation system based on BeiDou PPP-B2b service, which applies the aforementioned real-time coseismic displacement velocity calculation method based on BeiDou PPP-B2b service. The real-time coseismic displacement velocity calculation system based on BeiDou PPP-B2b service includes:

[0055] The precise ephemeris calculation module is used to combine broadcast ephemeris to complete the precise ephemeris calculation method based on the BeiDou-3 PPP-B2b service;

[0056] The displacement calculation method construction module constructs a real-time coseismic displacement calculation method, including GNSS observation data preprocessing and error correction methods, PPP observation and stochastic model establishment methods, constant acceleration dynamic model construction methods, and three-dimensional displacement, velocity and acceleration estimation methods of monitoring stations.

[0057] Another object of the present invention is to provide a computer device, the computer device including a memory and a processor, the memory storing a computer program, and when the computer program is executed by the processor, causing the processor to perform the steps of the real-time coseismic displacement velocity calculation method based on BeiDou PPP-B2b service.

[0058] Another objective of this invention is to provide a computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform the steps of the real-time coseismic displacement velocity calculation method based on BeiDou PPP-B2b service.

[0059] Another objective of this invention is to provide an information data processing terminal for implementing the aforementioned real-time coseismic displacement velocity calculation system based on BeiDou PPP-B2b service.

[0060] Based on the above technical solutions and the technical problems solved, the advantages and positive effects of the technical solution to be protected by this invention are as follows:

[0061] First, this invention discloses a real-time coseismic displacement and velocity calculation method based on BeiDou PPP-B2b service. Utilizing the orbit, clock bias, and code deviation corrections broadcast by PPP-B2b service, it calculates coseismic displacement and velocity in real time, solving the problems of inconvenient communication and low accuracy in coseismic displacement and velocity calculation during earthquakes. The key technical points are: acquiring PPP-B2b message data broadcast by BeiDou-3 GEO satellites; combining this with broadcast ephemeris data to recover precise clock bias and orbit based on PPP-B2b corrections; acquiring GNSS observation data; completing data preprocessing and error correction; constructing an ionospheric-neutral PPP observation equation; applying the PPP-B2b precise clock bias and orbit to the ionospheric-neutral PPP observation equation; weighting based on satellite elevation angles; and using Kalman filtering for parameter estimation based on the PPP observation model and the satellite elevation angle stochastic model. The main focus is on establishing a constant acceleration dynamic model during the filtering process to estimate the three-dimensional displacement, velocity, and acceleration of the monitoring station.

[0062] Secondly, based on the BeiDou PPP-B2b service and combined with a constant acceleration dynamics model, this invention provides a stable, reliable, and inexpensive real-time displacement, velocity, and acceleration calculation method for GNSS seismic monitoring stations without network communication in the near field of strong earthquakes. This provides observational data for earthquake monitoring and early warning research and can be widely applied in the field of satellite measurement.

[0063] Third, as supplementary evidence of the inventive step of the claims of this invention, it is also reflected in the following important aspects:

[0064] (1) The expected benefits and commercial value of the technical solution of the present invention after transformation are as follows: The present invention can provide accurate and reliable real-time coseismic displacement velocity extraction services, and can be applied to provincial earthquake bureaus, emergency search and rescue centers and other institutions in the future, which can generate significant economic benefits.

[0065] (2) The technical solution of this invention fills a technological gap in the domestic and international industry: Since its release in 2020, PPP-B2b technology has become a research hotspot in the field of geodesy, but it has not yet been applied in the field of earthquake monitoring. This invention combines PPP-B2b positioning technology and the VADASE method to overcome the impact of earthquakes on positioning accuracy, and successfully applies PPP-B2b technology to the field of earthquake monitoring, filling the technological gap of PPP-B2b technology in earthquake monitoring.

[0066] (3) Does the technical solution of this invention solve a technical problem that people have long desired to solve but have never been able to successfully address? my country is located in one of the world's two major earthquake belts and is one of the world's most earthquake-prone areas. Earthquake monitoring and early warning are crucial for earthquake prevention and disaster reduction. However, current earthquake monitoring methods cannot overcome the risk of network interruption during strong earthquakes and the defect of decreased monitoring accuracy during earthquakes. This results in the inability to provide accurate and reliable earthquake monitoring and early warning services, hindering earthquake emergency response and decision-making. This invention uses PPP-B2b technology to obtain real-time coseismic displacement and velocity, without relying on network communication, and combines the VADASE method to improve the calculation accuracy, thus solving two major pain points in the industry. Attached Figure Description

[0067] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the embodiments of the present invention will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0068] Figure 1 This is a flowchart of the real-time coseismic displacement and velocity calculation method provided in the embodiments of the present invention;

[0069] Figure 2 This is a comparison chart of coseismic displacement estimation results and PPP displacement results in actual earthquake events provided by embodiments of the present invention;

[0070] Figure 3 This is a comparison chart of coseismic velocity estimation results and VADASE velocity measurement results in international earthquake events provided by embodiments of the present invention. Detailed Implementation

[0071] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0072] To address the problems existing in the prior art, this invention provides a method and system for real-time coseismic displacement velocity calculation based on BeiDou PPP-B2b service. The invention will be described in detail below with reference to the accompanying drawings.

[0073] To enable those skilled in the art to fully understand how the present invention is specifically implemented, this section provides an explanatory description of the embodiments that expand upon the technical solutions of the claims.

[0074] Example 1

[0075] This invention provides a method for real-time coseismic displacement velocity calculation using the BeiDou-3 PPP-B2b service, including real-time PPP-VE for near-field seismic events, where PPP-VE is a precise single-point positioning and velocity measurement technology.

[0076] Steps to achieve real-time PPP-VE in the near field of an earthquake:

[0077] The first step is to combine broadcast ephemeris data to complete the precise ephemeris calculation method based on the BeiDou-3 PPP-B2b service. This mainly includes the calculation of precise orbits and clock errors. The specific steps are as follows:

[0078] 1) Express the precise orbital corrections in the satellite orbital coordinate system as precise orbital corrections in the geocentric-geocentric coordinate system:

[0079]

[0080] Among them, [δO r δO a δO c ] T These are the orbital corrections in the radial, tangential, and normal star-fixed (ACR) coordinate systems; [δO x δO y δO z ] T This is the orbital correction for the satellite within the Earth Center Earth Fixed (ECEF) reference frame.

[0081]

[0082] In the formula, r and To calculate satellite position and satellite velocity using broadcast ephemeris.

[0083] 2) Calculate the precision orbit based on PPP-B2b service:

[0084]

[0085] Among them, [XYZ] T These are the coordinates of the satellite within the Earth Center Earth Fixed (ECEF) reference frame.

[0086] 3) Calculate the precise clock bias based on PPP-B2b service by combining the satellite clock bias provided by the broadcast ephemeris.

[0087]

[0088] in, The precision clock error recovered from the PPP-B2b signal, where C0 is the correction for the clock error in the line-of-sight direction calculated using a third-order polynomial.

[0089] The second step involves constructing a real-time coseismic displacement calculation method, which mainly includes preprocessing and error correction of GNSS observation data, establishment of PPP observation and stochastic models, construction of constant acceleration dynamic models, and estimation methods for the three-dimensional displacement, velocity, and acceleration of monitoring stations. The specific steps are as follows:

[0090] 1) In the PPP positioning model, the commonly used GNSS primitive pseudorange and phase observations can be expressed as:

[0091]

[0092] In the formula, P i and L i The original code / phase observations; ρ represents the geometric distance from the receiver to the satellite; c represents the speed of light; δt r and δt s These represent the clock difference between the receiver and the satellite, respectively; T is the tropospheric delay; f i λ is the frequency value; I1 is the ionospheric delay of L1; i For wavelength, N i b represents integer ambiguity; r,i and These represent the receiver and satellite phase hardware offsets, respectively; B r,i and These represent the hardware deviations of the receiver and the satellite pseudorange, respectively. and These represent the unmodeled errors in pseudorange and carrier phase, including thermal noise and multipath error. Furthermore, errors such as unspecified relativistic effects, Sagnac effects, Shapiro delay, phase entanglement, and tidal effects have been corrected using existing models.

[0093] 2) In the PPP-VE positioning model, a dual-frequency ionospheric de-de-spheric combination is used to eliminate the influence of the first-order ionospheric delay. The established observation equation for the dual-frequency ionospheric de-de-spheric combination can be expressed as:

[0094]

[0095] In the formula,

[0096] 3) The dynamic model of PPP-VE during filtering calculations is generally based on Newton's second law, typically using dynamic position vectors, velocities, and accelerations to describe the vehicle's trajectory. In practical solutions, finite-order differentials are often used for approximation. The constant acceleration (CA) dynamic model can be expressed as:

[0097]

[0098] In the formula, k and k-1 represent the current time and the previous time, respectively; τ is the sampling interval; w represents the state noise matrix; q a q is the power density of acceleration. a Set to 0.01 m·s -5 / 2 I is a three-dimensional identity matrix.

[0099] 4) Due to the adoption of the precise orbit and precise satellite clock bias of BeiDou PPP-B2b, the satellite's orbital and clock errors can be considered eliminated, and the hardware deviations of the receiver and satellite pseudorange are absorbed by the receiver and satellite clock biases. Furthermore, tropospheric delay can be divided into dry delay and wet delay. The zenith tropospheric dry delay is generally corrected using the Saastamonien model, while the zenith tropospheric wet delay component needs to be estimated. Therefore, the linearized GNSS pseudorange and carrier phase for each epoch can be expressed as the observation equations under the CA model as follows:

[0100]

[0101]

[0102] In the formula, and To absorb the ambiguity caused by hardware deviations in receiver and satellite pseudorange; x and e represent the station position increment and its corresponding coefficient matrix, respectively; The receiver clock bias is the result of absorbing hardware deviations in both the receiver and satellite pseudorange; v and a represent the receiver's three-dimensional velocity and acceleration, respectively; zwd is the zenith wet delay; M w This is the projection function of the tropospheric wet component.

[0103] The above equation can be simplified to the observation equation in Kalman filtering:

[0104] L k =A k X k +v k ,v k ~N(0,R k )

[0105] In the formula, k represents the current time, L is the observation vector, A is the coefficient matrix of the observation equation, ν is the observation noise matrix, and R is the observation noise variance matrix, which can be expressed in simplified matrix form as follows:

[0106]

[0107] In the formula, and Let represent the observation noise of the carrier phase and pseudorange, respectively. In the above observation function model, the estimated state parameters of PPPVE can be expressed as:

[0108]

[0109] The PPP-VE model relates the system's velocity and acceleration to other parameters through a state transition matrix. Its state transition equation is as follows:

[0110] X k =F k,k-1 X k-1 +w k ,w k ~N(0,Q) w )

[0111] In the formula, k-1 represents the previous time step, w is the state noise matrix; F is the state transition matrix; Q w Let F and Q be the variance matrix of the state noise. w The specific expression is as follows:

[0112]

[0113] In the formula, m t and m z These represent the power densities of the receiver clock and the zenith wet delay, respectively. Set to 100m 2 ·s -1 , Set to 10 -9 m 2 ·s -1

[0114] Using the Kalman filter method to calculate the observation epoch, the coseismic displacement of the station can be obtained by integration from the following two equations:

[0115]

[0116] In the formula, s is the coseismic displacement value obtained by integration, and k0 is the initial epoch.

[0117] Compared with existing technologies, the method for real-time coseismic displacement and velocity calculation using BeiDou-3 PPP-B2b service provided in this invention has at least the following beneficial effects:

[0118] First, the method for real-time coseismic displacement and velocity calculation using BeiDou-3 PPP-B2b service provided in this embodiment of the invention can calculate high-precision displacement and velocity sequences in real time when the strong earthquake network is interrupted, and can be applied to the field of earthquake monitoring and early warning.

[0119] Secondly, the PPP-B2b service in the method for real-time coseismic displacement and velocity calculation using BeiDou-3 PPP-B2b service provided in this embodiment of the invention is publicly available free of charge. It can calculate coseismic displacement and velocity in real time without broadcasting RTS service through commercial satellite communication services, thus reducing the burden of GNSS station deployment.

[0120] Third, the method for real-time coseismic displacement and velocity calculation using BeiDou-3 PPP-B2b service provided in this embodiment of the invention directly uses a dynamic model to estimate the velocity and acceleration of the Earth's surface, and almost eliminates the drift problem that occurs in the coseismic displacement calculated by the VADASE method.

[0121] Example 2

[0122] The following examples illustrate the application scenarios of the GNSS coseismic displacement, velocity, and acceleration calculation method based on BeiDou-3 PPP-B2b service:

[0123] This invention, based on the PPP-B2b service of BeiDou-3, constructs a real-time coseismic displacement, velocity, and acceleration extraction method based on the PPP-VE model. To evaluate the performance of the real-time displacement, velocity, and acceleration calculation method based on the PPP-B2b service, VADASE velocity measurement results based on broadcast ephemeris were calculated for comparison. The post-processed PPP displacement results based on the final ephemeris product (GBM) of the multi-system system released by the GFZ center were used as reference values ​​to discuss the advantages and disadvantages of this invention (Table 2).

[0124] Table 2 Validation Scheme

[0125]

[0126] An application embodiment of the present invention provides a computer device, which includes a memory and a processor. The memory stores a computer program. When the computer program is executed by the processor, the processor performs the steps of the real-time coseismic displacement velocity calculation method based on BeiDou PPP-B2b service.

[0127] An application embodiment of the present invention provides a computer-readable storage medium storing a computer program. When the computer program is executed by a processor, the processor performs the steps of the real-time coseismic displacement velocity calculation method based on BeiDou PPP-B2b service.

[0128] An application embodiment of the present invention provides an information data processing terminal, which is used to implement the real-time coseismic displacement velocity calculation system based on BeiDou PPP-B2b service.

[0129] As a preferred embodiment, this invention discloses a real-time coseismic displacement and velocity calculation method based on BeiDou PPP-B2b service. Utilizing the orbit, clock bias, and code deviation corrections broadcast by PPP-B2b service, it calculates coseismic displacement and velocity in real time, solving the problems of inconvenient communication and low accuracy in coseismic displacement and velocity calculation during earthquakes. The key technical points are: acquiring PPP-B2b message data broadcast by BeiDou-3 GEO satellites; combining this with broadcast ephemeris data to recover precise clock bias and orbit based on PPP-B2b corrections; acquiring GNSS observation data; completing data preprocessing and error correction; constructing an ionospheric-neutral PPP observation equation; applying the PPP-B2b precise clock bias and orbit to the ionospheric-neutral PPP observation equation; weighting based on satellite elevation angles to construct a stochastic model; and using Kalman filtering for parameter estimation based on the PPP observation model and the satellite elevation angle stochastic model. The main focus is on establishing a constant acceleration dynamic model during the filtering process to estimate the three-dimensional displacement, velocity, and acceleration of the monitoring station. Based on the BeiDou PPP-B2b service and combined with a constant acceleration dynamics model, this invention provides a stable, reliable, and inexpensive real-time displacement, velocity, and acceleration calculation method for GNSS seismic monitoring stations without network communication in the near field of strong earthquakes. It provides observational data for earthquake monitoring and early warning research and can be widely applied in the field of satellite measurement.

[0130] The embodiments of the present invention have achieved some positive results during the research and development or use process, and have indeed great advantages compared with the prior art. The following content describes them in conjunction with the data, charts and other information of the experimental process.

[0131] 1) Static Experiment

[0132] The static experiment selected 10 IGS stations distributed in the surrounding areas of China. The experimental data used were dual-frequency data from both GPS and BDS systems on March 17, 2022 (76th day of the year), with a sampling rate of 1 second. During the calculation, the GNSS data for one day was divided into a calculation arc every 4 hours. Assuming the stations remained stationary, with 0 as the reference value, the velocity and displacement positioning obtained by the VADASE and PPP-VE methods were evaluated. The convergence criterion for PPP-VE was set as the displacement change between two consecutive epochs being less than 0.01 m, and remaining within 0.01 m. Since PPP-VE has a convergence time of several minutes, the velocity and displacement time series 15 minutes after the start of each calculation arc were selected for accuracy evaluation, with an evaluation time length of 15 minutes and a total of 900 epochs. Because the mainshock rupture time of an earthquake is usually less than 15 minutes, the velocity within each four-hour calculation arc was divided into five-minute intervals, and the displacement integration restarted every 15 minutes. Table 2 presents the RMS values ​​of velocity and displacement results for the VADASE and PPP-VE methods in four dimensions: north-south, east-west, vertical, and 3D.

[0133] Table 2. Statistical Table of Static Experiment Accuracy

[0134]

[0135] Experimental results: As can be seen from the table, the velocity measurement accuracy of the VADASE method and the PPP-VE method is comparable, and is better than 5 mm / s in the north-south, east-west and vertical directions. However, the displacement extraction accuracy of the PPP-VE method is significantly improved compared to the VADASE method, with improvements of 60%, 70% and 66% in the north-south, east-west and vertical directions, respectively.

[0136] 2) Actual earthquake events

[0137] According to the China Earthquake Networks Center, a magnitude 7.4 earthquake occurred in Madoi County, Golog Tibetan Autonomous Prefecture, Qinghai Province, China (34.59°N, 98.34°E) at 02:04:11 Beijing time on May 22, 2021 (18:04:13 UTC on May 21, 2021), with a focal depth of 17 kilometers. This earthquake was successfully recorded by 43 Qinghai Continuously Operating Reference Stations (QHCORS) and 16 stations of the China Crustal Movement Observation Network (CMONOC). This experiment selected high-frequency GNSS data from 58 stations, including Madoi (QHMD), for three hours before and after the earthquake. The coseismic displacement and velocity of all stations were extracted using the PPP-VE and VADASE methods. The displacement results calculated by post-processing PPP were used as reference values ​​for accuracy evaluation. Figure 2 and Figure 3 Coseismic displacement and velocity images extracted from typical QHMD stations are presented. It can be seen that the velocity sequences obtained by the VADASE method and the PPP-VE method are generally consistent in their overall trends, both capturing significant seismic signals. The difference between the two is less than 1 mm / s. The PPP-VE method shows a significant improvement in displacement extraction accuracy compared to the VADASE method, without exhibiting obvious drift. The displacement extraction accuracies in the north-south, east-west, and vertical directions are 2.61 cm, 1.52 cm, and 2.97 cm, respectively.

[0138] It should be noted that embodiments of the present invention can be implemented in hardware, software, or a combination of both. The hardware portion can be implemented using dedicated logic; the software portion can be stored in memory and executed by a suitable instruction execution system, such as a microprocessor or dedicated-design hardware. Those skilled in the art will understand that the above-described devices and methods can be implemented using computer-executable instructions and / or included in processor control code, for example, such code provided on a carrier medium such as a disk, CD, or DVD-ROM, a programmable memory such as read-only memory (firmware), or a data carrier such as an optical or electronic signal carrier. The devices and modules of the present invention can be implemented by hardware circuitry such as very large-scale integrated circuits or gate arrays, semiconductors such as logic chips, transistors, or programmable hardware devices such as field-programmable gate arrays, programmable logic devices, etc., or by software executed by various types of processors, or by a combination of the above-described hardware circuitry and software, such as firmware.

[0139] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications, equivalent substitutions, and improvements made by those skilled in the art within the scope of the technology disclosed in the present invention, and within the spirit and principles of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method for real-time coseismic displacement velocity calculation based on BeiDou PPP-B2b service, characterized in that, include: Step 1: Combine broadcast ephemeris to complete the precise ephemeris calculation method based on BeiDou-3 PPP-B2b service, including the calculation of precise orbit and clock bias; Step 2: Construct a real-time coseismic displacement calculation method, including preprocessing and error correction of GNSS observation data, establishment of PPP observation and stochastic model, construction of constant acceleration dynamic model, and estimation methods for three-dimensional displacement, velocity and acceleration of monitoring stations. In step one, in order to ensure the correlation between the information content broadcast by different information types, as well as the correlation between the correction number information and the broadcast ephemeris, the information is identified by four version numbers: Spatial State Expression Data Period Number (IOD SSR), Pseudo-Random Noise Code Data Period Number (IODP), Navigation Data Period Number (IODN), and Correction Data Identification and Data Period Number (IOD Corr), to facilitate matching and use. The orbit correction data in the PPP-B2b message has an epoch time interval of 48 seconds, and the broadcast orbit correction information includes the orbit correction vector. The satellite position correction vector can be calculated using the orbit correction values. After the IODN of the orbit correction information is successfully matched with the IODE in the wide-interpolation ephemeris navigation message, the precise orbit correction data in the satellite orbit coordinate system can be expressed as the precise orbit correction data in the geocentric-ground-fixed coordinate system. ; in, These are the orbital corrections in radial, tangential, and normal star-fixed coordinate systems; For the orbital correction of the satellite in the geocentric-solid reference frame; The unit vectors corresponding to the radial, tangential, and normal directions, respectively, can be represented as: ； In the formula, and To calculate satellite position and satellite velocity using broadcast ephemeris; Calculations are performed on precision orbits based on PPP-B2b services: ； in, The coordinates of the satellite within a geocentric, geofixed reference frame; In step one, the epoch time interval of the clock correction number in the PPP-B2b clock correction information is 6 seconds; after the clock correction information is successfully matched with the broadcast ephemeris, the clock correction parameters contained in the message can be used. Clock bias parameters obtained from broadcast ephemeris calculations Corrections are made to calculate the precision clock bias based on PPP-B2b services: ； in, Precision clock error recovered from PPP-B2b signal, The clock error correction in the line-of-sight direction is calculated using a third-order polynomial. When the track correction and clock error correction information in the PPP-B2b message are the same, and their IOD Corr can be matched with the IODE in the navigation message, the track and clock error correction can be used to correct errors. In step two, the dynamic model of PPP-VE during filtering calculation is based on Newton's second law, using dynamic position vectors, velocities, and accelerations to describe the trajectory of the carrier. In the actual solution, finite-order differentials are used for approximation, and the constant acceleration CA dynamic model can be expressed as: ； In the formula, These represent the current time and the previous time, respectively. The sampling interval; Represents the state noise matrix; The power density of acceleration, Set to 0.01 ; Because of the use of precise orbits and precise satellite clock biases in BeiDou PPP-B2b, satellite orbital and clock errors can be considered eliminated, while receiver and satellite pseudorange hardware biases are absorbed by receiver and satellite clock biases. Tropospheric delay can be divided into dry delay and wet delay. The zenith tropospheric dry delay is corrected using the Saastamonien model, while the zenith tropospheric wet delay component needs to be estimated. The linearized GNSS pseudorange and carrier phase for each epoch are expressed in the observation equations under the CA model as follows: ； In the formula, ; and To absorb the ambiguity caused by hardware deviations in the receiver and satellite pseudorange; and These represent the station location increments and their corresponding coefficient matrices, respectively. To absorb the receiver clock bias after the hardware deviation between the receiver and the satellite pseudorange; and These represent the receiver's three-dimensional velocity and acceleration, respectively. For zenith moisture delay; This is the projection function of the tropospheric wet component; The above equation can be simplified to the observation equation in Kalman filtering: ； In the formula, They represent the current time, For the observation vector, The coefficient matrix of the observation equation, For the observation noise matrix, The observation noise variance matrix can be represented by a simplified matrix as follows: ； In the formula, and Let represent the observation noise of the carrier phase and pseudorange, respectively. In the above observation function model, the estimated state parameters of PPPVE can be expressed as: ； In step two, the PPP-VE model uses a state transition matrix to link the system's velocity and acceleration with other parameters. Its state transition equation is as follows: ； In the formula, Indicates the previous moment, This is the state noise matrix; This is the state transition matrix; Let be the variance matrix of the state noise; and The specific expression is as follows: ； ； In the formula, It is a three-dimensional identity matrix; and These represent the power densities of the receiver clock and the zenith wet delay, respectively. Set to 100 , Set as ; Finally, the observation epochs are calculated using the Kalman filtering method, and the coseismic displacement of the station can be obtained by integration from the following equation: ； In the formula, The coseismic displacement value is obtained by integration. This is the initial epoch.

2. The real-time coseismic displacement velocity calculation method based on BeiDou PPP-B2b service as described in claim 1, characterized in that, In step two, the raw GNSS pseudorange and phase observations used in the PPP positioning model can be expressed as follows: ； In the formula, and These are the original code / phase observations; Indicates the geometric distance from the receiver to the satellite; Represents the speed of light; and These are the clock differences between the receiver and the satellite, respectively. For tropospheric delay; This is the frequency value; for Ionospheric delay; For wavelength, For integer ambiguity; and These represent the receiver and satellite phase hardware offsets, respectively. and These represent the hardware deviations of the receiver and the satellite pseudorange, respectively. and These represent the unmodeled errors in pseudorange and carrier phase, including thermal noise and multipath error, respectively. In addition, unspecified relativistic effects, Sagnac effects, Shapiro delay, phase entanglement, and tidal effects have been corrected using existing models. In the PPP-VE positioning model, a dual-frequency ionospheric de-de-spheric combination is used to eliminate the influence of the first-order ionospheric delay. The established observation equation for the dual-frequency ionospheric de-de-spheric combination can be expressed as: ； In the formula, ; ; ; , ; The weighting is performed using a satellite elevation angle stochastic model, expressed as follows: ； In the formula, For the calculated observation weight matrix, The observation weighting factor for each satellite. For satellite The satellite elevation angle.

3. A real-time coseismic displacement-velocity calculation system based on BeiDou PPP-B2b service, which applies the aforementioned real-time coseismic displacement-velocity calculation method, is characterized in that... The real-time coseismic displacement velocity calculation system based on BeiDou PPP-B2b service includes: The precise ephemeris calculation module is used to combine broadcast ephemeris to complete the precise ephemeris calculation method based on the BeiDou-3 PPP-B2b service; The displacement calculation method construction module constructs a real-time coseismic displacement calculation method, including GNSS observation data preprocessing and error correction methods, PPP observation and stochastic model establishment methods, constant acceleration dynamic model construction methods, and three-dimensional displacement, velocity and acceleration estimation methods of monitoring stations.

4. A computer device, comprising a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor performs the steps of the real-time coseismic displacement velocity calculation method based on BeiDou PPP-B2b service as described in claims 1-2.

5. A computer-readable storage medium storing a computer program, wherein when the computer program is executed by a processor, the processor performs the steps of the real-time coseismic displacement velocity calculation method based on BeiDou PPP-B2b service as described in claims 1-2.

6. An information data processing terminal, the information data processing terminal being used to implement the real-time coseismic displacement velocity calculation system based on BeiDou PPP-B2b service as described in claim 3.

Images