CN112612045B - GNSS earthquake earth surface displacement monitoring method considering multipath and homomorphic errors - Google Patents

GNSS earthquake earth surface displacement monitoring method considering multipath and homomorphic errors Download PDF

Info

Publication number
CN112612045B
CN112612045B CN202011358334.8A CN202011358334A CN112612045B CN 112612045 B CN112612045 B CN 112612045B CN 202011358334 A CN202011358334 A CN 202011358334A CN 112612045 B CN112612045 B CN 112612045B
Authority
CN
China
Prior art keywords
gnss
station
seismic
period
displacement
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202011358334.8A
Other languages
Chinese (zh)
Other versions
CN112612045A (en
Inventor
郑凯
刘克中
张小红
马杰
王克浩
陈默子
曾旭明
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan University of Technology WUT
Original Assignee
Wuhan University of Technology WUT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Wuhan University of Technology WUT filed Critical Wuhan University of Technology WUT
Priority to CN202011358334.8A priority Critical patent/CN112612045B/en
Publication of CN112612045A publication Critical patent/CN112612045A/en
Application granted granted Critical
Publication of CN112612045B publication Critical patent/CN112612045B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B15/00—Measuring arrangements characterised by the use of electromagnetic waves or particle radiation, e.g. by the use of microwaves, X-rays, gamma rays or electrons
    • G01B15/06—Measuring arrangements characterised by the use of electromagnetic waves or particle radiation, e.g. by the use of microwaves, X-rays, gamma rays or electrons for measuring the deformation in a solid
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/24—Acquisition or tracking or demodulation of signals transmitted by the system
    • G01S19/25—Acquisition or tracking or demodulation of signals transmitted by the system involving aiding data received from a cooperating element, e.g. assisted GPS
    • G01S19/256—Acquisition or tracking or demodulation of signals transmitted by the system involving aiding data received from a cooperating element, e.g. assisted GPS relating to timing, e.g. time of week, code phase, timing offset

Abstract

The invention provides a GNSS earthquake earth surface displacement monitoring method considering multipath and homomorphic errors. The GNSS survey stations are processed station by PPP to obtain displacement sequences in the east direction, the north direction and the vertical direction, the start and stop time of seismic waves at each survey station is detected, and the observation time interval is divided into two non-seismic time intervals and seismic time intervals. And calculating the spatial response and the principal component of the displacement sequence of the measuring station in the non-earthquake period by using a principal component analysis method. And calculating homomorphic errors in the displacement sequence of the measuring station in the earthquake-affected period by combining the spatial response of the measuring station in the non-earthquake-affected period and the principal component of the measuring station in the earthquake-affected period. And calculating homomorphic errors in the displacement sequence of the measuring station in the non-seismic period by combining the spatial response and the principal component of the measuring station in the non-seismic period. And (3) deducting homomorphic errors from the displacement sequence of the measuring station, and eliminating multi-path errors by adopting star daily filtering to finally obtain the earth surface displacement sequence with high accuracy in the east direction, the north direction and the vertical direction. The multi-path error between the stations has no correlation, and simultaneously considers the characteristic of heterogeneous distribution of homomorphic errors.

Description

GNSS earthquake earth surface displacement monitoring method considering multipath and homomorphic errors
Technical Field
The invention belongs to the field of GNSS data processing, and relates to a GNSS seismic earth surface displacement monitoring method considering multipath and homomorphic errors.
Background
Surface displacements caused by major earthquakes can be sensitive probes of frictional conditions at fault interfaces and of the rheology of the near crust and upper mantle. The GNSS precise point positioning technology (PPP) is used as a new earth surface displacement acquisition means and can assist the earthquake table network to monitor the earthquake. At present, research mainly focuses on data processing of tens of seconds and minutes during earthquake, and high-frequency GNSS is used for rapidly extracting same-earthquake displacement and earthquake wave waveforms so as to meet application scenes with high time response requirements, such as earthquake early warning, rapid fracture inversion, rapid earthquake source parameter determination and the like. Although the GNSS can achieve short-time (several tens of seconds) dynamic displacement extraction accuracy on the order of millimeters under good observation conditions, three-dimensional dynamic displacement extraction accuracy is maintained around 10cm in many cases.
Sometimes accompanied by slow surface deformation before and after a major earthquake occurs. These deformations, with time spans on the order of hours (or less), help us to better understand the process of release of the large seismic stresses. The slow sliding phenomenon of the earthquake front earth surface can be used for providing effective constraint conditions for earthquake magnitude estimation, tsunami prediction and the like. Meanwhile, the phenomenon of aftersliding after early earthquake may be obvious within hours after the occurrence of the main earthquake. The earth surface creeping motion and aftershock between the main shock and the aftershock are easily submerged in the noise of the displacement sequence due to small magnitude, so that the misjudgment is generated on certain geophysical phenomena. In addition, the isoseismal displacement is used as an important observation quantity for seismic magnitude estimation and fault sliding distribution inversion, and is easily polluted by displacement generated by aftersliding, aftershock and the like within minutes or hours after the earthquake. These application scenarios all require sub-centimeter or even millimeter ground surface displacement monitoring accuracy.
Multipath errors and homomorphic errors (CME) in GNSS observations are considered the two major unmodeled errors in PPP. At present, the following problems mainly exist in the processing method aiming at the two errors:
the existing GNSS multi-path error processing method mainly aims at the application scene of short baseline relative positioning, the influence of most homomorphic errors can be eliminated or weakened by adopting an inter-station difference method, and the method requires that a reference station is not interfered by multi-paths. In addition, the method obtains the relative displacement of the monitoring station relative to the reference station, and the reference station can also be influenced by seismic waves to generate violent motion when a large earthquake occurs, so that the estimated earth surface same-seismic displacement is not accurate enough.
The study aiming at homomorphic errors is mainly based on a time sequence of a GNSS observation network long-time space solution, and is weakened by adopting a spatial filtering mode generally. In a traditional stack filtering (stacking) method, homomorphic errors are considered to be consistent to all stations in the whole measuring area at the same time, and can be obtained by stack filtering of the displacement of the stations in the measuring area.
Principal Component Analysis (PCA) can extract homomorphic errors simultaneously, homomorphic errors are allowed to have inconsistency in space, input data is required to be subjected to averaging and detrending, and extraction accuracy of homomorphic errors is affected if displacement of a survey station affected by seismic waves is taken into consideration.
Aiming at the problems, the invention provides the GNSS earthquake earth surface displacement monitoring method considering the multipath and homomorphic errors, which effectively improves the high-frequency GNSS dynamic displacement monitoring precision and provides a new finding for researching the process of induction and rupture of the earthquake.
Disclosure of Invention
The invention provides a GNSS earthquake earth surface displacement monitoring method considering multipath and homomorphic errors, which can effectively improve the high-frequency GNSS dynamic displacement monitoring precision.
The invention provides a GNSS seismic earth surface displacement monitoring method considering multipath and homomorphic errors, which comprises the following steps:
step 1, PPP processing is carried out on N GNSS survey stations in a survey area one by one to obtain an east displacement sequence, a north displacement sequence and a vertical displacement sequence of the N GNSS survey stations, average energy of seismic wave signals carried in the displacement sequences of the GNSS survey stations is calculated by constructing an average energy function, the arrival time of seismic waves of each GNSS survey station is further detected, and the ending time of the seismic waves of each GNSS survey station is detected;
step 2, calculating a seismic wave arrival time reference value of the GNSS survey station by combining the arrival time of the seismic wave of the GNSS survey station, calculating a seismic wave ending time reference value of the GNSS survey station by combining the ending time of the seismic wave of the GNSS survey station, and deducting the co-seismic deformation in the east, north and vertical displacement sequence of the GNSS survey station after earthquake elimination;
step 3, respectively calculating the space response and the principal component of the east, north and vertical displacement sequence of the GNSS survey station in the non-earthquake period by using a principal component analysis method;
step 4, combining the reference value of the arrival time of the seismic waves of the GNSS survey stations and the reference value of the ending time of the seismic waves of the GNSS survey stations to determine a plurality of GNSS survey stations which are not influenced by the seismic waves, and respectively calculating the main components of the east, north and vertical displacement sequences of the survey stations within the time period influenced by the seismic waves by using the displacements of the GNSS survey stations which are not influenced by the seismic waves;
step 5, combining the spatial response of the east, north and vertical displacement sequences of the GNSS measurement station in the non-earthquake period calculated in the step 3 and the principal components of the east, north and vertical displacement sequences of the GNSS measurement station in the earthquake wave influenced period in the step 4, calculating homomorphic errors of the GNSS measurement station in the east, north and vertical displacement sequences including the earthquake-influenced station in the earthquake wave influenced period; combining the spatial response and the principal component of the east, north and vertical displacement sequences of the GNSS observation station in the non-earthquake period calculated in the step 3, and calculating homomorphic errors of the GNSS observation station in the east, north and vertical displacement sequences in the non-earthquake period;
step 6, subtracting homomorphic errors in the non-seismic period and homomorphic errors of the GNSS measurement station in the time period influenced by seismic waves from the east, north and vertical displacement sequences of the GNSS measurement station, and then carrying out star daily filtering to eliminate multipath errors so as to finally obtain east, north and vertical high-precision GNSS earth surface displacement sequences;
preferably, the east displacement sequence, the north displacement sequence, and the vertical displacement sequence of the N GNSS stations in step 1 are:
Figure BDA0002803280350000031
k∈[1,N]
wherein, the first and the second end of the pipe are connected with each other,
Figure BDA0002803280350000032
representing an east displacement sequence of the kth GNSS rover,
Figure BDA0002803280350000033
representing the sequence of the north displacement of the kth GNSS station,
Figure BDA0002803280350000034
representing a sequence of vertical displacements of a kth GNSS observation station, N representing the number of GNSS observation stations;
step 1, calculating the energy of the seismic wave signal carried in the displacement sequence of the GNSS survey station by constructing an average energy function, wherein the energy of the seismic wave signal is as follows:
Figure BDA0002803280350000035
wherein the content of the first and second substances,
Figure BDA0002803280350000036
for the kth GNSS station at tepThe speed of the east direction of the moment,
Figure BDA0002803280350000037
for the kth GNSS station at tepThe speed of the north-going direction at the moment,
Figure BDA0002803280350000038
for the kth GNSS station at tepVertical velocity of time, Velk(tep) Representing the kth GNSS station at tepThe integrated speed of the three directions of the time, delta t is the sampling interval of the observed data, m is the window length, Fk(tep) Representing the kth GNSS station at tepAverage energy of seismic signals carried in time-shifted sequences, F0For background noise, k ∈[1,N]；
The kth station being in the ith directionepVelocity of time of day
Figure BDA0002803280350000039
The difference between the epochs of the displacement sequence can be used for obtaining:
Figure BDA0002803280350000041
wherein the content of the first and second substances,
Figure BDA0002803280350000042
represents the displacement of the kth station in the ith direction, delta t is the sampling interval of observed data, and k belongs to [1, N ]]Wherein, i-1 represents east direction, i-2 represents north direction, and i-3 represents vertical direction;
step 1, detecting the arrival time of the seismic waves of the GNSS survey station and the ending time of the seismic waves of the GNSS survey station are as follows:
wherein, when Fk(tep) Greater than 3 times F0When it is, consider tepThe time is the arrival time t of the seismic wave of the kth GNSS survey stationk,arrival；
When F is presentk(tep) Less than 3 times F0When it is, consider tepThe moment is the ending moment t of the seismic wave of the kth measuring stationk,stop；
Preferably, in the step 2, the seismic deformation displacement in the sequence of east, north and vertical displacement of the post-seismic GNSS survey station is deducted:
selecting the earliest arrival time of the seismic waves in all the GNSS survey stations as a seismic wave arrival time reference value:
tarrival＝min(tk,arrival)
selecting the latest seismic wave ending time in all GNSS survey stations as a seismic wave ending time reference value:
tstop＝max(tk,stop)
will observe the arrival time t1The time interval between the arrival time of the seismic wave and the reference value is used as the seismic observationSection (2):
Tpre＝[t1,tarrival]
defining the time interval between the seismic wave arrival time reference value and the seismic wave ending time reference value as the time interval affected by the seismic waves:
Tmove＝[tarrival,tstop]
the reference value of the seismic wave ending time is changed to the observation ending time tendThe time period in between is taken as the post-earthquake observation time period:
Tafter＝[tstop,tend]
defining the pre-earthquake observation period and the post-earthquake observation period as non-earthquake periods:
Tpa＝Tpre∪Tafter
and deducting the deformation of the same earthquake from the east, north and vertical displacement sequence of the post-earthquake GNSS survey station:
Figure BDA0002803280350000051
wherein the content of the first and second substances,
Figure BDA0002803280350000052
represents the displacement sequence of the k GNSS observation station after the homoseismal deformation is deducted in the i direction,
Figure BDA0002803280350000053
representing the original displacement sequence of the kth GNSS station in the ith direction,
Figure BDA0002803280350000054
representing the homoseismal deformation of the kth GNSS station in the ith direction, k belongs to [1, N ∈]Wherein, i-1 represents east direction, i-2 represents north direction, and i-3 represents vertical direction;
preferably, the step 3 of calculating the spatial response and the principal component of the east, north and vertical displacement sequence of the non-seismic GNSS positioning station by using the principal component analysis method respectively comprises:
the calculation of the displacement component matrix of the ith direction of the GNSS survey station in the non-seismic period can be represented as follows:
Figure BDA0002803280350000055
wherein the content of the first and second substances,
Figure BDA0002803280350000056
representing the raw displacement matrix of the i-th direction GNSS survey station before the earthquake,
Figure BDA0002803280350000057
a displacement matrix C representing the displacement of the GNSS survey station in the ith direction after earthquake after deducting the deformation of the same earthquakepreIs the total number of observation epochs, C, in the pre-earthquake observation periodafterIs the total number of observation epochs in the observation period after earthquake, Cpa＝Cpre+CafterThe total number of observation epochs in a non-seismic period is N, the number of GNSS stations is N, each column of the matrix represents a displacement vector of the ith direction of each station, i is 1 to represent the east direction, i is 2 to represent the north direction, and i is 3 to represent the vertical direction;
the covariance matrix formed by the ith direction displacement component can be expressed as
Figure BDA0002803280350000058
Which can be decomposed into:
Figure BDA0002803280350000059
wherein the content of the first and second substances,
Figure BDA00028032803500000510
as a characteristic value
Figure BDA00028032803500000511
The formed matrix has characteristic values arranged from big to small,
Figure BDA00028032803500000512
representing the correlation of characteristic values of the ith direction in the non-seismic periodThe eigenvectors are mutually orthogonal, and any N-dimensional matrix can be represented by N groups of orthogonal bases, wherein N belongs to [1, N ∈]N is the number of GNSS stations, i ═ 1 indicates the east direction, i ═ 2 indicates the north direction, and i ═ 3 indicates the vertical direction;
selecting
Figure BDA00028032803500000513
Represents the displacement component of the ith direction of the kth GNSS station at the t-th time again by the first n groups of bases:
Figure BDA00028032803500000514
wherein, the first and the second end of the pipe are connected with each other,
Figure BDA00028032803500000515
representing the displacement matrix formed by the ith direction at the t moment in the non-seismic period,
Figure BDA0002803280350000061
represents the value of the jth principal component in the ith direction at the t moment in the non-seismic period,
Figure BDA0002803280350000062
matrix representing ith direction eigenvector in non-seismic period
Figure BDA0002803280350000063
The kth value in the jth feature vector, te ∈ Tpa，j∈[1,n]，k∈[1,N]N is the number of GNSS stations, i is 1 indicating the east direction, i is 2 indicating the north direction, and i is 3 indicating the vertical direction;
Figure BDA0002803280350000064
can be expressed as:
Figure BDA0002803280350000065
wherein the content of the first and second substances,
Figure BDA0002803280350000066
representing the displacement of the ith direction of the kth measuring station at the time t in the non-seismic period,
Figure BDA0002803280350000067
represents the jth principal component in the ith direction in the non-seismic period,
Figure BDA0002803280350000068
also called the spatial response corresponding to the jth principal component in the ith direction in the non-seismic period, ttep∈Tpa，j∈[1,n]The i-1 indicates an east direction, the i-2 indicates a north direction, and the i-3 indicates a vertical direction;
the jth modality may be represented by the jth principal component and the jth spatial response;
the first n modes represent homomorphic errors in the region;
preferably, in step 4, the GNSS stations that are not affected by the seismic waves are:
considering the average energy F of the seismic signals carried in the GNSS survey station displacement sequence at any time, as described in connection with step 1kAre all less than 3 times F0The GNSS survey station is not influenced by seismic waves;
preferably, in step 4, the displacements of the GNSS stations not affected by the seismic waves are used to calculate the principal components of the sequences of the eastern, northern and vertical displacements of the stations in the time period affected by the seismic waves:
the method comprises the steps of calculating the jth main component in the ith direction of the GNSS measuring station in the time period influenced by the seismic waves in the tth direction by using the Ns displacement of the GNSS measuring station which is not influenced by the seismic waves in the time period influenced by the seismic wavesnomValue of time:
Figure BDA0002803280350000069
wherein the content of the first and second substances,
Figure BDA00028032803500000610
representing the ith party in the time period affected by seismic wavesTo the jth principal component at the tnomThe value of the time of day is,
Figure BDA00028032803500000611
represents the jth principal component in the ith direction in the time period influenced by the seismic waves,
Figure BDA00028032803500000612
the kth measuring station which is not influenced by the seismic waves in the period influenced by the seismic waves is positioned at the t-th directionnomThe displacement at a moment in time is,
Figure BDA00028032803500000613
the kth value, t, of the spatial response corresponding to the jth mode in the ith direction in the time period influenced by the seismic wavesnom∈Tmove，j∈[1,Ns]，k∈[1,Ns]Ns is the number of stations which are not affected by the seismic waves in the time period affected by the seismic waves, i is 1 to indicate east direction, i is 2 to indicate north direction, and i is 3 to indicate vertical direction;
and (3) defining the time interval between the seismic wave arrival time reference value and the seismic wave ending time reference value as the time interval influenced by the seismic waves in combination with the step 2:
Tmove＝[tarrival,tstop]
preferably, the spatial response of the sequence of the east, north and vertical displacements of the GNSS positioning station in the non-seismic period calculated in step 5 and the principal components of the sequence of the east, north and vertical displacements of the GNSS positioning station in the seismic wave affected period in step 4 are combined, and the homomorphic errors of the GNSS positioning station in the east, north and vertical displacements of the GNSS positioning station in the seismic wave affected period including the seismic positioning station are calculated:
firstly, determining the mode participating in homomorphic error calculation by judging the ratio of the normalized space response to the characteristic value:
the spatial response corresponding to the jth modality in the ith direction in the non-seismic period is as follows:
Figure BDA0002803280350000071
the normalization was performed as:
Figure BDA0002803280350000072
wherein the content of the first and second substances,
Figure BDA0002803280350000073
expressed as the s-th value in the spatial response corresponding to the j-th modality in the i-th direction during the non-seismic period,
Figure BDA0002803280350000074
represents the maximum value of all elements in the space response corresponding to the jth mode in the ith direction in the non-seismic period, s belongs to [1, N ]]，j∈[1,N]N is the number of GNSS stations, i ═ 1 indicates the east direction, i ═ 2 indicates the north direction, and i ═ 3 indicates the vertical direction;
calculating the statistic corresponding to the jth spatial response in the ith direction as follows:
Figure BDA0002803280350000075
wherein the content of the first and second substances,
Figure BDA0002803280350000076
representing the j normalized spatial response vector in the ith direction during the non-seismic period
Figure BDA0002803280350000077
The number of stations tested, j ∈ [1, N ], greater than the threshold thre _ sr]N is the number of GNSS stations, i ═ 1 indicates the east direction, i ═ 2 indicates the north direction, and i ═ 3 indicates the vertical direction;
calculating the eigenvalue related to the jth space response vector in the ith direction in the non-earthquake period
Figure BDA0002803280350000081
The corresponding statistics are:
Figure BDA0002803280350000082
and if the statistic satisfies the following conditions, the j-th mode in the ith direction is considered to be used for calculating the homomorphic error of the displacement in the ith direction:
Figure BDA0002803280350000083
wherein, thre1(thre1 ∈ [ 50%, 100% ]) represents the threshold corresponding to the spatial response statistic, and thre2 represents the threshold thre2 ∈ [ 1%, 100% ] corresponding to the eigenvalue statistic;
the jth modality may be represented by the jth principal component and the jth spatial response;
the first n modes represent homomorphic errors in the region;
screening n0 modes according to the judgment condition, and calculating homomorphic errors of all the stations including the earthquake-affected station in the seismic wave affected time period:
Figure BDA0002803280350000084
wherein the content of the first and second substances,
Figure BDA0002803280350000085
representing homomorphic errors in a displacement sequence of a kth measuring station including the measuring station in the ith direction in the time period influenced by the seismic waves,
Figure BDA0002803280350000086
is the principal component corresponding to the jth mode in the ith direction in the time period influenced by the seismic waves,
Figure BDA0002803280350000087
for the k-th value in the spatial response corresponding to the j-th modality in the i-th direction calculated in step 2, j ∈ [1, n0]，k∈[1,N]N is the number of GNSS stations, i ═ 1 indicates the east direction, i ═ 2 indicates the north direction, and i ═ 3 indicates the vertical direction;
and 5, combining the spatial response and the principal component of the east, north and vertical displacement sequence of the GNSS observation station in the non-earthquake period calculated in the step 3, calculating homomorphic errors of the east, north and vertical displacement sequence of the GNSS observation station in the non-earthquake period:
according to the screened n0 modes, calculating homomorphic errors of the GNSS observation station in the non-earthquake period:
Figure BDA0002803280350000088
wherein the content of the first and second substances,
Figure BDA0002803280350000091
representing homomorphic errors in the displacement sequence of the kth station in the ith direction during the non-seismic period,
Figure BDA0002803280350000092
for the principal component corresponding to the jth mode in the ith direction in the non-seismic period calculated in step 3,
Figure BDA0002803280350000093
the k value j epsilon [1, n0 in the spatial response corresponding to the j mode in the i direction in the non-seismic period calculated in the step 3]，k∈[1,N]N is the number of GNSS stations, i ═ 1 indicates the east direction, i ═ 2 indicates the north direction, and i ═ 3 indicates the vertical direction;
preferably, the homomorphic errors of the sequences of east, north and vertical displacements of the GNSS stations including the seismometer station in the time period affected by the seismic waves, and the homomorphic errors of the sequences of east, north and vertical displacements of the GNSS stations in the non-seismal period, which are obtained in the step 5, are subtracted from the displacements of the N stations in the step 6:
Figure BDA0002803280350000094
Figure BDA0002803280350000095
wherein the content of the first and second substances,
Figure BDA0002803280350000096
representing the homomorphic error in the displacement sequence of the kth station in the ith direction in all observation times,
Figure BDA0002803280350000097
representing homomorphic errors in a displacement sequence of a kth measuring station including the measuring station in the ith direction in the time period influenced by the seismic waves,
Figure BDA0002803280350000098
representing the homomorphic error in the displacement sequence of the kth measuring station in the ith direction in the non-seismic period,
Figure BDA0002803280350000099
representing the original displacement sequence of the kth station in the ith direction, k is the [1, N ]]N is the number of GNSS stations, i ═ 1 indicates the east direction, i ═ 2 indicates the north direction, and i ═ 3 indicates the vertical direction;
step 6, the displacement sequence after homomorphic error is deducted
Figure BDA00028032803500000910
And (3) carrying out sidereal day filtering to finally obtain a high-precision position sequence:
Figure BDA00028032803500000911
wherein, the first and the second end of the pipe are connected with each other,
Figure BDA00028032803500000912
the multipath error of the kth GNSS station calculated by the sun-sun filtering in the ith direction is represented as:
Figure BDA00028032803500000913
wherein the content of the first and second substances,
Figure BDA0002803280350000101
represents the displacement sequence of the k testing station at d days after the homomorphic error is deducted in the i direction, k belongs to [1, N ∈]，d∈[1,D]N is the number of GNSS stations, D is the number of days involved in calculating multipath errors, i-1 indicates the east direction, i-2 indicates the north direction, and i-3 indicates the vertical direction.
Compared with the prior art, the invention has the following advantages and beneficial effects:
by adopting the PPP data processing method, the GNSS observation stations are not influenced by the displacement deviation of the reference station, the absolute coordinates relative to a global reference frame can be obtained, and meanwhile, the multipath errors among the observation stations have no correlation.
A principal component analysis method is introduced to separate homomorphic errors, so that the defect that the homomorphic errors are required to be uniformly distributed in the whole measuring area in the traditional stack filtering method is overcome, and the regional homomorphic errors with non-uniform characteristics in the region can be extracted.
The method for extracting the homomorphic errors of the seismological measurement station in the earthquake period is provided by combining the principal component analysis method, and the defect that the homomorphic errors are extracted to have deviation due to the influence of homomorphic displacement on the measurement station is avoided.
Drawings
FIG. 1: is a flow chart of the method of the present invention.
FIG. 2: the homomorphic error schematic diagram is solved by the principal component analysis method in the embodiment of the invention.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
The following describes an embodiment of the present invention with reference to fig. 1, and a GNSS seismic surface displacement monitoring method considering multipath and homomorphic errors includes the following steps:
step 1, PPP processing is carried out on N GNSS measurement stations in a measurement area one by one to obtain an east displacement sequence, a north displacement sequence and a vertical displacement sequence of the N GNSS measurement stations, average energy of seismic wave signals carried in the displacement sequences of the GNSS measurement stations is calculated by constructing an average energy function, the arrival time of seismic waves of each GNSS measurement station is further detected, and the ending time of the seismic waves of each GNSS measurement station is detected;
in the step 1, the east displacement sequence, the north displacement sequence and the vertical displacement sequence of the 200 GNSS stations are as follows:
Figure BDA0002803280350000102
k∈[1,N]
wherein the content of the first and second substances,
Figure BDA0002803280350000111
representing the east displacement sequence of the kth GNSS station,
Figure BDA0002803280350000112
representing the sequence of the north displacement of the kth GNSS station,
Figure BDA0002803280350000113
indicating a vertical displacement sequence of a kth GNSS station, wherein N is 200 and indicates the number of the GNSS stations;
the step 1 of calculating the energy of the seismic wave signal carried in the displacement sequence of the GNSS survey station by constructing an average energy function is as follows:
Figure BDA0002803280350000114
wherein the content of the first and second substances,
Figure BDA0002803280350000115
for the kth GNSS station at tepThe speed of the east direction of the time,
Figure BDA0002803280350000116
for the kth GNSS station at tepThe speed of the north-going direction at the moment,
Figure BDA0002803280350000117
for the kth GNSS station at tepVertical velocity of time, Velk(tep) Representing the kth GNSS station at tepThe integrated velocity in three directions at time, Δ t ═ 1s, m ═ 10, the window length, and F, is the sampling interval of the observed datak(tep) Representing the kth GNSS station at tepAverage energy of seismic signals carried in time-shifted sequences, F02cm as background noise, k ∈ [1,200 ]]；
The kth station being in the ith directionepVelocity of time of day
Figure BDA0002803280350000118
The difference between the epochs of the displacement sequence can be used for obtaining:
Figure BDA0002803280350000119
wherein the content of the first and second substances,
Figure BDA00028032803500001110
represents the displacement of the kth station in the ith direction, and Δ t ═ 1s is the sampling interval of observed data, k ∈ [1,200 ∈]Wherein, i-1 represents east direction, i-2 represents north direction, and i-3 represents vertical direction;
step 1, detecting the arrival time of the seismic waves of the GNSS survey station and the ending time of the seismic waves of the GNSS survey station are as follows:
wherein, when Fk(tep) Greater than 3 times F0When t is 2cm, t is consideredepThe time is the arrival time t of the seismic wave of the kth GNSS survey stationk,arrival；
When F is presentk(tep) Less than 3 times F0When t is 2cm, t is consideredepThe moment is the ending moment t of the seismic wave of the kth measuring stationk,stop；
Step 2, calculating a seismic wave arrival time reference value of the GNSS survey station by combining the arrival time of the seismic wave of the GNSS survey station, calculating a seismic wave ending time reference value of the GNSS survey station by combining the ending time of the seismic wave of the GNSS survey station, and deducting the co-seismic deformation in the east, north and vertical displacement sequence of the GNSS survey station after earthquake elimination;
deducting the same-seismic deformation displacement in the sequence of east, north and vertical displacements of the post-seismic GNSS survey station in the step 2:
selecting the earliest arrival time of the seismic waves in all the GNSS survey stations as a seismic wave arrival time reference value:
tarrival＝min(tk,arrival)
selecting the latest seismic wave ending time in all GNSS survey stations as a seismic wave ending time reference value:
tstop＝max(tk,stop)
will observe the arrival time t1The time interval between the arrival time reference value of the seismic wave is used as a pre-earthquake observation time interval:
Tpre＝[t1,tarrival]
defining the time interval between the seismic wave arrival time reference value and the seismic wave ending time reference value as the time interval influenced by the seismic waves:
Tmove＝[tarrival,tstop]
the reference value of the seismic wave ending time is changed to the observation ending time tendThe time period in between is taken as the post-earthquake observation time period:
Tafter＝[tstop,tend]
defining the pre-earthquake observation time interval and the post-earthquake observation time interval as non-earthquake time periods:
Tpa＝Tpre∪Tafter
deducting the homoseism deformation from the sequence of east, north and vertical displacement of the post-earthquake GNSS survey station:
Figure BDA0002803280350000121
wherein the content of the first and second substances,
Figure BDA0002803280350000122
represents the displacement sequence of the k GNSS survey station after deducting the homoseismal deformation in the ith direction,
Figure BDA0002803280350000123
representing the original displacement sequence of the kth GNSS station in the ith direction,
Figure BDA0002803280350000124
represents the homoseismal deformation of the kth GNSS survey station in the ith direction, and k belongs to [1,200 ]]Wherein, i-1 represents east direction, i-2 represents north direction, and i-3 represents vertical direction;
step 3, respectively calculating the space response and the principal component of the east, north and vertical displacement sequence of the GNSS survey station in the non-earthquake period by using a principal component analysis method;
referring to fig. 2, the step 3 of calculating the spatial response and the principal component of the east, north and vertical displacement sequence of the GNSS survey station in the non-seismic period by using the principal component analysis method respectively includes:
the calculation of the displacement component matrix of the ith direction of the GNSS survey station in the non-seismic period can be represented as follows:
Figure BDA0002803280350000131
wherein the content of the first and second substances,
Figure BDA0002803280350000132
representing the raw displacement matrix of the i-th direction GNSS survey station before the earthquake,
Figure BDA0002803280350000133
representing a displacement matrix of the GNSS survey station in the ith direction after the earthquake after the same-earthquake deformation is deducted, Cpre1000 is the total number of observation epochs in the observation period before earthquake, Cafter1000 is the total number of observation epochs in the post-earthquake observation period, Cpa＝Cpre+Cafter2000 is the total number of observation epochs in the non-seismic period, N200 is the number of GNSS stations, and each column of the matrix represents the station numberDisplacement vectors in i directions, i-1 represents an east direction, i-2 represents a north direction, and i-3 represents a vertical direction;
the covariance matrix formed by the ith direction displacement component can be expressed as
Figure BDA0002803280350000134
Which can be decomposed into:
Figure BDA0002803280350000135
wherein the content of the first and second substances,
Figure BDA0002803280350000136
as a characteristic value
Figure BDA0002803280350000137
The formed matrix has characteristic values arranged from big to small,
Figure BDA0002803280350000138
representing the matrix formed by the eigenvectors corresponding to the characteristic values in the ith direction in the non-seismic period, wherein the eigenvectors are mutually orthogonal, and any n-dimensional matrix can be represented by n groups of orthogonal bases, and n belongs to [1,200 ]]Wherein, i-1 represents east direction, i-2 represents north direction, and i-3 represents vertical direction;
selecting
Figure BDA0002803280350000139
The first n in (2) set of bases represents the displacement component of the kth GNSS station in the ith direction at the tth time again:
Figure BDA00028032803500001310
wherein the content of the first and second substances,
Figure BDA00028032803500001311
representing the displacement matrix formed by the ith direction at the t moment in the non-seismic period,
Figure BDA00028032803500001312
represents the value of the jth principal component in the ith direction at the t moment in the non-seismic period,
Figure BDA00028032803500001313
matrix representing ith direction eigenvector in non-seismic period
Figure BDA00028032803500001314
The kth value in the jth feature vector, te ∈ Tpa，j∈[1,2]，k∈[1,200]Wherein, i-1 represents east direction, i-2 represents north direction, and i-3 represents vertical direction;
Figure BDA0002803280350000141
can be expressed as:
Figure BDA0002803280350000142
wherein the content of the first and second substances,
Figure BDA0002803280350000143
representing the displacement of the ith direction of the kth measuring station at the time t in the non-seismic period,
Figure BDA0002803280350000144
represents the jth principal component in the ith direction in the non-seismic period,
Figure BDA0002803280350000145
also called the spatial response corresponding to the jth principal component in the ith direction in the non-seismic period, ttep∈Tpa，j∈[1,2]Wherein, i-1 represents east direction, i-2 represents north direction, and i-3 represents vertical direction;
the jth modality may be represented by the jth principal component and the jth spatial response;
the first n modes represent homomorphic errors in the region;
step 4, combining the reference value of the arrival time of the seismic waves of the GNSS measurement stations and the reference value of the ending time of the seismic waves of the GNSS measurement stations to determine a plurality of GNSS measurement stations which are not influenced by the seismic waves, and respectively calculating the principal components of the east, north and vertical displacement sequences of the measurement stations in the time period influenced by the seismic waves by using the displacements of the GNSS measurement stations which are not influenced by the seismic waves;
step 4, the GNSS survey station which is not influenced by the seismic waves is as follows:
considering the average energy F of the seismic signals carried in the GNSS survey station displacement sequence at any time, as described in connection with step 1kAre all less than 3 times F0When the distance is 2cm, the GNSS survey station is not influenced by seismic waves;
preferably, in step 4, the main components of the east, north and vertical displacement sequences of the measuring station in the time period affected by the seismic waves are respectively calculated by using the displacements of the GNSS measuring station not affected by the seismic waves:
using the displacement of Ns-50 GNSS measuring stations which are not influenced by the seismic waves in the time interval influenced by the seismic waves to calculate the jth principal component in the ith direction of the GNSS measuring station in the time interval influenced by the seismic waves in the tth directionnomThe value of the time:
Figure BDA0002803280350000146
wherein the content of the first and second substances,
Figure BDA0002803280350000147
representing the jth principal component in the ith direction in the t-th time interval affected by the seismic wavesnomThe value of the time of day is,
Figure BDA0002803280350000148
represents the jth principal component in the ith direction in the time period influenced by the seismic waves,
Figure BDA0002803280350000151
the kth measuring station which is not influenced by the seismic waves in the ith direction in the period of being influenced by the seismic waves is subjected to the tthnomThe displacement at a moment in time is,
Figure BDA0002803280350000152
when influenced by seismic wavesThe kth value, t, in the spatial response corresponding to the jth modality in the ith direction within the segmentnom∈Tmove，j∈[1,50]，k∈[1,50]Wherein, i-1 represents east direction, i-2 represents north direction, and i-3 represents vertical direction;
and (3) defining the time interval between the seismic wave arrival time reference value and the seismic wave ending time reference value as the time interval influenced by the seismic waves in combination with the step 2:
Tmove＝[tarrival,tstop]
step 5, combining the spatial response of the east, north and vertical displacement sequences of the GNSS measuring station in the non-earthquake period calculated in the step 3 and the main components of the east, north and vertical displacement sequences of the GNSS measuring station in the earthquake wave influenced period in the step 4, calculating homomorphic errors of the GNSS measuring station in the east, north and vertical displacement sequences including the earthquake-influenced station in the earthquake wave influenced period; combining the spatial response and the principal component of the east, north and vertical displacement sequences of the GNSS observation station in the non-earthquake period calculated in the step 3, and calculating homomorphic errors of the GNSS observation station in the east, north and vertical displacement sequences in the non-earthquake period;
combining the spatial response of the east, north and vertical displacement sequences of the GNSS measurement station in the non-seismic period calculated in the step 3 and the main components of the east, north and vertical displacement sequences of the GNSS measurement station in the seismic wave influenced period in the step 4, calculating homomorphic errors of the GNSS measurement station in the east, north and vertical displacement sequences including the seismic measurement station in the seismic wave influenced period:
firstly, determining the mode participating in homomorphic error calculation by judging the ratio of the normalized space response to the characteristic value:
the spatial response corresponding to the jth mode in the ith direction in the non-seismic period is as follows:
Figure BDA0002803280350000153
the normalization was performed as:
Figure BDA0002803280350000154
wherein the content of the first and second substances,
Figure BDA0002803280350000155
expressed as the s-th value in the spatial response corresponding to the j-th mode in the ith direction during the non-seismic period,
Figure BDA0002803280350000156
represents the maximum value of all elements in the spatial response corresponding to the jth mode in the ith direction in the non-seismic period, s is equal to [1,200 ]]，j∈[1,200]The i-1 indicates an east direction, the i-2 indicates a north direction, and the i-3 indicates a vertical direction;
calculating statistics corresponding to the jth spatial response in the ith direction as follows:
Figure BDA0002803280350000161
wherein, the first and the second end of the pipe are connected with each other,
Figure BDA0002803280350000162
representing the j normalized spatial response vector in the ith direction during the non-seismic period
Figure BDA0002803280350000163
More than 25% station number, j belongs to [1, N [ ]]N-200 is the GNSS station number, i-1 indicates the east direction, i-2 indicates the north direction, and i-3 indicates the vertical direction;
calculating the eigenvalue related to the jth space response vector in the ith direction in the non-earthquake period
Figure BDA0002803280350000164
The corresponding statistics are:
Figure BDA0002803280350000165
and if the statistic satisfies the following condition, the j mode in the ith direction is considered to be used for calculating the homomorphic error of the displacement in the ith direction:
Figure BDA0002803280350000166
wherein, thre 1-50% represents the threshold corresponding to the spatial response statistic, and thre 2-1% represents the threshold corresponding to the eigenvalue statistic;
the jth modality may be represented by the jth principal component and the jth spatial response;
the first n modes represent homomorphic errors in the region;
screening n 0-2 modes according to the judgment condition, and calculating homomorphic errors of all stations including the earthquake-affected station in the time period affected by the earthquake waves:
Figure BDA0002803280350000167
wherein the content of the first and second substances,
Figure BDA0002803280350000168
representing homomorphic errors in a displacement sequence of a kth measuring station including the measuring station in the ith direction in the time period influenced by the seismic waves,
Figure BDA0002803280350000169
is the principal component corresponding to the jth mode in the ith direction in the time period influenced by the seismic waves,
Figure BDA00028032803500001610
for the k-th value in the spatial response corresponding to the j-th mode in the i-th direction calculated in step 2, j ∈ [1,2 ]]，k∈[1,200]Wherein, i-1 represents east direction, i-2 represents north direction, and i-3 represents vertical direction;
and 5, combining the spatial response and the principal component of the east, north and vertical displacement sequence of the GNSS observation station in the non-earthquake period calculated in the step 3, calculating homomorphic errors of the east, north and vertical displacement sequence of the GNSS observation station in the non-earthquake period:
according to the screened n 0-2 modes, calculating homomorphic errors of the GNSS survey station in the non-earthquake period:
Figure BDA0002803280350000171
wherein the content of the first and second substances,
Figure BDA0002803280350000172
representing homomorphic errors in the displacement sequence of the kth station in the ith direction during the non-seismic period,
Figure BDA0002803280350000173
for the principal component corresponding to the jth mode in the ith direction in the non-seismic period calculated in step 3,
Figure BDA0002803280350000174
the k value j epsilon [1,2 ] in the space response corresponding to the j mode in the i direction in the non-seismic period calculated by the step 3]，k∈[1,200]Wherein, i-1 represents east direction, i-2 represents north direction, and i-3 represents vertical direction;
step 6, subtracting homomorphic errors in the non-earthquake period and homomorphic errors of the GNSS survey station in the time period influenced by earthquake waves from the east, north and vertical displacement sequences of the GNSS survey station, and then carrying out sidereal day filtering to eliminate multi-path errors so as to finally obtain east, north and vertical high-precision GNSS earth surface displacement sequences;
and 6, deducting homomorphic errors of east, north and vertical displacement sequences of the GNSS measuring station including the earthquake-affected station in the time period affected by the earthquake waves obtained in the step 5 from the displacements of the N measuring stations, and homomorphic errors of the east, north and vertical displacement sequences of the GNSS measuring station in the non-earthquake period:
Figure BDA0002803280350000175
Figure BDA0002803280350000176
wherein the content of the first and second substances,
Figure BDA0002803280350000177
representing homomorphic errors in the displacement sequence of the kth station in the ith direction for all observation times,
Figure BDA0002803280350000178
representing homomorphic errors in a displacement sequence of a kth measuring station including the measuring station in the ith direction in the time period influenced by the seismic waves,
Figure BDA0002803280350000179
representing homomorphic errors in the displacement sequence of the kth station in the ith direction during the non-seismic period,
Figure BDA00028032803500001710
representing the original displacement sequence of the kth station in the ith direction, k ∈ [1, N]N-200 is the GNSS station number, i-1 indicates the east direction, i-2 indicates the north direction, and i-3 indicates the vertical direction;
step 6, the displacement sequence after homomorphic error is deducted
Figure BDA0002803280350000181
And (3) carrying out sidereal day filtering to finally obtain a high-precision position sequence:
Figure BDA0002803280350000182
wherein the content of the first and second substances,
Figure BDA0002803280350000183
the multipath error of the kth GNSS station calculated by the sidereal day filtering in the ith direction is represented as:
Figure BDA0002803280350000184
wherein the content of the first and second substances,
Figure BDA0002803280350000185
represents the displacement sequence of the k station at d day after eliminating homomorphic error in the i direction, k is the [1,200 ]]，d∈[1,D]D-7 indicates the number of days involved in calculating the multipath error, i-1 indicates the east direction, i-2 indicates the north direction, and i-3 indicates the vertical direction.
It should be understood that the above description of the preferred embodiments is given for clarity and not for any purpose of limitation, and that various changes, substitutions and alterations can be made herein without departing from the spirit and scope of the invention as defined by the appended claims.

Claims (3)

1. A GNSS seismic earth surface displacement monitoring method considering multipath and homomorphic errors is characterized by comprising the following steps:
step 1, PPP processing is carried out on N GNSS measurement stations in a measurement area one by one to obtain an east displacement sequence, a north displacement sequence and a vertical displacement sequence of the N GNSS measurement stations, average energy of seismic wave signals carried in the displacement sequences of the GNSS measurement stations is calculated by constructing an average energy function, the arrival time of seismic waves of each GNSS measurement station is further detected, and the ending time of the seismic waves of each GNSS measurement station is detected;
step 2, calculating a seismic wave arrival time reference value of the GNSS survey station by combining the arrival time of the seismic wave of the GNSS survey station, calculating a seismic wave ending time reference value of the GNSS survey station by combining the ending time of the seismic wave of the GNSS survey station, and deducting the co-seismic deformation in the east, north and vertical displacement sequence of the GNSS survey station after earthquake elimination;
step 3, respectively calculating the space response and the principal component of the east, north and vertical displacement sequence of the GNSS survey station in the non-earthquake period by using a principal component analysis method;
step 4, combining the reference value of the arrival time of the seismic waves of the GNSS measurement stations and the reference value of the ending time of the seismic waves of the GNSS measurement stations to determine a plurality of GNSS measurement stations which are not influenced by the seismic waves, and respectively calculating the principal components of the east, north and vertical displacement sequences of the measurement stations in the time period influenced by the seismic waves by using the displacements of the GNSS measurement stations which are not influenced by the seismic waves;
step 5, combining the spatial response of the east, north and vertical displacement sequences of the GNSS measurement station in the non-earthquake period calculated in the step 3 and the principal components of the east, north and vertical displacement sequences of the GNSS measurement station in the earthquake wave influenced period in the step 4, calculating homomorphic errors of the GNSS measurement station in the east, north and vertical displacement sequences including the earthquake-influenced station in the earthquake wave influenced period; combining the spatial response and the principal component of the east, north and vertical displacement sequences of the GNSS observation station in the non-earthquake period calculated in the step 3, and calculating homomorphic errors of the GNSS observation station in the east, north and vertical displacement sequences in the non-earthquake period;
step 6, subtracting homomorphic errors in the non-earthquake period and homomorphic errors of the GNSS survey station in the time period influenced by earthquake waves from the east, north and vertical displacement sequences of the GNSS survey station, and then carrying out sidereal day filtering to eliminate multi-path errors so as to finally obtain east, north and vertical high-precision GNSS earth surface displacement sequences;
in the step 3, the spatial response and the principal component of the east, north and vertical displacement sequence of the GNSS survey station in the non-seismic period are respectively calculated by using a principal component analysis method as follows:
the calculation of the displacement component matrix of the ith direction of the GNSS survey station in the non-seismic period can be represented as follows:
Figure FDA0003673517160000011
wherein the content of the first and second substances,
Figure FDA0003673517160000012
a raw displacement matrix representing the i direction GNSS stations before the earthquake,
Figure FDA0003673517160000013
a displacement matrix C representing the displacement of the GNSS survey station in the ith direction after earthquake after deducting the deformation of the same earthquakepreIs the total number of observation epochs, C, in the pre-earthquake observation periodafterIs the total number of observation epochs in the observation period after earthquake, Cpa＝Cpre+CafterThe total number of observation epochs in a non-seismic period is N, the number of GNSS stations is N, each column of the matrix represents a displacement vector of the ith direction of each station, i is 1 to represent the east direction, i is 2 to represent the north direction, and i is 3 to represent the vertical direction;
the covariance matrix formed by the ith direction displacement component can be expressed as
Figure FDA0003673517160000021
It can be decomposed into:
Figure FDA0003673517160000022
wherein, the first and the second end of the pipe are connected with each other,
Figure FDA0003673517160000023
is a characteristic value
Figure FDA0003673517160000024
The formed matrix has characteristic values arranged from big to small,
Figure FDA0003673517160000025
representing a matrix formed by eigenvectors corresponding to characteristic values in the ith direction in a non-seismic period, wherein the eigenvectors are orthogonal to each other, and any N-dimensional matrix can be represented by N groups of orthogonal bases, and N belongs to [1, N ]]N is the number of GNSS stations, i ═ 1 indicates the east direction, i ═ 2 indicates the north direction, and i ═ 3 indicates the vertical direction;
selecting
Figure FDA0003673517160000026
Represents the displacement component of the ith direction of the kth GNSS station at the t-th time again by the first n groups of bases:
Figure FDA0003673517160000027
wherein the content of the first and second substances,
Figure FDA0003673517160000028
representing the displacement matrix formed by the ith direction at the t time in the non-seismic period,
Figure FDA0003673517160000029
represents the value of the jth principal component in the ith direction at the tth moment in the non-seismic period,
Figure FDA00036735171600000210
matrix representing ith direction eigenvector in non-seismic period
Figure FDA00036735171600000211
The kth value in the jth feature vector, te ∈ Tpa，j∈[1,n]，k∈[1,N]The i-1 indicates an east direction, the i-2 indicates a north direction, and the i-3 indicates a vertical direction;
Figure FDA00036735171600000212
can be expressed as:
Figure FDA00036735171600000213
wherein the content of the first and second substances,
Figure FDA00036735171600000214
represents the displacement of the ith direction of the kth measuring station in the non-seismic period at the time t,
Figure FDA00036735171600000215
represents the jth principal component in the ith direction in the non-seismic period,
Figure FDA00036735171600000216
also known as non-shockSpatial response t corresponding to jth principal component in ith direction in the periodtep∈Tpa，j∈[1,n]Wherein, i-1 represents east direction, i-2 represents north direction, and i-3 represents vertical direction;
the jth modality may be represented by the jth principal component and the jth spatial response;
the first n modes represent homomorphic errors in the region;
step 4, the GNSS survey station which is not influenced by the seismic waves is as follows:
considering the average energy F of the seismic signals carried in the GNSS survey station displacement sequence at any time, as described in connection with step 1kAre all less than 3 times F0The GNSS survey station is not influenced by seismic waves;
and 4, respectively calculating main components of east, north and vertical displacement sequences of the measuring station in the time period influenced by the seismic waves by using the displacement of the GNSS measuring station which is not influenced by the seismic waves:
calculating the jth principal component in the ith direction of the GNSS measuring station in the time period influenced by the seismic waves at the tth position by using the Ns displacement of the GNSS measuring station which is not influenced by the seismic waves in the time period influenced by the seismic wavesnomThe value of the time:
Figure FDA0003673517160000031
wherein the content of the first and second substances,
Figure FDA0003673517160000032
representing the j main component in the ith direction in the t time period influenced by the seismic wavesnomThe value of the time of day is,
Figure FDA0003673517160000033
represents the jth principal component of the ith direction in the time period influenced by the seismic waves,
Figure FDA0003673517160000034
the kth measuring station which is not influenced by the seismic waves in the period influenced by the seismic waves is positioned at the t-th directionnomDisplacement of time of day，
Figure FDA0003673517160000035
The kth value, t, of the spatial response corresponding to the jth mode in the ith direction in the time period influenced by the seismic wavesnom∈Tmove，j∈[1,Ns]，k∈[1,Ns]Ns is the number of stations which are not affected by the seismic waves in the time period affected by the seismic waves, i is 1 to indicate east, i is 2 to indicate north, and i is 3 to indicate vertical;
and (3) defining the time interval between the seismic wave arrival time reference value and the seismic wave ending time reference value as the time interval influenced by the seismic waves in combination with the step 2:
Tmove＝[tarrival,tstop]
combining the spatial response of the east, north and vertical displacement sequences of the GNSS measurement station in the non-seismic period calculated in the step 3 and the main components of the east, north and vertical displacement sequences of the GNSS measurement station in the seismic wave influenced period in the step 4, calculating homomorphic errors of the GNSS measurement station in the east, north and vertical displacement sequences including the seismic measurement station in the seismic wave influenced period:
firstly, determining the mode participating in homomorphic error calculation by judging the ratio of the normalized space response to the characteristic value:
the spatial response corresponding to the jth mode in the ith direction in the non-seismic period is as follows:
Figure FDA0003673517160000036
the normalization was performed as:
Figure FDA0003673517160000037
wherein the content of the first and second substances,
Figure FDA0003673517160000041
expressed as corresponding to the jth mode in the ith direction in the non-seismic periodThe s-th value in the spatial response,
Figure FDA0003673517160000042
represents the maximum value of all elements in the spatial response corresponding to the jth mode in the ith direction in the non-seismic period, and belongs to [1, N ∈]，j∈[1,N]N is the number of GNSS stations, i ═ 1 indicates the east direction, i ═ 2 indicates the north direction, and i ═ 3 indicates the vertical direction;
calculating statistics corresponding to the jth spatial response in the ith direction as follows:
Figure FDA0003673517160000043
wherein the content of the first and second substances,
Figure FDA0003673517160000044
representing the j normalized spatial response vector in the ith direction during the non-seismic period
Figure FDA0003673517160000045
The number of stations tested, j ∈ [1, N ], greater than the threshold thre _ sr]N is the number of GNSS stations, i is 1 indicating the east direction, i is 2 indicating the north direction, and i is 3 indicating the vertical direction;
computing eigenvalues related to the jth spatial response vector in the ith direction during the non-seismic period
Figure FDA0003673517160000046
The corresponding statistics are:
Figure FDA0003673517160000047
and if the statistic satisfies the following condition, the j mode in the ith direction is considered to be used for calculating the homomorphic error of the displacement in the ith direction:
Figure FDA0003673517160000048
wherein, thre1 represents the threshold corresponding to the spatial response statistic, thre1 belongs to [ 50%, 100% ], and thre2 represents the threshold thre2 corresponding to the eigenvalue statistic belongs to [ 1%, 100% ];
the jth modality may be represented by a jth principal component and a jth spatial response;
the first n modes represent homomorphic errors in the region;
screening n0 modes according to the judgment condition, and calculating homomorphic errors of all the stations including the earthquake-affected station in the seismic wave affected time period:
Figure FDA0003673517160000049
wherein the content of the first and second substances,
Figure FDA00036735171600000410
representing homomorphic errors in a displacement sequence of a kth measuring station including the measuring station in the ith direction in the time period influenced by the seismic waves,
Figure FDA00036735171600000411
is the principal component corresponding to the jth mode in the ith direction in the time period influenced by the seismic waves,
Figure FDA00036735171600000412
for the k-th value in the spatial response corresponding to the j-th mode in the i-th direction calculated in step 2, j ∈ [1, n0]，k∈[1,N]N is the number of GNSS stations, i ═ 1 indicates the east direction, i ═ 2 indicates the north direction, and i ═ 3 indicates the vertical direction;
and 5, combining the spatial response and the principal component of the east, north and vertical displacement sequence of the GNSS observation station in the non-earthquake period calculated in the step 3, calculating homomorphic errors of the east, north and vertical displacement sequence of the GNSS observation station in the non-earthquake period:
according to the screened n0 modes, calculating homomorphic errors of the GNSS survey station in the non-earthquake period:
Figure FDA0003673517160000051
wherein the content of the first and second substances,
Figure FDA0003673517160000052
representing homomorphic errors in the displacement sequence of the kth station in the ith direction during the non-seismic period,
Figure FDA0003673517160000053
for the principal component corresponding to the jth mode in the ith direction in the non-seismic period calculated in step 3,
Figure FDA0003673517160000054
the k value j epsilon [1, n0 in the spatial response corresponding to the j mode in the i direction in the non-seismic period calculated in the step 3]，k∈[1,N]N is the number of GNSS stations, i ═ 1 indicates the east direction, i ═ 2 indicates the north direction, and i ═ 3 indicates the vertical direction;
and 6, deducting homomorphic errors of the east, north and vertical displacement sequences of the GNSS measuring station including the earthquake-affected station in the time period affected by the earthquake waves and homomorphic errors of the east, north and vertical displacement sequences of the GNSS measuring station in the non-earthquake period, which are obtained in the step 5, from the displacements of the N measuring stations:
Figure FDA0003673517160000055
Figure FDA0003673517160000056
wherein the content of the first and second substances,
Figure FDA0003673517160000057
representing homomorphic errors in the displacement sequence of the kth station in the ith direction for all observation times,
Figure FDA0003673517160000058
representing homomorphic errors in a displacement sequence of a kth measuring station including the measuring station in the ith direction in the time period influenced by the seismic waves,
Figure FDA0003673517160000059
representing the homomorphic error in the displacement sequence of the kth measuring station in the ith direction in the non-seismic period,
Figure FDA00036735171600000510
representing the original displacement sequence of the kth station in the ith direction, k is the [1, N ]]N is the number of GNSS stations, i ═ 1 indicates the east direction, i ═ 2 indicates the north direction, and i ═ 3 indicates the vertical direction;
step 6, the displacement sequence after homomorphic error is deducted
Figure FDA00036735171600000511
And (3) carrying out star daily filtering to finally obtain a high-precision position sequence:
Figure FDA00036735171600000512
wherein the content of the first and second substances,
Figure FDA00036735171600000513
the multipath error of the kth GNSS station calculated by the sidereal day filtering in the ith direction is represented as:
Figure FDA00036735171600000514
wherein the content of the first and second substances,
Figure FDA00036735171600000515
represents the displacement sequence of the kth station on the d day after deducting homomorphic error in the ith direction, and k belongs to [1, N ]]，d∈[1,D]N is the number of GNSS stationsD indicates the number of days involved in calculating multipath errors, i 1 indicates the east direction, i 2 indicates the north direction, and i 3 indicates the vertical direction.
2. The method of claim 1 for GNSS seismic surface displacement monitoring accounting for multipath and homomorphic errors, wherein: the east displacement sequence, the north displacement sequence and the vertical displacement sequence of the N GNSS survey stations in the step 1 are as follows:
Figure FDA0003673517160000061
k∈[1,N]
wherein the content of the first and second substances,
Figure FDA0003673517160000062
representing the east displacement sequence of the kth GNSS station,
Figure FDA0003673517160000063
representing a sequence of north displacements for the kth GNSS station,
Figure FDA0003673517160000064
representing a vertical displacement sequence of a kth GNSS observation station, wherein N represents the number of the GNSS observation stations;
step 1, calculating the energy of the seismic wave signal carried in the displacement sequence of the GNSS survey station by constructing an average energy function, wherein the energy of the seismic wave signal is as follows:
Figure FDA0003673517160000065
wherein the content of the first and second substances,
Figure FDA0003673517160000066
for the kth GNSS station at tepThe speed of the east direction of the moment,
Figure FDA0003673517160000067
for the kth GNSS station at tepThe speed of the north-facing direction at the moment,
Figure FDA0003673517160000068
for the kth GNSS station at tepVertical velocity of the moment, Velk(tep) Representing the kth GNSS station at tepThe integrated speed of the three directions of the time, delta t is the sampling interval of the observed data, m is the window length, Fk(tep) Representing the kth GNSS station at tepAverage energy of seismic signals carried in time-shifted sequences, F0For background noise, k ∈ [1, N ]]；
The kth station being in the ith directionepVelocity of time of day
Figure FDA0003673517160000069
The difference between the epochs of the displacement sequence can be used for obtaining:
Figure FDA00036735171600000610
wherein the content of the first and second substances,
Figure FDA00036735171600000611
represents the displacement of the kth station in the ith direction, delta t is the sampling interval of observed data, and k belongs to [1, N ]]Wherein, i-1 represents east direction, i-2 represents north direction, and i-3 represents vertical direction;
step 1, detecting the arrival time of the seismic waves of the GNSS survey station and the ending time of the seismic waves of the GNSS survey station:
wherein, when Fk(tep) Greater than 3 times F0When it is, consider tepThe time is the arrival time t of the seismic wave of the kth GNSS survey stationk,arrival；
When F is presentk(tep) Less than 3 times F0When it is, consider tepThe moment is the ending moment t of the seismic wave of the kth measuring stationk,stop。
3. The method of claim 1 for GNSS seismic surface displacement monitoring accounting for multipath and homomorphic errors, wherein: deducting the same-seismic deformation displacement in the sequence of east, north and vertical displacements of the post-seismic GNSS survey station in the step 2:
selecting the earliest arrival time of the seismic waves in all the GNSS survey stations as a seismic wave arrival time reference value:
tarrival＝min(tk,arrival)
selecting the latest seismic wave ending time in all GNSS survey stations as a seismic wave ending time reference value:
tstop＝max(tk,stop)
will observe the arrival time t1The time interval between the arrival time reference value of the seismic wave is used as a pre-earthquake observation time interval:
Tpre＝[t1,tarrival]
defining the time interval between the seismic wave arrival time reference value and the seismic wave ending time reference value as the time interval affected by the seismic waves:
Tmove＝[tarrival,tstop]
the reference value of the seismic wave ending time is changed to the observation ending time tendThe time period in between is taken as the post-earthquake observation time period:
Tafter＝[tstop,tend]
defining the pre-earthquake observation period and the post-earthquake observation period as non-earthquake periods:
Tpa＝Tpre∪Tafter
and deducting the deformation of the same earthquake from the east, north and vertical displacement sequence of the post-earthquake GNSS survey station:
Figure FDA0003673517160000071
wherein the content of the first and second substances,
Figure FDA0003673517160000072
represents the displacement sequence of the k GNSS observation station after the homoseismal deformation is deducted in the i direction,
Figure FDA0003673517160000073
representing the original displacement sequence of the kth GNSS station in the ith direction,
Figure FDA0003673517160000074
representing the homoseismal deformation of the kth GNSS station in the ith direction, k belongs to [1, N ∈]Where, 1 denotes an east direction, 2 denotes a north direction, and 3 denotes a vertical direction.
CN202011358334.8A 2020-11-27 2020-11-27 GNSS earthquake earth surface displacement monitoring method considering multipath and homomorphic errors Active CN112612045B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011358334.8A CN112612045B (en) 2020-11-27 2020-11-27 GNSS earthquake earth surface displacement monitoring method considering multipath and homomorphic errors

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011358334.8A CN112612045B (en) 2020-11-27 2020-11-27 GNSS earthquake earth surface displacement monitoring method considering multipath and homomorphic errors

Publications (2)

Publication Number Publication Date
CN112612045A CN112612045A (en) 2021-04-06
CN112612045B true CN112612045B (en) 2022-07-19

Family

ID=75227965

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011358334.8A Active CN112612045B (en) 2020-11-27 2020-11-27 GNSS earthquake earth surface displacement monitoring method considering multipath and homomorphic errors

Country Status (1)

Country Link
CN (1) CN112612045B (en)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8874477B2 (en) * 2005-10-04 2014-10-28 Steven Mark Hoffberg Multifactorial optimization system and method
CN111596321A (en) * 2020-05-29 2020-08-28 武汉大学 Multi-GNSS multi-path error star day filtering method and system using non-difference correction

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR2827670A1 (en) * 2001-07-18 2003-01-24 Osmos Sa Seismic monitoring system uses embedded measurement line array
JP6779084B2 (en) * 2016-09-30 2020-11-04 株式会社パスコ Visualization device of ground surface displacement in the area of interest and visualization program of ground surface displacement in the area of interest
CN106569231A (en) * 2016-11-10 2017-04-19 中国地震局第监测中心 Method for determining co-seismic displacement by using single GNSS receiver

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8874477B2 (en) * 2005-10-04 2014-10-28 Steven Mark Hoffberg Multifactorial optimization system and method
CN111596321A (en) * 2020-05-29 2020-08-28 武汉大学 Multi-GNSS multi-path error star day filtering method and system using non-difference correction

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Capturing coseismic displacement in real time with mixed single and dual-frequency receivers: application to the 2018 Mw7.9 Alaska earthquake;Kai Zheng 等;《GPS Solutions》;20191231;第1-14页 *
单历元PPP分析日本Mw9...的我国东部沿海同震地表形变;沈飞 等;《武汉大学学报· 信息科学版》;20121130;第37卷(第11期);第1345-1347页 *
基于PPP动态定位技术的同震地表形变分析;方荣新等;《武汉大学学报(信息科学版)》;20091105(第11期);第1340-1343页 *
高频GPS测定的汶川Ms8.0级地震震时近场地表变形过程;殷海涛 等;《科学通报》;20101231;第55卷(第26期);第2621-2626页 *

Also Published As

Publication number Publication date
CN112612045A (en) 2021-04-06

Similar Documents

Publication Publication Date Title
Lacroix et al. Location of seismic signals associated with microearthquakes and rockfalls on the Séchilienne landslide, French Alps
EP1839074B1 (en) Method of seismic signal processing
Lois et al. A new automatic S-onset detection technique: Application in local earthquake data
Myers et al. Absolute locations of the North Korean nuclear tests based on differential seismic arrival times and InSAR
Kummerow Using the value of the crosscorrelation coefficient to locate microseismic events
Binder et al. Convolutional neural networks for automated microseismic detection in downhole distributed acoustic sensing data and comparison to a surface geophone array
US10920582B2 (en) Systems and methods to use triangulation through one sensor beamforming in downhole leak detection
Barbour et al. Coseismic strains on Plate Boundary Observatory borehole strainmeters in southern California
CN111487678B (en) Analysis method for determining high-resolution small multichannel seismic minimum offset distance and system delay
CN111487589B (en) Target drop point positioning method based on multi-source sensor network
Roux et al. Microseismic activity within a serac zone in an alpine glacier (Glacier d’Argentiere, Mont Blanc, France)
US20140142854A1 (en) Method for locating a microseismic event
KR101923166B1 (en) Method for correcting amplitude magnitude of seismic signal extracted from seismic ambient noise
Zeng et al. Turning a telecom fiber‐optic cable into an ultradense seismic array for rapid postearthquake response in an urban area
KR101914657B1 (en) Method for extracting phase and amplitude information of seismic signal from seismic ambient noise
CN112612045B (en) GNSS earthquake earth surface displacement monitoring method considering multipath and homomorphic errors
CN109425891B (en) Fracture imaging quality detection method based on geological model
Khadhraoui et al. Real-time detection and localization of microseismic events
Plourde et al. Earthquake depths, focal mechanisms, and stress in the Lower St. Lawrence Seismic Zone
Su et al. Visual identity-based earthquake ground displacement testing method
Orfanos et al. Automatic passive seismic data processing with no prior information: the contribution of surface wave tomography
Kitov et al. Detection, estimation of magnitude, and relative location of weak aftershocks using waveform cross-correlation: The earthquake of August 7, 2016, in the town of Mariupol
CN109490952B (en) Seismic coherence analysis method and system
Guo et al. Aftershocks of the 2012 M w 8.6 Wharton Basin Intraplate Earthquake in the Eastern Indian Ocean Revealed by Near‐Field Ocean‐Bottom Seismometers
Moore‐Driskell et al. Integration of arrival‐time datasets for consistent quality control: A case study of amphibious experiments along the Middle America Trench

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant