CN111352153B - Microseism interference positioning method based on instantaneous phase cross-correlation weighting - Google Patents
Microseism interference positioning method based on instantaneous phase cross-correlation weighting Download PDFInfo
- Publication number
- CN111352153B CN111352153B CN202010207285.1A CN202010207285A CN111352153B CN 111352153 B CN111352153 B CN 111352153B CN 202010207285 A CN202010207285 A CN 202010207285A CN 111352153 B CN111352153 B CN 111352153B
- Authority
- CN
- China
- Prior art keywords
- correlation
- waveform
- cross
- component
- time
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 33
- 238000003384 imaging method Methods 0.000 claims abstract description 56
- 239000000126 substance Substances 0.000 claims description 9
- 230000002123 temporal effect Effects 0.000 claims 1
- 238000004364 calculation method Methods 0.000 abstract description 15
- 230000007246 mechanism Effects 0.000 abstract description 5
- 230000008859 change Effects 0.000 abstract description 4
- 230000003252 repetitive effect Effects 0.000 abstract 1
- 238000010586 diagram Methods 0.000 description 14
- 238000001514 detection method Methods 0.000 description 5
- 238000012544 monitoring process Methods 0.000 description 5
- 238000012545 processing Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000005012 migration Effects 0.000 description 2
- 238000013508 migration Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000035699 permeability Effects 0.000 description 1
- 230000010287 polarization Effects 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/288—Event detection in seismic signals, e.g. microseismics
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/301—Analysis for determining seismic cross-sections or geostructures
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- Acoustics & Sound (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Business, Economics & Management (AREA)
- Emergency Management (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention discloses a microseism interference positioning method based on instantaneous phase cross-correlation weighting, which comprises the steps of obtaining N seismic waves of a microseism event through N detectors, and carrying out component on the seismic waves; combining the detectors one by one to form a non-repetitive detector pair; calculating a time domain cross-correlation waveform and an instantaneous phase cross-correlation waveform in each component; weighting each time domain cross-correlation waveform in each component and the corresponding instantaneous phase cross-correlation waveform to obtain a corresponding weighted cross-correlation waveform; in each component, offset superposition is carried out on the corresponding weighted cross-correlation waveform of each channel by utilizing the travel time difference of the theoretical direct wave to obtain an interference imaging section; and determining the position of the microseismic event according to the interference imaging section in the plurality of components multiplied to obtain a comprehensive interference imaging section. By the method, the problem of polarity change of the micro-seismic waveform caused by a complex seismic source mechanism is solved, and meanwhile, the influence of high-amplitude noise on time domain cross correlation and subsequent superposition calculation is reduced.
Description
Technical Field
The invention relates to the technical field of microseism monitoring and positioning, in particular to a microseism interference positioning method based on instantaneous phase cross-correlation weighting.
Background
With the increasing energy of the world, the emphasis of oil and gas exploration has gradually shifted from conventional oil and gas fields to unconventional oil and gas fields. The unconventional oil and gas fields have the characteristics of compactness, low permeability and the like, and are mainly explored and developed by adopting a hydraulic fracturing technology. The microseism monitoring technology is one of the most effective methods for realizing the real-time monitoring of the hydraulic fracturing effect. The basic flow is that when the fracturing operation is carried out in a reservoir through a fracturing well, the ground stress changes to form cracks, the seismic wave signals are received by a detector array and then can be used for inversion to obtain the position of a seismic source along with weak acoustic emission or seismic wave signals, and then the crack form is determined according to the distribution of the position of the seismic source, so that the construction design and the effect of the fracturing are optimized.
However, the traditional travel time inversion type microseism positioning algorithm needs to carry out first-arrival picking on effective microseism events, so that a larger picking error exists, and particularly for ground microseism monitoring with extremely low signal-to-noise ratio, the first-arrival picking on the microseism events can not be directly carried out generally. Accordingly, in recent years, some waveform stacking type automatic microseismic positioning methods that do not require first arrival pickup have been studied and developed. The positioning method may not obtain an accurate seismic source position because the influence of a complex seismic source mechanism on the amplitude and the phase of the waveform is not considered. Meanwhile, the method needs to search the origin time and the origin position of the seismic source at the same time, which limits the positioning accuracy and efficiency to a certain extent.
Disclosure of Invention
In order to overcome the defects in the prior art, the embodiment of the invention provides a microseism interference positioning method based on instantaneous phase cross-correlation weighting, which comprises the following steps:
acquiring N seismic waves of the microseism event through N detectors, and performing component on the seismic waves; combining N detectors one by one to form a detector pair with non-repeated H pairs, wherein H is N (N-1)/2, and N is a positive integer greater than or equal to 2; calculating H-channel time domain cross-correlation waveforms and H-channel instantaneous phase cross-correlation waveforms in each component; weighting each channel of time domain cross-correlation waveform in each component and the corresponding instantaneous phase cross-correlation waveform to obtain H channels of weighted cross-correlation waveforms in the corresponding components; in each component, carrying out offset superposition on the corresponding weighted cross-correlation waveform of each channel by using the travel time difference of the theoretical direct wave to obtain an interference imaging profile of each component; multiplying interference imaging sections in the multiple components to obtain a comprehensive interference imaging section; the location of the microseismic event is determined from the synthetic interferometric imaging profile.
In one embodiment, determining the location of the microseismic event based on the synthetic interferometric imaging profile includes determining the location of the microseismic event based on a maximum value of the synthetic interferometric imaging profile, the maximum value corresponding to the location of the microseismic event.
In one embodiment, the time-domain cross-correlation waveform is calculated as:
wherein v iskRepresenting a horizontal component vxOr the vertical component vz，
The waveform of the n-th micro-seismic event and the waveform of the m-th micro-seismic event are in vkA time-discrete sequence of time-domain cross-correlations over the components,
represents vkA time discrete sequence of the nth trace of the micro-seismic event waveform on the component,
represents vkThe time-discrete sequence of the mth trace of the microseismic event waveform on the component is passed through a shifted time-discrete sequence of length τ, TxaxRepresenting the discrete time length of the single trace microseismic event waveform before the time domain cross-correlation operation.
In one embodiment, the instantaneous phase cross-correlation waveform is calculated as:
wherein the content of the first and second substances,
is the cross-correlation of the instantaneous phase of the nth trace micro-seismic event waveform with the instantaneous phase of the mth trace micro-seismic event waveform at vkTime-discrete sequences over components, vkRepresenting a horizontal component vxOr the vertical component vz，
The instantaneous phase of the waveform representing the nth microseismic event is at vkA time-discrete sequence of the components is,
a time discrete sequence representing the instantaneous phase of the mth trace of the microseismic event waveform is shifted by a length τ in vkTime discrete sequences on components, TmaxRepresenting instantaneous phase cross-correlationThe discrete time length of the waveform of the previous single-track microseismic event, and q is an index for controlling the curve sharpness degree of the instantaneous phase cross-correlation curve from a high coherence position to a low coherence position.
In one embodiment, the weighted cross-correlation waveform is calculated as:
wherein the content of the first and second substances,
is shown at vkThe nth track and the mth track on the component are subjected to time-domain cross-correlation time-discrete sequences after instantaneous phase cross-correlation weighting.
In one embodiment, the travel time difference of the theoretical direct wave is the difference between the travel time of the nth theoretical direct wave and the travel time of the mth theoretical direct wave in the N seismic waves; establishing a speed model for the region to be detected, and obtaining the theoretical direct wave travel time from each point in the region to be detected to each detector according to the speed model; the theoretical direct wave travel time comprises theoretical direct longitudinal wave travel time and theoretical direct transverse wave travel time.
In one embodiment, in each component, offset stacking is performed on the corresponding weighted cross-correlation waveform of each channel by using the theoretical direct wave travel time difference, so as to obtain an interference imaging profile of each component, which specifically includes:
in the vertical component, offset superposition is carried out on the corresponding cross-correlation waveform after weighting by using the travel time difference of the theoretical direct longitudinal wave, so as to obtain an interference imaging section of the longitudinal wave in the vertical component; in the horizontal component, offset superposition is carried out on the corresponding cross-correlation waveform after weighting by using the travel time difference of the theoretical direct transverse wave, so as to obtain an interference imaging section of the transverse wave in the horizontal component.
In one embodiment, the seismic waves are subjected to components, and when the components are used for two-dimensional positioning, a vertical component and a horizontal component are obtained; when used for three-dimensional positioning, three components are obtained that are perpendicular to each other, the three components including: a vertical component, a first horizontal component, and a second horizontal component.
The embodiment of the invention has the advantages that: the problem of polarity change of micro-seismic waveforms caused by a complex seismic source mechanism is solved by multiplying any two instantaneous phase cross-correlation waveforms by corresponding time domain cross-correlation waveforms for weighting, meanwhile, the instantaneous phase cross-correlation curves are unbiased in amplitude, namely, the value of the instantaneous phase cross-correlation curves is not influenced by amplitude, the influence of high-amplitude noise on time domain cross-correlation and subsequent superposition calculation is reduced, the seismic source seismic origin time is decoupled from a micro-seismic positioning algorithm no matter in the time domain cross-correlation calculation or the instantaneous phase cross-correlation calculation, the influence of the seismic origin time on the micro-seismic positioning result is avoided, the positioning precision and the resolution are improved, and reliable basis and guidance can be provided for hydraulic fracture monitoring.
Drawings
FIG. 1 is a flow chart of a microseism interference positioning method based on instantaneous phase cross-correlation weighting according to an embodiment of the invention;
FIG. 2 is a flowchart of a method for calculating travel time of a theoretical direct arrival according to an embodiment of the present invention;
FIG. 3(a) is a schematic diagram illustrating components of the vibration velocity of seismic waves in the vertical direction according to an embodiment of the present invention;
FIG. 3(b) is a schematic diagram illustrating components of the seismic wave vibration velocity in the horizontal direction in the embodiment of the present invention;
FIG. 4 is a schematic diagram of a two-dimensional planar discrete model of a region to be measured according to an embodiment of the present invention;
FIG. 5(a) is a time-domain cross-correlation waveform diagram of a first trace waveform and other traces in the vertical direction according to an embodiment of the present invention;
FIG. 5(b) is a time-domain cross-correlation waveform diagram of the first trace waveform and other traces in the horizontal direction in the embodiment of the present invention;
FIG. 6(a) is a waveform diagram of instantaneous phase cross-correlation of a first waveform with other waveforms in the vertical direction according to an embodiment of the present invention;
FIG. 6(b) is a waveform diagram of instantaneous phase cross-correlation of a first waveform with other waveforms in the horizontal direction according to an embodiment of the present invention;
FIG. 7(a) is a weighted cross-correlation waveform of a first trace waveform and other trace waveforms in the vertical direction according to an embodiment of the present invention;
FIG. 7(b) is a weighted cross-correlation waveform of the first trace and other traces in the horizontal direction according to an embodiment of the present invention;
FIG. 8(a) is a two-dimensional microseismic source interferometric imaging profile of an embodiment of the present invention;
fig. 8(b) is a three-dimensional micro seismic source interferometric imaging diagram according to the embodiment of the invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the embodiments 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 drawings in the embodiments of the present invention, and it is obvious that the described embodiments are some, but not all, embodiments of the present invention. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
FIG. 1 is a flow chart of a microseismic interferometric positioning method based on instantaneous phase cross-correlation weighting according to an embodiment of the present invention, as shown in FIG. 1. The method comprises the following steps:
step S110: and acquiring the recording information of the microseism event of the area to be detected through the N detectors. Wherein the recorded information comprises N seismic waves corresponding to the N detectors.
And carrying out component on the N seismic waves. If the area to be measured is a two-dimensional plane area, performing two components on the seismic waves to obtain a vertical component and a horizontal component; and if the area to be measured is a three-dimensional area, performing three components on the seismic waves to obtain a vertical component, a first horizontal component and a second horizontal component, wherein the three components are vertical in pairs.
And combining N detectors one by one to obtain a detector pair with no repeated H pairs, wherein H is N (N-1)/2, and N is a positive integer greater than or equal to 2.
Step S120: and calculating H-channel time domain cross-correlation waveforms and H-channel instantaneous phase cross-correlation waveforms in each component.
Take the detection of a two-dimensional planar area as an example. Carrying out horizontal direction and vertical direction components on the seismic waves to obtain a horizontal component vxOr the vertical component vz。
At vzIn the method, time domain cross-correlation calculation and instantaneous phase cross-correlation calculation are respectively carried out on each detector pair to obtain vzH-channel time domain cross-correlation waveform and H-channel instantaneous phase cross-correlation waveform; and each cross-correlation waveform comprises the information of the direct transverse wave travel time difference and the direct longitudinal wave travel time difference. At vxIn the method, time domain cross-correlation calculation and instantaneous phase cross-correlation calculation are respectively carried out on each detector pair to obtain vxH-channel time domain cross-correlation waveform and H-channel instantaneous phase cross-correlation waveform; and each cross-correlation waveform comprises the information of the direct transverse wave travel time difference and the direct longitudinal wave travel time difference.
Specifically, the time-domain cross-correlation calculation formula is:
wherein v iskRepresenting a horizontal component vxOr the vertical component vz，
The waveform of the n-th micro-seismic event and the waveform of the m-th micro-seismic event are in vkA time-discrete sequence of time-domain cross-correlations over the components,
represents vkA time discrete sequence of the nth trace of the micro-seismic event waveform on branch x,
represents vkThe time-discrete sequence of the mth trace of the microseismic event waveform on the component is passed through a shifted time-discrete sequence of length tau, TmaxRepresenting the discrete time length of the single trace microseismic event waveform before the time domain cross-correlation operation.
The instantaneous phase cross-correlation calculation formula is as follows:
wherein the content of the first and second substances,
is the cross-correlation of the instantaneous phase of the nth trace micro-seismic event waveform with the instantaneous phase of the mth trace micro-seismic event waveform at vkTime-discrete sequences over components, vkRepresenting a horizontal component vxOr the vertical component vz，
The instantaneous phase of the waveform representing the nth microseismic event is at vkA time-discrete sequence of the components is,
a time discrete sequence representing the instantaneous phase of the mth trace of the microseismic event waveform is shifted by a length τ in vkTime discrete sequences on components, TmaxRepresenting the discrete time length of the waveform of a single trace microseismic event before the instantaneous phase cross-correlation operation, and q is an index controlling the sharpness of the instantaneous phase cross-correlation curve from a high coherence place to a low coherence place, namely, a higher q value can cause the attenuation of a signal with lower coherence to be larger, thereby being beneficial to better highlighting the coherent signal and adjusting the sensitivity of coherence measurement.
In the above formula
Obtained by the following formula:
wherein the content of the first and second substances,
is represented by vkTime discrete sequence of nth trace microseismic event waveform on component
And impulse response
Making Hilbert transform obtained by convolution on the time discrete sequence;
is represented by vkThe real part of the analytic sequence of the waveform time dispersion of the nth trace of the microseismic event on the component is vkTime discrete sequence of nth trace microseismic event waveform on component
Imaginary part at vkHilbert transform of nth trace micro-seismic event waveform time discrete sequence on component
Is represented by vkA time-discrete sequence of instantaneous phases in the components obtained by dividing the imaginary part of the time-discrete analytic sequence by the real part and by taking the inverse tangent value, vkRepresenting a horizontal component vxOr the vertical component vz. Meanwhile, in order to obtain the instantaneous phase varying between-pi and pi, the invention adopts an angle function to calculate the instantaneous phase in matlab.
The above calculation process is exemplified by a two-dimensional detection area. And if the detection area is a three-dimensional area, respectively calculating H-channel time domain cross-correlation waveforms and H-channel instantaneous phase cross-correlation waveforms of a first horizontal component and a second horizontal component which are perpendicular to each other in the horizontal components.
Step S130: and weighting each channel of time domain cross-correlation waveform in each component and the corresponding instantaneous phase cross-correlation waveform to obtain the H channel weighted cross-correlation waveform in the corresponding component.
Due to the influence of complex source mechanisms, the polarity of the detected waveform may be positive or negative for the same seismic phase (P-wave or S-wave). If the change of polarity is not considered, the direct offset superposition imaging of the time domain cross-correlation curve can make the positive value and the negative value mutually cancel each other, the peak value information in the waveform is lost, and the positioning result is influenced. Therefore, each instantaneous phase cross-correlation curve is multiplied by the corresponding time domain cross-correlation curve for weighting operation, and when the phase size and the polarity are consistent, the instantaneous phase cross-correlation value is a positive peak value; when the phases are consistent in size but opposite in polarity, the instantaneous phase cross-correlation value is a negative peak value; the absolute value of the cross-correlation value of other randomly varying instantaneous phases is smaller than the absolute value of the cross-correlation value when the phases are consistent. The method can not only achieve the purpose of weighting the time-domain cross-correlation curve, endow a high weight at the position with high correlation and endow a low weight at the position with low correlation, but also correct the polarity corresponding to the peak value in the time-domain cross-correlation curve and convert the polarity into a positive value, thereby avoiding the influence on superposition. And because the calculation of instantaneous phase cross-correlation is amplitude unbiased, the influence of high-amplitude noise on the positioning result can be effectively inhibited.
Specifically, the weighted time-domain cross-correlation calculation formula is:
wherein the content of the first and second substances,
is shown at vkThe nth track and the mth track on the component are subjected to time-domain cross-correlation time-discrete sequences after instantaneous phase cross-correlation weighting.
Step S140: in each component, the theoretical direct wave travel time difference is utilized to perform offset superposition on the corresponding weighted cross-correlation waveform of each channel, so as to obtain an interference imaging section of each component.
The theoretical direct wave travel time difference used here is calculated by establishing a speed model for the region to be measured. Specifically, the method for calculating the travel time of the theoretical direct arrival, as shown in fig. 2, includes the following steps:
step S210: and establishing a corresponding speed model according to the actual condition of the area to be measured.
Step S220: and solving an equation of an equation according to the speed model or adopting a ray tracing method to obtain the theoretical direct wave travel time of each point in the region to be detected reaching each detector.
The theoretical direct wave travel time obtained in step S220 includes theoretical direct longitudinal wave travel time and theoretical direct transverse wave travel time. Furthermore, a corresponding longitudinal wave travel time table and a transverse wave travel time table can be established for storing the theoretical direct wave travel time. When the specific theoretical direct wave travel time difference is calculated, the theoretical direct wave travel time difference can be calculated by storing the theoretical direct wave travel time in the travel time table by using the two detectors corresponding to each cross-correlation waveform.
Note that, in step S120, only the attention v is focused onzMedium direct longitudinal wave running time difference information and vxAnd (5) thinking of the travel time difference information of the intermediate direct transverse wave. P-wave at v due to the influence of polarization characteristics of P-wave and S-wavezStrong medium energy, S wave at vxThe medium energy is stronger. Therefore, in the examples of this application, at vxIn the time domain cross-correlation waveform and the instantaneous phase cross-correlation waveform, only the peak value corresponding to the S wave travel time difference is utilized to obtain an interference imaging section of the S wave; at vzThe time domain cross-correlation waveform and the instantaneous phase cross-correlation waveform of the P wave interference imaging section are obtained only by utilizing the peak value corresponding to the P wave travel time difference. And areThe final micro-seismic interference imaging section is obtained by multiplying the interference imaging sections of the two components, so that the influence of noise is inhibited, the artifact is reduced, and the energy in the imaging section is focused more.
Taking the example of detecting a two-dimensional planar area, we are at vxIn the method, the corresponding cross-correlation waveform after weighting is subjected to offset superposition by using the travel time difference of the theoretical direct transverse wave to obtain the correlation vxInterference imaging profile of
Wherein the content of the first and second substances,
when the direct transverse wave of the nth detector travels,
the time of the direct transverse wave of the mth wave detector.
At vzIn the method, the corresponding cross-correlation waveform after weighting is subjected to offset superposition by utilizing the travel time difference of the theoretical direct longitudinal wave to obtain the correlation vzInterference imaging profile of
Wherein the content of the first and second substances,
the time of the direct longitudinal wave of the nth detector,
is the direct longitudinal wave travel of the mth wave detectorThen (c) is performed.
V of seismic waves if the detection area is a three-dimensional areakThe components may include a vertical component vzA first horizontal component vx1And a second horizontal component vx2. Wherein, with respect to vzInterference imaging profile of
The calculation method is the same as that for detecting the two-dimensional plane area, and the description is omitted here. At vx1In the method, the corresponding cross-correlation waveform after weighting is subjected to offset superposition by using the travel time difference of the theoretical direct transverse wave to obtain the correlation vx1Interference imaging profile of
At vx2In the method, the corresponding cross-correlation waveform after weighting is subjected to offset superposition by using the travel time difference of the theoretical direct transverse wave to obtain the correlation vx2Interference imaging profile of
Step S150: multiplying interference imaging profiles in a plurality of components to obtain a comprehensive interference imaging profile
The location of the microseismic event is determined from the synthetic interferometric imaging profile.
And determining the position of the micro-seismic event according to the maximum value of the comprehensive interference imaging section, wherein the maximum value corresponds to the position of the micro-seismic event.
In a specific embodiment, the two-dimensional planar area to be measured shown in fig. 4 is discretized to obtain an area to be measured composed of a plurality of grid areas. The area is a two-dimensional plane rectangular area with the length of 500 meters and the depth of 500 meters, and the discretely divided grid is a rectangular area with the size of 2.5mx2.5 m. Setting the second place at the coordinate 0 point of the earth surface of the region to be measuredDetector M1Then, one detector is arranged at intervals of 10 meters until 51 detectors M are arranged at the position of 500 meters of coordinatesiI is 1, 2 … 51. By using a two-dimensional horizontal three-layer velocity model, a microseism event occurs at the position of position coordinates (250m ) in the area to be measured.
At the source of the micro-seismic event, a source of concentration force is employed that is excited in the horizontal direction. A waveform diagram as shown in fig. 3 is obtained. Fig. 3(a) shows a component vz of the seismic wave vibration velocity in the vertical direction z, and it can be seen that a wave starting between 100ms and 200ms is a P wave, the polarity of the P wave is symmetric with respect to the abscissa of the seismic source position, and the polarity of the P wave changes from positive to negative when viewed from left to right in the direction of the abscissa; the wave with the starting point between 200ms and 300ms is an S wave, the polarity of the S wave is also symmetrically distributed about the abscissa of the seismic source position, the polarity of the S wave changes from negative to positive when viewed from left to right in the direction of the abscissa, the energy of the P wave and the energy of the S wave are almost equal, and the wave detector (namely the 26 th wave detector M) is arranged right above the seismic source position26) The received P-wave and S-wave have the weakest waveform energy. FIG. 3(b) is a component v of the vibration velocity of the seismic wave in the horizontal direction xxIt can be seen that the wave starting between 100ms and 200ms is a P-wave, the polarity of the P-wave is positive, and the waveform energy received by the detector directly above the source location is weakest; the wave with the starting point between 200ms and 300ms is an S wave, the polarity of the S wave is unchanged and is always positive, and the energy of the S wave is obviously greater than that of the P wave.
In the case of a micro-seismic event,
there are a total of 1275 non-repeating detector pairs of 51 detectors.
Using the formula:
separately calculate vz1275 time-domain cross-correlation waveform sum v inx1275 time-domain cross-correlation waveforms in (a). As shown in FIG. 5, FIG. 5(a) shows vzTime domain cross-correlation waveform diagrams of the first channel waveform and other channels of waveforms; FIG. 5(b) shows vxThe time domain cross-correlation waveform diagram of the first trace waveform and other traces.
Using the formula:
separately calculate vz1275 instantaneous phase cross correlation waveform sum vx1275 instantaneous phase cross correlation waveforms in (a). As shown in FIG. 6, FIG. 6(a) shows vzInstantaneous phase cross-correlation waveform diagrams of the first waveform and other waveforms; FIG. 6(b) shows vxThe instantaneous phase cross-correlation waveform of the first trace and other traces.
Using the formula:
are each at vzAnd vxIn the method, each instantaneous phase cross-correlation waveform is multiplied by a corresponding time domain cross-correlation waveform to carry out weighting processing to obtain vz1275 weighted cross-correlation waveform sum vx1275 weighted cross-correlation waveforms in (a). As shown in FIG. 7, FIG. 7(a) shows vzA weighted cross-correlation waveform diagram of the first channel waveform and other channels of waveforms; FIG. 7(b) shows vxThe weighted cross-correlation waveform of the first trace and the other traces.
At vzMiddle, utilizing throughThe travel time difference of P wave travel time respectively carries out offset superposition on 1275 weighted cross-correlation waveforms to obtain an interference imaging section of the P wave
At vxRespectively carrying out offset superposition on 1275 weighted cross-correlation waveforms by using the travel time difference of the travel time of the direct S wave to obtain an interference imaging section of the S wave
V is to bexInterference imaging profile formed with respect to a transverse wave
And vzInterference imaging profile formed about longitudinal wave
Multiplying to form the final micro-seismic interference imaging section
As shown in fig. 8, wherein fig. 8(a) is a two-dimensional micro-seismic source interference imaging section view; fig. 8(b) is a three-dimensional micro seismic source interference imaging diagram.
In fig. 8(a), the intensity of energy is represented by gray scale, where white represents higher energy and black represents lower energy. The maximum of the location energy, located on the abscissa 250m and the ordinate 250m, represents the location as the source location, corresponding to the location where the micro-seismic event initially set in this embodiment occurred. From the results, the absolute error of the final imaging from the preset source position is shown to be (0m, 0m) (the difference between the position of the actual measured microseismic event and the position of the preset microseismic event: (250m, 250m), (250m )).
The three-dimensional microseismic source interferogram in fig. 8(b) may show that the source position (x-axis 250m, z-axis 250m) corresponds to a distinct energy peak, and that artifacts are well suppressed in directions other than perpendicular to the detector array direction. For artifacts perpendicular to the detector array direction, the placement of the geophone array relative to the source results in the artifacts extending perpendicular to the geophone array direction, since the source location is determined by the intersection of one of the H hyperbolas.
According to the microseism interference imaging method based on instantaneous phase cross-correlation weighting, the microseism event does not need to be picked up in a first arrival travel time mode, and meanwhile, the earthquake-initiating time of the microseism event does not need to be searched, so that the accuracy and the resolution of seismic source positioning are obviously improved; performing time domain cross-correlation operation and instantaneous phase cross-correlation operation according to the recorded information of the microseism event in the detector, performing multiplication weighting processing on the time domain cross-correlation of corresponding channels by using the instantaneous phase cross-correlation of each channel in a horizontal component and a vertical component, performing migration stacking processing on cross-correlation waveforms weighted by the horizontal component by using direct transverse wave to time difference to obtain interference imaging sections related to transverse waves, performing migration stacking processing on the cross-correlation waveforms weighted by the vertical component by using direct longitudinal wave to time difference to obtain interference imaging sections related to longitudinal waves, and multiplying the interference imaging sections of each component to obtain final interference imaging sections. Because the instantaneous phase cross-correlation is calculation with unbiased amplitude, the instantaneous phase cross-correlation is used for weighting the time domain cross-correlation, so that the influence of the change of the waveform polarity caused by different seismic source mechanisms on the positioning result can be overcome, the influence of high-amplitude noise on the positioning result can be inhibited, and the detection capability of the weak micro-seismic event is improved; meanwhile, the resolution and the precision of imaging can be effectively improved by properly increasing the value of the index q, particularly for a seismic source at a deep layer; subsequent superposition processing is more favorable for enabling the positioning algorithm to be suitable for low signal-to-noise ratio data; meanwhile, the horizontal component transverse wave interference imaging section and the vertical component longitudinal wave interference imaging section are utilized to facilitate energy focusing of imaging results, and positioning accuracy and resolution are improved.
The above embodiments are provided to further explain the objects, technical solutions and advantages of the present invention in detail, it should be understood that the above embodiments are merely exemplary embodiments of the present invention and are not intended to limit the scope of the present invention, and any modifications, equivalents, improvements and the like made within the spirit and principle of the present invention should be included in the scope of the present invention.
Claims (8)
1. A microseism interference positioning method based on instantaneous phase cross-correlation weighting is characterized by comprising the following steps:
acquiring N seismic waves of the microseism event through N detectors, and performing component on the seismic waves; combining N detectors one by one to form a detector pair with non-repeated H pairs, wherein H is N (N-1)/2, and N is a positive integer greater than or equal to 2;
calculating H-channel time domain cross-correlation waveforms and H-channel instantaneous phase cross-correlation waveforms in each component;
weighting each channel of time domain cross-correlation waveform in each component and the corresponding instantaneous phase cross-correlation waveform to obtain H channels of weighted cross-correlation waveforms in the corresponding components;
in each component, carrying out offset superposition on the corresponding weighted cross-correlation waveform of each channel by using the travel time difference of the theoretical direct wave to obtain an interference imaging profile of each component;
multiplying interference imaging sections in the multiple components to obtain a comprehensive interference imaging section; the location of the microseismic event is determined from the synthetic interferometric imaging profile.
2. The method of claim 1, wherein determining the location of the microseismic event based on the synthetic interferometric imaging profile comprises determining the location of the microseismic event based on a maximum value of the synthetic interferometric imaging profile, the maximum value corresponding to the location of the microseismic event.
3. The method of claim 1, wherein the time-domain cross-correlation waveform is calculated as:
wherein v iskRepresenting a horizontal component vxOr the vertical component vz，
The waveform of the n-th micro-seismic event and the waveform of the m-th micro-seismic event are in vkA time-discrete sequence of time-domain cross-correlations over the components,
represents vkA time discrete sequence of the nth trace of the micro-seismic event waveform on the component,
represents vkThe time-discrete sequence of the mth trace of the microseismic event waveform on the component is passed through a shifted time-discrete sequence of length tau, TmaxRepresenting the discrete time length of the single trace microseismic event waveform before the time domain cross-correlation operation.
4. The method of claim 1, wherein the instantaneous phase cross-correlation waveform is calculated as:
wherein the content of the first and second substances,
is the cross-correlation of the instantaneous phase of the nth trace micro-seismic event waveform with the instantaneous phase of the mth trace micro-seismic event waveform at vkTime-discrete sequences over components, vkTo representHorizontal component vxOr the vertical component vz，
The instantaneous phase of the waveform representing the nth microseismic event is at vkA time-discrete sequence of the components is,
a time discrete sequence representing the instantaneous phase of the mth trace of the microseismic event waveform is shifted by a length τ in vkTime discrete sequences on components, TmaxRepresents the discrete time length of the waveform of the single trace microseismic event before the instantaneous phase cross-correlation operation, and q is an index for controlling the curve sharpness degree of the instantaneous phase cross-correlation curve from a high coherence position to a low coherence position.
5. The method of claim 1, wherein the weighted cross-correlation waveform is calculated as:
wherein the content of the first and second substances,
is shown at vkThe nth and mth tracks on the component are time-discrete sequences of temporal cross-correlation weighted by instantaneous phase cross-correlation,
the waveform of the n-th micro-seismic event and the waveform of the m-th micro-seismic event are in vkA time-discrete sequence of time-domain cross-correlations over the components,
is the cross-correlation of the instantaneous phase of the nth trace micro-seismic event waveform with the instantaneous phase of the mth trace micro-seismic event waveform at vkTime discrete sequences over the components.
6. The method according to claim 1, wherein the travel time difference of the theoretical direct wave is a difference value between the travel time of the nth theoretical direct wave and the travel time of the mth theoretical direct wave in the N seismic waves; obtained by the following steps:
establishing a speed model for the area to be detected, and obtaining the theoretical direct wave travel time from each point in the area to be detected to each wave detector according to the speed model; the theoretical direct wave travel time comprises theoretical direct longitudinal wave travel time and theoretical direct transverse wave travel time.
7. The method according to claim 1, wherein in each component, offset stacking is performed on each corresponding weighted cross-correlation waveform by using a theoretical direct wave travel time difference, so as to obtain an interference imaging profile about each component, specifically:
in the vertical component, offset superposition is carried out on the corresponding cross-correlation waveform after weighting by using the travel time difference of the theoretical direct longitudinal wave, so as to obtain an interference imaging section of the longitudinal wave in the vertical component;
in the horizontal component, offset superposition is carried out on the corresponding cross-correlation waveform after weighting by using the travel time difference of the theoretical direct transverse wave, so as to obtain an interference imaging section of the transverse wave in the horizontal component.
8. The method of claim 1, wherein said subjecting seismic waves to components, when used for two-dimensional positioning, results in a vertical component and a horizontal component; when used for three-dimensional positioning, three components are obtained that are perpendicular to each other, the three components including: a vertical component, a first horizontal component, and a second horizontal component.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010207285.1A CN111352153B (en) | 2020-03-23 | 2020-03-23 | Microseism interference positioning method based on instantaneous phase cross-correlation weighting |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010207285.1A CN111352153B (en) | 2020-03-23 | 2020-03-23 | Microseism interference positioning method based on instantaneous phase cross-correlation weighting |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111352153A CN111352153A (en) | 2020-06-30 |
| CN111352153B true CN111352153B (en) | 2021-07-30 |
Family
ID=71196242
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010207285.1A Active CN111352153B (en) | 2020-03-23 | 2020-03-23 | Microseism interference positioning method based on instantaneous phase cross-correlation weighting |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111352153B (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113325467B (en) * | 2021-06-08 | 2023-10-24 | 中煤科工集团西安研究院有限公司 | Microseism focus positioning method based on slot wave frequency dispersion characteristics |
| CN113721296A (en) * | 2021-09-22 | 2021-11-30 | 应急管理部国家自然灾害防治研究院 | Method and device for processing remote seismic data |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103605151A (en) * | 2013-11-20 | 2014-02-26 | 中北大学 | Distributed group wave shallow-layer slight shock positioning method based on phase measuring |
| CN104765064A (en) * | 2015-03-25 | 2015-07-08 | 中国科学院声学研究所 | Microseism interference imaging method |
| CN105093317A (en) * | 2014-05-14 | 2015-11-25 | 中国石油化工股份有限公司 | Ground array type micro seismic data independent component separation denoising method |
| CN105403918A (en) * | 2015-12-09 | 2016-03-16 | 中国科学院地质与地球物理研究所 | Three-component microseism data effective event identification method and system |
| US10338246B1 (en) * | 2015-08-31 | 2019-07-02 | Seismic Innovations | Method and system for microseismic event wavefront estimation |
-
2020
- 2020-03-23 CN CN202010207285.1A patent/CN111352153B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103605151A (en) * | 2013-11-20 | 2014-02-26 | 中北大学 | Distributed group wave shallow-layer slight shock positioning method based on phase measuring |
| CN105093317A (en) * | 2014-05-14 | 2015-11-25 | 中国石油化工股份有限公司 | Ground array type micro seismic data independent component separation denoising method |
| CN104765064A (en) * | 2015-03-25 | 2015-07-08 | 中国科学院声学研究所 | Microseism interference imaging method |
| US10338246B1 (en) * | 2015-08-31 | 2019-07-02 | Seismic Innovations | Method and system for microseismic event wavefront estimation |
| CN105403918A (en) * | 2015-12-09 | 2016-03-16 | 中国科学院地质与地球物理研究所 | Three-component microseism data effective event identification method and system |
Non-Patent Citations (1)
| Title |
|---|
| "微地震地面监测资料反射波恢复研究";陈国金 等;《石油物探》;20200131;第59卷(第1期);第60-70页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111352153A (en) | 2020-06-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110133715B (en) | Microseism seismic source positioning method based on first-arrival time difference and waveform superposition | |
| CN103091710B (en) | Reverse time migration imaging method and device | |
| CN109669212B (en) | Seismic data processing method, stratum quality factor estimation method and device | |
| CN105182408A (en) | Manufacturing method and device for synthesizing earthquake record | |
| CN109856679B (en) | Method and system for imaging elastic wave Gaussian beam offset of anisotropic medium | |
| CN111352153B (en) | Microseism interference positioning method based on instantaneous phase cross-correlation weighting | |
| CN108089226B (en) | A kind of micro-seismic event automatic identifying method based on energy supposition between road | |
| Pan et al. | Estimating S-wave velocities from 3D 9-component shallow seismic data using local Rayleigh-wave dispersion curves–A field study | |
| US11754744B2 (en) | Surface wave prospecting method for jointly extracting Rayleigh wave frequency dispersion characteristics by seismoelectric field | |
| CN103926623A (en) | Method for suppressing reverse time migration low frequency noise | |
| CN110389377B (en) | Microseism offset imaging positioning method based on waveform cross-correlation coefficient multiplication | |
| Hayashi et al. | CMP spatial autocorrelation analysis of multichannel passive surface-wave data | |
| CN104570116A (en) | Geological marker bed-based time difference analyzing and correcting method | |
| CN102053260B (en) | Method for acquiring azimuth velocity of primary wave and method for processing earthquake data | |
| Luo et al. | Generation of a pseudo-2D shear-wave velocity section by inversion of a series of 1D dispersion curves | |
| CN112034520A (en) | Anisotropic medium dynamic focusing beam offset imaging method and system | |
| CN114415234A (en) | Method for determining shallow surface transverse wave velocity based on active source surface wave frequency dispersion and H/V | |
| Zhang et al. | Automated microseismic event location by amplitude stacking and semblance | |
| CN111999770B (en) | TTI medium conversion PS wave precise beam offset imaging method and system | |
| CN109143345B (en) | Quality factor Q nonlinear inversion method and system based on simulated annealing | |
| CN108693560B (en) | Scattered wave imaging method and system based on cross-correlation channel | |
| CN110780341B (en) | Anisotropic seismic imaging method | |
| CN110967751B (en) | Positioning method of micro-seismic event based on ground shallow well monitoring and storage medium | |
| CN106338760B (en) | The relief surface offset method of error compensation | |
| CN114114407B (en) | Surface wave and direct transverse wave suppression processing method for seismic wave detection |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |