CN110927779A - Fault constraint tomography inversion method and inversion system - Google Patents
Fault constraint tomography inversion method and inversion system Download PDFInfo
- Publication number
- CN110927779A CN110927779A CN201811094492.XA CN201811094492A CN110927779A CN 110927779 A CN110927779 A CN 110927779A CN 201811094492 A CN201811094492 A CN 201811094492A CN 110927779 A CN110927779 A CN 110927779A
- Authority
- CN
- China
- Prior art keywords
- fault
- constrained
- preconditioned
- geological
- tomographic
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 75
- 238000003325 tomography Methods 0.000 title description 20
- 238000003384 imaging method Methods 0.000 claims abstract description 47
- 230000005012 migration Effects 0.000 claims abstract description 24
- 238000013508 migration Methods 0.000 claims abstract description 24
- 238000004587 chromatography analysis Methods 0.000 claims abstract description 20
- 239000011159 matrix material Substances 0.000 claims description 26
- 230000006870 function Effects 0.000 claims description 12
- 238000009499 grossing Methods 0.000 claims description 10
- 238000002598 diffusion tensor imaging Methods 0.000 claims description 9
- 238000000354 decomposition reaction Methods 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 claims description 6
- 238000011161 development Methods 0.000 abstract description 10
- 230000008569 process Effects 0.000 abstract description 5
- 238000012545 processing Methods 0.000 abstract description 3
- 238000010586 diagram Methods 0.000 description 7
- 230000008859 change Effects 0.000 description 6
- 238000010276 construction Methods 0.000 description 6
- 238000001914 filtration Methods 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 230000009545 invasion Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention discloses a fault constraint chromatography inversion method and an inversion system, wherein the method comprises the following steps: step 1: calculating subsurface geological dip and fault properties of the migration profile based on the seismic migration image; step 2: constructing a structural constraint preconditioned chromatographic equation set based on the underground geological inclination and fault attributes; and step 3: and solving the structural constraint preconditioned chromatographic equation set to obtain an underground velocity model and an offset profile. The fault constraint chromatographic inversion method realizes the automation of the whole chromatographic inversion iteration process, does not need manual intervention, reduces the manual workload, shortens the processing period and improves the working efficiency; the velocity model is constrained by introducing underground geological dip angle and fault attributes, so that the velocity model obtained by inversion has higher precision in a complex structural area, particularly a fault development area, and meets the requirement of subsequent offset imaging.
Description
Technical Field
The invention belongs to the field of seismic imaging and inversion, and particularly relates to a fault constraint tomography inversion method and an inversion system.
Background
The method is limited by the influence of factors such as an illumination aperture, an acquisition signal-to-noise ratio and the like, and the resolution of the traditional ground seismic imaging is lack of medium wave number band information. That is, the long wavelength and low frequency components of the subsurface velocity can be better recovered by conventional ray tomography, the high wavenumber reflection coefficient can be characterized by offset imaging, and the high precision velocity model of the medium wavenumber band is difficult to be accurately inverted.
With the breakthrough progress of the two-width one-height (broadband, wide azimuth and high density) seismic acquisition and chromatographic inversion method, the resolution of the traditional seismic imaging is widened. The two-width one-height seismic acquisition expands the seismic migration imaging resolution to a low wavenumber range, and the chromatographic inversion expands the velocity inversion resolution to a medium wavenumber range and a high wavenumber range.
In recent years, the high-resolution chromatographic inversion method has been greatly improved, and the inversion resolution limit of the conventional chromatographic method is improved from 2-3Hz to 6-8 Hz. The main ideas of these works are to improve the quality of tomographic automatic picking and to introduce geological formation constraints in the inversion process. The method for implementing geological structure constraint mainly comprises two types, wherein one type introduces a geological horizon interface and a geological structure in velocity model parameterization, and the method needs to extract horizon information of a seismic migration image and has higher precision requirement on an automatic seismic interpretation algorithm; another method is to construct a tomographic preconditioner containing geological structure information for guiding the velocity update direction, so that the algorithm implementation is more flexible and the offset image quality is relatively lower, and the structure constraint accuracy depends on the accuracy of the preconditioner. However, the method for constructing the chromatography preconditioner containing the geological structure information is to manually pick up the information of the underground geological horizon and the fault on the offset profile, and the method is time-consuming and labor-consuming in the actual operation process and even can not meet the efficiency requirement of exploration and development.
Therefore, there is a particular need for a fault-constrained tomographic inversion method that can automatically pick up subsurface geological horizons and fault information.
Disclosure of Invention
The invention aims to provide a fault constraint tomography inversion method and an inversion system which can automatically pick up underground geological position and fault information and have high accuracy.
In order to achieve the above object, the present invention provides a fault-confined tomographic inversion method, including: step 1: calculating subsurface geological dip and fault properties of the migration profile based on the seismic migration image; step 2: constructing a structural constraint preconditioned chromatographic equation set based on the underground geological inclination angle and the fault attribute; and step 3: and solving the structural constraint preconditioned chromatographic equation set to obtain an underground velocity model and an offset profile.
Preferably, the subsurface geologic dip of the offset section is calculated by a structure tensor algorithm.
Preferably, the fault attribute is calculated according to the following formula:
wherein, s (i)x,it) As a fault attribute, ixAnd jxAs a transverse coordinate variable, itAs a longitudinal coordinate variable, sn(ix,it) And sd(ix,it) Respectively the numerator and denominator of the fault attribute expression, g (i)x,it) For use in calculating tomographic properties (i)x,it) Amplitude value of p (i)x,it) For seismic imaging in (i)x,it) The tangential slope of the structure, MxThe lateral window length of the fault attribute is calculated.
Preferably, the structural constraint preconditioned chromatography equation set is:
STLTLSu＝STLTτ (2)
wherein, L is a ray chromatography kernel function, Δ m ═ Su, Δ m and Su are velocity model updating quantities, τ is the difference between the forward calculation seismic wave propagation travel time and the received data travel time, and S is a preconditioner.
Preferably, the preconditioner S is a smoothing operator containing geological structure information, S being:
S＝(I+DTGD)-1(3)
where I is the unit matrix and D is the gradient operator
DTIs the transpose matrix of D
G is the structure tensor operator.
Preferably, the structural constraint preconditioned set of tomographic equations is solved by a least squares QR decomposition method.
The invention also provides a fault constraint tomography inversion system, which comprises: a memory storing computer-executable instructions; a processor executing computer executable instructions in the memory to perform the steps of: step 1: calculating subsurface geological dip and fault properties of the migration profile based on the seismic migration image; step 2: constructing a structural constraint preconditioned chromatographic equation set based on the underground geological inclination angle and the fault attribute; and step 3: and solving the structural constraint preconditioned chromatographic equation set to obtain an underground velocity model and an offset profile.
Preferably, the subsurface geologic dip of the offset section is calculated by a structure tensor algorithm.
Preferably, the fault attribute is calculated according to the following formula:
wherein, s (i)x,it) As a fault attribute, ixAnd jxAs a transverse coordinate variable, itAs a longitudinal coordinate variable, sn(ix,it) And sd(ix,it) Respectively the numerator and denominator of the fault attribute expression, g (i)x,it) For use in calculating tomographic properties (i)x,it) Amplitude value of p (i)x,it) For seismic imaging in (i)x,it) The tangential slope of the structure, MxFor calculating lateral direction of fault attributeThe window is long.
Preferably, the structural constraint preconditioned chromatography equation set is:
STLTLSu＝STLTτ (2)
wherein, L is a ray chromatography kernel function, Δ m ═ Su, Δ m and Su are velocity model updating quantities, τ is the difference between the forward calculation seismic wave propagation time and the received data propagation time, S is a preconditioner, and S ═ I + DTGD)-1I is the unit matrix and D is the gradient operator
DTIs the transpose matrix of D
G is the structure tensor operator.
The invention has the beneficial effects that: according to the fault constraint tomography inversion method and system, the underground geological inclination angle and the fault attribute of the migration profile are calculated, the structural constraint preconditioned tomography equation set is constructed based on the underground geological inclination angle and the fault attribute, the structural constraint preconditioned tomography equation set is solved, and the underground speed model and the migration profile are obtained, so that the automation of the whole tomography inversion iteration process is realized, manual intervention is not needed, the manual workload is reduced, the processing period is shortened, and the working efficiency is improved; the velocity model is constrained by introducing underground geological dip angle and fault attributes, so that the velocity model obtained by inversion has higher precision in a complex structural area, particularly a fault development area, and meets the requirement of subsequent offset imaging.
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 throughout.
FIG. 1 shows a flow diagram of a fault-confined tomographic inversion method according to an embodiment of the invention.
FIG. 2 shows a raw seismic imaging profile of a fault-confined tomographic inversion method according to one embodiment of the present invention.
FIG. 3 shows a fault attribute map for a fault-constrained tomographic inversion method according to an embodiment of the invention.
FIG. 4 shows a seismic imaging profile of a fault-confined tomographic inversion method according to an embodiment of the present invention.
FIG. 5a shows a high steep profile seismic image without consideration of fault-constrained velocity modeling for a fault-constrained tomographic inversion method according to an embodiment of the invention.
FIG. 5b shows a high steep profile seismic image modeled with consideration of fault constraint velocities for a fault constraint tomographic inversion method according to an embodiment of the invention.
FIG. 5c shows a high steep fault-bound velocity model diagram of a fault-bound tomographic inversion method according to an embodiment of the invention.
FIG. 6a shows a small-scale geological intrusion seismic imaging plot without consideration of fault-constrained velocity modeling for a fault-constrained tomographic inversion method according to an embodiment of the invention.
FIG. 6b shows a small-scale geological intrusion seismic imaging plot for fault-constrained velocity modeling with consideration of a fault-constrained tomographic inversion method according to an embodiment of the invention.
FIG. 6c shows a small-scale geological invaded body fault-constrained velocity model plot of a fault-constrained tomographic inversion method according to one embodiment of the present invention.
Detailed Description
Preferred embodiments of the present invention will be described in more detail below. While the following describes preferred embodiments of the present invention, it should be understood that the present invention may be embodied in various forms and should not be limited by 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.
The fault constraint tomography inversion method comprises the following steps: step 1: calculating subsurface geological dip and fault properties of the migration profile based on the seismic migration image; step 2: constructing a structural constraint preconditioned chromatographic equation set based on the underground geological inclination and fault attributes; and step 3: and solving the structural constraint preconditioned chromatographic equation set to obtain an underground velocity model and an offset profile.
Specifically, calculating an underground geological inclination angle and fault attributes of the migration profile, constructing a structural constraint preconditioned chromatographic equation set based on the underground geological inclination angle and the fault attributes, and solving the structural constraint preconditioned chromatographic equation set to obtain an underground speed model and the migration profile.
According to the fault constraint tomography inversion method of the exemplary embodiment, the automation of the whole tomography inversion iteration process is realized, manual intervention is not needed, the manual workload is reduced, the processing period is shortened, and the working efficiency is improved; the velocity model is constrained by introducing underground geological dip angle and fault attributes, so that the velocity model obtained by inversion has higher precision in a complex structural area, particularly a fault development area, and meets the requirement of subsequent offset imaging.
Preferably, the subsurface geologic dip of the offset section is calculated by a structure tensor algorithm.
Specifically, the local underground geological inclination angle information of the offset profile is calculated by adopting a structure tensor algorithm. And setting I as a two-dimensional seismic image, wherein a structure tensor representing spatial direction information in the two-dimensional image I is defined by an image gradient value, the structure tensor represents the change direction of a region and the variation along the change direction, and seismic stratum textures and fault textures are determined by the variation relation of azimuth information of local points. The introduction of Gaussian functions (Gaussian) blurs the local detail, so that the structure tensor highlights the complexity of the signal in the region. For two-dimensional images, the structure tensor is a 2 x 2 matrix:
wherein G is the structure tensor, GxAnd gyThe gradients of the seismic image in the horizontal and vertical directions,<and- > is two-dimensional Gaussian smooth filtering.
For a semi-positive definite matrix G, eigenvalues and eigenvectors may be obtained by solving | G- λ I | ═ 0:
wherein v is1Is the first feature vector, v2Is a second eigenvector, λ1For maximum eigenvalues, the tensor energy is in the first eigentensor direction v1Energy of (a), λ2For minimum eigenvalues, the tensor energy is in the second eigentensor direction v2Energy of (a)1-λ2)/λ1And the linearity reflects the consistency of local directions.
The geological dip angle information is hidden in the structure tensor G, the characteristic vector describes the directionality of the local linear structure of the image, and the characteristic vector v is specific to each point of the image1Normal to the main structural direction of the image, eigenvector v2Parallel to the main structure direction of the image, and the local inclination angle direction of the point is a vector v2In the direction of (a). The structure tensor operator G ideally contains local structural features of the subsurface geologic structure.
Preferably, the fault attribute is calculated according to the following formula:
wherein, s (i)x,it) As a fault attribute, ixAnd jxAs a transverse coordinate variable, itAs a longitudinal coordinate variable, sn(ix,it) And sd(ix,it) Respectively the numerator and denominator of the fault attribute expression, g (i)x,it) For use in calculating tomographic properties (i)x,it) Amplitude value of p (i)x,it) For seismic imaging in (i)x,it) The tangential slope of the structure, MxThe lateral window length of the fault attribute is calculated.
Specifically, the fault attribute is calculated by a method of calculating the similarity coefficient along the construction direction, and the specific formula is as follows:
wherein, s (i)x,it) As a fault attribute, ixAnd jxAs a transverse coordinate variable, itAs a longitudinal coordinate variable, sn(ix,it) And sd(ix,it) Respectively the numerator and denominator of the fault attribute expression, g (i)x,it) For use in calculating tomographic properties (i)x,it) Amplitude value of p (i)x,it) For seismic imaging in (i)x,it) The tangential slope of the structure, MxThe lateral window length of the fault attribute is calculated.
As a preferred scheme, the structural constraint preconditioned chromatographic equation set is as follows:
STLTLSu＝STLTτ (2)
wherein, L is a ray chromatography kernel function, Δ m ═ Su, Δ m and Su are velocity model updating quantities, τ is the difference between the forward calculation seismic wave propagation travel time and the received data travel time, and S is a preconditioner.
Preferably, the preconditioner S is a smoothing operator containing geological structure information, S being:
S＝(I+DTGD)-1(3)
where I is the unit matrix and D is the gradient operator
DTIs the transpose matrix of D
G is the structure tensor operator.
Specifically, the conventional chromatography equation set can be expressed as follows:
LΔm＝τ (6)
wherein, L is a ray chromatography kernel function, which is introduced in the field of seismic chromatography inversion, and the invention is not developed in this way. Δ m is the velocity model update quantity, τ is the difference between the forward computed seismic wave travel time and the received data travel time.
Considering the model preconditions, i.e., Δ m — Su, the construction of the system of constrained preconditioned tomographic equations can be expressed as:
STLTLSu＝STLTτ (2)
the pre-condition operator S is a smooth operator containing geological structure information, the equation is constructed as a chromatographic equation of geological structure constraint pre-conditions, and the corresponding solution is a smooth solution after the pre-conditions.
After the model parameterization of the underground medium, the basic geological rule of the underground medium is not changed, so that certain relation necessarily exists between the parameters. The accuracy of data measurement in tomography is determined by the inclination of the reflecting surface on the offset profile and the accuracy of the co-imaging to the depth of the concentrated pick-up image. Therefore, the inclination angle information of the reflecting surface and the distribution rule of scattering points in depth in the geological structure provide a feasible way for geological smoothing, and the prior information is not relied on. Adding the geological structure information into a constructed preconditioner S, wherein the preconditioner S is expressed as follows:
S＝(I+DTGD)-1(3)
where I is the unit matrix and D is the gradient operator
DTIs the transpose matrix of D
G is the structure tensor operator.
The structure tensor operator G implies geological dip angle information, and the eigenvalue lambda of the structure tensor operator G can be obtained by decomposing the eigenvalue of the operator G1And λ2Setting small positive real number in the area with fault attribute satisfying formula (7) to ensure that no chromatographic smoothing is applied in the fault development area to maintain the accuracy of chromatographic inversion at the fault boundary,
s(ix,it)＜σ (7)
wherein, σ is a threshold value for judging whether a fault exists at a certain point on the seismic imaging section, is a positive real number distributed from 0 to 1, and can be adjusted according to the signal-to-noise ratio of the seismic imaging section. Generally speaking, for the case that the signal-to-noise ratio is relatively high and the fault distance is relatively obvious, the threshold value is usually selected to be larger to ensure the accuracy, otherwise, the threshold value is relatively smaller to ensure the stability.
As a preferred scheme, the structural constraint preconditioned chromatographic equation set is solved by a least square QR decomposition method.
Specifically, the structural constraint preconditioned tomographic equation set is solved by a least square QR decomposition method, and the method is an iterative method and can efficiently solve a large-scale sparse matrix in the least square sense.
The fault constraint tomography inversion system comprises: a memory storing computer-executable instructions; a processor executing computer executable instructions in the memory to perform the steps of: step 1: calculating subsurface geological dip and fault properties of the migration profile based on the seismic migration image; step 2: constructing a structural constraint preconditioned chromatographic equation set based on the underground geological inclination and fault attributes; and step 3: and solving the structural constraint preconditioned chromatographic equation set to obtain an underground velocity model and an offset profile.
Specifically, calculating an underground geological inclination angle and fault attributes of the migration profile, constructing a structural constraint preconditioned chromatographic equation set based on the underground geological inclination angle and the fault attributes, and solving the structural constraint preconditioned chromatographic equation set to obtain an underground speed model and the migration profile.
Preferably, the subsurface geologic dip of the offset section is calculated by a structure tensor algorithm.
Specifically, the local underground geological inclination angle information of the offset profile is calculated by adopting a structure tensor algorithm. And setting I as a two-dimensional seismic image, wherein a structure tensor representing spatial direction information in the two-dimensional image I is defined by an image gradient value, the structure tensor represents the change direction of a region and the variation along the change direction, and seismic stratum textures and fault textures are determined by the variation relation of azimuth information of local points. The introduction of Gaussian functions (Gaussian) blurs the local detail, so that the structure tensor highlights the complexity of the signal in the region. For two-dimensional images, the structure tensor is a 2 x 2 matrix:
wherein G is the structure tensor, GxAnd gyThe gradients of the seismic image in the horizontal and vertical directions,<and- > is two-dimensional Gaussian smooth filtering.
For a semi-positive definite matrix G, eigenvalues and eigenvectors may be obtained by solving | G- λ I | ═ 0:
wherein v is1Is the first feature vector, v2Is a second eigenvector, λ1For maximum eigenvalues, the tensor energy is in the first eigentensor direction v1Energy of (a), λ2For minimum eigenvalues, the tensor energy is in the second eigentensor direction v2Energy of (a)1-λ2)/λ1And the linearity reflects the consistency of local directions.
The geological dip angle information is hidden in the structure tensor G, the characteristic vector describes the directionality of the local linear structure of the image, and the characteristic vector v is specific to each point of the image1Normal to the main structural direction of the image, eigenvector v2Parallel to the main structural direction of the image, the local of the pointThe direction of the inclination angle is the vector v2In the direction of (a). The structure tensor operator G ideally contains local structural features of the subsurface geologic structure.
Preferably, the fault attribute is calculated according to the following formula:
wherein, s (i)x,it) As a fault attribute, ixAnd jxAs a transverse coordinate variable, itAs a longitudinal coordinate variable, sn(ix,it) And sd(ix,it) Respectively the numerator and denominator of the fault attribute expression, g (i)x,it) For use in calculating tomographic properties (i)x,it) Amplitude value of p (i)x,it) For seismic imaging in (i)x,it) The tangential slope of the structure, MxThe lateral window length of the fault attribute is calculated.
Specifically, the fault attribute is calculated by a method of calculating the similarity coefficient along the construction direction, and the specific formula is as follows:
wherein, s (i)x,it) As a fault attribute, ixAnd jxAs a transverse coordinate variable, itAs a longitudinal coordinate variable, sn(ix,it) And sd(ix,it) Respectively the numerator and denominator of the fault attribute expression, g (i)x,it) For use in calculating tomographic properties (i)x,it) Amplitude value of p (i)x,it) For seismic imaging in (i)x,it) The tangential slope of the structure, MxThe lateral window length of the fault attribute is calculated.
As a preferred scheme, the structural constraint preconditioned chromatographic equation set is as follows:
STLTLSu＝STLTτ (2)
wherein, L is a ray chromatography kernel function, Δ m ═ Su, Δ m and Su are velocity model updating quantities, τ is the difference between the forward calculation seismic wave propagation travel time and the received data travel time, and S is a preconditioner.
Preferably, the preconditioner S is a smoothing operator containing geological structure information, S being:
S＝(I+DTGD)-1(3)
where I is the unit matrix and D is the gradient operator
DTIs the transpose matrix of D
G is the structure tensor operator.
Specifically, the conventional chromatography equation set can be expressed as follows:
LΔm＝τ (6)
wherein, L is a ray chromatography kernel function, which is introduced in the field of seismic chromatography inversion, and the invention is not developed in this way. Δ m is the velocity model update quantity, τ is the difference between the forward computed seismic wave travel time and the received data travel time.
Considering the model preconditions, i.e., Δ m — Su, the construction of the system of constrained preconditioned tomographic equations can be expressed as:
STLTLSu＝STLTτ (2)
the pre-condition operator S is a smooth operator containing geological structure information, the equation is constructed as a chromatographic equation of geological structure constraint pre-conditions, and the corresponding solution is a smooth solution after the pre-conditions.
After the model parameterization of the underground medium, the basic geological rule of the underground medium is not changed, so that certain relation necessarily exists between the parameters. The accuracy of data measurement in tomography is determined by the inclination of the reflecting surface on the offset profile and the accuracy of the co-imaging to the depth of the concentrated pick-up image. Therefore, the inclination angle information of the reflecting surface and the distribution rule of scattering points in depth in the geological structure provide a feasible way for geological smoothing, and the prior information is not relied on. Adding the geological structure information into a constructed preconditioner S, wherein the preconditioner S is expressed as follows:
S＝(I+DTGD)-1(3)
where I is the unit matrix and D is the gradient operator
DTIs the transpose matrix of D
G is the structure tensor operator.
The structure tensor operator G implies geological dip angle information, and the eigenvalue lambda of the structure tensor operator G can be obtained by decomposing the eigenvalue of the operator G1And λ2Setting small positive real number in the area with fault attribute satisfying formula (7) to ensure that no chromatographic smoothing is applied in the fault development area to maintain the accuracy of chromatographic inversion at the fault boundary,
s(ix,it)＜σ (7)
wherein, σ is a threshold value for judging whether a fault exists at a certain point on the seismic imaging section, is a positive real number distributed from 0 to 1, and can be adjusted according to the signal-to-noise ratio of the seismic imaging section. Generally speaking, for the case that the signal-to-noise ratio is relatively high and the fault distance is relatively obvious, the threshold value is usually selected to be larger to ensure the accuracy, otherwise, the threshold value is relatively smaller to ensure the stability.
As a preferred scheme, the structural constraint preconditioned chromatographic equation set is solved by a least square QR decomposition method.
Specifically, the structural constraint preconditioned tomographic equation set is solved by a least square QR decomposition method, and the method is an iterative method and can efficiently solve a large-scale sparse matrix in the least square sense.
Examples
FIG. 1 shows a flow diagram of a fault-confined tomographic inversion method according to an embodiment of the invention. FIG. 2 shows an original seismic imaging section of a fault-confined tomographic inversion method according to an embodiment of the present invention. FIG. 3 shows a fault attribute map for a fault-constrained tomographic inversion method according to an embodiment of the invention. FIG. 4 shows a seismic imaging profile of a fault-confined tomographic inversion method according to an embodiment of the present invention. FIG. 5a shows a high steep profile seismic image without consideration of fault-constrained velocity modeling for a fault-constrained tomographic inversion method according to an embodiment of the invention. FIG. 5b shows a high steep profile seismic image modeled with consideration of fault constraint velocities for a fault constraint tomographic inversion method according to an embodiment of the invention. FIG. 5c shows a high steep fault-bound velocity model diagram of a fault-bound tomographic inversion method according to an embodiment of the invention. FIG. 6a shows a small-scale geological intrusion seismic imaging plot without consideration of fault-constrained velocity modeling for a fault-constrained tomographic inversion method according to an embodiment of the invention. FIG. 6b shows a small-scale geological intrusion seismic imaging plot for fault-constrained velocity modeling with consideration of a fault-constrained tomographic inversion method according to an embodiment of the invention. FIG. 6c shows a small-scale geological invaded body fault-constrained velocity model plot of a fault-constrained tomographic inversion method according to one embodiment of the present invention.
As shown in fig. 1, the fault-constrained tomography inversion method includes:
step 1: calculating subsurface geological dip and fault properties of the migration profile based on the seismic migration image;
specifically, the local underground geological inclination angle information of the offset profile is calculated by adopting a structure tensor algorithm. And setting I as a two-dimensional seismic image, wherein a structure tensor representing spatial direction information in the two-dimensional image I is defined by an image gradient value, the structure tensor represents the change direction of a region and the variation along the change direction, and seismic stratum textures and fault textures are determined by the variation relation of azimuth information of local points. The introduction of Gaussian functions (Gaussian) blurs the local detail, so that the structure tensor highlights the complexity of the signal in the region. For two-dimensional images, the structure tensor is a 2 x 2 matrix:
wherein G is the structure tensor, GxAnd gyThe gradients of the seismic image in the horizontal and vertical directions,<and- > is two-dimensional Gaussian smooth filtering.
For a semi-positive definite matrix G, eigenvalues and eigenvectors may be obtained by solving | G- λ I | ═ 0:
wherein v is1Is the first feature vector, v2Is a second eigenvector, λ1For maximum eigenvalues, the tensor energy is in the first eigentensor direction v1Energy of (a), λ2For minimum eigenvalues, the tensor energy is in the second eigentensor direction v2Energy of (a)1-λ2)/λ1And the linearity reflects the consistency of local directions.
The geological dip angle information is hidden in the structure tensor G, the characteristic vector describes the directionality of the local linear structure of the image, and the characteristic vector v is specific to each point of the image1Normal to the main structural direction of the image, eigenvector v2Parallel to the main structure direction of the image, and the local inclination angle direction of the point is a vector v2In the direction of (a). The structure tensor operator G ideally contains local structural features of the subsurface geologic structure.
The fault attribute is calculated by adopting a method for calculating the similarity coefficient along the construction direction, and the specific formula is as follows:
wherein, s (i)x,it) As a fault attribute, ixAnd jxAs a transverse coordinate variable, itAs a longitudinal coordinate variable, sn(ix,it) And sd(ix,it) Respectively the numerator and denominator of the fault attribute expression, g (i)x,it) For use in calculating tomographic properties (i)x,it) Amplitude value of p (i)x,it) For seismic imaging in (i)x,it) The tangential slope of the structure, MxThe lateral window length of the fault attribute is calculated.
Step 2: constructing a structural constraint preconditioned chromatographic equation set based on the underground geological inclination and fault attributes;
specifically, the conventional chromatography equation set can be expressed as follows:
LΔm＝τ (6)
wherein, L is a ray chromatography kernel function, which is introduced in the field of seismic chromatography inversion, and the invention is not developed in this way. Δ m is the velocity model update quantity, τ is the difference between the forward computed seismic wave travel time and the received data travel time.
Considering the model preconditions, i.e., Δ m — Su, the construction of the system of constrained preconditioned tomographic equations can be expressed as:
STLTLSu＝STLTτ (2)
the pre-condition operator S is a smooth operator containing geological structure information, the equation is constructed as a chromatographic equation of geological structure constraint pre-conditions, and the corresponding solution is a smooth solution after the pre-conditions.
After the model parameterization of the underground medium, the basic geological rule of the underground medium is not changed, so that certain relation necessarily exists between the parameters. The accuracy of data measurement in tomography is determined by the inclination of the reflecting surface on the offset profile and the accuracy of the co-imaging to the depth of the concentrated pick-up image. Therefore, the inclination angle information of the reflecting surface and the distribution rule of scattering points in depth in the geological structure provide a feasible way for geological smoothing, and the prior information is not relied on. Adding the geological structure information into a constructed preconditioner S, wherein the preconditioner S is expressed as follows:
S＝(I+DTGD)-1(3)
where I is the unit matrix and D is the gradient operator
DTIs the transpose matrix of D
G is the structure tensor operator.
The structure tensor operator G implies geological dip angle information, and the eigenvalue lambda of the structure tensor operator G can be obtained by decomposing the eigenvalue of the operator G1And λ2Setting small positive real number in the area with fault attribute satisfying formula (7) to ensure that no chromatographic smoothing is applied in the fault development area to maintain the accuracy of chromatographic inversion at the fault boundary,
s(ix,it)＜σ (7)
wherein, σ is a threshold value for judging whether a fault exists at a certain point on the seismic imaging section, is a positive real number distributed from 0 to 1, and can be adjusted according to the signal-to-noise ratio of the seismic imaging section. Generally speaking, for the case that the signal-to-noise ratio is relatively high and the fault distance is relatively obvious, the threshold value is usually selected to be larger to ensure the accuracy, otherwise, the threshold value is relatively smaller to ensure the stability.
And step 3: solving a structural constraint preconditioned chromatographic equation set to obtain an underground speed model and an offset profile;
wherein, the structural constraint preconditioned chromatographic equation set is solved by a least square QR decomposition method.
FIG. 2 shows a raw seismic imaging profile of a fault-confined tomographic inversion method according to one embodiment of the present invention.
As shown in fig. 2, the original seismic imaging section is used to calculate the structure tensor and extract the fault attribute information, which is used as the basis for subsequently constructing the subsequent tomographic preconditioner.
FIG. 3 shows a fault attribute map for a fault-constrained tomographic inversion method according to an embodiment of the invention.
As shown in fig. 3, a computed fault attribute map is obtained by taking an original seismic imaging section as an input, and it can be seen that the fault attribute has higher correlation with a fault development area of the seismic section shown in fig. 2, and shows better stability and adaptability to noise. The fault attributes are distributed between 0 and 1, with the greater the attribute value the greater the likelihood of developing a fault.
FIG. 4 shows a seismic imaging profile of a fault-confined tomographic inversion method according to an embodiment of the present invention.
As shown in FIG. 4, the tomography preconditioner is constructed by utilizing the fault extraction information, and the preconditioner is applied to the original seismic imaging section, so that the processed image has higher continuity as a whole, the section is effectively strengthened, and the tomography preconditioner has higher precision in the fault development area.
FIG. 5a shows a high steep profile seismic image without consideration of fault-constrained velocity modeling for a fault-constrained tomographic inversion method according to an embodiment of the invention. FIG. 5b shows a high steep profile seismic image modeled with consideration of fault constraint velocities for a fault constraint tomographic inversion method according to an embodiment of the invention. FIG. 5c shows a high steep fault-bound velocity model diagram of a fault-bound tomographic inversion method according to an embodiment of the invention.
As shown in fig. 5a, 5b and 5c, fig. 5a and 5b are seismic imaging section diagrams corresponding to modeling considering fault constraint velocity without considering fault approximation, respectively, and it can be seen from the diagrams that the seismic imaging section considering fault constraint velocity modeling shows higher imaging signal-to-noise ratio, and fig. 5c shows that the velocity model obtained by inversion of the invention is consistent with the seismic fault characteristics in a high steep section development area, and shows very high inversion resolution.
FIG. 6a shows a small-scale geological intrusion seismic imaging plot without consideration of fault-constrained velocity modeling for a fault-constrained tomographic inversion method according to an embodiment of the invention. FIG. 6b shows a small-scale geological intrusion seismic imaging plot for fault-constrained velocity modeling with consideration of a fault-constrained tomographic inversion method according to an embodiment of the invention. FIG. 6c shows a small-scale geological invaded body fault-constrained velocity model plot of a fault-constrained tomographic inversion method according to one embodiment of the present invention.
As shown in fig. 6a, 6b and 6c, the quality of the seismic imaging profile modeled by considering fault constraint velocity in fig. 6b is obviously better than that of the seismic imaging profile modeled by considering fault constraint velocity in fig. 6a, and the velocity model in fig. 6c well recovers a small-scale geological invasion body, corresponding to the seismic imaging profile.
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 (10)
1. A fault-constrained tomographic inversion method, comprising:
step 1: calculating subsurface geological dip and fault properties of the migration profile based on the seismic migration image;
step 2: constructing a structural constraint preconditioned chromatographic equation set based on the underground geological inclination angle and the fault attribute;
and step 3: and solving the structural constraint preconditioned chromatographic equation set to obtain an underground velocity model and an offset profile.
2. The fault-constrained tomographic inversion method of claim 1, wherein the subsurface geological dip of the offset section is calculated by a structure tensor algorithm.
3. The fault-constrained tomographic inversion method of claim 1, wherein the fault attribute is calculated according to the following formula:
wherein, s (i)x,it) As a fault attribute, ixAnd jxAs a transverse coordinate variable, itAs a longitudinal coordinate variable，sn(ix,it) And sd(ix,it) Respectively the numerator and denominator of the fault attribute expression, g (i)x,it) For use in calculating tomographic properties (i)x,it) Amplitude value of p (i)x,it) For seismic imaging in (i)x,it) The tangential slope of the structure, MxThe lateral window length of the fault attribute is calculated.
4. The fault-constrained tomographic inversion method of claim 1, wherein the structural-constrained preconditioned tomographic equations are:
STLTLSu＝STLTτ (2)
wherein, L is a ray chromatography kernel function, Δ m ═ Su, Δ m and Su are velocity model updating quantities, τ is the difference between the forward calculation seismic wave propagation travel time and the received data travel time, and S is a preconditioner.
5. The fault-constrained tomographic inversion method of claim 4, wherein the preconditioner S is a smoothing operator containing geological structure information, S being:
S＝(I+DTGD)-1(3)
where I is the unit matrix and D is the gradient operator
DTIs the transpose matrix of D
G is the structure tensor operator.
6. The fault-constrained tomographic inversion method of claim 1, wherein the set of construction-constrained preconditioned tomographic equations is solved by a least-squares (QR) decomposition method.
7. A fault-constrained tomographic inversion system, comprising: a memory storing computer-executable instructions;
a processor executing computer executable instructions in the memory to perform the steps of:
step 1: calculating subsurface geological dip and fault properties of the migration profile based on the seismic migration image;
step 2: constructing a structural constraint preconditioned chromatographic equation set based on the underground geological inclination angle and the fault attribute;
and step 3: and solving the structural constraint preconditioned chromatographic equation set to obtain an underground velocity model and an offset profile.
8. The fault-constrained tomographic inversion system of claim 7, wherein subsurface geological dips of the offset section are calculated by a structure tensor algorithm.
9. The fault-constrained tomographic inversion system of claim 7, wherein the fault attribute is calculated according to the following formula:
wherein, s (i)x,it) As a fault attribute, ixAnd jxAs a transverse coordinate variable, itAs a longitudinal coordinate variable, sn(ix,it) And sd(ix,it) Respectively the numerator and denominator of the fault attribute expression, g (i)x,it) For use in calculating tomographic properties (i)x,it) Amplitude value of p (i)x,it) For seismic imaging in (i)x,it) The tangential slope of the structure, MxThe lateral window length of the fault attribute is calculated.
10. The fault-constrained tomographic inversion system of claim 7, wherein the construction-constrained preconditioned tomographic equations are:
STLTLSu＝STLTτ (2)
wherein, L is a ray chromatography kernel function, Δ m ═ Su, Δ m and Su are velocity model updating quantities, τ is the difference between the forward calculation seismic wave propagation time and the received data propagation time, S is a preconditioner, and S ═ I + DTGD)-1I is the unit matrix and D is the gradient operator
DTIs the transpose matrix of D
G is the structure tensor operator.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811094492.XA CN110927779B (en) | 2018-09-19 | 2018-09-19 | Fault constraint tomography inversion method and inversion system |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811094492.XA CN110927779B (en) | 2018-09-19 | 2018-09-19 | Fault constraint tomography inversion method and inversion system |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110927779A true CN110927779A (en) | 2020-03-27 |
| CN110927779B CN110927779B (en) | 2021-09-17 |
Family
ID=69856043
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811094492.XA Active CN110927779B (en) | 2018-09-19 | 2018-09-19 | Fault constraint tomography inversion method and inversion system |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110927779B (en) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111830561A (en) * | 2020-07-24 | 2020-10-27 | 中国科学技术大学 | Method for constructing fault three-dimensional structure based on seismic distribution characteristics |
| CN112083144A (en) * | 2020-09-01 | 2020-12-15 | 中国科学院地质与地球物理研究所 | Fault on-off prediction method and device, computer equipment and storage medium |
| CN113589385A (en) * | 2021-08-11 | 2021-11-02 | 成都理工大学 | Reservoir characteristic inversion method based on seismic scattering wave field analysis |
| CN113589375A (en) * | 2020-04-30 | 2021-11-02 | 中国石油化工股份有限公司 | VSP layer velocity inversion method based on inclined layer constraint travel time calculation |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20070032955A1 (en) * | 2005-08-08 | 2007-02-08 | Williams Kenneth E | Method and system for combining seismic data and basin modeling |
| US20110267921A1 (en) * | 2010-04-30 | 2011-11-03 | Schlumberger Technology Corporation | Multicomponent seismic inversion of vsp data |
| CN103135136A (en) * | 2011-11-25 | 2013-06-05 | 中国石油化工股份有限公司 | Automatic fault interpretation device for three-dimensional seismic data body |
| CN103514630A (en) * | 2013-10-16 | 2014-01-15 | 北京石油化工学院 | Fault structure three-dimensional modeling method |
| CN108072892A (en) * | 2016-11-09 | 2018-05-25 | 中国石油化工股份有限公司 | A kind of geological structure constraint chromatography conversion method of automation |
-
2018
- 2018-09-19 CN CN201811094492.XA patent/CN110927779B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20070032955A1 (en) * | 2005-08-08 | 2007-02-08 | Williams Kenneth E | Method and system for combining seismic data and basin modeling |
| US20110267921A1 (en) * | 2010-04-30 | 2011-11-03 | Schlumberger Technology Corporation | Multicomponent seismic inversion of vsp data |
| GB2490278A (en) * | 2010-04-30 | 2012-10-24 | Schlumberger Holdings | Multicomponent seismic inversion of VSP data |
| CN103135136A (en) * | 2011-11-25 | 2013-06-05 | 中国石油化工股份有限公司 | Automatic fault interpretation device for three-dimensional seismic data body |
| CN103514630A (en) * | 2013-10-16 | 2014-01-15 | 北京石油化工学院 | Fault structure three-dimensional modeling method |
| CN108072892A (en) * | 2016-11-09 | 2018-05-25 | 中国石油化工股份有限公司 | A kind of geological structure constraint chromatography conversion method of automation |
Non-Patent Citations (1)
| Title |
|---|
| 彭海龙等: "基于地质构造约束的3D速度建模方法在琼东南盆地深水复杂断块区域成像中的应用", 《物探与化探》 * |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113589375A (en) * | 2020-04-30 | 2021-11-02 | 中国石油化工股份有限公司 | VSP layer velocity inversion method based on inclined layer constraint travel time calculation |
| CN111830561A (en) * | 2020-07-24 | 2020-10-27 | 中国科学技术大学 | Method for constructing fault three-dimensional structure based on seismic distribution characteristics |
| CN111830561B (en) * | 2020-07-24 | 2022-09-06 | 中国科学技术大学 | Method for constructing fault three-dimensional structure based on seismic distribution characteristics |
| CN112083144A (en) * | 2020-09-01 | 2020-12-15 | 中国科学院地质与地球物理研究所 | Fault on-off prediction method and device, computer equipment and storage medium |
| CN113589385A (en) * | 2021-08-11 | 2021-11-02 | 成都理工大学 | Reservoir characteristic inversion method based on seismic scattering wave field analysis |
| CN113589385B (en) * | 2021-08-11 | 2023-08-04 | 成都理工大学 | Reservoir characteristic inversion method based on seismic scattered wave field analysis |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110927779B (en) | 2021-09-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110927779B (en) | Fault constraint tomography inversion method and inversion system | |
| CN108072892B (en) | Automatic geological structure constraint chromatography inversion method | |
| CA2683618C (en) | Inverse-vector method for smoothing dips and azimuths | |
| CN107688201B (en) | RBM-based seismic prestack signal clustering method | |
| US20110115787A1 (en) | Visulation of geologic features using data representations thereof | |
| US20120004849A1 (en) | Efficient windowed radon transform | |
| US9022129B2 (en) | Tracking geologic object and detecting geologic anomalies in exploration seismic data volume | |
| CN111290021A (en) | Carbonate rock fracture-cave enhanced identification method based on gradient structure tensor spectrum decomposition | |
| US11733413B2 (en) | Method and system for super resolution least-squares reverse time migration | |
| AU2013200609A1 (en) | Visulation of geologic features using data representations thereof | |
| Li et al. | Extended full waveform inversion with matching filter | |
| Zhou et al. | Stochastic structure-constrained image-guided inversion of geophysical data | |
| CN112731520A (en) | Full waveform inversion method and system based on structure tensor diffusion filtering | |
| US11965999B2 (en) | Method and system for processing seismic images to obtain seismic horizon surfaces for a geological formation | |
| CN109655910A (en) | The two-parameter full waveform inversion method of Ground Penetrating Radar based on phasing | |
| Cheng et al. | Multiscale fracture prediction technique via deep learning, seismic gradient disorder, and aberrance: Applied to tight sandstone reservoirs in the Hutubi block, southern Junggar Basin | |
| Cullison et al. | An image-guided method for automatically picking common-image-point gathers | |
| Zhang | Ensemble methods of data assimilation in porous media flow for non-Gaussian prior probability density | |
| Zand et al. | Least-squares reverse time migration with shifted total variation regularization | |
| Chen et al. | Augmented deep learning workflow for robust fault prediction over multiple tectonic regimes | |
| CN110927780B (en) | Geological stratum position constrained small-scale geologic body velocity modeling method and system | |
| US20240111067A1 (en) | Faulted seismic horizon mapping | |
| US11630226B2 (en) | Formation evaluation based on seismic horizon mapping with multi-scale optimization | |
| US20240052734A1 (en) | Machine learning framework for sweep efficiency quantification | |
| US20240069236A1 (en) | Method and system for processing seismic images to obtain an rgt image by using a decimated seismic dip image |
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 |