CN110749925A - Broadband reverse time migration imaging processing method - Google Patents
Broadband reverse time migration imaging processing method Download PDFInfo
- Publication number
- CN110749925A CN110749925A CN201810819600.9A CN201810819600A CN110749925A CN 110749925 A CN110749925 A CN 110749925A CN 201810819600 A CN201810819600 A CN 201810819600A CN 110749925 A CN110749925 A CN 110749925A
- Authority
- CN
- China
- Prior art keywords
- wave field
- wave
- shot
- frequency
- decomposition
- 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.)
- Pending
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 141
- 230000005012 migration Effects 0.000 title claims abstract description 99
- 238000013508 migration Methods 0.000 title claims abstract description 99
- 238000003672 processing method Methods 0.000 title claims abstract description 18
- 238000001514 detection method Methods 0.000 claims abstract description 64
- 238000004458 analytical method Methods 0.000 claims abstract description 32
- 239000000523 sample Substances 0.000 claims abstract description 23
- 238000000354 decomposition reaction Methods 0.000 claims description 91
- 238000000034 method Methods 0.000 claims description 25
- 230000008569 process Effects 0.000 claims description 10
- 230000008859 change Effects 0.000 claims description 7
- 230000009286 beneficial effect Effects 0.000 abstract description 3
- 230000001737 promoting effect Effects 0.000 abstract 1
- 238000001914 filtration Methods 0.000 description 6
- 238000001228 spectrum Methods 0.000 description 6
- 230000005540 biological transmission Effects 0.000 description 4
- 238000005070 sampling Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 230000008030 elimination Effects 0.000 description 2
- 238000003379 elimination reaction Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000013213 extrapolation Methods 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000003860 storage Methods 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- 230000001133 acceleration Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000013501 data transformation Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000001629 suppression 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/38—Seismology; Seismic or acoustic prospecting or detecting specially adapted for water-covered areas
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Oceanography (AREA)
- Engineering & Computer Science (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- Remote Sensing (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Image Analysis (AREA)
Abstract
The invention relates to a broadband reverse time migration imaging processing method, which comprises the following steps: step 1: obtaining a shot point wave field and a wave detection point wave field at a certain moment; step 2: constructing a shot point analysis wave field and a wave detection point analysis wave field at the moment; and step 3: constructing a broadband reverse time migration imaging condition; and 4, step 4: applying the shot point analysis wave field and the demodulator probe analysis wave field to a broadband reverse time migration imaging condition to obtain a broadband reverse time migration imaging result at the moment; and 5: obtaining broadband reverse time migration imaging results of all propagation moments; the invention eliminates the influence of low-frequency migration false image in reverse time migration on imaging frequency band, protects the low-frequency information of effective information number, is beneficial to improving the imaging quality of high and steep fracture structure, improves the accuracy and resolution of imaging section and is beneficial to promoting high-accuracy earthquake imaging.
Description
Technical Field
The invention belongs to the technical field of seismic data prestack depth migration imaging, and particularly relates to a broadband reverse time migration imaging processing method.
Background
Compared with Kirchhoff's law (Kirchhoff) migration, the frequency band of the traditional reverse time migration imaging section is much narrower, and the requirement of broadband seismic data processing cannot be met. In particular, the marine data frequency band is generally wide, so when the conventional reverse time migration method is applied to the marine data, it is difficult to obtain high-resolution imaging results, and the imaging quality of a high-steep fracture structure is poor.
Therefore, a broad-band reverse time shift imaging technology needs to be researched, and especially, low-frequency information in an imaging section is protected from being damaged. The double-pass wave equation-based reverse time migration method performs zero-delay cross-correlation imaging on a shot point forward-transmission wave field and a demodulator probe reverse-transmission wave field along the time direction, and because a double-pass wave field propagation operator does not distinguish the wave field propagation directions, the wave fields with the same shot point and demodulator probe propagation directions are also imaged, which is the root cause of low-frequency migration noise generation on a reverse time migration propagation path.
Scholars at home and abroad carry out intensive research on the problem, for example, Baysal and the like adopt a reflection-free equation to reduce back reflection waves; loewenthal et al propose to smooth the model slowness by using a window function longer than the wavelength length to reduce back reflection; yoon et al use poynting vector imaging conditions to suppress offset noise; zhang et al propose a post-stack laplacian filtering method to suppress offset noise. The method has certain applicable conditions, wherein the Laplace filtering method after the superposition is simple to implement and high in calculation efficiency, and is widely applied to practical application. The traditional reverse time migration method for constructing the decomposition of the up wave field and the down wave field is an effective means for suppressing the low-frequency noise in the reverse time migration and protecting the low-frequency information in the data, but is generally realized in a frequency wave number domain or needs to solve a multiple fluctuation equation, which has great challenge for the practical application of the three-dimensional reverse time migration, because the transformation of mass data from a time space domain to the frequency wave number domain needs to continuously store a forward wave field and a backward wave field along a time direction and perform Fourier transformation related to the time direction, and the time dimension is in the slowest dimension of the data in the forward and backward extrapolation process of the reverse time migration wave field, so even if the acceleration is performed by utilizing CPU/GPU cooperative computing, the computing efficiency is still low, the storage requirement is very high, and the production requirement of the three-dimensional practical data cannot be met.
Disclosure of Invention
In order to solve the problems, the invention provides a broadband reverse time migration imaging processing method which can solve the problems of low resolution, narrow frequency band, low frequency information loss and poor imaging quality of high and steep section waves when marine seismic data is subjected to high-precision imaging.
The invention provides a broadband reverse time migration imaging processing method, which comprises the following steps:
step 1: obtaining a shot point wave field and a wave detection point wave field at a certain moment;
step 2: constructing a shot point analysis wave field and a wave detection point analysis wave field at the moment;
and step 3: constructing a broadband reverse time migration imaging condition;
and 4, step 4: applying the shot point analysis wave field and the demodulator probe analysis wave field to a broadband reverse time migration imaging condition to obtain a broadband reverse time migration imaging result at the moment;
and 5: and obtaining the broadband reverse time migration imaging results of all the propagation moments.
In one embodiment, in step 1, the shot wave field and the geophone wave field at a certain time are extrapolated by using a wave equation to obtain a shot wave field s (t, x) and a geophone wave field r (t, x) at the time; wherein, the wave equation expression is:
wherein x ═ x, y, z represents a spatial position; v (x) represents the offset velocity field at spatial position x; t represents time; p (t, x) represents the pressure field.
In one embodiment, the method specifically comprises the following steps:
step 2.1: performing Hilbert transform on the shot point wave field s (t, x) in a certain dimension direction to obtain sH(t, x); wherein s isH(t,x)=HT(s(t,x));
Step 2.2: subjecting the wave field r (t, x) of the wave detection point to Hilbert transform in a certain dimension direction to obtain rH(t,x),rH(t,x)=HT(r(t,x))。
In one embodiment, said step 2.2 is followed by a step 2.3: will sH(t, x) and rH(t, x) substituting into analytic wave field expression, and respectively constructing shot analytic wave fieldSum-point resolved wavefieldWherein, the shot point analysis wave field expression and the wave detection point analysis wave field expression are respectively as follows:
where i represents an imaginary symbol.
In an embodiment, the constructing the wideband reverse time shift imaging condition in step 3 further includes the following steps:
step 3.1: constructing a wave field decomposition related imaging condition of a frequency domain;
step 3.2: constructing a shot decomposition wave field of a frequency wave number domain and a wave detection point decomposition wave field of the frequency wave number domain;
step 3.3: and (3) applying the inverse Fourier change result of the shot decomposition wave field of the frequency wave number domain and the inverse Fourier change result of the wave detection point decomposition wave field of the frequency wave number domain in the step (3.2) to the wave field decomposition related imaging condition of the frequency domain in the step (3.1) to obtain the broadband reverse time migration imaging condition of time-space domain wave field implicit decomposition.
In an embodiment, the step 3.1 further specifically includes the following steps:
step 3.11: the zero-delay correlation imaging conditions of the shot point wave field and the wave detection point wave field are as follows:
wherein, TmaxExtending the maximum time for the wavefield;
step 3.12: respectively carrying out up-down decomposition on the shot point wave field and the wave detection point wave field as follows:
s(t,x)=su(t,x)+sd(t,x)
r(t,x)=ru(t,x)+rd(t,x)
wherein s isd(t, x) is the down-going wavefield at the shot, su(t, x) is the shot upgoing wavefield; r isd(t, x) is the wave field descending from the detection point, ru(t, x) is a demodulator probe upgoing wave field;
step 3.13: respectively bringing the up-and-down wave field of the shot point wave field and the up-and-down wave field of the wave detection point wave field into a zero-delay cross-correlation imaging condition to obtain a wave field decomposition correlation imaging condition as follows:
wherein, I1(x) The correlation imaging result of the shot point downgoing wave field and the wave detection point upgoing wave field is obtained; i is2(x) For correlation of the up-going wave field of the shot with the down-going wave field of the demodulator probe, I3(x) The correlation imaging result of the shot point descending wave field and the wave detection point descending wave field is obtained; i is4(x) The correlation imaging result of the shot point upgoing wave field and the demodulator probe upgoing wave field is obtained;
step 3.14: in the reverse time migration imaging process, selecting an imaging result contributing to imaging information as a final imaging result, and then eliminating the wave field decomposition related imaging conditions after low-frequency migration noise:
step 3.15: the wave field decomposition related imaging condition after eliminating the low-frequency offset noise is popularized to a frequency domain, and the expression is as follows:
wherein,
up-going wavefield s for shotuTransposing the result of the one-dimensional Fourier transform of (t, x),for the down-going wavefield of the shot sdTransposing the one-dimensional Fourier transform result of (t, x);
for the detection point up-going wave field ruTransposing the result of the one-dimensional Fourier transform of (t, x),for the detector point down-going wave field rdTransposing the one-dimensional Fourier transform result of (t, x); "+" indicates the complex conjugate.
In one embodiment, the step 3.2 further specifically includes the following steps:
step 3.21: in the frequency wave number domain, respectively carrying out up-going and down-going wave field decomposition on a shot point wave field and a wave detection point wave field;
step 3.22: defining generalized uplink and downlink shot point decomposed wave fields and generalized uplink and downlink demodulator probe decomposed wave fields;
step 3.23: and (5) combining the step 3.21 and the step 3.22 to obtain a shot decomposition wave field of a frequency wave number domain and a demodulator probe decomposition wave field of the frequency wave number domain.
In one embodiment, in step 3.21, the shot wavefield up-and-down wavefield decomposition expression is:
wherein w is frequency; k is a radical ofzIs wave number, wkzIs the product of frequency and wave number; sd(w,kz) For the down-going wave field, s, of the frequency-wavenumber-domain shot-wave fieldu(w,kz) Up-going wavefields, s (w, k), for the frequency-wavenumber-domain shot wavefieldz) Fourier forward transform of shot wavefields in the time direction and the depth direction;
the decomposition expression of the up-going wave field and the down-going wave field of the wave detection point wave field is as follows:
wherein r isd(w,kz) For the frequency-wavenumber domain, the detection point, the wavefield, ru(w,kz) For the frequency-wavenumber domain, the detection-point wavefield, the upgoing wavefield, r (w, k)z) The Fourier forward transform of the wave field at the wave detection point along the time direction and the depth direction is carried out.
In step 3.22, the generalized uplink and downlink shot decomposition wave field expression is as follows:
wherein s is+(w,kz) And s-(w,kz) The wave fields are decomposed by generalized uplink and downlink shot points;
the generalized up-down detection point decomposition wave field expression is as follows:
wherein r is+(w,kz) And r-(w,kz) The wave fields are decomposed by generalized uplink and downlink detection points.
In one embodiment, in step 3.23, the up-and-down wavefield decomposition expression of the shot point wavefield in step 3.21 is combined with the generalized up-and-down shot point decomposition wavefield expression in step 3.22 to obtain a shot point decomposition wavefield in the frequency wave number domain, where the expression is:
similarly, the analysis point wave field up-down going wave field decomposition expression in step 3.1 and the generalized analysis point wave field up-down going wave field decomposition expression in step 3.2 are combined to obtain the analysis point decomposition wave field in the frequency wave number domain, and the expression is:
from the shot decomposed wavefield in the frequency-wavenumber domain and the demodulator decomposed wavefield in the frequency-wavenumber domain, the direction of the generalized decomposed wavefield is implicitly present, which depends on the sign of the current frequency w.
In one embodiment, the inverse fourier transform result of the shot decomposition wavefield in the frequency wavenumber domain and the inverse fourier transform result of the demodulator decomposition wavefield in the frequency wavenumber domain in step 3.23 are applied to the wavefield decomposition related imaging conditions in the frequency domain in step 3.15 to obtain a broadband reverse time migration imaging condition:
wherein Re represents a real part.
The invention provides a broadband reverse time migration imaging condition aiming at the problem of narrow reverse time migration imaging frequency band, realizes implicit decomposition of a wave field efficiently by a fast algorithm in the time-space domain imaging process, avoids two wave field extrapolation processes required by frequency wave number domain data transformation and wave field display decomposition, eliminates the influence of a low frequency migration false image in reverse time migration on an imaging frequency band, protects the low frequency information of an effective information number, is beneficial to improving the imaging quality of a high and steep fracture structure, and improves the imaging section precision and resolution.
Drawings
The invention will be described in more detail hereinafter on the basis of embodiments and with reference to the accompanying drawings. Wherein:
FIG. 1 is a flow chart of a broadband reverse time migration imaging processing method according to the present invention;
FIG. 2 is a model of the velocity of a multi-layer horizontal reflective interface of the present invention;
FIGS. 3 a-3 c are conventional reverse time migration wavefield snapshots and single shot migration results; wherein, fig. 3a is a shot wavefield when t is 1.7s, and fig. 3b is a geophone wavefield when t is 1.7 s; FIG. 3c is a single shot migration result;
FIGS. 4a to 4c are wave field snapshots and single shot migration results after decomposition of the broadband reverse time imaging condition; wherein, fig. 4a shows the down-going wave (left) and the up-going wave (right) of the shot wave field when t is 1.7s, and fig. 4b shows the down-going wave (left) and the up-going wave (right) of the inspection wave field when t is 1.7 s; FIG. 4c is a single shot migration result;
FIGS. 5 a-5 c are comparisons of real data line-inversion offset results; wherein, FIG. 5a shows the conventional reverse time migration result; FIG. 5b shows the reverse time migration result of the present invention; FIG. 5c is the Laplace filtering result after imaging of FIG. 5 a;
FIGS. 6a to 6c are two reverse time migration result comparisons of actual data measurement lines; wherein, FIG. 6a shows the conventional reverse time migration result; FIG. 6b shows the reverse time migration result of the present invention; fig. 6c shows the laplacian filtering results after imaging in fig. 6 a.
In the drawings like parts are provided with the same reference numerals. The figures are not drawn to scale.
Detailed Description
The invention will be further explained with reference to the drawings. Therefore, the realization process of how to apply the technical means to solve the technical problems and achieve the technical effect can be fully understood and implemented. It should be noted that the technical features mentioned in the embodiments can be combined in any way as long as no conflict exists. It is intended that the invention not be limited to the particular embodiments disclosed, but that the invention will include all embodiments falling within the scope of the appended claims.
The invention provides a broadband reverse time migration imaging processing method, which comprises the following steps:
step 1: obtaining a shot point wave field and a wave detection point wave field at a certain moment;
step 2: constructing a shot point analysis wave field and a wave detection point analysis wave field at the moment;
and step 3: constructing a broadband reverse time migration imaging condition;
and 4, step 4: applying the shot point analysis wave field and the demodulator probe analysis wave field to a broadband reverse time migration imaging condition to obtain a broadband reverse time migration imaging result at the moment;
and 5: and obtaining the broadband reverse time migration imaging results of all the propagation moments.
Preferably, in step 1, in a constant-density acoustic medium, the shot wave field and the geophone wave field at a certain time are extrapolated by a wave equation to obtain a shot wave field s (t, x) and a geophone wave field r (t, x) at the time; wherein, the wave equation expression is:
wherein x ═ x, y, z represents a spatial position; v (x) represents the offset velocity field at spatial position x; t represents time; p (t, x) represents the pressure field.
Preferably, the method specifically comprises the following steps:
step 2.1: performing Hilbert transform on the shot point wave field s (t, x) in a certain dimension direction to obtain sH(t, x); wherein s isH(t,x)=HT(s(t,x));
Step 2.2: subjecting the wave field r (t, x) of the wave detection point to Hilbert transform in a certain dimension direction to obtain rH(t,x),rH(t,x)=HT(r(t,x))。
Preferably, step 2.2 is followed by step 2.3: will sH(t, x) and rH(t, x) substituting into analytic wave field expression, and respectively constructing shot analytic wave fieldSum-point resolved wavefieldWherein, the shot point analysis wave field expression and the wave detection point analysis wave field expression are respectively as follows:
where i represents an imaginary symbol.
Preferably, the constructing of the broadband reverse time shift imaging condition in step 3 further specifically includes the following steps:
step 3.1: constructing a wave field decomposition related imaging condition of a frequency domain;
step 3.2: constructing a shot decomposition wave field of a frequency wave number domain and a wave detection point decomposition wave field of the frequency wave number domain;
step 3.3: and (3) applying the inverse Fourier change result of the shot decomposition wave field of the frequency wave number domain and the inverse Fourier change result of the wave detection point decomposition wave field of the frequency wave number domain in the step (3.2) to the wave field decomposition related imaging condition of the frequency domain in the step (3.1) to obtain the broadband reverse time migration imaging condition of time-space domain wave field implicit decomposition.
Preferably, in step 3.1, the method further specifically comprises the following steps:
step 3.11: and extracting a reverse time migration imaging result of the underground interface by using the shot wave field s (t, x) and the wave detection point wave field r (t, x) under a zero delay correlation imaging condition, wherein the zero delay correlation imaging condition of the shot wave field and the wave detection point wave field is as follows:
wherein, TmaxExtending the maximum time for the wavefield;
the zero-delay correlation imaging condition can generate serious strong amplitude and low-frequency offset noise when the propagation directions in the shot point wave field and the wave detection point wave field are the same, so that the reverse time offset imaging quality is reduced.
Therefore, the suppression of low-frequency offset noise can be realized by constructing the reverse time offset imaging condition of wave field propagation direction decomposition, and simultaneously, low-frequency information is protected;
step 3.12: respectively carrying out up-down decomposition on the shot point wave field and the wave detection point wave field as follows:
s(t,x)=su(t,x)+sd(t,x)
r(t,x)=ru(t,x)+rd(t,x)
wherein s isd(t, x) is the down-going wavefield at the shot, su(t, x) is the shot pointAn up-going wavefield; r isd(t, x) is the wave field descending from the detection point, ru(t, x) is a demodulator probe upgoing wave field;
step 3.13: respectively bringing the up-and-down wave field of the shot point wave field and the up-and-down wave field of the wave detection point wave field into a zero-delay cross-correlation imaging condition to obtain a wave field decomposition correlation imaging condition as follows:
wherein, I1(x) The correlation imaging result of the shot point downgoing wave field and the wave detection point upgoing wave field is obtained; i is2(x) For correlation of the up-going wave field of the shot with the down-going wave field of the demodulator probe, I3(x) The correlation imaging result of the shot point descending wave field and the wave detection point descending wave field is obtained; i is4(x) The correlation imaging result of the shot point upgoing wave field and the demodulator probe upgoing wave field is obtained;
I1(x) And I2(x) Reflecting the imaging information of the underground interface for the wave field imaging result with the opposite propagation direction; i is3(x) And I4(x) The wave field imaging result with the same propagation direction does not contribute to the underground structure imaging result and is the root cause of the low-frequency offset noise.
Step 3.14: in the reverse time migration imaging process, selecting an imaging result contributing to imaging information as a final imaging result, and then eliminating the wave field decomposition related imaging conditions after low-frequency migration noise:
step 3.15: the wave field decomposition related imaging condition after eliminating the low-frequency offset noise is popularized to a frequency domain, and the expression is as follows:
wherein,
up-going wavefield s for shotuTransposing the result of the one-dimensional Fourier transform of (t, x),for the down-going wavefield of the shot sdTransposing the one-dimensional Fourier transform result of (t, x);
for the detection point up-going wave field ruTransposing the result of the one-dimensional Fourier transform of (t, x),for the detector point down-going wave field rdTransposing the one-dimensional Fourier transform result of (t, x); "+" indicates the complex conjugate.
The wave field decomposition related imaging conditions show that the key point for obtaining an ideal reverse time migration imaging result is to obtain wave fields with opposite propagation directions in a shot point wave field and a wave detection point wave field and perform cross correlation along the propagation time.
Therefore, the broadband reverse time migration cross-correlation imaging condition is provided, implicit fast decomposition of a wave field is achieved, interference of low-frequency migration noise on a reverse time migration imaging result is eliminated, low-frequency information in effective signals is protected, and a reverse time migration imaging result with a relatively broadband is obtained.
Preferably, said step 3.2, in turn, specifically comprises the steps of:
step 3.21: in the frequency wave number domain, respectively carrying out up-going and down-going wave field decomposition on a shot point wave field and a wave detection point wave field;
step 3.22: defining generalized uplink and downlink shot point decomposed wave fields and generalized uplink and downlink demodulator probe decomposed wave fields;
step 3.23: and (5) combining the step 3.21 and the step 3.22 to obtain a shot decomposition wave field of a frequency wave number domain and a demodulator probe decomposition wave field of the frequency wave number domain.
Preferably, in step 3.21, the shot wavefield up-and-down wavefield decomposition expression is:
wherein w is frequency; k is a radical ofzIs wave number, wkzIs the product of frequency and wave number; sd(w,kz) For the down-going wave field, s, of the frequency-wavenumber-domain shot-wave fieldu(w,kz) Up-going wavefields, s (w, k), for the frequency-wavenumber-domain shot wavefieldz) Fourier forward transform of shot wavefields in the time direction and the depth direction;
the decomposition expression of the up-going wave field and the down-going wave field of the wave detection point wave field is as follows:
wherein r isd(w,kz) For the frequency-wavenumber domain, the detection point, the wavefield, ru(w,kz) For the frequency-wavenumber domain, the detection-point wavefield, the upgoing wavefield, r (w, k)z) The Fourier forward transform of the wave field at the wave detection point along the time direction and the depth direction is carried out.
In step 3.22, in order to implement efficient implicit decomposition of the shot point wave field and the wave detection point wave field, a generalized uplink and downlink shot point decomposition wave field expression is defined as follows:
wherein s is+(w,kz) And s-(w,kz) Are all made ofDecomposing a wave field for generalized uplink and downlink shot points;
defining a generalized uplink and downlink demodulator probe decomposition wave field expression as follows:
wherein r is+(w,kz) And r-(w,kz) The wave fields are decomposed by generalized uplink and downlink detection points.
Preferably, in step 3.23, the up-and-down wavefield decomposition expression of the shot point wavefield in step 3.21 and the generalized up-and-down shot point decomposition wavefield expression in step 3.22 are combined to obtain a shot point decomposition wavefield in the frequency wave number domain, where the expression is:
similarly, the analysis point wave field up-down going wave field decomposition expression in step 3.1 and the generalized analysis point wave field up-down going wave field decomposition expression in step 3.2 are combined to obtain the analysis point decomposition wave field in the frequency wave number domain, and the expression is:
from the shot decomposed wavefield in the frequency-wavenumber domain and the demodulator decomposed wavefield in the frequency-wavenumber domain, the direction of the generalized decomposed wavefield is implicitly present, which depends on the sign of the current frequency w.
Preferably, the result s of the inverse Fourier transform of the shot decomposition wavefield in the frequency wavenumber domain of step 3.23 is usedd(t,x),suResults r of inverse Fourier transform of the (t, x) and demodulator-decomposed wavefield in the frequency wavenumber domaind(t,x),ru(t, x) is applied to the wave field decomposition correlation imaging condition of the frequency domain after the low frequency offset noise elimination in step 3.15, at this time, no matter whether the sign of the current frequency is positive or negative, the broadband reverse time offset imaging condition of the constant time-space domain wave field implicit decomposition holds:
wherein Re represents a real part.
Based on the above analysis, to obtain a generalized implicit decomposed wavefield s+(t, z) it is critical to ensure s+The wave number component of (t, z) is constantly greater than zero. In the digital signal, the analytic signal has a property that the spectrum (or wavenumber spectrum) contains only positive-valued components, and the positive-valued spectrum (or wavenumber spectrum) thereof has a 2-fold relationship with the spectrum (or wavenumber spectrum) of the original signal.
For a comparative description of the imaging process of the conventional reverse time migration method and the imaging process of the method, a multilayer horizontal reflection interface velocity model shown in fig. 2 is constructed, the velocity is 2000m/s, 2500m/s and 3000m/s respectively, the model size is 500 longitudinal sampling points, 400 transverse sampling points, the sampling interval is 10m, and fig. 3 and 4 are snapshot and single-shot migration results after wave field decomposition obtained by forward transmission at the shot point end and backward transmission at the detection end when the propagation time is 1.7s under the conventional reverse time migration and application of broadband imaging conditions, respectively. FIG. 3a and FIG. 3b show that the forward and backward wave fields include both up and down going waves, and the single shot migration result is shown in FIG. 3c, which has strong low frequency migration noise interference; fig. 4 is a shot point end forward transmission wave field snapshot and a detection end backward transmission wave field snapshot at the same time, an up-going wave field and a down-going wave field are separated, and a single shot migration result is compared, as can be seen from fig. 3c and fig. 4c, the new imaging condition implicitly only carries out selective imaging on the wave fields with opposite propagation directions, the elimination of low-frequency migration noise is realized in the migration imaging process, and the imaging quality of an effective wave field is ensured.
The technology of the invention is applied to actual data, and application processing of different measuring lines is completed. Fig. 5 and fig. 6 are comparison results of reverse time migration of different lines in the same probe region. Fig. 5a and fig. 6a are both conventional reverse time migration results, fig. 5b and fig. 6b are reverse time migration results obtained by the present invention, and fig. 5c and fig. 6c are laplacian filtering results after imaging of fig. 5a and fig. 6a, respectively. The conventional reverse time migration profile has low wave number migration noise interference, which affects the imaging quality of the imaging profile, especially the shallow and middle layers, and affects the frequency band distribution of the profile, so that the resolution of the profile is reduced. The method suppresses the low-frequency offset artifact well, protects the low-frequency information in the data, and improves the imaging resolution (compare fig. 5a and 5b or fig. 6a and 6 b); due to the post-stack laplacian filtering method, the whole data frequency band moves to a high wave number, low-frequency information in the data is lost, the continuity of the in-phase axis is deteriorated, and the signal-to-noise ratio is also reduced to a certain extent (compare fig. 5b and 5c or fig. 6b and 6 c). In addition, the conventional reverse time migration results are about 1210 at cdp and the depth ranges from 3.5km to 4km, as shown in fig. 6a, so that a sharp fracture structure is obvious, the reverse time migration results obtained by the method are well maintained as shown in fig. 6b, and the sharp fracture structure is blurred due to the frequency band change and the loss of low-frequency effective signals as shown in fig. 6 c.
The application results of the theoretical model and the actual data show that: the broadband reverse time migration imaging condition provided by the invention realizes the rapid implicit decomposition of a wave field in a time space domain, eliminates the influence of a low-frequency migration false image in reverse time migration on an imaging frequency band, protects the low-frequency information of an effective information number, and improves the accuracy and the resolution of an imaging section.
Compared with the prior art, the invention has the advantages that: compared with the conventional reverse time migration algorithm, the method does not need extra storage space, only needs to pay the calculation cost of constructing and analyzing the wave field, increases the calculated amount to be less than 10%, has important practical significance for the three-dimensional reverse time migration practical production application facing mass data, reserves low-frequency information in the obtained imaging result frequency band, and can improve the imaging quality of high and steep structures in the broadband data.
While the present invention has been described with reference to the preferred embodiments as above, the description is only for the convenience of understanding the present invention and is not intended to limit the present invention. It will be understood by those skilled in the art that various changes in form and details may be made therein without departing from the spirit and scope of the invention as defined by the appended claims.
Claims (10)
1. A broadband reverse time migration imaging processing method is characterized by comprising the following steps:
step 1: obtaining a shot point wave field and a wave detection point wave field at a certain moment;
step 2: constructing a shot point analysis wave field and a wave detection point analysis wave field at the moment;
and step 3: constructing a broadband reverse time migration imaging condition;
and 4, step 4: applying the shot point analysis wave field and the demodulator probe analysis wave field to a broadband reverse time migration imaging condition to obtain a broadband reverse time migration imaging result at the moment;
and 5: and obtaining the broadband reverse time migration imaging results of all the propagation moments.
2. The broadband reverse time migration imaging processing method according to claim 1, wherein in step 1, the shot wave field and the geophone wave field at a certain time are extrapolated by a wave equation to obtain the shot wave field s (t, x) and the geophone wave field r (t, x) at the time; wherein, the wave equation expression is:
wherein x ═ x, y, z represents a spatial position; v (x) represents the offset velocity field at spatial position x; t represents time; p (t, x) represents the pressure field.
3. The broadband reverse time migration imaging processing method according to claim 2, wherein the step 2 specifically includes the following steps:
step 2.1: performing Hilbert transform on the shot point wave field s (t, x) in a certain dimension direction to obtain sH(t, x); wherein s isH(t,x)=HT(s(t,x));
Step 2.2: subjecting the wave field r (t, x) of the wave detection point to Hilbert transform in a certain dimension direction to obtain rH(t,x),rH(t,x)=HT(r(t,x))。
4. The broadband reverse time migration imaging processing method according to claim 3, wherein said step 2.2 is further followed by a step 2.3: will sH(t, x) and rH(t, x) substituting into analytic wave field expression, and respectively constructing shot analytic wave fieldSum-point resolved wavefieldWherein, the shot point analysis wave field expression and the wave detection point analysis wave field expression are respectively as follows:
where i represents an imaginary symbol.
5. The broadband reverse time migration imaging processing method according to claim 4, wherein the constructing of the broadband reverse time migration imaging condition in step 3 further comprises the following steps:
step 3.1: constructing a wave field decomposition related imaging condition of a frequency domain;
step 3.2: constructing a shot decomposition wave field of a frequency wave number domain and a wave detection point decomposition wave field of the frequency wave number domain;
step 3.3: and (3) applying the inverse Fourier change result of the shot decomposition wave field of the frequency wave number domain and the inverse Fourier change result of the wave detection point decomposition wave field of the frequency wave number domain in the step (3.2) to the wave field decomposition related imaging condition of the frequency domain in the step (3.1) to obtain the broadband reverse time migration imaging condition of time-space domain wave field implicit decomposition.
6. The broadband reverse time migration imaging processing method according to claim 5, wherein the step 3.1 further comprises the following steps:
step 3.11: the zero-delay correlation imaging conditions of the shot point wave field and the wave detection point wave field are as follows:
wherein, TmaxExtending the maximum time for the wavefield;
step 3.12: respectively carrying out up-down decomposition on the shot point wave field and the wave detection point wave field as follows:
s(t,x)=su(t,x)+sd(t,x)
r(t,x)=ru(t,x)+rd(t,x)
wherein s isd(t, x) is the down-going wavefield at the shot, su(t, x) is the shot upgoing wavefield; r isd(t, x) is the wave field descending from the detection point, ru(t, x) is a demodulator probe upgoing wave field;
step 3.13: respectively bringing the up-and-down wave field of the shot point wave field and the up-and-down wave field of the wave detection point wave field into a zero-delay cross-correlation imaging condition to obtain a wave field decomposition correlation imaging condition as follows:
wherein, I1(x) The correlation imaging result of the shot point downgoing wave field and the wave detection point upgoing wave field is obtained;I2(x) For correlation of the up-going wave field of the shot with the down-going wave field of the demodulator probe, I3(x) The correlation imaging result of the shot point descending wave field and the wave detection point descending wave field is obtained; i is4(x) The correlation imaging result of the shot point upgoing wave field and the demodulator probe upgoing wave field is obtained;
step 3.14: in the reverse time migration imaging process, selecting an imaging result contributing to imaging information as a final imaging result, and then eliminating the wave field decomposition related imaging conditions after low-frequency migration noise:
step 3.15: the wave field decomposition related imaging condition after eliminating the low-frequency offset noise is popularized to a frequency domain, and the expression is as follows:
wherein,
up-going wavefield s for shotuTransposing the result of the one-dimensional Fourier transform of (t, x),for the down-going wavefield of the shot sdTransposing the one-dimensional Fourier transform result of (t, x);
7. The broadband reverse time migration imaging processing method according to claim 6, wherein said step 3.2 further comprises the following steps:
step 3.21: in the frequency wave number domain, respectively carrying out up-going and down-going wave field decomposition on a shot point wave field and a wave detection point wave field;
step 3.22: defining generalized uplink and downlink shot point decomposed wave fields and generalized uplink and downlink demodulator probe decomposed wave fields;
step 3.23: and (5) combining the step 3.21 and the step 3.22 to obtain a shot decomposition wave field of a frequency wave number domain and a demodulator probe decomposition wave field of the frequency wave number domain.
8. The broadband reverse time migration imaging processing method according to claim 7, wherein in step 3.21, the down-going and up-going wavefield decomposition expression of the shot wavefield is:
wherein w is frequency; k is a radical ofzIs wave number, wkzIs the product of frequency and wave number; sd(w,kz) For the down-going wave field, s, of the frequency-wavenumber-domain shot-wave fieldu(w,kz) Up-going wavefields, s (w, k), for the frequency-wavenumber-domain shot wavefieldz) Fourier forward transform of shot wavefields in the time direction and the depth direction;
the decomposition expression of the up-going wave field and the down-going wave field of the wave detection point wave field is as follows:
wherein r isd(w,kz) For the frequency-wavenumber domain, the detection point, the wavefield, ru(w,kz) For the frequency-wavenumber domain, the detection-point wavefield, the upgoing wavefield, r (w, k)z) Fourier forward transform of a wave field at a wave detection point along a time direction and a depth direction;
in step 3.22, the generalized uplink and downlink shot decomposition wave field expression is as follows:
wherein s is+(w,kz) And s- (w, k)z) The wave fields are decomposed by generalized uplink and downlink shot points;
the generalized up-down detection point decomposition wave field expression is as follows:
wherein r is+(w,kz) And r- (w, k)z) The wave fields are decomposed by generalized uplink and downlink detection points.
9. The processing method of broadband reverse time migration imaging according to claim 8, wherein in step 3.23, the up-and-down going shot point wavefield decomposition expression in step 3.21 is combined with the generalized up-and-down going shot point decomposition wavefield expression in step 3.22 to obtain the shot point decomposition wavefield in the frequency wavenumber domain, and the expression is:
similarly, the analysis point wave field up-down going wave field decomposition expression in step 3.1 and the generalized analysis point wave field up-down going wave field decomposition expression in step 3.2 are combined to obtain the analysis point decomposition wave field in the frequency wave number domain, and the expression is:
from the shot decomposed wavefield in the frequency-wavenumber domain and the demodulator decomposed wavefield in the frequency-wavenumber domain, the direction of the generalized decomposed wavefield is implicitly present, which depends on the sign of the current frequency w.
10. The broadband reverse time migration imaging processing method according to claim 9, wherein the inverse fourier transform result of the shot decomposition wave field in the frequency wavenumber domain and the inverse fourier transform result of the demodulator decomposition wave field in the frequency wavenumber domain in step 3.23 are applied to the wave field decomposition related imaging conditions in the frequency domain in step 3.15 to obtain the broadband reverse time migration imaging conditions:
wherein Re represents a real part.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810819600.9A CN110749925A (en) | 2018-07-24 | 2018-07-24 | Broadband reverse time migration imaging processing method |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810819600.9A CN110749925A (en) | 2018-07-24 | 2018-07-24 | Broadband reverse time migration imaging processing method |
Publications (1)
Publication Number | Publication Date |
---|---|
CN110749925A true CN110749925A (en) | 2020-02-04 |
Family
ID=69275485
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810819600.9A Pending CN110749925A (en) | 2018-07-24 | 2018-07-24 | Broadband reverse time migration imaging processing method |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110749925A (en) |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112904368A (en) * | 2021-01-25 | 2021-06-04 | 中国科学院西安光学精密机械研究所 | Non-visual field three-dimensional reconstruction method and system based on analytic signal and compensation reference function |
CN114285393A (en) * | 2021-12-22 | 2022-04-05 | 成都理工大学 | Method and device for eliminating reverse time migration noise |
CN116299691A (en) * | 2022-11-25 | 2023-06-23 | 西南石油大学 | F-K domain-based passive source surface wave imaging method and data screening method |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20120275268A1 (en) * | 2011-04-28 | 2012-11-01 | Cggveritas Services Sa | Device and method for extrapolating specular energy of reverse time migration three dimensional angle gathers |
CN104570114A (en) * | 2013-10-12 | 2015-04-29 | 中国石油化工股份有限公司 | Reverse time migration noise suppression method based on wave field decomposition |
US20170192118A1 (en) * | 2016-01-05 | 2017-07-06 | Schlumerger Technology Corporation | Amplitude Inversion on Partitioned Depth Image Gathers Using Point Spread Functions |
CN107153216A (en) * | 2017-07-05 | 2017-09-12 | 中国科学院地质与地球物理研究所 | Determine method, device and the computer-readable storage medium of the Poynting vector of seismic wave field |
CN108072895A (en) * | 2016-11-09 | 2018-05-25 | 中国石油化工股份有限公司 | A kind of anisotropy pre-Stack Reverse optimization method based on GPU |
-
2018
- 2018-07-24 CN CN201810819600.9A patent/CN110749925A/en active Pending
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20120275268A1 (en) * | 2011-04-28 | 2012-11-01 | Cggveritas Services Sa | Device and method for extrapolating specular energy of reverse time migration three dimensional angle gathers |
CN104570114A (en) * | 2013-10-12 | 2015-04-29 | 中国石油化工股份有限公司 | Reverse time migration noise suppression method based on wave field decomposition |
US20170192118A1 (en) * | 2016-01-05 | 2017-07-06 | Schlumerger Technology Corporation | Amplitude Inversion on Partitioned Depth Image Gathers Using Point Spread Functions |
CN108072895A (en) * | 2016-11-09 | 2018-05-25 | 中国石油化工股份有限公司 | A kind of anisotropy pre-Stack Reverse optimization method based on GPU |
CN107153216A (en) * | 2017-07-05 | 2017-09-12 | 中国科学院地质与地球物理研究所 | Determine method, device and the computer-readable storage medium of the Poynting vector of seismic wave field |
Non-Patent Citations (4)
Title |
---|
FAQI LIU ET AL.: "An effective imaging condition for reverse-time migration using wavefield decomposition", 《GEOPHYSICS》 * |
TONG W. FEI ET AL.: "Removing false images in reverse time migration:The concept of de-primary", 《GEOPHYSICS》 * |
徐兴荣等: "基于波场分离理论的逆时偏移成像条件研究及应用", 《地球物理学进展》 * |
王一博等: "基于 Hilbert变换的全波场分离逆时偏移成像", 《地球物理学报》 * |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112904368A (en) * | 2021-01-25 | 2021-06-04 | 中国科学院西安光学精密机械研究所 | Non-visual field three-dimensional reconstruction method and system based on analytic signal and compensation reference function |
CN112904368B (en) * | 2021-01-25 | 2023-09-29 | 中国科学院西安光学精密机械研究所 | Non-visual field three-dimensional reconstruction method and system based on analytic signal and compensation reference function |
CN114285393A (en) * | 2021-12-22 | 2022-04-05 | 成都理工大学 | Method and device for eliminating reverse time migration noise |
CN114285393B (en) * | 2021-12-22 | 2022-08-02 | 成都理工大学 | Method and device for eliminating reverse time migration noise |
CN116299691A (en) * | 2022-11-25 | 2023-06-23 | 西南石油大学 | F-K domain-based passive source surface wave imaging method and data screening method |
CN116299691B (en) * | 2022-11-25 | 2023-10-27 | 西南石油大学 | F-K domain-based passive source surface wave imaging method and data screening method |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN106970419B (en) | A kind of non-homogeneous bent wave 3D seismic data method for reconstructing based on linear Bregman algorithm | |
CN106249291B (en) | A kind of high precision seismic data re-establishing method based on the non-homogeneous warp wavelet of two dimension | |
Marfurt et al. | 3-D broad-band estimates of reflector dip and amplitude | |
US9625593B2 (en) | Seismic data processing | |
CN101614826B (en) | Method and device for realizing binning homogenization in three-dimensional seismic data processing | |
CN106597539B (en) | For the bent wave zone Radon converter noise drawing method of Huangtuyuan area | |
US20120275267A1 (en) | Seismic Data Processing | |
CN102156296A (en) | Elastic reverse time migration imaging method by combining seismic multi-component | |
US9250341B2 (en) | Device and method for extrapolating specular energy of reverse time migration three dimensional angle gathers | |
CN111505718B (en) | High-resolution underground structure amplitude-preserving imaging method | |
CN110749925A (en) | Broadband reverse time migration imaging processing method | |
CN110456417B (en) | Seismic data multiple suppression method | |
US20200393582A1 (en) | Method and system for generating geophysical data | |
CN104459770B (en) | A kind of method for regularizing high-dimensional seismic data | |
Biondi et al. | Transformation of 3-D prestack data by azimuth moveout (AMO) | |
CN103558636B (en) | A kind of method that gathers footprint decay from post-stack seismic data | |
CN108562937A (en) | A kind of seismic imaging method | |
CN111257938A (en) | Time-lapse seismic virtual source wave field reconstruction method and system based on wavelet cross-correlation | |
CN113866821B (en) | Passive source interference offset imaging method and system based on illumination direction constraint | |
Wu et al. | High resolution beam forming for 3D common offset Kirchhoff beam migration | |
CN111562616B (en) | Method and device for suppressing scattered noise of seismic data | |
Wang et al. | Ground roll wave suppression based on wavelet frequency division and radial trace transform | |
CN102998702A (en) | Amplitude-preserving plane wave prestack depth migration method | |
CN107589456B (en) | Method and device for acquiring seismic data and computer readable storage medium | |
CN112213775A (en) | Fidelity frequency-boosting method for high-coverage-frequency pre-stack seismic data |
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 | ||
RJ01 | Rejection of invention patent application after publication |
Application publication date: 20200204 |
|
RJ01 | Rejection of invention patent application after publication |