CN110857998B - Elastic reverse time migration method and system based on lowrank finite difference - Google Patents

Elastic reverse time migration method and system based on lowrank finite difference Download PDF

Info

Publication number
CN110857998B
CN110857998B CN201810967971.1A CN201810967971A CN110857998B CN 110857998 B CN110857998 B CN 110857998B CN 201810967971 A CN201810967971 A CN 201810967971A CN 110857998 B CN110857998 B CN 110857998B
Authority
CN
China
Prior art keywords
wave
lowrank
finite difference
representing
vector
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
Application number
CN201810967971.1A
Other languages
Chinese (zh)
Other versions
CN110857998A (en
Inventor
韩冬
李博
许璐
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Original Assignee
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China Petroleum and Chemical Corp, Sinopec Geophysical Research Institute filed Critical China Petroleum and Chemical Corp
Priority to CN201810967971.1A priority Critical patent/CN110857998B/en
Publication of CN110857998A publication Critical patent/CN110857998A/en
Application granted granted Critical
Publication of CN110857998B publication Critical patent/CN110857998B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
    • 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/362—Effecting static or dynamic corrections; Stacking
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/50—Corrections or adjustments related to wave propagation
    • G01V2210/51—Migration
    • G01V2210/512—Pre-stack
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/70—Other details related to processing
    • G01V2210/74—Visualisation of seismic data

Abstract

Disclosed are a method and a system for elastic reverse time migration based on lowrank finite difference, the method and the system comprising: constructing a lowrank finite difference coefficient, realizing forward and reverse elastic wave field continuation based on a discretization wave equation of the constructed lowrank finite difference coefficient, and obtaining PP and PS imaging sections by utilizing the obtained forward and reverse elastic wave fields and scalar product imaging conditions. The invention constructs the finite difference operator with both time and space by utilizing the plane wave theory and the lowrank decomposition method, can meet the high-precision elastic wave field continuation in the time and space directions, fully utilizes multi-component seismic data, can obtain vector longitudinal waves and transverse waves in the continuation process, keeps the vector characteristic of an underground wave field, avoids the problems of different wave type crosstalk and converted wave polarity correction in the conventional elastic reverse time migration in the imaging process, and improves the precision of the elastic reverse time migration imaging.

Description

Elastic reverse time migration method and system based on lowrank finite difference
Technical Field
The invention relates to the technical field of seismic data prestack depth migration imaging, in particular to an elastic reverse time migration method and system based on lowrank finite difference.
Background
Elastic reverse time migration is a migration imaging technique developed for multi-component seismic data, which makes full use of the vector characteristics of elastic wave fields and can provide longitudinal wave (P-wave) and converted transverse wave (S-wave) imaging results. In the elastic reverse time migration process, elastic wave field continuation is one of the key technologies, and the precision of the elastic wave field continuation is directly related to the imaging quality. The conventional finite difference method is a commonly used wave field continuation method. In order to meet the high-precision wave field continuation, when the grid space is dispersed by the finite difference method, the selected sampling interval is usually about 8-10 times smaller than the wavelength of the actual seismic wavelet, and is generally oversampling; meanwhile, in order to satisfy the stability condition of Courant-Friedrch-Lewy (CFL), the time sampling interval also needs to take a smaller value, thereby causing the increase of the calculation amount of the algorithm. On the other hand, although the pseudo-spectral operator can be viewed as a spatial extreme of the finite difference operator, it is difficult to accommodate the wavefield prolongation of inhomogeneous media. With the development of exploration technology, the field attaches more and more importance to the wave field continuation precision. The conventional finite difference method is difficult to meet the requirement of high-precision elastic reverse time migration, and therefore, it is necessary to develop an elastic reverse time migration method and system based on lowrank finite difference.
The information disclosed in this background section is only for enhancement of understanding of the general background of the invention and should not be taken as an acknowledgement or any form of suggestion that this information forms the prior art already known to a person skilled in the art.
Disclosure of Invention
The invention provides an elastic reverse time migration method and system based on lowrank finite difference, which can construct a high-precision elastic wave field prolongation operator with both time and space in a time-space domain by means of a plane wave field decomposition method, a lowrank decomposition method, a Fourier transform method and other methods, and can obtain a high-precision elastic vector wave field with longitudinal and transverse wave decoupling based on the operator, thereby improving the elastic reverse time migration quality.
According to an aspect of the present invention, an elastic reverse time migration method based on lowrank finite difference is provided. The method may include:
1) construction of lowrank finite difference coefficient
Figure BDA0001775439950000021
Where i is 1,3, j is 1,3, P, S represents P wave, S wave, x (x, z) represents a space coordinate vector, and k (k) represents a space coordinate vectorx,kz) The wave number vector is represented by a vector of wave numbers,
Figure BDA0001775439950000022
is a matrix of coefficients, B (ξ)lK) is a matrix of vector coefficients,
Figure BDA0001775439950000023
is used as an index set, and the index set,
Figure BDA0001775439950000024
the number of points on the 1 st and 3 rd Cartesian coordinate components is represented, L represents a summation variable, and L represents a summation final value;
2) obtaining a discretization wave equation based on the lowrank finite difference coefficient constructed in the step 1) to realize forward and reverse elastic wave field continuation;
3) and (3) acquiring PP and PS imaging sections by utilizing the forward and reverse elastic vector wave fields acquired in the step 2) and utilizing scalar product imaging conditions.
Preferably, in step 1), the lowrank finite difference coefficient is constructed by:
1-1) determining an elastic wave expression propagating in a non-homogeneous isotropic medium:
Figure BDA0001775439950000025
wherein u ═ u (u)x,uz) Representing a displacement vector, ux,uzRespectively representing the displacement components in the x-direction and the z-direction in a rectangular coordinate system, t representing time, cp(x) And cs(x) Respectively representing a P-wave velocity and an S-wave velocity, and x ═ x, z represents a space coordinate vector;
1-2) carrying out Fourier transform on the formula (2) to obtain:
Figure BDA0001775439950000026
Figure BDA0001775439950000027
wherein k is (k)x,kz) Representing wave number vector, kxAnd kzRespectively representing the wave number variables in the x-direction and the z-direction in a rectangular coordinate system,
Figure BDA0001775439950000031
representing displacement vectors in the wavenumber domain, cp、csP-wave and S-wave velocities in the wavenumber domain, respectively;
1-3) obtaining a decoupled elastic wave recursive time integral equation based on the formula (3), wherein the expression is as follows:
Figure BDA0001775439950000032
wherein, UpRepresenting P-wave displacement vector, UsRepresents the displacement vector of the S-wave,
Figure BDA0001775439950000033
to normalize the wavenumber vector, Δ t represents the time interval,
Figure BDA0001775439950000034
in order to recur the time-integration operator,
Figure BDA0001775439950000035
or
Figure BDA0001775439950000036
Respectively representing wave number variables in the x direction and the z direction in a rectangular coordinate system, wherein i is 1,3, and j is 1, 3;
1-4) the recursive time integration operator in step 1-3) is represented in the form of a matrix as follows:
Figure BDA0001775439950000037
wherein, cp,s(x) Representing P-waves or S-wavesThe speed of the motor is controlled by the speed of the motor,
Figure BDA0001775439950000038
a P wave or S wave recursion time integral operator;
1-5) decomposing the formula (5) to obtain:
Figure BDA0001775439950000039
wherein m and n both represent a matrix
Figure BDA00017754399500000310
The rank of (c) is determined,
Figure BDA00017754399500000311
in the form of a spatial correlation matrix, the correlation matrix,
Figure BDA00017754399500000312
is a wave number correlation matrix, amnIs an intermediate matrix;
1-6) decomposing the formula (6) to obtain
Figure BDA00017754399500000313
Wherein the content of the first and second substances,
Figure BDA00017754399500000314
is a coefficient matrix;
1-7) defining a coefficient matrix
Figure BDA00017754399500000315
Wherein M represents the final value of summation, N represents the final value of summation, and is substituted into formula (6) to finally obtain lowrank finite difference coefficient
Figure BDA00017754399500000316
Preferably, the discretization wave equation of the lowrank finite difference coefficient is obtained by substituting the formula (1) into the formula (4) in the step 2), so as to realize forward and reverse elastic wave field continuation, wherein the expression is as follows:
Figure BDA0001775439950000041
wherein the content of the first and second substances,
Figure BDA0001775439950000042
uxand uzRepresenting the x-and z-direction wave field components, u, in the spatial domainpx、upz、usxAnd uszRepresenting P-wave or S-wave field components in the x-direction and z-direction respectively in the spatial domain,
Figure BDA0001775439950000043
respectively representing the P-wave correlation x-direction and the second-order lowrank finite difference coefficient for the z-direction,
Figure BDA0001775439950000044
is the second-order lowrank finite difference coefficient of P wave relative to x and z directions,
Figure BDA0001775439950000045
respectively representing the second-order lowrank finite difference coefficients of the relevant x and z directions of the S wave,
Figure BDA0001775439950000046
the coefficients are second-order lowrank finite difference coefficients of S-wave related in x and z directions.
Preferably, in step 3), the PP and PS imaging profiles are obtained by the following formula:
Figure BDA0001775439950000047
wherein, IppRepresenting PP common image point, IpsRepresenting a PS common imaging point, Ss(x, t) is the S-wave source wavefield, Sp(x, t) is the P-wave source wavefield, Rp(x, t) is the P-wave detection wave field, Rs(x, t) is the S-wave detection wavefield.
Preferably, the vector coefficient matrix B (x)lThe expression of k) is:
Figure BDA0001775439950000048
wherein the content of the first and second substances,
Figure BDA0001775439950000049
Figure BDA00017754399500000410
represents the number of points, k, on the 1,3 th Cartesian coordinate componentjIs the jth wavenumber component, Δ xjIs the spatial sampling interval.
Preferably, the coefficient matrix
Figure BDA00017754399500000411
The expression of (a) is:
Figure BDA00017754399500000412
wherein the content of the first and second substances,
Figure BDA0001775439950000051
representing the pseudo-inverse of the matrix.
According to another aspect of the present invention, a system for resilient reverse time migration based on lowrank finite difference is proposed, on which a computer program is stored, wherein the program, when executed by a processor, performs the steps of:
step 1: construction of lowrank finite difference coefficients
Figure BDA0001775439950000052
Where i is 1,3, j is 1,3, P, S represents P wave, S wave, x (x, z) represents a space coordinate vector, and k (k) represents a space coordinate vectorx,kz) The wave number vector is represented by a vector of wave numbers,
Figure BDA0001775439950000053
is a matrix of difference coefficients, B (xi)lK) is a vector of trigonometric functions,
Figure BDA0001775439950000054
is used as an index set, and the index set,
Figure BDA0001775439950000055
the number of points on the 1 st and 3 rd Cartesian coordinate components is represented, L represents a summation variable, and L represents a summation final value;
step 2: obtaining a discretization wave equation based on the lowrank finite difference coefficient constructed in the step 1 to realize the continuation of the forward and reverse elastic wave fields;
and step 3: and (3) acquiring PP and PS imaging sections by utilizing the forward and reverse elastic vector wave fields acquired in the step (2) and utilizing scalar product imaging conditions.
Preferably, in step 1, the lowrank finite difference coefficient is constructed by:
step 1-1: determining an elastic wave expression propagating in a non-uniform isotropic medium:
Figure BDA0001775439950000056
wherein u ═ u (u)x,uz) Representing a displacement vector, ux,uzRespectively representing the displacement components in the x-direction and the z-direction in a rectangular coordinate system, t representing time, cp(x) And cs(x) Respectively representing a P-wave velocity and an S-wave velocity, and x ═ x, z represents a space coordinate vector;
step 1-2: performing Fourier transform on the formula (2) to obtain:
Figure BDA0001775439950000057
wherein k is (k)x,kz) Representing wave number vector, kxAnd kzRespectively representing the wave number variables in the x-direction and the z-direction in a rectangular coordinate system,
Figure BDA0001775439950000058
representing displacement vectors in the wavenumber domain, cp、csP-wave and S-wave velocities in the wavenumber domain, respectively;
step 1-3: obtaining a decoupled elastic wave recursive time integral equation based on the formula (3), wherein the expression is as follows:
Figure BDA0001775439950000061
wherein, UpRepresenting P-wave displacement vector, UsRepresents the displacement vector of the S-wave,
Figure BDA0001775439950000062
to normalize the wavenumber vector, Δ t represents the time interval,
Figure BDA0001775439950000063
in order to recur the time-integration operator,
Figure BDA0001775439950000064
or
Figure BDA0001775439950000065
Respectively representing wave number variables in the x direction and the z direction in a rectangular coordinate system, wherein i is 1,3, and j is 1, 3;
step 1-4: the recursive time integration operator in steps 1-3 is represented in the form of a matrix as follows:
Figure BDA0001775439950000066
wherein, cp,s(x) Representing the velocity of the P-wave or S-wave,
Figure BDA0001775439950000067
a P wave or S wave recursion time integral operator;
step 1-5: decomposing equation (5) to obtain:
Figure BDA0001775439950000068
wherein m and n both represent a matrix
Figure BDA0001775439950000069
The rank of (c) is determined,
Figure BDA00017754399500000610
in the form of a spatial correlation matrix, the correlation matrix,
Figure BDA00017754399500000611
is a wave number correlation matrix, amnIs an intermediate matrix;
step 1-6: decomposing the formula (6) to obtain
Figure BDA00017754399500000612
Wherein the content of the first and second substances,
Figure BDA00017754399500000613
is a coefficient matrix;
step 1-7: defining a matrix of coefficients
Figure BDA00017754399500000614
Wherein M represents the final value of summation, N represents the final value of summation, and is substituted into formula (6) to finally obtain lowrank finite difference coefficient
Figure BDA00017754399500000615
Preferably, in step 2, the discretization wave equation of the lowrank finite difference coefficient is obtained by substituting the formula (1) into the formula (4), so as to realize forward and reverse elastic wave field continuation, wherein the expression is as follows:
Figure BDA0001775439950000071
wherein the content of the first and second substances,
Figure BDA0001775439950000072
uxand uzRepresenting the x-and z-direction wave field components, u, in the spatial domainpx、upz、usxAnd uszRepresenting P-wave or S-wave field components in the x-direction and z-direction respectively in the spatial domain,
Figure BDA0001775439950000073
second-order lowrank finite difference respectively representing P-wave correlation in x-direction and z-directionThe coefficients of which are such that,
Figure BDA0001775439950000074
is the second-order lowrank finite difference coefficient of P wave relative to x and z directions,
Figure BDA0001775439950000075
respectively representing the second-order lowrank finite difference coefficients of the relevant x and z directions of the S wave,
Figure BDA0001775439950000076
the coefficients are second-order lowrank finite difference coefficients of S-wave related in x and z directions.
Preferably, in step 3, the PP and PS imaging profiles are obtained by the following formula:
Figure BDA0001775439950000077
wherein, IppRepresenting PP common image point, IpsRepresenting a PS common imaging point, Ss(x, t) is the S-wave source wavefield, Sp(x, t) is the P-wave source wavefield, Rp(x, t) is the P-wave detection wave field, Rs(x, t) is the S-wave detection wavefield.
The present invention has other features and advantages which will be apparent from or are set forth in detail in the accompanying drawings and the following detailed description, which are incorporated herein, and which together serve to explain certain principles of the invention.
Drawings
The above and other objects, features and advantages of the present invention will become more apparent by describing in more detail exemplary embodiments thereof with reference to the attached drawings, in which like reference numerals generally represent like parts.
FIG. 1 is a flow chart illustrating the steps of a lowrank finite difference based elastic reverse time migration method according to the present invention;
FIG. 2(a) is a P-wave field snapshot obtained by a finite difference method;
FIG. 2(b) is a snapshot of the S-wave field obtained by the finite difference method;
FIG. 2(c) is a P-wave field snapshot obtained by lowrank finite difference method;
FIG. 2(d) is a S-wave field snapshot obtained by lowrank finite difference method;
FIG. 3 is a vertical slice model and parameters;
FIG. 4(a) is a P-wave horizontal component wave field snapshot obtained by the lowrank finite difference method;
FIG. 4(b) is a S-wave horizontal component wave field snapshot obtained by the lowrank finite difference method;
FIG. 4(c) is a P-wave vertical component wave field snapshot obtained by lowrank finite difference method;
FIG. 4(d) is a S-wave vertical component wave field snapshot obtained by lowrank finite difference method;
FIG. 5(a) is a horizontal component synthetic seismic record obtained by the lowrank finite difference method;
FIG. 5(b) is a vertical component synthetic seismic record obtained by the lowrank finite difference method;
FIG. 6(a) is a cross section of a vertical slice model multi-shot stacked PP imaging;
FIG. 6(b) is a vertical slice model multi-shot stacked PS imaging profile;
FIG. 7 is a salt dome model;
FIG. 8(a) is a salt dome model multi-shot stack PP imaging profile;
FIG. 8(b) is a salt dome model multiple shot stack PS imaging profile.
Detailed Description
The invention will be described in more detail below with reference to the accompanying drawings. While the preferred embodiments of the present invention are shown in the drawings, it should be understood that the present invention 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 invention to those skilled in the art.
Fig. 1 shows a flow chart of the steps of a lowrank finite difference-based elastic reverse time migration method according to the present invention.
In this embodiment, an elastic reverse time migration method based on lowrank finite difference according to the present invention may include:
step 101, constructing lowrank finite difference coefficient
Figure BDA0001775439950000091
Where i is 1,3, j is 1,3, P, S represents P wave, S wave, x (x, z) represents a space coordinate vector, and k (k) represents a space coordinate vectorx,kz) The wave number vector is represented by a vector of wave numbers,
Figure BDA0001775439950000092
is a matrix of difference coefficients, B (xi)lK) is a vector of trigonometric functions,
Figure BDA0001775439950000093
is used as an index set, and the index set,
Figure BDA0001775439950000094
the number of points on the 1 st and 3 rd Cartesian coordinate components is represented, L represents a summation variable, and L represents a summation final value;
in an exemplary embodiment, the lowrank finite difference coefficient is constructed by:
in the time-space domain, the expression of an elastic wave propagating in an infinitely large, non-uniform isotropic medium is determined:
Figure BDA0001775439950000095
wherein u ═ u (u)x,uz) Representing a displacement vector, ux,uzRespectively representing the displacement components in the x-direction and the z-direction in a rectangular coordinate system, t representing time, cp(x) And cs(x) Respectively, P-wave velocity and S-wave velocity, x ═ x, z denotes a space coordinate vector,
Figure BDA0001775439950000096
and
Figure BDA0001775439950000097
gradient, divergence and rotation operators respectively;
performing Fourier transform on the formula (2) to obtain:
Figure BDA0001775439950000098
wherein k is (k)x,kz) Representing wave number vector, kxAnd kzRespectively representing the wave number variables in the x-direction and the z-direction in a rectangular coordinate system,
Figure BDA0001775439950000099
representing displacement vectors in the wavenumber domain, cp、csP-wave and S-wave velocities in the wavenumber domain, respectively;
according to the wave field decomposition theory and the plane wave theory, a decoupled elastic wave recursive time integral equation in a space-wave number domain is obtained based on a formula (3), and the expression is as follows:
Figure BDA0001775439950000101
wherein, UpRepresenting P-wave displacement vector, UsRepresents the displacement vector of the S-wave,
Figure BDA0001775439950000102
to normalize the wavenumber vector, Δ t represents the time interval,
Figure BDA0001775439950000103
in order to recur the time-integration operator,
Figure BDA0001775439950000104
or
Figure BDA0001775439950000105
Respectively representing wave number variables in the x direction and the z direction in a rectangular coordinate system, wherein i is 1,3, and j is 1, 3;
the recursive time integration operator is expressed in the form of a matrix as follows:
Figure BDA0001775439950000106
wherein, cp,s(x) Representing the velocity of the P-wave or S-wave,
Figure BDA0001775439950000107
a P wave or S wave recursion time integral operator;
specifically, because the calculation amount of the elastic recursion time integral operator is large, the complexity is the same as the number N of discrete grid pointsxThe square order of the wave field is equivalent, so that the huge calculation amount is not beneficial to the elastic wave field continuation, therefore, the recursive time integral operator is approximated by using the lowrank decomposition method to obtain the following matrix form:
Figure BDA0001775439950000108
wherein m and n both represent a matrix
Figure BDA0001775439950000109
The rank of (c) is determined,
Figure BDA00017754399500001010
in the form of a spatial correlation matrix, the correlation matrix,
Figure BDA00017754399500001011
is a wave number correlation matrix, amnIs an intermediate matrix; the complexity is the same as NxLinear correlation, matrix
Figure BDA00017754399500001012
Only the wave number is related, and the formula (6) is decomposed to obtain
Figure BDA00017754399500001013
Figure BDA00017754399500001014
Wherein the content of the first and second substances,
Figure BDA00017754399500001015
is a coefficient matrix;
specifically, vector coefficient matrix B (x)lThe expression of k) is:
Figure BDA00017754399500001016
wherein the content of the first and second substances,
Figure BDA00017754399500001017
is used as an index set, and the index set,
Figure BDA00017754399500001018
represents the number of points, k, on the 1,3 th Cartesian coordinate componentjIs the jth wavenumber component, Δ xjFor the spatial sampling interval, j is 1, and 3 corresponds to the x and z directions in cartesian coordinates.
In particular, a coefficient matrix
Figure BDA00017754399500001019
The expression of (a) is:
Figure BDA0001775439950000111
wherein the content of the first and second substances,
Figure BDA0001775439950000112
representing the pseudo-inverse of the matrix.
Defining a matrix of coefficients
Figure BDA0001775439950000113
Wherein M represents the final value of summation, N represents the final value of summation, and is substituted into formula (6) to finally obtain lowrank finite difference coefficient
Figure BDA0001775439950000114
The obtained lowrank finite difference coefficient depends on P-wave and S-wave velocity parameters, and can also be called an elastic lowrank finite difference operator.
Specifically, we use formula (1) to construct the lowrank finite difference coefficients of P-wave correlation and S-wave correlation, respectively.
102, obtaining a discretization wave equation based on the lowrank finite difference coefficient constructed in the step 101 to realize forward and reverse elastic wave field continuation;
specifically, the formula (1) is substituted into the formula (3) to obtain the discretization wave equation based on the lowrank finite difference coefficient, and the expression is as follows:
Figure BDA0001775439950000115
wherein the content of the first and second substances,
Figure BDA0001775439950000116
uxand uzRepresenting the x-and z-direction wave field components, u, in the spatial domainpx、upz、usxAnd uszRepresenting P-wave or S-wave field components in the x-direction and z-direction respectively in the spatial domain,
Figure BDA0001775439950000117
respectively representing the second-order lowrank finite difference coefficients of the P-wave in the x and z directions,
Figure BDA0001775439950000118
is the second-order lowrank finite difference coefficient of P wave relative to x and z directions,
Figure BDA0001775439950000119
second-order lowrank finite difference coefficients of the relevant x and z directions of the S wave respectively,
Figure BDA00017754399500001110
the coefficients are second-order lowrank finite difference coefficients of S-wave related in x and z directions.
Specifically, a vector elastic wave field of high-precision longitudinal and transverse wave decoupling at each moment is respectively constructed by using a formula (8), a velocity model and seismic data.
Step 103, using the forward and backward elastic wave fields obtained in step 102, and using scalar product imaging conditions to obtain PP and PS imaging profiles.
In one example, for a source detection field, an imaging process is realized by means of vector inner product and normalization. The method can ensure the precision of the P/S wave vector wave field, avoid the phenomenon of P/S wave crosstalk and the phenomenon of polarity reversal of a converted wave section in the imaging process, and effectively improve the imaging precision.
Specifically, PP and PS imaging profiles were obtained by the following formula:
Figure BDA0001775439950000121
wherein, IppRepresenting PP common image point, IpsRepresenting a PS common imaging point, Ss(x, t) is the S-wave source wavefield, Sp(x, t) is the P-wave source wavefield, Rp(x, t) is the P-wave detection wave field, Rs(x, t) is the S-wave detection wavefield.
According to the method, the high-precision mixed domain elastic continuation operator which is time-space compatible is constructed in the space-wave number domain, and the high-precision elastic wave field continuation operator in the time-space domain is formed by means of the lowrank decomposition method and the Fourier transform method, so that the precision of the elastic wave field continuation is ensured. The method can obtain high-precision longitudinal and transverse wave vector wave fields, is favorable for improving the imaging quality, and has important significance for promoting high-precision imaging.
Application example
To facilitate understanding of the solution of the embodiments of the present invention and the effects thereof, a specific application example is given below. It will be understood by those skilled in the art that this example is merely for the purpose of facilitating an understanding of the present invention and that any specific details thereof are not intended to limit the invention in any way.
As shown in fig. 2(a), fig. 2(b) is a P-wave horizontal component wave field snapshot obtained by a conventional finite difference method, fig. 2(c) is a P-wave horizontal component wave field snapshot obtained by a lowrank finite difference method, fig. 2(d) is a S-wave horizontal component wave field snapshot obtained by a lowrank finite difference method, which is the same as the conditions of time order 2 and space order 4, it can be seen by comparison that a wave field snapshot obtained by the conventional finite difference method suffers from a serious dispersion effect and cannot be directly applied to an elastic reverse time migration process, and a wave field snapshot obtained by the elastic lowrank finite difference method is hardly affected by the dispersion effect, thereby showing that the elastic lowrank finite difference has effectiveness in suppressing the dispersion.
And performing elastic wave field continuation on the vertical fault model by using an elastic lowrank finite difference method with the same order number and applying an elastic reverse time migration process. FIG. 3 shows the geometric structure and model parameters of a vertical fault model, wherein a fault interface divides the whole model medium into two parts, the longitudinal wave velocity is 2800m/s, the transverse wave velocity is 1800m/s, the transverse wave velocity is 3400m/s and the transverse wave velocity is 2200 m/s. The sampling interval in the horizontal direction and the sampling interval in the vertical direction are both 10m, a Gaussian wavelet with a larger time sampling interval of 1.5ms and a main frequency of 25Hz is used as a seismic source wavelet, and a seismic source point is located at (3000,800) m. Elastic wave field continuation is performed by using an elastic lowrank finite difference method, and a P-wave horizontal component wave field snapshot as shown in fig. 4(a), an S-wave horizontal component wave field snapshot as shown in fig. 4(b), a P-wave vertical component wave field snapshot as shown in fig. 4(c), and an S-wave vertical component wave field snapshot as shown in fig. 4(d) can be obtained, and further, a horizontal component synthetic seismic record as shown in fig. 5(a) and a vertical component synthetic seismic record as shown in fig. 5(b) can be obtained. The specific forms of P waves and S waves are clear and visible when the P waves and the S waves are transmitted in underground media, and special wave forms such as converted waves and prism waves are clear and depicted. As shown in fig. 6(a), the vertical slice model multi-shot superposition PP imaging section and fig. 6(b) the vertical slice model multi-shot superposition PS imaging section show that elastic reverse time migration based on lowrank finite difference method can clearly depict the slice, the slice morphology position expression is accurate, and the converted wave section has no polarity inversion phenomenon. Therefore, the elastic lowrank finite difference method can be suitable for elastic wave field continuation of complex media with large time step length and meets the requirement of elastic reverse time migration.
Fig. 7 shows an SEG/EAGE salt dome model, which mainly has the difficulty that the structures under salt are difficult to focus and image, elastic wave field extension is performed by using an elastic lowrank finite difference method, and a final imaging process is realized by using scalar product imaging conditions, so that a salt dome model multi-shot stack PP imaging section shown in fig. 8(a) and a salt dome model multi-shot stack PS imaging section shown in fig. 8(b) can be obtained, and the structures on the side wall of the salt dome and under the salt dome are clearly carved.
In conclusion, the invention utilizes the elastic lowrank finite difference operator to carry out forward and reverse elastic wave field continuation, directly obtains the high-precision longitudinal and transverse wave (P/S wave) source vector wave field in the continuation process, does not need to apply an additional wave field separation operator, simultaneously obtains the imaging result without polarity correction after applying a scalar product imaging condition, and solves a series of problems that the wave field continuation precision is low when the conventional finite difference method is mainly utilized in the conventional elastic reverse time migration, the scalar P wave and S wave fields are obtained by additionally utilizing divergence and rotation operators, the polarity correction processing is needed to be carried out on a converted wave imaging section in the imaging process, and the like. The elastic reverse time migration method based on lowrank finite difference can meet the requirement of high-precision elastic wave migration imaging.
It will be appreciated by persons skilled in the art that the above description of embodiments of the invention is intended only to illustrate the benefits of embodiments of the invention and is not intended to limit embodiments of the invention to any examples given.
Having described embodiments of the present invention, the foregoing description is intended to be exemplary, not exhaustive, and not limited to the embodiments disclosed. Many modifications and variations will be apparent to those of ordinary skill in the art without departing from the scope and spirit of the described embodiments.

Claims (8)

1. An elastic reverse time migration method based on lowrank finite difference is characterized by comprising the following steps:
1) construction of lowrank finite difference coefficients
Figure FDA0003190026090000011
Wherein, i is 1,3, j is 1,3, P, S represents P wave, S wave, x is (x, z) representsSpace coordinate vector, k ═ k (k)x,kz) The wave number vector is represented by a vector of wave numbers,
Figure FDA0003190026090000012
is a matrix of difference coefficients, B (xi)lK) is a vector of trigonometric functions,
Figure FDA0003190026090000013
is used as an index set, and the index set,
Figure FDA0003190026090000014
the number of points on the 1 st and 3 rd Cartesian coordinate components is represented, L represents a summation variable, and L represents a summation final value;
2) obtaining a discretization wave equation based on the lowrank finite difference coefficient constructed in the step 1) to realize forward and reverse elastic wave field continuation;
3) acquiring PP and PS imaging sections by utilizing the forward and reverse elastic vector wave fields acquired in the step 2) and scalar product imaging conditions;
in step 3), PP and PS imaging profiles are obtained by the following formula:
Figure FDA0003190026090000015
wherein, IppRepresenting PP common image point, IpsRepresenting a PS common imaging point, Ss(x, t) is the S-wave source wavefield, Sp(x, t) is the P-wave source wavefield, Rp(x, t) is the P-wave detection wave field, Rs(x, t) is the S-wave detection wavefield.
2. A lowrank finite difference-based elastic reverse time migration method according to claim 1, wherein in step 1), the lowrank finite difference coefficient is constructed by:
1-1) determining an elastic wave expression propagating in a non-homogeneous isotropic medium:
Figure FDA0003190026090000021
wherein u ═ u (u)x,uz) Representing a displacement vector, ux,uzRespectively representing the displacement components in the x-direction and the z-direction in a rectangular coordinate system, t representing time, cp(x) And cs(x) Respectively representing a P-wave velocity and an S-wave velocity, and x ═ x, z represents a space coordinate vector;
1-2) carrying out Fourier transform on the formula (2) to obtain:
Figure FDA0003190026090000022
Figure FDA0003190026090000023
wherein k is (k)x,kz) Representing wave number vector, kxAnd kzRespectively representing the wave number variables in the x-direction and the z-direction in a rectangular coordinate system,
Figure FDA0003190026090000024
representing displacement vectors in the wavenumber domain, cp、csP-wave and S-wave velocities in the wavenumber domain, respectively;
1-3) obtaining a decoupled elastic wave recursive time integral equation based on the formula (3), wherein the expression is as follows:
Figure FDA0003190026090000025
wherein, UpRepresenting P-wave displacement vector, UsRepresents the displacement vector of the S-wave,
Figure FDA0003190026090000026
to normalize the wavenumber vector, Δ t represents the time interval,
Figure FDA0003190026090000027
in order to recur the time-integration operator,
Figure FDA0003190026090000028
or
Figure FDA0003190026090000029
Respectively representing wave number variables in the x direction and the z direction in a rectangular coordinate system, wherein i is 1,3, and j is 1, 3;
1-4) the recursive time integration operator in step 1-3) is represented in the form of a matrix as follows:
Figure FDA00031900260900000210
wherein, cp,s(x) Representing the velocity of the P-wave or S-wave,
Figure FDA00031900260900000211
a P wave or S wave recursion time integral operator;
1-5) decomposing the formula (5) to obtain:
Figure FDA00031900260900000212
wherein m and n both represent a matrix
Figure FDA00031900260900000213
The rank of (c) is determined,
Figure FDA00031900260900000214
for the spatial correlation matrix, the correlation matrix is,
Figure FDA00031900260900000215
is a wave number correlation matrix, amnIs an intermediate matrix;
1-6) decomposing the formula (6) to obtain
Figure FDA00031900260900000216
Wherein the content of the first and second substances,
Figure FDA00031900260900000217
is a coefficient matrix;
1-7) defining a coefficient matrix
Figure FDA0003190026090000031
Wherein M represents the final value of summation, N represents the final value of summation, and is substituted into formula (6) to finally obtain lowrank finite difference coefficient
Figure FDA0003190026090000032
3. A method for elastic reverse time migration based on lowrank finite difference according to claim 2, wherein the discretization wave equation of lowrank finite difference coefficient is obtained by substituting formula (1) into formula (4) in step 2), so as to realize forward and reverse elastic wave field continuation, and the expression is:
Figure FDA0003190026090000033
wherein the content of the first and second substances,
Figure FDA0003190026090000034
uxand uzRepresenting the x-and z-direction wave field components, u, in the spatial domainpx、upz、usxAnd uszRepresenting P-wave or S-wave field components in the x-direction and z-direction respectively in the spatial domain,
Figure FDA0003190026090000035
respectively representing the second-order lowrank finite difference coefficients of the P-wave correlation in the x direction and the z direction,
Figure FDA0003190026090000036
is the second-order lowrank finite difference coefficient of P wave relative to x and z directions,
Figure FDA0003190026090000037
second-order lowrank finite difference system respectively representing X direction and Z direction of S wave correlationThe number of the first and second groups is,
Figure FDA0003190026090000038
the coefficients are second-order lowrank finite difference coefficients of S-wave related in x and z directions.
4. A lowrank finite difference-based elastic reverse time migration method according to claim 2, wherein the vector coefficient matrix B (ξ)lThe expression of k) is:
Figure FDA0003190026090000039
wherein the content of the first and second substances,
Figure FDA0003190026090000041
represents the number of points, k, on the 1,3 th Cartesian coordinate component1Is the 1 st wave number component, k3Is the 3 rd wavenumber component, Δ x1、Δx2Is the spatial sampling interval.
5. A lowrank finite difference-based elastic reverse time migration method according to claim 2, wherein the coefficient matrix
Figure FDA0003190026090000042
The expression of (a) is:
Figure FDA0003190026090000043
wherein the content of the first and second substances,
Figure FDA0003190026090000044
representing the pseudo-inverse of the matrix.
6. A lowrank finite difference based elastic reverse time migration system having a computer program stored thereon, wherein said program when executed by a processor performs the steps of:
step 1: construction of lowrank finite difference coefficients
Figure FDA0003190026090000045
Where i is 1,3, j is 1,3, P, S represents P wave, S wave, x (x, z) represents a space coordinate vector, and k (k) represents a space coordinate vectorx,kz) The wave number vector is represented by a vector of wave numbers,
Figure FDA0003190026090000046
is a matrix of difference coefficients, B (xi)lK) is a vector of trigonometric functions,
Figure FDA0003190026090000047
is used as an index set, and the index set,
Figure FDA0003190026090000048
the number of points on the 1 st and 3 rd Cartesian coordinate components is represented, L represents a summation variable, and L represents a summation final value;
step 2: obtaining a discretization wave equation based on the lowrank finite difference coefficient constructed in the step 1 to realize the continuation of the forward and reverse elastic wave fields;
and step 3: acquiring PP and PS imaging sections by utilizing the forward and reverse elastic vector wave fields acquired in the step 2 and utilizing scalar product imaging conditions;
in step 3), PP and PS imaging profiles are obtained by the following formula:
Figure FDA0003190026090000049
wherein, IppRepresenting PP common image point, IpsRepresenting a PS common imaging point, Ss(x, t) is the S-wave source wavefield, Sp(x, t) is the P-wave source wavefield, Rp(x, t) is the P-wave detection wave field, Rs(x, t) is the S-wave detection wave field
7. A lowrank finite difference-based elastic reverse time migration system according to claim 6, wherein in step 1, the lowrank finite difference coefficients are constructed by:
step 1-1: determining an expression for an elastic wave propagating in a non-homogeneous isotropic medium:
Figure FDA0003190026090000051
wherein u ═ u (u)x,uz) Representing a displacement vector, ux,uzRespectively representing the displacement components in the x-direction and the z-direction in a rectangular coordinate system, t representing time, cp(x) And cs(x) Respectively representing a P-wave velocity and an S-wave velocity, and x ═ x, z represents a space coordinate vector;
step 1-2: performing Fourier transform on the formula (2) to obtain:
Figure FDA0003190026090000052
wherein k is (k)x,kz) Representing wave number vector, kxAnd kzRespectively representing the wave number variables in the x-direction and the z-direction in a rectangular coordinate system,
Figure FDA0003190026090000053
representing displacement vectors in the wavenumber domain, cp、csP-wave and S-wave velocities in the wavenumber domain, respectively;
step 1-3: obtaining a decoupled elastic wave recursive time integral equation based on the formula (3), wherein the expression is as follows:
Figure FDA0003190026090000054
wherein, UpRepresenting P-wave displacement vector, UsRepresents the displacement vector of the S-wave,
Figure FDA0003190026090000055
to normalize the wavenumber vector, Δ t represents the time interval,
Figure FDA0003190026090000056
in order to recur the time-integration operator,
Figure FDA0003190026090000057
or
Figure FDA0003190026090000058
Respectively representing wave number variables in the x direction and the z direction in a rectangular coordinate system, wherein i is 1,3, and j is 1, 3;
step 1-4: the recursive time integration operator in steps 1-3 is represented in the form of a matrix as follows:
Figure FDA0003190026090000061
wherein, cp,s(x) Representing the velocity of the P-wave or S-wave,
Figure FDA0003190026090000062
a P wave or S wave recursion time integral operator;
step 1-5: decomposing equation (5) to obtain:
Figure FDA0003190026090000063
wherein m and n both represent a matrix
Figure FDA0003190026090000064
The rank of (c) is determined,
Figure FDA0003190026090000065
in the form of a spatial correlation matrix, the correlation matrix,
Figure FDA0003190026090000066
is a wave number correlation matrix, amnIs an intermediate matrix;
step 1-6: decomposing the formula (6) to obtain
Figure FDA0003190026090000067
Wherein the content of the first and second substances,
Figure FDA0003190026090000068
is a coefficient matrix;
step 1-7: defining a matrix of coefficients
Figure FDA0003190026090000069
Wherein M represents the final value of summation, N represents the final value of summation, and is substituted into formula (6) to finally obtain lowrank finite difference coefficient
Figure FDA00031900260900000610
8. A lowrank finite difference-based elastic reverse time migration system according to claim 6, wherein the forward and reverse elastic wave field prolongation is realized by substituting formula (1) into formula (4) to obtain the discretized wave equation of lowrank finite difference coefficient in step 2, which is expressed as:
Figure FDA00031900260900000611
wherein the content of the first and second substances,
Figure FDA00031900260900000612
uxand uzRepresenting the x-and z-direction wave field components, u, in the spatial domainpx、upz、usxAnd uszRepresenting P-wave or S-wave field components in the x-direction and z-direction respectively in the spatial domain,
Figure FDA00031900260900000613
respectively representing the second-order lowrank finite difference coefficients of the P-wave correlation in the x direction and the z direction,
Figure FDA0003190026090000071
is the second-order lowrank finite difference coefficient of P wave relative to x and z directions,
Figure FDA0003190026090000072
respectively representing the second-order lowrank finite difference coefficients of the relevant x and z directions of the S wave,
Figure FDA0003190026090000073
the coefficients are second-order lowrank finite difference coefficients of S-wave related in x and z directions.
CN201810967971.1A 2018-08-23 2018-08-23 Elastic reverse time migration method and system based on lowrank finite difference Active CN110857998B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810967971.1A CN110857998B (en) 2018-08-23 2018-08-23 Elastic reverse time migration method and system based on lowrank finite difference

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810967971.1A CN110857998B (en) 2018-08-23 2018-08-23 Elastic reverse time migration method and system based on lowrank finite difference

Publications (2)

Publication Number Publication Date
CN110857998A CN110857998A (en) 2020-03-03
CN110857998B true CN110857998B (en) 2021-11-05

Family

ID=69635189

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810967971.1A Active CN110857998B (en) 2018-08-23 2018-08-23 Elastic reverse time migration method and system based on lowrank finite difference

Country Status (1)

Country Link
CN (1) CN110857998B (en)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111158047B (en) * 2020-03-04 2021-05-11 中国石油大学(北京) Three-dimensional elastic wave field vector decomposition method, device and computer storage medium
CN112327359B (en) * 2020-10-14 2022-06-14 山东省科学院海洋仪器仪表研究所 Elastic reverse time migration method based on imaging energy flow vector
CN112904426B (en) * 2021-03-27 2022-09-30 中国石油大学(华东) Decoupling elastic wave reverse time migration method, system and application
CN116088049B (en) * 2023-04-07 2023-06-20 清华大学 Least square inverse time migration seismic imaging method and device based on wavelet transformation

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104122585A (en) * 2014-08-08 2014-10-29 中国石油大学(华东) Seismic forward modeling method based on elastic wave field vector decomposition and low-rank decomposition
CN104459773A (en) * 2014-08-08 2015-03-25 中国石油天然气集团公司 Unconditionally stable seismic wave field continuation method based on staggered grid Lowrank decomposition

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015104637A2 (en) * 2014-01-13 2015-07-16 Cgg Services Sa Device and method for deghosting seismic data using sparse tau-p inversion
WO2018004789A1 (en) * 2016-06-28 2018-01-04 Exxonbil Upstream Research Company Reverse time migration in anisotropic media with stable attenuation compensation

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104122585A (en) * 2014-08-08 2014-10-29 中国石油大学(华东) Seismic forward modeling method based on elastic wave field vector decomposition and low-rank decomposition
CN104459773A (en) * 2014-08-08 2015-03-25 中国石油天然气集团公司 Unconditionally stable seismic wave field continuation method based on staggered grid Lowrank decomposition

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Elastic wave propagation on a staggered-grid using lowrank finite-difference method;Dong Han,等;《SEG》;20161231;第3900-3904页 *
交错网格Lowrank有限差分及其在逆时偏移中的应用;方刚，等;《中国石油大学学报(自然科学版)》;20140430;第45-49页 *
基于lowrank分解的地震正演模拟与逆时偏移方法研究;方刚;《中国博士学位论文全文数据库》;20160715;第9、23-26、93-96页 *
基于Lowrank近似的TTI介质一步法qP波正演模拟;朱浩然，等;《CT理论与应用研究 》;20170228;第9-18页 *

Also Published As

Publication number Publication date
CN110857998A (en) 2020-03-03

Similar Documents

Publication Publication Date Title
CN110857998B (en) Elastic reverse time migration method and system based on lowrank finite difference
Yu et al. Monte Carlo data-driven tight frame for seismic data recovery
Chen et al. Iterative deblending of simultaneous-source seismic data using seislet-domain shaping regularization
Zhao et al. Strain Green’s tensors, reciprocity, and their applications to seismic source and structure studies
CN111221037B (en) Decoupling elastic reverse time migration imaging method and device
EP2689273A2 (en) System and method for seismic data modeling and migration
GB2548555A (en) A method of redatuming geophysical data
Yang et al. Seislet-based morphological component analysis using scale-dependent exponential shrinkage
Zand et al. Consensus optimization of total variation–based reverse time migration
Li et al. An lp-space matching pursuit algorithm and its application to robust seismic data denoising via time-domain radon transforms
Shragge Acoustic wave propagation in tilted transversely isotropic media: Incorporating topography
Wang et al. Robust nonstationary local slope estimation
Feng et al. Estimating primaries by sparse inversion of the 3D curvelet transform and the L1-norm constraint
Fan et al. Multisource least-squares reverse-time migration with structure-oriented filtering
Liu et al. Seismic data interpolation using generalised velocity‐dependent seislet transform
CN113866826A (en) Mixed domain seismic migration hessian matrix estimation method
CN111045084B (en) Multi-wave self-adaptive subtraction method based on prediction feature extraction
CN109725346B (en) Time-space high-order VTI rectangular grid finite difference method and device
Raknes et al. Challenges and solutions for performing 3D time-domain elastic full-waveform inversion
Kazemi et al. Block row recursive least-squares migration
Shi et al. Random Noise Attenuation of Common Offset Gathers Using Iteratively Reweighted $\ell_ {2, 1} $ Norm Minimization
Gholami et al. Simultaneous constraining of model and data smoothness for regularization of geophysical inverse problems
CN113341460B (en) Circulating minimization seismic data reconstruction method based on continuous operator splitting
CN112444849A (en) Elastic reverse time migration imaging method based on staggered grid low-rank finite difference
Burnett et al. Reversible Stolt migration

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