Dense array earthquake real-time monitoring and risk analysis method and system

Abstract

The application discloses a kind of intensive array seismic real-time monitoring and risk analysis method and system, belong to energy development production field.The method includes: the original seismic data of real-time acquisition is carried out format conversion and segmentation processing;And based on the seismic data after pre-processing, in turn, seismic phase picking, earthquake association, earthquake relocation and magnitude calculation, finally by two kinds of evaluation methods based on the four-stage results of the above earthquake risk quantification evaluation, finally the evaluation result is visualized and fed back to user end.The application realizes the automation from data acquisition, real-time processing to risk analysis, significantly improves the timeliness of induced seismic monitoring and risk assessment reliability.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of energy development and production, and particularly relates to a method and system for real-time monitoring and risk analysis of intensive array earthquakes. BACKGROUND

[0002] Exploration and development of shale gas is of great significance to energy strategy. "Hydraulic fracturing" (referred to as "fracturing") is the main means of unconventional shale gas exploitation. However, during drilling, large-scale loss of drilling fluid may occur when faults are encountered, and stress disturbance may occur around the well during fracturing and other operations, which may induce seismic activity. Although the earthquake is usually small in size, it may still have an impact on the surrounding environment. In addition, once a medium-strong induced earthquake occurs, not only may it damage nearby infrastructure, affect production, and cause economic losses, but it may also cause casualties and arouse community concerns about seismic risk. How to monitor seismic risk during production and adjust production plans in a timely manner to avoid potentially destructive induced earthquakes has become an increasingly important issue for the oil and gas industry and related departments.

[0003] In shale gas production and development, a large amount of fluid is injected into the underground reservoir, changing the stress state of the underground space and disrupting the stability of the fault, thereby inducing or triggering the occurrence of an earthquake. According to statistical data, in the past two decades, there have been more than 110 destructive earthquakes caused by water injection operations. Therefore, before and after the development of shale gas production, accurate monitoring methods and scientific production guidance are needed to avoid large earthquakes.

[0004] Currently, the main problems in the monitoring task of induced earthquakes can be summarized as follows:

[0005] First, there is a lag in the processing of collected data and the output of real-time monitoring results, making it difficult to grasp the seismic activity situation in the monitoring area in the first time, resulting in insufficient earthquake risk warning and avoidance capability, and thus failing to achieve timely analysis and warning of seismic risk.

[0006] Second, the current induced earthquake processing flow has significant instability, and the monitoring scheme is difficult to be directly reused in different regions and different equipment systems, and the difference in flow adaptability further reduces the reliability of the final result.

[0007] Third, there is a lack of sufficient auxiliary information to support risk analysis in real-time seismic risk monitoring. Currently, real-time risk analysis relies only on seismic frequency and magnitude for qualitative evaluation, making it difficult to achieve objective and quantitative judgment of seismic risk. SUMMARY

[0008] In view of the short board existing in the induced earthquake monitoring, the present application aims to solve two core problems: one is to improve the real-time performance of induced earthquake monitoring and the reliability of monitoring results; the other is to solve the single index of earthquake risk evaluation, so as to realize the quantitative evaluation of risk. In order to solve the above problems, the present application provides a dense array earthquake real-time monitoring and risk analysis method and system.

[0009] In order to achieve the above-mentioned application purposes, the present application specifically adopts the following technical solutions:

[0010] In the first aspect, the present application provides a dense array earthquake real-time monitoring and risk analysis method, which comprises the following steps:

[0011] S1. Real-time acquisition of original seismic data continuously monitored by a dense array, and format conversion and segmentation processing of the original seismic data to form preprocessed seismic data;

[0012] S2. Inputting the preprocessed seismic data into a machine learning-based seismic phase picking model to output the seismic phase picking result corresponding to the seismic data;

[0013] S3. Associating discrete seismic phases in the seismic phase picking result into a specific seismic event by using a fast seismic correlation and positioning method, outputting the seismic correlation result, and returning to step S1 directly when not associated with a seismic event;

[0014] S4. Positioning the seismic correlation result by a two-step positioning method, obtaining an accurate spatial position of the seismic source through independent inversion of the seismic plane position and the depth position, and outputting the final positioning result;

[0015] S5. Reading the seismic event information from the final positioning result, calculating the magnitude of each seismic station, and taking the median of the magnitudes of all seismic stations as the magnitude result of the seismic event;

[0016] S6. Based on the final positioning result and the magnitude result of the seismic event, using a first evaluation method based on real-time dynamic b value and a second evaluation method based on the distribution of P wave velocity and S wave velocity ratio for seismic risk quantitative evaluation; in the first evaluation method, the b value is negatively correlated with the seismic event risk, and the value is cross-validated by the maximum likelihood estimation method and the b+ method; in the second evaluation method, high-precision seismic event distribution is obtained by the double-difference positioning method, P wave travel time difference and S wave travel time difference are calculated by the waveform cross-correlation method, and the ratio of P wave velocity to S wave velocity is estimated, and the estimated value of the ratio of P wave velocity to S wave velocity is positively correlated with the fluid influence.

[0017] On the basis of the above scheme, each step can be realized in the following preferred specific manner.

[0018] As a preferred embodiment of the first aspect, in step S1, the original seismic data is converted into SAC format, and the seismic data after format conversion is divided into data records with a length of 1 minute as preprocessed seismic data; when detecting a seismic data transmission interruption, firstly, the difference between the seismic data acquisition time and the real world time is detected: if the difference is less than a preset time length threshold, it is determined that the current seismic station is in real-time processing state, and no backtracking processing is performed; if the difference is greater than or equal to the time length threshold, the previous interruption time is read from the record log, the seismic data corresponding to the previous interruption time is taken as a historical data node, and backtracking processing is performed from the historical data node until the difference is less than the time length threshold.

[0019] As a preferred embodiment of the first aspect, in step S2, the phase picking model uses PhaseNet, and PhaseNet outputs the phase picking results with a probability higher than a preset phase reservation threshold in a CSV file format, wherein the phase picking results include the seismic station ID, the seismic data acquisition time, the number of phase picking sampling points, the phase arrival time, the phase picking probability value, the amplitude value at the picked phase, the phase type, and the seismic data storage path; the phase type is a P-wave phase or an S-wave phase.

[0020] As a preferred embodiment of the first aspect, in the phase association process of step S3, in addition to reading the phase picking results of the current preprocessed seismic data, the phase picking results of the last Ttravel seconds in the previous preprocessed seismic data are also read synchronously, wherein Ttravel represents the maximum value of the seismic wave travel time calculated based on the farthest position of the possible occurrence of the seismic event and the seismic wave velocity.

[0021] As a preferred embodiment of the first aspect, in step S4, the specific process of positioning using the two-step positioning method is as follows:

[0022] S41. Convert the phase association result file into a data input format suitable for the Hypoinverse algorithm or the VELEST algorithm;

[0023] S42. After format conversion, use the Hypoinverse algorithm or the VELEST algorithm to preliminarily invert the planar position of the seismic event source;

[0024] S43. According to a preset epicenter distance threshold, only the phase association results within the epicenter distance threshold are selected and read, and the Hypoinverse algorithm or the VELEST algorithm is used to perform positioning calculation based on the selected phase association results to obtain the depth position of the seismic event source, and the planar position and the depth position together constitute the spatial position of the source.

[0025] As a preferred embodiment of the first aspect, in step S6, the process of using the second evaluation method for seismic risk quantitative evaluation is as follows:

[0026] S61. Extract the origin time and spatial location of each seismic event from the final positioning result as the earthquake positioning data, obtain the longitude, latitude and station name of the seismic station as the seismic station data, obtain the earthquake catalog containing all seismic events, obtain the travel time difference data by screening the seismic event pairs in the earthquake catalog, and store the three kinds of data in the form of a dictionary;

[0027] S62. For each seismic event pair, screen the stations that record P-wave travel time difference and S-wave travel time difference at the same time to obtain a station list and corresponding P-wave travel time difference and S-wave travel time difference, and perform mean removal operation on the two kinds of travel time difference respectively to obtain the relative change amount of P-wave actual travel time difference and the relative change amount of S-wave actual travel time difference;

[0028] S63. The random sample consensus algorithm is used to coarsely screen the relative change amounts of the two kinds of actual travel time difference, and the coarsely screened relative change amounts of the two kinds of actual travel time difference are recorded as the first change amount and the second change amount respectively; a data point is formed by the first change amount and the second change amount, the perpendicular distance from the data point to the known straight line and the mean of all distances are calculated, the distance mean is multiplied by a preset hyperparameter to obtain a distance retention threshold, and the data points with a distance less than the distance retention threshold are retained; wherein, the known straight line is a corresponding straight line with a theoretical P-wave velocity to theoretical S-wave velocity ratio equal to 1.732;

[0029] S64. The P-wave velocity and S-wave velocity of the seismic event are preset, the straight line distance between the hypocenter and the station is obtained based on the final positioning result and the station coordinates, and the P-wave theoretical travel time is obtained by dividing the straight line distance by the P-wave velocity, and the S-wave theoretical travel time is obtained by dividing the straight line distance by the S-wave velocity;

[0030] S65. The P-wave theoretical travel time difference and S-wave theoretical travel time difference of each seismic event pair are calculated according to the P-wave theoretical travel time and S-wave theoretical travel time of the seismic event, and mean removal operation is performed on the two kinds of theoretical travel time difference respectively to obtain the relative change amount of P-wave theoretical travel time difference and the relative change amount of S-wave theoretical travel time difference;

[0031] S66. Least squares method is used to fit based on the screened data points to obtain the preliminary estimate value of the P-wave velocity to S-wave velocity ratio and the preliminary estimate value of the intercept, and the total least squares method is combined with the two preliminary estimate values and the screened data points to re-fit in a constantly iterative manner, and the difference between the theoretical travel time difference relative change amount and the actual travel time difference relative change amount is calculated as a residual during the iteration process, and the final estimate value of the P-wave velocity to S-wave velocity ratio is obtained after the iteration is completed.

[0032] As a preferred form of the first aspect, in S61, the travel time difference data is obtained as follows: a pair of seismic events is formed by two different seismic events in the earthquake catalog, the travel time difference of each pair of seismic events is calculated to form a travel time difference file, target seismic events are selected from the earthquake location data according to the latitude and longitude of the target region to form a target event seismic ID list, and a pair of seismic events that satisfies the following conditions is selected: both seismic events of the pair of seismic events are in the target event seismic ID list, the pair of seismic events has at least S seismic stations simultaneously recording P-wave travel time difference and S-wave travel time difference, and the travel time difference of the selected pair of seismic events is used as the travel time difference data; wherein S is a preset seismic array quantity threshold.

[0033] In a second aspect, the present application provides a dense array seismic real-time monitoring and risk analysis system, which comprises:

[0034] a data acquisition module configured to acquire original seismic data within a user-set geographical space range and time range;

[0035] a result acquisition module configured to process the original seismic data according to the dense array seismic real-time monitoring and risk analysis method of any one of the first aspect, calculate the b value, and estimate the in-situ ratio of P-wave velocity to S-wave velocity of the source region;

[0036] a visualization module configured to send a warning information or a structured risk assessment report to the user terminal in a preset manner when the calculated b value deviates from the expected value or exceeds the preset risk threshold.

[0037] In a third aspect, the present application provides a computer readable storage medium, wherein the storage medium stores a computer program, and when the computer program is executed by a processor, the dense array seismic real-time monitoring and risk analysis method of any one of the first aspect is realized.

[0038] In a fourth aspect, the present application provides a computer electronic device comprising a memory and a processor;

[0039] the memory is configured to store a computer program;

[0040] the processor is configured to realize the dense array seismic real-time monitoring and risk analysis method of any one of the first aspect when the computer program is executed.

[0041] Compared with the prior art, the present application has the following beneficial effects:

[0042] The method constructs a set of induced earthquake real-time dynamic monitoring processes suitable for the oil and gas field production field, improves the traditional artificial and step-by-step processing non-continuous monitoring mode, and achieves the effect of full-process automation and real-time operation, which significantly improves the timeliness of induced earthquake monitoring.

[0043] The method first realizes efficient acquisition of real-time monitoring results: using the Internet of Things data acquisition equipment to upload seismic data to the server through the 4G network, combining the Internet remote data download function, breaking the geographical limitations of data processing, realizing real-time processing and analysis of data.

[0044] Secondly, the method adopts the PhaseNet method based on machine learning for phase picking, improves the traditional manual picking or automatic picking method based on fixed threshold, and achieves the effect of greatly improving the picking efficiency and accuracy, and adapting to the real-time processing demand of continuous data stream.

[0045] In the earthquake association link, the cross-period data reading mechanism (reading the last Ttravel seconds of the phase in the previous minute data) is introduced, the defect that the phase may be missed due to data segmentation cutting is improved, and the effect of ensuring all phases associated with the earthquake event and improving the completeness of earthquake event detection is achieved.

[0046] In the positioning link, the "two-step positioning method" (first inversion of the plane position, and then optimization of the depth based on the epicenter distance threshold) is proposed, the problem of poor depth parameter accuracy in the conventional positioning algorithm is improved, and the effect of obtaining high-precision plane position and depth position respectively and improving the overall accuracy of the source spatial positioning is achieved.

[0047] In the magnitude calculation, the median of the single station magnitude of each station is used as the final magnitude, the algorithm of directly taking the average value is improved, and the effect of effectively suppressing the interference of abnormal values of individual stations and improving the robustness of the magnitude result is achieved.

[0048] In the earthquake risk quantitative evaluation, the cross-validation mechanism based on real-time dynamic b value and b+ method is introduced, the limitation of single dependence on maximum likelihood estimation method is improved, the effect of reducing the probability of earthquake risk misjudgment and improving the reliability of risk evaluation result is achieved; and in the analysis of the ratio of P wave velocity and S wave velocity, the multi-round abnormal value suppression strategy of "coarse screening + fine screening" is adopted, the limitation of single denoising method is improved, and the effect of effectively suppressing noise and abnormal values in the observation travel time difference and improving the calculation accuracy of the ratio is achieved.

[0049] In addition, the real-time dynamic display system matched with the process enables the user to grasp the latest earthquake distribution and earthquake risk index at any time and any place, significantly simplifies the communication process, improves the communication efficiency, and provides a solid guarantee for the safety of energy development. BRIEF DESCRIPTION OF DRAWINGS

[0050] Figure 1 is a flowchart of the method of the present application;

[0051] Figure 2 is a system block diagram of the present application;

[0052] Figure 3 is a real-time earthquake monitoring array distribution map provided by the present embodiment;

[0053] Figure 4 is a real-time dynamic earthquake event distribution map provided by the present embodiment;

[0054] Figure 5 is a seismic event dynamic b value calculation result map provided by the present embodiment;

[0055] Figure 6 is a schematic diagram of a computer electronic device according to the present application. DETAILED DESCRIPTION

[0056] In order to make the above-mentioned purposes, features and advantages of the present application more obvious and easy to understand, the specific embodiments of the present application will be described in detail below with reference to the drawings. In the following description, a large number of specific details are set forth in order to facilitate a full understanding of the present application. However, the present application can be implemented in many other ways different from those described herein, and those skilled in the art can make similar improvements without departing from the spirit of the present application, so the present application is not limited to the specific embodiments disclosed below. The technical features in each embodiment of the present application can be combined accordingly without conflict.

[0057] In the description of the present application, it should be understood that the terms "first", "second" are only used for distinguishing description purposes and cannot be understood as indicating or implying relative importance or implicitly indicating the number of the technical features indicated. Therefore, the features defined with "first", "second" can explicitly or implicitly include at least one of the features.

[0058] As shown in Figure 1 In a preferred implementation of the present application, the above-mentioned dense array real-time earthquake monitoring and risk analysis method includes the following S1-S6 steps. The specific implementation process will be described in detail below.

[0059] I. Real-time acquisition and processing of earthquake data stream

[0060] S1. Real-time acquisition of original seismic data continuously monitored by dense station array, and format conversion and segmentation of original seismic data to form preprocessed seismic data.

[0061] It should be noted that in step S1 of the present application, the original seismic data is converted into SAC format, and the seismic data after format conversion is cut into 1-minute-long data records as preprocessed seismic data; when a transmission interruption of seismic data is detected, first detect the difference between the seismic data acquisition time and the real world time: if the difference is less than a preset time length threshold, it is determined that the current seismic station is in real-time processing state, and no backtracking processing is performed; if the difference is greater than or equal to the time length threshold, the previous interruption time is read from the log, the seismic data corresponding to the previous interruption time is taken as a historical data node, and the backtracking processing is started from the historical data node until the difference is less than the time length threshold.

[0062] In step S1 of the present embodiment, the real-time continuous original seismic data is first converted into SAC seismic data format, and is simultaneously cut into 1-minute-long data records. The advantage of using SAC format is that each component of the seismic station corresponds to an independent data file, which facilitates subsequent data reading operations and can quickly determine whether the seismic station is in normal recording state; and the data length of 1 minute can balance the efficiency and real-time performance while ensuring data processing efficiency.

[0063] Considering the volatility of 4G data transmission, possible failure of monitoring equipment, and possible computer crash, etc., these problems may cause data transmission delay or even interruption. Therefore, the present application specially designs a seismic data detection algorithm: when an unexpected situation causes the seismic data transmission to stop or the program to be unexpectedly interrupted, after restarting the real-time seismic processing program, the difference between the seismic data acquisition time and the real world time will be first detected, and the above-mentioned time length threshold will be preset in advance. In the present embodiment, the threshold is set to 1 minute. When the difference is less than 1 minute, it is determined that the current seismic station is in "real-time processing state"; when the difference is greater than or equal to 1 minute, the previous interruption time of the program is read from the log, the seismic data corresponding to the previous interruption time is taken as a historical data node, and the backtracking processing is started from the historical data node until the difference between the seismic data acquisition time and the real world time is less than 1 minute. Thereafter, the program will automatically detect the time difference every 5 seconds, and if the difference is again greater than 1 minute, the next round of backtracking processing will be started, and through the cyclic mechanism, the continuity of real-time data processing is ensured, and a stable real-time processing process is realized.

[0064] II. Real-time phase picking based on machine learning

[0065] S2. input the preprocessed seismic data into the machine learning-based phase picking model, and output the phase picking result corresponding to the seismic data.

[0066] It should be noted that in step S2 of the present application, the phase picking model adopts PhaseNet, which outputs the phase picking result with a probability higher than a preset phase retention threshold in the form of a CSV file. The phase picking result includes the seismic station ID, the seismic data acquisition time, the number of phase picking sampling points, the phase arrival time, the phase picking probability value, the amplitude value at the picked phase, the phase type, and the seismic data storage path. The phase type is P-wave phase or S-wave phase.

[0067] In step S2 of the present application, the above-mentioned phase picking model can adopt PhaseNet, or other machine learning-based phase picking methods, and is not limited to PhaseNet. In this embodiment, PhaseNet based on machine learning is selected as the phase picking tool to replace the traditional manual picking mode, which greatly improves the phase picking efficiency and effectively guarantees the picking accuracy.

[0068] To ensure the time sequence consistency of phase picking and data processing, the program automatically records the seismic data acquisition time when the 1-minute data record processing is completed, and passes it to the phase picking model in the form of a parameter. After receiving the corresponding parameter, the model reads the seismic data in the matching time period and performs phase picking operation.

[0069] PhaseNet takes the preprocessed seismic data (three-component seismic waveform) as input. In this embodiment, the phase retention threshold is set to 0.6. PhaseNet automatically selects the phase picking result with a probability higher than the threshold and saves it as a file with the default file name "Picks.csv". Considering that the output file of PhaseNet is named "Picks.csv" by default each time, to avoid data loss caused by subsequent file overwriting, the program copies the current "Picks.csv" file to the specified save directory after phase picking is completed, and renames the file according to the current seismic data acquisition time (for example, "20250101_010100.csv"). The CSV file records the phase picking result of the corresponding time seismic data, including the following core fields: seismic station ID, seismic data acquisition time, number of phase picking sampling points, phase arrival time, phase picking probability value, amplitude value at the picked phase, phase type (P-wave phase or S-wave phase), and seismic data storage path.

[0070] III. Real-time seismic association

[0071] S3. The discrete earthquakes are associated into a specific earthquake event by using a rapid earthquake association and location method, and an earthquake association result is outputted, and the step S1 is directly returned when not associated with an earthquake event.

[0072] It should be noted that in the earthquake association process of step S3 of the present application, in addition to reading the phase picking result of the current preprocessed seismic data, the phase picking result of the last Ttravel seconds in the previous preprocessed seismic data is also read synchronously, where Ttravel represents the maximum seismic wave travel time calculated based on the farthest possible location of the earthquake event and the seismic wave velocity.

[0073] In step S3 of the present application, earthquake association is a core process of matching and associating the seismic phases recorded by different seismic stations to determine whether these seismic phases are derived from the same earthquake event, thereby laying a foundation for determining key parameters such as earthquake location.

[0074] After obtaining the phase picking result by PhaseNet or other machine learning-based phase picking method, the present application uses a rapid earthquake association and location method (REAL) to associate these discrete earthquakes into a specific earthquake event. Consistent with the logic of step S2, after the phase picking process is completed, the same time parameter is applied in the earthquake association process, and the REAL algorithm also automatically reads the phase picking result of the corresponding time period according to the time information in the parameter.

[0075] It should be noted that in the real-time processing flow, 1 minute of data records are read each time, so the seismic phases may be distributed across time periods, i.e., part of the seismic phases fall in the previous minute of data records, and the other part of the seismic phases fall in the current minute of data records, which may easily lead to missing of seismic phases. To solve this problem, the present application first calculates the maximum seismic wave travel time Ttravel by dividing the farthest possible location of the earthquake event by the seismic wave velocity, and when reading the data, the phase picking result of the current preprocessed seismic data and the phase picking result of the last Ttravel seconds in the previous preprocessed seismic data are read synchronously, so as to ensure that all seismic phases of the earthquake event can be associated and information loss caused by time period cutting is avoided.

[0076] Since the present flow processes continuous seismic data, and earthquake events occur intermittently, in most cases, the current 1 minute of data records does not contain an earthquake event. For this purpose, the present application designs a trigger-type execution logic: if the earthquake event is successfully associated, the earthquake association result phase.txt is automatically outputted so as to continue the subsequent process; if the earthquake event is not successfully associated, the step S1 is directly returned, and the processing cycle of the next round of seismic data is entered.

[0077] In addition, after the calculation is completed by the REAL method, the method can not only realize the effective association of seismic phases, but also preliminarily determine the focal position of the seismic event. Although the preliminary positioning result has certain error, it can meet the input data requirement of the subsequent relocation process. The core advantage of using the REAL method is that the calculation efficiency of real-time processing is ensured, and the relatively reliable seismic event positioning result with high precision is quickly obtained, so as to realize the balance between efficiency and positioning precision.

[0078] Four, real-time positioning of seismic events

[0079] S4. The associated results of the seismic phases are positioned by the two-step positioning method, the precise focal spatial position is obtained through independent inversion of the seismic plane position and the depth position, and the final positioning result is output.

[0080] It should be noted that in S4 of the present application, the specific process of positioning by the two-step positioning method is as follows:

[0081] S41. The associated results of the seismic phases are converted into a data input format suitable for the Hypoinverse algorithm or the VELEST algorithm.

[0082] In S41 of the present embodiment, when it is detected that the phase.txt file exists, the real-time positioning process is started. The first step of the process is data format adaptation: that is, the phase.txt file is converted into a data input format suitable for the Hypoinverse algorithm or the VELEST algorithm.

[0083] Further, after the format conversion is completed, the original phase.txt file is deleted in the present embodiment, so as to avoid misjudgment caused by the next iteration processing. In addition, the phase.txt can be saved as a copy file named by the processing time according to actual needs, so as to be used in subsequent other scenes.

[0084] S42. After the format conversion is completed, the plane position of the seismic event focal is preliminarily inverted by using the Hypoinverse algorithm or the VELEST algorithm.

[0085] S43. According to the preset epicenter distance threshold, only the associated results of the seismic phases within the epicenter distance threshold are selected and read, the Hypoinverse algorithm or the VELEST algorithm is used to reposition the selected associated results of the seismic phases, the depth position of the seismic event focal is obtained, and the plane position and the depth position together constitute the spatial position of the focal.

[0086] In the embodiment S42-S43, the application innovatively proposes a "two-step positioning method", that is, through independent inversion of the seismic plane position and the depth position, the precise determination of the spatial position of the seismic source is realized, and the specific process is as follows:

[0087] First step: preliminary inversion of plane position

[0088] Read all the earthquake associated results, use Hypoinverse algorithm or VELEST algorithm for preliminary inversion, and get the high-precision plane position of the earthquake source, that is, the longitude and latitude of the earthquake source.

[0089] Second step: depth information optimization calculation

[0090] According to the preset epicenter distance threshold, only the earthquake associated results in the threshold range are selected and read (that is, the earthquake associated results with too far epicenter distance and large interference on depth calculation are removed), and the Hypoinverse algorithm or VELEST algorithm is used for positioning calculation again, and the depth position with higher precision, that is, the depth of the earthquake source, can be obtained.

[0091] Therefore, the application obtains high-precision plane position and depth position through "two-step positioning", overcomes the problem that the conventional earthquake positioning algorithm can output high-precision plane position (longitude and latitude), but the depth result of the earthquake event may have large error due to too far epicenter distance of some stations, and finally the accurate spatial position of the earthquake event can be determined, and the final positioning result file phase.pha is generated, further improving the accuracy of the earthquake source positioning.

[0092] Five, real-time magnitude update of earthquake event

[0093] S5. After reading the earthquake event information from the final positioning result and calculating the magnitude of each seismic station, the median of the magnitudes of all seismic stations is taken as the magnitude result of the earthquake event.

[0094] In step S5 of the application, the magnitude of the earthquake is the most critical parameter in earthquake monitoring, which determines whether the construction plan can continue or needs to be suspended. Therefore, an accurate magnitude calculation method is particularly important. In this embodiment, first, the relevant information of the earthquake event is read from the final positioning result file phase.pha, including the longitude, latitude, depth, occurrence time of the earthquake, and the S-wave phase arrival time of the corresponding seismic station. When starting to calculate the magnitude, this embodiment will intercept the seismic data of the corresponding seismic station according to the occurrence time of the earthquake event and the S-wave phase arrival time. Since the amplitude physical quantity and unit recorded by different types of stations are different, this paper takes the common velocity type station as an example, whose unit is millivolt (mV), to illustrate the calculation process of the magnitude.

[0095] Firstly, the embodiment can exclude the influence of abnormal records and perform high-pass filtering on the original seismic data at 1 Hz to remove the interference of ultra-low frequency noise. Then, the amplitude unit of the filtered seismic data is converted from millivolt (mV) to volt (V), and the seismic data after unit conversion is subjected to spike removal and mean value processing. Then, the instrument response in the processed seismic data is removed, the velocity record is converted into a displacement record to measure the displacement of the ground surface when the seismic event occurs, the maximum displacement value of the horizontal component of the seismic station is measured, the distance between the seismic station and the seismic source is calculated according to the spatial position of the seismic source and the maximum displacement value, the magnitude of the corresponding seismic station is obtained by using the magnitude calculation formula (for reference GB 17740-2017 "Provisions for Seismic Magnitude" or the magnitude calculation formula corrected according to the region), and the median of the magnitudes of all seismic stations is taken as the magnitude result of the seismic event , wherein, represents the magnitude of the i-th seismic station, represents the seismic station ID.

[0096] Therefore, the present application designs a complete real-time monitoring process of dynamic induced seismic data by S1-S5, which is the core of the present application. The purpose of the process is to timely invert the occurrence time, source location and magnitude of the seismic event. Specifically, the method of the present application first identifies all possible seismic phases from the continuously monitored seismic data waveform, then performs correlation calculation on the detected seismic phases, thereby obtaining the true seismic phase record corresponding to the seismic event, and further obtaining the high-precision positioning result and updating the magnitude of the seismic event through positioning calculation.

[0097] It should be noted that the automatic picking of seismic phases, seismic correlation and seismic positioning method involved in the above process are not the only scheme - existing research has proposed various similar technologies, and in the real-time monitoring process of the present application, additional operation steps can also be added according to actual needs, such as seismic data filtering and denoising, seismic event positioning error quantification analysis, etc., to further improve the applicability and accuracy of the process.

[0098] Six, seismic event risk analysis

[0099] ​​S6. Based on the final positioning result and the magnitude result of the seismic event, a first evaluation method based on real-time dynamic b value and a second evaluation method based on the distribution of the ratio of P wave velocity to S wave velocity are used for seismic risk quantitative evaluation; in the first evaluation method, the b value is negatively correlated with the seismic event risk, and the value is cross-validated by the maximum likelihood estimation method and the b+ method; in the second evaluation method, the high-precision seismic event distribution is obtained by the double-difference positioning method, the waveform cross-correlation method is used to calculate the P wave travel time difference and the S wave travel time difference, and the ratio of P wave velocity to S wave velocity is estimated, and the estimated value of the ratio of P wave velocity to S wave velocity is positively correlated with the fluid effect.

[0100] It should be noted that obtaining real-time and accurate seismic positioning results is the core basis for preventing induced seismic risk. After obtaining the accurate earthquake catalog through the above process, the quantitative evaluation of seismic risk needs to be further carried out. Therefore, the application provides two evaluation methods for quantitative evaluation of seismic risk, which are a first evaluation method based on real-time dynamic b value and a second evaluation method based on the distribution of the ratio of P wave velocity to S wave velocity. The two evaluation methods will be described in detail below.

[0101] It should be noted that in step S6 of the application, the seismic b value is a core parameter for describing the proportional relationship of earthquakes in seismic activity, which is derived from the classic Gutenberg-Richter law (G-R law), and is a key indicator for seismic activity analysis, regional tectonic stress evaluation and seismic risk prediction. It is widely used to judge the change of seismic activity and potential rupture mechanism. The core of the G-R law is that the higher the magnitude, the lower the frequency of earthquake occurrence, and its mathematical expression is:

[0102]

[0103] Among them, is the total number of earthquakes (i.e. "cumulative frequency"); is the activity parameter (positively correlated with the total number of earthquakes in a certain region, the greater the value, the more frequent the seismic activity in the region); is the b value parameter; is the magnitude result of the seismic event.

[0104] Based on the real-time updated earthquake catalog, the b value change calculated in a specific spatial range and different time windows is the dynamic b value; at the same time, the reciprocal of the b value (i.e. ) is taken, when the dynamic b value rises, it can correspond to represent higher seismic risk, so as to realize the quantitative goal of induced seismic risk.

[0105] Considering that the single maximum likelihood estimation method for calculating b value is susceptible to algorithm limitations, incomplete seismic data and the number of earthquakes, in order to further improve the robustness of b value calculation and avoid misjudgment of singular values on the analysis results, the application introduces two b value estimation methods, maximum likelihood estimation method and b+ method, and verifies the results of the two methods with each other, and then quantifies the risk index of induced earthquakes, effectively reduces the misjudgment probability of earthquake risk, that is, b value and earthquake risk are negatively correlated, when b value is small, the earthquake risk is high, and vice versa, when b value is large, the earthquake risk is low. The maximum likelihood estimation method is a mature calculation method, and the b+ method is a recently disclosed method. In order to better understand the application by those skilled in the art, the principles of the two methods are introduced as follows.

[0106] 1) Maximum likelihood method for b value estimation

[0107] Since the monitoring range of dense stations is large, the range of the target area is first specified, including the minimum longitude, maximum longitude, minimum latitude and maximum latitude, then the earthquake event information corresponding to the range is read, the earthquake event information includes the occurrence time of the earthquake event and the magnitude of the earthquake event, and then the magnitude list Mags is sorted according to the occurrence time of the earthquake event. The magnitude list stores the magnitudes of multiple earthquake events.

[0108] After sorting the magnitude list, the embodiment converts it into an array and selects the earthquake events with a magnitude greater than the "modified cutoff magnitude", to ensure that the complete magnitude data near and above the cutoff magnitude (i.e. the completeness magnitude ) is included in the analysis. Then the number of effective earthquake events after screening is counted. If the number of effective earthquake events is insufficient, return NaN; if the number of effective earthquake events meets the statistical requirements , estimate the b value according to the following formula, and then calculate the activity parameter according to the total number of earthquakes and the estimated b value:

[0109]

[0110] Wherein, represents the estimated b value; represents a natural constant; represents a logarithmic operation; represents the average magnitude above the completeness magnitude; represents the modified cutoff magnitude; represents the magnitude interval width, which is set to 0.1 in this embodiment.

[0111] After estimating b value and calculating the activity parameter, two kinds of uncertainty estimation results are calculated: one is Aki uncertainty based on sample size, that is, the estimated b value is divided by the square root of the number of effective seismic events, and the calculated result is the uncertainty estimation result of the activity parameter ; the second is Shibolt uncertainty, which combines the deviation square sum of magnitude and average magnitude, the number of effective seismic events and other parameters to calculate, and the calculated result is the uncertainty estimation result of b value . Finally, the a value, b value and two kinds of uncertainty estimation results are returned. The specific calculation method is as follows:

[0112]

[0113]

[0114] Wherein, Mags represents the magnitude of the i th seismic event.

[0115] In the above process, the completeness magnitude is estimated by MBS method. The estimation of b value is based on the completeness magnitude in the earthquake catalog, that is, the minimum magnitude that can be completely recorded in the earthquake catalog. First, the magnitude distribution in the magnitude list Mags is counted, and the number of seismic events and the total number of earthquakes in each magnitude interval are obtained. For each magnitude interval, the corresponding b value and uncertainty are calculated with the completeness magnitude as the cutoff magnitude. The sliding window method is used to calculate the average b value of each cutoff magnitude and several magnitude intervals after it, and to judge the stability of b value (that is, whether the difference between b value and its average value is within the preset uncertainty range). If there is a stable b value interval, the first magnitude interval that meets the stability condition and is greater than the minimum value in the magnitude list is taken as the completeness magnitude, otherwise the magnitude interval with the most seismic events is taken as the completeness magnitude.

[0116] 2) b+ method for b value estimation

[0117] ​When measuring b value in active aftershock sequences, the short-term incompleteness of earthquake catalog and the data windowing process have significant effects on the results, and the biases introduced by these two factors are consistent with the target signal (true b value change), which can seriously interfere with the reliability of the results. To solve this problem, Vander Elst (2021) developed a new b value estimator, "b-positive" (hereinafter referred to as "b+ estimator"). This estimator has two major advantages: first, it is not sensitive to the instantaneous fluctuations of earthquake catalog completeness, and second, it does not rely on the traditional "data windowing" process, which can avoid the errors caused by window division. The core principle is to take the "positive difference subset" of the magnitude difference between consecutive earthquake events as the basis for analysis. This positive difference subset conforms to the double exponential (Laplace) distribution characteristics, and its corresponding b value is completely consistent with the b value of the earthquake magnitude distribution itself. Through this design, the b+ estimator greatly improves the robustness of continuous b value measurement in two key scenarios: on the one hand, it can adapt to the real-time b value monitoring needs during active seismic sequences, and on the other hand, it can also be applied to historical earthquake catalog analysis with unknown or fluctuating completeness over time, effectively solving the applicability limitations of traditional methods in complex scenarios.

[0118] In this embodiment, first input the earthquake time series (datetimes) and magnitude list (Mags) formed by the occurrence time of earthquake events, and set the magnitude difference threshold DIFF to 0.02. Also need to convert the magnitude list to an array. By calculating the second difference of each occurrence time relative to the start time of the earthquake time series, the time axis (hour_list) is obtained, and the magnitude difference (magdif) of adjacent earthquake events is calculated, and 0 is added at the beginning of the earthquake time series to align the index. Then, select the earthquake events with a magnitude difference greater than or equal to the magnitude difference threshold, count the number of these earthquake events (evt_num), and record the occurrence time (cata) and magnitude (mag) of each selected earthquake event. Dynamically set the sliding window size Np according to the size of evt_num (the more earthquake events, the larger Np, Np range is 45-100) to ensure statistical stability. For each selected earthquake event, when its index is greater than or equal to Np (to ensure enough data in the window), take the Np magnitude difference data before the earthquake event, and calculate the b value of the current window. At the same time, use the bootstrap sampling method (randomly replace the data in the window, extract Np-10 earthquake events each time, repeat Np times), and take the standard deviation of the sampled b value as the uncertainty of the current b value. Finally, extract the b value of all valid sliding windows to provide a basis for subsequent analysis. This process realizes the quantification of b value and its reliability in the magnitude difference sequence through dynamic windowing and statistical sampling.

[0119] In addition, it is necessary to point out that in the above two methods, the b value estimation is dynamically defined, that is, the number of selected seismic events is selected when calculating the b value. The reason for setting this parameter as a variable parameter is that the number of seismic event distribution in different regions is significantly different. If the same value is used, the number of regions with fewer events will result in sparse calculation results, and for the same reason, the number of events will result in unstable calculation results. Therefore, according to the number of seismic events in different regions or different time periods, the b value window number parameter is dynamically set, that is, according to the number of different seismic events, the b value window number is dynamically set to 60-120.

[0120] The most common mechanism of induced seismic risk is that the change of fault plane pore pressure caused by large fluid injection will reduce the normal stress on both sides of the fault plane when the pore pressure rises, so that the fault will lose stability under the shear stress of the stress field, and then induce an earthquake through rapid sliding. The matching relationship between the migration path of the injected fluid, the accumulation position and the spatial distribution of the fault is the key basis for judging the induced seismic risk. Because the fluid injected into the underground reservoir will change the physical properties of the medium, and The ratio of P wave velocity and S wave velocity has high sensitivity to the difference between dry pores and saturated pores, and the spatial and temporal variation characteristics of this parameter can be calculated to study the migration and accumulation of fluid, and then achieve the purpose of analyzing the induced seismic risk. The present application first obtains the high-precision distribution of seismic events by the double-difference positioning method (HypoDD); then calculates the P wave and S wave travel time difference between different seismic event pairs by waveform cross-correlation method; finally, the spatial and temporal distribution characteristics of are obtained by a specific calculation formula.

[0121] Specifically, the process of using the second evaluation method to quantitatively evaluate the seismic risk is as follows:

[0122] S61. Extract the occurrence time and source space position of each seismic event from the final positioning result as the seismic positioning data, obtain the longitude, latitude and seismic station name of the seismic station as the seismic station data, obtain the earthquake catalog containing all seismic events, obtain the travel time difference data by screening the seismic event pairs in the earthquake catalog, and store the three kinds of data in the form of a dictionary.

[0123] In the S61 of the present application, the travel time difference data is obtained in the following manner: a pair of earthquake events is formed by two different earthquake events in the earthquake catalog, the travel time difference of each earthquake event pair is calculated to form a travel time difference file, the target earthquake events are selected from the earthquake location data according to the latitude and longitude of the target area to form a target event earthquake ID list, the earthquake event pairs that meet the following conditions are selected: both earthquake events of the earthquake event pair are in the target event earthquake ID list, the earthquake event pair has at least S seismic stations that simultaneously record P-wave travel time difference and S-wave travel time difference, and the travel time difference of the selected earthquake event pair is used as the travel time difference data; wherein S is a preset seismic array quantity threshold.

[0124] In the S61 of the present embodiment, three types of core data required for seismological analysis need to be read and converted into a dictionary structure (for fast indexing) to provide a data basis for subsequent analysis. The first type of data is earthquake location data, which needs to be read from the final location result, and the unique ID of each earthquake event and the latitude, longitude and depth of the source are extracted to be stored in the form of a dictionary to realize fast query of the location of the earthquake event. The second type of data is seismic station data, which needs to be read from the seismic station parameter file (such as station.dat), and the longitude, latitude and seismic station name (sta) of the seismic station are extracted to be stored in the form of a dictionary for subsequent source-seismic station distance calculation. The third type of data is travel time difference data, which is generated based on the selected earthquake event pair, and is stored in the form of a dictionary to generate a temporary file (dtcc_tmp.txt) to avoid repeated calculation and improve efficiency.

[0125] In the S61 of the present embodiment, for the target area, the process of selecting target earthquake events is as follows: first, read the complete information (ID, latitude, longitude, depth, magnitude, and occurrence time) of the earthquake event according to the final location result file type (.pha or.reloc); then perform the first round of spatial screening, which is based on the latitude, longitude and depth range of the target area to retain the earthquake events located within the target area; then perform the second round of time screening, which is based on the preset earthquake event start time evt_start_time and earthquake event end time evt_end_time to retain the earthquake events within a specific time period. After two rounds of screening, all target earthquake events that meet the time and space requirements are obtained.

[0126] S62. For each earthquake event pair, select the stations that simultaneously record P-wave travel time difference and S-wave travel time difference to obtain the station list and the corresponding P-wave travel time difference and S-wave travel time difference, and perform mean removal operation on the two types of travel time difference respectively to obtain the P-wave actual travel time difference relative variation and the S-wave actual travel time difference relative variation. ​​​

[0127] In this embodiment S62, since the travel time difference is easily affected by station noise and seismic phase identification errors, multiple rounds of preprocessing are required to improve data quality and lay the foundation for subsequent fitting analysis. Specifically, for each earthquake event pair, stations that "simultaneously record P-wave and S-wave travel time differences" are selected (i.e., both station name + "P" and station name + "S" exist in the travel time difference dictionary, thus obtaining a list of stations and their corresponding P-wave travel times). S-wave travel time difference This ensures the integrity of data from individual stations. Then, the filtered P-wave travel time difference is analyzed. Time difference with S wave Calculate the mean travel time difference for each wave, and obtain the relative change in each type of travel time difference by subtracting the mean travel time difference from the observed travel time difference. This is the relative change in the actual travel time difference of the P wave. The relative change in the actual travel time difference of the S-wave By using the travel time difference to remove the mean, the interference of the overall system offset on the fitting results can be eliminated, highlighting the relative differences in travel time difference (which are directly related to the changes in underground medium velocity).

[0128] S63. A random sample consensus algorithm is used to coarsely screen the two relative changes in actual travel time. The two relative changes in actual travel time after coarse screening are recorded as the first change and the second change, respectively. A data point is formed by the first change and the second change. The vertical distance from the data point to the known straight line and the mean of all distances are calculated. The mean distance is multiplied by a preset hyperparameter to obtain the distance retention threshold. Data points with distances less than the distance retention threshold are retained. The known straight line is the straight line corresponding to the ratio of the theoretical P-wave velocity to the theoretical S-wave velocity of the underground medium being equal to 1.732.

[0129] In this embodiment S63, a two-step method of "coarse screening + fine screening" is used to eliminate abnormal relative changes in actual travel time. Specifically, the Random Sample Consensus (RANSAC) algorithm is used for coarse screening, and this embodiment is based on sklearn.linear_model.RANSAC Regressor. The linear model iteratively filters inlier points (non-outliers) to remove extreme outliers caused by station malfunctions or misjudgments of seismic phases. A distance threshold is used for fine-tuning; this embodiment calculates the data points. arrive For the perpendicular distance to the corresponding line, data points with a distance less than 0.8 times the mean distance are retained to further suppress residual noise. Among these, and These are the theoretical longitudinal wave velocity and the theoretical transverse wave velocity of the underground medium, respectively.

[0130] Further, the embodiment S63 also performs global clipping, i.e. all the effective data points of the earthquake event pairs are summarized, and the extreme abnormal values are clipped according to a fixed threshold (such as 0.4s), to ensure the consistency of the input fitting data.

[0131] S64. The P-wave velocity and S-wave velocity of the preset earthquake event are obtained based on the final positioning result and the coordinates of the seismic station to obtain the straight-line distance between the hypocenter and the station, and then the P-wave theoretical travel time is obtained by dividing the straight-line distance by the P-wave velocity, and the S-wave theoretical travel time is obtained by dividing the straight-line distance by the S-wave velocity.

[0132] In the embodiment S64, a preset average velocity model (such as P-wave velocity 5450m / s, S-wave velocity = P-wave velocity / 1.732) is adopted, based on the final positioning result of the earthquake event and the coordinates of the seismic station, the horizontal distance between the hypocenter and the seismic station is calculated by the Haversine formula, and the straight-line distance (oblique distance) between the hypocenter and the station is obtained combined with the depth position of the hypocenter. The P-wave theoretical travel time and the S-wave theoretical travel time of a single earthquake event are calculated by “oblique distance / velocity”.

[0133] S65. The P-wave theoretical travel time difference and the S-wave theoretical travel time difference of each earthquake event pair are calculated according to the P-wave theoretical travel time and the S-wave theoretical travel time of the earthquake event , and the P-wave theoretical travel time difference relative change and the S-wave theoretical travel time difference relative change are obtained by de-meaning the two kinds of theoretical travel time difference respectively.

[0134] In the embodiment S65, the P-wave and S-wave theoretical travel time difference and of the two earthquake events of the earthquake event pair are calculated respectively, and the P-wave and S-wave theoretical travel time difference relative change and are obtained after de-meaning, which are used as the reference benchmark of the observation data and provide the basis for the reliability test of the observation data.

[0135] S66. The least square method is used to fit the screened data points to obtain the preliminary estimated value of the ratio of P-wave velocity to S-wave velocity and the preliminary estimated value of the intercept, and the total least squares method is used to re-fit the two preliminary estimated values and the screened data points in a constantly iterative manner, and the difference between the theoretical travel time difference relative change and the actual travel time difference relative change is calculated as the residual during the iteration process, and the final estimated value of the ratio of P-wave velocity to S-wave velocity is obtained after the iteration is completed.

[0136] In the embodiment S65, the P-wave velocity and the S-wave velocity ​The core logic for solving the ratio is to fit a... The slope of the fitted line is... ratio, This is the intercept. Considering... and Both methods involve measurement errors; ordinary least squares only considers... Errors can introduce biases. Therefore, this embodiment adopts a two-layer strategy of "linear fitting of least squares method + iterative of overall least squares method" to balance efficiency and accuracy.

[0137] Specifically, the least squares method is implemented based on the numpy.polyfit function. A preliminary linear fit was obtained quickly. Preliminary estimates of the data and the intercept are used for rapid assessment of data quality (e.g., ...). The value must be greater than 0.5 to exclude obviously unreasonable fitting results.

[0138] Secondly, initialize with the least squares fitting result or the user-given value (initial_guess). Ratio and intercept; calculate the ratio to the fitted line for each data point. The signed vertical distance; calculate the objective function (sum of squared residuals) and gradient (for each case). , (partial derivatives); adaptive step-size gradient descent is used to adjust... and The initial step size is 1.0. If the new residual is greater than the original residual, the step size is halved until the residual decreases or the step size is less than 1e-12. and When the changes are all less than the convergence threshold (tol=1e-5) or the maximum number of iterations (max_iter=1e5) is reached, the iteration stops; the final output is the final fitted slope (Vp / Vs ratio), intercept, and residual sum of squares (reflecting the goodness of fit).

[0139] In this embodiment S65, the final estimated value of the ratio of P-wave velocity to S-wave velocity is related to the risk of earthquake events. Whether it is positive or negative depends on the value of b (which reflects the slip properties of the fault). The final estimated value is positively correlated with the influence of fluid. When the fluid increases, the risk of fault instability increases.

[0140] It should also be noted that the method for real-time earthquake monitoring and risk analysis using dense arrays 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 system for real-time earthquake monitoring and risk analysis using dense arrays, corresponding to the method provided in the above embodiments, such as...Figure 2 as shown, comprising:

[0141] a data acquisition module configured to acquire raw seismic data within a user-set geospatial range and time range;

[0142] a result acquisition module configured to process the raw seismic data according to the dense seismic array real-time monitoring and risk analysis method described in the above embodiments, calculate the b value and estimate the in-situ P-wave velocity to S-wave velocity ratio of the source region;

[0143] a visualization module configured to send a warning message or a structured risk assessment report to the user terminal in a preset manner when the calculated b value deviates from the expected value or exceeds the preset risk threshold.

[0144] It should be noted that when an earthquake event with a magnitude exceeding 2.5 is identified, or the 1 / b value rises significantly, reaching a value of 1.5 or more, it is determined that there is a possibility of an increase in seismic risk. Warning information can be sent to the user terminal through message pushing, email or SMS, and a risk assessment report can also be transmitted to the user terminal. The report can be presented in PDF or HTML format, which can include a time series graph of the calculated b value over time, a time series graph of the P-wave velocity to S-wave velocity ratio over time, a spatial distribution cloud map of the P-wave velocity to S-wave velocity ratio, etc. The report can also include a comprehensive risk assessment map integrating earthquake event distribution, b value risk indicators and abnormal regions. The form of feedback can be selected and adjusted by those skilled in the art according to actual needs, so it is not limited in the present application.

[0145] The present application will now be described through a specific example to demonstrate the application effect of the dense seismic array real-time monitoring and risk analysis method described in S1-S6 of the above embodiments on a specific data set, so as to facilitate understanding of the essence of the present application.

[0146] Embodiment

[0147] The steps of this embodiment are the same as the dense seismic array real-time monitoring and risk analysis method shown in the above S1-S6 steps, which will not be repeated here. This embodiment mainly demonstrates some specific parameter settings and implementation results. For ease of description, the method shown in S1-S6 steps is referred to as the present application method, denoted as RISEMAP.

[0148] The present application embodiment takes the practice of deploying a dense seismic array on a high-risk hydraulic fracturing platform in a certain region, conducting real-time seismic monitoring and dynamic seismic risk assessment as an example.

[0149] On a certain high-risk shale gas development platform in the above region, a real-time seismic monitoring array consisting of 55 seismic instruments is deployed, as shown in Figure 3The seismic observation equipment of the array has real-time data transmission capability, can realize real-time monitoring of seismic activity in the region, and the monitoring range focuses on covering the horizontal well trajectory and the surrounding fault distribution area. In order to ensure the monitoring accuracy, the seismic station layout adopts a differentiated spacing design: in the core area of the well trajectory, the spacing between seismic stations is set to 0.5 to 1 km; in the peripheral area around the well, the spacing between seismic stations is adjusted to 1 to 2 km, so as to ensure that small and micro seismic events can be accurately captured and monitored. Figure 3 In the figure, the black solid line is the shale gas development well trajectory, the red solid line is the known natural fault distribution near the platform, and the dark green triangle marks the position of the seismic station layout.

[0150] The data transmission is based on a high-efficiency and stable data transmission link built on 4G data transmission technology, with "cloud hub" as the core, realizing reliable flow of seismic data from the acquisition end to the user end. In the field monitoring scene, the instrument deployed has a built-in SIM card as the basis for data transmission, which quickly pushes the real-time collected seismic data to the nearby base station in a wireless communication manner, forming data storage in the cloud. Then, with the help of the high-speed channel of the Internet, the data is seamlessly connected to the data server, completing the migration from the physical terminal to the virtual cloud. The cloud server plays the role of "intelligent housekeeper" in this transmission process.

[0151] Then, real-time data stream is picked up by PhaseNet, and the arrival time key information of P wave, S wave and other seismic phases is extracted from the seismic waveform. Through 500Hz high sampling rate, the picking accuracy is ensured, and the reliability of subsequent seismic analysis is improved. The REAL algorithm is used to complete the association of seismic phases: a three-dimensional grid is set at the reference position, the number of seismic phases of each grid point is counted combined with the travel time list, the grid point with the most seismic phases and the smallest travel time residual is selected as the initial position, and after simulated annealing optimization, the unreliable arrival time is excluded, and the effective seismic event is determined according to the seismic phase quantity threshold. Based on the REAL preliminary positioning result, the Hypoinverse, VELEST and other absolute positioning methods are used to quickly invert the source position, and at the same time, combined with the regional seismic wave attenuation law, the appropriate magnitude calculation method is selected to update the magnitude of the seismic event in real time, such as Figure 4 as shown in the figure.

[0152] Based on the real-time obtained source position and magnitude, RISEMAP real-time inverts the seismic source mechanism, and synchronously uploads the results to the real-time dynamic display system. Users can independently set the geographical space range and event time period in the system, view the dynamic distribution of earthquakes in the corresponding interval, and calculate the dynamic b value change, a key risk indicator, such as Figure 5 as shown in the figure, to provide data support for seismic risk assessment.

[0153] It can be understood that the dense seismic array real-time monitoring and risk analysis method described in S1-S6 can be realized by a computer program.

[0154] Similarly, based on the same inventive concept, another preferred embodiment of the present application also provides a computer electronic device corresponding to the dense seismic array real-time monitoring and risk analysis method provided by the above-mentioned embodiments, as shown in the figure, which comprises a memory and a processor. Figure 6

[0155] The memory is used to store a computer program.

[0156] The processor is used to realize the dense seismic array real-time monitoring and risk analysis method in the above-mentioned embodiments when the computer program is executed.

[0157] In addition, the logical instructions in the above-mentioned memory can be realized in the form of a software functional unit and sold or used as an independent product, which can be stored in a computer readable storage medium. Based on such understanding, the technical solutions of the present application or the part of the prior art that essentially contributes or the part of the technical solutions can be embodied in the form of a software product, which is stored in a storage medium and includes a plurality of instructions for making a computer device (which can be a personal computer, a server, or a network device, etc.) execute all or part of the steps of the method described in the embodiments of the present application.

[0158] Therefore, based on the same inventive concept, another preferred embodiment of the present application also provides a computer readable storage medium corresponding to the dense seismic array real-time monitoring and risk analysis method provided by the above-mentioned embodiments, and the storage medium stores a computer program, and when the computer program is executed by a processor, the dense seismic array real-time monitoring and risk analysis method in the above-mentioned embodiments can be realized.

[0159] It can be understood that the above-mentioned storage medium can include a random access memory (RAM) and can also include a non-volatile memory (NVM), such as at least one disk memory. Meanwhile, the storage medium can also be a U disk, a mobile hard disk, a magnetic disk or an optical disk, and various media that can store program codes.

[0160] ​It can be understood that the processor described above can be a general processor, including a central processing unit (CPU), a network processor (NP), etc.; can also be a digital signal processor (DSP), an application specific integrated circuit (ASIC), a field-programmable gate array (FPGA) or other programmable logic device, a discrete gate or transistor logic device, a discrete hardware component.

[0161] In addition, it should be noted that those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working process of the system described above can refer to the corresponding process in the foregoing method embodiments, which will not be described here. In the embodiments provided in the present application, the division of steps or modules in the system and method is only a logical functional division, and there can be another division mode in actual implementation, for example, multiple modules or steps can be combined or integrated together, or a module or step can be split.

[0162] The above-described embodiments are only a preferred scheme of the present application, and are not intended to limit the present application. Those skilled in the art can make various changes and modifications without departing from the spirit and scope of the present application. Therefore, any technical scheme obtained by equivalent replacement or equivalent transformation falls within the protection scope of the present application.

Claims

1. A method for dense array seismic real-time monitoring and risk analysis, characterized in that, The method comprises the following steps: S1. Real-time acquisition of original seismic data continuously monitored by a dense station array, and format conversion and segmentation of the original seismic data to form preprocessed seismic data; S2. Inputting the preprocessed seismic data into a machine learning-based phase picking model to output the phase picking results corresponding to the seismic data; S3. Using a fast seismic correlation and positioning method to correlate discrete seismic phases in the phase picking results into a specific seismic event, outputting the seismic correlation results, and directly returning to step S1 when not correlated to a seismic event; S4. Positioning the seismic correlation results by a two-step positioning method to obtain an accurate spatial position of the seismic source through independent inversion of the seismic plane position and depth position, and outputting the final positioning results; S5. Reading the seismic event information from the final positioning results, calculating the magnitude of each seismic station, and taking the median of the magnitudes of all seismic stations as the magnitude result of the seismic event; S6. Based on the final positioning results and the magnitude result of the seismic event, using a first evaluation method based on real-time dynamic b value and a second evaluation method based on the distribution of P wave velocity and S wave velocity ratio for seismic risk quantification evaluation; in the first evaluation method, the b value is negatively correlated with the risk of the seismic event, and the b value is cross-validated by the maximum likelihood estimation method and the b+ method; in the second evaluation method, high-precision seismic event distribution is obtained by a double-difference positioning method, P wave travel time difference and S wave travel time difference are calculated by a waveform cross-correlation method, and the ratio of P wave velocity to S wave velocity is estimated, and the estimated value of the ratio of P wave velocity to S wave velocity is positively correlated with the fluid effect.

2. The dense array seismic real-time monitoring and risk analysis method of claim 1, wherein, In step S1, the original seismic data is converted into SAC format, and the seismic data after format conversion is cut into 1-minute-long data records as preprocessed seismic data; When detecting a transmission interruption of seismic data, first detect the difference between the seismic data acquisition time and the real world time: if the difference is less than a preset time length threshold, it is determined that the current seismic station is in real-time processing state, and no backtracking processing is performed; if the difference is greater than or equal to the time length threshold, read the previous interruption time from the record log, take the seismic data corresponding to the previous interruption time as a historical data node, and start backtracking processing from the historical data node until the difference is less than the time length threshold.

3. The dense array seismic real-time monitoring and risk analysis method of claim 1, wherein, In step S2, the phase picking model uses PhaseNet, which outputs the phase picking results with a probability higher than a preset phase retention threshold in CSV file format. The phase picking results include seismic station ID, seismic data acquisition time, number of phase picking sampling points, phase arrival time, phase picking probability value, amplitude value at the picked phase, phase type, and seismic data storage path. The phase type is P wave phase or S wave phase.

4. The dense array seismic real-time monitoring and risk analysis method of claim 1, wherein, In the seismic correlation process of step S3, in addition to reading the phase picking results of the current preprocessed seismic data, the phase picking results of the last Ttravel seconds in the previous preprocessed seismic data are also read synchronously, where Ttravel represents the maximum seismic wave travel time calculated based on the maximum distance of possible seismic event occurrence and seismic wave velocity.

5. The dense array seismic real-time monitoring and risk analysis method of claim 1, wherein, In step S4, the specific process of positioning by using the two-step positioning method is as follows: S41. Convert the earthquake association result file into a data input format suitable for the Hypoinverse algorithm or the VELEST algorithm; S42. After the format conversion is completed, the plane position of the earthquake source of the seismic event is obtained by using the Hypoinverse algorithm or the VELEST algorithm for preliminary inversion; S43. According to a preset epicenter distance threshold, only the earthquake association results within the epicenter distance threshold range are selected and read, and the Hypoinverse algorithm or the VELEST algorithm is used to perform positioning calculation based on the selected earthquake association results to obtain the depth position of the earthquake source of the seismic event, and the plane position and the depth position together constitute the spatial position of the earthquake source.

6. The dense array seismic real-time monitoring and risk analysis method of claim 1, wherein, In step S6, the process of quantitatively evaluating the earthquake risk by using the second evaluation method is as follows: S61. Extract the earthquake occurrence time and the spatial position of the earthquake source of each seismic event from the final positioning result as the earthquake positioning data, obtain the longitude, latitude and seismic station name of the seismic station as the seismic station data, obtain the earthquake catalog containing all seismic events, obtain the travel time difference data by screening the seismic event pairs in the earthquake catalog, and store the three kinds of data in the form of a dictionary; S62. For each seismic event pair, select the stations that record the P-wave travel time difference and the S-wave travel time difference at the same time to obtain the station list and the corresponding P-wave travel time difference and S-wave travel time difference, and perform mean removal operation on the two kinds of travel time differences respectively to obtain the P-wave actual travel time difference relative variation and the S-wave actual travel time difference relative variation; S63. The random sample consensus algorithm is used to coarsely screen the two kinds of actual travel time difference relative variations, and the coarsely screened two kinds of actual travel time difference relative variations are respectively denoted as the first variation and the second variation; a data point is formed by the first variation and the second variation, the perpendicular distance from the data point to the known straight line and the mean value of all distances are calculated, the distance mean value is multiplied by a preset hyperparameter to obtain a distance retention threshold, and the data points with a distance less than the distance retention threshold are retained; wherein, the known straight line is a corresponding straight line with a theoretical P-wave velocity to theoretical S-wave velocity ratio equal to 1.732; S64. The P-wave velocity and the S-wave velocity of the seismic event are preset, the straight line distance between the earthquake source and the station is obtained based on the final positioning result and the station coordinates, the P-wave theoretical travel time is obtained by dividing the straight line distance by the P-wave velocity, and the S-wave theoretical travel time is obtained by dividing the straight line distance by the S-wave velocity; S65. The P-wave theoretical travel time difference and the S-wave theoretical travel time difference of each seismic event pair are calculated according to the P-wave theoretical travel time and the S-wave theoretical travel time of the seismic event, and mean removal operation is performed on the two kinds of theoretical travel time differences respectively to obtain the P-wave theoretical travel time difference relative variation and the S-wave theoretical travel time difference relative variation; S66. The least squares method is used to fit the screened data points to obtain a preliminary estimate of the P-wave velocity to S-wave velocity ratio and a preliminary estimate of the intercept. The total least squares method is used to re-fit the two preliminary estimates and the screened data points in an iterative manner. The difference between the relative change in the theoretical travel time difference and the relative change in the actual travel time difference is calculated as a residual during the iteration process. The final estimate of the P-wave velocity to S-wave velocity ratio is obtained after the iteration process is completed.

7. The dense array seismic real-time monitoring and risk analysis method of claim 6, wherein, In S61, the travel time difference data is obtained as follows: a pair of seismic events is formed by two different seismic events in the earthquake catalog, the travel time difference of each seismic event pair is calculated to form a travel time difference file, the target seismic event ID list is formed by screening the target region latitude and longitude from the earthquake location data, and the seismic event pairs that meet the following conditions are screened: both seismic events of the seismic event pair are in the target event ID list, and the seismic event pair has at least S seismic stations that record P-wave travel time difference and S-wave travel time difference at the same time. The travel time difference of the screened seismic event pair is used as the travel time difference data; wherein S is a preset seismic array quantity threshold.

8. A dense array seismic real-time monitoring and risk analysis system, characterized in that, The method comprises: a data acquisition module for acquiring original seismic data within a user-set geographic space range and time range; a result acquisition module for processing the original seismic data according to the dense array seismic real-time monitoring and risk analysis method of any one of claims 1-7, calculating the b value, and estimating the in-situ P-wave velocity to S-wave velocity ratio of the source region; a visualization module for sending a warning message or a structured risk assessment report to the user end in a preset manner when the calculated b value deviates from the expected value or exceeds the preset risk threshold.

9. A computer-readable storage medium, characterized in that, The storage medium stores a computer program, and when the computer program is executed by the processor, the dense array seismic real-time monitoring and risk analysis method of any one of claims 1-7 is implemented.

10. A computer electronic device, comprising: The device comprises a memory and a processor; The memory is used to store a computer program; The processor is used to implement the dense array seismic real-time monitoring and risk analysis method of any one of claims 1-7 when the computer program is executed.

Images