Microseismic interference positioning method, system and device based on time shift correlation coefficient, and medium

Abstract

The invention discloses a microseism interference positioning method, system and device based on a time shift correlation coefficient, and a medium, and the method comprises the steps: constructing a DTW time shift correlation coefficient through a dynamic time warping algorithm, multiplying the constructed time shift correlation coefficient with an imaging operator, superposing all time shift correlation gathers at a same grid point, obtaining a single grid point imaging section, and carrying out the imaging of the single grid point. And all single grid point imaging profiles are further superposed, a final imaging effect representing the position of the seismic source is obtained, the position with the maximum imaging value is the seismic source, and seismic source positioning is realized. According to the method, the phenomenon that the imaging resolution of a cross-correlation imaging operator is reduced due to the influence of a wave propagation effect, noise and the like can be solved, and the positioning effect of the method under the influence of noise pollution is improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of earthquake detection, and in particular relates to a microseismic interferometric positioning method, system, equipment and medium based on a time-shift correlation coefficient. Background Art

[0002] In the process of oil and gas extraction, hydraulic fracturing technology is usually used to increase the production of oil and gas wells. The fracturing operation causes underground rocks to break and radiate energy outward, causing the rocks to move. These movements are similar to natural earthquakes along the fault plane, but with a smaller magnitude, so they are called "microearthquakes." A large number of microseismic events can provide enough seismic data for seismic processing and interpretation.

[0003] Microseismic monitoring technology uses microseismic events induced by underground rock fractures to monitor reservoirs and describe underground fractures. Unlike traditional seismic exploration, the location, intensity and time of the earthquake source in microseismic monitoring are unknown, so source location is the basis of microseismic monitoring. Traditional microseismic positioning methods use the picked arrival time and a given velocity model to perform ray-based linear travel time inversion to determine the source location. However, due to the small magnitude of microseismic events and low signal-to-noise ratio, it is difficult to accurately pick arrival times and locate the source. Compared with traditional travel time inversion methods, waveform-based positioning methods do not require phase picking operations and can detect and locate more events using low signal-to-noise ratio data.

[0004] Waveform-based positioning methods can be roughly divided into reverse time migration imaging and waveform stacking methods that mainly use the first arrival phase. Reverse time migration imaging requires wave field extrapolation, which is time-consuming and the imaging conditions used have a great impact on the source positioning results. The waveform stacking method initially mainly used the diffraction stacking operator to focus or back-project the energy along the travel time curve onto the spatial grid points. In recent years, this type of method has been successfully used to locate natural earthquakes and induced seismic activities related to mining operations, geothermal extraction, and oil and gas reservoirs. However, the diffraction stacking method needs to search for the time of earthquake occurrence when imaging, and the error of the earthquake occurrence time will also affect the final positioning accuracy.

[0005] The cross-correlation stacking source location method, which belongs to the interferometric imaging method, is a relatively new waveform stacking location method. This method uses the travel time difference information extracted by the cross-correlation operator for stacking imaging, which can avoid searching for the time of earthquake occurrence, and weighting different travel time difference information according to the characteristics of different seismic phases can also improve the positioning accuracy. Under the influence of factors such as focal mechanism, medium change and environmental noise, microseismic waveforms will produce different time differences at different times at different geophones. The overall time shift calculated by cross-correlation cannot achieve good alignment of the waveform of microseismic events, which in turn affects the positioning accuracy. At present, most methods use characteristic functions to transform the original waveform. The characteristic function can effectively avoid the influence of the first arrival waveform change and reduce the sensitivity to velocity errors, but it will sacrifice the spatial resolution of imaging and is easily affected by noise. Therefore, it is necessary to find a method that can achieve precise alignment between changing waveforms and use it to accurately locate microseisms. Summary of the invention

[0006] The purpose of the present invention is to provide a microseismic interferometric positioning method, system, device and medium based on time-shift correlation coefficient to solve the problems existing in the above-mentioned prior art.

[0007] To achieve the above object, the present invention provides a microseismic interferometry positioning method based on time-shift correlation coefficient, comprising:

[0008] Obtain velocity models and microseismic records of the monitoring area;

[0009] Dividing the velocity model into a plurality of grid points based on a preset grid size, performing travel time calculation on each of the grid points, and obtaining a travel time table containing theoretical travel time data;

[0010] Perform dynamic time warping calculation on every two seismic records in the microseismic records to obtain the time shift of the corresponding sampling points between every two microseismic records, pick up the first arrival time of the seismic event, and build a time shift sequence of each sampling point relative to the first arrival time based on the time shift of the first arrival time, calculate the correlation coefficient between seismic traces based on the time shift sequence, and obtain a time shift correlation trace set containing actual observed travel time difference information;

[0011] Acquire the imaging section corresponding to each grid point based on the time-shift correlation gather and the travel time table;

[0012] A final imaging map representing the earthquake source position is constructed based on the imaging profiles corresponding to the grid points, and the earthquake source position is determined based on the final imaging map.

[0013] Preferably, the calculating travel time of each grid point specifically includes:

[0014] The theoretical travel time data from each grid point to each detector is calculated, and a travel time table is constructed based on the calculated theoretical travel time data.

[0015] Preferably, the time shift sequence of each sampling point relative to the first arrival time is constructed, and the specific calculation formula is:

[0016] τ TC (t i )=τ t (t i )-τ t (t fb ),t i =t first ,t max

[0017] In the formula, i is the corresponding sampling point of the sequence, t i is the time corresponding to the sampling point, t first to max is the time range corresponding to the sampling point, τ t (t i ) is the calculated time shift corresponding to each sampling point, τ t (t fb ) is the first arrival time of microseismic signal t fb The corresponding time shift, τ TC The calculated time t for each sampling point i The relative time shift.

[0018] Preferably, the correlation coefficient between seismic traces is calculated based on the time-shift sequence to obtain a time-shift correlation trace gather containing actual observed travel time difference information. The specific calculation formula is:

[0019]

[0020] Where i and j are detector channel numbers, τ is the time shift between the corresponding sampling points of the two channels, and τ TC is the relative time shift, u(t i ,i) is the time series corresponding to detector i, u(t i +τ+τ TC (t i ), j) is the relative time shift sequence calculated by detector j relative to detector i, t i is the time corresponding to the sampling point, t first to max is the time range of the sampling point, C TC is the correlation coefficient between the two traces based on DTW.

[0021] Preferably, acquiring the imaging profile corresponding to each grid point based on the time-shift correlation gather and the travel time table specifically includes:

[0022] The time-shift correlation gathers are multiplied by the interferometric imaging operator containing the theoretical travel time data to perform interferometric imaging, and all the time-shift correlation gathers are superimposed at the same grid point to obtain a single grid point imaging profile.

[0023] Preferably, determining the earthquake source location based on the final imaging image specifically includes:

[0024] By determining the position corresponding to the maximum imaging value in the final imaging image, the earthquake source is obtained and the earthquake source positioning is achieved.

[0025] A microseismic interferometry positioning system based on time-shift correlation coefficient, comprising:

[0026] Data acquisition module, used to obtain velocity models and microseismic records of the monitoring area;

[0027] A travel time table calculation module, used to divide the velocity model into a plurality of grid points according to a preset grid size, perform travel time calculation on each of the grid points, and obtain a travel time table containing theoretical travel time data;

[0028] A time-shift correlation coefficient calculation module is used to perform dynamic time warping calculation on every two seismic records in the microseismic records, obtain the time shift of the corresponding sampling points between every two microseismic records, pick up the first arrival time of the seismic event, and use the time shift of the first arrival time as a reference to construct a time shift sequence of each sampling point relative to the first arrival time, calculate the correlation coefficient between seismic traces based on the time shift sequence, and obtain a time-shift correlation trace set containing actual observed travel time difference information;

[0029] A positioning module is used to obtain the imaging section corresponding to each grid point according to the time-shift correlation gather and the travel time table; construct a final imaging map representing the source position based on the imaging section corresponding to each grid point, and determine the source position based on the final imaging map.

[0030] An electronic device comprises a memory and a processor, wherein the memory is used to store a computer program, and the processor runs the computer program to enable the electronic device to execute the microseismic interference positioning method based on time-shift correlation coefficient.

[0031] A computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the microseismic interferometry positioning method based on time-shift correlation coefficient is implemented.

[0032] The technical effects of the present invention are:

[0033] The present invention uses a dynamic time warping (DTW) algorithm to construct a DTW time-shift correlation coefficient, multiplies the constructed time-shift correlation coefficient with an imaging operator, and superimposes all time-shift correlation gathers at the same grid point to obtain a single grid point imaging profile, and further superimposes all single grid point imaging profiles to obtain a final imaging effect that characterizes the source position. The position with the largest imaging value is the source, and the source location is achieved. The present invention constructs a time-shift correlation coefficient by using the DTW method, multiplies the time-shift correlation coefficient with an imaging operator, and superimposes to achieve imaging of the source position. The imaging effect is more focused and has better results than the results obtained by the traditional interference imaging method based on cross-correlation. The present invention can solve the phenomenon of reduced imaging resolution of the cross-correlation imaging operator caused by wave propagation effects and noise, and improve the positioning effect of the method under the influence of noise pollution. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0035] The drawings constituting a part of the present application are used to provide a further understanding of the present application. The illustrative embodiments and descriptions of the present application are used to explain the present application and do not constitute an improper limitation on the present application. In the drawings:

[0036] Figure 1 It is a schematic diagram of the implementation process in the embodiment of the present invention;

[0037] Figure 2 A schematic diagram of the time-shift correlation coefficient and the traditional mutual correlation coefficient proposed and constructed in an embodiment of the present invention;

[0038] Figure 3 Schematic diagram of the velocity model used in the forward modeling in the embodiment of the present invention;

[0039] Figure 4 A schematic diagram of noisy seismic data generated by adding -12dB Gaussian noise to original seismic data in an embodiment of the present invention;

[0040] Figure 5 The filtered seismic data obtained by performing 5-80Hz bandpass filtering on the data containing -12dB Gaussian noise in the embodiment of the present invention;

[0041] Figure 6 It is a schematic diagram of positioning results obtained by using the traditional cross-correlation interference positioning method and the interference positioning method based on the time-shift correlation coefficient proposed in the present invention respectively using the filtered data in an embodiment of the present invention. DETAILED DESCRIPTION

[0042] Various exemplary embodiments of the present invention will now be described in detail. This detailed description should not be considered as limiting the present invention, but should be understood as a more detailed description of certain aspects, features, and embodiments of the present invention.

[0043] It should be understood that the terms described in the present invention are only for describing special embodiments and are not intended to limit the present invention. In addition, for the numerical range in the present invention, it should be understood that each intermediate value between the upper and lower limits of the scope is also specifically disclosed. Each smaller range between the intermediate value in any stated value or stated range and any other stated value or intermediate value in the described range is also included in the present invention. The upper and lower limits of these smaller ranges can be independently included or excluded in the scope.

[0044] It will be apparent to those skilled in the art that various modifications and variations may be made to the specific embodiments of the present invention description without departing from the scope or spirit of the present invention. Other embodiments derived from the present invention description will be apparent to those skilled in the art. The present application description and examples are exemplary only.

[0045] The words “include,” “including,” “have,” “contain,” etc. used in this article are open-ended terms, meaning including but not limited to.

[0046] It should be noted that, in the absence of conflict, the embodiments and features in the embodiments of the present application can be combined with each other. The present application will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

[0047] Embodiment 1

[0048] like Figure 1 - Figure 6 As shown, in this embodiment, a microseismic interferometric positioning method based on time-shift correlation coefficient is provided, including: obtaining a velocity model and microseismic records of a monitoring area; dividing the velocity model into a number of grid points based on a preset grid size, performing travel time calculations on each of the grid points, and obtaining a travel time table containing theoretical travel time data; performing dynamic time regularization calculations on every two seismic records in the microseismic records, obtaining the time shift of the corresponding sampling points between every two microseismic records, picking up the first arrival time of the earthquake event, and constructing a time shift sequence of each sampling point relative to the first arrival time based on the time shift of the first arrival time, calculating the correlation coefficient between seismic traces based on the time shift sequence, and obtaining a time-shift correlation track set containing actual observed travel time difference information; obtaining an imaging profile corresponding to each grid point based on the time-shift correlation track set and the travel time table; constructing a final imaging map representing the source position based on the imaging profile corresponding to each of the grid points, and determining the source position based on the final imaging map.

[0049] This embodiment uses a dynamic time warping (DTW) algorithm to construct a DTW time-shift correlation coefficient, multiplies the constructed time-shift correlation coefficient with an imaging operator, and superimposes all time-shift correlation traces at the same grid point to obtain a single grid point imaging profile, and further superimposes all single grid point imaging profiles to obtain a final imaging effect that characterizes the source position. The position with the largest imaging value is the source, and the source location is achieved. This embodiment uses the DTW method to construct a time-shift correlation coefficient, and uses the time-shift correlation coefficient to multiply and superimpose the imaging operator to achieve imaging of the source position. The imaging effect is more focused and better than the result obtained by the traditional interference imaging method based on cross-correlation. This embodiment can solve the phenomenon of reduced imaging resolution of the cross-correlation imaging operator caused by wave propagation effects and noise, and improve the positioning effect of the method under the influence of noise pollution.

[0050] The specific technical solutions adopted in this embodiment are as follows:

[0051] Step 1: Calculate the travel time table. Divide the velocity model of the monitored area into several grid points according to a certain grid size, regard each grid point as a potential earthquake source, calculate the theoretical travel time from each grid point to each detector, and obtain the corresponding travel time table. Optionally, in order to efficiently calculate the travel time, this embodiment uses the open source Python package fteikpy to quickly solve the Eikonal equation. Compared with the prior art, the use of fteikpy significantly reduces the computational complexity while ensuring the accuracy of the solution;

[0052] Step 2: Calculate the time-shift correlation coefficient. Perform DTW calculation between every two seismic traces. The calculation result of the DTW algorithm is the time shift of the corresponding sampling points between every two seismic records. From this, the time shift τ of each sampling point of the waveform corresponding to the microseismic event can be obtained. t . Pick the first arrival time t of the corresponding earthquake event first , with the first arrival time t picked up first The corresponding time shift is τ fb , construct the waveform about τ fb Optionally, in order to pick up the first arrival time of the corresponding earthquake event, this embodiment uses the EventPicker Python package. Through the EventPicker package, the first arrival time t can be accurately picked up. fb On the basis of picking up the first arrival time and taking the time shift corresponding to the first arrival time as the benchmark, the relative time shift sequence of the waveform with respect to time is constructed as follows:

[0053] τ TC (t i )=τ t (ti )-τ t (t fb ),t i =t first ,t max

[0054] In the formula, i is the corresponding sampling point of the sequence, t i is the time corresponding to the sampling point, t first to max For the range to which it belongs, τ t (t i ) is the calculated time shift corresponding to each sampling point, τ t (t fb ) is the first arrival time of microseismic signal t fb The corresponding time shift, τ TC The calculated time t for each sampling point i The relative time shift.

[0055] Calculate the correlation coefficient between the relative time-shifted seismic trace and another seismic trace to obtain the time-shift correlation trace set containing the actual observed travel time difference information:

[0056]

[0057] Where i and j are detector channel numbers, τ is the time shift between the corresponding sampling points of the two channels, and τ TC is the relative time shift, u(t i ,i) is the time series corresponding to detector i, u(t i +τ+τ TC (t i ), j) is the relative time shift sequence calculated by detector j relative to detector i, t i is the time corresponding to the sampling point, t first to max is the time range of the sampling point, C TC is the correlation coefficient between the two traces based on DTW.

[0058] Step 3: Get a single-gather imaging profile. Multiply the time-shift correlation gathers with the interferometric imaging operator containing theoretical travel time difference information to perform interferometric imaging, and superimpose all time-shift correlation gathers at the same grid point to obtain a single-grid point imaging profile:

[0059]

[0060] Where x is the source position vector, δ is the Dirac Delta function, τ is the time, N is the number of geophones, and τ i,x , τ j,x is the theoretical travel time from source x to detectors i and j, S TC(x,i) is the imaging profile of a single grid point. The interferometric imaging operator uses the travel time difference from the source to the receiver pair {i,j}.

[0061] Step 4: Superimpose all single grid point imaging profiles to obtain the final imaging effect that represents the source location. The location with the largest imaging value is the source, thus achieving source location.

[0062]

[0063] Where N is the number of detectors, S TC (x) is the final imaging result.

[0064] Compared with the traditional interference imaging method based on the mutual correlation coefficient, the time-shift correlation coefficient is constructed in Formula 2 of this embodiment. The comparison between the time-shift correlation coefficient and the mutual correlation coefficient is as follows: Figure 2 shown. Figure 6 The comparison between the microseismic interferometric imaging positioning effect based on the time-shift correlation coefficient of this embodiment and the conventional interferometric imaging positioning effect based on the mutual correlation coefficient shows that the imaging effect based on the time-shift correlation coefficient is greatly improved.

[0065] A microseismic interferometry positioning system based on time-shift correlation coefficient, comprising:

[0066] Data acquisition module, used to obtain velocity models and microseismic records of the monitoring area;

[0067] A travel time table calculation module, used to divide the velocity model into a plurality of grid points according to a preset grid size, perform travel time calculation on each of the grid points, and obtain a travel time table containing theoretical travel time data;

[0068] A time-shift correlation coefficient calculation module is used to perform dynamic time warping calculation on every two seismic records in the microseismic records, obtain the time shift of the corresponding sampling points between every two microseismic records, pick up the first arrival time of the seismic event, and use the time shift of the first arrival time as a reference to construct a time shift sequence of each sampling point relative to the first arrival time, calculate the correlation coefficient between seismic traces based on the time shift sequence, and obtain a time-shift correlation trace set containing actual observed travel time difference information;

[0069] A positioning module is used to obtain the imaging section corresponding to each grid point according to the time-shift correlation gather and the travel time table; construct a final imaging map representing the source position based on the imaging section corresponding to each grid point, and determine the source position based on the final imaging map.

[0070] An electronic device comprises a memory and a processor, wherein the memory is used to store a computer program, and the processor runs the computer program to enable the electronic device to execute the microseismic interference positioning method based on time-shift correlation coefficient.

[0071] A computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the microseismic interferometry positioning method based on time-shift correlation coefficient is implemented.

[0072] The above is only a preferred specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed in the present application should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be based on the protection scope of the claims.

Claims

1. A microseismic interferometry positioning method based on time-shift correlation coefficient, characterized in that: include: Obtain velocity models and microseismic records of the monitoring area; Dividing the velocity model into a plurality of grid points based on a preset grid size, performing travel time calculation on each of the grid points, and obtaining a travel time table containing theoretical travel time data; Perform dynamic time warping calculation on every two seismic records in the microseismic records to obtain the time shift of the corresponding sampling points between every two microseismic records, pick up the first arrival time of the seismic event, and build a time shift sequence of each sampling point relative to the first arrival time based on the time shift of the first arrival time, calculate the correlation coefficient between seismic traces based on the time shift sequence, and obtain a time shift correlation trace set containing actual observed travel time difference information; Acquire the imaging section corresponding to each grid point based on the time-shift correlation gather and the travel time table; A final imaging map representing the earthquake source position is constructed based on the imaging profiles corresponding to the grid points, and the earthquake source position is determined based on the final imaging map.

2. The microseismic interferometry positioning method based on time-shift correlation coefficient according to claim 1, characterized in that: The calculation of travel time for each grid point specifically includes: The theoretical travel time data from each grid point to each detector is calculated, and a travel time table is constructed based on the calculated theoretical travel time data.

3. The microseismic interferometry positioning method based on time-shift correlation coefficient according to claim 1, characterized in that: The specific calculation formula for constructing the time shift sequence of each sampling point relative to the first arrival time is: τ TC (t i )=τ t (t i )-τ t (t fb ),t i =t first ,t max In the formula, i is the corresponding sampling point of the sequence, t i is the time corresponding to the sampling point, t first to max is the time range corresponding to the sampling point, τ t (t i ) is the calculated time shift corresponding to each sampling point, τ t (t fb ) is the first arrival time of microseismic signal t fb The corresponding time shift, τ TC The calculated time t for each sampling point i The relative time shift.

4. The microseismic interferometry positioning method based on time-shift correlation coefficient according to claim 1, characterized in that: The correlation coefficient between seismic traces is calculated based on the time-shift sequence to obtain a time-shift correlation trace set containing actual observed travel time difference information. The specific calculation formula is: Where i and j are detector channel numbers, τ is the time shift between the corresponding sampling points of the two channels, and τ TC is the relative time shift, u(t i ,i) is the time series corresponding to detector i, u(t i +τ+τ TC (t i ), j) is the relative time shift sequence calculated by detector j relative to detector i, t i is the time corresponding to the sampling point, t first to max is the time range of the sampling point, C TC is the correlation coefficient between the two traces based on DTW.

5. The microseismic interferometry positioning method based on time-shift correlation coefficient according to claim 1, characterized in that: The obtaining of the imaging profile corresponding to each grid point based on the time-shift correlation gather and the travel time table specifically includes: The time-shift correlation gathers are multiplied by the interferometric imaging operator containing the theoretical travel time data to perform interferometric imaging, and all the time-shift correlation gathers are superimposed at the same grid point to obtain a single grid point imaging profile.

6. The microseismic interferometry positioning method based on time-shift correlation coefficient according to claim 1, characterized in that: Determining the earthquake source location based on the final imaging image specifically includes: By determining the position corresponding to the maximum imaging value in the final imaging image, the earthquake source is obtained and the earthquake source positioning is achieved.

7. A microseismic interferometry positioning system based on time-shift correlation coefficient, characterized in that: include: Data acquisition module, used to obtain velocity models and microseismic records of the monitoring area; A travel time table calculation module, used to divide the velocity model into a plurality of grid points according to a preset grid size, perform travel time calculation on each of the grid points, and obtain a travel time table containing theoretical travel time data; A time-shift correlation coefficient calculation module is used to perform dynamic time warping calculation on every two seismic records in the microseismic records, obtain the time shift of the corresponding sampling points between every two microseismic records, pick up the first arrival time of the seismic event, and use the time shift of the first arrival time as a reference to construct a time shift sequence of each sampling point relative to the first arrival time, calculate the correlation coefficient between seismic traces based on the time shift sequence, and obtain a time-shift correlation trace set containing actual observed travel time difference information; A positioning module is used to obtain the imaging section corresponding to each grid point according to the time-shift correlation gather and the travel time table; construct a final imaging map representing the source position based on the imaging section corresponding to each grid point, and determine the source position based on the final imaging map.

8. An electronic device, characterized in that: It comprises a memory and a processor, wherein the memory is used to store a computer program, and the processor runs the computer program to enable the electronic device to execute a microseismic interference positioning method based on a time-shift correlation coefficient according to any one of claims 1-6.

9. A computer-readable storage medium, characterized in that: It stores a computer program, which, when executed by a processor, implements a microseismic interferometry positioning method based on time-shift correlation coefficient as described in any one of claims 1 to 6.

Images