A distributed acoustic sensing seismic monitoring data stream processing method and system

Abstract

The application discloses a kind of distributed acoustic sensing earthquake monitoring data stream processing method and system, belong to earthquake monitoring and data processing technical field.The method of the application includes: after the original earthquake monitoring data is preprocessed, it is sliced using sliding processing time window, and the data segment suitable for real-time detection is obtained;Design adaptive short-time average / long-time average detection algorithm, detect the data segment, and dynamically adjust the adaptive parameter in the algorithm according to the data characteristics;Finally, the initial candidate list of earthquake event is executed multi-channel fusion detection, to generate the final detection list of earthquake event.The application is based on distributed stream data processing architecture, combined with the innovative advantage of adaptive earthquake detection algorithm, can improve the accuracy of earthquake detection while ensuring the real-time of detection result, can provide efficient and reliable technical support for large-scale earthquake monitoring task.

Description

Technical Field

[0001] This invention belongs to the field of earthquake monitoring and data processing technology, and in particular relates to a distributed acoustic sensing earthquake monitoring data stream processing method and system. Background Technology

[0002] Earthquake monitoring technology, as an important branch of geophysics, has evolved from traditional seismograph networks to modern digital monitoring systems. Traditional earthquake monitoring mainly relies on distributed seismograph networks, using the STA / LTA (short-time averaging / long-time averaging) algorithm for earthquake event detection. This method performs well in high signal-to-noise ratio (SNR) scenarios and with the arrival of impactful P-waves and S-waves, but its sensitivity is low in low SNR, overlapping events, and complex noise environments. In recent years, with the development of machine learning technology, deep learning methods have shown superior performance in earthquake detection, but they face challenges such as dataset differences and the diversity of evaluation methods when deployed in heterogeneous environments.

[0003] Distributed acoustic sensing (DAS) technology, as an emerging earthquake monitoring method, offers a revolutionary solution to the limitations of traditional monitoring technologies. DAS technology can convert ordinary optical fibers into a continuous array of seismic sensors, detecting vibration signals along the fiber optic cable through the Rayleigh scattering effect of laser pulses. A single fiber can provide thousands of virtual sensor channels, achieving a spatial resolution of up to 1 meter and a sampling frequency exceeding 1000 Hz. Compared to traditional seismograph networks, DAS technology has significant advantages: First, DAS systems utilize existing telecommunications fiber optic infrastructure, eliminating the need for additional power supply and maintenance equipment in the field, significantly reducing long-term operating costs; second, DAS can achieve dense sampling over a kilometer-scale range, providing unprecedented spatial resolution for near-surface monitoring; furthermore, DAS systems have been successfully applied to remote earthquake detection, such as detecting the Fiji Islands earthquake within a city fiber optic network more than 9000 kilometers from the epicenter, demonstrating its reliability in complex noisy environments.

[0004] However, DAS technology still faces significant technical challenges in practical applications, primarily stemming from the mismatch between the massive amounts of data it generates and existing processing technologies. The amount of data generated by DAS systems grows exponentially at high sampling rates. A typical DAS array contains thousands of channels with sampling frequencies exceeding 1000Hz, typically generating terabytes of data daily. Traditional data processing systems, employing single-machine or small-scale cluster architectures, lack horizontal scalability and struggle to meet real-time processing demands. Furthermore, traditional batch processing models suffer from processing delays ranging from minutes to hours, severely impacting the response speed of earthquake early warning systems, which require issuing warnings within 4 seconds of P-wave arrival. In addition, traditional STA / LTA earthquake detection algorithms use fixed parameter configurations, failing to adapt to varying geological environments and noise conditions. This results in low detection accuracy and high false alarm and false negative rates in complex geological environments. Modern earthquake monitoring systems need to maintain detection sensitivity while avoiding false alarms caused by noise sources such as industrial activities and equipment malfunctions. Summary of the Invention

[0005] The purpose of this invention is to provide a distributed acoustic sensing seismic monitoring data stream processing method and system. It aims to solve key technical problems in the prior art, such as insufficient real-time processing capability of large-scale DAS data, poor horizontal scalability of distributed systems, excessively long real-time seismic event detection response time (unable to achieve second-level detection response), and poor environmental adaptability of fixed parameter detection algorithms, by using adaptive short-time averaging / long-time averaging detection algorithms and multi-channel fusion technology.

[0006] To achieve the above-mentioned objectives, the present invention specifically adopts the following technical solution:

[0007] In a first aspect, the present invention provides a method for processing distributed acoustic sensing seismic monitoring data streams, comprising the following steps:

[0008] S1. Receive raw seismic monitoring data collected by distributed acoustic sensors through the data source, and use a distributed acoustic sensor data parser to parse and preprocess the raw seismic monitoring data;

[0009] S2. Apply a sliding processing time window to slice the preprocessed seismic monitoring data, and use the seismic monitoring data in each window as a data segment suitable for real-time detection;

[0010] S3. In the seismic detection processor, an adaptive short-time average / long-time average detection algorithm is executed on the data segments, and the adaptive parameters in the algorithm are dynamically adjusted according to the characteristics of the seismic monitoring data during the detection process, and an initial candidate list of seismic events is output.

[0011] S4. Perform multi-channel fusion detection on the initial candidate list of earthquake events to generate the final detection list of earthquake events.

[0012] Based on the above scheme, each step can be implemented in the following preferred manner.

[0013] As a preferred embodiment of the first aspect above, in step S3, the adaptive parameters include short-term average window length, long-term average window length, and detection trigger threshold.

[0014] The short-time average window length is obtained by limiting the calculated first base value to between the first minimum value and the first maximum value through a constraint function. The first maximum value is twice the sampling rate. The first base value is obtained by multiplying the preset first constant, the first intermediate term, and the second intermediate term. The first intermediate term is the ratio of the sampling rate to the first component. The first component is the maximum value between the main frequency and 1. The second intermediate term is obtained by adding the second component to 1. The second component is obtained by multiplying the signal variability index and the second constant.

[0015] The long-term average window length is obtained by limiting the calculated second basic value to between the second minimum and the second maximum value through a constraint function. The second minimum value is an integer multiple of the short-term average window length, and the second maximum value is an integer multiple of the sampling rate. The second basic value is obtained by multiplying the third intermediate term, the fourth intermediate term, and the fifth intermediate term. The third intermediate term is an integer multiple of the short-term average window length, the fourth intermediate term is obtained by adding the third component to 1, the third component is a multiple of the noise level, and the fifth intermediate term is the difference between the third constant and the signal stationarity index.

[0016] The detection trigger threshold is obtained by limiting the calculated third basic value to between the third minimum and the third maximum value through a constraint function. Both the third minimum and the third maximum value are constants. The third basic value is obtained by multiplying the fourth constant, the sixth intermediate term, the seventh intermediate term, and the eighth intermediate term. The sixth intermediate term is obtained by adding the fourth component to 1. The fourth component is a multiple of the noise level. The seventh intermediate term is calculated based on the signal-to-noise ratio. The eighth intermediate term is obtained by adding the fifth component to 1. The fifth component is a multiple of the signal variability index.

[0017] Furthermore, the short-time average window length The calculation method is as follows:

[0018]

[0019]

[0020] in, Represents the limiting function; This is the first basic value; It is the first minimum value; It is the first maximum value; It is the first constant; It is the first intermediate term; Indicates the sampling rate; Indicates the clock speed; This indicates taking the maximum value; Indicates the first component; It is the second intermediate term; Indicates the second component; It is the second constant; This is an indicator of signal variability.

[0021] Furthermore, the long-term average window length The calculation method is as follows:

[0022]

[0023]

[0024] in, The second basic value, the second minimum value is The second maximum value is The third intermediate term is ; It is the fourth intermediate term; The third component, Indicates noise level; The fifth intermediate term, It is the third constant; It is an indicator of signal stability.

[0025] Furthermore, detect trigger thresholds The calculation methods are as follows:

[0026]

[0027]

[0028] in, The third basic value, the third minimum value is The third maximum value is ; It is the fourth constant; The sixth intermediate term, It is the fourth component; Indicates the seventh middle term; Indicates the signal-to-noise ratio; The eighth intermediate term, This represents the fifth component.

[0029] As a preferred embodiment of the first aspect mentioned above, the noise level, frequency distribution of earthquake monitoring data, signal variability index, and signal stationarity index are all obtained by analyzing the earthquake monitoring data using a signal feature analyzer.

[0030] As a preferred embodiment of the first aspect above, in step S3, for each channel's data segment, the average signal energy within a short time window and the average signal energy within a long time window are calculated, and the ratio of the two average signal energy values ​​is taken as the STA / LTA ratio. When the STA / LTA ratio is greater than the adaptive detection threshold, it is considered that a seismic event has occurred in that channel and it is added to the initial candidate list of seismic events.

[0031] As a preferred embodiment of the first aspect above, the adaptive detection threshold is obtained by multiplying the calculated detection trigger threshold by a correction term, wherein the correction term is... The result is obtained by adding the ninth intermediate term, which is the product of the sixth component and the preset adjustment sensitivity factor. The sixth component is the result of taking the logarithm of the seventh component. The seventh component is the ratio of two signal-to-noise ratios. The first signal-to-noise ratio is the real-time signal-to-noise ratio calculated from the data segment of the current channel, and the second signal-to-noise ratio is the preset reference signal-to-noise ratio.

[0032] Furthermore, adaptive detection threshold The calculation formula is as follows:

[0033]

[0034] in, This is a preset sensitivity adjustment factor; This represents the real-time signal-to-noise ratio calculated from the data segment of the current channel; This indicates the preset reference signal-to-noise ratio.

[0035] As a preferred embodiment of the first aspect above, in step S3, the adaptive short-time average / long-time average detection algorithm employs a parameter caching mechanism: maintaining a parameter cache mapping table to store the adaptive parameters for each channel; setting a cache expiration time, and when the cache expiration time is reached, recalculating the adaptive parameters and updating the parameter cache mapping table.

[0036] As a preferred embodiment of the first aspect mentioned above, the specific process of performing multi-channel fusion detection on the initial candidate list in step S4 is as follows:

[0037] S41. Extract the earthquake events detected by each channel from the initial candidate list of earthquake events and use them as initial earthquake events. Each initial earthquake event contains six key pieces of information, namely, initial earthquake event ID, channel ID, earthquake start time, earthquake end time, STA / LTA ratio, and confidence level.

[0038] S42. When two initial earthquake events simultaneously satisfy both the temporal correlation constraint and the spatial correlation constraint, the two initial earthquake events are added to the candidate earthquake event list; wherein, the temporal correlation constraint is that the trigger time difference between the two initial earthquake events is less than or equal to the maximum allowed trigger time difference, and the trigger time difference is the difference in the earthquake start time of the two initial earthquake events; the spatial correlation constraint is that the spatial correlation between the two initial earthquake events is greater than the preset spatial correlation threshold.

[0039] S43. For each candidate earthquake event in the candidate earthquake event list, if the same candidate earthquake event is detected by more than Nc channels, then the candidate earthquake event is added to the final detection list as a fused earthquake event. The earthquake start time and earthquake end time of the fused earthquake event are the earliest earthquake start time and the latest earthquake end time of all corresponding candidate earthquake events, respectively. The STA / LTA ratio of the fused earthquake event is the maximum STA / LTA ratio of all corresponding candidate earthquake events. The confidence of the fused earthquake event is obtained by fusing the confidence of all corresponding candidate earthquake events. Wherein, Nc is a preset channel number threshold.

[0040] As a preferred embodiment of the first aspect above, in step S42, the spatial correlation between the two initial seismic events is obtained by multiplying the tenth intermediate term, the time overlap factor, and the distance weighting factor.

[0041] The tenth intermediate term is obtained by dividing the minimum STA / LTA ratio by the maximum STA / LTA ratio between the two initial seismic events;

[0042] The time overlap factor is obtained by dividing the eleventh intermediate term and the twelfth intermediate term. The eleventh intermediate term is the maximum value between 0 and the eighth component. The eighth component is the difference between the minimum earthquake end time and the maximum earthquake start time in two initial earthquake events. The twelfth intermediate term is the minimum value between the ninth component and the tenth component. The ninth component is the difference between the earthquake end time and the earthquake start time of one initial earthquake event. The tenth component is the difference between the earthquake end time and the earthquake start time of another initial earthquake event.

[0043] The distance weighting factor is the maximum value between 0 and the eleventh component. The eleventh component is the difference between 1 and the twelfth component. The twelfth component is obtained by dividing the geographic spatial distance between the channels where the two initial seismic events are located by a preset spatial distance threshold.

[0044] Furthermore, the spatial correlation between the two initial seismic events Calculate using the following formula:

[0045]

[0046]

[0047]

[0048] in, These represent the initial earthquake events. and the initial earthquake event The STA / LTA ratio; Indicates the time overlap factor; Indicates the distance weighting factor; and These are the initial earthquake events. The start and end times of the earthquake; and These are the initial earthquake events. The start and end times of the earthquake; Indicates the initial earthquake event Location and initial earthquake event The geographical distance between the passageways; This indicates the preset spatial distance threshold.

[0049] As a preferred embodiment of the first aspect above, in step S43, the confidence level of the fused seismic event is obtained by dividing the weighted sum of the confidence levels of all candidate seismic events by the sum of all fusion weights; the fusion weight of the candidate seismic event is obtained by multiplying the STA / LTA ratio of the candidate seismic event by the reliability coefficient of the channel and then dividing by a preset normalized benchmark threshold; the reliability coefficient is obtained by dividing the historical effective detection count of the channel by the corrected historical total detection count of the channel, and the corrected historical total detection count of the channel is obtained by adding the historical total detection count of the channel and a smoothing term.

[0050] Furthermore, the confidence levels of fused seismic events The calculation uses the following formula:

[0051]

[0052]

[0053]

[0054] in, This indicates the number of candidate seismic events participating in the fusion; Indicates the first The fusion weights of candidate seismic events; Indicates the first Confidence level of each candidate earthquake event; Indicates the first The STA / LTA ratio of candidate earthquake events; The normalized benchmark threshold represents the significance of candidate seismic events. In this embodiment, the normalized benchmark threshold is taken as 2 to 3 times the adaptive detection threshold. Indicates the reliability coefficient of the channel; This indicates the number of historical valid detections for the channel; This indicates the total number of historical tests conducted on the channel; This indicates a smoothing term to prevent the denominator from being zero.

[0055] In a second aspect, the present invention provides a distributed acoustic sensing earthquake monitoring data stream processing system, comprising:

[0056] The data acquisition module is used to receive raw seismic monitoring data collected by distributed acoustic sensors through a data source;

[0057] The result acquisition module is used to process the raw earthquake monitoring data according to the distributed acoustic sensing earthquake monitoring data stream processing method described in any of the above-mentioned first aspects, and obtain the final detection list of earthquake events.

[0058] Compared with the prior art, the present invention has the following advantages:

[0059] This invention provides an end-to-end real-time processing system architecture for distributed acoustic sensing (DAS) seismic monitoring data. This architecture implements a distributed streaming pipeline from data access via KafkaSource, preprocessing by a distributed acoustic sensor data parser (DASDataParser), to detection within a sliding processing time window. In typical scenarios with thousands of channels and sampling rates greater than or equal to 1 kHz, the overall processing latency is reduced by 2-3 orders of magnitude compared to traditional batch processing solutions. While parallelism scales linearly with the number of channels, throughput increases nearly linearly, maintaining exactly-once consistency and fault tolerance, adapting to the needs of large-scale sensor array deployments. This invention is based on Flink's event-time semantics and checkpointing mechanism, and is engineered from standardized data sources and operator interfaces. It provides horizontally scalable processing capabilities while ensuring real-time performance and consistency, adapting to different DAS deployment models and operational scenarios. Compared with existing production systems that use frameworks such as AQMS / Earthworm for data access and event processing, it can overcome bottlenecks in scalability, end-to-end consistency, and processing latency under conditions of high channel density and ultra-large data volume.

[0060] This invention provides a seismic event detection method combining adaptive parameter optimization and multi-channel spatial correlation analysis. It employs an adaptive short-time / long-time averaging detection algorithm and a signal feature analyzer to achieve parameter adaptation and multi-scale calculation. Multi-channel fusion detection performs spatiotemporal correlation constraints, seismic event fusion, and confidence calculation, supplemented by parameter caching and memory pool optimization to reduce computational and memory overhead. Compared to fixed-parameter short-time / long-time averaging detection algorithms, the detection accuracy is improved by approximately 15-25%, and the false alarm rate is reduced by approximately 30-40%. Simultaneously, multi-channel collaboration effectively suppresses local noise interference, improving detection stability and repeatability. To a certain extent, it overcomes the problems of traditional short-time / long-time averaging detection algorithms in DAS data, which are affected by noise and coupling uncertainties, making threshold setting difficult and prone to false positives and false negatives. Attached Figure Description

[0061] Figure 1 This is a flowchart of the steps of the method of the present invention;

[0062] Figure 2 This is a schematic diagram illustrating the principle of the adaptive short-time average / long-time average detection algorithm in this invention;

[0063] Figure 3 This is a schematic diagram of the multi-channel fusion detection process in this invention;

[0064] Figure 4 This is a flowchart of the earthquake monitoring data analysis and processing in this invention;

[0065] Figure 5 This is a spatiotemporal distribution diagram of the DAS strain wave field after flow cytometry preprocessing provided in an embodiment of the present invention.

[0066] Figure 6 This is a schematic diagram of the residual distribution of automatic seismic phase picking provided in an embodiment of the present invention;

[0067] Figure 7 This is a schematic diagram of the data stream processing results of distributed acoustic sensing earthquake monitoring provided in an embodiment of the present invention;

[0068] Figure 8 The graph shows the end-to-end processing latency performance test results of this invention under different data input scales;

[0069] Figure 9 This is a graph showing the test results of the horizontal scalability of the present invention, where throughput varies with parallelism configuration;

[0070] Figure 10 This is a graph showing the changes in memory resource consumption during a 72-hour long-term stability test of the present invention.

[0071] Figure 11 This is a system block diagram of the present invention. Detailed Implementation

[0072] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Many specific details are set forth in the following description to provide a thorough understanding of the present invention. However, the present invention can be practiced in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below. Technical features in the various embodiments of the present invention can be combined accordingly without mutual conflict.

[0073] In the description of this invention, it should be understood that the terms "first" and "second" are used only for descriptive purposes and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Therefore, a feature defined with "first" and "second" may explicitly or implicitly include at least one of those features.

[0074] To fully leverage the advantages of DAS technology in earthquake monitoring and overcome existing technological limitations, efforts should be focused on addressing issues such as insufficient real-time processing capabilities for massive amounts of data, poor adaptability of traditional detection algorithms, and low efficiency in multi-channel data fusion. This invention provides a distributed acoustic sensing earthquake monitoring data stream processing method. This method, by constructing an efficient distributed stream processing architecture, achieves intelligent real-time processing of massive DAS data and accurate earthquake event detection.

[0075] like Figure 1 As shown, in a preferred embodiment of the present invention, the above-mentioned distributed acoustic sensing seismic monitoring data stream processing method includes the following steps S1 to S4. The specific implementation process of each step will be described in detail below.

[0076] S1. Receive raw seismic monitoring data collected by distributed acoustic sensors (DAS) through the data source, and use a distributed acoustic sensor data parser (DAS Data Parser) to parse and preprocess the raw seismic monitoring data.

[0077] It should be noted that in step S1 of this invention, the DAS data stream is received in real time through a message queue component (Kafka). This component supports a multi-partition parallel consumption mechanism and can handle data traffic of up to millions of records per second, ensuring high-throughput data processing capabilities. The message queue cluster is configured with 3 nodes, each topic has 12 partitions, and the replication factor is set to 2 to ensure data reliability and system fault tolerance. Data preprocessing uses a data mapper for parsing. This mapper is implemented based on the MapFunction of the distributed stream processing framework and has state management and fault recovery capabilities. The data format standardization process converts the raw seismic monitoring data into a standard seismic record format, which includes fields such as channel identifier, timestamp, geographic coordinates, strain rate values, and signal quality indicators, ensuring data structure consistency and ease of subsequent processing. The timestamp correction module corrects the seismic monitoring data timestamps based on the Network Time Protocol (NTP) and the system clock, using a linear interpolation method to handle clock deviations, ensuring time accuracy at the millisecond level and eliminating data timing errors caused by clock drift. The channel mapping module establishes the mapping relationship between physical channels and logical channels, supporting unified management of multi-fiber networks. This module maintains a dynamically updated channel configuration table, which can automatically identify newly added or invalid channels and adjust the processing strategy accordingly.

[0078] Furthermore, this embodiment also partitions the preprocessed seismic monitoring data based on channel identifiers, allocating seismic monitoring data from the same channel to the same processing nodes (computing nodes) for easier storage and to ensure the continuity of time-series computation. Specifically, logical partitioning is achieved through the KeyBy operator of the Flink stream processing framework: using the "channel ID" as the partition key, the unordered or mixed seismic monitoring data stream output from the distributed acoustic sensor data parser is reorganized into logical sub-streams isolated by channel. This operation not only partitions by a single channel ID but also maintains a shared state using "group ID" as the key, ensuring that continuous sampling points from the same physical fiber location are strictly routed to the same parallel processing task slot. This satisfies the requirements of the Short-Time Average / Long-Time Average Detection Algorithm (STA / LTAAlgorithm) for the temporal continuity and state locality of data from the same channel, ensuring the locality principle of data processing and improving computational efficiency. This design allows subsequent processing nodes to access not only historical data of the current channel, but also data of its neighboring channels with low latency, providing a good data foundation for subsequent time series analysis and pattern recognition.

[0079] S2. Apply sliding processing time windows to slice the preprocessed seismic monitoring data, and treat the seismic monitoring data in each window as a data segment suitable for real-time detection.

[0080] It should be noted that step S2 of this invention aims to convert a continuous, unbounded data stream into discrete, bounded data blocks to adapt to the windowed computation requirements of the Adaptive Short-Time Average / Long-Time Average Detection Algorithm (STA / LTA Algorithm). Specifically, this invention configures a sliding processing time window with a fixed length (WindowSize, e.g., 10 seconds) and a sliding step (e.g., 1 second). As time progresses, the window slides forward continuously, and each slide triggers an aggregation computation of the data within the window.

[0081] In this embodiment, the specific process of slicing the preprocessed seismic monitoring data using a sliding processing time window is as follows: This embodiment adopts an optimization strategy based on memory management to reduce the overhead caused by high-frequency window triggering. Unlike the traditional Flink window mechanism where each window stores data independently, this embodiment implements a custom window mechanism of "incremental aggregation + circular buffer". Specifically, the process window function maintains a circular buffer based on time indexing. Newly arrived data is only written to the head of the circular buffer, while expired old data is removed from the tail of the circular buffer. When a window is triggered, the method of this invention does not need to scan all seismic records within the window, but directly reuses the existing data in the circular buffer and only performs incremental calculations on the newly added data. This design reduces the time complexity of window operations from O(N) to O(1) (where N is the window size), greatly improving the processing throughput in high sampling rate (e.g., above 1kHz) scenarios.

[0082] S3. In the Seismic Detection Processor, an adaptive short-time average / long-time average detection algorithm is executed on the data segments. During the detection process, the adaptive parameters in the algorithm are dynamically adjusted according to the characteristics of the seismic monitoring data, and an initial candidate list of seismic events is output.

[0083] It should be noted that in step S3 of the present invention, the adaptive parameters include short-term average window length, long-term average window length, and detection trigger threshold.

[0084] The short-time average window length is obtained by limiting the calculated first base value to between the first minimum value and the first maximum value through a constraint function. The first maximum value is twice the sampling rate. The first base value is obtained by multiplying the preset first constant, the first intermediate term, and the second intermediate term. The first intermediate term is the ratio of the sampling rate to the first component. The first component is the maximum value between the main frequency and 1. The second intermediate term is obtained by adding the second component to 1. The second component is obtained by multiplying the signal variability index and the second constant.

[0085] The long-term average window length is obtained by limiting the calculated second basic value to between the second minimum and the second maximum value through a constraint function. The second minimum value is an integer multiple of the short-term average window length, and the second maximum value is an integer multiple of the sampling rate. The second basic value is obtained by multiplying the third intermediate term, the fourth intermediate term, and the fifth intermediate term. The third intermediate term is an integer multiple of the short-term average window length, the fourth intermediate term is obtained by adding the third component to 1, the third component is a multiple of the noise level, and the fifth intermediate term is the difference between the third constant and the signal stationarity index.

[0086] The detection trigger threshold is obtained by limiting the calculated third basic value to between the third minimum and the third maximum value through a constraint function. Both the third minimum and the third maximum value are constants. The third basic value is obtained by multiplying the fourth constant, the sixth intermediate term, the seventh intermediate term, and the eighth intermediate term. The sixth intermediate term is obtained by adding the fourth component to 1. The fourth component is a multiple of the noise level. The seventh intermediate term is calculated based on the signal-to-noise ratio. The eighth intermediate term is obtained by adding the fifth component to 1. The fifth component is a multiple of the signal variability index.

[0087] Furthermore, in this embodiment, the aforementioned short-time average window length The calculation method is as follows:

[0088]

[0089]

[0090] in, This represents a constraint function, used to restrict the variables of the function. Limited to and between, and These are preset limit values; This is the first basic value; As the first minimum value, in this embodiment, ; It is the first maximum value; As the first constant, in this embodiment, ; It is the first intermediate term; Indicates the sampling rate; Indicates the clock speed; This indicates taking the maximum value; Indicates the first component; It is the second intermediate term; Indicates the second component; As the second constant, in this embodiment, ; This is an indicator of signal variability.

[0091] Furthermore, in this embodiment, the aforementioned long-term average window length The calculation method is as follows:

[0092]

[0093]

[0094] in, The second basic value, the second minimum value is The second maximum value is The third intermediate term is ; It is the fourth intermediate term; The third component, Indicates noise level; The fifth intermediate term, As the third constant, in this embodiment, ; This is an indicator of signal stability.

[0095] Furthermore, in this embodiment, a trigger threshold is detected. The calculation methods are as follows:

[0096]

[0097]

[0098] in, The third basic value, the third minimum value is The third maximum value is ; As the fourth constant, in this embodiment, ; The sixth intermediate term, It is the fourth component; Indicates the seventh middle term; Indicates the signal-to-noise ratio; The eighth intermediate term, This represents the fifth component.

[0099] In this embodiment, the noise level, frequency distribution of earthquake monitoring data, signal variability index, and signal stationarity index are all obtained by analyzing the earthquake monitoring data using a signal characteristics analyzer. Furthermore, the noise level, signal variability index, and signal stationarity index are all normalized, with values ​​ranging from [0, 1]. This signal characteristics analyzer can quantify the original earthquake monitoring data and extract key characteristic indicators of the seismic waves. For example, a wavelet transform-based signal characteristics analyzer first performs wavelet multi-scale decomposition on the original earthquake monitoring data to obtain sub-band signals in different frequency bands. Then, it calculates the energy in the high-frequency sub-bands to estimate the noise level, obtains the frequency distribution characteristics by analyzing the energy distribution of each sub-band, calculates the variance of the energy in a specific sub-band as a signal variability index, and evaluates the signal stationarity by comparing the stability of sub-band statistics over different time periods, thus obtaining the signal stationarity index. This systematically completes the extraction of all four characteristic indicators.

[0100] It should be noted that traditional short-time average / long-time average detection algorithms (STA / LTA algorithms) often suffer from high false alarm rates or missed detections when processing DAS data due to the non-stationarity of background noise. To address this technical problem, this invention introduces a dynamic parameter adjustment mechanism based on the STA / LTA detection algorithm. Specifically, it introduces three adaptive parameters: short-time average window length, long-time average window length, and detection trigger threshold, and designs an adaptive STA / LTA algorithm that can dynamically update the short-time average window length, long-time average window length, and detection trigger threshold based on the real-time calculated signal-to-noise ratio (SNR) and background noise variance. When the background noise is stable, the adaptive STA / LTA detection algorithm automatically extends the LTA window to obtain a more stable baseline; when a sudden high-frequency interference is detected, it automatically increases the detection trigger threshold to suppress false alarms. In addition, this step also integrates vectorized computation technology, which uses the SIMD instruction set of modern CPUs to compute STA / LTA values ​​at multiple time points in parallel, significantly improving computing performance.

[0101] like Figure 2As shown, the adaptive short-time average / long-time average detection algorithm comprises five core steps: data segment input, STA calculation, LTA calculation, adaptive detection threshold calculation, and seismic event detection. These steps are executed sequentially to achieve high-precision seismic event identification. Specifically, for each channel's data segment, the average signal energy within the short-time window (STA) and the average signal energy within the long-time window (LTA) are calculated, and the ratio of these two average signal energy values ​​is taken as the STA / LTA ratio. When the STA / LTA ratio is greater than the adaptive detection threshold, a seismic event is considered to have occurred in that channel, and the channel is added to the initial candidate list of seismic events (Detection Event).

[0102] Furthermore, the aforementioned adaptive detection threshold The calculated detection trigger threshold is obtained by multiplying it by a correction term, where the correction term is... The result is obtained by adding the ninth intermediate term, which is the product of the sixth component and the preset adjustment sensitivity factor. The sixth component is the result of taking the logarithm of the seventh component. The seventh component is the ratio of two signal-to-noise ratios. The first signal-to-noise ratio is the real-time signal-to-noise ratio calculated from the data segment of the current channel, and the second signal-to-noise ratio is the preset reference signal-to-noise ratio.

[0103] In this embodiment, the classic short-time average / long-time average detection algorithm (STA / LTA) based on energy characteristics is used as the basic detection operator to obtain the STA / LTA ratio. The specific calculation formula is as follows:

[0104]

[0105] in, express The STA / LTA ratio at time 1; express The average signal energy within a short time window; express The average signal energy within a long time window; Indicates the short-time window length (number of sampling points); Indicates the length of the long-term window (number of sampling points); Indicates the short window The signal of the moment; Indicates the long-term window The signal of the moment; express The amplitude is used to calculate the instantaneous energy at the current moment; express The amplitude is used to calculate the average energy of the background noise; , These are the time offset index variables within the short-time window and the long-time window, respectively. Specifically, in this embodiment, the short-time window length... Setting the sampling point to 500 points, corresponding to 0.5 seconds of data, this length is sufficient to capture the initial motion characteristics of seismic P-waves without excessively smoothing transient signal changes; long-term window length Setting the sampling points to 30,000, corresponding to 30 seconds of data, this length can effectively estimate the background noise level while avoiding including too many seismic events.

[0106] After obtaining the STA / LTA ratio, a seismic event trigger decision is executed. To overcome the limitations of a fixed threshold in non-stationary noise environments, this embodiment designs an adaptive parameter adjustment strategy, which can dynamically calculate the adaptive detection threshold based on the real-time signal-to-noise ratio. To adapt to different environmental noise conditions, when the STA / LTA ratio exceeds an adaptive detection threshold, a seismic event is considered to have occurred. The formula for calculating this threshold is as follows:

[0107]

[0108] in, This is a preset sensitivity adjustment factor; This represents the real-time signal-to-noise ratio calculated from the data segment of the current channel; This represents a preset reference signal-to-noise ratio. Specifically, in this embodiment, the adaptive detection threshold... It will dynamically adjust based on real-time signal-to-noise ratio conditions; adjust the sensitivity factor. Setting it to 0.3 controls the sensitivity of adaptive adjustment; real-time signal-to-noise ratio. The reference signal-to-noise ratio was calculated by analyzing the power spectral density of the data segment. It is then set to 10dB as a reference standard for standardization.

[0109] It should be noted that in step S3 of this invention, the adaptive short-time averaging / long-time averaging detection algorithm employs a parameter caching mechanism: maintaining a parameter cache mapping table to store the adaptive parameters for each channel; setting a cache expiration time, and when the cache expiration time is reached, recalculating the adaptive parameters and updating the parameter cache mapping table. In this embodiment, the cache expiration time is set to 300,000 milliseconds, meaning that the adaptive parameters are updated every 300,000 milliseconds, and the parameter cache mapping table is also updated.

[0110] Additionally, it should be noted that, to support efficient adaptive detection, this embodiment employs a memory-based circular buffer design for the sliding window buffer. This structure is specifically optimized for the continuous sliding characteristics of seismic data streams, providing extremely low-latency data access while ensuring efficient memory usage. Specifically, the size of the circular buffer is set according to the maximum long-term window length, employing a fixed memory allocation strategy to avoid the overhead of frequent garbage collection by the Java Virtual Machine. As the window slides forward over time, newly arriving seismic data directly overwrites the earliest time-slice data in the buffer. This first-in-first-out update mechanism ensures that the data used for STA / LTA calculations always remains within the latest time window, and the memory operation complexity is only O(1), thus meeting the real-time requirements of large-scale channel parallel detection. A recursive calculation method is used to avoid repeatedly calculating historical data, significantly improving computational efficiency. Furthermore, the adaptive short-term average / long-term average detection algorithm of this invention adopts a parallel processing architecture, allocating independent processing threads to each channel group, and achieving efficient concurrent processing through thread pool management and task scheduling mechanisms.

[0111] S4. Perform multi-channel fusion detection on the initial candidate list of earthquake events to generate the final detection list of earthquake events.

[0112] It should be noted that step S4 of this invention implements a multi-channel fusion detection and seismic event confirmation mechanism. This mechanism is a key link in seismic detection. Its core function lies in leveraging the advantage of dense DAS deployment. By analyzing the correlation between alarms from different channels in terms of spatial distance, temporal consistency, and amplitude correlation, it can effectively integrate STA / LTA detection results from different channels, fusing multiple related single-channel seismic events into a high-confidence fused seismic event. For isolated single-channel seismic events, this invention further marks them as random noise (such as passing vehicles, human knocking, etc.) and filters them out. This multi-channel collaborative mechanism effectively solves the problem of high false alarm rate in traditional single-channel methods under high-noise environments, and can improve the accuracy and reliability of seismic event identification.

[0113] like Figure 3 As shown, multi-channel fusion detection comprises three core sub-steps: single-channel candidate seismic event collection, spatiotemporal correlation verification, and seismic event confirmation. These sub-steps work collaboratively to achieve high-confidence seismic event identification. The specific process of performing multi-channel fusion detection on the initial candidate list is as follows:

[0114] S41. Extract the earthquake events detected by each channel from the initial candidate list of earthquake events and use them as initial earthquake events. Each initial earthquake event contains six key pieces of information: initial earthquake event ID, channel ID, earthquake start time, earthquake end time, STA / LTA ratio, and confidence level.

[0115] S42. When two initial earthquake events simultaneously satisfy both the temporal correlation constraint and the spatial correlation constraint, the two initial earthquake events are added to the candidate earthquake event list; wherein, the temporal correlation constraint is that the trigger time difference between the two initial earthquake events is less than or equal to the maximum allowed trigger time difference, and the trigger time difference is the difference in the earthquake start time of the two initial earthquake events; the spatial correlation constraint is that the spatial correlation between the two initial earthquake events is greater than the preset spatial correlation threshold.

[0116] In this embodiment S42, the time correlation constraint is used to determine its rationality by analyzing the distribution of seismic events within a time window. Specifically, the following judgment criteria are adopted:

[0117]

[0118] in, Indicates the trigger time difference; , These represent the initial earthquake events. and the initial earthquake event The start time of the earthquake; This indicates the maximum allowed trigger time difference.

[0119] In this embodiment S42, the spatial correlation constraint is designed based on the geographical distance between channels and the correlation of seismic events. Specifically, this embodiment maintains a channel distance mapping table (channelDistances) to record the geographical spatial distance between each channel, and sets the spatial correlation threshold to 0.7. The spatial correlation of two initial seismic events... It is obtained by multiplying the tenth intermediate term, the time overlap factor, and the distance weight factor.

[0120] The tenth intermediate term is obtained by dividing the minimum STA / LTA ratio by the maximum STA / LTA ratio between the two initial seismic events;

[0121] The time overlap factor is obtained by dividing the eleventh intermediate term and the twelfth intermediate term. The eleventh intermediate term is the maximum value between 0 and the eighth component. The eighth component is the difference between the minimum earthquake end time and the maximum earthquake start time in two initial earthquake events. The twelfth intermediate term is the minimum value between the ninth component and the tenth component. The ninth component is the difference between the earthquake end time and the earthquake start time of one initial earthquake event. The tenth component is the difference between the earthquake end time and the earthquake start time of another initial earthquake event.

[0122] The distance weighting factor is the maximum value between 0 and the eleventh component. The eleventh component is the difference between 1 and the twelfth component. The twelfth component is obtained by dividing the geographic spatial distance between the channels where the two initial seismic events are located by a preset spatial distance threshold.

[0123] Calculate using the following formula:

[0124]

[0125]

[0126]

[0127] in, These represent the initial earthquake events. and the initial earthquake event The STA / LTA ratio; Indicates the time overlap factor; Indicates the distance weighting factor; and These are the initial earthquake events. The start and end times of the earthquake; and These are the initial earthquake events. The start and end times of the earthquake; Indicates the initial earthquake event Location and initial earthquake event The geographical distance between the passageways; This indicates the preset spatial distance threshold.

[0128] S43. For each candidate earthquake event in the candidate earthquake event list, if the same candidate earthquake event is detected by more than Nc channels, then the candidate earthquake event is added to the final detection list as a fused earthquake event. The earthquake start time and earthquake end time of the fused earthquake event are the earliest earthquake start time and the latest earthquake end time of all corresponding candidate earthquake events, respectively. The STA / LTA ratio of the fused earthquake event is the maximum STA / LTA ratio of all corresponding candidate earthquake events. The confidence of the fused earthquake event is obtained by fusing the confidence of all corresponding candidate earthquake events. Wherein, Nc is a preset channel number threshold.

[0129] Furthermore, the confidence level of the fused seismic events is obtained by dividing the weighted sum of the confidence levels of all candidate seismic events by the sum of all fusion weights; the fusion weight of a candidate seismic event is obtained by multiplying the STA / LTA ratio of the candidate seismic event by the reliability coefficient of the channel and then dividing by a preset normalized benchmark threshold; the reliability coefficient is obtained by dividing the historical effective detection count of the channel by the corrected total historical detection count of the channel, and the corrected total historical detection count of the channel is obtained by adding the total historical detection count of the channel to a smoothing term.

[0130] In this embodiment S43, relevant candidate seismic events are found using the "find Correlated Events" method. Fusion of these candidate seismic events is only performed when the number of candidate seismic events reaches a preset channel number threshold (default 3). During fusion, for each fused seismic event, the trigger count is automatically incremented, and the fused seismic event type is marked as EARTHQUAKE. The confidence level of the fused seismic events is also considered. The calculation uses the following formula:

[0131]

[0132]

[0133]

[0134] in, This indicates the number of candidate seismic events participating in the fusion; Indicates the first The fusion weights of candidate seismic events; Indicates the first Confidence level of each candidate earthquake event; Indicates the first The STA / LTA ratio of candidate earthquake events; The normalized benchmark threshold represents the significance of candidate seismic events. In this embodiment, the normalized benchmark threshold is taken as 2 to 3 times the adaptive detection threshold. Indicates the reliability coefficient of the channel; This indicates the number of historical valid detections for the channel; This indicates the total number of historical tests conducted on the channel; This represents a smoothing term to prevent the denominator from being zero. Each execution of the S4 process is considered a detection, and generating a merged seismic event is considered a valid detection.

[0135] In this embodiment, multi-channel fusion detection employs the performMultiChannelFusion method, supporting a high-confidence single-channel seismic event retention mechanism. When the confidence level of a single-channel seismic event exceeds a threshold of 1.2 times, it is still retained. Furthermore, a concurrent processing architecture is adopted, using ConcurrentHashMap to ensure thread safety, controlling processing latency to the millisecond level, meeting the performance requirements of real-time detection. Seismic event confirmation results are monitored for performance using the DetectionStatistics class, including key indicators such as total number of detections, effective number of detections, fusion detections, effective detection rate, and fusion rate. Further, this embodiment employs a memory pool management mechanism, reducing memory allocation overhead through object reuse and improving processing performance.

[0136] To better demonstrate the specific implementation and technical effects of the present invention, the distributed acoustic sensing seismic monitoring data stream processing method shown in steps S1 to S4 of the above preferred implementation is applied to a specific example below.

[0137] Example

[0138] The specific implementation process of the distributed acoustic sensing seismic monitoring data stream processing method used in this embodiment is as described above and will not be repeated here. A brief description of the implementation process of this embodiment is provided below. Figure 4 As shown, the data stream processing method provided by this invention is a closed-loop real-time monitoring process. The process begins with connecting to the data source and sequentially proceeds through raw seismic monitoring data parsing, data validity verification, channel grouping, sliding window processing, amplitude data conversion, adaptive short-time / long-time average detection algorithm calculation, and multi-channel fusion detection. Finally, it outputs the fused seismic event and determines whether to continue detection. If the data is invalid or an anomaly occurs, this invention automatically records an error log and terminates the current processing flow, ensuring the stable operation of the method.

[0139] In this embodiment, a region with a complex fault network and a significant regional stress field is selected as the target region. This region contains multiple major fault zones and has a rich history of seismic activity. Its complex geological structure and frequent seismic activity provide a sufficient data foundation and practical application scenario for verifying the effectiveness of the method of the present invention for DAS seismic monitoring.

[0140] This embodiment utilizes DAS monitoring data deployed along existing fiber optic lines, with a fiber length of approximately 55 km and a typical channel spacing of about 8-10 m. 4096 channels were selected for processing. The Ridgecrest region DAS data used in this embodiment originates from the Southern California Earthquake Data Center (SCEDC), in SEG-Y standard format, with a sampling frequency of 250 Hz, conforming to SCEDC's official data specifications. The time base is uniformly UTC, and all channel data undergoes verification and integrity checks before entering the processing flow. The fiber optic cables are laid along existing communication infrastructure, covering major fault zones and seismically active areas within the target region, ensuring effective capture of seismic wave signal propagation characteristics.

[0141] To verify the effectiveness of the method of the present invention, the following ablation experiments were conducted in this embodiment. 1) The short-time averaging / long-time averaging detection algorithm adopts the traditional scheme as the basis for performance comparison; 2) The adaptive short-time averaging / long-time averaging detection algorithm, within the framework of the present invention, only enables the adaptive short-time averaging / long-time averaging algorithm without performing multi-channel fusion detection; 3) The "multi-channel fusion detection", within the framework of the present invention, only enables the multi-channel fusion detection process without using the adaptive short-time averaging / long-time averaging detection algorithm; 4) The method of the present invention is the complete scheme executed according to the aforementioned S1~S4 processes. The phase pickup performance of the ablation experiment is comprehensively evaluated using multiple indicators, including detection accuracy, false alarm rate, and time accuracy assessment. The time accuracy assessment uses the root mean square error (RMSE). ) method. Among them, detection accuracy The formula for calculating the root mean square error is as follows:

[0142]

[0143]

[0144] in, This indicates the number of true positives, i.e., the number of correctly detected seismic phases; This indicates the number of false positives, i.e., the number of false alarms. This indicates the number of false negatives, i.e., the number of missed seismic phases. This represents the total number of earthquake events used in the assessment; Indicates the first The detection time for each earthquake event; Indicates the first The actual time of each earthquake event is manually annotated.

[0145] As shown in Table 1, the adaptive short-time averaging / long-time averaging detection algorithm improves the detection rate and reduces the false alarm rate in complex noisy environments by dynamically adjusting parameters; multi-channel fusion detection further reduces the false alarm rate by utilizing the spatial correlation of the original seismic monitoring data to eliminate isolated noise interference. This invention combines two methods to further improve the detection rate and reduce the false alarm rate. Due to the increased computational load, the latency of the method in this invention is increased to some extent, but it is still kept within an acceptable range for practical applications.

[0146] Table 1. Ablation Experiment Results

[0147]

[0148] Figure 5 The two-dimensional wavefield morphology of the original seismic monitoring data after preprocessing steps such as data analysis, invalid data filtering, and bandpass filtering is shown. Figure 5 The horizontal axis represents the channel number of the distributed acoustic sensor (DAS) array, the vertical axis represents time (seconds), and the color scale on the right represents the normalized strain amplitude intensity. Figure 5 The spatially continuous tilted fringe features presented in the image correspond to the P-wave and S-wave wavefront signals propagating along the optical fiber; the red and cyan dashed lines superimposed on the wavefield background are the theoretically calculated P-wave reference path and S-wave reference path, respectively, used to indicate the expected wavefront position for signal feature analysis.

[0149] Figure 6 This is used to illustrate the distribution of deviations between the detection results and theoretical values ​​of the present invention. The horizontal axis represents the channel number of the distributed acoustic sensor (DAS) array, and the vertical axis represents the residual value (in seconds), which is the difference between the automatically picked-up time and the theoretically calculated time. Figure 6 In the diagram, the discrete data points located near the zero baseline reflect the detection bias of each channel. The red circles represent P-wave residuals, and the blue squares represent S-wave residuals. Figure 6 The root mean square value indicated in the legend represents the statistical error level of the final earthquake event detection results output by the method of this invention after performing multi-channel spatial correlation analysis.

[0150] Figure 7The invention demonstrates the automatic detection and phase picking results of earthquake events using the method of this invention. The horizontal axis represents the channel number of the distributed acoustic sensor (DAS) array, corresponding to the spatial position along the optical fiber; the vertical axis represents the relative time (seconds) of the earthquake occurrence. The red solid line and the blue solid line represent the theoretical P-wave arrival time curve and the theoretical S-wave arrival time curve calculated based on the source parameters, respectively. The red circular markers and blue square markers distributed along the curves represent the automatically picked P-wave and S-wave event points generated after adaptive short-time averaging / long-time averaging detection algorithms and multi-channel fusion detection processing, respectively. Figure 7 The values ​​shown in the legend represent the total number of seismic phases that the present invention effectively identified and fused within the current time window.

[0151] Figure 8 In the graph, the horizontal axis represents the number of concurrently accessing DAS channels, covering a test range from 1,000 to 10,000 channels; the vertical axis represents the average end-to-end latency of the invention from receiving earthquake monitoring data packets to outputting earthquake detection results, in milliseconds (ms). The test curves show that as the number of channels increases linearly, the latency of the invention remains stable within a low range of 120ms to 180ms, without an exponential increase. This indicates that the sliding window and parallel computing mechanism employed in this invention can effectively handle large-scale data concurrency, ensuring millisecond-level real-time response capability for earthquake monitoring.

[0152] Figure 9 In the graph, the horizontal axis represents the parallelism of the stream processing job configuration, i.e., the total number of threads on the computing nodes; the vertical axis represents the maximum stable data throughput, expressed in megabytes per second (MB / s). The results show a significant linear positive correlation between throughput and parallelism. When the parallelism is increased from 8 to 32, the throughput of this invention increases by nearly four times, verifying the excellent horizontal scalability of the distributed architecture and its ability to flexibly adapt to DAS monitoring arrays of different sizes by increasing computing resources.

[0153] Figure 10 In the graph, the horizontal axis represents the continuous running time in hours; the vertical axis represents the real-time utilization of the Java Virtual Machine (JVM) heap memory. Figure 10 The curves exhibit regular sawtooth-like fluctuations, clearly reflecting the periodic creation and garbage collection process of memory objects. Crucially, the baseline for memory usage remains horizontal throughout the entire test period, without any drift over time. This demonstrates that the memory pool management and object reuse strategies in this invention effectively eliminate the risk of memory leaks and ensure the stable operation of this invention in long-term unattended scenarios.

[0154] It should also be noted that the distributed acoustic sensing seismic monitoring data stream processing method in the above embodiments can essentially be executed by a computer program or module. Therefore, similarly, based on the same inventive concept, another preferred embodiment of the present invention also provides a distributed acoustic sensing seismic monitoring data stream processing system corresponding to the distributed acoustic sensing seismic monitoring data stream processing method provided in the above embodiments, such as... Figure 11 As shown, it includes:

[0155] The data acquisition module is used to receive raw seismic monitoring data collected by distributed acoustic sensors through a data source;

[0156] The result acquisition module is used to process the raw earthquake monitoring data according to the distributed acoustic sensing earthquake monitoring data stream processing method described in the above embodiments to obtain the final detection list of earthquake events.

[0157] It is understood that the distributed acoustic sensing seismic monitoring data stream processing method described in S1-S4 above can essentially be implemented by a computer program. Therefore, based on the same inventive concept, another preferred embodiment of the present invention also provides a computer program product corresponding to the distributed acoustic sensing seismic monitoring data stream processing method provided in the above embodiments, which includes a computer program / instructions. When the computer program / instructions are executed by a processor, they can implement the distributed acoustic sensing seismic monitoring data stream processing method as described in the above embodiments.

[0158] It should also be noted that those skilled in the art will understand that, for the sake of convenience and brevity, the specific working process of the system described above can be referred to the corresponding process in the foregoing method embodiments, and will not be repeated here. In the embodiments provided in this application, the division of steps or modules in the system and method is merely a logical functional division, and there may be other division methods in actual implementation. For example, multiple modules or steps may be combined or integrated together, and a module or step may also be split.

[0159] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the invention. Those skilled in the art can make various changes and modifications without departing from the spirit and scope of the invention. Therefore, all technical solutions obtained through equivalent substitution or transformation fall within the protection scope of the present invention.

Claims

1. A method for processing distributed acoustic sensing seismic monitoring data streams, characterized in that, Includes the following steps: S1. Receive raw seismic monitoring data collected by distributed acoustic sensors through the data source, and use a distributed acoustic sensor data parser to parse and preprocess the raw seismic monitoring data; S2. Apply a sliding processing time window to slice the preprocessed seismic monitoring data, and use the seismic monitoring data in each window as a data segment suitable for real-time detection; S3. In the seismic detection processor, an adaptive short-time average / long-time average detection algorithm is executed on the data segments, and the adaptive parameters in the algorithm are dynamically adjusted according to the characteristics of the seismic monitoring data during the detection process, and an initial candidate list of seismic events is output. S4. Perform multi-channel fusion detection on the initial candidate list of earthquake events to generate the final detection list of earthquake events; In step S3, the adaptive parameters include short-term average window length, long-term average window length, and detection trigger threshold; In step S3, the short-time average window length is obtained by limiting the calculated first base value to between the first minimum value and the first maximum value through a constraint function. The first maximum value is twice the sampling rate. The first base value is obtained by multiplying the preset first constant, the first intermediate term, and the second intermediate term. The first intermediate term is the ratio of the sampling rate to the first component. The first component is the maximum value between the main frequency and 1. The second intermediate term is obtained by adding the second component to 1. The second component is obtained by multiplying the signal variability index and the second constant. The long-term average window length is obtained by limiting the calculated second base value to between the second minimum and the second maximum value through a constraint function. The second minimum value is an integer multiple of the short-term average window length, and the second maximum value is an integer multiple of the sampling rate. The second fundamental value is obtained by multiplying the third, fourth, and fifth intermediate terms. The third intermediate term is an integer multiple of the short-time average window length. The fourth intermediate term is obtained by adding the third component to 1. The third component is a multiple of the noise level. The fifth intermediate term is the difference between the third constant and the signal stationarity index. The detection trigger threshold is obtained by limiting the calculated third basic value to between the third minimum value and the third maximum value through a constraint function, where the third minimum value and the third maximum value are both constants; The third fundamental value is obtained by multiplying the fourth constant, the sixth intermediate term, the seventh intermediate term, and the eighth intermediate term. The sixth intermediate term is obtained by adding the fourth component to 1. The fourth component is a multiple of the noise level. The seventh intermediate term is calculated based on the signal-to-noise ratio. The eighth intermediate term is obtained by adding the fifth component to 1. The fifth component is a multiple of the signal variability index.

2. The distributed acoustic sensing seismic monitoring data stream processing method as described in claim 1, characterized in that, Noise levels, frequency distribution of earthquake monitoring data, signal variability indicators, and signal stationarity indicators were all obtained by analyzing earthquake monitoring data using a signal feature analyzer.

3. The distributed acoustic sensing seismic monitoring data stream processing method as described in claim 1, characterized in that, In step S3, for each channel's data segment, the average signal energy within the short-time window and the average signal energy within the long-time window are calculated, and the ratio of the two average signal energy values ​​is taken as the STA / LTA ratio. When the STA / LTA ratio is greater than the adaptive detection threshold, it is considered that a seismic event has occurred in that channel and it is added to the initial candidate list of seismic events.

4. The distributed acoustic sensing seismic monitoring data stream processing method as described in claim 3, characterized in that, The adaptive detection threshold is obtained by multiplying the calculated detection trigger threshold by a correction term, where the correction term is... The result is obtained by adding the ninth intermediate term, which is the product of the sixth component and the preset adjustment sensitivity factor. The sixth component is the result of taking the logarithm of the seventh component. The seventh component is the ratio of two signal-to-noise ratios. The first signal-to-noise ratio is the real-time signal-to-noise ratio calculated from the data segment of the current channel, and the second signal-to-noise ratio is the preset reference signal-to-noise ratio.

5. The distributed acoustic sensing seismic monitoring data stream processing method as described in claim 1, characterized in that, In step S3, the adaptive short-time average / long-time average detection algorithm adopts a parameter caching mechanism: maintaining a parameter cache mapping table to store the adaptive parameters of each channel; setting a cache expiration time, and when the cache expiration time is reached, recalculating the adaptive parameters and updating the parameter cache mapping table.

6. The distributed acoustic sensing seismic monitoring data stream processing method as described in claim 1, characterized in that, In step S4, the specific process of performing multi-channel fusion detection on the initial candidate list is as follows: S41. Extract the earthquake events detected by each channel from the initial candidate list of earthquake events and use them as initial earthquake events. Each initial earthquake event contains six key pieces of information, namely, initial earthquake event ID, channel ID, earthquake start time, earthquake end time, STA / LTA ratio, and confidence level. S42. When two initial earthquake events simultaneously satisfy both the temporal correlation constraint and the spatial correlation constraint, the two initial earthquake events are added to the candidate earthquake event list; wherein, the temporal correlation constraint is that the trigger time difference between the two initial earthquake events is less than or equal to the maximum allowed trigger time difference, and the trigger time difference is the difference in the earthquake start time of the two initial earthquake events; the spatial correlation constraint is that the spatial correlation between the two initial earthquake events is greater than the preset spatial correlation threshold. S43. For each candidate earthquake event in the candidate earthquake event list, if the same candidate earthquake event is detected by more than Nc channels, then the candidate earthquake event is added to the final detection list as a fused earthquake event. The earthquake start time and earthquake end time of the fused earthquake event are the earliest earthquake start time and the latest earthquake end time of all corresponding candidate earthquake events, respectively. The STA / LTA ratio of the fused earthquake event is the maximum STA / LTA ratio of all corresponding candidate earthquake events. The confidence of the fused earthquake event is obtained by fusing the confidence of all corresponding candidate earthquake events. Wherein, Nc is a preset channel number threshold.

7. The distributed acoustic sensing seismic monitoring data stream processing method as described in claim 6, characterized in that, In step S42, the spatial correlation between the two initial seismic events is obtained by multiplying the tenth intermediate term, the time overlap factor, and the distance weighting factor. The tenth intermediate term is obtained by dividing the minimum STA / LTA ratio by the maximum STA / LTA ratio between the two initial seismic events; The time overlap factor is obtained by dividing the eleventh intermediate term and the twelfth intermediate term. The eleventh intermediate term is the maximum value between 0 and the eighth component. The eighth component is the difference between the minimum earthquake end time and the maximum earthquake start time in two initial earthquake events. The twelfth intermediate term is the minimum value between the ninth component and the tenth component. The ninth component is the difference between the earthquake end time and the earthquake start time of one initial earthquake event. The tenth component is the difference between the earthquake end time and the earthquake start time of another initial earthquake event. The distance weighting factor is the maximum value between 0 and the eleventh component. The eleventh component is the difference between 1 and the twelfth component. The twelfth component is obtained by dividing the geographic spatial distance between the channels where the two initial seismic events are located by a preset spatial distance threshold.

8. The distributed acoustic sensing seismic monitoring data stream processing method as described in claim 6, characterized in that, In step S43, the confidence level of the fused seismic events is obtained by dividing the weighted sum of the confidence levels of all candidate seismic events by the sum of all fusion weights. The fusion weight of a candidate seismic event is obtained by multiplying the STA / LTA ratio of the candidate seismic event by the reliability coefficient of the channel and then dividing by a preset normalized benchmark threshold. The reliability coefficient is obtained by dividing the number of historical valid detections of the channel by the corrected total number of historical detections of the channel. The corrected total number of historical detections of the channel is obtained by adding the total number of historical detections of the channel to a smoothing term.

9. A distributed acoustic sensing seismic monitoring data stream processing system, characterized in that, include: The data acquisition module is used to receive raw seismic monitoring data collected by distributed acoustic sensors through a data source; The result acquisition module is used to process the original earthquake monitoring data according to the distributed acoustic sensing earthquake monitoring data stream processing method according to any one of claims 1 to 8, and obtain the final detection list of earthquake events.

Images