GNSS aseismic slip signal extraction method based on template matching and linear velocity correction

Abstract

The present invention discloses a GNSS aseismic slip signal extraction method based on template matching and linear rate correction, which relates to the field of geodetic surveying and geophysical monitoring technology, including: decomposing the mixed noise in the GNSS / InSAR time series and extracting the denoised residual signal; using masking technology to separate transient and non-transient signal segments, dynamically correcting the linear trend error based on the variance-weighted average rate of the non-transient data before and after the transient segment, and achieving millimeter-level amplitude extraction of transient signals through iterative optimization; constructing an exponential decay and cosine-type transient signal template library, calculating the full sequence cross-correlation by sliding the template along the residual sequence, and accurately identifying the signal amplitude, duration, and spatial distribution. The present invention breaks through the bottleneck of linear rate interference and multi-signal superposition technology, provides high-precision data support for fault slip behavior, aseismic slip identification, and earthquake risk assessment, and significantly enhances the engineering application value of GNSS / InSAR technology in geophysical research and disaster warning.

Description

Technical Field

[0001] The present invention relates to the technical field of geodesy and geophysical monitoring, and in particular to a GNSS aseismic slip signal extraction method based on template matching and linear rate correction. Background Art

[0002] Transient aseismic slip encompasses a variety of crustal deformation phenomena, including slow slip events (SSEs), fault creep, and post-seismic afterslip. Its precise identification and accurate extraction of characteristic parameters are of vital scientific significance and application value for understanding the evolution of locked fault states, revealing diverse fault slip behaviors, searching for potential earthquake precursor signals, and conducting reliable earthquake risk assessments. These transient events play an integral role in the complete cycle of earthquake energy accumulation and release, significantly altering the stress state along the fault and surrounding areas, and thus influencing the development and occurrence of future earthquakes. Currently, Global Navigation Satellite System (GNSS) and Interferometric Synthetic Aperture Radar (InSAR) are the primary geodetic techniques for monitoring these subtle surface deformations and extracting transient aseismic slip signals. GNSS technology provides time series of station coordinates with high temporal resolution (days or even higher), accurately capturing slow surface displacements at the millimeter level. InSAR technology, with its high spatial resolution and wide coverage, provides a unique means for capturing the spatial distribution of transient deformation fields. However, accurately extracting the target transient signal from these geodetic time series faces numerous challenges. The original time series are often contaminated by long-term tectonic trends, seasonal loading effects, co-seismic and post-seismic deformation, instrument or reference frame errors (such as common-mode error), and various random noise factors. In standard data processing, linear trend correction is a common step to remove the long-term tectonic background. However, this process, particularly when using global or static linear models, often inappropriately absorbs some of the transient signal energy, significantly underestimating the amplitude of the extracted transient aseismic slip signal and potentially distorting its true shape and duration. This signal amplitude bias introduced by improper linear trend correction can have a series of negative consequences on subsequent geophysical interpretations. For example, underestimated slip amplitude directly leads to inaccurate estimates of strain energy and seismic moment released during aseismic slip, which in turn affects the accurate assessment of fault energy distribution (i.e., the ratio of energy released during aseismic to seismic slip). Furthermore, inaccurate transient signal parameters can mislead assessments of fault coupling, the potential for future major earthquakes, and the seismic hazard of specific regions. Internationally, scholars have noted the complexity of signal separation in geodetic time series processing. For example, in studying GNSS time series, some studies have focused on improving noise models (e.g., accounting for colored noise) to obtain more accurate rate estimates. Furthermore, some studies have employed iterative methods to optimize between slow slip event detection and linear trend correction during interseismic periods, aiming to obtain more reliable background rates and transient signals.For example, using low-frequency earthquakes (LFEs) and tectonic tremor as a guide, previously unidentified small slow-slip events were extracted from a noisy background, revealing that interseismic periods are not entirely quiescent. These studies highlight the importance of accurately separating long-term trends from transient signals. However, most existing methods focus on optimizing the linear trend estimation itself or fitting multiple signal components simultaneously through complex models, lacking direct attention to and systematic correction for the inherent bias in the amplitude of identified transient signals caused by standard linear detrending operations. In summary, despite continuous advancements in geodetic data processing technology, significant technical bottlenecks remain in the accurate extraction of transient anesthetic slip signals, particularly in effectively overcoming and correcting the underestimation of signal amplitudes caused by conventional linear trend correction. This bias directly affects the accurate interpretation of fault slip behavior and energy distribution, limiting our in-depth understanding of earthquake physics and the accuracy of earthquake risk assessment. To address this issue, we propose a GNSS anesthetic slip signal extraction method based on template matching and linear velocity correction. Summary of the Invention

[0003] In order to solve the above technical problems, a GNSS aseismic slip signal extraction method based on template matching and linear rate correction is provided. This technical solution solves the above-mentioned deficiencies and limitations in the existing technology for identifying transient aseismic slip signals based on GNSS or InSAR time series.

[0004] In order to achieve the above objects, the technical solution adopted by the present invention is:

[0005] The GNSS aseismic slip signal extraction method based on template matching and linear rate correction includes the following steps:

[0006] S1. Obtain the original GNSS or InSAR time series observation data in the fault area and decompose them into seasonal signals, co-seismic step terms, post-seismic deformation terms, initial linear trend terms, common mode error terms, and residual signals using a mixed noise model.

[0007] S2. Preliminarily dividing the transient signal segment into a transient signal segment and a non-transient signal segment based on the residual signal; for each transient signal segment, performing the following steps: calculating a weighted average modified rate based on the linear rate and variance of the preceding and following non-transient signal segments; if only the preceding or following segment exists, directly using the rate of that segment; iteratively subtracting the difference between the modified rate and the initial rate from the transient signal segment data to optimize the residual signal;

[0008] S3. Construct a transient signal template library, which includes an exponential decay template and a cosine template; calculate the normalized cross-correlation function between the optimized residual signal and the template: the numerator is the sum of the point-by-point products of the two signals within the template length window, and the denominator is the square root of the product of the sum of the squares of the two signals; identify the amplitude, start time and duration of the transient signal by searching for the maximum value of the cross-correlation function, and determine it as a valid signal when the normalized cross-correlation coefficient is greater than 0.7, and compare it with the recognition result of the original residual signal.

[0009] Preferably, template matching improves the signal-to-noise ratio through normalized cross-correlation analysis, which is specifically manifested as follows: for transient signals with an amplitude less than three times the noise level, the peak value of the cross-correlation coefficient is increased by forty to sixty percent compared with conventional correlation analysis; the template library sets a weight coefficient, the weight of the exponential decay template is 0.6, and the weight of the cosine template is 0.4, and the final recognition parameter takes the weighted comprehensive value of the two template results; the cross-correlation calculation adopts a frequency domain acceleration method, and the time domain convolution operation is converted into a frequency domain product operation through Fourier transform, and the calculation efficiency is improved by an order of magnitude of the data volume N multiplied by the logarithm of N.

[0010] Preferably, the seasonal signal is composed of harmonic components of an annual cycle and a semi-annual cycle; the co-seismic step term is obtained by accumulating the step displacements of all known earthquake events; the post-seismic deformation term is obtained by accumulating the exponential decay displacements of significant post-seismic deformation events; the initial linear trend term is expressed as the product of the initial estimated long-term linear rate and time; the common-mode error term is extracted using principal component analysis or independent component analysis; and the residual signal is the remainder after subtracting the above five items from the original data.

[0011] Preferably, in step S1, the seasonal signal contains only two harmonic components, namely, an annual cycle and a semi-annual cycle; the harmonic amplitude is estimated by a Lomb-Scargle periodogram algorithm; the time variable is measured in Julian years, with a minimum time resolution of 0.001 years;

[0012] The earthquake event time of the coseismic step term is derived from the earthquake catalog database, and the time error is controlled within 10 seconds; the step amplitude is calculated by the difference of the displacement mean in the three-day window before and after the event, and the influence of the common mode error is deducted in the calculation;

[0013] The amplitude coefficient of the post-seismic deformation term is determined by fitting the logarithmic time function, and the decay constant is obtained by searching the residual variance minimization criterion.

[0014] Preferably, in step S1, the rules for introducing the post-seismic deformation term are as follows: for earthquakes with a moment magnitude of 7.0 or above, the post-seismic deformation term is forcibly retained; for earthquakes with a moment magnitude of 6.5 or below, it is introduced only when the deformation amount one year after the earthquake exceeds twice the noise level; the co-seismic step amplitude and the post-seismic deformation amplitude coefficient of the same earthquake event are jointly solved, and the constraint condition is that the total displacement at the time point one year after the earthquake is equal to the co-seismic step amplitude plus the post-seismic deformation amplitude multiplied by one minus the inverse of the negative attenuation constant of the natural constant e.

[0015] Preferably, in step S1, when the common-mode error term is extracted using principal component analysis, the principal components with a cumulative contribution rate of eigenvalues ​​of not less than 85% are retained; when independent component analysis is used, the FastICA algorithm is used to optimize non-Gaussianity; before extraction, all time series of the measuring station network need to be standardized and preprocessed to remove the linear trends specific to each station; the extraction results need to pass the spatial uniformity test, and the test standard is that the inter-station difference in the annual change of the common-mode error of all measuring stations is less than 0.5 mm.

[0016] Preferably, in step S2, the transient signal segmentation adopts a displacement change rate threshold method, and when the displacement change in a continuous time window exceeds three times the standard deviation of the background noise, it is marked as a transient signal segment; the weighted average correction rate is calculated by: dividing the linear rate of the first segment by the variance of the first segment, adding the linear rate of the second segment divided by the variance of the second segment, and then dividing by the sum of the two inverses of the variance; the variance is obtained by least squares linear fitting of non-transient data before and after the transient signal segment; the iterative operation is specifically: at each time point in the transient signal segment, subtracting the difference between the product of the correction rate and time and the initial linear trend term from the residual value to generate a new residual sequence, and repeating the segmentation and rate calculation until the change in the correction rate between two adjacent times is less than 0.01 mm per year;

[0017] The least squares linear fitting adopts the weighted least squares method, and the observation weight is inversely proportional to the noise variance of the time series; the noise variance is obtained by calculating the residual autocorrelation function through the autoregressive moving average model, and the model order is fixed to 1st-order autoregressive and 1st-order moving average; the iteration stopping condition is supplemented with a dynamic detection mechanism: the standard deviation of the correction rate change in three consecutive iterations is less than 0.005 mm per year; the transient signal segment boundary adjustment rule is: if the second-order derivative value at the endpoint of the corrected transient signal segment exceeds the threshold of 0.1 mm per day squared, the mask range is expanded until the derivative value is continuously smooth.

[0018] Preferably, in step S3, the exponential decay template is defined as an amplitude parameter multiplied by a negative exponential function with a starting time as a reference and a duration as a decay constant, wherein the amplitude parameter ranges from 0.5 mm to 20 mm, and the duration parameter ranges from 10 days to 365 days, and is used to describe the fault slip relaxation process;

[0019] The cosine template has an amplitude parameter range of 1 mm to 15 mm and a duration parameter range of 30 days to 180 days, and is used to describe periodic slow slip events;

[0020] During the template matching process, the template is slid along the time axis with a step length of 1 day, and the duration is traversed within the preset range at intervals of 5 days; the normalized cross-correlation coefficient is calculated for each sliding position.

[0021] Preferably, in step S3, the denominator calculation of the normalized cross-correlation function includes: calculating the cumulative sum of the square values ​​of the corrected residual signal of all data points in the template length time window, calculating the cumulative sum of the square values ​​of the template signal of all data points in the same window, multiplying the two cumulative sums and taking the square root; the numerator calculation is the cumulative sum of the point-by-point products of the corrected residual signal and the template signal in the same window; the specific process of parameter identification is: moving the template with a sliding step on the complete time axis, recording the sliding time point corresponding to the maximum value of the cross-correlation coefficient as the signal start time, when the maximum value is greater than 0.7, the template duration at this position is used as the identification duration, and the signal amplitude is calculated by the amplitude ratio relationship between the template and the residual signal; the identification result must meet the physical constraints that the interval between adjacent transient events is greater than the duration and the spatial distribution conforms to the fault direction.

[0022] Preferably, in step S3, the comparison of the recognition results includes: recording the transient signal amplitude values ​​identified based on the optimized residual and the original residual respectively; when the absolute value of the relative deviation of the amplitudes of the two exceeds 20%, triggering the secondary iterative optimization of S2; the parameter consistency verification conditions are: the absolute value of the start time deviation is less than 3 days, the relative deviation of the duration is less than 10%, and the relative deviation of the amplitude is less than 15%; finally, the recognition result using the optimized residual is output, and the amplitude deviation percentage of the original residual recognition result is marked in the result file.

[0023] Compared with the prior art, the present invention has the following beneficial effects:

[0024] The GNSS aseismic slip signal extraction method proposed in the present invention performs fine processing on geodetic time series such as GNSS / InSAR, separates transient and non-transient signal segments using masking technology, and dynamically corrects the linear trend error based on the variance-weighted average rate of non-transient data before and after the transient segment, and accurately extracts the transient signal amplitude through iterative optimization; combines the preset transient signal template library for sliding cross-correlation analysis, and can selectively incorporate fault geometry and multi-station joint constraints to achieve accurate identification of transient aseismic slip signal amplitude and duration, thereby more realistically reflecting the actual aseismic slip physical process of the fault, and providing quantitative It provides a highly reliable data basis for analyzing fault energy distribution, accurately depicting the spatiotemporal evolution of fault coupling states, and deeply understanding the stress loading and release mechanisms. It can effectively reduce geophysical interpretation deviations caused by inaccurate signal parameter extraction, and is of great value in improving the application efficiency and scientific output of existing GNSS, InSAR and other observation data in the fields of fault activity monitoring, slow earthquake phenomenon research, earthquake precursor exploration, and dynamic assessment of earthquake hazards. It has also actively promoted the development of related geophysical research and disaster warning technologies towards higher precision and greater intelligence, and has significant scientific significance and broad application prospects. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Figure 1 is a flow chart of the method of the present invention;

[0026] Figure 2 This is a graph of synthetic time series data;

[0027] Figure 3 is the residual linear rate correction map;

[0028] Figure 4 The SSE amplitude and correlation coefficient of the transient signal for template matching Figure 1 ;

[0029] Figure 5 The SSE amplitude and correlation coefficient of the transient signal for template matching Figure 2 ;

[0030] Figure 6 The SSE amplitude and correlation coefficient of the transient signal for template matching Figure 3 ;

[0031] Figure 7 The SSE amplitude and correlation coefficient of the transient signal for template matching Figure 4 . DETAILED DESCRIPTION

[0032] The following description is intended to disclose the present invention so that those skilled in the art can implement the present invention. The preferred embodiments described below are merely examples, and those skilled in the art may conceive of other obvious variations.

[0033] The method proposed in the present invention mainly includes the following three core steps, which have been summarized in the "Summary of the Invention" section. Here, we will provide a more detailed explanation based on considerations in actual operation. The present invention will be further explained below with reference to the accompanying drawings and examples.

[0034] GNSS / InSAR time series data acquisition and decomposition:

[0035] This step aims to obtain and preprocess time series data containing surface deformation information and decompose it into multiple known signal components and target residual signals.

[0036] First, obtain GNSS or InSAR time series data related to fault activity in the study area. GNSS data typically originate from a network of continuously operating reference stations, providing high-resolution three-dimensional coordinate time series. InSAR data, on the other hand, can be acquired from SAR satellites and processed using time-series InSAR processing techniques to generate deformation time series.

[0037] The original observation time series data obtained It is usually a mixture of multiple signals and noise. According to step S1 of "Summary of the Invention", the following mixed noise model is used to decompose Y(t):

[0038] ;

[0039] Among them, each item is time The specific definition and processing method of the function are as follows:

[0040] Seasonal signal S(t): Mainly caused by hydrological loads, atmospheric pressure changes, etc., usually showing annual and semi-annual cycle characteristics. Its mathematical expression is:

[0041] ;

[0042] in, is the value of the seasonal signal at time t. Is the summation symbol, indicating the harmonic order To the highest order Accumulate the subsequent items (for annual and semi-annual periods, Usually 2 is selected). and They are The amplitudes of the sine and cosine components of the first harmonic. and are the sine and cosine function parts of the kth order harmonic, respectively, where Usually in years, is the ratio of pi. These amplitude parameters can be estimated and removed from the time series using the least squares fitting method.

[0043] coseismic step term : The instantaneous displacement caused by an earthquake. It is defined as:

[0044] ;

[0045] in, It's time The cumulative effect of the coseismic step at represents the accumulation of all known co-seismic events j. is a step function (Heaviside function), when (the moment when the jth co-seismic event occurs), its value is 1, otherwise it is 0. It is The displacement step amplitude caused by a coseismic event. and It can be obtained from earthquake catalogs or preliminary data analysis and fixed or estimated in the model.

[0046] Post-seismic deformation term P(t): A major earthquake is usually followed by post-seismic deformation for several months to several years, which is often described by a logarithmic or exponential decay model:

[0047] ;

[0048] in, It's time Cumulative effect of post-earthquake deformation. Represents all earthquake events that produce significant post-seismic deformation Accumulate. It is The amplitude coefficient of the post-seismic deformation caused by an earthquake event. is the decay constant of the post-seismic deformation caused by the event, indicating how quickly the post-seismic deformation decays. These parameters are usually estimated and removed using nonlinear fitting methods for earthquake events known to produce significant post-seismic deformation.

[0049] Initial linear trend term L(t): represents the long-term tectonic movement background rate. It is defined as:

[0050] ;

[0051] is the initial linear trend value at time t, where is the initial estimate of the long-term linear rate. This is a preliminary, global rate, typically derived through least squares fitting of the entire time series (or a selected "quiet" period). It is this global or static linear trend removal that can lead to an underestimation of the subsequent transient signal amplitude, which is the problem addressed by step S2 of the present invention.

[0052] Common mode error term E(t): This term is prevalent in GNSS regional networks and is caused by reference frame instability, satellite orbit errors, etc. It can be extracted and subtracted from the time series data of the station network using spatial filtering methods such as principal component analysis (PCA) or independent component analysis (ICA).

[0053] Residual signal R(t): The signal obtained by removing all the above known signals and error terms from the original time series Y(t), that is:

[0054] ;

[0055] The residual signal R(t) mainly contains the target transient aseismic slip signal and unmodeled random noise, and is the core object of subsequent analysis.

[0056] Dynamic linear rate correction of transient signal segments in residual signal

[0057] This step is one of the key innovations of the present invention, and is intended to correct the initial linear trend term L(t) in step S1 to remove the possible suppression effect on the transient signal amplitude.

[0058] Preliminary determination and masking of transient signal segments:

[0059] Based on the residual signal R(t) obtained in step S1, it is first necessary to preliminarily identify time periods that may contain transient anechoic slip events. This can be achieved through visual inspection, setting simple amplitude or rate change thresholds, or referring to a catalog of known slow slip events (if relevant information is available for the study area). Once transient signal segments are identified, they are distinguished from non-transient signal segments using masking techniques.

[0060] Construction and application of segmented dynamic correction model:

[0061] For each identified transient signal segment, the core idea is to use the weighted deformation rates of the preceding and following non-transient segments adjacent to the transient segment to define a more accurate correction rate .

[0062] Calculate the weighted average rate of the non-transient segments before and after the transient signal segment :

[0063] ;

[0064] in, 、 are the linear rates obtained by performing linear least squares fitting on the residual signal R(t) of the non-transient segment before and after the transient segment, respectively; 、 are the variances of the two rate estimates, used as weights. If there are only non-transient segments before or after the transient segment (for example, the transient event occurs at the beginning or end of the time series), it degenerates to or .

[0065] this Represents the expected rate of correction during this time period in the absence of transient events.

[0066] Iterative correction and dynamic optimization:

[0067] Use the calculated To adjust the baseline of the transient signal segment. Specifically, for the time point within the transient signal segment, subtract the value of The difference between the local linear trend defined by the initial global rate and the linear trend defined by the initial global rate is equivalent to correcting the original residual signal R(t) to obtain the corrected residual signal .

[0068] The advantage of this "dynamic" correction is that it provides a personalized background rate correction for each transient event, preventing the global linear trend from over-absorbently absorbing the local transient signal, thereby more accurately recovering the true amplitude of the transient signal. This is a direct solution to the problem of signal amplitude underestimation described in the background.

[0069] Parameter extraction of transient aseismic slip signal based on template matching

[0070] After obtaining the residual signal after dynamic linear rate correction , a template matching method is used to identify and parameterize the transient aseismic slip signal.

[0071] Transient signal template library construction:

[0072] A set of mathematical functions that can describe the typical transient aseismic slip signal morphology is predefined as templates. Commonly used templates include:

[0073] Exponential decay The template is defined as:

[0074] ;

[0075] is the template value of the exponentially decaying transient signal at time t. is the amplitude of the transient signal, is the start time of the transient signal, is the duration of the transient signal.

[0076] Cosine template Defined as:

[0077] ;

[0078] is the template value of the cosine transient signal at time t. The template library can contain multiple function forms to adapt to different types of transient signals.

[0079] Normalized cross-correlation analysis:

[0080] Compare each template function in the template library with the corrected residual signal Perform normalized cross-correlation calculation. Normalized cross-correlation function Defined as:

[0081] ;

[0082] is the normalized cross-correlation coefficient in time The value at , its range is usually [-1,1], the length of the transient signal template is .molecular is the corrected residual signal With template signal The inner product of and Represents signals and the time-shifted template The energy within this time window. Through this normalization, the similarity of the two signal waveforms can be measured without being affected by their absolute amplitudes. According to the normalized time series cross-correlation function, the maximum correlation coefficient is identified to identify the amplitude A and start time of the transient aseismic slip signal. , duration D, and compare the results obtained by normalizing the correlation coefficient with the original residual signal. The template and its parameters with the highest correlation and reasonable physical meaning are used as the features of the identified transient aseismic slip signal.

[0083] A key validation step is to apply this method to the corrected residuals The obtained results are compared with those obtained by directly applying them to the original residual R(t) without correction in step S2. This comparison allows to quantify the contribution of the dynamic linear rate correction in step S2 to the improvement of the accuracy of transient signal parameter extraction (especially amplitude A).

[0084] Simulation experiment

[0085] To verify the effectiveness and superiority of the method proposed in the present invention, the following simulation experiment was designed and conducted. This embodiment aims to demonstrate the performance of the method under known signal parameters, especially its ability to accurately recover the transient signal amplitude.

[0086] Reference Figure 2 The concept shown generates a synthetic time series Y(t) that simulates a typical GNSS station displacement time series and contains the following components:

[0087] Real transient aseismic slip signals (SSEs): embed one or more transient signals with known parameters. In this simulation, three SSE events with different characteristics are embedded, and their real parameters (amplitude , start time , duration )exist Figure 2 The winning bid was awarded.

[0088] Long-term linear signal: Set a known signal with different long-term linear rates and different noise levels to examine the model's sensitivity to different noise and long-term linear rates. Figure 4 The middle bars represent different combinations of noise level and long-term linear rate.

[0089] Seasonal signal: Contains sine waves with an annual cycle (amplitude 10 mm) and a semi-annual cycle (amplitude 5 mm), with set amplitude and phase.

[0090] Noise signal: Add Gaussian white noise or colored noise (such as a combination of flicker noise and white noise) to simulate the random errors in real observation data.

[0091] Synthetic time series Y(t): Add all the above components to get the final simulated observation data.

[0092] The three-step method proposed in this invention is applied to the synthetic time series Y(t) generated above ( Figure 2 ), in the initial residual signal R(t), transient signal segments are preliminarily identified based on the known approximate occurrence period of SSEs. For each transient signal segment, a weighted average rate is calculated using the data of the preceding and / or following non-transient segments. A linear rate correction is then performed on the corresponding transient signal segment to obtain the corrected residual signal. Figure 3 The correction process is intuitively demonstrated: different line segments represent R(t), the local correction trend estimated based on the non-transient segment, and the scattered points represent the transient aseismic slip signal after the dynamic linear rate correction. Then, the initial residual signal R(t) (i.e., the residual without dynamic rate correction) and the corrected residual signal are respectively (i.e. the residual after dynamic rate correction) Apply template matching method. Find the best start time by normalized cross correlation scanning and duration D, and fit the amplitude A. Figure 4 Shows the relationship between R(t) and Comparison of the SSE amplitude recovery percentage (i.e., the ratio of the difference between the estimated amplitude and the true value to the true value) obtained after template matching. The results show that the amplitude loss ratio of the transient signal estimated based on the corrected residual is much smaller than the amplitude loss ratio without correction. Figure 5 、 Figure 6 、 Figure 7 The three transient signal segments show the The detailed process of template matching: the solid lines are the normalized correlation coefficient and the cumulative correlation coefficient, and the dotted lines are Sequence, the shaded area indicates the start and duration of the real signal.

[0093] Result analysis:

[0094] Amplitude estimation accuracy improvement: from Figure 4 It can be seen that for the three simulated transient events SSE1, SSE2, and SSE3, the amplitude loss ratio extracted by performing template matching directly on the residual R(t) without dynamic rate correction can reach 80%. This confirms the problem that the traditional linear detrending method described in the background art will lead to the underestimated amplitude of transient signals. In contrast, under different noise levels and long-term linear rates, after adopting the dynamic linear rate correction in step S2 of the present invention, the corrected residual The extracted amplitude is very close to the true value (the amplitude loss ratio is mostly within 10%), which fully demonstrates the significant advantage of the method of the present invention in restoring the true amplitude of transient signals.

[0095] Effectiveness of template matching: Figure 5 、 Figure 6 、 Figure 7 It shows that for the residual sequence after dynamic rate correction ,The template signal is in good agreement with the identified transient signal segment, and the normalized correlation coefficient reaches its peak during the period when the real signal occurs, indicating that the template matching algorithm can effectively identify and locate transient events from background noise.

[0096] Overall Performance: Simulation experimental results demonstrate that the proposed "template matching aseismic slip signal extraction method based on linear velocity correction" can effectively separate and accurately extract various parameters of transient aseismic slip signals from complex synthetic time series. In particular, the dynamic linear velocity correction in step S2 significantly overcomes the signal amplitude underestimation problem caused by improper linear trend removal in conventional processing.

[0097] These results directly support the advantages of the present invention described in the “Beneficial Effects” section: it can more realistically reflect the actual aseismic slip physical process of the fault, provide a more reliable data basis for subsequent geophysical interpretation (such as fault energy distribution and coupling state assessment), and reduce the deviation caused by inaccurate signal parameter extraction.

[0098] In summary, through the description of the above detailed steps and the verification of simulation experiments, the present invention provides an effective method for extracting transient aseismic slip signals. This method can significantly improve the extraction accuracy of parameters such as transient signal amplitude by introducing a dynamic linear rate correction mechanism for transient signal segments and combining it with template matching technology. The simulation experimental results clearly demonstrate the improvement of this method over the traditional processing flow, especially its superiority in correcting the underestimated signal amplitude. The process of the present invention is clear and does not rely on complex prior models. It has important practical value and scientific significance for accurately extracting and analyzing weak transient crustal deformation signals from geodetic time series such as GNSS and InSAR, and helps to deepen the understanding of fault activity and earthquake physical processes. The feasibility and reliability of this method have been confirmed by simulation data.

[0099] The above shows and describes the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The above embodiments and descriptions merely illustrate the principles of the present invention. Various changes and modifications may be made to the present invention without departing from the spirit and scope of the present invention. Such changes and modifications are intended to fall within the scope of the present invention. The scope of protection claimed by the present invention is defined by the appended claims and their equivalents.

Claims

1. A GNSS aseismic slip signal extraction method based on template matching and linear rate correction, characterized in that: The following steps are involved: S1. Obtain the original GNSS or InSAR time series observation data in the fault area and decompose them into seasonal signals, co-seismic step terms, post-seismic deformation terms, initial linear trend terms, common mode error terms, and residual signals using a mixed noise model. S2. Preliminarily dividing the transient signal segment and the non-transient signal segment based on the residual signal; For each transient signal segment, the following steps are performed: Calculate the weighted average correction rate based on the linear rate and variance of the preceding and following non-transient signal segments; If only the preceding or following segment exists, use the rate of that segment directly; Iteratively subtract the correction term obtained by multiplying the correction rate by time from the entire time series data to optimize the residual signal; S3, constructing a transient signal template library, wherein the template library includes an exponential decay template and a cosine template; The normalized cross-correlation function of the optimized residual signal and the template is calculated: the numerator is the sum of the point-by-point products of the two signals within the template length window, and the denominator is the square root of the product of the sum of the squares of the two signals. The amplitude, start time and duration of the transient signal are identified by searching for the maximum value of the cross-correlation function. When the normalized cross-correlation coefficient is greater than 0.7, it is determined to be a valid signal and compared with the recognition result of the original residual signal.

2. The GNSS aseismic slip signal extraction method based on template matching and linear velocity correction according to claim 1, characterized in that: Template matching improves the signal-to-noise ratio through normalized cross-correlation analysis. Specifically, for transient signals with amplitudes less than the noise level, the peak value of the cross-correlation coefficient is 40 to 60 percent higher than that of conventional correlation analysis. The template library sets weight coefficients, with the exponential decay template having a weight of 0.4 and the cosine template having a weight of 0.

6. The final recognition parameter is the weighted combined value of the two template results. The cross-correlation calculation adopts the frequency domain acceleration method, which converts the time domain convolution operation into the frequency domain product operation through Fourier transform. The calculation efficiency is improved by an order of magnitude of the data volume N multiplied by the logarithm of N.

3. The GNSS aseismic slip signal extraction method based on template matching and linear velocity correction according to claim 1, characterized in that: The seasonal signal is composed of harmonic components of annual and semi-annual cycles; the co-seismic step term is obtained by accumulating the step displacements of all known earthquake events; the post-seismic deformation term is obtained by accumulating the exponentially decaying displacements of significant post-seismic deformation events; the initial linear trend term is expressed as the product of the initial estimated long-term linear rate and time; the common-mode error term is extracted using principal component analysis or independent component analysis; the above items can be used to obtain the optimal prediction model through least squares fitting, and finally the residual signal is the optimal prediction model of the original data minus the above five items.

4. The GNSS aseismic slip signal extraction method based on template matching and linear velocity correction according to claim 1, characterized in that: In step S1, the seasonal signal contains only two harmonic components, namely, an annual cycle and a semi-annual cycle; the harmonic amplitudes are estimated using a Lomb-Scargle periodogram algorithm; the time variable is measured in days, and the minimum time resolution has no lower limit; The earthquake event time of the coseismic step term is derived from the earthquake catalog database or instrument maintenance catalog, and the time error is controlled within 10 seconds; the step amplitude is calculated by the difference of the displacement mean in the three-day window before and after the event time, and the influence of the common mode error is deducted in the calculation; The amplitude coefficient of the post-seismic deformation term is determined by fitting the logarithmic time function, and the decay constant is obtained by searching the residual variance minimization criterion.

5. The GNSS aseismic slip signal extraction method based on template matching and linear velocity correction according to claim 1, characterized in that: In step S1, the rules for introducing the post-seismic deformation term are as follows: for earthquakes with a moment magnitude of 7.0 or above, the post-seismic deformation term is forcibly retained; for earthquakes with a moment magnitude of 6.5 or below, it is introduced only when the deformation one year after the earthquake exceeds twice the noise level; the co-seismic step amplitude and the post-seismic deformation amplitude coefficient of the same earthquake event are jointly solved, and the constraint condition is that the total displacement at the time point one year after the earthquake is equal to the co-seismic step amplitude plus the post-seismic deformation amplitude multiplied by one minus the inverse of the negative attenuation constant of the natural constant e.

6. The GNSS aseismic slip signal extraction method based on template matching and linear velocity correction according to claim 1, characterized in that: In step S1, when the common-mode error term is extracted using principal component analysis, the principal components with a cumulative contribution rate of eigenvalues ​​of not less than 85% are retained; when independent component analysis is used, the FastICA algorithm is used to optimize non-Gaussianity; before extraction, all time series of the measuring station network need to be standardized and preprocessed to remove the linear trends specific to each station; the extracted results need to pass the spatial uniformity test, and the test standard is that the inter-station difference in the annual change of the common-mode error of all measuring stations is less than 0.5 mm.

7. The GNSS aseismic slip signal extraction method based on template matching and linear velocity correction according to claim 1, characterized in that: In step S2, the transient signal segmentation adopts the displacement change rate threshold method. When the displacement change in the continuous time window exceeds three times the standard deviation of the background noise, it is marked as a transient signal segment. The weighted average correction rate is calculated by dividing the linear rate of the first segment by the variance of the first segment, adding the linear rate of the second segment divided by the variance of the second segment, and then dividing by the sum of the inverses of the two variances. The variance is obtained by least squares linear fitting of the non-transient data before and after the transient signal segment. The iterative operation is specifically as follows: at each time point in the transient signal segment, subtract the difference between the product of the correction rate and time and the initial linear trend term from the residual value to generate a new residual sequence, and repeat the segmentation and rate calculation until the change in the correction rate between two adjacent times is less than 0.01 mm / year. The least squares linear fitting adopts the weighted least squares method, and the observation weight is inversely proportional to the variance of the time series noise; The noise variance is obtained by calculating the residual autocorrelation function using an autoregressive moving average model, and the model order is fixed at 1st-order autoregressive and 1st-order moving average; the iteration stopping condition is supplemented by a dynamic detection mechanism: the standard deviation of the change in the correction rate in three consecutive iterations is less than 0.005 mm per year; the transient signal segment boundary adjustment rule is: if the second-order derivative value at the endpoint of the corrected transient signal segment exceeds the threshold of 0.1 mm per day squared, the mask range is expanded until the derivative value is continuously smooth.

8. The GNSS aseismic slip signal extraction method based on template matching and linear velocity correction according to claim 1, characterized in that: In step S3, the exponential decay template is defined as an amplitude parameter multiplied by a negative exponential function with a starting time as a reference and a duration as a decay constant, wherein the amplitude parameter ranges from 1 mm to 20 mm and the duration parameter ranges from 10 days to 1000 days, and is used to describe the fault slip relaxation process; The cosine template has an amplitude parameter range of 1 mm to 20 mm and a duration parameter range of 10 days to 1000 days, and is used to describe periodic slow slip events; During the template matching process, the template is slid along the time axis with a step length of 1 day, and the duration is traversed within the preset range at intervals of 5 days; the normalized cross-correlation coefficient is calculated for each sliding position.

9. The GNSS aseismic slip signal extraction method based on template matching and linear velocity correction according to claim 1, characterized in that: In step S3, the denominator calculation of the normalized cross-correlation function includes: calculating the cumulative sum of the square values ​​of all data points of the corrected residual signal in the template length time window, calculating the cumulative sum of the square values ​​of all data points of the template signal in the same window, multiplying the two cumulative sums and taking the square root; the numerator calculation is the cumulative sum of the point-by-point products of the corrected residual signal and the template signal in the same window; the specific process of parameter identification is: moving the template with a sliding step on the complete time axis, recording the sliding time point corresponding to the maximum value of the cross-correlation coefficient as the signal start time, when the maximum value is greater than 0.7, the template duration at this position is used as the identification duration, and the signal amplitude is calculated by the amplitude ratio relationship between the template and the residual signal; the identification result must meet the physical constraints that the interval between adjacent transient events is greater than the duration and the spatial distribution conforms to the fault direction.

10. The GNSS aseismic slip signal extraction method based on template matching and linear velocity correction according to claim 1, characterized in that: In step S3, the comparison of the recognition results includes: recording the transient signal amplitude values ​​based on the optimized residual and the original residual respectively; when the absolute value of the relative deviation of the amplitudes of the two exceeds 20%, triggering the secondary iterative optimization of S2; the parameter consistency verification conditions are: the absolute value of the start time deviation is less than 20 days, the relative deviation of the duration is less than 10%, and the relative deviation of the amplitude is less than 15%; finally, the recognition result using the optimized residual is output, and the amplitude deviation percentage of the original residual recognition result is marked in the result file.

Images