CN107632321A - A kind of offset imaging method - Google Patents
A kind of offset imaging method Download PDFInfo
- Publication number
- CN107632321A CN107632321A CN201710780533.XA CN201710780533A CN107632321A CN 107632321 A CN107632321 A CN 107632321A CN 201710780533 A CN201710780533 A CN 201710780533A CN 107632321 A CN107632321 A CN 107632321A
- Authority
- CN
- China
- Prior art keywords
- dreamlet
- local
- detector
- domain
- propagation
- 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.)
- Granted
Links
Landscapes
- Medicines Containing Antibodies Or Antigens For Use As Internal Diagnostic Agents (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
A kind of offset imaging method, by geological data u (t, xs,xr) dreamlet Orthogonal Decompositions are carried out, obtain the dreamlet coefficients of wave field；In each depth, first using the common small echo electron gun of dreamlet operators continuation, the wave field after continuation is settled into minor beam wave detector altogether using dreamlet operators again；Wave field transformation to frequency domain is subjected to local dip correction, then using image-forming condition, extraction zero-offset zero moment wave field amplitude is as picture value.Dreamlet propagation operators utilize temporal localization in the present invention, in continuation wave field, if dreamlet coefficients get over zero-acrross ing moment after propagation and reach the negative time, disregarded by casting out, therefore hash can be automatically removed during wave field extrapolation, without doing extra process；Wave field depth continuation of the present invention is completed in real number field, and complex operation is avoided in communication process, improves computational efficiency.
Description
Technical Field
The invention belongs to the technical field of seismic exploration, and relates to an offset imaging method, which mainly utilizes a space-time localization one-way wave operator to realize settlement method offset of an observation system.
Background
The proposal of various localized decomposition methods is a revolutionary breakthrough in the field of modern signal processing, and is mainly applied to a time-frequency analysis technology for analyzing non-stationary time signals at first for describing the localized properties of the signals in time and frequency. When relating to the propagation and imaging of seismic waves, localization has different physical meanings, mainly in space and in the direction of propagation. The localized decomposition can be independent of the wave propagation process, and only the wave field is decomposed to extract the local angle information of the wave field, such as the illumination analysis and aperture compensation technique based on the two-way wave, or the reflection angle gather extraction technique. Another very attractive application of spatial localization is for wave propagation, enabling high-precision propagation and imaging of seismic waves in laterally rapidly strongly varying media, such as the Beamlet migration (Beamlet migration) technique.
The beamlet migration is a double-domain migration method based on a single-pass wave fluctuation equation, the propagation process comprises Local wavenumber domain background propagation and spatial domain Local disturbance correction, and a beamlet migration algorithm based on a Gabor-Daubechies frame, a Local Cosine base (Local Cosine base) and a Local Harmonic base (Local Harmonic base) is mainly developed. The spatially localized beamlet propagation operator can adapt to transverse severe changes in velocity, so that local velocity disturbances are significantly reduced, and good seismic wave propagation accuracy can be maintained through simple local phase shift correction. When complex speed medium environment imaging is involved, such as the condition that high-speed salt exists, good precision can be obtained. However, after transforming the time-space domain wave field into the frequency domain by the wavelet beam migration, the method only performs localized decomposition along the space dimension, lacks the localization in time, and cannot flexibly adapt to the non-stationary change of the seismic data in time and space, so the effect of sparse representation of the wave field is not ideal.
The Dreamlet shift method adds time-frequency localization on the basis of the spatial localization of the beamlets, and is a fully localized shift method in time and space. The Dreamlet atom is composed of tensor products of time and space localized atoms, can be respectively realized by flexibly using different local decomposition methods, and mainly realizes a local index frame (time) -local cosine base (space) with the redundancy of 2 and an orthogonal local cosine base (time) -local cosine base (space) Dreamlet offset method at present. For wavefield decomposition, orthogonal dreamlet decomposition is equivalent to multi-dimensional local cosine transform, such as Wang and Wu (2000) applying semi-adaptive multi-dimensional local cosine transform (time-adaptive or space-adaptive) to two-dimensional post-stack seismic data compression and testing migration effects at different compression rates. The result shows that the two-dimensional local cosine transform can still keep the effective components of the seismic signals under the high compression rate, but the compression is transparent relative to the transmission, namely the compressed data is reconstructed and then subjected to migration imaging. The Dreamlet migration simultaneously realizes space-time local decomposition transformation on a wave field and a single-pass wave propagation operator, and the space-time local decomposition coefficient of the seismic wave field is directly propagated, so that not only can seismic data be effectively represented, but also the space-time local single-pass wave Dreamlet propagation operator can flexibly and variously operate the wave field in a time correlation manner in the migration process.
The one-way wave migration method is a forward algorithm, the wave field is always stepped and extended along a specific direction (such as the depth direction), and the energy in the wave field correspondingly flows along a fixed direction of a time axis (such as a zero time point). At each depth, imaging conditions are applied to the extended wavefield to obtain reflection intensity information, returning the reflected energy to its spatial location. In the one-way wave migration method, the energy that has been "parked" no longer carries the reflection information of the deeper structure, and does not need to propagate further downward. However, in all frequency domain migration algorithms, due to the continuous phase shift operation and the periodicity of the Fourier transform in the frequency domain, even after the energy reaches the zero time of the seismic wavefield, it collapses back to the other end of the time axis and continues to propagate down the wavefield. The Dreamlet propagation operator uses temporal localization and when extending the wavefield, if the Dreamlet coefficient is propagated past zero to a negative time, it is discarded. Therefore, with the increase of the depth, the signals of the outgoing wave field are more and more, the effective length of the time shaft is shorter and shorter, the wave field coefficient needing to be extended is correspondingly reduced, and the migration efficiency can be improved. In the shot domain dreamlet migration, the source and detector wavefields are extended to all locations in space, respectively, and the reflection value at each point is derived from the source and detector wavefields. In general, the reflection energy will not occur at zero time when it reaches its reflection location, and even if the dreamlet operator discards all the signals that arrive at negative time, a large amount of useless data will stay in the wavefield. While the seismic travel time can be calculated to further discard the unwanted data, a rough estimate of travel time does not completely discard the unwanted data, and calculating the travel time accurately adds additional computation.
Disclosure of Invention
In view of the above problems, an object of the present invention is to provide a migration imaging method, in which a spatio-temporal localized propagation operator is combined with an observation system settlement concept, a wave field dreamlet coefficient crossing zero is automatically discarded in a dreamlet migration operation during wave field continuation, all wave field values arriving at zero are extracted as reflection energy and cannot be continuously propagated, the extended wave field only contains reflection information of an observation system substructure, and the number of dreamlet coefficients required to be propagated is automatically reduced without additional processing.
In order to realize the purpose, the invention adopts the following technical scheme:
an offset imaging method for obtaining seismic data u (t, x) s ,x r ) Performing dreamlet orthogonal decomposition to obtain a dreamlet coefficient of a wave field; at each depth, firstly, a dreamlet operator is used for extending a common wavelet beam source, and a dreamlet operator is used for subsiding the common wavelet beam detector for the wave field after extension; and transforming the wave field to a frequency domain for local disturbance correction, and then applying an imaging condition to extract the wave field amplitude at the zero offset and zero moment as an image value.
The invention is further improved in that the method specifically comprises the following steps:
1) Constructing an orthogonal dreamlet basis function:
in the spatial domain, the starting position isA space width ofAnd a local cosine base b with a local wave number index m mn (x) The expression is as follows:
in the formula B n (x) Is shown in the intervalBell function with smooth upper surface and tight support, L n Represents the width of the nth window, M = 0.., M-1, M represents the number of wave numbers of the local cosine basis;
expression b of local cosine base in wavenumber domain mn (ξ) is as follows:
wherein B is n (xi) is the clock function B n (x) Fourier transform of (2), and represents the local wavenumber;
similar in the time domain, the local cosine base is started from the timeLength of timeAnd local frequencyRepresents the time domain b of ij (t) and frequency domain b ij The expressions (ω) are respectively:
wherein B is j (ω) is B j (t) Fourier transform and has
Separately decomposing the source and detector wavefields, dreamlet atoms, at two-dimensional shot domain migrationFrom the time domain b ij (t) and spatial domain b mn (x) The tensor product of the local cosine basis constitutes:
the branch of the Dreamlet atom is the Cartesian product of the time and space bell-type functions, the horizontal lines on the letters represent local variables with respect to the time window, in equation (5)Which is indicative of a local time of day,the local frequency is represented by a local frequency,indicating a local position andrepresenting a local wavenumber;
2) Reading all received seismic data into a memory and according to a common source point x s And a common detector point x r Spread, denoted as seismic data u (t, x) s ,x r ) Wherein t represents a recording time;
seismic data u (t, x) s ,x r ) Dreamlet atom of (1)From the source beamlet b mn (x s ) Detector beamlet b bq (x r ) And time domain b ij (t) a tensor product;
wherein the local frequencySource local wavenumberLocal wave number of detector The spatial location of the source beamlets is indicated,representing the detector beamlet spatial position;
wave field u (t, x) of each layer s ,x r ) By dreamlet atomIs represented as follows:
in order to localize the parameter index,<,> represents the inner product operation, defined as follows:
is thatIn a simplified representation of the same, the reference numbers,representing the source beamlet b mn (x s ) Detector beamlet b bq (x r ) And time frequencyA lower wavefield dreamlet coefficient; in the process of migration, only coefficient values higher than a certain threshold value are propagated, and small dreamlet coefficients are abandoned, namely, a compression propagation process is realized; based on the local disturbance theory, all the coefficients are extended to the next layer at the speed in the spatial local window;
3) Dreamlet domain observation system sedimentation method offset imaging
In the space-time domain, the expression of the double square root operator of the source-detector domain is
Wherein v (x) s ) And v (x) r ) Velocity at the source and detector locations, respectively; the dual square root operator extends both the source and detector down simultaneously, and both propagate down sequentially, equation (9) is equivalently split into two single square root equations:
equations (10) and (11) have the same form, for seismic data u (t, x) r ；x s ) Respectively extending a wave field of the detector and a source wave field downwards; deducing a propagation operator of the source-end dreamlet observation system by the formula (10), deducing a propagation operator of the detector dreamlet by the formula (11), and propagating the source-end dreamlet observation system in the case of offsetThe operator and the detector dreamlet propagation operator respectively act on the common source small beam and the common detector small beam of the data in sequence, the wave field is transformed to the frequency domain for local disturbance correction, and then the wave field amplitude at the zero offset and the zero moment is extracted as an image value by applying the imaging condition.
A further improvement of the present invention is that in step (1) all seismic data received is edge-extended with a zero value to enlarge the computed aperture.
The invention is further improved in that the concrete process of deriving the dreamlet observation system propagation operator by the formula (10) is as follows:
solution of formula (10) z+Δz (t,x r ；x s ) Is composed of
Δ z is the depth continuation step, where
P S The single square root propagation operator is separated into propagation operators P under local uniform velocity according to local disturbance theory S0 And local phase screen correction P S1 :
And is provided with
WhereinFor detector beamlet spatial positionLocal constant velocity reference velocity inside;
after passing through the medium, the dreamlet atoms evolve and follow the wave equation; after propagation, a single dreamlet atom can diffuse into other space-time windows of the phase space; the dreamlet atoms after propagation distortion are decomposed again to form a dreamlet propagation matrix
Wherein the subscripts with a left-falling stroke are dreamlet coefficient parameters propagated to the depth of a new layer;
corresponding to dreamlet propagation matrix atoms in observation system settling offset imagingIs composed of
dreamlet propagation matrix atomRepresenting the coupling coefficient weight between different dreamlet atoms before and after propagation;
after formula (6) is substituted for formula (17), there are
WhereinThe propagation operator is a single square root dreamlet propagation operator and is used for extending a wave field of the detector;
equation (18) is equivalent to the partial cosine base due to its orthogonality
WhereinRepresenting a source beamlet energy normalization coefficient;
the formula (19) shows that the dreamlet extension operator settled by the observation system is completely the same as the shot domain dreamlet operator, and the wave field of the common wavelet beam source is only diffused among the dreamlet coefficients of the detector after being propagated, and can not be propagated to other wavelet beam sources;
the analytical solution of the wave equation in the frequency wave number domain and the local cosine base expression of the frequency and wave number domain are substituted into an expression (19) to obtain
The dreamlet propagation operator is obtained by analytic derivation of a complex-domain frequency wave number domain one-way wave fluctuation equation, and finally only a real number part needs to be stored.
The invention is further improved in that the threshold is set to 10-4 times of the maximum value of the absolute value of the dreamlet coefficient corresponding to each layer of depth, namely, the dreamlet coefficient smaller than the threshold is ignored and does not participate in the propagation.
Compared with the prior art, the invention has the following beneficial effects:
(1) According to the method, orthogonal dreamlet decomposition is adopted for data, the total number of coefficients after decomposition is the same as the original data amount, but the dreamlet coefficients are distributed on a certain frequency wave number in a concentrated manner due to the energy compression capacity of local cosine transform, so that the pre-stack seismic data can be expressed in a sparse manner;
(2) According to the invention, orthogonal dreamlet decomposition is adopted for the double square root propagation operator, so that the gun domain dreamlet propagation operator can be applied to realize settlement offset of the observation system;
(3) According to the method, the Dreamlet propagation operator utilizes the localization on time, and when a wave field is extended, if the Dreamlet coefficient crosses zero time and reaches negative time after propagation, the Dreamlet coefficient is discarded and ignored, so that useless data can be automatically removed in the wave field extension process, and additional processing is not needed;
(4) The wave field depth continuation is completed in a real number domain, so that complex operation is avoided in the transmission process, and the calculation efficiency is improved.
Drawings
FIG. 1 is a schematic diagram of arrangement of settlement method prestack data of a two-dimensional Marmousi model observation system. The horizontal coordinate is the position coordinate of the detector, the vertical axis is the position coordinate of the source point, the zero value is used for expanding the edge and widening the calculation aperture, and the middle part is the position of the earth surface for receiving the seismic data.
Fig. 2 shows the co-wavelet source low wavelet number dreamlet coefficients in the surface received data dreamlet decomposition coefficients under the source-detector coordinates of the two-dimensional Marmousi model.
Fig. 3 is an enlarged view of the box numbered 1 in fig. 2.
Fig. 4 is an enlarged view of the box numbered 2 in fig. 2.
Fig. 5 is an enlarged view of the box numbered 3 in fig. 2.
Fig. 6 shows the co-wavelet source intermediate beamlet wavenumber dreamlet coefficients in the surface received data dreamlet decomposition coefficients under the source-detector coordinates of the two-dimensional Marmousi model.
Fig. 7 is an enlarged view of the box numbered 1 in fig. 6.
Fig. 8 is an enlarged view of the box numbered 2 in fig. 6.
Fig. 9 is an enlarged view of the box numbered 3 in fig. 6.
Fig. 10 shows the co-wavelet source high beamlet wave number dreamlet coefficients in the surface received data dreamlet decomposition coefficients under the source-detector coordinates of the two-dimensional Marmousi model.
Fig. 11 is an enlarged view of the box numbered 1 in fig. 10.
Fig. 12 is an enlarged view of fig. 10 at box 2.
Fig. 13 is an enlarged view of the box numbered 3 in fig. 10.
FIG. 14 shows co-beamlet detector gather low beamlet wave number dreamlet coefficients in surface receive data dreamlet decomposition coefficients at two-dimensional Marmousi model source-detector coordinates.
Fig. 15 is an enlarged view of the box numbered 1 in fig. 14.
Fig. 16 is an enlarged view of fig. 14 at box 2.
Fig. 17 is an enlarged view of the box numbered 3 in fig. 14.
FIG. 18 shows the co-beamlet detector gather median beamlet wave number dreamlet coefficients in the surface receive data dreamlet decomposition coefficients at the source-detector coordinates of the two-dimensional Marmousi model.
Fig. 19 is an enlarged view of the box numbered 1 in fig. 18.
Fig. 20 is an enlarged view of fig. 18 at box 2.
Fig. 21 is an enlarged view of the box numbered 3 in fig. 18.
FIG. 22 shows co-beamlet detector gather high beamlet wave number dreamlet coefficients in surface receive data dreamlet decomposition coefficients at source-detector coordinates in a two-dimensional Marmousi model.
Fig. 23 is an enlarged view of the box numbered 1 in fig. 22.
Fig. 24 is an enlarged view of fig. 22 at box 2.
Fig. 25 is an enlarged view of the box numbered 3 in fig. 22.
FIG. 26 is a two-dimensional SEG/EAGE post-stack seismic data and its dreamlet decomposition. Wherein, (a) is post-surface-stack data; (b) spatio-temporal window partitioning of seismic data for Dreamlet transform; (c) Dreamlet transform coefficients; (d) displaying the seismic signals adjacent to the window for amplification; and (e) displaying the dreamlet coefficients of the adjacent window for enlargement.
FIG. 27 is a plot of wavefields at different depths in two-dimensional SEG/EAGE post-stack data. The seismic map is placed in superposition with the velocity model, and the zero time of the seismic map is located at the corresponding depth position. Wherein, (a), (c) and (e) are extended wave fields of Dreamlet operators; and (b), (d) and (f) are Beamlet operator continuation wave fields.
FIG. 28 is the two-dimensional SEG/EAGE post-stack migration imaging results. Wherein, (a) is a Dreamlet imaging result; (b) The length of a time axis is 5 seconds for the result of the Beamlet offset imaging; (c) The time axis length zero padding is extended to 8.192 seconds for the Beamlet offset imaging results.
FIG. 29 shows the two-dimensional Marmousi model pre-stack depth imaging results. Wherein, (a) is a Marmousi offset velocity model; (b) is the Dreamlet shot field offset; and (c) is Dreamlet observation system sedimentation method offset.
FIG. 30 shows the migration of a co-point source gather wavefield by dreamlet observation system settlement: wherein (a) is the earth's surface; (b) 1.0km depth; (c) 2.0km depth; and (d) 3.0km depth.
Fig. 31 is a schematic diagram of the change of the sedimentation method offset dreamlet coefficient of the two-dimensional Marmousi model observation system with the depth.
Detailed Description
The invention is described in further detail below with reference to the accompanying drawings:
1) Constructing orthogonal dreamlet basis functions
The local cosine transform is a localized discrete cosine transform whose basis function is the product of a smooth, tightly supported clock function and a cosine function. The local cosine transform has many properties similar to windowed fourier transform, and fast fourier transform can be used to implement fast algorithms.
In the spatial domain, the starting position isA space width ofAnd a local cosine base b having a local wave number index of M (M = 0.., M-1, M represents the number of wave numbers of the local cosine base) mn (x) The expression is as follows:
in the formula B n (x) Is shown in the intervalBell function with smooth upper support, L n The representation shows the width of the nth window.
Expression b of local cosine base in wavenumber domain mn (ξ) is as follows:
wherein B is n (xi) is the clock function B n (x) Fourier transform of (2), and the local wavenumber is indicated.
Similar to in the time domain, the local cosine base may be started from the timeLength of timeAnd local frequencyRepresents the time domain b of ij (t) and frequency domain b ij The expressions (ω) are respectively:
wherein B is j (ω) is B j (t) Fourier transformAnd is provided with
Separately decomposing the source and detector wavefields, dreamlet atoms, at two-dimensional shot domain migrationBy the time domain b ij (t) and spatial domain b mn (x) The tensor product of the local cosine basis constitutes:
the set of Dreamlet atoms is the Cartesian product of the time and space bell-type functions (Cartesian product). The horizontal lines on the letters in the present invention represent local variables with respect to the time-space window, as in formula (5)The local time is represented by the time of the local,the local frequency is represented by a local frequency,indicating a local position andthe local wavenumber is indicated.
2) Multidimensional prestack seismic data dreamlet decomposition
The wave field continuation of each step of the traditional observation system settlement method is divided into two parts: firstly, extending all wave fields of the detector downwards for each point source gather; second, for each detector gather, all source wavefields are extended downward. For wavefield decomposition, the dreamlet decomposition continues along the source, detector and time records, which is a fully localized decomposition of the wavefield. The source is decomposed into a beamlet source and the detector is decomposed into a beamlet detector. After decomposition, the wavefield propagation also changes from common point source, common point detector to run with common beamlet source and common beamlet detector gathers.
And (3) respectively imaging the single cannons by shot domain pre-stack depth migration, superposing migration results of the cannons to form a final migration section, and independently processing data of the cannons. However, the seismic wave propagation and migration imaging are a linear process, and all the received shot data can be propagated at the same time, namely the observation system settlement method migration. Due to the multiple coverage of the earth surface by multiple guns and multiple offset distances, the earth surface prestack seismic records in the two-dimensional space are three-dimensional continuous data volumes (5-dimensional data volumes in the three-dimensional case). Reading all received seismic data into a memory before beginning of settlement method migration of an observation system, and reading all received seismic data into the memory according to a common source point x s And a common detector point x r Permutation, denoted u (t, x) s ,x r ) Where t represents a recording time. FIG. 1 is a schematic diagram of a two-dimensional Marmousi model prestack seismic trace data arrangement. The horizontal sampling interval of the model is 12.5m, the total number of the shots is 240, and the shot distance is 25m. The time sampling interval was 0.004 seconds, the recording length was 3 seconds. In fig. 1, the abscissa is the detector position coordinate, the ordinate is the source point position coordinate, the time axis is perpendicular to the source-detector coordinate plane, and similar to shot domain calculation, all seismic data received are subjected to zero value edge expansion to expand the calculation aperture.
In observation system sedimentation method migration under a source-detector coordinate system, complete localization decomposition of seismic data (three-dimensional data volume) u (t, x) is realized s ,x r ) Dreamlet atom of (1)From the source beamlet b mn (x s ) Detector beamlet b bq (x r ) And time domain b ij (t) a tensor product;
wherein the local frequencyLocal wavenumber of sourceAnd detector local wavenumber The spatial location of the source beamlets is indicated,representing the detector beamlet spatial position.
Wave field u (t, x) of each layer s ,x r ) Can be composed of dreamlet atomsIs represented as follows:
in order to localize the parameter index,<,> represents the inner product operation, defined as follows:
is thatIn a simplified representation of the above-described embodiments,representing the source beamlet b mn (x s ) Detector beamlet b bq (x r ) And time frequencyThe wavefield dreamlet coefficient of. FIGS. 2-25 show a common beamlet source for two-dimensional Marmousi model surface receiving dataAnd-beamlet detectorGather dreamlet coefficient. The space window and the time window are both fixed in length and respectively comprise 16 sampling points, the upper left corner in each time window is low-frequency low wave number, and the lower right corner in each time window is high-frequency high wave number. In each co-beamlet trace set, three spatial and temporal windows are randomly selected to magnify the display coefficient distribution details. Due to the fact that orthogonal dreamlet decomposition is adopted, the total number of coefficients after decomposition is the same as the original data amount, but due to the energy compression capacity of local cosine transformation, the dreamlet coefficients are distributed on a certain frequency wave number in a concentrated mode. At offset, only coefficient values above a certain threshold are propagated, while small dreamlet coefficients are discarded, i.e. a compact propagation process is implemented. Based on the local perturbation theory, all coefficients continue to the next layer at a speed within the spatial local window.
The dreamlet transformation in the observation system settlement method is complete localized decomposition of a wave field, and local cosine basis transformation is respectively carried out on a source point space coordinate, a detector space coordinate and a time axis. The dreamlet coefficient of the decomposed seismic data is concentrated in a low-frequency and low-wave-number part and is sparse representation of the seismic data.
3) Dreamlet domain observation system sedimentation method offset imaging
In the space-time domain, the source-detector domain dual square root operator (DSR) is expressed as
Wherein v (x) s ) And v (x) r ) Respectively the velocity at the source and detector locations. The dual square root operator extends both the source and detector down simultaneously, but they can also be propagated down in sequence, respectively, then equation (9) can be equivalently split into two Simple Square Root (SSR) equations:
equations (10) and (11) have the same form, and are applied to seismic data u (t, x) r ；x s ) The detector wave field (common source) and the source wave field (common detector) are extended downwards respectively. In the invention, equation (10) is taken as an example to derive the dreamlet observation system propagation operator.
Solution of formula (10) z+Δz (t,x r ；x s ) Is composed of
Δ z is the depth continuation step, where
Is a formal single square root one-way operator, and cannot be solved directly in the time-space domain, v (x) s ) And when the transverse direction changes, the solution cannot be converted into the wave number domain.
According to the local perturbation theory, the single square root propagation operator (13) can be separated into a propagation operator P under local uniform velocity S0 And local phase screen correction P S1 :
And is provided with
WhereinFor the detector beamlet spatial positionLocal constant velocity reference velocities within.
After passing through the medium, the evolution of dreamlet atoms must follow the wave equation. After propagation, a single dreamlet atom will diffuse into the other spatio-temporal window of the phase space. The dreamlet atoms after propagation distortion are decomposed again to form a dreamlet propagation matrix
Where the primed subscripts are the dreamlet coefficient parameters that propagate to the new layer depth.
Thus corresponding to dreamlet propagation matrix atoms in observation system settling offset imagingIs composed of
It represents the coupling coefficient weights between different dreamlet atoms before and after propagation.
After the expression of dreamlet, namely the expression (6) is substituted into the expression (17), the expression
WhereinIs a square root dreamlet propagation operator used to extend the detector wavefield.
Equation (18) is equivalent to the partial cosine base due to its orthogonality
WhereinRepresenting the source beamlet energy normalization coefficient.
It can be seen from equation (19) that the dreamlet extension operator of the observation system settlement is completely the same as the shot domain dreamlet operator, and the wave field of the common wavelet beam source is only diffused among the detector dreamlet coefficients after being propagated, and is not propagated to other wavelet beam sources.
The analytic solution of the wave equation in the frequency wave number domain and the local cosine base expression of the frequency and wave number domain are substituted into an expression (19) to obtain
Although the dreamlet propagation operator is obtained by analytic derivation of a complex domain frequency wave number domain one-way wave equation, only the real number part needs to be stored finally.
The dreamlet propagation operator of the common detector set out by the formula (11) is completely the same as the source end, when the offset occurs, the dreamlet operator (formula (20)) of the source end is respectively and sequentially acted on the common source small beam and the small beam of the common detector of the data, a wave field is converted into a frequency domain for local disturbance correction, then, the imaging condition is applied, and the amplitude of the wave field at the zero offset and the zero moment is extracted as an image value.
In the invention, only coefficient values higher than a certain threshold value are transmitted during the offset, and small dreamlet coefficients are discarded, namely, a compression transmission process is realized. In the invention, the threshold value of each depth is selected to be set as the depth dre of the layerMaximum value of the absolute value of the coefficient of the amlet is 10 -4 The double, i.e., dreamlet coefficients less than the threshold are negligible and do not participate in propagation, thereby speeding up the offset speed.
Numerical simulation result
Two-dimensional SEG/EAGE post-stack depth migration
Post-stack seismic data migration is a special case of observation system sedimentation migration. The post-stack seismic data is the result of the pre-stack data after stacking processing, and can be regarded as seismic records acquired when the source and the detector are at the same position, and is also called as a zero offset profile. The post-stack data is of lower dimension and can be more easily used to characterize the dreamlet offset.
The post-stack recorded profile can also be understood by another theoretical model, i.e. the explosion reflector model. The model assumes that all the detectors are still on the surface, but the sources are distributed on various reflecting surfaces in the ground and are excited at zero time with the intensity proportional to the reflection coefficient of the reflecting layer, the propagation process ignores the multiple reflection and transmission effects of all interfaces, and the energy propagated from the reflecting surfaces to the surface is recorded. The reflection record received at the actual surface is a two-way trip from the surface to the reflective layer and back to the surface from the reflective layer. Therefore, in the explosion reflecting surface model, although it is necessary to convert the two-way travel profile into the one-way travel profile, the offset may be performed simply by setting the corresponding post-stack offset velocity to half the actual velocity.
FIG. 26 (a) is a two-dimensional SEG/EAGE synthetic post-stack seismic record with 1408 sample-space points, 12.5m (40 ft) sampling intervals, 5 seconds per record time length, and 8ms sampling rate. The grid shown in fig. 26 (b) is a spatiotemporal window partition of the data by dreamlet transform, and the length of the temporal and spatial windows is also 16 sample points. Fig. 26 (c) is a plot of dreamlet decomposition coefficients against input data (fig. 26 (a)), with horizontal axis representing space-wavenumber coordinates and vertical axis representing time-frequency, visually showing the sparseness of the wavefield in the dreamlet domain. Also for details, the seismogram of the adjacent window (fig. 26 (d)) and its corresponding dreamlet coefficient (fig. 26 (e)) are displayed enlarged.
The seismic migration imaging process is divided into two steps, wavefield reconstruction and application of imaging conditions. Wave field reconstruction obtains the wave field of the underground target area from the data received by the earth surface, and then reflection information is obtained by applying imaging conditions. FIG. 27 shows the wavefield resulting from the dreamlet and beamlet propagation operator depth prolongation with the zero time wavefield value placed at the corresponding depth on the velocity model. It is clear from the seismic map that the reflected energy received at the earth's surface, after wavefield continuation, appears at the zero time of the reflected spatial location wavefield. Therefore, the reflection information of the stratum can be obtained by extracting the zero-time amplitude of the extended wave field based on the imaging condition of the post-stack data of the explosion reflecting surface model. These energies that are relocated to spatial locations no longer carry reflected information from deeper formations and should not be propagated further after the application of the finished image conditions. However, in the frequency domain wavefield continuation wavefield, due to the finite time length fourier transform periodicity and the phase shift operation that persists in the frequency domain shift, the energy after reaching zero is equivalent to wrapping around to the other end of the wavefield time axis in the time-space domain and periodically, staying in the wavefield and continuing downward (fig. 27 (b), fig. 27 (d), fig. 27 (f)). The Dreamlet propagation operator has localization in both space and time, extending the wavefield by propagating the Dreamlet coefficients in each time-space window. For each space-time window, the dreamlet propagation operator only propagates each coefficient to the time window where the coefficient is relatively smaller, and when the dreamlet coefficient reaches the negative time after propagating beyond the zero moment, the dreamlet operator discards the atom. Thus, the wavefield energy arriving at time zero is discarded after further propagation, and only useful signals for imaging deeper structures remain in the wavefield. As the depth increases, the effective signal length also decreases (fig. 27 (a), 27 (c), and 27 (e)). Comparing the deepest dreamlet and the beamlet extended wave fields (fig. 27 (e), fig. 27 (f)), it can be clearly seen intuitively that a large amount of energy intensive but not contributing to the imaging is propagated during the extension of the frequency domain method.
In the post-stack data migration, the zero time of the wave field obtained by extending the depth of each layer is extracted to form a final imaging result, and the final imaging result is equivalent to the summation of the wave fields of all frequencies in a frequency domain. Fig. 28 (a) shows the dreamlet imaging result, and fig. 28 (b) shows the frequency domain beamlet imaging result. Careful comparison of the two sub-salt structure images shows that the dreamlet imaging result is less noisy, and even some false reflection structures appear in the frequency domain beamlet shift result due to the influence of the folding signal, as shown by the white arrows in fig. 28 (b). Of course, the effect of the periodic folding signal can also be reduced by adding the time axis length by zero padding at the end of each seismic record (Yilmaz, 2001), as shown in fig. 28 (c). However, zero padding can encrypt frequency domain samples, increase the number of frequencies to be transmitted, and thus increase the amount of calculation.
The two-dimensional Marmousi observation system settlement method prestack depth migration:
the two-dimensional Marmousi prestack migration velocity model is provided with 737 and 750 sampling points in the transverse direction and the depth direction respectively, and the sampling intervals are 12.5m and 4m respectively. The maximum and minimum speeds are 1500m/s and 5500m/s respectively, and 50 speeds are selected as reference speeds of background propagation in the test, and the sampling is uniform from the minimum speed to the maximum speed. Fig. 29 (a) is a Marmousi offset velocity model. Fig. 29 (b) shows the dreamlet observation system sedimentation imaging results. By way of comparison, fig. 29 (c) shows shot domain dreamlet offset imaging results. Comparing the offset results in fig. 29, it can be seen that the dreamlet observation system sedimentation method imaging result has the same quality as the shot domain dreamlet imaging result, and especially the model deep structure has higher resolution.
FIG. 30 is a common point source wavefield at different depths in a dreamlet observation system sedimentation migration. In shot domain migration, the extended detector wavefield is correlated with the source wavefield propagating downward from a virtual source placed on the surface to obtain reflection information. Typically, the zero time does not occur when the reflected energy reaches its depth. To improve computational efficiency, the present invention estimates travel time by maximum and minimum velocities in the model to cut off some of the time window energy that has been used to image shallow structures. And the observation system settlement method offset imaging condition is a zero-moment zero-offset wave field value, and the zero-moment wave field energy used for imaging can be automatically discarded in the subsequent transmission without additional processing.
Fig. 31 shows that the number of dreamlet coefficients of the two-dimensional Marmousi model changes with the depth when the model is shifted, and generally, the dreamlet coefficients decrease with the increase of the depth. For the observation system of the Marmousi model tested at this time, a gap exists between the source and the detector, so that in the initial propagation stage, the wave field can be diffused to generate equivalent sources which are distributed in the whole observation system, but after the wave field is settled to a certain depth, the dreamlet coefficient is rapidly reduced. In the test, the number of dreamlet coefficients on the earth surface is 14 times of the number of coefficients at the deepest layer, and the offset efficiency at the deepest layer is higher than that at the first layer by one order of magnitude.
Claims (5)
1. An offset imaging method characterized by subjecting seismic data u (t, x) s ,x r ) Performing dreamlet orthogonal decomposition to obtain a dreamlet coefficient of a wave field; at each depth, firstly, a dreamlet operator is used for extending a common wavelet beam source, and a dreamlet operator is used for subsiding the common wavelet beam detector for the wave field after extension; and transforming the wave field to a frequency domain for local disturbance correction, and then applying an imaging condition to extract the wave field amplitude at the zero offset and zero moment as an image value.
2. The offset imaging method according to claim 1, comprising the steps of:
1) Constructing an orthogonal dreamlet basis function:
in the spatial domain, the starting position isA space width ofAnd a local cosine base b with a local wave number index m mn (x) The expression is as follows:
in the formula B n (x) Is shown in the intervalBell function with smooth upper support, L n Denotes the width of the nth window, M = 0.., M-1, M denotes the number of wave numbers of the local cosine basis;
expression b of local cosine base in wave number domain mn (ξ) is as follows:
wherein B is n (xi) is the clock function B n (x) Fourier transform of (A), and representing a local wavenumber;
similar in time domain, the local cosine base starts from the timeLength of timeAnd local frequencyRepresents the time domain b of ij (t) and frequency domain b ij The expressions (ω) are:
wherein B is j (ω) is B j (t) Fourier transform of (t) and having
Separately decomposing the source and detector wavefields, dreamlet atoms, at two-dimensional shot domain migrationFrom the time domain b ij (t) and spatial domain b mn (x) The tensor product of the local cosine basis constitutes:
the branch of the Dreamlet atom is the Cartesian product of the time and space bell-type functions, the horizontal lines on the letters represent local variables with respect to the time window, in equation (5)The local time is represented by the time of the local,the local frequency is represented by a local frequency,indicating a local position andrepresenting a local wavenumber;
2) Reading all received seismic data into a memory and according to a common source point x s And a common detector point x r Spread, denoted as seismic data u (t, x) s ,x r ) Wherein t represents a recording time;
seismic data u (t, x) s ,x r ) Dreamlet atom of (1)Wavelet from sourceBundle b mn (x s ) Detector beamlet b bq (x r ) And time domain b ij (t) a tensor product;
wherein the local frequencySource local wavenumberLocal wave number of detector The spatial location of the source beamlets is indicated,representing the detector beamlet spatial position;
wave field u (t, x) of each layer s ,x r ) By dreamlet atomIs represented as follows:
in order to localize the parameter index,<,> represents the inner product operation, defined as follows:
is thatIn a simplified representation of the same, the reference numbers,representative source beamlet b mn (x s ) Detector beamlet b bq (x r ) And time frequencyA lower wave field dreamlet coefficient; in the process of offset, only coefficient values higher than a certain threshold value are transmitted, and small dreamlet coefficients are discarded, namely, a compression transmission process is realized; based on the local disturbance theory, all the coefficients are extended to the next layer at the speed in the spatial local window;
3) Dreamlet domain observation system sedimentation method offset imaging
In the space-time domain, the expression of the source-detector domain dual square root operator is
Wherein v (x) s ) And v (x) r ) Velocity at the source and detector locations, respectively; the dual square root operator extends the source and detector simultaneously downward, and both propagate downward in sequence, equation (9) splits equivalently into two single square root equations:
equations (10) and (11) have the same form, and are applied to seismic data u (t, x) r ；x s ) Is divided intoRespectively extending a wave field and a source wave field of the detector downwards; a source-end dreamlet observation system propagation operator is deduced by a formula (10), a detector dreamlet propagation operator is deduced by a formula (11), the source-end dreamlet observation system propagation operator and the detector dreamlet propagation operator respectively and sequentially act on a common-source small beam and a common detector small beam of data in the process of shifting, a wave field is converted into a frequency domain for local disturbance correction, and then the amplitude of the wave field at zero offset and zero time is extracted as an image value by applying an imaging condition.
3. An offset imaging method according to claim 1 wherein all seismic data received in step (1) is marginalized to enlarge the computational aperture.
4. The offset imaging method as claimed in claim 1, wherein the specific process of deriving the dreamlet observation system propagation operator from equation (10) is:
solution u of formula (10) z+Δz (t,x r ；x s ) Is composed of
Δ z is the depth continuation step, where
P S The single square root propagation operator is separated into propagation operators P under local uniform velocity according to local disturbance theory S0 And local phase screen correction P S1 :
And is provided with
WhereinFor the detector beamlet spatial positionLocal constant speed reference speed within;
after passing through the medium, the dreamlet atoms evolve and follow the wave equation; after propagation, a single dreamlet atom can diffuse into other space-time windows of the phase space; decomposing the dreamlet atoms after propagation distortion again to form a dreamlet propagation matrix
Wherein the subscripts with a skimming are dreamlet coefficient parameters propagated to the depth of a new layer;
corresponding to dreamlet propagation matrix atoms in observation system settlement offset imagingIs composed of
dreamlet propagation matrix atomRepresenting the coupling coefficient weight between different dreamlet atoms before and after propagation;
after the formula (6) is substituted for the formula (17), there are
WhereinThe propagation operator is a single square root dreamlet propagation operator and is used for extending a wave field of the detector;
equation (18) is equivalent to the orthogonality of the local cosine basis
WhereinRepresenting a source beamlet energy normalization coefficient;
the formula (19) shows that the dreamlet extension operator settled in the observation system is completely the same as the shot domain dreamlet operator, and the wave field of the common wavelet beam source is only diffused among the dreamlet coefficients of the detector after being propagated, and cannot be propagated to other wavelet beam sources;
the analytic solution of the wave equation in the frequency wave number domain and the local cosine base expression of the frequency and wave number domain are substituted into an expression (19) to obtain
The dreamlet propagation operator is obtained by analytic derivation of a complex domain frequency wave number domain one-way wave fluctuation equation, and finally only a real number part needs to be stored.
5. The offset imaging method as claimed in claim 1, wherein the threshold is set to 10 of the maximum of the absolute value of the dreamlet coefficient corresponding to each layer depth -4 Double, i.e., dreamlet coefficients less than the threshold are ignored and do not participate in the propagation.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710780533.XA CN107632321B (en) | 2017-09-01 | 2017-09-01 | A kind of offset imaging method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710780533.XA CN107632321B (en) | 2017-09-01 | 2017-09-01 | A kind of offset imaging method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107632321A true CN107632321A (en) | 2018-01-26 |
| CN107632321B CN107632321B (en) | 2018-12-07 |
Family
ID=61100264
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710780533.XA Active CN107632321B (en) | 2017-09-01 | 2017-09-01 | A kind of offset imaging method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107632321B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110988869A (en) * | 2019-11-22 | 2020-04-10 | 中国科学院电子学研究所 | Imaging method and device based on MIMO array |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102890290A (en) * | 2012-09-25 | 2013-01-23 | 中国石油天然气股份有限公司 | Prestack depth migration method under condition of undulating surface |
| CN106646595A (en) * | 2016-10-09 | 2017-05-10 | 电子科技大学 | Earthquake data compression method based on tensor adaptive rank truncation |
-
2017
- 2017-09-01 CN CN201710780533.XA patent/CN107632321B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102890290A (en) * | 2012-09-25 | 2013-01-23 | 中国石油天然气股份有限公司 | Prestack depth migration method under condition of undulating surface |
| CN106646595A (en) * | 2016-10-09 | 2017-05-10 | 电子科技大学 | Earthquake data compression method based on tensor adaptive rank truncation |
Non-Patent Citations (2)
| Title |
|---|
| 吴帮玉 等: "基于局部余弦基小波束的观测系统沉降法叠前深度偏移", 《地球物理学报》 * |
| 吴帮玉 等: "时-空局域化地震波传播方法：Dreamlet叠前深度偏移", 《地球物理学报》 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110988869A (en) * | 2019-11-22 | 2020-04-10 | 中国科学院电子学研究所 | Imaging method and device based on MIMO array |
| CN110988869B (en) * | 2019-11-22 | 2022-02-11 | 中国科学院电子学研究所 | Imaging method and device based on MIMO array |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107632321B (en) | 2018-12-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| KR101549388B1 (en) | Prestack elastic generalized-screen migration method for seismic multicomponent data | |
| Chen et al. | Iterative deblending of simultaneous-source seismic data using seislet-domain shaping regularization | |
| Chattopadhyay et al. | Imaging conditions for prestack reverse-time migration | |
| KR101548976B1 (en) | Estimation of soil properties using waveforms of seismic surface waves | |
| US10705233B2 (en) | Determining a component of a wave field | |
| CN109669212B (en) | Seismic data processing method, stratum quality factor estimation method and device | |
| CN106896409B (en) | A kind of varying depth cable ghost reflection drawing method based on wave equation boundary values inverting | |
| CN110579795B (en) | Joint velocity inversion method based on passive source seismic waveform and reverse-time imaging thereof | |
| CN107884829A (en) | A kind of method for combining compacting shallow sea OBC Multiple Attenuation in Seismic Data | |
| Stoffa et al. | Plane-wave depth migration | |
| Hayashi et al. | CMP spatial autocorrelation analysis of multichannel passive surface-wave data | |
| CN111856577B (en) | Method for reducing calculation amount of reverse-time migration earth surface offset gather | |
| Yin et al. | Improving horizontal resolution of high-frequency surface-wave methods using travel-time tomography | |
| CN113885079A (en) | Elastic wave field decoupling-based high-precision multi-azimuth reverse-time seismic source imaging method | |
| CN113805237A (en) | Method and system for offset land cross-spread seismic using compressed sensing models | |
| US11143774B2 (en) | Method and system for separating blended seismic data | |
| CN111999764B (en) | Method for constructing least square reverse time migration under salt based on time-frequency domain objective function | |
| Kazei et al. | Amplitude-based DAS logging: Turning DAS VSP amplitudes into subsurface elastic properties | |
| CN107632321A (en) | A kind of offset imaging method | |
| Hu et al. | An iterative focal denoising strategy for passive seismic data | |
| EA037490B1 (en) | Method of marine seismic acquisition | |
| CN112505782B (en) | Interference reference plane reconstruction method and system for radiation mode correction in four-dimensional earthquake | |
| Yang* et al. | Prestack depth migration method using the time-space Gaussian beam | |
| Wu et al. | Attenuation compensation in multicomponent Gaussian beam prestack depth migration | |
| Chen et al. | Research on vertical cable seismic interferometry imaging |
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 |