Method for identifying low-quality data of strong motion records

Abstract

The present invention relates to a method for identifying low-quality data of strong vibration records, and belongs to the technical field of earthquake monitoring data processing. The method solves the problems of low efficiency and high missed detection rate of manual identification. The key points of the technical solution include: first, obtaining the three-component acceleration time history data of the strong vibration record; then constructing a multi-dimensional dynamic quality assessment index system including data integrity, signal-to-noise ratio dynamic range, waveform anomaly index, frequency response offset and polarity consistency coefficient; based on the system, the data is scored for quality in real time, and when the score is lower than the dynamic threshold, it is marked as low-quality data; finally, an adaptive classification algorithm is used to identify the abnormal type of the marked data, including at least one of waveform distortion, baseline drift, polarity reversal and site response anomaly. The method is mainly used to improve the reliability and engineering application value of strong vibration observation data.

Description

Technical Field

[0001] The present invention relates to the field of seismic data processing, and more particularly to a method for identifying low-quality data of strong earthquake records. Background Art

[0002] Strong motion records play a vital role in earthquake engineering research and applications, but their data quality is susceptible to a variety of factors. Currently, the identification of low-quality records relies primarily on manual judgment, which is inefficient and highly subjective. Operators must individually examine factors such as waveform morphology and spectral characteristics, which is time-consuming and difficult to process with massive amounts of data. Manual identification accuracy is particularly low for records containing complex anomalies or hidden defects.

[0003] Existing technologies have limitations in identifying some low-quality types. For example, while waveform anomalies can be detected through visual observation, there is a lack of quantitative assessment standards, and different personnel may have significant differences in their judgment of the degree of anomaly. Existing methods often fail to effectively identify non-waveform anomalies, such as sensor polarity reversal or abnormal site response. Polarity reversal is often caused by incorrect sensor wiring, but traditional technical methods have difficulty automatically detecting spatial directional deviations. Site response anomalies are affected by local geological conditions, and there is no mature automated identification solution.

[0004] Some studies have attempted to use fixed thresholds for automated screening, but this approach lacks adaptability. Strong motion records are influenced by source characteristics, propagation paths, and site conditions, making a single threshold difficult to adapt to diverse scenarios. The signal-to-noise ratio distributions between near-field and far-field records differ significantly, making fixed thresholds prone to misjudgment. Furthermore, quality assessment involves multi-dimensional indicators, and the lack of a dynamic adjustment mechanism for weighting these indicators makes it difficult to accurately reflect data reliability in complex environments.

[0005] Another shortcoming of the existing recognition process is its weak correction and verification processes. After identifying low-quality data, most methods simply mark or eliminate it, lacking targeted remediation measures. For example, correcting baseline drift requires dynamic adjustment of filter parameters based on signal characteristics, but existing technologies lack a closed-loop optimization mechanism. Furthermore, the confidence assessment system for automated recognition results is still underdeveloped, making it difficult to determine when manual review is necessary, impacting overall processing efficiency.

[0006] These issues have resulted in a large number of strong motion records not being effectively utilized, limiting their application value in areas such as engineering seismic design and earthquake early warning. Developing more efficient and comprehensive methods for identifying low-quality data has become a key requirement for improving the availability of strong motion observation data. Summary of the Invention

[0007] An object of the present invention is to solve at least the above problems and to provide at least the advantages which will be described hereinafter.

[0008] In order to achieve these objects and other advantages according to the present invention, a method for identifying low-quality data of a vibration recording is provided, comprising the following steps:

[0009] Obtain three-component acceleration time history data of strong vibration records;

[0010] Construct a multi-dimensional dynamic quality assessment index system, including data integrity, signal-to-noise ratio dynamic range, waveform anomaly index, frequency response deviation and polarity consistency coefficient;

[0011] Performing real-time quality scoring on the acceleration time history data based on the indicator system, and marking it as low-quality data if the score is lower than a dynamic threshold;

[0012] An adaptive classification algorithm is used to identify abnormal types of the marked low-quality data, where the abnormal types include at least one of waveform distortion, baseline drift, polarity reversal and field response abnormality.

[0013] This invention establishes a systematic framework for identifying low-quality data. It uses a multidimensional dynamic quality assessment system to quantitatively score three-component acceleration time histories, overcoming the limitations of traditional reliance on manual judgment. A dynamic threshold mechanism adapts to different earthquake scenarios, significantly improving the accuracy of low-quality data labeling. An adaptive classification algorithm enables intelligent identification of anomaly types, providing a precise basis for subsequent data repair or removal, significantly improving the usability of strong motion records.

[0014] Preferably, the waveform anomaly index is calculated by:

[0015] Detect the waveform amplitude mutation rate. If the amplitude change rate of consecutive sampling points exceeds the preset mutation threshold, it is determined to be a mutation abnormality;

[0016] Analyze the smoothness of the waveform envelope. If the local curvature of the envelope exceeds the curvature threshold, it is determined to be a non-smooth anomaly.

[0017] Identify periodic noise interference and extract the noise energy ratio of the characteristic frequency through Fourier transform. If the ratio exceeds the noise threshold, it is determined to be periodic interference.

[0018] This paper proposes a three-dimensional waveform anomaly index calculation, effectively covering the main types of waveform distortion. Amplitude mutation rate detection captures transient interference, envelope curvature analysis identifies slowly varying distortion, and periodic noise energy ratio quantifies persistent interference. A weighted combination of these three factors forms a comprehensive assessment. This method can accurately locate anomaly sources such as instrument failures and environmental interference, reducing the waveform misjudgment rate by approximately 18%.

[0019] Preferably, the adaptive classification algorithm comprises:

[0020] Time domain analysis: Detects data integrity loss and baseline drift through a sliding window. Baseline drift is determined when the slope of the linear trend of the displacement time history exceeds the drift threshold.

[0021] Frequency domain analysis: Calculate the root mean square error between the frequency response and the standard station response curve. If the error exceeds the offset threshold, the frequency response is determined to be offset.

[0022] Polarization analysis: Calculates the azimuth deviation of the epicenter based on the eigenvectors of the three-component covariance matrix. If the deviation exceeds the azimuth tolerance, it is considered a polarity reversal.

[0023] Site response analysis: The H / V spectrum ratio coefficient of variation is used to identify site response anomalies. If the coefficient of variation exceeds the stability threshold, it is judged as abnormal.

[0024] This invention achieves fine-grained anomaly classification through four-dimensional parallel analysis. Time-domain sliding window detection detects data loss and baseline drift in real time. Frequency-domain root mean square error assesses device performance degradation. Polarization analysis ensures spatial orientation accuracy. H / V spectrum coefficient of variation monitors site response stability. Multi-dimensional cross-validation increases the accuracy of complex anomaly identification to 92%.

[0025] Preferably, the polarization analysis specifically includes:

[0026] Extract the three-component acceleration time history within the P wave initial motion window;

[0027] Calculate the eigenvalues ​​and principal eigenvector directions of the covariance matrix within the window;

[0028] Polarity reversal is determined based on the angle deviation between the main eigenvector and the measured azimuth of the station. If the angle deviation exceeds 30°, the correction process is triggered.

[0029] This invention standardizes the polarization analysis process, limiting the P-wave onset window to ensure optimal signal-to-noise ratio. Eigendecomposition of the covariance matrix accurately extracts the primary vibration direction. A 30° azimuth deviation threshold, verified through large-scale field measurements, balances sensitivity and false alarm rate. This solution achieves a success rate of over 95% for automatic polarity reversal correction, addressing inherent issues such as sensor wiring errors.

[0030] Preferably, it also includes:

[0031] Establish an indicator weight optimization model: Use a lightweight neural network to dynamically adjust the weight distribution of multi-dimensional indicators;

[0032] The neural network takes the waveform time domain characteristics and frequency domain characteristics as input, outputs the weight coefficients of each indicator, and optimizes the weight distribution strategy through back propagation.

[0033] This method employs a dynamic weighting mechanism, using a lightweight neural network to adaptively adjust indicator weights based on waveform characteristics. Near-field earthquakes enhance waveform anomaly indices, while far-field events prioritize signal-to-noise ratio assessment. Backpropagation optimization enables continuous evolution of the weighting strategy, resulting in a 23% improvement in quality score consistency with expert evaluation compared to a fixed-weight model.

[0034] Preferably, the neural network training step includes:

[0035] Construct a mixed sample set: including synthetic data simulating environmental noise and measured low-quality recorded data;

[0036] Feature enhancement: Expand sample diversity through time-frequency domain data enhancement technology;

[0037] Transfer learning: Adapt the pre-trained temporal feature extraction model to the weight distribution task, freezing the underlying convolutional layers and fine-tuning the fully connected layers.

[0038] This method uses synthetic data to simulate complex noisy environments and measured data to provide real-world examples. Time-frequency domain enhancement improves model generalization, and transfer learning reuses pretrained feature extractors to reduce data requirements. After 200 rounds of training, the weight allocation error is kept within 0.08, and the computational efficiency meets real-time processing requirements.

[0039] Preferably, a dynamic correction process is also included:

[0040] For labeled low-quality data, filter banks are selected based on the anomaly type;

[0041] Waveform distortion uses an adaptive bandpass filter, and the passband range is dynamically adjusted according to the signal's main frequency;

[0042] The baseline drift was corrected for the displacement time course using piecewise cubic spline fitting;

[0043] Polarity inversion automatically flips the direction of channel data using polarization analysis results.

[0044] This invention establishes a directional correction mechanism for anomaly types. An adaptive bandpass filter dynamically tracks the main frequency to eliminate waveform distortion; piecewise spline fitting accurately removes baseline drift trends; and a polarity flip function provides one-click correction for directional errors. This process has increased the availability of low-quality data from 54% to 89%, providing a reliable data foundation for earthquake engineering.

[0045] Preferably, the dynamic correction process further includes:

[0046] Iterative optimization mechanism: Taking the signal-to-noise ratio improvement rate as the objective function, the filter parameters are adjusted by the gradient descent method until the output signal meets the signal-to-noise ratio threshold or reaches the maximum number of iterations.

[0047] This invention incorporates an iterative optimization process, using the signal-to-noise ratio improvement rate as the objective function. Using a gradient descent method to dynamically adjust filter parameters such as the cutoff frequency, the signal-to-noise ratio can be improved by 133% (e.g., from 15dB to 35dB) after 3-5 iterations. A maximum number of iterations is limited to avoid infinite loops and ensure efficient processing.

[0048] Preferably, a confidence assessment is also included:

[0049] A confidence model is constructed based on the multi-dimensional indicator score deviation, historical identification accuracy and station reliability;

[0050] Output confidence score. If the score is lower than the reliability threshold, the manual review process is triggered.

[0051] This invention introduces a confidence quantification model that calculates a confidence score based on current indicator deviation, historical accuracy, and device status. This threshold trigger mechanism accurately selects cases requiring review. In practice, this reduces unnecessary manual intervention by 68%, focusing resources on true-positive anomaly verification.

[0052] Preferably, the manual review process includes:

[0053] Generate visual reports: mark the waveform abnormality location, algorithm judgment basis and correction suggestions;

[0054] An interactive interface is provided to support manual correction of labels and feedback to the training set for optimizing the recognition model.

[0055] This invention establishes a closed loop of human-machine collaboration, with visual reports that intuitively display the algorithm's decision logic; an interactive interface supports rapid label correction; and human feedback data optimizes the model in real time. A continuous learning mechanism has reduced the system's false alarm rate by 23.5% within six months, creating an intelligent evolutionary capability that becomes more accurate with use.

[0056] Other advantages, objectives and features of the present invention will be reflected in part from the following description and will be understood by those skilled in the art through study and practice of the present invention. DETAILED DESCRIPTION

[0057] The present invention is further described in detail below with reference to the embodiments so that those skilled in the art can implement the invention with reference to the description.

[0058] The present invention provides a method for identifying low-quality data of vibration recordings, comprising the following steps:

[0059] Obtain three-component acceleration time history data of strong vibration records;

[0060] Construct a multi-dimensional dynamic quality assessment index system, including data integrity, signal-to-noise ratio dynamic range, waveform anomaly index, frequency response deviation and polarity consistency coefficient;

[0061] Performing real-time quality scoring on the acceleration time history data based on the indicator system, and marking it as low-quality data if the score is lower than a dynamic threshold;

[0062] An adaptive classification algorithm is used to identify abnormal types of the marked low-quality data, where the abnormal types include at least one of waveform distortion, baseline drift, polarity reversal and field response abnormality.

[0063] Specifically, first, obtain three-component acceleration time history data for the strong motion event to be assessed. This data is typically derived from seismic network records and contains time-varying acceleration series data for two mutually perpendicular horizontal components (e.g., north-south and east-west) and one vertical component (UD). The sampling rate and range must comply with relevant observation specifications.

[0064] Next, a multidimensional dynamic quality assessment index system was constructed, which comprehensively examines multiple key quality dimensions of the data. The data integrity index is used to detect whether there are missing data points or invalid values ​​(such as NaN) in the record. The signal-to-noise ratio dynamic range index evaluates the energy ratio of the effective component of the signal (usually referring to ground motion) to the background noise and examines how this ratio changes over the entire recording period. For example, a dynamic range of at least 80 dB is required to ensure that weak seismic phases can be identified. The waveform anomaly index specifically captures waveform distortions. It is calculated by detecting abnormal amplitude mutations (such as the rate of change between consecutive sampling points exceeding a preset mutation threshold by 1.5 times), analyzing the smoothness of the waveform envelope (such as local curvature exceeding a curvature threshold of 0.05), and identifying periodic noise interference (such as detecting, through Fourier transform, that the noise energy proportion of a specific frequency exceeds a noise threshold of 5%). The frequency response deviation index measures the degree of difference between the frequency response characteristics of the measured record and the standard response curve calibrated by the station. The root mean square error (RMSE) between the two is usually calculated. The polarity consistency coefficient is based on the spatial relationship of the three-component data and evaluates whether the direction recorded by each component is consistent with the expected polarization characteristics of seismic wave propagation.

[0065] Based on the aforementioned indicator system, the input three-component acceleration time history data is scored in real time. Each indicator generates a sub-score based on its calculation rules and preset quality standards (e.g., data integrity must be 100% and the signal-to-noise ratio dynamic range must be >60dB). Each sub-score is then combined using weighting or normalization methods (e.g., normalization to a range of 0-1) to create an overall quality score. The system sets a dynamic threshold (which may be adaptively adjusted based on background noise levels or event scale). If the calculated overall quality score falls below this dynamic threshold, the entire strong motion record or the problematic data segment is marked as "low-quality data."

[0066] Finally, for data segments marked as low-quality, an adaptive classification algorithm is used to identify the specific anomaly type. Based on the data's characteristics, this algorithm intelligently determines the primary anomaly type or types of anomalies contributing to the low quality. Common anomaly types include waveform distortion (distortion of the signal morphology, such as clipping or oscillation), baseline drift (a non-physical, slow shift of the recorded signal's zero line), polarity reversal (the recording direction of one or more components is opposite to physical reality), and site response anomalies (the amplification effect of the local site beneath the station on seismic waves significantly deviates from the expected model). The algorithm outputs the specific anomaly type identified, providing a basis for subsequent data correction or exclusion decisions.

[0067] Furthermore, the waveform anomaly index is calculated in the following manner:

[0068] Detect the waveform amplitude mutation rate. If the amplitude change rate of consecutive sampling points exceeds the preset mutation threshold, it is determined to be a mutation abnormality;

[0069] Analyze the smoothness of the waveform envelope. If the local curvature of the envelope exceeds the curvature threshold, it is determined to be a non-smooth anomaly.

[0070] Identify periodic noise interference and extract the noise energy ratio of the characteristic frequency through Fourier transform. If the ratio exceeds the noise threshold, it is determined to be periodic interference.

[0071] Specifically, the waveform amplitude mutation rate is first detected. For the acceleration time history data, continuous sampling points are scanned in units of a preset time window (such as 0.5 seconds). The ratio of the absolute value of the amplitude difference between adjacent sampling points to the amplitude of the point is calculated and defined as the mutation rate δ:

[0072]

[0073] Among them A i and A i+1 Representing the acceleration amplitudes at the i-th and i+1-th sampling points, respectively, this formula calculates the strength of the relative amplitude change. The system scans a continuous sequence of sampling points within a preset time window (e.g., 0.5 seconds). If the average δ value calculated for several consecutive sampling points (e.g., three points) exceeds the preset mutation threshold (typically 1.5), an amplitude mutation anomaly is determined to have occurred within that time window. This is often caused by transient instrument failures or external transient interference, resulting in sharp glitches on the waveform.

[0074] Envelope smoothness analysis aims to evaluate the regularity of the overall waveform profile. First, the envelope E(t) of the acceleration time history data A(t) is extracted using the Hilbert transform. Then, the local curvature C of the envelope is calculated within the sliding window:

[0075]

[0076] Here, dE / dt and d²E / dt² are the first and second derivatives of the envelope E(t) with respect to time, respectively. The physical meaning of curvature C is the measure of the degree of curvature of the envelope (the inverse of the curvature radius). Calculate the maximum curvature value C of the envelope E(t) within the window max , if C max If the curvature exceeds the preset threshold (recommended value is 0.05), the window is considered to have a non-smooth anomaly. This anomaly usually manifests as an irregular jagged or stepped envelope, which may be caused by continuous mechanical vibration of the station or poor instrument contact.

[0077] The goal of identifying periodic noise interference is to detect whether there is persistent interference of a specific frequency in the recording. The acceleration time history data A(t) is Fourier transformed to obtain its spectrum F(ω). A specific characteristic frequency band Ω is selected (for example, common power frequency interference may be at 50±2 Hz or 100±5 Hz), and the ratio η of the noise energy in this characteristic frequency band to the total signal energy is calculated:

[0078]

[0079] The numerator calculates the integral of the signal energy within the characteristic frequency band Ω, while the denominator calculates the integral of the total signal energy within the entire analysis frequency band. If the calculated value η exceeds the preset noise threshold (typically 5%), significant periodic noise interference is detected. A typical application scenario involves 50Hz or 60Hz electromagnetic interference generated by high-voltage transmission lines near a station, which will produce a noticeable peak at the corresponding frequency on the spectrum graph.

[0080] The final waveform anomaly index is calculated by statistically analyzing the frequency and severity of the three types of anomalies in the analysis records, performing a weighted comprehensive calculation, and normalizing the results to a score range of 0 to 1. The closer the score is to 1, the more severe the waveform morphology abnormality.

[0081] Furthermore, the adaptive classification algorithm includes:

[0082] Time domain analysis: Detects data integrity loss and baseline drift through a sliding window. Baseline drift is determined when the slope of the linear trend of the displacement time history exceeds the drift threshold.

[0083] Frequency domain analysis: Calculate the root mean square error between the frequency response and the standard station response curve. If the error exceeds the offset threshold, the frequency response is determined to be offset.

[0084] Polarization analysis: Calculates the azimuth deviation of the epicenter based on the eigenvectors of the three-component covariance matrix. If the deviation exceeds the azimuth tolerance, it is considered a polarity reversal.

[0085] Site response analysis: The H / V spectrum ratio coefficient of variation is used to identify site response anomalies. If the coefficient of variation exceeds the stability threshold, it is judged as abnormal.

[0086] Specifically, the adaptive classification algorithm includes four core analysis modules for accurately identifying the abnormal types of low-quality data. The time domain analysis module uses a sliding window technique to traverse the acceleration time history data, and the window length is usually set to 1 second. In the data integrity test, the proportion ρ of invalid values ​​(such as NaN or zero values) in each window is counted. When ρ exceeds the missing threshold (such as 5%), it is marked as a data missing anomaly. For baseline drift detection, the acceleration time history must first be integrated twice to obtain the displacement time history D(t), and the displacement time history must be linearly fitted to obtain the slope k. If the absolute value of the fitting slope |k| exceeds the drift threshold (typical value 1×10 -4 m / s 2 ), then it is determined that there is baseline drift. For example, when the power supply voltage of the seismometer is unstable, the displacement time history will show a continuous upward trend.

[0087] The frequency domain analysis module first obtains the standard station response curve S(ω), which is provided by the station calibration file. The frequency response R(ω) of the measured acceleration data is calculated, and the root mean square error of S(ω) in the effective frequency band (e.g., 0.1-35Hz) is calculated as follows:

[0088]

[0089] Where N is the number of frequency points. Frequency response deviation is determined when ε exceeds a deviation threshold (e.g., 3dB). This is often caused by sensor sensitivity loss, resulting in a decrease in high-frequency response.

[0090] The polarization analysis module focuses on the 0.5-2 second time window after the first arrival of the P wave. The three-component acceleration data is extracted to form a data matrix X, and its covariance matrix is ​​calculated. , n is the number of sampling points in the window. Solve the eigenvalues ​​of the covariance matrix λ1≥λ2≥λ3 and the corresponding eigenvector v1. The azimuth angle of the main eigenvector v1 is θ = arctan(v 1,y / v 1,x ), the deviation from the station's measured azimuth θ0 is Δθ = |θ – θ0|. If Δθ > 30°, polarity reversal is determined. For example, incorrect sensor wiring can cause vertical component phase reversal.

[0091] The site response analysis module calculates the Fourier spectrum ratio H / V(f) of the horizontal and vertical components. The complete record is divided into M sub-windows (M ≥ 10 is recommended) and the spectrum ratio curve H / V of each sub-window is calculated. m (f). At the characteristic frequency point f j Calculate the coefficient of variation at:

[0092]

[0093] Where σ and μ represent the standard deviation and mean respectively. If CV(f j )>0.3 (stability threshold), the site response is judged to be abnormal. For example, when soil liquefaction occurs under the station, the spectrum ratio curve will fluctuate violently.

[0094] These four analysis modules run in parallel and output the abnormality type discrimination results. For example, a station record shows |k|=5×10 -4 m / s² (baseline drift) and Δθ = 45° (polarity reversal), it is comprehensively determined that the record has a complex abnormality.

[0095] Furthermore, the polarization analysis specifically includes:

[0096] Extract the three-component acceleration time history within the P wave initial motion window;

[0097] Calculate the eigenvalues ​​and principal eigenvector directions of the covariance matrix within the window;

[0098] Polarity reversal is determined based on the angle deviation between the main eigenvector and the measured azimuth of the station. If the angle deviation exceeds 30°, the correction process is triggered.

[0099] Specifically, the three-component acceleration time history data are accurately extracted within a specific time window after the first arrival of the P wave. This window length is typically set to 0.5 to 2 seconds. The time of the first arrival of the P wave is automatically determined using the ratio of the short-term average to the long-term average combined with the AIC criterion. The three-component data within the window form a three-dimensional data matrix X with dimensions n × 3 (n is the number of sampling points in the window, e.g., 50–200 sampling points at a 100 Hz sampling rate).

[0100] Then calculate the covariance matrix of the data matrix, its mathematical expression is:

[0101]

[0102] The covariance matrix is ​​a 3×3 real symmetric matrix, with the diagonal elements representing the variance of each component and the off-diagonal elements representing the covariance between components. Eigenvalue decomposition is used to determine the eigenvalues ​​λ1, λ2, and λ3 of this matrix (arranged in descending order: λ1 ≥ λ2 ≥ λ3) and their corresponding eigenvectors v1, v2, and v3. The eigenvector v1 corresponding to the largest eigenvalue λ1 is the dominant vibration direction of the signal, reflecting the dominant energy direction of seismic wave propagation.

[0103] The azimuth of the epicenter is calculated based on the main eigenvector v1, and the calculation formula of its horizontal projection azimuth θ is:

[0104]

[0105] Where v1x and v 1y Represent the horizontal components of v1 in the east-west (EW) and north-south (NS) directions respectively. Comparing the calculated azimuth θ with the epicenter azimuth θ0 measured at the station (determined by the station coordinates and the epicenter position), the azimuth deviation is obtained:

[0106]

[0107] When Δθ exceeds 30°, a polarity reversal anomaly is detected. For example, during a certain earthquake, the measured epicenter azimuth θ0 was 120°, but polarization analysis calculated θ=75°. In this case, Δθ=45°, which is greater than 30°, and the system determines a horizontal component polarity anomaly. After triggering the correction process, the EW and NS component data are automatically multiplied by -1 to reverse the polarity. After the reversal, the azimuth deviation is recalculated to within 5°, verifying the effectiveness of the correction. This process can resolve data reversal issues caused by incorrect sensor wiring.

[0108] Furthermore, it also includes:

[0109] Establish an indicator weight optimization model: Use a lightweight neural network to dynamically adjust the weight distribution of multi-dimensional indicators;

[0110] The neural network takes the waveform time domain characteristics and frequency domain characteristics as input, outputs the weight coefficients of each indicator, and optimizes the weight distribution strategy through back propagation.

[0111] Specifically, the core of the model is to dynamically adjust the weight distribution of multi-dimensional quality assessment indicators through a lightweight neural network to solve the problem that traditional fixed weights cannot adapt to different earthquake event characteristics and station environments. The input layer of the neural network receives the time domain features and frequency domain feature vectors extracted from the original three-component acceleration time history. The time domain features include statistical parameters such as the mean μ, standard deviation σ, skewness γ, kurtosis κ of the acceleration time history, as well as the time series change rate of the waveform anomaly index; the frequency domain features include the main frequency f after Fourier transform. m , bandwidth B, spectral entropy H s Key parameters such as , all input features need to be normalized by Z-score.

[0112] The neural network adopts a three-layer fully connected structure: input layer (feature dimension d=32), hidden layer (number of neurons h=16, ReLU activation function) and output layer (number of neurons equal to the number of indicators k=5, Softmax activation function). The output layer generates the weight coefficient vector W=(w1,w2,w3,w4,w5) of each quality indicator (data integrity I1, signal-to-noise ratio I2, waveform anomaly index I3, frequency response deviation I4, polarity consistency I5), satisfying ∑w i = 1. The final quality score is calculated by the weighted formula Q = ∑(w i·S i ) is calculated, where S i is the normalized score of the i-th indicator.

[0113] The training process uses the adaptive momentum optimization algorithm (Adam), and the loss function is defined as:

[0114]

[0115] Where N is the batch sample size, λ=0.01 is the regularization coefficient, Score the prediction quality of the i-th sample, The true quality score of the i-th sample is scored. During training, the network parameters θ are dynamically adjusted through the back propagation algorithm. The mean absolute error (MAE) of the validation set is calculated after each round of iteration. When the MAE drops less than 10 for 5 consecutive rounds, the training set is ranked as the best. -4 This dynamic weighting mechanism significantly improves the adaptability of quality assessment. For example, it automatically increases the weight of the waveform anomaly index (typical value w3≈0.35) in near-field seismic records, while strengthening the weight of the signal-to-noise ratio indicator (w2≈0.40) in far-field records.

[0116] Furthermore, the training step of the neural network includes:

[0117] Construct a mixed sample set: including synthetic data simulating environmental noise and measured low-quality recorded data;

[0118] Feature enhancement: Expand sample diversity through time-frequency domain data enhancement technology;

[0119] Transfer learning: Adapt the pre-trained temporal feature extraction model to the weight distribution task, freezing the underlying convolutional layers and fine-tuning the fully connected layers.

[0120] Specifically, the training starts with constructing a mixed sample set, which contains two types of data sources: one is the strong earthquake records synthesized by the physical model, using random source parameters (magnitude M w =4.0~8.0, epicenter distance R=10~200km) superimposed on an environmental noise model (spectral characteristics such as traffic vibration and industrial interference); second, measured low-quality records obtained from the National Strong Motion Network Center, covering typical abnormal cases such as waveform distortion and baseline drift. The total sample set is approximately 100,000, with a ratio of 7:3 between synthetic data and measured data. Each sample is annotated with a true quality score Q true (The average value is obtained by independent scoring by 5 experts).

[0121] During the feature enhancement phase, multidimensional data augmentation techniques are used to improve model generalization. In the time domain, random time stretching (scaling factor α∈[0.8,1.2]) and amplitude scaling (scaling factor β∈[0.7,1.3]) are applied to the acceleration time histories to simulate waveform distortion due to different propagation paths. In the frequency domain, a controllable frequency offset (offset Δf ≤ 5% of the Nyquist frequency) and band-limited noise injection (with a signal-to-noise ratio dynamic range of 30-50dB) are introduced through Fourier transform. For special scenarios such as periodic interference, a generative adversarial network (GAN) is used to synthesize enhanced samples with specific noise spectra (e.g., 50Hz power frequency interference). The resulting sample size is expanded to 1.5 times the original data, significantly improving the ability to identify rare anomalous patterns.

[0122] The transfer learning phase uses a pre-trained time series feature extraction model as the infrastructure. Specifically, the InceptionTime model, pre-trained on the UCR time series classification dataset, is used. Its underlying layer contains four convolutional modules (with kernel sizes of 10, 20, 40, and 80 sampling points, respectively), which can capture multi-scale waveform features. During migration, the weight parameters of the three underlying convolutional layers are frozen, and only the last convolutional layer and the newly added fully connected layer (dimensional configuration: input 1024 → hidden layer 256 → output layer 5) are fine-tuned. Fine-tuning uses a layered learning rate strategy: the base layer learning rate η base =10 -5 , fine-tuning layer η fine =10 -3 , the optimizer uses Nesterov momentum SGD (momentum coefficient μ=0.9). After each round of training, the weighted score error of the validation set is calculated:

[0123]

[0124] Where K is the number of samples in the validation set, Score the prediction quality of the k-th sample in the validation set, Score the true quality of the kth sample in the validation set. val The early stopping mechanism was triggered when the probability dropped by less than 0.001 for three consecutive epochs. After 200 rounds of training, the model achieved a mean absolute error of 0.08 on the test set, and the consistency of weight distribution with expert decision-making reached 92%.

[0125] Furthermore, it also includes a dynamic correction process:

[0126] For labeled low-quality data, filter banks are selected based on the anomaly type;

[0127] Waveform distortion uses an adaptive bandpass filter, and the passband range is dynamically adjusted according to the signal's main frequency;

[0128] The baseline drift was corrected for the displacement time course using piecewise cubic spline fitting;

[0129] Polarity inversion automatically flips the direction of channel data using polarization analysis results.

[0130] Specifically, the dynamic correction process calls a specific correction module according to the identified abnormality type. For waveform distortion abnormalities, an adaptive bandpass filter is used for processing: First, the main frequency of the acceleration time history is calculated by the Welch power spectrum estimation method. f m (the frequency corresponding to the maximum spectrum value is taken), the filter passband range is dynamically set to [0.6 f m ,1.4 f m ], with a stopband attenuation of no less than 40dB. For example, when the main frequency of a recording is detected to be 5Hz, a 4th-order Butterworth filter with a passband of 3-7Hz is automatically generated, effectively eliminating high-frequency noise interference while retaining the main frequency energy of the seismic motion.

[0131] For baseline drift anomalies, displacement time history reconstruction is performed: the original acceleration time history is double numerically integrated to obtain the displacement time history D ( t ), and a piecewise cubic spline function is used to fit its trend term. Specifically, the interval is divided into 1 second nodes, and a cubic polynomial is constructed in each interval:

[0132]

[0133] The coefficient a j , b j , c j , d j Solve by boundary conditions (displacement and velocity continuity at nodes). D ( t ) minus the trend term S ( t ), eliminating slow zero drift caused by instrument temperature drift (typical drift is reduced from 15mm to 0.5mm).

[0134] For polarity reversal anomalies, the data direction is automatically flipped according to the polarization analysis results: when the system determines that the polarity of a component is wrong (such as the azimuth deviation of the EW component exceeds the limit), the full time course data of the component is multiplied by -1. After correction, the eigenvector azimuth of the covariance matrix is ​​recalculated. i new , verify that it is consistent with the measured azimuth i Deviation Δ from 0 i new <5°, ensuring the validity of the correction. For example, the original data of the NS component of a station results in Δ i =152°, after flipping Δ i new =3°.

[0135] Furthermore, the dynamic correction process also includes:

[0136] Iterative optimization mechanism: With the signal-to-noise ratio improvement rate as the objective function, the filter parameters are adjusted by the gradient descent method until the output signal meets the signal-to-noise ratio threshold or the maximum number of iterations is reached.

[0137] Specifically, the iterative optimization mechanism is based on the signal-to-noise ratio improvement rate As the objective function, SNR new is the original signal-to-noise ratio, SNR new To obtain the new signal-to-noise ratio, the filter parameters are adjusted using the gradient descent method. Specifically, the cutoff frequency of the bandpass filter is f c Calculating gradients , update the parameters along the gradient direction until or ≥0.3 or the maximum number of iterations reaches 20. Experiments show that the signal-to-noise ratio of a typical waveform distortion record can be improved from the original 15dB to over 35dB after 3-5 iterations.

[0138] Furthermore, confidence assessment is also included:

[0139] A confidence model is constructed based on the multi-dimensional indicator score deviation, historical identification accuracy and station reliability;

[0140] Output confidence score. If the score is lower than the reliability threshold, the manual review process is triggered.

[0141] Specifically, the confidence evaluation mechanism is used to quantify the reliability of the quality identification results. The core of the confidence model consists of three parameters: the multi-dimensional index score deviation reflects the anomaly detection consistency of the current record, and the actual score S of each index is calculated. i and the historical average score μ i The normalized difference:

[0142]

[0143] Among them, k is the number of indicators (taken as 5), σ i is the standard deviation of the indicator in the historical data set. For example, if the frequency response deviation I4 of a station is μ4=0.85 and σ4=0.12, and the new record S4=0.52, the deviation contribution is (0.33 / 0.12)=2.75 times the standard deviation.

[0144] The historical recognition accuracy parameter A is calculated based on time-decay weighting, with more recent recognition results given a higher weight:

[0145]

[0146] Where M is the number of historical records, t m Indicates the number of months from the present time for the mth record, with a decay factor of α=0.9, A m A1 is a binary accuracy flag for manual verification (correct is 1, incorrect is 0). If a station has three identifications in the past six months: correct two months ago (A1=1), incorrect four years ago (A2=0), and correct six months ago (A3=1), then A=(0.9 2 ×1+0.9 4 ×0+0.9 6 ×1) / (0.81+0.656+0.531)=0.61.

[0147] The station reliability parameter R is dynamically evaluated through the equipment status log, which includes three types of deduction items: basic failure rate R0=0.3×(1-failure number / 10), calibration overdue rate R c =0.4×(1-number of months exceeded / 24), major accident record R a =0.3×(1-number of accidents / 5). If a station experiences two failures within two years, is overdue for calibration by 8 months, and has no major accidents, then R=0.3×(1-0.2)+0.4×(1-0.33)+0.3×1=0.77.

[0148] The final confidence score integrates three parameters:

[0149]

[0150] When the C value falls below the reliability threshold of 0.65, manual review is triggered. For example, if the recognition result is D = 0.35 (moderate deviation), A = 0.61 (recent low accuracy), and R = 0.77 (good device condition), the calculated C = 0.4 × 0.65 + 0.3 × 0.61 + 0.3 × 0.77 = 0.678, and no review is required. However, if D rises to 0.55 (serious deviation), the calculated C = 0.4 × 0.45 + 0.3 × 0.61 + 0.3 × 0.77 = 0.582, and the system automatically transfers the result to manual processing.

[0151] Once the manual review process is initiated, the system generates a visual report: abnormal sections are marked on the waveform time-history graph (e.g., periods of baseline drift highlighted in red), along with the algorithm's judgment rationale (e.g., "displacement time-history slope 0.00015 > threshold 0.0001") and correction suggestions (segmented spline fitting is recommended). The interactive interface supports manual corrections: engineers can adjust the anomaly type label (e.g., changing a record misclassified as "waveform distortion" to "periodic interference"). These changes are fed back to the training set in real time for model optimization. The typical feedback cycle is weekly updates to model parameters to continuously improve recognition accuracy.

[0152] Furthermore, the manual review process includes:

[0153] Generate visual reports: mark the waveform abnormality location, algorithm judgment basis and correction suggestions;

[0154] An interactive interface is provided to support manual correction of labels and feedback to the training set for optimizing the recognition model.

[0155] Specifically, when the confidence assessment model determines that manual intervention is required, the system automatically generates a structured visual report consisting of three core sections: waveform anomaly location markers, an explanation of the algorithm's judgment basis, and correction recommendations. In the waveform display area, the original three-component acceleration time history data is plotted as colored curves (typically blue for vertical components and red / green for horizontal components). Anomalous sections are highlighted with highlighted blocks—for example, a red semi-transparent overlay indicates baseline drift anomaly at 12.5-13.2 seconds on the time axis, while the abnormal energy peak at 50 Hz is circled in the spectrum. The judgment basis section uses natural language to explain the algorithm's decision logic, such as "EW component azimuth deviation Δθ = 42° > 30° threshold" or "H / V spectral ratio coefficient of variation CV (2.5Hz) = 0.38 > 0.3 threshold." Correction recommendations provide specific technical solutions, such as "suggesting segmented cubic spline fitting correction for baseline drift" and "applying a 50 Hz notch filter for periodic interference."

[0156] The interactive interface design supports engineers to manually correct the system's judgment results. The original waveform and algorithm mark are displayed on the left side of the interface, and the operation panel is provided on the right: including anomaly type checkbox (waveform distortion, site anomaly, etc. can be reselected), polarity flip button (supports clicking the "reverse" icon for a single component) and parameter adjustment slider (such as the signal-to-noise ratio threshold range of 0-60dB). After the engineer completes the correction, click the "Confirm Feedback" button to trigger the data storage process. The correction record is saved in a structured format, including the original data ID, manual correction label, modified parameter value (such as changing the automatically identified polarity reversal to a site response anomaly) and the operator's signature. The storage format is, for example:

[0157] <record id="EQ20240710_0856">

[0158] <Original_Label>Polarity_Inversion< / Original_Label>

[0159] <Corrected_Label>Site_Response_Anomaly< / Corrected_Label>

[0160] <Modified_Parameters>

[0161] <H_V_Threshold>0.35→0.40< / H_V_Threshold>

[0162] < / Modified_Parameters>

[0163] <operator> Engineer_Zhang< / operator>

[0164] < / record>

[0165] The corrected data is synchronized to the training set update queue in real time. The system triggers a model optimization cycle every time 100 new samples are accumulated. The incremental learning strategy is used during optimization: the weights of the neural network base layer are frozen, and only the fully connected parameters of the output layer are fine-tuned. The loss function is fine-tuned to add manual correction weight coefficients:

[0166]

[0167] in βi The weight is set according to the engineer's professional title (1.0 for senior engineer and 0.7 for intermediate engineer). is the index weight vector of the newly generated i-th sample, W i human is the manually adjusted indicator weight vector for the i-th sample. After 200 iterations, the model's accuracy on the validation set improved by approximately 12%. For example, the false positive rate for site anomalies at a certain station dropped from 18.7% to 6.2%. Historical data shows that six months of continuous manual feedback can reduce the system's overall false positive rate by 23.5%.

[0168] Although the embodiments of the present invention have been disclosed above, they are not limited to the applications listed in the description and implementation methods. They can be fully applied to various fields suitable for the present invention. For those familiar with the art, additional modifications can be easily implemented. Therefore, without departing from the general concept defined by the claims and the scope of equivalents, the present invention is not limited to the specific details and embodiments shown and described herein.

Claims

1. A method for identifying low-quality data of strong earthquake records, characterized in that: The following steps are involved: Obtain three-component acceleration time history data of strong vibration records; Construct a multi-dimensional dynamic quality assessment index system, including data integrity, signal-to-noise ratio dynamic range, waveform anomaly index, frequency response deviation and polarity consistency coefficient; Performing real-time quality scoring on the acceleration time history data based on the indicator system, and marking it as low-quality data if the score is lower than a dynamic threshold; Adopting an adaptive classification algorithm to identify abnormal types of the marked low-quality data, wherein the abnormal types include at least one of waveform distortion, baseline drift, polarity reversal and field response abnormality; The waveform anomaly index is calculated as follows: Detect the waveform amplitude mutation rate. If the amplitude change rate of consecutive sampling points exceeds the preset mutation threshold, it is determined to be a mutation abnormality; Analyze the smoothness of the waveform envelope. If the local curvature of the envelope exceeds the curvature threshold, it is determined to be a non-smooth anomaly. Identify periodic noise interference and extract the noise energy ratio of the characteristic frequency through Fourier transform. If the ratio exceeds the noise threshold, it is determined to be periodic interference; The adaptive classification algorithm includes: Time domain analysis: Detects data integrity loss and baseline drift through a sliding window. Baseline drift is determined when the slope of the linear trend of the displacement time history exceeds the drift threshold. Frequency domain analysis: Calculate the root mean square error between the frequency response and the standard station response curve. If the error exceeds the offset threshold, the frequency response is determined to be offset. Polarization analysis: Calculates the azimuth deviation of the epicenter based on the eigenvectors of the three-component covariance matrix. If the deviation exceeds the azimuth tolerance, it is considered a polarity reversal. Site response analysis: The H / V spectrum ratio coefficient of variation is used to identify site response anomalies. If the coefficient of variation exceeds the stability threshold, it is judged as abnormal.

2. The method for identifying low-quality data of strong motion recordings according to claim 1, wherein: The polarization analysis specifically includes: Extract the three-component acceleration time history within the P wave initial motion window; Calculate the eigenvalues ​​and principal eigenvector directions of the covariance matrix within the window; Polarity reversal is determined based on the angle deviation between the main eigenvector and the measured azimuth of the station. If the angle deviation exceeds 30°, the correction process is triggered.

3. The method for identifying low-quality data of strong motion recordings according to claim 1, wherein: Also includes: Establish an indicator weight optimization model: Use a lightweight neural network to dynamically adjust the weight distribution of multi-dimensional indicators; The neural network takes the waveform time domain characteristics and frequency domain characteristics as input, outputs the weight coefficients of each indicator, and optimizes the weight distribution strategy through back propagation.

4. The method for identifying low-quality data of strong motion recordings according to claim 3, wherein: The training steps of the neural network include: Construct a mixed sample set: including synthetic data simulating environmental noise and measured low-quality recorded data; Feature enhancement: Expand sample diversity through time-frequency domain data enhancement technology; Transfer learning: Adapt the pre-trained temporal feature extraction model to the weight distribution task, freezing the underlying convolutional layers and fine-tuning the fully connected layers.

5. The method for identifying low-quality data of strong motion recordings according to claim 1, wherein: Also includes dynamic correction process: For labeled low-quality data, filter banks are selected based on the anomaly type; Waveform distortion uses an adaptive bandpass filter, and the passband range is dynamically adjusted according to the signal's main frequency; The baseline drift was corrected for the displacement time course using piecewise cubic spline fitting; Polarity inversion automatically flips the direction of channel data using polarization analysis results.

6. The method for identifying low-quality data of strong motion recordings according to claim 5, wherein: The dynamic correction process also includes: Iterative optimization mechanism: Taking the signal-to-noise ratio improvement rate as the objective function, the filter parameters are adjusted by the gradient descent method until the output signal meets the signal-to-noise ratio threshold or reaches the maximum number of iterations.

7. The method for identifying low-quality data of strong motion recordings according to claim 1, wherein: Also includes confidence assessment: A confidence model is constructed based on the multi-dimensional indicator score deviation, historical identification accuracy and station reliability; Output confidence score. If the score is lower than the reliability threshold, the manual review process is triggered.

8. The method for identifying low-quality data of strong motion records according to claim 7, wherein: The manual review process includes: Generate visual reports: mark the waveform abnormality location, algorithm judgment basis and correction suggestions; An interactive interface is provided to support manual correction of labels and feedback to the training set for optimizing the recognition model.

Images