Broadband seismograph performance online evaluation method based on spatial-temporal characteristics of multiple stations

Abstract

The invention discloses a broadband seismograph performance online evaluation method based on spatial-temporal characteristics of multiple stations, and aims to realize real-time monitoring and accurate evaluation of the performance state of a seismograph. The method comprises the following steps: carrying out stationarity test on seismic background noise data of a plurality of stations, removing unstable data, carrying out preprocessing, dividing the processed data into continuous data segments by adopting a sliding window technology, and providing a time sequence data stream for subsequent analysis; multiple time-frequency features are extracted through continuous wavelet transform, Hilbert-Huang transform and third-order statistics, information entropy quantification features and principal component analysis dimension reduction feature processing are combined, and dimension reduction features serve as seismograph performance evaluation values; and training the processed feature sequence by using an attention-based space-time diagram convolutional network, capturing the correlation between time and space, evaluating the performance states of the plurality of seismographs in real time, calculating a performance anomaly index, and performing accurate performance evaluation. Experimental results show that the limitation that real-time monitoring cannot be achieved through a traditional seismograph performance evaluation method is solved, and the accuracy and the real-time performance of broadband seismograph performance evaluation are improved.

Description

Technical Field

[0001] The present invention relates to the field of seismograph performance evaluation, and in particular to an online performance evaluation method for a broadband seismograph based on spatiotemporal characteristics of multiple stations. Background Art

[0002] Broadband seismographs are high-precision instruments capable of capturing broadband dynamic variations in seismic signals. They are widely used in seismology, geophysics, and other related fields to monitor and record data such as seismic waveforms, crustal movement, and seismic background noise. To ensure the accuracy of broadband seismographs in recording seismic background noise, real-time monitoring is required during operation. If anomalies, such as amplitude or phase anomalies, occur, it is necessary to determine whether the abnormality is due to the performance of the broadband seismograph or to external natural factors. If a broadband seismograph does exhibit abnormal performance, prompt calibration and maintenance are required to ensure data accuracy.

[0003] However, currently, the most widely used methods for evaluating seismograph performance rely on on-site evaluation and assessment of seismic waveform data. This on-site evaluation method exhibits a certain degree of latency; when a seismograph's performance anomaly occurs, it can take days or even weeks to detect the problem, making it impossible to monitor instrument performance in real time. Furthermore, this method analyzes data from a single station, ignoring the spatiotemporal correlations between multiple stations. For example, data quality at a particular station may degrade due to instrument aging, environmental interference, or other factors. However, in a single-station evaluation, it is impossible to compare the performance of the data with that of surrounding stations. This single-perspective analysis can mask instrument performance anomalies. Seismic waveform data evaluation primarily focuses on surface and body waves. Surface-wave-based evaluation requires calculating the travel time of the background noise cross-correlation function over a long period of time, which is time-consuming. Body-wave-based evaluation relies on the occurrence of seismic events and lacks the persistence of data. Therefore, to address the shortcomings of existing seismograph performance evaluation methods, a method for online performance monitoring that integrates data from multiple stations is needed to promptly detect changes in seismograph performance, ensure the accuracy of earthquake observations, and mitigate the risks associated with instrument performance degradation. Summary of the Invention

[0004] In view of the shortcomings of the background technology, the purpose of the present invention is to propose an online evaluation method for broadband seismograph performance based on the spatiotemporal characteristics of multiple stations, so as to solve the limitations of current seismograph performance evaluation methods that cannot be monitored in real time and the performance of a single station, thereby improving the accuracy of broadband seismograph performance evaluation.

[0005] The technical solution adopted by the present invention includes the following steps:

[0006] S1, obtain seismic background noise data from multiple stations and perform preprocessing;

[0007] S2, feature extraction of preprocessed data;

[0008] S3, construction of spatiotemporal graph convolutional network based on attention mechanism;

[0009] S4, training, validation and test data set division, model training and validation;

[0010] S5, model testing, evaluates seismograph performance results.

[0011] The specific step S1 is:

[0012] S2.1, perform ADF test on the background noise data collected by multiple broadband seismometers to eliminate unstable sequences.

[0013] S2.2, the stationary data are subjected to band-pass filtering to remove the mean, trend and period of 10 to 50 s to eliminate useless noise interference in the seismic background noise.

[0014] S2.3, uses sliding windows to divide the preprocessed data into multiple data streams to achieve data enhancement.

[0015] The specific steps of step S2 are:

[0016] S3.1, continuous wavelet transform (CWT), Hilbert-Huang transform (HHT) and third-order statistics (TOC) are used to extract the time-frequency and statistical features from the seismic background noise.

[0017] CWT is used to extract features from seismic background noise data, which is determined according to the following formula:

[0018]

[0019] Where u is the translation factor that controls the position of the wavelet window in the time domain; v is the proportional factor that controls the size of the wavelet window and its position in the frequency domain; is the wavelet basis function.

[0020] HHT is used to extract features from seismic background noise data, which is determined according to the following formula:

[0021]

[0022] Where τ is the integral variable; t is time; by constructing the analytical signal At this time, the seismic background noise data is determined according to the following formula:

[0023]

[0024] Where ai (t) is the instantaneous amplitude, ω i (t) is the instantaneous frequency,

[0025] TOC is used to extract features from seismic background noise data, which is determined according to the following formula:

[0026] C 3x (τ1,τ2)=E[x(n)x(n+τ1)x(n+τ2)]

[0027] Where E is the mathematical expectation; τ1 and τ2 are time delays.

[0028] S3.2, information entropy is used to quantitatively describe the key information of seismic background noise data after CWT, HHT and TOC feature extraction. For a set of probabilities x i The information entropy of a discrete random variable X is determined according to the following formula:

[0029]

[0030] In the formula, p(x i ) is every value x i Probability of occurrence.

[0031] S3.3, perform principal component analysis PCA dimensionality reduction on the calculated feature information entropy. For an n-dimensional feature vector X, the data matrix is Where m is the number of samples, n is the dimension of the information entropy feature, and each feature in X is standardized according to the following formula:

[0032] Y=X-μ

[0033] The covariance matrix C is then calculated using the standardized data Y, which is determined according to the following formula:

[0034]

[0035] For the covariance matrix C, the eigenvalues and eigenvectors are determined according to the following formula:

[0036] Cv i =λ i v i

[0037] Where λ is the eigenvalue; v is the corresponding eigenvector; arrange the eigenvalues of the covariance matrix and select the largest eigenvalue λ1; the eigenvector v1 corresponding to this eigenvalue is the first principal component; project the data onto the direction of the first principal component to obtain data reduced to 1 dimension, which is determined according to the following formula:

[0038] XPCA =Yv1

[0039] Where, X PCA It is the feature after being reduced to one dimension, X PCA Normalize to [0,1], set X PCA is the performance evaluation value δ, and the instrument performance is evaluated according to the change of δ at the station.

[0040] The specific steps of step S3 are:

[0041] Furthermore, the structure of the attention-based spatiotemporal graph convolutional network (ASTGCN) in step S3 is introduced in detail:

[0042] S4.1, construct the input layer. The input layer includes a data input layer, which inputs the data set into the network. The input layer receives a data set with a shape of N×C×T h Input data, where N represents the number of stations; C represents the performance index of each station, which is the performance evaluation value here; T h Represents the historical time step, the input data contains the past T h The input data is then processed by the position embedding layer, which includes temporal position embedding and spatial position embedding.

[0043] S4.2, build the encoder. The network contains multiple encoder layers. Each encoder layer consists of a multi-head self-attention module, a graph convolutional network module, and a normalization and residual connection module:

[0044] A. Multi-head self-attention module, which replaces some projection operations with 1D convolution operations to enable it to perceive local changes in the data. Each head processes the input data in parallel and calculates the attention weights in different subspaces. The results are finally spliced and the output is generated through linear projection. The multi-station performance evaluation value δ is input into the encoder to extract temporal feature information. The long-term temporal dependencies of the data are learned through the attention mechanism to capture the temporal relationship in the features. For a single attention mechanism matrix, it is determined according to the following formula:

[0045]

[0046] Where Q, K, and V represent matrices consisting of multiple query, key, and value vectors. Each query, key, and value vector is represented by a d-dimensional vector. In the attention spatiotemporal data, a multi-head attention mechanism is used to extract temporal features, which are determined according to the following formula:

[0047] MultiHead(Q,K,V)=Concat(head1,...,head h )W O

[0048] head i =Attention(QW i Q ,KW i K ,VW i V )

[0049] Where h is the number of heads, W is i Q 、W i K and W i V are the projection matrices for Q, K, and V.

[0050] B. Graph Convolutional Network (GCN) module, which uses the self-attention mechanism to calculate the spatial correlation between nodes. For each time step, the GCN module calculates the spatial correlation of the adjacency matrix and updates the node features through graph convolution to capture the spatial correlation between stations. In GCN, two spatial data structures, nodes and edges, are used to represent spatial information. For the broadband seismograph in the present invention, the station is represented as a node structure in the graph, and the distance relationship between the stations is represented as an edge structure in the graph. GCN is used to learn the spatial correlation information between multiple stations, and the spatial distribution characteristics of the seismograph are encoded into the model together with the spatial characteristics of the neighboring stations. Each layer of GCN is specifically determined according to the following formula:

[0051]

[0052] Where H (k) represents the output of the kth layer; H (0) Indicates the input value N is the number of stations; F represents the feature dimension of each station; σ represents the nonlinear activation function, and ReLU is used as the activation function; It is represented as an adjacency matrix between multiple stations; is the degree matrix; W (k) is the weight matrix of the kth layer, where each element in the adjacency matrix is determined according to the following formula:

[0053]

[0054] Where α represents a non-zero parameter, set to 0.003; d represents the distance between stations, which is set to 0 for distances greater than 150 km and is constructed according to the formula for distances less than 150 km. The distance between stations is determined according to the Haversine distance formula, specifically according to the following formula:

[0055]

[0056] Where r is the radius of the Earth; φ and λ are the latitude and longitude of the two stations.

[0057] The azimuth between stations is determined according to the following formula:

[0058]

[0059] C. Layer Normalization and Residual Connection Module,Layer normalization and residual connection are used in the encoder to make the model,training more stable and avoid the gradient vanishing problem in,deep networks.

[0060] S4.3, build the decoder. The decoder structure consists of four modules: masked multi-head self-attention module, multi-head self-attention module, graph convolution module, and layer normalization and residual connection module:

[0061] A. Masked Multi-Head Self-Attention Module. Unlike the multi-head self-attention module in the encoder, the masked multi-head self-attention module in the decoder uses causal convolution instead of 1D convolution. Through the causal convolution operation, the decoder can only extract past input information and cannot obtain information about subsequent input sequences, thus ensuring the rationality of the prediction process.

[0062] B. Multi-head self-attention module. This module uses the encoder input and the masked multi-head self-attention module input to capture temporal relationships. Each head processes the input data in parallel and calculates attention weights in different subspaces. The results are finally concatenated and linearly projected to generate the output.

[0063] C. Graph Convolution Module,The decoder also contains a graph convolution module, which is used to adjust the spatial correlation,between the node features generated by the decoder.

[0064] D. Layer Normalization and Residual Connection Module,Layer normalization and residual connection are used in the decoder to avoid the problem of gradient disappearance in the network.

[0065] S4.4, construct the output layer, and the final output maps the high-dimensional features to the target output space through linear projection to generate the current performance evaluation value matrix The shape of the output matrix is N×C×T p , that is, evaluate the current T p The performance evaluation value of each time step.

[0066] The specific step S4 is:

[0067] S5.1, input the training set and validation set into ASTGCN for training, use the performance evaluation value of the previous n segments of each station to evaluate the performance evaluation value of the current segment, and realize real-time monitoring and status evaluation of the seismometer. According to the determination coefficient R 2 The performance of the model is evaluated by the mean absolute error (MAE), which is determined by the following formula:

[0068]

[0069] Where R 2 The interpretability of the model error is measured by calculating the complementary number of the ratio between the sum of squares of the errors between the model estimate and the true value and the variance of the true value. 2 The value of is in the interval [0,1]. The larger the value, the better the performance. MAE evaluates the performance of the model by calculating the average of the absolute errors between the model's estimated value and the true value. The smaller the value, the smaller the deviation between the model's estimated result and the true value, that is, the better the model performance.

[0070] S5.2. Use the data in the test set to evaluate the performance of the model and calculate the R of the test set data. 2 and MAE. Then calculate the abnormal deviation of the abnormal station and calculate the mean δ of the performance evaluation value of the abnormal station in the window f , and calculate the normal station performance evaluation value δ h and calculate the mean It is determined according to the following formula:

[0071]

[0072] Where m is the number of normal stations; δ h,i is the performance evaluation value of the i-th normal station in the time window, and the mean of the performance evaluation value of the abnormal station and the performance evaluation value of the normal station is calculated. The absolute difference between the two is used to obtain the performance deviation value β of the abnormal station, which is determined according to the following formula:

[0073]

[0074] The performance deviation β ranges from [0, +∞]. The abnormal threshold coefficient of β is set to 0.5. When β>0.5, performance abnormality occurs.

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

[0076] (1) The present invention combines wavelet transform, Hilbert-Huang transform and third-order cumulant feature extraction method to comprehensively analyze seismic background noise data from multiple features to improve the accuracy and robustness of seismograph performance evaluation.

[0077] (2) The present invention jointly evaluates the seismic background noise data of multiple stations through a spatiotemporal graph convolutional network based on the attention mechanism, effectively eliminating the influence of local environmental interference on the performance data of a single station, thereby ensuring the stability and reliability of the performance evaluation results.

[0078] (3) The present invention calculates the performance deviation value of the station, realizes the accurate quantification of the performance of the seismograph, and provides a scientific basis for the real-time monitoring and maintenance of the seismograph in practical applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0079] In order to make the purpose, technical solutions and advantages of the invention more clear, the present invention will be further described in detail below with reference to the accompanying drawings, in which:

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

[0081] Figure 2 A station diagram related to an embodiment of the present invention;

[0082] Figure 3 This is a signal diagram of the seismic background noise in the embodiment after being tested for stationarity, removed from the mean, trended, and filtered;

[0083] Figure 4 This is a feature extraction diagram of earthquake background noise after processing in the embodiment;

[0084] Figure 5 This is the ASTGCN architecture diagram;

[0085] Figure 6 Loss curve diagram for training set and validation set;

[0086] Figure 7 This is a regression diagram of the performance evaluation values of 9 stations under normal conditions in the embodiment;

[0087] Figure 8 This is a regression diagram of the performance evaluation values of the nine stations under the abnormal condition of the O17K station in the embodiment;

[0088] Figure 9 : is a performance deviation curve diagram of the O17K station in the abnormal time window in the embodiment;

[0089] Figure 10 1 is a graph showing the cross-correlation function of one day's data and two waveforms after removing different instrument responses at the O17K station in the embodiment. DETAILED DESCRIPTION

[0090] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0091] like Figure 1 As shown, the present invention discloses an online evaluation method for broadband seismograph performance based on the spatiotemporal characteristics of multiple stations, comprising the following steps:

[0092] S1, obtain seismic background noise data from multiple stations and perform preprocessing;

[0093] S2, feature extraction of preprocessed data;

[0094] S3, construction of spatiotemporal graph convolutional network based on attention mechanism;

[0095] S4, training, validation and test data set division, model training and validation;

[0096] S5, model testing, evaluates seismograph performance results.

[0097] The specific step S1 is:

[0098] S2.1, perform ADF test on the background noise data collected by multiple broadband seismometers to eliminate unstable sequences.

[0099] S2.2, the stationary data are subjected to band-pass filtering to remove the mean, trend and period of 10 to 50 s to eliminate useless noise interference in the seismic background noise.

[0100] S2.3, uses sliding windows to divide the preprocessed data into multiple data streams to achieve data enhancement.

[0101] The specific steps of step S2 are:

[0102] S3.1, continuous wavelet transform (CWT), Hilbert-Huang transform (HHT) and third-order statistics (TOC) are used to extract the time-frequency and statistical features from the seismic background noise.

[0103] CWT is used to extract features from seismic background noise data, which is determined according to the following formula:

[0104]

[0105] Where u is the translation factor that controls the position of the wavelet window in the time domain; v is the proportional factor that controls the size of the wavelet window and its position in the frequency domain; is the wavelet basis function.

[0106] HHT is used to extract features from seismic background noise data, which is determined according to the following formula:

[0107]

[0108] Where τ is the integral variable; t is time; by constructing the analytical signal At this time, the seismic background noise data is determined according to the following formula:

[0109]

[0110] Where a i (t) is the instantaneous amplitude, ω i (t) is the instantaneous frequency,

[0111] TOC is used to extract features from seismic background noise data, which is determined according to the following formula:

[0112] C 3x (τ1,τ2)=E[x(n)x(n+τ1)x(n+τ2)]

[0113] Where E is the mathematical expectation; τ1 and τ2 are time delays.

[0114] S3.2, information entropy is used to quantitatively describe the key information of seismic background noise data after CWT, HHT and TOC feature extraction. For a set of probabilities x i The information entropy of a discrete random variable X is determined according to the following formula:

[0115]

[0116] In the formula, p(x i ) is every value x i Probability of occurrence.

[0117] S3.3, perform principal component analysis PCA dimensionality reduction on the calculated feature information entropy. For an n-dimensional feature vector X, the data matrix is Where m is the number of samples, n is the dimension of the information entropy feature, and each feature in X is standardized according to the following formula:

[0118] Y=X-μ

[0119] The covariance matrix C is then calculated using the standardized data Y, which is determined according to the following formula:

[0120]

[0121] For the covariance matrix C, the eigenvalues and eigenvectors are determined according to the following formula:

[0122] Cv i =λi v i

[0123] Where λ is the eigenvalue; v is the corresponding eigenvector; arrange the eigenvalues of the covariance matrix and select the largest eigenvalue λ1; the eigenvector v1 corresponding to this eigenvalue is the first principal component; project the data onto the direction of the first principal component to obtain data reduced to 1 dimension, which is determined according to the following formula:

[0124] X PCA =Yv1

[0125] Where, X PCA It is the feature after being reduced to one dimension, X PCA Normalize to [0,1], set X PCA is the performance evaluation value δ, and the instrument performance is evaluated according to the change of δ at the station.

[0126] The specific steps of step S3 are:

[0127] Furthermore, the structure of the attention-based spatiotemporal graph convolutional network (ASTGCN) in step S3 is introduced in detail:

[0128] S4.1, construct the input layer. The input layer includes a data input layer, which inputs the data set into the network. The input layer receives a data set with a shape of N×C×T h Input data, where N represents the number of stations; C represents the performance index of each station, which is the performance evaluation value here; T h Represents the historical time step, the input data contains the past T h The input data is then processed by the position embedding layer, which includes temporal position embedding and spatial position embedding.

[0129] S4.2, build the encoder. The network contains multiple encoder layers. Each encoder layer consists of a multi-head self-attention module, a graph convolutional network module, and a normalization and residual connection module:

[0130] A. Multi-head self-attention module, which replaces some projection operations with 1D convolution operations to enable it to perceive local changes in the data. Each head processes the input data in parallel and calculates the attention weights in different subspaces. The results are finally spliced and the output is generated through linear projection. The multi-station performance evaluation value δ is input into the encoder to extract temporal feature information. The long-term temporal dependencies of the data are learned through the attention mechanism to capture the temporal relationship in the features. For a single attention mechanism matrix, it is determined according to the following formula:

[0131]

[0132] Where Q, K, and V represent matrices consisting of multiple query, key, and value vectors. Each query, key, and value vector is represented by a d-dimensional vector. In the attention spatiotemporal data, a multi-head attention mechanism is used to extract temporal features, which are determined according to the following formula:

[0133] MultiHead(Q,K,V)=Concat(head1,...,head h )W O

[0134] head i =Attention(QW i Q ,KW i K ,VW i V )

[0135] Where h is the number of heads, W is i Q 、W i K and W i V are the projection matrices for Q, K, and V.

[0136] B. Graph Convolutional Network (GCN) module, which uses the self-attention mechanism to calculate the spatial correlation between nodes. For each time step, the GCN module calculates the spatial correlation of the adjacency matrix and updates the node features through graph convolution to capture the spatial correlation between stations. In GCN, two spatial data structures, nodes and edges, are used to represent spatial information. For the broadband seismograph in the present invention, the station is represented as a node structure in the graph, and the distance relationship between the stations is represented as an edge structure in the graph. GCN is used to learn the spatial correlation information between multiple stations, and the spatial distribution characteristics of the seismograph are encoded into the model together with the spatial characteristics of the neighboring stations. Each layer of GCN is specifically determined according to the following formula:

[0137]

[0138] Where H (k) represents the output of the kth layer; H (0) Indicates the input value N is the number of stations; F represents the feature dimension of each station; σ represents the nonlinear activation function, and ReLU is used as the activation function; It is represented as an adjacency matrix between multiple stations; is the degree matrix; W (k) is the weight matrix of the kth layer, where each element in the adjacency matrix is determined according to the following formula:

[0139]

[0140] Where α represents a non-zero parameter, set to 0.003; d represents the distance between stations, which is set to 0 for distances greater than 150 km and is constructed according to the formula for distances less than 150 km. The distance between stations is determined according to the Haversine distance formula, specifically according to the following formula:

[0141]

[0142] Where r is the radius of the Earth; φ and λ are the latitude and longitude of the two stations.

[0143] The azimuth between stations is determined according to the following formula:

[0144]

[0145] C. Layer Normalization and Residual Connection Module,Layer normalization and residual connection are used in the encoder to make the model,training more stable and avoid the gradient vanishing problem in,deep networks.

[0146] S4.3, build the decoder. The decoder structure consists of four modules: masked multi-head self-attention module, multi-head self-attention module, graph convolution module, and layer normalization and residual connection module:

[0147] A. Masked Multi-Head Self-Attention Module. Unlike the multi-head self-attention module in the encoder, the masked multi-head self-attention module in the decoder uses causal convolution instead of 1D convolution. Through the causal convolution operation, the decoder can only extract past input information and cannot obtain information about subsequent input sequences, thus ensuring the rationality of the prediction process.

[0148] B. Multi-head self-attention module. This module uses the encoder input and the masked multi-head self-attention module input to capture temporal relationships. Each head processes the input data in parallel and calculates attention weights in different subspaces. The results are finally concatenated and linearly projected to generate the output.

[0149] C. Graph Convolution Module,The decoder also contains a graph convolution module, which is used to adjust the spatial correlation,between the node features generated by the decoder.

[0150] D. Layer Normalization and Residual Connection Module,Layer normalization and residual connection are used in the decoder to avoid the problem of gradient disappearance in the network.

[0151] S4.4, construct the output layer, and the final output maps the high-dimensional features to the target output space through linear projection to generate the current performance evaluation value matrix The shape of the output matrix is N×C×T p , that is, evaluate the current T p The performance evaluation value of each time step.

[0152] The specific step S4 is:

[0153] S5.1, input the training set and validation set into ASTGCN for training, use the performance evaluation value of the previous n segments of each station to evaluate the performance evaluation value of the current segment, and realize real-time monitoring and status evaluation of the seismometer. According to the determination coefficient R 2 The performance of the model is evaluated by the mean absolute error (MAE), which is determined by the following formula:

[0154]

[0155] Where R 2 The interpretability of the model error is measured by calculating the complementary number of the ratio between the sum of squares of the errors between the model estimate and the true value and the variance of the true value. 2 The value of is in the interval [0,1]. The larger the value, the better the performance. MAE evaluates the performance of the model by calculating the average of the absolute errors between the model's estimated value and the true value. The smaller the value, the smaller the deviation between the model's estimated result and the true value, that is, the better the model performance.

[0156] S5.2. Use the data in the test set to evaluate the performance of the model and calculate the R of the test set data. 2 and MAE. Then calculate the abnormal deviation of the abnormal station and calculate the mean δ of the performance evaluation value of the abnormal station in the window f , and calculate the normal station performance evaluation value δ h and calculate the mean It is determined according to the following formula:

[0157]

[0158] Where m is the number of normal stations; δ h,i is the performance evaluation value of the i-th normal station in the time window, and the mean of the performance evaluation value of the abnormal station and the performance evaluation value of the normal station is calculated. The absolute difference between the two is used to obtain the performance deviation value β of the abnormal station, which is determined according to the following formula:

[0159]

[0160] The performance deviation β ranges from [0, +∞]. The abnormal threshold coefficient of β is set to 0.5. When β>0.5, performance abnormality occurs.

[0161] The following is a more detailed description of the online performance evaluation method of a broadband seismograph based on spatiotemporal characteristics of multiple stations based on embodiments of the present invention.

[0162] The specific process of earthquake background noise preprocessing is as follows: First, collect Figure 2 Seismic background noise data collected by 9 stations in the China TA network. Non-stationary signals were removed by ADF test, and the mean and trend were removed. Bandpass filtering was then performed for 10 to 50 seconds. The data was then divided into multiple data streams using a sliding window. The sliding window parameters were set, with a window width of two hours and a 50% overlap between windows. Figure 3 This is a data stream of the O17K station in the embodiment. Each station generates a total of 7653 data streams, and the total number of samples is 9×7653.

[0163] The specific process of seismic background noise feature extraction is as follows: First, feature extraction is performed on the data stream, and feature information is extracted using CWT, HHT, and TOC. Then, information entropy processing is performed on the extracted features, where the information entropy of wavelet coefficients greater than 0.0015 is extracted as feature 1; the information entropy of Hilbert amplitude greater than 0.0082 is extracted as feature 2; the information entropy of TOC greater than 0 is extracted as feature 3, and the information entropy of TOC less than 0 is extracted as feature 4. The feature extraction process is as follows: Figure 4 Finally, PCA is used to reduce the dimensionality of these four features to one dimension, which is used as the input features for evaluating the performance of the instrument.

[0164] The specific process of ASTGCN training is as follows: First, the dataset is divided into training set, validation set and test set in a ratio of 8:1:1. The architecture diagram of ASTGCN is as follows Figure 5 As shown in Figure 2. ASTGCN is used to train the training set and validation set. The window size is set to 30 time periods, and the performance evaluation value of each window is gradually evaluated. The evaluation length is 60 time periods. The number of iterations is set to 300 rounds, the learning rate is 0.0001, and the batch size of each round is 32. The loss curves of the training set and validation set are shown in Figure 2. Figure 6 As shown, Figure 6 The training set in (a) reaches convergence at the 50th round. Figure 6 The validation set in (b) reaches convergence at round 150. The performance index results of the proposed method on the validation data are shown in Table 1.

[0165] Table 1 Performance index results of the method proposed in the present invention

[0166]

[0167] Table 1 shows the R of the model on the validation data. 2 The score result value is 0.996, which is close to 1, indicating that the model can effectively perform regression fitting; the MAE score result is 0.0003, indicating that the mean absolute error of the model is 0.0003, which meets the accuracy requirements of the model.

[0168] The ASTGCN station performance evaluation process involves conducting an online performance evaluation of nine stations in the test set using ASTGCN. The model uses the station's historical performance evaluation values from the previous 30 periods to assess the performance of the current period. Subsequently, as the time window continuously slides, the model gradually uses the historical performance evaluation values from the previous 30 periods to assess the performance of the current period, thus enabling real-time online evaluation of seismometer performance. Figure 7 The evaluation results of 60 performance evaluation values of 9 stations are presented.

[0169] In the 30th time period, the instrument response sensitivity and extreme point of the O17K station were modified. The sensitivity was modified from 5.09222E+8 to 6.26628E+8, and the extreme point was modified from -0.148600+0.148600i to -0.113100+0.00000i. In order to specifically measure the time deviation after the change in the instrument performance of the O17K station, a cross-correlation operation was performed on the background noise data of the normal day of O17K and the background noise data of the abnormal day. Figure 8 As shown in (c), the maximum value of the cross-correlation function of the two waveforms deviates from one sample point at time zero. Since the sampling frequency of the data is 40Hz, this indicates that the data of the O17K station has a time deviation of 0.025s. Figure 9 In the data, starting from the 30th time period, the regression effect of the O17K station deteriorated significantly, and the model could not be effectively fitted, indicating that the seismic instrument had a performance anomaly. At this time, the performance evaluation value of the O17K station gradually increased from the 30th time period, reaching a peak of about 0.87%. Then, as the window slid, the abnormal data tended to stabilize, and the performance evaluation value began to slowly decline to a steady state; while the difference between the true value and the evaluation value of the performance evaluation value of other stations was small. According to the performance deviation formula, the sliding window size is set to 5, the window sliding step is 1, a total of 26 window values are calculated, and the performance deviation of O17K is calculated within these windows. Figure 10 As shown in the figure, the performance deviations of the O17K station are all greater than 0.5, and the trend of the performance deviations is consistent with the trend of the performance evaluation value. This indicates that performance anomalies occurred in the 30th time period and timely calibration is required.

[0170] In summary, this paper proposes an online performance evaluation method for broadband seismographs based on the spatiotemporal characteristics of multiple stations. This method can comprehensively assess the performance of multiple stations and promptly detect the 0.025s time deviation generated by the station data in the embodiment, effectively improving the accuracy and real-time performance evaluation.

Claims

1. A method for online performance evaluation of broadband seismographs based on spatiotemporal characteristics of multiple stations, characterized by: The method comprises the following steps: S1, obtain seismic background noise data from multiple stations and perform preprocessing; S2, feature extraction of preprocessed data; S3, construction of spatiotemporal graph convolutional network based on attention mechanism; S4, training, validation and test data set division, model training and validation; S5, model testing, evaluates seismograph performance results. The specific step S1 is: S2.1, perform ADF test on the background noise data collected by multiple broadband seismometers to eliminate unstable sequences. S2.2, the stationary data are subjected to band-pass filtering to remove the mean, trend and period of 10 to 50 s to eliminate useless noise interference in the seismic background noise. S2.3, uses sliding windows to divide the preprocessed data into multiple data streams to achieve data enhancement. The specific step S2 is: S3.1, continuous wavelet transform (CWT), Hilbert-Huang transform (HHT) and third-order statistics (TOC) are used to extract the time-frequency and statistical features from the seismic background noise. CWT is used to extract features from seismic background noise data, which is determined according to the following formula: Where u is the translation factor that controls the position of the wavelet window in the time domain; v is the proportional factor that controls the size of the wavelet window and its position in the frequency domain; is the wavelet basis function. HHT is used to extract features from seismic background noise data, which is determined according to the following formula: Where τ is the integral variable; t is time; by constructing the analytical signal At this time, the seismic background noise data is determined according to the following formula: Where a i (t) is the instantaneous amplitude, ω i (t) is the instantaneous frequency, TOC is used to extract features from seismic background noise data, which is determined according to the following formula: C 3x (τ1,τ2)=E[x(n)x(n+τ1)x(n+τ2)] Where E is the mathematical expectation; τ1 and τ2 are time delays. S3.2, information entropy is used to quantitatively describe the key information of seismic background noise data after CWT, HHT and TOC feature extraction. For a set of probabilities x i The information entropy of a discrete random variable X is determined according to the following formula: In the formula, p(x i ) is every value x i Probability of occurrence. S3.3, perform principal component analysis PCA dimensionality reduction on the calculated feature information entropy. For an n-dimensional feature vector X, the data matrix is Where m is the number of samples, n is the dimension of the information entropy feature, and each feature in X is standardized according to the following formula: Y=X-μ The covariance matrix C is then calculated using the standardized data Y, which is determined according to the following formula: For the covariance matrix C, the eigenvalues and eigenvectors are determined according to the following formula: Cv i =λ i v i Where λ is the eigenvalue; v is the corresponding eigenvector; arrange the eigenvalues of the covariance matrix and select the largest eigenvalue λ1; the eigenvector v1 corresponding to this eigenvalue is the first principal component; project the data onto the direction of the first principal component to obtain data reduced to 1 dimension, which is determined according to the following formula: X PCA =Yv1 Where, X PCA It is the feature after reducing to one dimension, X PCA Normalize to [0,1], set X PCA is the performance evaluation value δ, and the instrument performance is evaluated according to the change of δ at the station. The specific steps of step S3 are: Furthermore, the structure of the attention-based spatiotemporal graph convolutional network (ASTGCN) in step S3 is introduced in detail: S4.1, construct the input layer. The input layer includes a data input layer, which inputs the data set into the network. The input layer receives a data set with a shape of N×C×T h Input data of , where N represents the number of stations; C represents the performance index of each station, which is the performance evaluation value here; T h Represents the historical time step, the input data contains the past T h The input data is then processed by the position embedding layer, which includes temporal position embedding and spatial position embedding. S4.2, build the encoder. The network contains multiple encoder layers. Each encoder layer consists of a multi-head self-attention module, a graph convolutional network module, and a normalization and residual connection module: A. Multi-head self-attention module. This module replaces some projection operations with 1D convolution operations, enabling it to perceive local changes in the data. Each head processes the input data in parallel and calculates the attention weights in different subspaces. The results are finally concatenated and the output is generated through linear projection. The multi-station performance evaluation value δ is input into the encoder to extract temporal feature information. The attention mechanism learns the long-term temporal dependencies of the data and captures the temporal relationships in the features. For a single attention mechanism matrix, it is determined according to the following formula: Where Q, K, and V represent matrices consisting of multiple query, key, and value vectors. Each query, key, and value vector is represented by a d-dimensional vector. In the attention spatiotemporal data, a multi-head attention mechanism is used to extract temporal features, which is determined by the following formula: MultiHead(Q,K,V)=Concat(head1,...,head h )W O head i =Attention(QW i Q ,KW i K ,VW i V ) Where h is the number of heads, W is i Q 、W i K and W i V are the projection matrices for Q, K, and V. B. Graph Convolutional Network (GCN) module, which uses the self-attention mechanism to calculate the spatial correlation between nodes. For each time step, the GCN module calculates the spatial correlation of the adjacency matrix and updates the node features through graph convolution to capture the spatial correlation between stations. In GCN, two spatial data structures, nodes and edges, are used to represent spatial information. For the broadband seismograph in the present invention, the station is represented as a node structure in the graph, and the distance relationship between the stations is represented as an edge structure in the graph. GCN is used to learn the spatial correlation information between multiple stations, and the spatial distribution characteristics of the seismograph are encoded into the model together with the spatial characteristics of the neighboring stations. Each layer of GCN is specifically determined according to the following formula: Where H (k) represents the output of the kth layer; H (0) Indicates the input value N is the number of stations; F represents the feature dimension of each station; σ represents the nonlinear activation function, and ReLU is used as the activation function; It is represented as an adjacency matrix between multiple stations; is the degree matrix; W (k) is the weight matrix of the kth layer, where each element in the adjacency matrix is determined according to the following formula: Where α represents a non-zero parameter, set to 0.003; d represents the distance between stations, which is set to 0 for distances greater than 150 km and is constructed according to the formula for distances less than 150 km. The distance between stations is determined according to the Haversine distance formula, specifically according to the following formula: Where r is the radius of the Earth; φ and λ are the latitude and longitude of the two stations. The azimuth between stations is determined according to the following formula: C. Layer Normalization and Residual Connection Module,Layer normalization and residual connection are used in the encoder to make the model,training more stable and avoid the gradient vanishing problem in,deep networks. S4.3, build the decoder. The decoder structure consists of four modules: masked multi-head self-attention module, multi-head self-attention module, graph convolution module, and layer normalization and residual connection module: A. Masked Multi-Head Self-Attention Module. Unlike the multi-head self-attention module in the encoder, the masked multi-head self-attention module in the decoder uses causal convolution instead of 1D convolution. Through the causal convolution operation, the decoder can only extract past input information and cannot obtain information about subsequent input sequences, thus ensuring the rationality of the prediction process. B. Multi-head self-attention module. This module uses the encoder input and the masked multi-head self-attention module input to capture temporal relationships. Each head processes the input data in parallel and calculates attention weights in different subspaces. The results are finally concatenated and linearly projected to generate the output. C. Graph Convolution Module,The decoder also contains a graph convolution module, which is used to adjust the spatial correlation,between the node features generated by the decoder. D. Layer Normalization and Residual Connection Module,Layer normalization and residual connection are used in the decoder to avoid the problem of gradient disappearance in the network. S4.4, construct the output layer, and the final output maps the high-dimensional features to the target output space through linear projection to generate the current performance evaluation value matrix The shape of the output matrix is N×C×T p , that is, evaluate the current T p The performance evaluation value of each time step. The specific step S4 is: S5.1, input the training set and validation set into ASTGCN for training, use the performance evaluation value of the previous n segments of each station to evaluate the performance evaluation value of the current segment, and realize real-time monitoring and status evaluation of the seismometer. According to the determination coefficient R 2 The performance of the model is evaluated by the mean absolute error (MAE), which is determined by the following formula: Where R 2 The interpretability of the model error is measured by calculating the complementary number of the ratio between the sum of squares of the errors between the model estimate and the true value and the variance of the true value. 2 The value of is in the interval [0,1]. The larger the value, the better the performance. MAE evaluates the performance of the model by calculating the average of the absolute errors between the model's estimated value and the true value. The smaller the value, the smaller the deviation between the model's estimated result and the true value, that is, the better the model performance. S5.

2. Use the data in the test set to evaluate the performance of the model and calculate the R of the test set data. 2 and MAE. Then calculate the abnormal deviation of the abnormal station and calculate the mean δ of the performance evaluation value of the abnormal station in the window f , and calculate the normal station performance evaluation value δ h and calculate the mean It is determined according to the following formula: Where m is the number of normal stations; δ h,i is the performance evaluation value of the i-th normal station in the time window, and the mean of the performance evaluation value of the abnormal station and the performance evaluation value of the normal station is calculated. The absolute difference between the two is used to obtain the performance deviation value β of the abnormal station, which is determined according to the following formula: The performance deviation β ranges from [0, +∞]. The abnormal threshold coefficient of β is set to 0.

5. When β>0.5, performance abnormality occurs.

Images