A hydraulic fracturing microseismic event identification method for real-time ground monitoring
Patent Information
- Application Number
- CN202410129602.0
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2024-01-30
- Publication Date
- 2026-09-25
- Estimated Expiration
- 2044-01-30
AI Technical Summary
[0003]传统的基于震相拾取的方法针对单台设备,容易受噪声干扰,低信噪比情况下很难准确识别初至信号,导致定位误差较大
[0042]本发明通过若干个目标点对压裂段区域进行覆盖,找到叠加能量最强的目标点,避免叠加能量小于事件识别的阈值,对地面监测中获取的信噪比较低的数据,仍能提供准确的识别结果,避免遗漏微弱的事件;同时相比传统的偏移叠加定位方法遍历所有网格及发震时刻的方式,本发明的计算效率大幅提高,适用于地面微震实时监测的事件识别。
Smart Images

Figure CN117950031B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of hydraulic fracturing ground microseismic monitoring technology, specifically to a method for identifying hydraulic fracturing microseismic events for real-time ground monitoring. Background Technology
[0002] Hydraulic fracturing is an effective technique for enhancing oil and gas production and developing low-permeability oilfields. During hydraulic fracturing, seismic waves are generated along with the fracturing of rocks. By deploying seismographs on the ground for microseismic monitoring, analyzing the received waveforms, identifying microseismic events, and determining parameters such as the source location and focal mechanism, the spatial distribution characteristics of fractures can be obtained in near real-time, providing crucial information for evaluating the fracturing effect.
[0003] Traditional phase-based methods rely on single devices and are susceptible to noise interference. In low signal-to-noise ratio conditions, they struggle to accurately identify the first arrival signal, leading to significant positioning errors. Existing offset-stacking positioning methods traverse all grids and earthquake occurrence times, calculating the superimposed energy of each grid cell in the 3D mesh for each occurrence time. This results in high computational complexity and makes them unsuitable for real-time monitoring.
[0004] Therefore, it is evident that developing a microseismic event identification method that is applicable to low signal-to-noise ratio data and has high computational efficiency is crucial for real-time ground monitoring of hydraulic fracturing. Summary of the Invention
[0005] To address the aforementioned technical problems, this invention provides a method for identifying hydraulic fracturing microseismic events for real-time ground monitoring, specifically comprising the following steps:
[0006] S1. Calculate the spatial position of the midpoint of the fracturing section based on the well trajectory. Based on the fracture extension range, select several target points at the farthest extension range on both sides and at equal intervals in the middle, and calculate the spatial position of each target point.
[0007] S2. Based on the velocity model, target point location, and seismograph location, calculate the travel time from each target point to the seismograph.
[0008] S3. Receive real-time data streams through the microseismic monitoring system, put the data into a buffer, and then extract segmented data from the real-time data buffer;
[0009] S4. Preprocess the segmented data;
[0010] S5. Use the travel time to offset and superimpose the preprocessed segmented data to obtain the event recognition waveform;
[0011] S6. Determine whether there is a valid event based on the event recognition waveform. If there is a valid event in the event recognition waveform, locate the event location and update the target point list, and update step S2. If there is no valid event in the event recognition waveform, return to step S3.
[0012] Furthermore, the well trajectory in S1 consists of four columns of data denoted as (x, y, z, md), and the midpoint of the fracturing section is denoted as target point 1. The target point 1 is calculated by finding two points A and B near the fracturing section in the fourth column of well depth md, using linear interpolation to obtain the position coordinates of target point 1, and denoting them as P1.
[0013] Furthermore, calculate two target points perpendicular to the fracturing section and reaching the maximum fracture extension distance, denoted as P2 and P3. P2 and P3 have the same depth as P1, and the fracture extension distance is denoted as Len. The angle between the horizontal projection of the well trajectory curve of the fracturing section at point P1 and the due north direction is denoted as radians θ. Then, the horizontal positions of P2 and P3 satisfy the formula:
[0014] P2.x = P1.x - cos(θ) * Len
[0015] P2.y=P1.y+sin(θ)*Len
[0016] P3.x = P1.x + cos(θ) * Len
[0017] P3.y=P1.y-sin(θ)*Len
[0018] Furthermore, between positions P2 and P3, a target point is taken every 50m, and the nth target point is denoted as Pn, resulting in target point P4-Pn. The position of Pn satisfies the formula:
[0019]
[0020]
[0021] Furthermore, in S2, ray tracing is used to calculate the travel time from each target point to all seismographs to obtain the travel time table corresponding to each target point. Each seismograph includes a seismic trace signal with one vertical component and two horizontal components.
[0022] Furthermore, in S4, the segmented data undergoes preprocessing, including mean removal, pre-whitening, and bandpass filtering of the seismic trace signals of the vertical and horizontal components, respectively. The pre-whitening process for the segmented data includes:
[0023] The seismic trace signal is subjected to a fast Fourier transform to calculate the amplitude at each frequency point, and the minimum amplitude values are removed and smoothed.
[0024] Divide the amplitude spectrum of the seismic trace signal by the amplitude value and then perform an inverse Fourier transform to obtain the pre-whitened data.
[0025] Furthermore, the P-wave characteristic waveform CF is calculated from the preprocessed vertical component data. 1 The S-wave characteristic waveform CF is calculated from the preprocessed horizontal component data. 2 and the characteristic waveform CF 1 and CF 2 The following processes are performed: Equalization and normalization are applied.
[0026] The P-wave characteristic waveform CF 1 The calculation formula is:
[0027]
[0028] Among them, X i = (X1, X2, ..., X p ) represents the preprocessed signal, m represents the number of sampling points in the long time window, and n represents the number of sampling points in the short time window;
[0029] The S-wave characteristic waveform CF 2 The calculation formula is:
[0030] X(j)=x(j)+iH{x(j)}
[0031] Y(j)=y(j)+iH{y(j)}
[0032]
[0033] Let the eigenvalues of the covariance matrix M(j) be λ1 and λ2, and λ1 ≥ λ2, then CF 2 (j)=λ1(j) 2 ,
[0034] Where X(j) and Y(j) are the analytical forms of the seismic trace signals of the two horizontal components, and H is the Hilbert transform;
[0035] For the characteristic waveform CF 1 and CF 2 The calculation formula for the balancing process is as follows:
[0036]
[0037] Here, we assume that the characteristic waveform data with N sampling points is divided into multiple time windows, and each time window has 2M+1 sampling points, a j For the characteristic waveform data before processing, a j ′ represents the processed characteristic waveform data.
[0038] Furthermore, the offset superposition formula in S5 is as follows:
[0039]
[0040] Where S(t) is the event recognition waveform, t is the sampling time, M is the number of stations superimposed, and T is the time of sampling. n The travel time from the target point to each seismograph.
[0041] Compared with the prior art, the present invention has the following beneficial effects:
[0042] This invention covers the fracturing section area with several target points to find the target point with the strongest superimposed energy, avoiding situations where the superimposed energy is less than the threshold for event identification. It can still provide accurate identification results for data with low signal-to-noise ratio obtained from ground monitoring, avoiding the omission of weak events. At the same time, compared with the traditional offset superposition positioning method that traverses all grids and seismic times, the computational efficiency of this invention is greatly improved, making it suitable for event identification in real-time ground microseismic monitoring. Attached Figure Description
[0043] Figure 1 This is a flowchart of the present invention;
[0044] Figure 2 The waveform of the seismic data before preprocessing in this invention;
[0045] Figure 3 This is the preprocessed seismic data waveform of the present invention;
[0046] Figure 4 This refers to the three-component data of a single seismograph in this invention, as well as the characteristic waveforms of P-waves and S-waves.
[0047] Figure 5 These are the waveforms for identifying P-wave and S-wave events in this invention. Detailed Implementation
[0048] To make the technical solutions and effects of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments.
[0049] like Figure 1 As shown, this invention aims to provide a method for identifying hydraulic fracturing microseismic events that can be accurately identified under low signal-to-noise ratio conditions and used for real-time ground monitoring. The method specifically includes the following steps:
[0050] Step 1: Calculate the spatial location of the midpoint of the fracturing section (i.e., target point 1) based on the well trajectory. According to the designed fracture extension range, select several target points 2-n at the farthest extension points on both sides and at equal intervals in the middle, and calculate their spatial locations. Specifically:
[0051] (1) The well trajectory consists of four columns of data (x, y, z, md). The target point 1 is calculated by finding two points A and B near the fracturing section in the fourth column of well depth md, and using linear interpolation to obtain the position coordinates of the target point 1, which is denoted as P1.
[0052] (2) Calculate the two target points perpendicular to the fracturing section to reach the maximum fracture extension distance, denoted as P2 and P3. The depths of P2 and P3 are the same as those of P1. The designed fracture extension distance is denoted as Len. The angle between the horizontal projection of the well trajectory curve of the fracturing section at point P1 and the due north direction is radians θ. Then the horizontal positions of P2 and P3 satisfy the formula:
[0053] P2.x = P1.x - cos(θ) * Len
[0054] P2.y=P1.y+sin(θ)*Len
[0055] P3.x = P1.x + cos(θ) * Len
[0056] P3.y=P1.y-sin(θ)*Len
[0057] Between P2 and P3, a target point is selected every 50m, denoted as Pn for the nth target point, resulting in target point P4-Pn. The position of Pn then satisfies the formula:
[0058]
[0059]
[0060] Step 2: For each target point, based on the velocity model, target point location, and seismograph location (each seismograph has one vertical component and two horizontal component seismic traces), use ray tracing to calculate the travel time from the target point to all seismographs, obtaining the travel time table corresponding to each target point.
[0061] Step 3: Receive real-time data streams through the microseismic monitoring system and put the data into a buffer. If the buffer contains data for the current time period, extract segmented data from the buffer.
[0062] Step 4: Preprocess the segmented data, including mean removal, pre-whitening, and bandpass filtering of the vertical and horizontal components respectively. Specifically, the pre-whitening calculation method includes performing an FFT (Fast Fourier Transform) on the seismic trace signal, calculating the amplitude at each frequency point, removing and smoothing the minimum amplitude values, dividing the amplitude spectrum by the aforementioned amplitude, and then performing an IFFT (Inverse Fast Fourier Transform) to obtain the pre-whitened signal. (Refer to...) Figure 2 For the data waveform before preprocessing, Figure 3 This is the waveform of the preprocessed data.
[0063] Step 5, as follows Figure 4 The three-component data and characteristic waveforms of P-wave and S-wave from a single seismograph are shown. The characteristic waveform CF of the P-wave is calculated from the preprocessed vertical component data. 1 The S-wave characteristic waveform CF is calculated from the preprocessed horizontal component. 2 ; for the characteristic waveform CF 1 and CF 2 Perform balancing and normalization. Specifically:
[0064] (1) P-wave characteristic waveform CF 1 The calculation formula is:
[0065]
[0066] Among them, X i = (X1, X2, ..., X p ) represents the preprocessed signal, m represents the number of sampling points in the long time window, and n represents the number of sampling points in the short time window;
[0067] (2) S-wave characteristic waveform CF 2 The calculation formula is:
[0068] X(j)=x(j)+iH{x(j)}
[0069] Y(j)=y(j)+iH{y(j)}
[0070] ^ represents conjugate, and the covariance matrix M(j) is:
[0071]
[0072] Let the eigenvalues of the covariance matrix M(j) be λ1 and λ2, and λ1 ≥ λ2, then CF 2 (j)=λ1(j) 2 Where X(j) and Y(j) are the analytical forms of the seismic trace signals of the two horizontal components, and H is the Hilbert transform;
[0073] (3) For the characteristic waveform CF1 and CF 2 The calculation formula for the balancing process is as follows:
[0074]
[0075] Here, we assume that the characteristic waveform data with N sampling points is divided into multiple time windows, and each time window has 2M+1 sampling points, a j For the characteristic waveform data before processing, a j ′ represents the processed characteristic waveform data.
[0076] Step 6: For each target point, use the corresponding theoretical travel time to calculate the CF. 1 After offset superposition, the event recognition waveform Stack1 is obtained, and the S-wave characteristic waveform CF is obtained. 2 After offset stacking, the event recognition waveform Stack2 is obtained. For example... Figure 5 The P-wave and S-wave event identification waveforms Stack1 and Stack2 are shown in the figure. The vertical dashed line in the figure indicates the location of the maximum value picked up. The offset stacking formula in this step is:
[0077]
[0078] Where S(t) is the event recognition waveform, t is the sampling time, M is the number of stations superimposed, and T is the time of sampling. n The travel time from the target point to each seismograph.
[0079] Step 7: In the sliding window, find the maximum values Max1 and Max2 of Stack1 and Stack2. If Max1 or Max2 in the window exceeds the set threshold, then the location is determined by the subsequent location algorithm; otherwise, proceed to the next window. If no event is detected, continue to step 3.
[0080] Step 8: Based on the positioning results, if the distance between the positioning result and the nearest target point is greater than 100m, store the positioning result point in the target point list and calculate the travel time from that point to each seismograph. Continue to Step 3.
[0081] Since fracturing typically focuses on microseismic events near the fracturing section, and the travel times of nearby microseismic events are similar, overlaying the travel time data of the closest target point can enhance the superimposed energy. By covering the fracturing section area with several target points, the target point with the strongest superimposed energy can be found, preventing the superimposed energy from falling below the event identification threshold. Even with low signal-to-noise ratio data acquired from ground monitoring, this method still provides accurate identification results compared to traditional seismic phase identification methods, avoiding the omission of weak events. Furthermore, compared to traditional migration overlay positioning methods that traverse all grids and seismic times, the computational efficiency is significantly improved, making it suitable for event identification in real-time ground microseismic monitoring.
[0082] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.
Claims
1. A method for identifying hydraulic fracturing microseismic events in real-time ground monitoring, characterized in that, Specifically, the following steps are included: S1. Calculate the spatial position of the midpoint of the fracturing section based on the well trajectory. Based on the fracture extension range, select several target points at the farthest extension range on both sides and at equal intervals in the middle, and calculate the spatial position of each target point. S2. Based on the velocity model, target point location, and seismograph location, calculate the travel time from each target point to the seismograph. S3. Receive real-time data streams through the microseismic monitoring system, put the data into a buffer, and then extract segmented data from the real-time data buffer; S4. Preprocess the segmented data; S5. Use the travel time to offset and superimpose the preprocessed segmented data to obtain the event recognition waveform; S6. Determine whether there is a valid event based on the event recognition waveform. If there is a valid event in the event recognition waveform, locate the event location and update the target point list, and update step S2. If there is no valid event in the event recognition waveform, return to step S3.
2. The method for identifying hydraulic fracturing microseismic events for real-time ground monitoring according to claim 1, characterized in that, The well trajectory in S1 consists of four columns of data, denoted as follows: The midpoint of the fracturing section is denoted as target point 1. Target point 1 is calculated by finding two points A and B near the fracturing section at the fourth well depth md, using linear interpolation to obtain the coordinates of target point 1, and recording them as follows: .
3. The method for identifying hydraulic fracturing microseismic events for real-time ground monitoring according to claim 2, characterized in that, Calculate two target points perpendicular to the fracturing section that reach the maximum fracture extension distance, and denote them as follows: and ,and , and The cracks are at the same depth, and the distance they extend is denoted as... The well trajectory curve of the fracturing section is in The angle between the horizontal projection of a point and the direction of true north is denoted in radians. ,but and The horizontal position satisfies the formula: 。 4. The method for identifying hydraulic fracturing microseismic events for real-time ground monitoring according to claim 3, characterized in that, In the and In the middle of the position, take a target point every 50m, and record it as the first... The target points are , obtain the target point - The The position satisfies the formula: 。 5. The method for identifying hydraulic fracturing microseismic events for real-time ground monitoring according to claim 1, characterized in that, In step S2, ray tracing is used to calculate the travel time from each target point to all seismographs to obtain the travel time table corresponding to each target point. Each seismograph includes a seismic trace signal with one vertical component and two horizontal components.
6. The method for identifying hydraulic fracturing microseismic events for real-time ground monitoring according to claim 5, characterized in that, S4 involves preprocessing the segmented data, including mean removal, pre-whitening, and bandpass filtering of the seismic trace signals for the vertical and horizontal components, respectively. The pre-whitening process for the segmented data includes: The seismic trace signal is subjected to a fast Fourier transform to calculate the amplitude at each frequency point, and the minimum amplitude values are removed and smoothed. Divide the amplitude spectrum of the seismic trace signal by the amplitude value and then perform an inverse Fourier transform to obtain the pre-whitened data.
7. The method for identifying hydraulic fracturing microseismic events for real-time ground monitoring according to claim 6, characterized in that, Calculate the P-wave characteristic waveform from the preprocessed vertical component data. S-wave characteristic waveforms were calculated from the preprocessed horizontal component data. and the characteristic waveform and The following processes are performed: Equalization and normalization are applied. The P-wave characteristic waveform The calculation formula is: in, =( , ,…, ) The preprocessed signal This represents the number of sampling points within the long time window. This represents the number of sampling points within the short time window; The S-wave characteristic waveform The calculation formula is: Let the covariance matrix be... eigenvalues and ,and ≥ ,but , in, and This is the analytical form of the seismic trace signal for the two horizontal components. For Hilbert transform; For the characteristic waveform and The calculation formula for the balancing process is as follows: Among them, let it be that it has The characteristic waveform data of each sampling point is divided into multiple time windows, and each time window contains... One sampling point, For the characteristic waveform data before processing, This is the processed characteristic waveform data.
8. The method for identifying hydraulic fracturing microseismic events for real-time ground monitoring according to claim 1, characterized in that, The offset superposition formula in S5 is: in, To identify waveforms for events, Sampling time, This refers to the number of stations being superimposed. The travel time from the target point to each seismograph.
Citation Information
Patent Citations
Microseism interference imaging method
CN104765064A
Micro-seismic offset superposition positioning method and device, electronic equipment and readable storage medium
CN115951403A