Real-time reconstruction method of broadband coseismic velocity waveform based on fusion data of cascaded filtering

Abstract

The application provides a fusion data broadband coseismic velocity waveform real-time reconstruction method based on cascade filtering, which uses cascade filtering to fuse high-frequency GNSS original Doppler data, high-frequency GNSS carrier phase data and accelerometer data to realize broadband coseismic velocity waveform real-time reconstruction, and the realization process comprises the following steps: high-frequency GNSS / accelerometer co-located station real-time acquisition of high-frequency GNSS observation data, broadcast ephemeris and original acceleration data, and sending to a ground computing center; the ground computing center real-time calculation of high-frequency GNSS TDCP average velocity and high-frequency GNSS RD instantaneous velocity; the ground computing center zero offset correction of the real-time received original acceleration data to obtain zero offset correction acceleration data; and the ground computing center real-time reconstruction of the broadband coseismic velocity waveform of the high-frequency GNSS / accelerometer co-located station according to the high-frequency GNSS TDCP average velocity, the high-frequency GNSS RD instantaneous velocity and the zero offset correction acceleration data by using cascade filtering.

Description

Technical Field

[0001] This invention relates to the field of geological disaster monitoring technology, and in particular to a method for real-time reconstruction of broadband coseismic velocity waveforms based on cascaded filtering and fusing high-frequency GNSS and accelerometer data. Background Technology

[0002] High-frequency GNSS can measure high-precision low-frequency seismic signals, even permanent coseismic displacements, but its data sampling rate is low, typically only 1–5 Hz. Accelerometers can obtain seismic acceleration information with high sampling rates (generally as high as 100–200 Hz), but the integrated velocity or displacement sequences usually exhibit nonlinear trend drift, making it difficult to obtain accurate low-frequency seismic information. The two sensors have excellent complementary characteristics; by fusing high-frequency GNSS and accelerometer data, broadband seismic waveform reconstruction can be achieved.

[0003] Currently, the high-frequency GNSS data processing portion of high-frequency GNSS / accelerometer fusion algorithms generally includes precise single-point positioning, relative positioning, and carrier phase epoch differential velocimetry. Precise single-point positioning requires a relatively long ambiguity convergence process to obtain centimeter-level displacement, while relative positioning accuracy heavily depends on the reference station. Stable reference stations are often difficult to find in the near field of a large earthquake, and the accuracy of relative positioning results decreases with increasing baseline length. High-frequency GNSS carrier phase epoch differential velocimetry avoids ambiguity convergence or reconvergence processes by differentiating carrier phase observations from adjacent epochs; it can achieve real-time continuous reconstruction of coseismic velocity and displacement waveforms using only real-time available broadcast ephemeris and single-station operation. However, the above three high-frequency GNSS data processing algorithms, because they only use carrier phase and pseudorange observations, can only reconstruct the coseismic average velocity waveform (essentially the displacement increment between epochs). In strong-motion deformation monitoring applications such as earthquakes, there is a significant difference between average velocity and instantaneous velocity; therefore, further exploration is needed to discover methods that combine high-frequency GNSS and accelerometers to capture broadband coseismic velocity waveforms. Summary of the Invention

[0004] In static scenarios, based on 1Hz high-frequency GNSS observation data, the carrier phase epoch differential technique (TDCP) can obtain average velocities with accuracy on the order of mm / s, while the original Doppler instantaneous velocity method (RD) can only measure instantaneous velocities with accuracy on the order of cm / s. Compared with TDCP, although RD has greater velocity measurement noise, the estimated velocity waveform has a higher signal-to-noise ratio in strong vibration scenarios. This invention discloses a real-time reconstruction method for broadband coseismic velocity waveforms based on cascaded filtering and fusing high-frequency GNSS and accelerometer data. This invention adaptively selects the optimal velocity from TDCP average velocity and RD instantaneous velocity using a priori threshold filtering to obtain a high signal-to-noise ratio high-frequency GNSS coseismic velocity waveform. Furthermore, this invention also uses dual-speed Kalman filtering to fuse high-frequency GNSS velocity with accelerometer data, broadening the frequency band of the coseismic velocity waveform. In summary, by using cascaded filtering to fuse high-frequency GNSS TDCP average velocity, RD instantaneous velocity, and accelerometer data, this invention can reconstruct high signal-to-noise ratio broadband coseismic velocity waveforms in real time, providing reliable ground motion parameters for the rapid production of seismic intensity maps.

[0005] This invention provides a method for real-time reconstruction of broadband coseismic velocity waveforms based on cascaded filtering. It utilizes cascaded filtering to fuse raw high-frequency GNSS Doppler data, high-frequency GNSS carrier phase data, and accelerometer data to achieve real-time reconstruction of broadband coseismic velocity waveforms. The implementation process is as follows.

[0006] The high-frequency GNSS / accelerometer co-located station collects high-frequency GNSS observation data, broadcast ephemeris and raw acceleration data in real time, and sends the collected high-frequency GNSS observation data, broadcast ephemeris and raw acceleration data to the ground computing center;

[0007] The ground computing center calculates the average velocity of high-frequency GNSS TDCP and the instantaneous velocity of high-frequency GNSS RD in real time based on the high-frequency GNSS observation data received in real time and the broadcast ephemeris.

[0008] The ground computing center performs pre-earthquake bias correction on the raw acceleration data received in real time to obtain bias-corrected acceleration data.

[0009] The ground computing center reconstructs the broadband coseismic velocity waveform of the high-frequency GNSS / accelerometer co-located station in real time using cascaded filtering based on the high-frequency GNSS TDCP average velocity, the high-frequency GNSS RD instantaneous velocity, and the zero-bias correction acceleration data.

[0010] Furthermore, the ground computing center calculates the average velocity of the high-frequency GNSS TDCP and the instantaneous velocity of the high-frequency GNSS RD in real time based on the high-frequency GNSS observation data and the broadcast ephemeris, including the following processing:

[0011] The high-frequency GNSS observation data includes high-frequency GNSS pseudorange observations, high-frequency GNSS carrier phase observations, and high-frequency GNSS raw Doppler observations;

[0012] The ground computing center calculates the approximate coordinates of the high-frequency GNSS / accelerometer co-located station using standard single-point positioning based on the high-frequency GNSS pseudorange observations and the broadcast ephemeris.

[0013] The ground computing center calculates the average velocity of the high-frequency GNSS TDCP of the high-frequency GNSS / accelerometer co-located station using the high-frequency GNSS carrier phase observation value, the broadcast ephemeris and the approximate coordinates, and the high-frequency GNSS carrier phase epoch difference technology.

[0014] The ground computing center calculates the high-frequency GNSS RD instantaneous velocity of the high-frequency GNSS / accelerometer co-located station using the high-frequency GNSS raw Doppler observations, the broadcast ephemeris, and the approximate coordinates, employing the high-frequency GNSS raw Doppler instantaneous velocity method.

[0015] Furthermore, the ground computing center, based on the high-frequency GNSS TDCP average velocity, the high-frequency GNSS RD instantaneous velocity, and the zero-bias correction acceleration data, uses cascaded filtering to reconstruct the broadband coseismic velocity waveform of the high-frequency GNSS / accelerometer co-located station in real time, including the following processing.

[0016] The cascaded filtering includes a prior threshold filter and a dual-speed Kalman filter;

[0017] The ground computing center adaptively selects the optimal high-frequency GNSS velocity using prior threshold filtering based on the average velocity of the high-frequency GNSS TDCP and the instantaneous velocity of the high-frequency GNSS RD.

[0018] The ground computing center reconstructs the broadband coseismic velocity waveform of the high-frequency GNSS / accelerometer co-located station in real time using dual-speed Kalman filtering, based on the optimal high-frequency GNSS velocity and the zero-bias correction acceleration data.

[0019] Furthermore, the ground computing center adaptively selects the optimal high-frequency GNSS velocity based on the average velocity of the high-frequency GNSS TDCP and the instantaneous velocity of the high-frequency GNSS RD using prior threshold filtering, as follows:

[0020] The ground computing center calculates the high-frequency GNSS TDCP / RD cross-difference sequence based on the high-frequency GNSS TDCP average velocity sequence and the high-frequency GNSS RD instantaneous velocity sequence, as follows:

[0021]

[0022] In the formula, This represents the cross-difference between high-frequency GNSS TDCP / RD at epoch k. Represents the average velocity of high-frequency GNSS TDCP at epoch k; This represents the instantaneous velocity of the high-frequency GNSS RD at epoch k; abs() indicates taking the absolute value.

[0023] The ground computing center calculates the root mean square error of the high-frequency GNSS TDCP / RD cross-difference sequence before the earthquake event based on the high-frequency GNSS TDCP / RD cross-difference sequence as follows:

[0024]

[0025] In the formula, rms() represents the root mean square error before the seismic event calculated from the cross-difference sequence of high-frequency GNSS TDCP / RD from epoch k0 to epoch k1; k2 represents the P-wave arrival time of the high-frequency GNSS / accelerometer co-located station.

[0026] The ground computing center adaptively selects the optimal high-frequency GNSS velocity from the high-frequency GNSS TDCP average velocity and the high-frequency GNSS RD instantaneous velocity by setting a prior threshold based on the root mean square error before the earthquake event and the cross-difference sequence of the high-frequency GNSS TDCP / RD.

[0027] Moreover, the values ​​of k0 and k1 are k1 = k2 and k0 = k1 - 300s.

[0028] Furthermore, the adaptive selection of the optimal high-frequency GNSS velocity from the high-frequency GNSS TDCP average velocity and the high-frequency GNSS RD instantaneous velocity is achieved as follows:

[0029]

[0030] In the formula, This represents the optimal high-frequency GNSS velocity adaptively selected using prior threshold filtering; ratio is the prior ratio threshold.

[0031] Furthermore, the ratio value is set to 6.

[0032] Furthermore, the ground computing center, based on the optimal high-frequency GNSS velocity and the zero-bias correction acceleration data, uses dual-speed Kalman filtering to reconstruct the broadband co-seismic velocity waveform of the high-frequency GNSS / accelerometer co-located station in real time. The dual-speed Kalman filtering fusion method employed is as follows.

[0033] Let the state variable x of the dual-speed Kalman filter be... k For three-dimensional velocity v k ,

[0034] x k =[v e,k v n,k v u,k ] T

[0035] In the formula, x k v is the state variable for epoch k; e,k v n,k and v u,k These are the east-west, north-south, and vertical components of the velocity at epoch k, respectively; T Indicates the transpose operation;

[0036] The dual-speed Kalman filter system model is

[0037] x k =A k-1,k x k-1 +B k-1,k u k-1,k +w k-1,k ,w k-1,k ~N(0,Q) k-1,k )

[0038] A k-1.k =I 3×3 B k-1,k =τ a I 3×3

[0039] u k-1,k =(a k-1 +a k ) / 2

[0040] Q k-1,k =q a τ a I 3×3

[0041] In the formula, A k-1,k Let x be the state transition matrix from epoch k-1 to epoch k. k-1 For k-1 epoch state variables; u k-1,k For system input; a k-1 and a kThe zero-bias correction accelerations for epochs k-1 and k, respectively; B k-1,k Input matrix for the system; w k-1,k The system noise, w, follows a normal distribution. k-1,k ~N(0,Q) k-1,k ), Q k-1,k τ is the variance-covariance matrix of the system noise; a q represents the sampling interval of the accelerometer. a For the noise variance of the zero-bias correction acceleration, I 3×3 It is a 3×3 identity matrix;

[0042] The dual-speed Kalman filter observation model is

[0043]

[0044] H k =I k

[0045]

[0046] In the formula, For high-frequency GNSS, the optimal speed is ε. k The optimal velocity noise for high-frequency GNSS follows a normal distribution ε k ~(0,R k ), R k H is the variance-covariance matrix of the velocity noise; k Design a matrix for the observation model; and These are the noise variances of the east-west, north-south, and vertical components of the optimal velocity for high-frequency GNSS, respectively.

[0047] The dual-speed Kalman filter consists of two parts: time update and measurement update. The time update frequency is consistent with the accelerometer sampling interval, while the measurement update frequency is synchronized with the high-frequency GNSS sampling interval.

[0048] The time update formula is shown below.

[0049]

[0050] The measurement update formula is shown below.

[0051]

[0052] In the formula, () -1 P represents the matrix inversion operation; k For x k The variance and covariance matrix; (-) and (+) represent the time update and measurement update results, respectively.

[0053] This invention discloses a real-time reconstruction method for broadband coseismic velocity waveforms based on cascaded filtering and fusing high-frequency GNSS and accelerometer data. The invention utilizes prior threshold filtering to adaptively select the optimal velocity from the TDCP average velocity and RD instantaneous velocity to obtain a high signal-to-noise ratio (SNR) high-frequency GNSS coseismic velocity waveform. Furthermore, the invention employs dual-speed Kalman filtering to fuse high-frequency GNSS velocity with accelerometer data, broadening the frequency band of the coseismic velocity waveform. In summary, this invention, through cascaded filtering, fuses high-frequency GNSS TDCP average velocity, RD instantaneous velocity, and accelerometer data, enabling real-time reconstruction of high SNR broadband coseismic velocity waveforms, providing reliable ground motion parameters for the rapid production of seismic intensity maps. The solution is simple and convenient to implement, highly practical and universal, solving the problems of low practicality and inconvenience in actual application of related technologies, improving user experience, and possessing significant market value. Attached Figure Description

[0054] Figure 1 This is a flowchart illustrating the overall process of a real-time reconstruction method for broadband coseismic velocity waveforms based on cascaded filtering and fusing high-frequency GNSS and accelerometer data, according to an embodiment of the present invention.

[0055] Figure 2 This is a flowchart illustrating the calculation of average velocity and instantaneous velocity of high-frequency GNSS TDCP based on high-frequency GNSS data according to an embodiment of the present invention.

[0056] Figure 3 This is a flowchart illustrating an embodiment of the present invention that utilizes cascaded filtering to fuse high-frequency GNSS TDCP average velocity, RD instantaneous velocity, and accelerometer data. Detailed Implementation

[0057] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings and embodiments. Although exemplary embodiments of the present invention are shown in the drawings, it should be understood that the present invention can be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that the present invention can be thoroughly understood and its scope can be fully conveyed to those skilled in the art.

[0058] It will be understood by those skilled in the art that, unless otherwise defined, all terms used herein (including technical and scientific terms) have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. It should also be understood that terms such as those defined in general dictionaries should be understood to have the meaning consistent with their meaning in the context of the prior art, and should not be interpreted in an idealized or overly formal sense unless specifically defined.

[0059] In static scenarios, based on 1Hz high-frequency GNSS observation data, the carrier phase epoch differential technique (TDCP) can obtain average velocities with accuracy on the order of mm / s, while the original Doppler instantaneous velocity method (RD) can only measure instantaneous velocities with accuracy on the order of cm / s. Compared with TDCP, although RD has greater velocity measurement noise, the estimated velocity waveform has a higher signal-to-noise ratio in strong vibration scenarios. This invention discloses a real-time reconstruction method for broadband coseismic velocity waveforms based on cascaded filtering and fusing high-frequency GNSS and accelerometer data. This invention adaptively selects the optimal velocity from TDCP average velocity and RD instantaneous velocity using a priori threshold filtering to obtain a high signal-to-noise ratio high-frequency GNSS coseismic velocity waveform. Furthermore, this invention also uses dual-speed Kalman filtering to fuse high-frequency GNSS velocity with accelerometer data, broadening the frequency band of the coseismic velocity waveform. In summary, by using cascaded filtering to fuse high-frequency GNSS TDCP average velocity, RD instantaneous velocity, and accelerometer data, this invention can reconstruct high signal-to-noise ratio broadband coseismic velocity waveforms in real time, providing reliable ground motion parameters for the rapid production of seismic intensity maps.

[0060] Figure 1 This is a flowchart illustrating the overall process of a real-time reconstruction method for broadband coseismic velocity waveforms based on cascaded filtering and fusing high-frequency GNSS and accelerometer data, provided in an embodiment of the present invention. (Refer to...) Figure 1 The specific steps of the embodiment are as follows:

[0061] S11. The high-frequency GNSS / accelerometer co-located station collects high-frequency GNSS observation data, broadcast ephemeris and raw acceleration data in real time, and sends the collected high-frequency GNSS observation data, broadcast ephemeris and raw acceleration data to the ground computing center;

[0062] The aforementioned co-located high-frequency GNSS / accelerometer stations refer to a pair of high-frequency GNSS observation stations and accelerometer stations with a station spacing of less than an empirical threshold; wherein the high-frequency GNSS / accelerometer station spacing refers to the spherical distance between the high-frequency GNSS observation station and the accelerometer station, and the specific formula for calculating the spherical distance is as follows:

[0063] sepdist=R earth ×arccos(sin(B gnss )×sin(B acc )+cos(B gnss )×cos(B acc )×cos(L gnss -L acc )) (1)

[0064] In the formula, sepdist represents the spherical distance between the high-frequency GNSS observation station and the accelerometer station, i.e., the high-frequency GNSS / accelerometer station spacing; R earth B represents the average radius of the Earth, which is generally taken as 6371 km; gnss and L gnss B represents the latitude and longitude coordinates of the high-frequency GNSS observation station, respectively; acc and L acc These represent the latitude and longitude coordinates of the accelerometer station, respectively.

[0065] Furthermore, in S11, the determination of the co-located high-frequency GNSS / accelerometer stations, i.e., the selection of the empirical threshold, should follow the following criteria:

[0066] Extensive data studies have shown that the seismic waveforms recorded by high-frequency GNSS observation stations and accelerometer stations with a station spacing of less than 4 km have good consistency. Therefore, the empirical threshold can be set to 4 km. At the same time, since high-frequency seismic signals in the near field of an earthquake are more prone to rapid attenuation, the empirical threshold for selecting co-located stations should be smaller the closer they are to the source region.

[0067] S12. The ground computing center calculates the average velocity of high-frequency GNSS TDCP and the instantaneous velocity of high-frequency GNSS RD in real time based on the high-frequency GNSS observation data received in real time and the broadcast ephemeris.

[0068] Further, in S12, the ground computing center calculates the average velocity of the high-frequency GNSS TDCP and the instantaneous velocity of the high-frequency GNSS RD in real time based on the high-frequency GNSS observation data received in real time and the broadcast ephemeris, such as... Figure 2 As shown, this can be achieved through the following specific steps:

[0069] S121: The high-frequency GNSS observation data includes high-frequency GNSS pseudorange observations, high-frequency GNSS carrier phase observations, and high-frequency GNSS raw Doppler observations;

[0070] S122: The ground computing center calculates the approximate coordinates of the high-frequency GNSS / accelerometer co-located station using standard single-point positioning based on the high-frequency GNSS pseudorange observations and the broadcast ephemeris.

[0071] Standard single-point positioning is existing technology and will not be described in detail in this invention.

[0072] S123: The ground computing center calculates the average velocity of the high-frequency GNSS TDCP of the high-frequency GNSS / accelerometer co-located station using the high-frequency GNSS carrier phase observation value, the broadcast ephemeris and the approximate coordinates, and the high-frequency GNSS carrier phase epoch difference technology.

[0073] Among them, the high-frequency GNSS carrier phase epoch differential technology is existing technology and will not be described in detail in this invention.

[0074] S124: The ground computing center calculates the high-frequency GNSS RD instantaneous velocity of the high-frequency GNSS / accelerometer co-located station using the high-frequency GNSS raw Doppler observations, the broadcast ephemeris, and the approximate coordinates, employing the high-frequency GNSS raw Doppler instantaneous velocity method.

[0075] The high-frequency GNSS raw Doppler instantaneous velocity method is existing technology and will not be described in detail in this invention.

[0076] S13. The ground computing center performs pre-earthquake bias correction on the raw acceleration data received in real time to obtain bias-corrected acceleration data.

[0077] S14. The ground computing center reconstructs the broadband coseismic velocity waveform of the high-frequency GNSS / accelerometer co-located station in real time using cascaded filtering based on the high-frequency GNSS TDCP average velocity, the high-frequency GNSS RD instantaneous velocity, and the zero-bias correction acceleration data.

[0078] Further, in S14, the ground computing center, based on the high-frequency GNSS TDCP average velocity, the high-frequency GNSS RD instantaneous velocity, and the zero-bias correction acceleration data, uses cascaded filtering to reconstruct the broadband coseismic velocity waveform of the high-frequency GNSS / accelerometer co-located station in real time, such as... Figure 3 As shown, this can be achieved through the following specific steps:

[0079] S141: Set up a cascaded filter, which includes a prior threshold filter and a dual-speed Kalman filter;

[0080] S142: The ground computing center adaptively selects the optimal high-frequency GNSS velocity based on the instantaneous velocity of the high-frequency GNSS RD and the average velocity of the high-frequency GNSS TDCP using prior threshold filtering;

[0081] Furthermore, in step S142, the ground computing center adaptively selects the optimal high-frequency GNSS velocity based on the instantaneous velocity of the high-frequency GNSS RD and the average velocity of the high-frequency GNSS TDCP using prior threshold filtering, through the following specific process:

[0082] Step 1: The ground computing center calculates the high-frequency GNSS TDCP / RD cross-difference sequence based on the high-frequency GNSS TDCP average velocity sequence and the high-frequency GNSS RD instantaneous velocity sequence. The specific formula is as follows:

[0083]

[0084] In the formula, This represents the high-frequency GNSS TDCP / RD cross difference at epoch k, i.e., the cross difference between the high-frequency GNSS TDCP average velocity and the high-frequency GNSS RD instantaneous velocity. Represents the average velocity of high-frequency GNSS TDCP at epoch k; This represents the instantaneous velocity of the high-frequency GNSS RD at epoch k; abs() indicates taking the absolute value.

[0085] Step 2: The ground computing center calculates the root mean square error of the high-frequency GNSS TDCP / RD cross-difference sequence before the earthquake event based on the high-frequency GNSS TDCP / RD cross-difference sequence. The specific formula is as follows:

[0086]

[0087] In the formula, The root mean square error (RMSE) before the seismic event is calculated from the cross-difference sequence of the high-frequency GNSS TDCP / RD from epoch k0 to epoch k1; rms() represents the calculation of the RMSE; k2 represents the P-wave arrival time of the aforementioned high-frequency GNSS / accelerometer co-located station. It is worth noting that in practical applications, the preferred values ​​for k0 and k1 are: k1 = k2, k0 = k1 - 300s, where s represents seconds.

[0088] Step 3: Based on the root mean square error before the earthquake event and the cross-difference sequence of the high-frequency GNSS TDCP / RD, the ground computing center adaptively selects the optimal high-frequency GNSS velocity from the average velocity of the high-frequency GNSS TDCP and the instantaneous velocity of the high-frequency GNSS RD by setting a prior threshold. The specific formula is as follows:

[0089]

[0090] In the formula, This represents the optimal high-frequency GNSS speed adaptively selected using prior threshold filtering; ratio is the prior ratio threshold, and extensive data testing shows that a ratio value of 6 is preferred. It should be noted that this is merely providing a feasible empirical value for ratio to facilitate a deeper and more thorough understanding of the various process details and features of this invention during implementation, and is not intended to limit the invention.

[0091] S143: The ground computing center reconstructs the broadband coseismic velocity waveform of the high-frequency GNSS / accelerometer co-located station in real time using dual-speed Kalman filtering based on the optimal high-frequency GNSS velocity and the zero-bias correction acceleration data.

[0092] Furthermore, in step S143, the ground computing center, based on the optimal high-frequency GNSS velocity and the zero-bias correction acceleration data, uses dual-speed Kalman filtering to reconstruct the broadband co-seismic velocity waveform of the high-frequency GNSS / accelerometer co-located station in real time. The specific dual-speed Kalman filtering fusion method is implemented as follows:

[0093] Let the state variable x of the dual-speed Kalman filter be... k For three-dimensional velocity v k ,Right now

[0094] x k =[v e,k v n,k v u,k ] T (5)

[0095] In the formula, x k v is the state variable for epoch k; e,k v n,k and v u,k These are the east-west, north-south, and vertical components of the velocity at epoch k, respectively; T This indicates the transpose operation.

[0096] The dual-speed Kalman filter system model is

[0097] x k =A k-1,k x k-1 +B k-1,k u k-1,k +w k-1,k ,w k-1,k ~N(0,Q) k-1,k (6)

[0098] A k-1.k =I 3×3 B k-1,k =τ a I 3×3

[0099] u k-1,k =(a k-1 +a k ) / 2 (7)

[0100] Q k-1,k =q a τa I 3×3

[0101] In the formula, A k-1,k Let x be the state transition matrix from epoch k-1 to epoch k. k-1 For k-1 epoch state variables; u k-1,k For system input; a k-1 and a k The zero-bias correction accelerations for epochs k-1 and k, respectively; B k-1,k Input matrix for the system; w k-1,k The system noise, w, follows a normal distribution. k-1,k ~N(0,Q) k-1,k ), Q k-1,k τ is the variance-covariance matrix of the system noise; a q represents the sampling interval of the accelerometer. a For the noise variance of the zero-bias correction acceleration, I 3×3 It is a 3×3 identity matrix.

[0102] The dual-speed Kalman filter observation model is

[0103]

[0104] In the formula, For high-frequency GNSS, the optimal speed is ε. k The optimal velocity noise for high-frequency GNSS follows a normal distribution ε k ~(0,R k ), R k H is the variance-covariance matrix of the velocity noise; k Design a matrix for the observation model; and These represent the noise variances of the east-west, north-south, and vertical components of the optimal velocity for high-frequency GNSS.

[0105] The dual-speed Kalman filter consists of two parts: time update and measurement update. The time update frequency should be consistent with the accelerometer sampling interval, while the measurement update frequency should be synchronized with the high-frequency GNSS sampling interval. The formulas for time update and measurement update are shown below:

[0106]

[0107]

[0108] Formula (9) is the time update formula, and formula (10) is the measurement update formula. In the formulas, () -1 P represents the matrix inversion operation; k For x kThe variance and covariance matrix; (-) and (+) represent the time update and measurement update results, respectively.

[0109] For the sake of simplicity, the method embodiments are described as a series of actions. However, those skilled in the art should understand that the embodiments of the present invention are not limited to the described order of actions, because according to the embodiments of the present invention, some steps can be performed in other orders or simultaneously. Furthermore, those skilled in the art should also understand that the embodiments described in the specification are preferred embodiments, and the actions involved are not necessarily essential to the embodiments of the present invention.

[0110] In specific implementation, the method proposed in the technical solution of this invention can be automatically executed by those skilled in the art using computer software technology. System devices for implementing the method, such as computer-readable storage media storing the corresponding computer program of the technical solution of this invention and computer equipment including the computer program running the corresponding computer program, should also be within the protection scope of this invention.

[0111] In some possible embodiments, a real-time reconstruction system for broadband coseismic velocity waveforms based on cascaded filtering is provided, including a processor and a memory. The memory is used to store program instructions, and the processor is used to call the stored instructions in the memory to execute the real-time reconstruction method for broadband coseismic velocity waveforms based on cascaded filtering as described above.

[0112] In some possible embodiments, a real-time reconstruction system for broadband coseismic velocity waveforms based on cascaded filtering is provided, including a readable storage medium storing a computer program. When the computer program is executed, it implements the real-time reconstruction method for broadband coseismic velocity waveforms based on cascaded filtering as described above.

[0113] The specific embodiments described herein are merely illustrative of the spirit of the invention. Those skilled in the art to which this invention pertains may make various modifications or additions to the described specific embodiments or use similar methods to substitute them, without departing from the spirit of the invention or exceeding the scope defined by the appended claims.

Claims

1. A method for real-time reconstruction of broadband co-oscillation velocity waveforms based on cascaded filtering data, characterized in that: Real-time reconstruction of broadband coseismic velocity waveforms is achieved by fusing cascaded filtering with raw high-frequency GNSS Doppler data, high-frequency GNSS carrier phase data, and accelerometer data. The process is as follows. The high-frequency GNSS / accelerometer co-located station collects high-frequency GNSS observation data, broadcast ephemeris and raw acceleration data in real time, and sends the collected high-frequency GNSS observation data, broadcast ephemeris and raw acceleration data to the ground computing center; The ground computing center calculates the average velocity of high-frequency GNSS TDCP and the instantaneous velocity of high-frequency GNSS RD in real time based on the high-frequency GNSS observation data received in real time and the broadcast ephemeris. The ground computing center performs pre-earthquake bias correction on the raw acceleration data received in real time to obtain bias-corrected acceleration data. The ground computing center reconstructs the broadband coseismic velocity waveform of the high-frequency GNSS / accelerometer co-located station in real time using cascaded filtering based on the high-frequency GNSS TDCP average velocity, the high-frequency GNSS RD instantaneous velocity, and the zero-bias correction acceleration data. This includes the following processing: The cascaded filtering includes a prior threshold filter and a dual-speed Kalman filter; The ground computing center adaptively selects the optimal high-frequency GNSS velocity using prior threshold filtering based on the average velocity of the high-frequency GNSS TDCP and the instantaneous velocity of the high-frequency GNSS RD. The ground computing center reconstructs the broadband coseismic velocity waveform of the high-frequency GNSS / accelerometer co-located station in real time using dual-speed Kalman filtering based on the optimal high-frequency GNSS velocity and the zero-bias correction acceleration data. The ground computing center adaptively selects the optimal high-frequency GNSS velocity based on the average velocity of the high-frequency GNSS TDCP and the instantaneous velocity of the high-frequency GNSS RD using a priori threshold filtering. The implementation method is as follows: The ground computing center calculates the high-frequency GNSS TDCP / RD cross-difference sequence based on the high-frequency GNSS TDCP average velocity sequence and the high-frequency GNSS RD instantaneous velocity sequence, as follows: In the formula, express High-frequency GNSS TDCP / RD cross-difference at epochs express High-frequency GNSS TDCP average velocity at each epoch; express Instantaneous velocity of high-frequency GNSS RD at an epoch; Indicates taking the absolute value; The ground computing center calculates the root mean square error of the high-frequency GNSS TDCP / RD cross-difference sequence before the earthquake event based on the high-frequency GNSS TDCP / RD cross-difference sequence as follows: In the formula, Indicates by Era to Root mean square error before seismic events calculated from high-frequency GNSS TDCP / RD cross-difference sequences at epochs; This indicates the calculation of the root mean square error. This indicates the arrival time of the P-wave at the aforementioned high-frequency GNSS / accelerometer co-located station; The ground computing center adaptively selects the optimal high-frequency GNSS velocity from the high-frequency GNSS TDCP average velocity and the high-frequency GNSS RD instantaneous velocity by setting a prior threshold based on the root mean square error before the earthquake event and the cross-difference sequence of the high-frequency GNSS TDCP / RD.

2. The method for real-time reconstruction of broadband co-oscillation velocity waveforms based on cascaded filtering according to claim 1, characterized in that: The ground computing center calculates the average velocity of high-frequency GNSS TDCP and the instantaneous velocity of high-frequency GNSS RD in real time based on the high-frequency GNSS observation data and the broadcast ephemeris received in real time, including the following processing. The high-frequency GNSS observation data includes high-frequency GNSS pseudorange observations, high-frequency GNSS carrier phase observations, and high-frequency GNSS raw Doppler observations; The ground computing center calculates the approximate coordinates of the high-frequency GNSS / accelerometer co-located station using standard single-point positioning based on the high-frequency GNSS pseudorange observations and the broadcast ephemeris. The ground computing center calculates the average velocity of the high-frequency GNSS TDCP of the high-frequency GNSS / accelerometer co-located station using the high-frequency GNSS carrier phase observation value, the broadcast ephemeris and the approximate coordinates, and the high-frequency GNSS carrier phase epoch difference technology. The ground computing center calculates the high-frequency GNSSRD instantaneous velocity of the high-frequency GNSS / accelerometer co-located station using the high-frequency GNSS raw Doppler observations, the broadcast ephemeris, and the approximate coordinates, employing the high-frequency GNSS raw Doppler instantaneous velocity method.

3. The method for real-time reconstruction of broadband co-oscillation velocity waveforms based on cascaded filtering according to claim 1, characterized in that: and The value can be , .

4. The method for real-time reconstruction of broadband co-oscillation velocity waveforms based on cascaded filtering according to claim 1, characterized in that: The adaptive selection of the optimal high-frequency GNSS velocity from the high-frequency GNSS TDCP average velocity and the high-frequency GNSS RD instantaneous velocity is achieved as follows: In the formula, This indicates the optimal high-frequency GNSS velocity adaptively selected using prior threshold filtering; This is the prior ratio threshold.

5. The method for real-time reconstruction of broadband co-oscillation velocity waveforms based on cascaded filtering according to claim 4, characterized in that: The value is 6.

6. The method for real-time reconstruction of broadband co-oscillation velocity waveforms based on cascaded filtering according to claim 1, 2, 3, 4, or 5, characterized in that: The ground computing center reconstructs the broadband coseismic velocity waveform of the high-frequency GNSS / accelerometer co-located station in real time using dual-speed Kalman filtering based on the optimal high-frequency GNSS velocity and the zero-bias correction acceleration data. The dual-speed Kalman filtering fusion method is as follows. Let the state variables of the dual-speed Kalman filter be... For three-dimensional velocity , In the formula, for epoch state variables; , and They are respectively The east-west, north-south, and vertical components of epochal velocity; Indicates the transpose operation; The dual-speed Kalman filter system model is In the formula, From Era to The state transition matrix of an epoch. for -1 epoch state quantity; For system input; and They are respectively and Zero-biased correction acceleration at epochs; Input matrix to the system; The system noise follows a normal distribution. , Let be the variance-covariance matrix of the system noise; The sampling interval of the accelerometer; The noise variance of the zero-bias correction acceleration. It is a 3×3 identity matrix; The dual-speed Kalman filter observation model is In the formula, The optimal speed for high-frequency GNSS; The optimal velocity noise for high-frequency GNSS follows a normal distribution. , Let be the variance-covariance matrix of the velocity noise; Design a matrix for the observation model; , and These are the noise variances of the east-west, north-south, and vertical components of the optimal velocity for high-frequency GNSS, respectively. The dual-speed Kalman filter consists of two parts: time update and measurement update. The time update frequency is consistent with the accelerometer sampling interval, while the measurement update frequency is synchronized with the high-frequency GNSS sampling interval. The time update formula is shown below. The measurement update formula is shown below. In the formula, This represents the matrix inversion operation; for The variance-covariance matrix; and These represent the results of time updates and measurement updates, respectively.

Images