CN115469362B - Energy flow density vector calculation method in seismic exploration - Google Patents
Energy flow density vector calculation method in seismic exploration Download PDFInfo
- Publication number
- CN115469362B CN115469362B CN202211122422.7A CN202211122422A CN115469362B CN 115469362 B CN115469362 B CN 115469362B CN 202211122422 A CN202211122422 A CN 202211122422A CN 115469362 B CN115469362 B CN 115469362B
- Authority
- CN
- China
- Prior art keywords
- field
- wave
- wavefield
- wave field
- point
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
- G01V1/30—Analysis
- G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
Abstract
The application relates to the field of seismic wave oil and gas exploration, and discloses an energy flow density vector calculation method in seismic exploration, which comprises the following steps of: forward continuation is carried out on the wave field based on the wave equation, and a displacement field of the seismic source wave field is obtained; gradient operation is carried out on the displacement field of the seismic source wave field stored in the step one to obtain a stress field of the seismic source wave field, and derivative operation of propagation time is carried out on the displacement field of the seismic source wave field stored in the step one to obtain a velocity field of the seismic source wave field; scanning the velocity field of the seismic source wave field in the second step on the propagation time, and storing the stress field and the velocity field of each space point when the amplitude mode value of the velocity field is maximum; obtaining energy flow density vectors of the wave field according to an energy flow density calculation method; the application scans wave field amplitude values of each space point underground, then obtains energy flow density vectors with higher accuracy, and can obtain high-precision angle domain common imaging point gathers based on the energy flow density vectors.
Description
Technical Field
The application relates to the field of seismic wave oil and gas exploration, in particular to an energy flow density vector calculation method in seismic exploration.
Background
Prestack depth migration imaging is a seismic data processing technique used to image subsurface geologic structures, aid in hydrocarbon exploration and production, and although migration imaging is one of the most widely used techniques in the industry today, one of the important research directions in the industry today is the construction of subsurface velocity fields because it is limited by the accuracy of the input velocity fields. With the continuous development of the seismic exploration technology, the velocity modeling method is gradually changed from a conventional interactive NMO velocity analysis method to an offset velocity analysis method, and among a plurality of velocity modeling methods, a tomography method based on angle domain common imaging point gathers (angle gathers) is widely applied. The angle gather calculation method based on the direction vector method is gradually applied in industry due to simple theory and easy realization.
There are various direction vector methods for calculating the propagation angle of the seismic wave and thus the angle gather, and among them, calculation methods based on energy density vectors (poyning vectors) (Yoon k and Marfurt k, 2006, reverse-time migration using the Poynting vector, expression geodesics, 37 (1): 102-107) are theoretically clear and gradually applied. The inventors of the present application found that: because the seismic waves are complex in underground propagation, the problem of wave field aliasing at different moments of imaging points in the same space is easy to occur, so that the energy density calculation is not accurate enough, and the final speed modeling precision is affected.
Disclosure of Invention
The application aims to provide an energy flow density vector calculation method in seismic exploration, which is used for scanning wave field amplitude values of each space point underground and then obtaining an energy flow density vector with higher accuracy, so that a high-accuracy angle domain common imaging point gather can be obtained based on the energy flow density vector.
In order to achieve the above object, the present application provides a method for calculating energy flow density vectors in seismic exploration, comprising the steps of:
step one, forward continuation is carried out on wave fields based on wave equation, and displacement fields of the source wave fields are obtained;
step two, gradient operation is carried out on the displacement field of the seismic source wave field stored in the step one to obtain a stress field of the seismic source wave field, and derivative operation of propagation time is carried out on the displacement field of the seismic source wave field stored in the step one to obtain a velocity field of the seismic source wave field;
scanning the velocity field of the seismic source wave field in the second step on the propagation time, and storing the stress field and the velocity field of each space point when the amplitude modulus of the velocity field is maximum;
and step four, obtaining the energy flow density vector of the wave field according to the energy flow density calculation method.
As a preferred embodiment of the present application, the following verification step is further included,
step five, calculating the incidence angle of the wave field according to the wave field propagation direction calculation formula by using the stratum dip angle and the energy flow density vector of the seismic source wave field obtained in the step four;
step six, reverse continuation is carried out on the wave field based on the wave equation, and a displacement field of the detected wave field is obtained;
and step seven, calculating an angle domain common imaging point gather by using the seismic source wave field and detection wave field displacement field, the incident angle of the wave field and an angle domain common imaging point gather calculation formula, and comparing the incident angle in the theoretical incident angle and the angle domain common imaging point gather so as to verify the calculation accuracy of the energy flow density vector calculated in the step four.
In the first step, the velocity field of the longitudinal wave medium is based on a second-order acoustic wave equation to obtain displacement fields with different forward extension times and store the displacement fields P of the source wave field with different propagation times of each spatial point S (t,x)，
Calculating the displacement field P of the source wave field by adopting the following second-order acoustic wave equation S (t,x)：
Wherein x=x, z represents the coordinate of a space point in a cartesian rectangular coordinate system, x is the coordinate in the horizontal direction, and z is the coordinate in the vertical direction; t is the travel time of the wavefield; p is the displacement field of the wave field, and S is the source wave field; v (V) P (x) Is the longitudinal wave velocity at the spatial point x.
In a preferred embodiment of the present application, in the second step,
the stress field of the seismic source wave field is calculated by adopting a displacement field gradient operation formula of the seismic source wave field:
in the middle of For the derivative of the displacement field at the spatial point x in the x-direction, +.>Is the derivative of the displacement field at the spatial point x in the z-direction; />Is the source wavefield stress field at spatial point x; x=x, z represents the coordinate of a space point in a cartesian rectangular coordinate system, x is the coordinate in the horizontal direction, and z is the coordinate in the vertical direction; p is the displacement field of the wave field, and S is the source wave field; t is the travel time of the wavefield; />Is a gradient operator.
The derivative of the displacement field of the source wavefield with respect to the travel time t of the wavefield is calculated using the formula:
in the method, in the process of the application,the derivative of the travel time wave field of the seismic source wave field when the displacement field mode value is maximum, P is the displacement field of the wave field, and the superscript S is the seismic source wave field; t is the travel time of the wave field, δt is the sampling interval of the travel time of the wave field;
in the third step, as a preferable mode of the present application,
calculating the maximum amplitude module value of each underground space point on all propagation time by adopting the following displacement field maximum amplitude module value formula, and simultaneously storing the propagation time of the maximum amplitude module value:
in the method, in the process of the application,is the maximum value of the displacement field mode at the space point x; x=x, z represents the coordinate of a space point in a cartesian rectangular coordinate system, x is the coordinate in the horizontal direction, and z is the coordinate in the vertical direction; i represents modulo operation, max represents the maximum calculation, T max For maximum travel time, Δt is the sampling interval of the travel time of the wavefield.
In a preferred embodiment of the present application, in the fourth step,
the energy flow density vector in the seismic wavefield is calculated using the formula:
where PV represents the energy flow density vector in the seismic wavefield, v represents the particle velocity of vibration of the wavefield propagation, τ represents the stress;
the energy flow density vector in the source wavefield is calculated using the following formula:
in the method, in the process of the application,for the component of the energy flow density vector of the source wavefield in the x-direction at the spatial point x, +.>A component along the z-direction of the fluence density vector of the source wavefield at the spatial point x; />The derivative in the x direction of the maximum value of the displacement field mode at the spatial point x; />The derivative in the z direction of the maximum value of the displacement field mode at the spatial point x; />A derivative of the maximum value of the displacement field mode at the spatial point x with respect to the propagation time t; x=x, z denotes cartesian rectangular coordinatesThe coordinate of a space point in the system, x is the coordinate in the horizontal direction, and z is the coordinate in the vertical direction; p is the displacement field of the wave field, and S is the source wave field; t is the travel time of the wavefield.
In a preferred embodiment of the present application, in the fifth step,
the angle of incidence of the wavefield to the achieved image point is calculated using the following wavefield propagation direction calculation formula:
where α (x) represents the angle of incidence of the wavefield to the point x of the achieved image;for the component of the energy flow density vector of the source wavefield in the x-direction at the spatial point x, +.>A component along the z-direction of the fluence density vector of the source wavefield at the spatial point x; θ (x) is the formation dip.
In the sixth aspect of the present application, in the step,
taking the seismic data recorded on the ground as a boundary value condition, performing reverse time prolongation on a longitudinal wave medium velocity field based on a second-order acoustic wave equation, and storing displacement fields P of detection wave fields with different propagation times of each space point R (t,x)，
Determination of the displacement field P of a detected wavefield using the second-order acoustic wave equation R (t,x)：
Wherein x=x, z represents the coordinate of a space point in a cartesian rectangular coordinate system, x is the coordinate in the horizontal direction, and z is the coordinate in the vertical direction; t is the travel time of the wavefield; p is the displacement field of the wave field, and the superscript R is the detected wave field; v (V) P (x) To be located at a space pointLongitudinal wave velocity at x.
In the preferred embodiment of the present application, in the seventh step,
calculating an angle domain common imaging point gather by adopting the following angle domain common imaging point gather calculation formula:
wherein I is ADCIGs Representing the angle domain co-imaging point gather, x represents the coordinate of a space point in a Cartesian coordinate system, T is the propagation time of the wave field, T max Is the maximum propagation time; θ (x) is the angle of incidence of the wave field to the point where the image point is x, θ k For angles of incidence in the angular domain common imaging point gather, δ represents a δ function; p is the displacement field of the wavefield, the superscript S represents the source wavefield, the superscript R represents the detected wavefield, and x represents the convolution operator.
In the preferred embodiment of the present application, in the seventh step,
and selecting an angle domain common imaging point gather at the horizontal position of 1.75km, calculating a theoretical incident angle at the horizontal position of 1.75km through Snell's law and a triangular relation, extracting the angle domain common imaging point gather at the horizontal position of 1.75km through an angle domain common imaging point gather calculation formula, and comparing the theoretical incident angle with the incident angle in the angle domain common imaging point gather obtained through the angle domain common imaging point gather calculation formula so as to verify the accuracy of energy flow density vector calculation in the fourth step.
Compared with the prior art, the energy flow density vector calculation method in the seismic exploration has the beneficial effects that:
according to the application, the wave field amplitude value of each space point underground is scanned, then the energy flow density vector with higher accuracy is obtained, and the high-accuracy angle domain common imaging point gather can be obtained on the basis of the energy flow density vector with high accuracy.
Drawings
In order to more clearly illustrate the technical solution of the embodiments of the present application, the drawings of the embodiments will be briefly described below.
FIG. 1 is a flow chart of steps provided by an embodiment of the present application;
FIG. 2 is a schematic diagram of an embodiment of the present application providing a longitudinal wave medium velocity field model constructed from horizontal and sloped interfaces;
FIG. 3 is a schematic illustration of the dip angle of the reflective interface formation obtained by scanning the longitudinal wave medium velocity field of FIG. 1;
FIG. 4 is a schematic diagram of a source wavefield displacement field when forward extending to 750 ms;
FIG. 5 is a schematic diagram of a velocity field with maximum source wavefield amplitude value when forward extended to 750 ms;
FIG. 6a is a schematic diagram of the component of the fluence vector in the x-direction when forward extended to 750 ms;
FIG. 6b is a schematic diagram of the component of the fluence vector in the z-direction when forward extension is 750 ms;
FIG. 7 is a schematic representation of the imaging angle for each spatial point obtained from the calculated fluence vectors and formation dip;
FIG. 8 is a schematic diagram of the detected wavefield displacement field when the inverse extension is 750 ms;
FIG. 9a is a schematic diagram of a longitudinal wave medium velocity field labeled shot, detector, ray path, and theoretical calculated wave field incident angle;
FIG. 9b is a schematic diagram of an angle gather labeled with the migration calculated wave field incident angle and the theoretical calculated wave field incident angle.
Detailed Description
The following describes in further detail the embodiments of the present application with reference to the drawings and examples. The following examples are illustrative of the application and are not intended to limit the scope of the application.
In the description of the present application, it should be understood that the directions or positional relationships indicated by the terms "upper", "lower", "left", "right", "top", "bottom", etc. are based on the directions or positional relationships shown in the drawings, are merely for convenience of describing the present application and simplifying the description, and do not indicate or imply that the devices or elements referred to must have a specific orientation, be constructed and operated in a specific orientation, and thus should not be construed as limiting the present application. It should be understood that the terms "first," "second," and the like are used herein to describe various information, but such information should not be limited to these terms, which are used merely to distinguish one type of information from another. For example, a "first" message may also be referred to as a "second" message, and similarly, a "second" message may also be referred to as a "first" message, without departing from the scope of the application.
As shown in fig. 1 to 9b, a method for calculating energy flow density vectors in seismic exploration according to a preferred embodiment of the present application includes the steps of:
step one, forward continuation is carried out on wave fields based on wave equation, and a displacement field of a seismic source wave field is obtained (S1); in the present embodiment, taking the longitudinal wave medium velocity field model shown in fig. 2 as an example, based on the second-order acoustic wave equation, the displacement fields of different forward extension times are obtained and the displacement field P of the source wave field of different propagation times at each spatial point is stored S (t,x)，
Calculating the displacement field P of the source wave field by adopting the following second-order acoustic wave equation S (t,x)：
Wherein x=x, z represents the coordinate of a space point in a cartesian rectangular coordinate system, x is the coordinate in the horizontal direction, and z is the coordinate in the vertical direction; t is the travel time of the wavefield; p is the displacement field of the wave field, and S is the source wave field; v (V) P (x) Is the longitudinal wave velocity at the spatial point x. The displacement field of the seismic source wave field shown in fig. 4 can be obtained based on the second-order acoustic wave equation at the moment of forward extension to 750 ms;
step two, gradient operation is carried out on the displacement field of the seismic source wave field stored in the step one to obtain a stress field of the seismic source wave field, and derivative operation of propagation time is carried out on the displacement field of the seismic source wave field stored in the step one to obtain a velocity field of the seismic source wave field (S2);
specifically, the stress field of the seismic source wave field is calculated by adopting a displacement field gradient operation formula of the seismic source wave field:
in the method, in the process of the application, for the derivative of the displacement field at the spatial point x in the x-direction, +.>Is the derivative of the displacement field at the spatial point x in the z-direction; />Is the source wavefield stress field at spatial point x; x=x, z represents the coordinate of a space point in a cartesian rectangular coordinate system, x is the coordinate in the horizontal direction, and z is the coordinate in the vertical direction; p is the displacement field of the wave field, and S is the source wave field; t is the travel time of the wavefield; />Is a gradient operator;
the derivative of the displacement field of the source wavefield with respect to the travel time t of the wavefield is calculated using the formula:
in the method, in the process of the application,the derivative of the travel time of the wave field when the displacement field mode value is maximum, P is the displacement field of the wave field, and the superscript S is the source wave field; t is the travel time of the wavefield.
Scanning the velocity field of the seismic source wave field in the second step on the propagation time, and storing the stress field and the velocity field of each space point when the amplitude modulus of the velocity field is maximum (S3);
specifically, the maximum amplitude modulus of each underground space point on all propagation times is calculated by adopting the following displacement field maximum amplitude modulus formula, and the propagation time t with the maximum amplitude modulus is saved max ：
In the method, in the process of the application,is the maximum value of the displacement field mode at the space point x; x=x, z represents the coordinate of a space point in a cartesian rectangular coordinate system, x is the coordinate in the horizontal direction, and z is the coordinate in the vertical direction; i represents modulo operation, max represents maximum operation, T max For the maximum travel time, t is the travel time of the wavefield and Δt is the sampling interval of the travel time of the wavefield. In this embodiment, this enables to obtain a point of underground space when the amplitude mode maximum occurs when forward continuation reaches 750ms (as shown in fig. 5).
Step four, obtaining energy flow density vectors of wave fields according to an energy flow density calculation method;
specifically, the energy flow density vector in the seismic wavefield is calculated using the following formula:
where PV represents the energy flow density vector in the seismic wavefield, v represents the particle velocity of vibration of the wavefield propagation, τ represents the stress;
the energy flow density vector in the source wavefield is calculated using the following formula:
in the method, in the process of the application,for the component of the energy flow density vector of the source wavefield in the x-direction at the spatial point x, +.>A component along the z-direction of the fluence density vector of the source wavefield at the spatial point x; />The derivative in the x direction of the maximum value of the displacement field mode at the spatial point x; />The derivative in the z direction of the maximum value of the displacement field mode at the spatial point x; />A derivative of the maximum value of the displacement field mode at the spatial point x with respect to the propagation time t; x=x, z represents the coordinate of a space point in a cartesian rectangular coordinate system, x is the coordinate in the horizontal direction, and z is the coordinate in the vertical direction; p is the displacement field of the wave field, and S is the source wave field; t is the travel time of the wavefield. In the present embodiment, it is thus possible to obtain a component of the energy density vector in the x direction at a spatial point when the mode maximum occurs when the forward extension is 750ms (as shown in fig. 6 a), and a component of the energy density vector in the z direction at a spatial point when the mode maximum occurs when the forward extension is 750ms (as shown in fig. 6 b).
Illustratively, the energy-flow density vector computation in the seismic survey further includes the step of verifying that the energy-flow density vector computation of the wavefield is accurate,
step five, calculating the incidence angle of the wave field according to a wave field propagation direction calculation formula by using the known stratum dip angle (shown in figure 3) and the energy flow density vector of the seismic source wave field obtained in step four (S4); specifically, in the fifth step,
the angle of incidence of the wavefield to the achieved image point is calculated using the following wavefield propagation direction calculation formula:
where α (x) represents the angle of incidence of the wavefield to the point x of the achieved image;for the component of the energy flow density vector of the source wavefield in the x-direction at the spatial point x, +.>A component along the z-direction of the fluence density vector of the source wavefield at the spatial point x; θ (x) is the formation dip. This enables the angle of incidence at each imaging point x over the whole space to be obtained (as shown in fig. 7).
Step six, reverse continuation is carried out on the wave field based on the wave equation, and a displacement field of the detected wave field is obtained (S5); in this embodiment, taking the longitudinal wave medium velocity field model shown in fig. 2 as an example, using the seismic data recorded on the ground as a boundary value condition, performing reverse time extension based on a second-order acoustic wave equation, and storing the displacement field P of the detected wave field with different propagation times at each spatial point R (t,x)，
Determination of the displacement field P of a detected wavefield using the second-order acoustic wave equation R (t,x)：
Wherein x=x, z represents the coordinate of a space point in a cartesian rectangular coordinate system, x is the coordinate in the horizontal direction, and z is the coordinate in the vertical direction; t is the travel time of the wavefield; p is the displacement field of the wave field, and the superscript R is the detected wave field; v (V) P (x) Is the longitudinal wave velocity at the spatial point x. This can obtain a second-order acoustic wave equation based on the reverse continuation to 750msTo the displacement field of the detected wavefield as shown in fig. 8.
And step seven, calculating an angle domain common imaging point gather by using the seismic source wave field and detection wave field displacement field, the incident angle of the wave field and an angle domain common imaging point gather calculation formula, and comparing the theoretical incident angle with the incident angle in the angle domain common imaging point gather to verify the calculation accuracy of the energy flow density vector calculated in the step four (S6). Specifically, the following angle domain common imaging point gather calculation formula is adopted to calculate the angle domain common imaging point gather:
wherein I is ADCIGs Representing the angle domain co-imaging point gather, x represents the coordinate of a space point in a Cartesian coordinate system, T is the propagation time of the wave field, T max Is the maximum propagation time; θ (x) is the angle of incidence of the wave field to the point where the image point is x, θ k For angles of incidence in the angular domain common imaging point gather, δ represents a δ function; p is the displacement field of the wavefield, the superscript S represents the source wavefield, the superscript R represents the detected wavefield, and x represents the convolution operator.
Illustratively, in the seventh step, an angle domain co-imaging point gather at a horizontal position of 1.75km is selected, a theoretical incident angle at the horizontal position of 1.75km is calculated by Snell's law and trigonometric relation (as shown in fig. 9 a), the angle domain co-imaging point gather at the horizontal position of 1.75km is extracted by an angle domain co-imaging point gather calculation formula, and the theoretical incident angle is compared with the incident angle in the obtained angle domain co-imaging point gather in the angle domain co-imaging point gather calculation formula to verify the accuracy of the energy flow density vector calculation in the fourth step. And when the theoretical incidence angle is the same as the incidence angle in the calculated angle domain common imaging point trace set, the energy flow density vector is accurately calculated. In the seventh step, the angle domain common imaging point gather at any position of the horizontal position can be selected according to the actual situation, and then the subsequent test is performed.
According to the application, the wave field amplitude value of each space point underground is scanned, then the energy flow density vector with higher accuracy is obtained, and the high-accuracy angle domain common imaging point gather can be obtained on the basis of the energy flow density vector with high accuracy.
The application also provides an energy flow density vector calculation device in seismic exploration, and a seismic source wave field displacement field acquisition module which is used for forward continuation of the wave field based on a wave equation to acquire a displacement field of the seismic source wave field;
the stress field and velocity field acquisition module is used for carrying out gradient operation and derivative operation on propagation time on the stored displacement field to acquire a stress field and a velocity field of a seismic source wave field respectively;
the system comprises a space point maximum amplitude mode value stress field and a speed field acquisition module, wherein the space point maximum amplitude mode value stress field and speed field acquisition module is used for scanning the seismic source wave field speed field in the propagation time and storing the stress field and the speed field when the amplitude mode value of each space point speed field is maximum;
the energy flow density vector acquisition module is used for acquiring energy flow density vectors of the wave field according to an energy flow density calculation method;
the wave field incidence angle acquisition module is used for calculating the incidence angle of the wave field according to a wave field propagation direction calculation formula by using the stratum dip angle and the energy flow density vector;
the wave field detection displacement field acquisition module is used for carrying out reverse continuation on the wave field based on the wave equation to obtain a displacement field of the wave field detection;
and the angle domain common imaging point gather acquisition module is used for calculating an angle domain common imaging point gather by utilizing displacement fields of the seismic source wave field and the detection wave field, the incidence angle of the wave field and an angle domain common imaging point gather calculation formula, and comparing the theoretical incidence angle with the incidence angle in the angle domain common imaging point gather.
In addition, the application also provides energy flow density vector computing equipment in seismic exploration, which comprises a processor and a memory, wherein the memory is used for storing at least one executable instruction, and the executable instruction enables the processor to execute the operation corresponding to the energy flow density vector computing method, and the specific operation steps are as follows:
forward continuation is carried out on the wave field based on the wave equation, and a displacement field of the seismic source wave field is obtained;
performing gradient operation and derivative operation on the travel time with the saved displacement field to obtain a stress field and a velocity field of a seismic source wave field respectively;
scanning the seismic source wave field velocity field in the propagation time, and storing a stress field and a velocity field when the amplitude modulus of each spatial point velocity field is maximum;
obtaining energy flow density vectors of the wave field according to an energy flow density calculation method;
calculating the incidence angle of the wave field according to a wave field propagation direction calculation formula by using the stratum dip angle and the energy flow density vector;
carrying out reverse continuation on the wave field based on the wave equation to obtain a displacement field of the detected wave field;
and calculating an angle domain common imaging point gather by using the seismic source wave field and the detection wave field displacement field, the incidence angle of the wave field and an angle domain common imaging point gather calculation formula, and comparing the angle domain common imaging point gather with a theoretical value incidence angle.
In the description of the present application, it should be noted that, unless explicitly specified and limited otherwise, the terms "connected," "connected," and "connected" are to be construed broadly, and may be either fixedly connected, detachably connected, or integrally connected, for example; can be mechanically or electrically connected; can be directly connected or indirectly connected through an intermediate medium, and can be communication between two elements. The specific meaning of the above terms in the present application will be understood in specific cases by those of ordinary skill in the art.
The foregoing is merely a preferred embodiment of the present application, and it should be noted that modifications and substitutions can be made by those skilled in the art without departing from the technical principles of the present application, and these modifications and substitutions should also be considered as being within the scope of the present application.
Claims (10)
1. The energy flow density vector calculation method in the seismic exploration is characterized by comprising the following steps of:
step one, forward continuation is carried out on wave fields based on wave equation, and displacement fields of the source wave fields are obtained;
step two, gradient operation is carried out on the displacement field of the seismic source wave field stored in the step one to obtain a stress field of the seismic source wave field, and derivative operation of propagation time is carried out on the displacement field of the seismic source wave field stored in the step one to obtain a velocity field of the seismic source wave field;
scanning the velocity field of the seismic source wave field in the second step on the propagation time, and storing the stress field and the velocity field of each space point when the amplitude modulus of the velocity field is maximum;
and step four, obtaining the energy flow density vector of the wave field according to the energy flow density calculation method.
2. The method of energy density vector computation in seismic exploration according to claim 1, further comprising the step of validating,
step five, calculating the incidence angle of the wave field according to the wave field propagation direction calculation formula by using the stratum dip angle and the energy flow density vector of the seismic source wave field obtained in the step four;
step six, reverse continuation is carried out on the wave field based on the wave equation, and a displacement field of the detected wave field is obtained;
and step seven, calculating an angle domain common imaging point gather by using the seismic source wave field and detection wave field displacement field, the incident angle of the wave field and an angle domain common imaging point gather calculation formula, and comparing the incident angle in the theoretical incident angle and the angle domain common imaging point gather so as to verify the calculation accuracy of the energy flow density vector calculated in the step four.
3. The method of energy density vector calculation in seismic exploration according to claim 2, wherein in said step one, the longitudinal wave medium velocity field is based on a second order acoustic wave equation to obtain displacement fields of different forward extension times and preserve the displacement field P of the source wavefield of different propagation times for each spatial point S (t,x)，
Calculating the displacement field P of the source wave field by adopting the following second-order acoustic wave equation S (t,x)：
Wherein x=x, z represents the coordinate of a space point in a cartesian rectangular coordinate system, x is the coordinate in the horizontal direction, and z is the coordinate in the vertical direction; t is the travel time of the wavefield; p is the displacement field of the wave field, and S is the source wave field; v (V) P (x) Is the longitudinal wave velocity at the spatial point x.
4. The method for energy density vector calculation in seismic exploration according to claim 3, wherein in said step two,
the stress field of the seismic source wave field is calculated by adopting a displacement field gradient operation formula of the seismic source wave field:
in the middle ofFor the derivative of the displacement field at the spatial point x in the x-direction, +.>Is the derivative of the displacement field at the spatial point x in the z-direction; />Is the source wavefield stress field at spatial point x; x=x, z represents the coordinate of a space point in a cartesian rectangular coordinate system, x is the coordinate in the horizontal direction, and z is the coordinate in the vertical direction; p is the displacement field of the wave field, and S is the source wave field; t is the travel time of the wavefield; />For gradient calculationA seed;
the derivative of the displacement field of the source wavefield with respect to the travel time t of the wavefield is calculated using the formula:
in the method, in the process of the application,the derivative of the travel time of the wave field when the displacement field mode value is maximum, P is the displacement field of the wave field, and the superscript S is the source wave field; t is the travel time of the wavefield.
5. The method of energy density vector computation in seismic exploration according to claim 4, wherein in said step three,
calculating the maximum amplitude modulus of each underground space point on all propagation times by adopting the following displacement field maximum amplitude modulus formula, and simultaneously storing the propagation time t with the maximum amplitude modulus max ：
In the method, in the process of the application,is the maximum value of the displacement field mode at the space point x; x=x, z represents the coordinate of a space point in a cartesian rectangular coordinate system, x is the coordinate in the horizontal direction, and z is the coordinate in the vertical direction; i represents modulo operation, max represents maximum operation, T max For maximum travel time, Δt is the sampling interval of the travel time of the wavefield.
6. The method of energy density vector computation in seismic exploration according to claim 5, wherein in said step four,
the energy flow density vector in the seismic wavefield is calculated using the formula:
where PV represents the energy flow density vector in the seismic wavefield, v represents the particle velocity of vibration of the wavefield propagation, τ represents the stress;
the energy flow density vector in the source wavefield is calculated using the following formula:
in the method, in the process of the application,for the component of the energy flow density vector of the source wavefield in the x-direction at the spatial point x, +.>A component along the z-direction of the fluence density vector of the source wavefield at the spatial point x; />The derivative in the x direction of the maximum value of the displacement field mode at the spatial point x; />The derivative in the z direction of the maximum value of the displacement field mode at the spatial point x; />A derivative of the maximum value of the displacement field mode at the spatial point x with respect to the propagation time t; x=x, z represents the coordinate of a space point in a cartesian rectangular coordinate system, x is the coordinate in the horizontal direction, and z is the coordinate in the vertical direction; p is the displacement field of the wave field, and S is the source wave field; t is the propagation of the wave fieldTime.
7. The method of energy density vector computation in seismic exploration according to claim 6, wherein in said step five,
the angle of incidence of the wavefield to the achieved image point is calculated using the following wavefield propagation direction calculation formula:
where α (x) represents the angle of incidence of the wavefield to the point x of the achieved image;for the component of the energy flow density vector of the source wavefield in the x-direction at the spatial point x, +.>A component along the z-direction of the fluence density vector of the source wavefield at the spatial point x; θ (x) is the formation dip.
8. The method of energy density vector computation in seismic exploration according to claim 7, wherein in said step six,
taking the seismic data recorded on the ground as a boundary value condition, performing reverse time prolongation on a longitudinal wave medium velocity field based on a second-order acoustic wave equation, and storing displacement fields P of detection wave fields with different propagation times of each space point R (t,x)，
Determination of the displacement field P of a detected wavefield using the second-order acoustic wave equation R (t,x)：
Wherein x=x, z represents the coordinate of a space point in a Cartesian coordinate system, x is the coordinate in the horizontal direction, and z is the coordinate in the vertical directionMarking; t is the travel time of the wavefield; p is the displacement field of the wave field, and the superscript R is the detected wave field; v (V) P (x) Is the longitudinal wave velocity at the spatial point x.
9. The method of energy density vector computation in seismic exploration according to claim 8, wherein in said step seven,
calculating an angle domain common imaging point gather by adopting the following angle domain common imaging point gather calculation formula:
wherein I is ADCIGs Representing the angle domain co-imaging point gather, x represents the coordinate of a space point in a Cartesian coordinate system, T is the propagation time of the wave field, T max Is the maximum propagation time; θ (x) is the angle of incidence of the wave field to the point where the image point is x, θ k For angles of incidence in the angular domain common imaging point gather, δ represents a δ function; p is the displacement field of the wavefield, the superscript S represents the source wavefield, the superscript R represents the detected wavefield, and x represents the convolution operator.
10. The method of energy density vector computation in seismic exploration according to claim 9, wherein in said step seven,
and selecting an angle domain common imaging point gather at the horizontal position of 1.75km, calculating a theoretical incident angle at the horizontal position of 1.75km through Snel l's law and a triangular relation, extracting the angle domain common imaging point gather at the horizontal position of 1.75km through an angle domain common imaging point gather calculation formula, and comparing the theoretical incident angle with the incident angle in the angle domain common imaging point gather obtained through the angle domain common imaging point gather calculation formula so as to verify the accuracy of energy flow density vector calculation in the fourth step.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211122422.7A CN115469362B (en) | 2022-09-15 | 2022-09-15 | Energy flow density vector calculation method in seismic exploration |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211122422.7A CN115469362B (en) | 2022-09-15 | 2022-09-15 | Energy flow density vector calculation method in seismic exploration |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115469362A CN115469362A (en) | 2022-12-13 |
| CN115469362B true CN115469362B (en) | 2023-10-10 |
Family
ID=84371313
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211122422.7A Active CN115469362B (en) | 2022-09-15 | 2022-09-15 | Energy flow density vector calculation method in seismic exploration |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115469362B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116381783B (en) * | 2023-03-30 | 2023-10-03 | 东北石油大学 | Stable and efficient Potentilla vector extraction method and system |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102636809A (en) * | 2012-03-27 | 2012-08-15 | 中国科学院地质与地球物理研究所 | Method for generating spreading angle domain common image point gathers |
| CN110618459A (en) * | 2019-11-01 | 2019-12-27 | 中国科学院地质与地球物理研究所 | Seismic data processing method and device |
| CN110945383A (en) * | 2017-04-06 | 2020-03-31 | 沙特阿拉伯石油公司 | Generating co-imaging gathers using wave field separation |
| CN112327359A (en) * | 2020-10-14 | 2021-02-05 | 山东省科学院海洋仪器仪表研究所 | Elastic reverse time migration method based on imaging energy flow vector |
| CN112462427A (en) * | 2020-11-13 | 2021-03-09 | 中国石油大学(华东) | Multi-component seismic data amplitude-preserving angle domain common imaging point gather extraction method and system |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8457899B2 (en) * | 2007-12-14 | 2013-06-04 | Shell Oil Company | Method of processing data obtained from seismic prospecting |
| US8619498B2 (en) * | 2010-09-24 | 2013-12-31 | CGGVeritas Services (U.S.) Inc. | Device and method for calculating 3D angle gathers from reverse time migration |
-
2022
- 2022-09-15 CN CN202211122422.7A patent/CN115469362B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102636809A (en) * | 2012-03-27 | 2012-08-15 | 中国科学院地质与地球物理研究所 | Method for generating spreading angle domain common image point gathers |
| CN110945383A (en) * | 2017-04-06 | 2020-03-31 | 沙特阿拉伯石油公司 | Generating co-imaging gathers using wave field separation |
| CN110618459A (en) * | 2019-11-01 | 2019-12-27 | 中国科学院地质与地球物理研究所 | Seismic data processing method and device |
| CN112327359A (en) * | 2020-10-14 | 2021-02-05 | 山东省科学院海洋仪器仪表研究所 | Elastic reverse time migration method based on imaging energy flow vector |
| CN112462427A (en) * | 2020-11-13 | 2021-03-09 | 中国石油大学(华东) | Multi-component seismic data amplitude-preserving angle domain common imaging point gather extraction method and system |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115469362A (en) | 2022-12-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Wu et al. | Directional illumination analysis using beamlet decomposition and propagation | |
| US5999488A (en) | Method and apparatus for migration by finite differences | |
| EP2409176B1 (en) | Seismic imaging systems and methods employing a fast target-oriented illumination calculation | |
| US9632192B2 (en) | Method of processing seismic data by providing surface offset common image gathers | |
| US7460437B2 (en) | Seismic data processing method and system for migration of seismic signals incorporating azimuthal variations in the velocity | |
| CN107817523B (en) | The analysis method and device of diffracted wave migration velocity | |
| CN105549081A (en) | Anisotropic medium common shot domain Gaussian beam migration imaging method | |
| CN102841376A (en) | Retrieval method for chromatography speed based on undulating surface | |
| CN106226818A (en) | Seismic data processing technique and device | |
| Hua et al. | Parsimonious 2D prestack Kirchhoff depth migration | |
| CN115469362B (en) | Energy flow density vector calculation method in seismic exploration | |
| CN105629299A (en) | Travel-time table and angle table acquisition method for angle domain prestack depth migration and imaging method | |
| CN107656308B (en) | A kind of common scattering point pre-stack time migration imaging method based on time depth scanning | |
| Toxopeus et al. | Simulating migrated and inverted seismic data by filtering a geologic model | |
| CN111999770B (en) | TTI medium conversion PS wave precise beam offset imaging method and system | |
| CN108693560B (en) | Scattered wave imaging method and system based on cross-correlation channel | |
| CN102798888B (en) | Method for calculating velocity ratio of longitudinal wave to transverse wave by using non-zero wellhead distance data | |
| CN114647003B (en) | Method and equipment for calculating residual time difference of converted wave angle domain common imaging point gather | |
| RU2705519C2 (en) | Method of producing migrated seismic images of geologic media based on 2d seismic survey data | |
| Pavlis | Imaging the earth with passive seismic arrays | |
| CN111999769B (en) | Complex surface anisotropy multicomponent seismic data prestack depth migration method | |
| CN113075734B (en) | Residual curvature spectrum calculation method and device based on signal-to-noise ratio constraint | |
| CN112379431B (en) | PS wave seismic data migration imaging method and system under complex surface condition | |
| RU2747628C1 (en) | Method for determining slope angle of reflective borders according to cdpm-2d data | |
| CN113391351B (en) | Method for extracting mine collection area structure based on passive source seismic wave field analysis |
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 |