CN117908095A

CN117908095A - Seismic wave imaging method and device

Info

Publication number: CN117908095A
Application number: CN202211233509.1A
Authority: CN
Inventors: 徐基祥
Original assignee: China National Petroleum Corp
Current assignee: China National Petroleum Corp
Priority date: 2022-10-10
Filing date: 2022-10-10
Publication date: 2024-04-19

Abstract

The invention discloses a seismic wave imaging method and device. The method comprises the steps of calculating an illumination dip angle according to an incident ray parameter and a scattered ray parameter of a scattering point aiming at the scattering point at each imaging point position; determining the common imaging aperture of the reflected wave and the scattered wave according to the maximum and minimum illumination dip angles, and establishing an illumination dip angle gather according to the imaging aperture and the illumination dip angle of a scattering point on the imaging channel for each imaging channel; performing generalized radon transformation of three dimensions of depth, radon parameters and stratum dip angle on the illumination dip angle gather to obtain a radon model with scattered wave and reflected wave separated; and respectively imaging the scattered wave and the reflected wave in the radon model to obtain a scattered wave imaging channel and a reflected wave imaging channel, and further obtaining a scattered wave imaging data body and a reflected wave imaging data body. The method can utilize weak scattered waves in the seismic waves and retain strong reflected waves, so that high-definition imaging of the seismic waves is realized.

Description

Seismic wave imaging method and device

Technical Field

The invention relates to the technical field of seismic data processing, in particular to a seismic wave imaging method and device.

Background

The seismic exploration is a method for deducing the morphology and the rock stratum property of an underground structure by utilizing the difference of elasticity and density of an underground medium and observing and analyzing the response of the underground elastic medium to manually excited seismic waves, has become the most important means for petroleum and natural gas resource exploration, and is widely applied to the aspects of coal field and engineering geological exploration, regional geological research, crust research and the like. At present, the seismic exploration mainly adopts reflected wave seismic exploration and consists of three parts of seismic data acquisition, processing and interpretation, wherein seismic data migration imaging represented by prestack depth migration is a key core link of seismic data processing.

The primary task of seismic exploration is to improve the capability of seismic data to distinguish underground structures and lithology, and the main object of reflected wave seismic exploration is a large underground sleeve (large scale) layered and transversely spread stratum, so that the method plays an irreplaceable important role in the conventional oil and gas exploration stage. Along with expansion of the oil and gas exploration field, the method changes from the construction oil and gas reservoir exploration to lithologic oil and gas reservoir exploration, changes from the middle shallow exploration to the middle-deep exploration, and extends from the sedimentary rock oil and gas reservoir exploration to the carbonate rock and igneous rock oil and gas reservoir exploration, and even the tight oil and gas reservoir exploration, so that higher requirements are provided for the seismic exploration precision. For non-structural hydrocarbon reservoirs such as carbonate, igneous rock, small fragments and lithology, the common feature of them is small-scale geologic bodies or irregular surfaces, which are difficult to form effective reflected waves, but generate a large number of scattered waves. The imaging definition of the underground geologic body can be improved by utilizing scattered waves, and the method has important practical significance for lithologic hydrocarbon reservoir exploration.

Near the wavelength scale of the seismic wave, velocity and density inhomogeneities of the subsurface medium produce scattered waves. Unlike conventional reflected waves, reflected wave energy is concentrated in a certain direction, and scattered wave energy diverges to the periphery, so that scattered wave energy is weaker than reflected wave energy, which makes scattered wave imaging much more difficult than reflected wave imaging. Broadly, seismic scattered waves include reflected waves that are the result of directional coherence enhancement of the scattered waves. Seismic data collected in the field are rearranged to obtain a common center point gather, namely a time-offset domain gather. On the common-center point gather, both the reflected wave in-phase axis and the scattered wave in-phase axis are in hyperbolic shapes, but the curvature of the scattered wave in-phase axis is obviously larger than that of the reflected wave in-phase axis. The conventional reflected wave seismic wave exploration regards reflected waves as effective signals, so that the reflected waves are leveled on the same phase axis through dynamic correction and then are horizontally overlapped, and a conventional reflected wave overlapped section is obtained. Obviously, on the common-center-point gather after dynamic correction, the reflection wave in-phase axis is leveled, and the scattering wave in-phase axis is generally curved and appears as an upward-convex in-phase axis. Due to the weak scattered wave energy, it is often suppressed as noise. The horizontal superposition profile is the result of the early stage of reflected wave exploration in the fifth sixties of the last century, and is generally suitable for oil and gas exploration of basins with simple structures such as horizontal stratum and the like, such as Daqing oilfield of Songliao basin. As the complexity of constructing hydrocarbon reservoir exploration increases, the horizontal stacking profile has failed to meet the exploration requirements, and seismologists have developed offset techniques. The essence of this technique is to return the oblique layer to the position, also for reflected waves. The development of a time shift from post-stack to pre-stack depth shift solves the more complex structural imaging problems, such as imaging of Bohai Bay basin segments and imaging of Tarim basin kuda front Liu Chongduan with high and steep structures. The pre-stack depth migration may generate a common-imaging-point gather, most commonly a depth-scatter-angle-domain common-imaging-point gather. Like the dynamically corrected common-center point gathers, the reflected wave phase axis is also straight and the scattered wave phase axis is also often convex up on the depth-scatter angle domain common-image point gathers. After the depth-scattering angle domain common imaging point gathers are horizontally overlapped, reflected waves are still imaged, and scattered waves are suppressed. Until the large-scale structure oil and gas reservoir is found almost, the application value of the reflected wave reaches the peak. The application value of scattered waves is only important when entering a small-scale lithologic oil and gas exploration stage.

In the early stage of the development of scattered wave imaging technology, even for data before imaging in the imaging process, a weighting factor is set to suppress reflected waves and strengthen the scattered waves, and finally the reflected waves are too strong to take effect. The primary task is now to protect the weakly scattered wave, which is no longer regarded as noise, but rather to use it as an effective signal. Therefore, the same phase axis of the scattered wave is leveled on the gather, and the scattered wave imaging can be ensured through horizontal superposition. Of course, the scattered wave event is flat and the reflected wave event appears as a concave event. In the context of complex construction, the same phase axis of the scattered wave can only be leveled on the depth-illumination tilt domain common-image point gathers. This is only the first step, which exploits the difference in radiation direction of the scattered wave and the reflected wave. Their energy differences have not been exploited. If depth-illumination dip domain common imaging point gathers are directly superimposed, the strong energy reflected waves will produce significant imaging noise, albeit with enhanced scattered waves. In reflected wave imaging, one limits the offset aperture to solve the problem of arcing such as imaging noise, and in fact, the problem of limited aperture has not been solved until now because it is determined by the inclination of the subsurface reflecting layer. Not only can weak scattered waves be utilized, but also strong reflected waves can be reserved, so that high-definition imaging of the seismic waves is realized, and the method is a technical problem to be solved urgently at present.

Disclosure of Invention

In view of the foregoing, the present invention has been made to provide a seismic wave imaging method and apparatus that overcomes or at least partially solves the foregoing problems, and that can realize high-definition imaging of seismic waves using both weak scattered waves and strong reflected waves in the seismic waves.

In a first aspect, an embodiment of the present invention provides a seismic wave imaging method, including:

An illumination inclination calculation step, namely calculating an illumination inclination of a scattering point according to an incident ray parameter and a scattered ray parameter of the scattering point aiming at the scattering point at each imaging point position of a research area;

An illumination dip angle gather establishing step, namely determining an imaging aperture common to reflected waves and scattered waves in the seismic waves according to the maximum and minimum illumination dip angles, and establishing an illumination dip angle gather according to the imaging aperture and the illumination dip angle of each scattering point on each imaging channel aiming at each imaging channel in a research area;

A Lato transformation step, namely, carrying out generalized Lato transformation of three dimensions of depth, lato parameters and stratum dip angle aiming at each illumination dip angle gather to obtain a Lato model with scattered wave and reflected wave separated;

imaging the scattered wave and the reflected wave in each radon model to obtain a scattered wave imaging channel and a reflected wave imaging channel;

And an imaging data body synthesizing step, wherein the scattered wave imaging data body is synthesized by the scattered wave imaging channels, and the reflected wave imaging data body is synthesized by the reflected wave imaging channels.

In a second aspect, an embodiment of the present invention provides a seismic imaging apparatus, including:

The illumination inclination angle calculation module is used for calculating the illumination inclination angle of the scattering points according to the incident ray parameters and the scattered ray parameters of the scattering points aiming at the scattering points at each imaging point position of the research area;

The illumination dip angle gather establishing module is used for determining the common imaging aperture of the reflection wave and the scattering wave in the seismic wave according to the maximum and minimum illumination dip angles, and establishing an illumination dip angle gather according to the imaging aperture and the illumination dip angle of each scattering point on the imaging channel aiming at each imaging channel of the research area;

The Lato transformation module is used for carrying out generalized Lato transformation of three dimensions of depth, lato parameters and stratum dip angle aiming at each illumination dip angle gather to obtain a Lato model with scattered wave and reflected wave separated;

The imaging channel imaging module is used for respectively imaging the scattered wave and the reflected wave in each radon model to obtain a scattered wave imaging channel and a reflected wave imaging channel;

and the imaging data body synthesis module is used for synthesizing the scattered wave imaging data body by the scattered wave imaging channels and synthesizing the reflected wave imaging data body by the reflected wave imaging channels.

In a third aspect, embodiments of the present invention provide a computer storage medium having stored therein computer executable instructions that when executed by a processor implement the above-described seismic wave imaging method.

In a fourth aspect, embodiments of the present disclosure provide a server, including: the system comprises a memory, a processor and a computer program stored in the memory and capable of running on the processor, wherein the processor realizes the seismic wave imaging method when executing the program.

The technical scheme provided by the embodiment of the invention has the beneficial effects that at least:

(1) According to the seismic wave imaging method provided by the embodiment of the invention, each imaging point is assumed to have one scattering point, and the illumination dip angle of each scattering point is calculated, so that an illumination dip angle gather is established; for each illumination dip angle gather, generalized radon transformation of three dimensions of depth, radon parameters and stratum dip angle is carried out, so that a radon model with scattered wave and reflected wave separated is obtained; the scattered wave and the reflected wave are separated, so that the two images are respectively formed. The method strengthens scattered wave and focused reflected wave, and realizes high-definition imaging of seismic wave.

(2) According to the seismic wave imaging method provided by the embodiment of the invention, the common imaging aperture of the reflected wave and the scattered wave in the seismic wave is determined according to the maximum and minimum illumination dip angles, so that the problem that the limited aperture required by the reflected wave imaging is difficult to determine is solved, and the problem that the calculation amount of the scattered wave imaging by using an infinite aperture is too large is solved.

Additional features and advantages of the invention will be set forth in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. The objectives and other advantages of the invention will be realized and attained by the structure particularly pointed out in the written description and claims thereof as well as the appended drawings.

The technical scheme of the invention is further described in detail through the drawings and the embodiments.

Drawings

The accompanying drawings are included to provide a further understanding of the invention and are incorporated in and constitute a part of this specification, illustrate the invention and together with the embodiments of the invention, serve to explain the invention. In the drawings:

FIG. 1 is a flow chart of a seismic imaging method according to a first embodiment of the invention;

FIG. 2 is a flowchart showing the implementation of step S11 in FIG. 1;

FIG. 3 is a flowchart showing a specific implementation of step S12 in FIG. 1;

FIG. 4 is a flowchart of a seismic imaging method according to a second embodiment of the invention;

FIG. 5 is a schematic view of illumination tilt gathers according to a second embodiment of the present invention;

FIG. 6 is a conventional imaging, reflected wave imaging trace and scattered wave imaging trace of the gather of FIG. 5;

FIG. 7 is a conventional imaging profile in a second embodiment of the invention;

FIG. 8 is a reflected wave imaging profile in a second embodiment of the invention;

FIG. 9 is a cross-section of scattered wave imaging in accordance with a second embodiment of the invention;

FIG. 10 is an imaging section of a scattered wave and reflected wave synthesized in a 10:1 ratio in accordance with a second embodiment of the present invention;

FIG. 11 is a schematic diagram of a seismic imaging apparatus according to an embodiment of the invention.

Detailed Description

Exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be embodied in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

It is to be understood that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. In addition, for numerical ranges in this disclosure, it is understood that each intermediate value between the upper and lower limits of the ranges is also specifically disclosed. Every smaller range between any stated value or stated range, and any other stated value or intermediate value within the stated range, is also encompassed within the invention. The upper and lower limits of these smaller ranges may independently be included or excluded in the range.

Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs. Although only preferred methods and materials are described herein, any methods and materials similar or equivalent to those described herein can be used in the practice or testing of the present invention. All documents mentioned in this specification are incorporated by reference for the purpose of disclosing and describing the methods and/or materials associated with the documents. In case of conflict with any incorporated document, the present specification will control.

In order to solve the problem that the scattered wave and the reflected wave cannot be simultaneously utilized in the seismic wave imaging in the prior art, the embodiment of the invention provides a seismic wave imaging method and device, which can utilize the weak scattered wave in the seismic wave and retain the strong reflected wave to realize the high-definition imaging of the seismic wave.

Example 1

The first embodiment of the invention provides a seismic wave imaging method, which is suitable for the fields of three lithology hydrocarbon reservoir exploration of clastic rock, carbonate rock and igneous rock, and also suitable for the fields of mineral exploration, engineering exploration and the like, and the flow is shown in a figure 1, and comprises the following steps:

step S11: and aiming at the scattering points at the positions of each imaging point of the research area, calculating the illumination inclination angles of the scattering points according to the incident ray parameters and the scattered ray parameters of the scattering points.

Referring to fig. 2, the method specifically includes the following steps:

Step S111: and according to a seismic wave velocity model of the research area, using a program function equation to respectively calculate ray paths from the shot point and the receiving point to each imaging point, and obtaining an incident ray parameter and a scattered ray parameter of a scattering point at the position of the imaging point.

Assuming scattering points at the positions of each imaging point in advance, under the condition that a velocity macro model is known, using a program function equation, the ray paths from the seismic data shot point and the receiving point to each imaging point (scattering point) underground are calculated.

The ray from the shot point to the scattering point is an incident ray, the downward inclination angle (the included angle with the plumb line) is alpha _SM, and the azimuth angle (the included angle with the x direction) is beta _SM. The rays from the receiving point to the scattering point are scattered rays (opposite to the actual scattered rays), and the downtilt angle is alpha _GM, and the azimuth angle is beta _GM.

Specifically, the x direction is located on the upper plane, and the specific direction can be flexibly set according to the requirement.

Step S112: and calculating the illumination inclination angle of the scattering point according to the incidence angle of the scattering point and the downward inclination angle and the azimuth angle of the scattering ray.

According to the incident ray of the scattering point and the downward inclination angle and azimuth angle of the scattering ray, the illumination inclination angle of the scattering point is calculated by the following calculation formula of the illumination inclination angle theta relative to the x direction:

in the above formula, α is a scattering angle, specifically, half of an included angle between an incident ray and a scattered ray, and the calculation formula is as follows:

The above formulas (1) and (2) are generic and apply to both two-dimensional and three-dimensional situations, even to steep and inverted configurations.

Step S113: and determining the positive and negative of the illumination inclination angle of the scattering point according to the inclination direction of the stratum where the scattering point is located.

The illumination inclination angle of the scattering point is positive when the stratum where the scattering point is located is inclined upwards (downwards left) relative to the x direction; when the stratum where the scattering points are located is declined right (tilted left upwards) relative to the x direction, the illumination tilt angle of the scattering points is negative.

Step S12: and determining the imaging aperture of the seismic waves, which is common to the reflected wave and the scattered wave, according to the maximum and the minimum illumination dip angles, and establishing an illumination dip angle gather according to the imaging aperture and the illumination dip angle of each scattering point on the imaging channel aiming at each imaging channel of the research area.

Referring to fig. 3, the method specifically includes the following steps:

step S121: and calculating a central ray for each scattering point by using a kinematic ray tracing method to obtain travel time from the shot point and the receiving point to the scattering point.

For each scattering point r= (x, y, z) in the ground, the central ray is calculated by a kinematic ray tracing method to obtain a shot point r _S and a receiving point r _G, which are respectively marked as tau _S＝(r_S, r) and tau _G＝(r_G, r when the travel from the shot point r to the point r is carried out, wherein r _S and r _G are both on the ground surface, the ground surface is represented by a symbol S, and the coordinate of the ground surface is represented as xi= (xi ₁,ξ₂).

Step S122: and calculating a paraxial ray by using a dynamic ray tracing method to obtain a geometric diffusion matrix, obtaining the amplitudes from shot points and receiving points to scattering points, obtaining a second partial derivative of the central ray coordinate during travel, and further obtaining Beylkin determinant h (r, ζ).

The amplitudes a _S and a _G from shot, receiver to scatter, are denoted as a _S＝A(r_S,r)、A_G＝A(r_G, r respectively).

Step S123: and taking the projection point of the imaging channel on the ground surface as a central point, and determining an integration surface according to the central point and the imaging aperture.

From the ground-measured seismic data D (t, r _S,r_G), where t is the recording time, the illumination dip gather at spatial location r= (x, y, z) can be calculated:

in the above formula (3), d (r, θ) is an illumination inclination angle gather, r is a scattering point, and θ is a variable illumination inclination angle; hit (r, θ) is the number of illumination times; The upper points represent the derivative over time; s is an integrating plane; ζ is the position of the point on the integration surface.

For weighting operators, in true amplitude/>Wherein A _S and A _G are amplitudes from the shot point and the receiving point to the scattering point, h (r, ζ) is Beylkin determinant,/>, respectivelyRepresenting the gradient operator sign, the scalar in the numerator is represented by taking its absolute value, and the vector in the denominator is represented by taking its modulus.

The illumination tilt gather d (r, θ) is a four-dimensional array, with the illumination tilt gather being a two-dimensional array d (z, θ) at each (x, y) location (imaging track). In the illumination dip angle gather d (z, θ), since the scattered wave diverges around, theoretically the integral surface S may take the entire acquisition surface of the probe region, and taking the scattered wave attenuation and the calculation amount into consideration, an integral surface S ₁ centered on (x, y) may be taken instead of S. This integration surface S ₁ is defined by the imaging aperture. Reflected wave imaging requires a finite aperture, while scattered wave imaging aperture may be infinite. Since the finite aperture of reflected wave imaging is centered at the intersection of the subsurface interface normal reflection line at the surface, this finite aperture is also related to the inclination of the reflection interface. When the seismic architecture is complex, the determination of a finite aperture has heretofore been an inefficient approach. Instead, the method does not require determining the imaging aperture of the reflected wave for each reflective interface in the subsurface, and only requires that the dip angle θ ₀ of all formations in the subsurface be contained within the illumination dip range, denoted as θ ₀∈[θ_min,θ_max, where θ _min and θ _max are the minimum and maximum illumination dip, respectively. In d (z, θ), the scattered wave phase axis is flattened and the reflected wave phase axis appears as a concave curve, the depth at which the apex of the curve lies is the depth of the reflecting interface, and the illumination tilt at the apex is the tilt of the reflecting interface.

The above formula (3) can be converted into formula (4):

in the above formula (4), d (z, θ) is an illumination inclination angle gather, z is a position variable in the vertical direction, and θ is a variable illumination inclination angle; hit (z, θ) is the number of illumination times; the points above D represent the derivation of time, t is the recording time of the seismic data, and r _S、r_G、τ_S and τ _G are the shot position, the receiving position, the travel time from shot to scattering and the travel time from receiving to scattering in sequence; s ₁ is an integration surface; ζ is the position of the point on the integration surface; In order for the weighting operator to be a weighted operator, Wherein A _S and A _G are amplitudes from a shot point and a receiving point to a scattering point, and h (z, ζ) is a Beylkin determinant.

Step S124: and establishing an illumination dip angle gather according to the integral surface of the imaging channel and the illumination dip angle of each scattering point on the imaging channel.

The illumination tilt angle gather is established by the above formula (4).

This step is not limited to the above scheme, and other schemes such as wave equation may be substituted.

Step S13: and aiming at each illumination dip angle gather, carrying out generalized radon transformation of three dimensions of depth, radon parameters and stratum dip angle to obtain a radon model with scattered wave and reflected wave separated.

The following generalized radon transform is performed on the illumination dip gather d (z, θ):

in the formula (5) of the present invention, For a three-dimensional radon model, χ represents the depth (km) of the three-dimensional radon model, p is the radon parameter (km),/>For stratum dip angle, generalized radon transform kernel function/>Can be expressed as:

The generalized radon transformation is directly calculated, and the obtained three-dimensional radon model has lower resolution and can generate obvious noise and false images. For this purpose the above generalized radon transform needs to be calculated by means of a sparse constraint inversion technique. This needs to start with the following generalized radon inverse transform:

to obtain a high resolution three-dimensional radon model, the following optimization problem needs to be solved:

In the formula (8), m is a three-dimensional radon model, d is an illumination dip angle domain gather, Fourier transform operator representing depth, L is defined by the kernel function/>Constructed generalized radon transform operator,/>Denote l ₁ norm,/>The square of the norm of l ₂ is represented, e the allowable residual. By least squares, this optimization problem can be transformed into the following system of linear equations:

In the formula (9) of the present invention, For the radon model, χ is depth, p is radon parameter,/>Is the formation dip angle; d (Z, θ) is the illumination tilt angle gather, Z is the position variable in the vertical direction, θ is the variable illumination tilt angle; w is the weight coefficient of the Ladong model; fourier inverse transform operator representing depth,/> A fourier transform operator representing depth; l is defined by the kernel function/>Constructed generalized radon transform operator,/>Is the complex conjugate transpose of L.

For each illumination dip angle gather, a radon model with scattered wave and reflected wave separated in three dimensions of depth, radon parameters and formation dip angle can be obtained by solving the above formula (9). Further, a radon model in which scattered waves are distributed in a region with a radon parameter of 0 and reflected waves are distributed in a region with a radon parameter of less than 0 is obtained.

Equation (9) may be solved using an iterative algorithm, for example, a preconditioned conjugate gradient algorithm PCG.

Step S14: and respectively imaging the scattered wave and the reflected wave in each radon model to obtain a scattered wave imaging channel and a reflected wave imaging channel.

In high resolution three-dimensional radon modelIn (c), the scattered wave is distributed over p=0, the reflected wave is focused on the formation dip position, and p= -z <0 of the reflected wave focus position. The scattered wave and the reflected wave are now distributed in different areas. Thus, the scattered wave and reflected wave imaging can be obtained by horizontal superposition in the divided regions. Because of the energy difference between the strongly reflected wave and the weakly scattered wave, they must be imaged independently to avoid flooding the weakly scattered wave. After independent imaging of scattered waves and reflected waves is obtained, the two imaging channels can be weighted and overlapped according to explanation requirements, and a final high-definition composite imaging channel can be obtained.

Step S15: a scattered wave imaging data volume is synthesized from the scattered wave imaging channels, and a reflected wave imaging data volume is synthesized from the reflected wave imaging channels.

It may also include obtaining a composite imaging data volume from the composite imaging modality.

According to the seismic wave imaging method provided by the embodiment of the invention, each imaging point is assumed to have one scattering point, and the illumination dip angle of each scattering point is calculated, so that an illumination dip angle gather is established; for each illumination dip angle gather, generalized radon transformation of three dimensions of depth, radon parameters and stratum dip angle is carried out, so that a radon model with scattered wave and reflected wave separated is obtained; the scattered wave and the reflected wave are separated, so that the two images are respectively formed. The method strengthens scattered wave and focused reflected wave, and realizes high-definition imaging of seismic wave.

The imaging aperture of the seismic waves, which is common to the reflected wave and the scattered wave, is determined according to the maximum and minimum illumination dip angles, so that the problem that the limited aperture required by the reflected wave imaging is difficult to determine is solved, and the problem that the calculation amount of the scattered wave imaging by utilizing an infinite aperture is too large is solved.

In order to meet the requirement of lithologic hydrocarbon reservoir exploration on the imaging precision of seismic waves, the scattered wave effect generated by small-scale geologic bodies must be exerted. As well as reflected waves generated by large scale structured interfaces, scattered waves are also effective signals. The embodiment of the invention aims to fully utilize scattered wave information, improve imaging precision of small-scale geologic bodies such as ancient river sand bodies, karst caves and the like, and simultaneously utilize reflected wave information imaging to construct a background so as to realize high-definition imaging of underground media.

Example two

The second embodiment of the present invention provides a specific implementation flow of a seismic wave imaging method, as shown in fig. 4, including the following steps:

Step S41: and calculating the illumination dip angle.

The velocity macro model in this embodiment has 772X1151×401 evenly distributed discrete points with voxel sizes of 25m×50m× 20m and depths from 0km to 8km. The shots and the receiving points are located on the ground surface, and the total of 29065 shots and 64160 receiving points (1 st shot 6282, receiving point 17030, 2 nd shot 7595, receiving point 15716, 3 rd shot 7564, receiving point 15726, 3 rd shot 7594, receiving point 15688). The subsurface imaging point has 500×60×1401 points, and its voxel size is 12.5×25×5m, and its depth is from 3km to 7km.

Step S42: an illumination dip angle gather is calculated.

In the embodiment of the invention, the (x, y) positions of the imaging areas are 500×60, namely 60 imaging lines, and each imaging line has 500 common center points (CDPs). That is, 30000 illumination tilt gathers were calculated in total. The input seismic data channel length is 7s, the sampling interval is 2ms, and the sampling point number is 3500. The offset apertures in the inline and crossline directions were 6000m and 4500m, respectively. Each of the output illumination inclination angle gathers has 21 channels, the illumination inclination angle is from-30 degrees to 30 degrees, and the illumination inclination angle increment is 3 degrees. The track length of each track in the illumination dip angle track set is 4km (from 3km to 7 km), the sampling interval is 5m, and the imaging point number is 801.

Referring to fig. 5, a schematic diagram of a calculated illumination tilt angle gather is shown.

Step S43: a generalized radon transform is computed.

In this embodiment, a preconditioned conjugate gradient algorithm (PCG) is used to calculate a high resolution three-dimensional radon modelThe maximum iteration number of the outer loop is 5 times, the iteration number of the inner loop is 10 times, and the allowable residual ratio parameter is 1%.

Step S44: and (5) superposition imaging.

In this embodiment, the three-dimensional radon model channels of p E < -0.2,0.2 > are horizontally overlapped to obtain a scattered wave imaging channel, and the three-dimensional radon model channels of p E < -7 > -3 > are horizontally overlapped to obtain a reflected wave imaging channel. Finally, according to 10: and 1, synthesizing the scattered wave and the reflected wave imaging channel to obtain a final imaging channel.

Step S45: imaging data volume synthesis.

And (3) calculating the illumination dip angle gathers at each (x, y) position in the steps S43 and S44 to obtain scattered wave imaging data bodies, reflected wave imaging data bodies and imaging data bodies synthesized by the scattered wave imaging data bodies and the reflected wave imaging data bodies.

In this embodiment, it is necessary to repeat the calculations S43 and S44 30000 times. The resulting imaged data volume has 60 imaged lines, each with 500 CDPs. Each imaging track has a track length of 4km (from 3km to 7 km), a sampling interval of 5m, and an imaging point number of 801.

Compared with the technology mentioned in the background art, the embodiment of the invention is characterized in that the scattered wave and the reflected wave are used as effective signals to image the underground medium, wherein the reflected wave is used for imaging a large-scale structure, the scattered wave is used for imaging a small-scale geologic body, and the important point is to protect the weak scattered wave, thereby realizing high-definition imaging. The embodiment of the invention has the advantages that scattered waves are changed into valuable, the imaging definition of a small-scale scatterer is improved, the signal to noise ratio and the resolution of a large-scale reflection interface are improved, and the lithologic oil and gas reservoir exploration is facilitated. FIG. 5 is a plot of illumination dip calculated in an embodiment of the invention, the plot depth ranging from 3km to 7km, illumination dip ranging from-30 to 30, depth sampling interval of 5m, and illumination dip sampling interval of 3. Fig. 5 contains multiple groups of concave-down homophase axes of reflected waves. The vertex of the same phase axis of the reflected wave with the depth of more than 4km is positioned at 0 degrees and corresponds to a horizontal stratum; the off-axis peak of the reflected wave near 5.3km in depth is at-6 deg., corresponding to the right declined formation. Carefully recognizing that there is a set of straight scattering wave in-phase axes around 5.9km in depth, corresponding to karst cave scatterers. Fig. 6 is a diagram of a conventional imaging trace (a), a reflected wave imaging trace (b), and a scattered wave imaging trace (c) corresponding to the trace set of fig. 5 in an embodiment of the present invention. Fig. 6a is an imaging trace obtained by horizontally superimposing the illumination inclination angle trace set of fig. 5, which is called a conventional imaging trace in which both reflected waves and scattered waves are imaged, but imaging noise is contained in reflected wave imaging. In fig. 6b is a reflected wave imaging trace, without imaging noise. In fig. 6 c is a scattered wave imaging trace, and the scattered wave imaging at a depth of 5.9km is clear. Comparing the three imaging channels in fig. 6, it can be found that the reflected wave and the scattered wave are imaged independently, and imaging noise caused by the imaging aperture is eliminated. FIG. 7 shows a conventional imaging profile calculated in this embodiment of the present invention, with common center point (CDP) numbers from 1502 to 2300, wherein circles indicate imaging of karst cave scatterers, and it can be seen that the profile contains a lot of noise, has low signal-to-noise ratio, poor continuity, and has complex waves below the non-integration surface with depth of 4.4km, which affects resolution. Fig. 8 is a reflected wave imaging profile calculated in an embodiment of the present invention. Comparing fig. 7 and 8, it is easy to see that the signal-to-noise ratio of the reflected wave imaging profile is obviously improved, imaging noise is not included, the continuity of the reflecting layer is obviously improved, the composite wave below the non-integrated surface with the depth of 4.4km is eliminated, and the resolution is improved. The top surface of the ancient buried hill, the multi-stage sedimentary body between the ancient buried hill and the non-integrated surface and the buried hill inner curtain fracture system are clearly imaged. FIG. 9 is a cross section of scattered wave imaging calculated in an embodiment of the invention that more clearly images subsurface small scale bodies, with clear imaging of both the top and bottom of the karst cave bodies shown by circles, without the appearance of scissors artifacts. FIG. 10 is an imaging profile synthesized in an embodiment of the invention, combining scattered wave imaging and scattered wave imaging at a 10:1 ratio. The invention can ensure the synthesizing effect just because the phase characteristics and the waveform characteristics of the reflected wave and the scattered wave are kept. As can be seen from fig. 10, both the large scale reflective interface and the small scale geologic volume are clearly imaged. This example demonstrates that the invention improves the imaging definition of small scale scatterers, the signal to noise ratio and resolution of large scale reflective interfaces.

Based on the inventive concept, the embodiment of the present invention further provides a seismic wave imaging apparatus, the structure of which is shown in fig. 11, including:

an illumination inclination calculation module 111, configured to calculate, for each scattering point at the imaging point position of the investigation region, an illumination inclination of the scattering point according to an incident ray parameter and a scattered ray parameter of the scattering point;

an illumination dip angle gather establishing module 112, configured to determine an imaging aperture common to the reflected wave and the scattered wave in the seismic wave according to the maximum and minimum illumination dip angles, and establish an illumination dip angle gather for each imaging lane of the investigation region according to the imaging aperture and the illumination dip angle of each scattering point on the imaging lane;

The radon transformation module 113 is configured to perform generalized radon transformation in three dimensions of depth, radon parameters and stratum inclination angle for each illumination inclination angle gather, so as to obtain a radon model with separated scattered wave and reflected wave;

An imaging channel imaging module 114, configured to image the scattered wave and the reflected wave in each radon model respectively, so as to obtain a scattered wave imaging channel and a reflected wave imaging channel;

an imaging data volume synthesis module 115 for synthesizing the scattered wave imaging data volume from the scattered wave imaging channels and synthesizing the reflected wave imaging data volume from the reflected wave imaging channels.

The specific manner in which the various modules perform the operations in the apparatus of the above embodiments have been described in detail in connection with the embodiments of the method, and will not be described in detail herein.

Based on the inventive concept, the embodiments of the present invention further provide a computer storage medium, in which computer executable instructions are stored, which when executed by a processor, implement the above-mentioned seismic wave imaging method.

Based on the inventive concept of the present invention, an embodiment of the present invention further provides a server, including: the system comprises a memory, a processor and a computer program stored in the memory and capable of running on the processor, wherein the processor realizes the seismic wave imaging method when executing the program.

Unless specifically stated otherwise, terms such as processing, computing, calculating, determining, displaying, or the like, may refer to an action and/or process of one or more processing or computing systems, or similar devices, that manipulates and transforms data represented as physical (e.g., electronic) quantities within the processing system's registers or memories into other data similarly represented as physical quantities within the processing system's memories, registers or other such information storage, transmission or display devices. Information and signals may be represented using any of a variety of different technologies and techniques. For example, data, instructions, commands, information, signals, bits, symbols, and chips that may be referenced throughout the above description may be represented by voltages, currents, electromagnetic waves, magnetic fields or particles, optical fields or particles, or any combination thereof.

It should be understood that the specific order or hierarchy of steps in the processes disclosed are examples of exemplary approaches. Based on design preferences, it is understood that the specific order or hierarchy of steps in the processes may be rearranged without departing from the scope of the present disclosure. The accompanying method claims present elements of the various steps in a sample order, and are not meant to be limited to the specific order or hierarchy presented.

In the foregoing detailed description, various features are grouped together in a single embodiment for the purpose of streamlining the disclosure. This method of disclosure is not to be interpreted as reflecting an intention that the claimed embodiments of the subject matter require more features than are expressly recited in each claim. Rather, as the following claims reflect, invention lies in less than all features of a single disclosed embodiment. Thus the following claims are hereby expressly incorporated into this detailed description, with each claim standing on its own as a separate preferred embodiment of this invention.

Those of skill would further appreciate that the various illustrative logical blocks, modules, circuits, and algorithm steps described in connection with the embodiments disclosed herein may be implemented as electronic hardware, computer software, or combinations of both. To clearly illustrate this interchangeability of hardware and software, various illustrative components, blocks, modules, circuits, and steps have been described above generally in terms of their functionality. Whether such functionality is implemented as hardware or software depends upon the particular application and design constraints imposed on the overall system. Skilled artisans may implement the described functionality in varying ways for each particular application, but such implementation decisions should not be interpreted as causing a departure from the scope of the present disclosure.

The steps of a method or algorithm described in connection with the embodiments disclosed herein may be embodied directly in hardware, in a software module executed by a processor, or in a combination of the two. A software module may reside in RAM memory, flash memory, ROM memory, EPROM memory, EEPROM memory, registers, hard disk, a removable disk, a CD-ROM, or any other form of storage medium known in the art. An exemplary storage medium is coupled to the processor such the processor can read information from, and write information to, the storage medium. In the alternative, the storage medium may be integral to the processor. The processor and the storage medium may reside in an ASIC. The ASIC may reside in a user terminal. The processor and the storage medium may reside as discrete components in a user terminal.

For a software implementation, the techniques described in this disclosure may be implemented with modules (e.g., procedures, functions, and so on) that perform the functions described herein. These software codes may be stored in memory units and executed by processors. The memory unit may be implemented within the processor or external to the processor, in which case it can be communicatively coupled to the processor via various means as is known in the art.

The foregoing description includes examples of one or more embodiments. It is, of course, not possible to describe every conceivable combination of components or methodologies for purposes of describing the aforementioned embodiments, but one of ordinary skill in the art may recognize that many further combinations and permutations of various embodiments are possible. Accordingly, the embodiments described herein are intended to embrace all such alterations, modifications and variations that fall within the scope of the appended claims. Furthermore, as used in the specification or claims, the term "comprising" is intended to be inclusive in a manner similar to the term "comprising," as interpreted when employed as a transitional word in a claim. Furthermore, any use of the term "or" in the specification of the claims is intended to mean "non-exclusive or".

Claims

1. A method of seismic imaging, comprising:

2. The method according to claim 1, wherein calculating the illumination tilt angle of the scattering point based on the incident ray parameter and the scattered ray parameter of the scattering point comprises:

according to the incidence angle and the declination angle of the scattered rays, the illumination inclination angle of the scattering point is calculated by the following formula:

In the above formula, θ is an illumination inclination angle, α _SM and β _SM are respectively a downinclination angle and an azimuth angle of an incident ray, α _GM and β _GM are respectively a downinclination angle and an azimuth angle of a scattered ray, and α is half of an included angle between the incident ray and the scattered ray:

3. the method as recited in claim 2, further comprising:

And determining the positive and negative of the illumination inclination angle of the scattering point according to the inclination direction of the stratum where the scattering point is located.

4. The method of claim 2, wherein the incident and scattered radiation parameters of the scattering point are determined by:

And according to a seismic wave velocity model of the research area, using a program function equation to respectively calculate ray paths from the shot point and the receiving point to each imaging point, and obtaining an incident ray parameter and a scattered ray parameter of a scattering point at the position of the imaging point.

5. The method of claim 1, wherein said creating an illumination dip gather based on the illumination dip of each scattering point on the imaging aperture and imaging tunnel, comprises:

taking a projection point of an imaging channel on the ground surface as a central point, and determining an integration surface according to the central point and the imaging aperture;

According to the integral surface of the imaging channel and the illumination dip angle of each scattering point on the imaging channel, an illumination dip angle channel set is established by the following formula:

In the above formula, d (z, theta) is an illumination inclination angle gather, z is a position variable in the vertical direction, and theta is a variable illumination inclination angle; hit (z, θ) is the number of illumination times; the points above D represent the derivation of time, t is the recording time of the seismic data, and r _S、r_G、τ_S and τ _G are the shot position, the receiving position, the travel time from shot to scattering and the travel time from receiving to scattering in sequence; s ₁ is an integration surface; ζ is the position of the point on the integration surface; for weighting operators,/> Wherein A _S and A _G are amplitudes from a shot point and a receiving point to a scattering point, and h (z, ζ) is a Beylkin determinant.

6. The method of claim 5, wherein prior to establishing the illumination dip gather, further comprising:

Calculating a central ray by using a kinematic ray tracing method aiming at each scattering point to obtain tau _S and tau _G when the shot point and the receiving point travel to the scattering point;

And calculating paraxial rays by using a dynamic ray tracing method to obtain a geometric diffusion matrix, obtaining amplitudes A _S and A _G from shot points and receiving points to scattering points, obtaining a second partial derivative of coordinates of the central rays when traveling, and further obtaining a Beylkin determinant h (z, ζ).

7. The method of claim 1, wherein for each illumination dip gather, the generalized radon transformation in three dimensions of depth, radon parameters and formation dip is performed to obtain a radon model with scattered and reflected wave separation, specifically comprising:

aiming at each illumination dip angle gather, the radon model with the separation of scattered waves and reflected waves in three dimensions of depth, radon parameters and formation dip angle is obtained by solving the following formula:

In the above-mentioned method, the step of, For the radon model, χ is depth, p is radon parameter,/>Is the formation dip angle; d (z, θ) is the illumination tilt angle gather, z is the position variable in the vertical direction, θ is the variable illumination tilt angle; w is the weight coefficient of the Ladong model; /(I)Fourier inverse transform operator representing depth,/>A fourier transform operator representing depth; l is defined by the kernel function/>Constructed generalized radon transform operator,/>Is the complex conjugate transpose of L.

8. The method of claim 7, wherein solving the following formula specifically comprises:

the following formula is solved using an iterative algorithm.

9. The method of claim 8, wherein said solving the following formula using an iterative algorithm comprises:

The following formula is solved using a preconditioned conjugate gradient algorithm PCG.

10. The method according to claim 1, wherein the radon model for obtaining the separation of the scattered wave and the reflected wave comprises:

and obtaining a radon model with scattered waves distributed in a region with a radon parameter of 0 and reflected waves distributed in a region with a radon parameter of less than 0.

11. The method of any one of claims 1 to 10, wherein the imaging step further comprises:

the scattered wave imaging channel and the reflected wave imaging channel are weighted and overlapped to obtain a synthetic imaging channel;

correspondingly, the imaging data volume synthesis step further comprises the following steps:

A synthetic imaging data volume is obtained from the synthetic imaging modality.

12. A seismic imaging apparatus, comprising:

13. A computer storage medium having stored therein computer executable instructions which when executed by a processor implement the seismic wave imaging method of any of claims 1 to 11.