CN115097520A - Frequency-dependent anisotropic shear wave splitting algorithm for multiple VSP data - Google Patents
Frequency-dependent anisotropic shear wave splitting algorithm for multiple VSP data Download PDFInfo
- Publication number
- CN115097520A CN115097520A CN202210650897.7A CN202210650897A CN115097520A CN 115097520 A CN115097520 A CN 115097520A CN 202210650897 A CN202210650897 A CN 202210650897A CN 115097520 A CN115097520 A CN 115097520A
- Authority
- CN
- China
- Prior art keywords
- frequency
- algorithm
- dependent
- transverse wave
- axis
- 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.)
- Pending
Links
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 45
- 230000001419 dependent effect Effects 0.000 title claims abstract description 44
- 239000013598 vector Substances 0.000 claims abstract description 29
- 238000000034 method Methods 0.000 claims abstract description 12
- 230000000694 effects Effects 0.000 claims abstract description 11
- 238000012545 processing Methods 0.000 claims abstract description 8
- 238000009795 derivation Methods 0.000 claims abstract description 5
- 238000007405 data analysis Methods 0.000 claims abstract description 4
- 238000004088 simulation Methods 0.000 claims abstract description 4
- 239000011159 matrix material Substances 0.000 claims description 26
- 239000012530 fluid Substances 0.000 claims description 12
- 239000011435 rock Substances 0.000 claims description 12
- 230000010287 polarization Effects 0.000 claims description 10
- 238000004458 analytical method Methods 0.000 claims description 8
- 239000002245 particle Substances 0.000 claims description 8
- 230000014509 gene expression Effects 0.000 claims description 7
- 238000006073 displacement reaction Methods 0.000 claims description 6
- 239000011148 porous material Substances 0.000 claims description 6
- 238000012937 correction Methods 0.000 claims description 3
- 230000007246 mechanism Effects 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 230000001131 transforming effect Effects 0.000 claims description 3
- 230000017105 transposition Effects 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 claims 1
- 238000001914 filtration Methods 0.000 abstract description 3
- 230000009466 transformation Effects 0.000 description 3
- 238000009825 accumulation Methods 0.000 description 2
- 230000015572 biosynthetic process Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000009977 dual effect Effects 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 230000004075 alteration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012795 verification Methods 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
- G01V1/30—Analysis
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- 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
- G01V1/282—Application of seismic models, synthetic seismograms
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 frequency-dependent anisotropic transverse wave splitting algorithm for multiple VSP data, which comprises the following steps: s1, obtaining a synthetic record of frequency-dependent anisotropy by adopting a numerical simulation method; s2, the existence of frequency-dependent anisotropy is proved through data analysis and based on a crack model; s3, deducing a multi-component VSP crack parameter inversion algorithm containing frequency-dependent anisotropic characteristics through processing and operation based on a vector convolution model; s4, inversion derivation is carried out on a direction angle theta of the fast transverse wave and time difference delta tau of the fast transverse wave and the slow transverse wave is calculated; and S5, analyzing the example effect of the algorithm. The algorithm is based on the crack model, proves the existence of frequency-dependent anisotropy, and meanwhile, based on the vector convolution model, through proper processing and operation, the multi-component VSP crack parameter inversion algorithm containing the frequency-dependent anisotropy characteristic can be deduced, so that the anisotropy parameters can be directly calculated for different frequencies in a frequency domain conveniently, and errors possibly brought by segmented filtering are eliminated.
Description
Technical Field
The invention belongs to the technical field of transverse wave splitting, and particularly relates to a frequency-dependent anisotropic transverse wave splitting algorithm for multiple VSP data.
Background
Seismic anisotropy, which has been the main phenomenon of shear wave splitting, has been widely applied to fractured rocks since the 80's 20 th century, and is considered as one of important means for obtaining internal physical properties of rocks in a formation, which have been difficult to obtain before, by direct measurement. In shear wave seismic exploration, one common application is the determination of fast shear waves (qS) using multi-component data 1 ) And its polarization direction with the split slow transverse wave (qS) 2 ) The time difference between them. The two parameters can directly give out the characteristic direction and strength of anisotropy, and further possibly give out detailed information such as stress in the oil reservoir, geometric characteristics of cracks and the like. In order to determine these two parameters, algorithms have emerged in recent years, which can be broadly divided into accumulation algorithms that calculate the accumulated average response between source and detector, such as Alford rotation, linear transformation, dual source accumulation, etc., and local algorithms that calculate the local average response between two detectors, such as propagation matrix, dual source transformation, etc. In addition, the method is based on a delaminating method between the two, the algorithms are all based on the convolution equation of wave propagation in the anisotropic medium, and people can select different algorithms according to different geological conditions. However, these conventional shear wave splitting algorithms, such as Alford rotation, linear transformation, etc., are performed in the time domain.
At present, cracks (fracure) and cracks (crack) are considered synonyms in seismic exploration physics, and we cannot distinguish microcracks from macrocracks. While both micro-scale and macro-scale fractures may be considered the primary cause of reservoir anisotropy, the latter is useful for oil production because the flow of fluids in the reservoir is dominated by the large-scale fluid cells. Therefore, it is likely that the quantitative characterization of fractures in the formation from seismic data is key to predicting flow field mobility and flow patterns. The existing equivalent medium theory under the long-wavelength assumption is independent of frequency, and all the equivalent medium theories have a parameter, namely crack density, related to the degree of anisotropy although the precision of each model is different. Usually the fracture density can be estimated by the time difference of shear wave splitting, but it is not sensitive to the size of the fracture and cannot be used to draw specific conclusions about whether the anisotropy of the rock of interest in practical applications is caused by micro-cracks or macro-cracks. This may be one of the main reasons why shear wave splitting techniques have not been widely used in the industry to date.
To solve this problem, models that take into account the actual fracture dimensions, effects due to fluid flow, etc. must be employed, which typically cause frequency-dependent anisotropy. Indeed, some recent practical data processing has found evidence of anisotropy as a function of frequency: the time delay between the fast and slow shear waves decreases with increasing frequency. Theoretical models related to the fracture model are also developed, such as a connected fracture and constant-diameter porous model proposed by Hudson et al, Tod populates the model to an oriented fracture, a fracture model with a combination of micro and macro scales proposed by Chapman, and the like. Common to these models is that the fracture scale is contained in the model and is frequency dependent.
However, the research in this area is just started internationally, and further development and verification of test data are needed, especially, no relevant transverse wave splitting algorithm exists at present about frequency anisotropy. The method used today is to filter the data piecewise and then analyze it using conventional methods. Since the time difference generated by transverse wave splitting is usually in the order of milliseconds, the artificial error generated by the analysis method is possibly larger than the divergence of the time difference, so that the analysis of the data is misled, and therefore, a frequency-dependent anisotropic transverse wave splitting algorithm of a large amount of VSP data is provided.
Disclosure of Invention
The present invention aims to provide a frequency-dependent anisotropic shear wave splitting algorithm for a large amount of VSP data to solve the problems mentioned in the background art.
In order to achieve the purpose, the invention provides the following technical scheme: a frequency dependent anisotropic shear wave splitting algorithm for multiple VSP data, comprising the steps of:
s1, obtaining a synthetic record of frequency-dependent anisotropy by adopting a numerical simulation method;
s2, the existence of frequency-dependent anisotropy is proved through data analysis based on a crack model;
s3, deducing a multi-component VSP crack parameter inversion algorithm containing frequency-dependent anisotropic characteristics through processing and operation based on a vector convolution model;
s4, inversion derivation is carried out on a direction angle theta of the fast transverse wave and time difference delta tau of the fast transverse wave and the slow transverse wave is calculated;
and S5, analyzing the example effect of the algorithm.
Preferably, the fracture model in S2 is used to theoretically illustrate that anisotropy is frequency dependent, and the fracture model has two different sizes, micro size and macro size: particle size and fracture size, the fluid in the pore space flows following a jet mechanism during propagation of the elastic wave, and the equivalent elastic stiffness tensor in porous rocks containing fractures is:
C＝C 0 -Φ p C 1 -ε c C 2 -ε f C 3 (1)
wherein, C 0 Is the elastic constant of isotropic surrounding rock, C 1 、C 2 And C 3 Correction terms generated when the pores, the micro-cracks and the cracks act independently are respectively; phi p Is the porosity of the surrounding rock, epsilon c Is the density of microcracks,. epsilon. f Is the crack density.
Preferably, the particle size of the fluid flow and relaxation time τ m In connection with the inflow and outflow of fluid in the fracture with a lower characteristic frequency or a longer relaxation time tau f Related, the correlation formula is as follows:
wherein af is the fracture radius, a m Is the particle size.
Preferably, the formula for the far-field vector wavefield in isotropic media generated by a point source in the vector convolution model in S3 is as follows:
d(t)＝Aδ(t-τ)*B(t)*s(t)/r (3)
wherein the seismic source is in a homogeneous medium and the vector s (t) represents the direction of polarization; τ represents the travel time (r/V) of the source wave with origin at the source s (t); r represents the total propagation distance, V represents the wave velocity; b (t) represents an attenuation factor dependent on the quality factor Q; a represents a constant amplitude factor, and represents convolution in the time domain;
the polarization vector of the seismic source is:
wherein s is X (t) and s Y (t) projections of the seismic source in the X-axis and Y-axis directions, respectively;
vectors s for equivalent motion produced by the seismic source in the directions of the f-axis and the s-axis, orthogonal according to the coordinate system p (t) represents, i.e.:
s p (t)＝R(θ)s(t)； (5)
where the subscript p indicates that the physical quantity is transformed into a propagation coordinate system consisting of the f-axis and the s-axis, R (θ) is a rotation matrix, i.e.:
the propagation of the transverse wave along the f-axis and the s-axis is respectively expressed by formula (3), and the displacement vector recorded in the propagation coordinate system is obtained as follows:
d p (t)＝Λ(t)*[R(θ)s(t)] (7)
wherein the operator Λ (t) is:
in the time domain, λ f (t) and lambda s (t) convolution with the source waveform s (t) yields fast and slow shear amplitude attenuation and time shift, respectively, as:
for isotropic media, λ f (t)＝λ s (t)＝λ(t)。
Preferably, said displacement vector d recorded in the propagation coordinate system p (t) transforming back to the coordinate system defined by the detector, resulting in a record of the detector components in the X-axis and Y-axis directions as:
d(t)＝R T (θ)d p (t)＝[R T (θ)Λ(t)R(θ)]*s(t) (10)
expressed as independent variables:
d(t)＝[R T (θ)Λ(t；t f ,t s ,a f ,a s )R(θ)]*s(t)+n(t) (11)
wherein, the superscript T represents the matrix transposition, and n (T) is an added noise item to represent the influence of random noise or local correlation;
a common factor d is extracted from Λ (t) of equation (11) based on the time independence of R (θ) 0 (t)＝a f b(t)δ(t-t f ) And/r, obtaining:
d(t)＝d 0 (t；t f ,a f )*[R T (θ)Λ(t；Δτ；Δa)R(θ)]*s(t)+n(t) (12)
wherein, Δ τ ═ t s -t f A is the time difference between the fast and slow shear waves, Δ a ═ a s /a f Is the relative amplitude between the fast and slow transverse waves.
Preferably, the specific steps of inversely deriving the fast transverse wave direction angle θ and calculating the time difference Δ τ between the fast transverse wave and the slow transverse wave in S4 are as follows:
setting the vibration directions of two orthogonal transverse wave seismic sources as the X-axis direction and the Y-axis direction respectively, recording the transverse wave generated by each seismic source by using a detector with two horizontal components to obtain a four-component data matrix, wherein the four-component data matrix comprises: the XX component and the XY component are respectively records which are excited by a seismic source along the X direction and then received along the X direction and the Y direction; the YX and YY components are respectively records which are excited by a seismic source along the Y direction and then received along the X direction and the Y direction, and the data matrix component form is as follows:
ignoring the noise term, the detector recordings for a set of orthogonal sources are:
without loss of generality, writing them in a matrix form of a dyadic, one can obtain:
D(t)＝[d X (t)|d Y (t)]＝[R T (θ)Λ(t)R(θ)]*[s X (t)|s Y (t)] (15)
if the intensity and waveform of the two seismic sources are the same, then:
thus, d (t) ═ R T (θ)Λ(t)R(θ)I*s(t) (17)
Where I is the identity matrix and s (t) is the general scalar source time function.
Preferably, the R (theta) is an orthogonal matrix, and R (theta) is multiplied on the left side and R (theta) is multiplied on the right side of the formula (17) T (θ), obtaining:
Λ(t)*s(t)＝R(θ)D(t)R T (θ) (18)
FFT is carried out on two sides of the above formula, and the expression in the frequency domain is as follows:
Λ(ω)*s(ω)＝R(θ)D(ω)R T (θ) (19)
the written component is in the form of:
wherein:
and has s f ＝λ f (ω).s(ω)，s s ＝λ s S (ω) with λ in the frequency domain f ＝A f exp(-iωt f )，λ s ＝A s exp(-iωt s )；
The formula (20) is simplified:
wherein
Preferably, the fast transverse wave direction angle θ obtains a θ value with the minimum energy on the off-diagonal component by using a least square method, and the energy function on the off-diagonal component is:
where n is the number of sampling points in the frequency domain, and the first derivative of the extreme point satisfying the function E (θ) is zero, that is:
this gives:
where k is an arbitrary integer, and θ has four solutions spaced by pi/4 in the interval from 0 to pi, corresponding to k being 0,1,2, and 3, respectively, and the minimum point of the function E (θ) satisfies that the second derivative of the function is positive, that is:
the obtained polarization angle is substituted back to the formula (20), and the time difference of the fast and slow transverse waves is calculated as follows:
where real represents the real part of the computation.
Preferably, the example effect analysis of the algorithm in S5 is tested using a set of 2 × 2 four-component VSP synthetic records obtained by computing a multilayer model composed of four anisotropic layers with different forms and parameters using anisieis simulation software.
Preferably, the expression defining the anisotropy factor r in the example effect analysis of the algorithm in S5 is as follows:
wherein, V s1 For fast transverse wave velocity, V s2 Is a slow shear wave velocity.
Compared with the prior art, the invention has the beneficial effects that:
the algorithm is based on the crack model, proves the existence of frequency-dependent anisotropy, and meanwhile, based on the vector convolution model, through proper processing and operation, the multi-component VSP crack parameter inversion algorithm containing the frequency-dependent anisotropy characteristic can be deduced, so that the anisotropy parameters can be directly calculated for different frequencies in a frequency domain conveniently, and errors possibly brought by segmented filtering are eliminated.
Drawings
FIG. 1 is a flow chart of the present invention;
FIG. 2 is a schematic diagram of a coordinate system of a vector convolution model according to the present invention;
FIG. 3 is a diagram of a four component data matrix according to the present invention.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
Referring to fig. 1, the present invention provides a technical solution: a frequency dependent anisotropic shear wave splitting algorithm for multiple VSP data, comprising the steps of:
s1, obtaining a synthetic record of frequency-dependent anisotropy by adopting a numerical simulation method;
s2, the existence of frequency-dependent anisotropy is proved through data analysis and based on a crack model;
s3, deducing a multi-component VSP crack parameter inversion algorithm containing frequency-dependent anisotropic characteristics through processing and operation based on a vector convolution model;
s4, inversion derivation is carried out on a direction angle theta of the fast transverse wave and time difference delta tau of the fast transverse wave and the slow transverse wave is calculated;
and S5, analyzing the example effect of the algorithm.
In this embodiment, preferably, the fracture model in S2 is used to theoretically illustrate that the anisotropy is frequency-dependent, and the fracture model has two different sizes, micro size and macro size: particle size and fracture size, the fluid in the pore space flows according to a jet mechanism during propagation of the elastic wave, and the equivalent elastic stiffness tensor in a fractured porous rock is:
C＝C 0 -Φ p C 1 -ε c C 2 -ε f C 3 (1)
wherein, C 0 Is the elastic constant of isotropic surrounding rock, C 1 、C 2 And C 3 The correction terms are generated when the pores, the microcracks and the cracks act independently; phi p Is the porosity of the surrounding rock, epsilon c Is the density of microcracks,. epsilon. f Is the crack density.
In this embodiment, it is preferable that the particle-sized fluid flow and relaxation time τ m In connection with the inflow and outflow of fluid in the fracture with a lower characteristic frequency or a longer relaxation time tau f Correlation, the correlation formula is as follows:
wherein, a f Is the radius of the crack, a m Is the particle size.
In this embodiment, it is preferable that the formula for the far-field vector wave field in the isotropic medium generated by the point source in the vector convolution model in S3 is as follows:
d(t)＝Aδ(t-τ)*B(t)*s(t)/r (3)
wherein the seismic source is in a homogeneous medium and the vector s (t) represents the direction of polarization; τ represents the travel time (r/V) of the source wave with origin at the source s (t); r represents the total propagation distance, and V represents the wave velocity; b (t) represents a decay factor dependent on the quality factor Q; a represents a constant amplitude factor, and represents convolution in the time domain;
the polarization vector of the seismic source is:
wherein s is X (t) and s Y (t) projections of the seismic source in the X-axis and Y-axis directions, respectively;
vectors s for equivalent motion produced by the seismic source in the directions of the f-axis and the s-axis, orthogonal according to the coordinate system p (t) represents, i.e.:
s p (t)＝R(θ)s(t)； (5)
where the subscript p indicates that the physical quantity is transformed into a propagation coordinate system consisting of f-axis and s-axis, R (θ) is a rotation matrix, i.e.:
the propagation of the transverse wave along the f-axis and the s-axis is respectively expressed by formula (3), and the displacement vector recorded in the propagation coordinate system is obtained as follows:
d p (t)＝Λ(t)*[R(θ)s(t)] (7)
wherein the operator Λ (t) is:
in the time domain, λ f (t) and lambda s (t) convolution of the source waveform s (t) with the source waveform s (t) produces fast shear, slow shear amplitude attenuation, and time shift, respectively, as:
for isotropic media, λ f (t)＝λ s (t)＝λ(t)。
In this embodiment, it is preferable that the displacement vector d recorded in the propagation coordinate system p (t) transforming back to the coordinate system defined by the detector, resulting in a record of the detector components in the X-axis and Y-axis directions as:
d(t)＝R T (θ)d p (t)＝[R T (θ)Λ(t)R(θ)]*s(t) (10)
expressed as independent variables:
d(t)＝[R T (θ)Λ(t；t f ,t s ,a f ,a s )R(θ)]*s(t)+n(t) (11)
wherein, the superscript T represents the matrix transposition, and n (T) is an added noise item to represent the influence of random noise or local correlation;
a common factor d is extracted from Λ (t) of equation (11) based on the time independence of R (θ) 0 (t)＝a f b(t)δ(t-t f ) And/r, obtaining:
d(t)＝d 0 (t；t f ,a f )*[R T (θ)Λ(t；Δτ；Δa)R(θ)]*s(t)+n(t) (12)
wherein, Δ τ ═ t s -t f A is the time difference between the fast and slow shear waves, Δ a ═ a s /a f Is the relative amplitude between the fast and slow transverse waves.
In this embodiment, preferably, the specific steps of inversely deriving the fast transverse wave direction angle θ and calculating the time difference Δ τ between the fast transverse wave and the slow transverse wave in S4 are as follows:
setting the vibration directions of two orthogonal transverse wave seismic sources as the X-axis direction and the Y-axis direction respectively, recording the transverse wave generated by each seismic source by using a detector with two horizontal components to obtain a four-component data matrix, wherein the four-component data matrix comprises: the XX component and the XY component are respectively records which are excited by a seismic source along the X direction and then received along the X direction and the Y direction; the YX and YY components are respectively records which are excited by a seismic source along the Y direction and then received along the X direction and the Y direction, and the data matrix component form is as follows:
ignoring the noise term, the detector recordings for a set of orthogonal sources are:
without loss of generality, writing them in a matrix form of a dyadic, one can obtain:
D(t)＝[d X (t)|d Y (t)]＝[R T (θ)Λ(t)R(θ)]*[s X (t)|s Y (t)] (15)
if the intensity and waveform of the two seismic sources are the same, then:
thus, D (t) R T (θ)Λ(t)R(θ)I*s(t) (17)
Where I is the identity matrix and s (t) is the general scalar source time function.
In this embodiment, preferably, R (θ) is an orthogonal matrix, and R (θ) is multiplied to the left and R (R) is multiplied to the right of formula (17) T (θ), obtaining:
Λ(t)*s(t)＝R(θ)D(t)R T (θ) (18)
FFT is carried out on two sides of the above formula, and the expression in the frequency domain is as follows:
Λ(ω)*s(ω)＝R(θ)D(ω)R T (θ) (19)
the written component is in the form of:
wherein:
and has s f ＝λ f (ω).s(ω)，s s ＝λ s S (ω) with λ in the frequency domain f ＝A f exp(-iωt f )，λ s ＝A s exp(-iωt s )；
The formula (20) is simplified:
wherein
In this embodiment, preferably, the fast transverse wave direction angle θ obtains a θ value with the minimum energy on the off-diagonal component by using a least square method, and the energy function on the off-diagonal component is:
where n is the number of sampling points in the frequency domain, and the first derivative of the extreme point satisfying the function E (θ) is zero, that is:
this gives:
where k is an arbitrary integer, and θ has four solutions spaced by pi/4 in the interval from 0 to pi, which correspond to k being 0,1,2, and 3, respectively, and the minimum point of the function E (θ) satisfies that the second derivative of the function is positive, that is:
the obtained polarization angle is substituted back to the formula (20), and the time difference of the fast and slow transverse waves is calculated as follows:
where real represents the real part of the computation.
In this embodiment, it is preferable that the example effect analysis of the algorithm in S5 is performed by using a set of 2 × 2 four-component VSP synthetic records, where the synthetic records are obtained by calculating a multilayer model with anisieis simulation software, and the multilayer model is composed of four anisotropic layers with different forms and different parameters.
In this embodiment, preferably, the expression of the anisotropy factor r defined in the example effect analysis of the algorithm in S5 is:
wherein, V s1 For fast transverse wave velocity, V s2 Is the slow shear wave velocity.
The principle and the advantages of the invention are as follows: the algorithm is based on the crack model, proves the existence of frequency-dependent anisotropy, and meanwhile, based on the vector convolution model, through proper processing and operation, the multi-component VSP crack parameter inversion algorithm containing the frequency-dependent anisotropy characteristic can be deduced, so that the anisotropy parameters can be directly calculated for different frequencies in a frequency domain conveniently, and errors possibly brought by segmented filtering are eliminated.
Although embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions and alterations can be made in these embodiments without departing from the principles and spirit of the invention, the scope of which is defined in the appended claims and their equivalents.
Claims (10)
1. A frequency-dependent anisotropic shear wave splitting algorithm of multiple VSP data is characterized in that: the method comprises the following steps:
s1, obtaining a synthetic record of frequency-dependent anisotropy by adopting a numerical simulation method;
s2, the existence of frequency-dependent anisotropy is proved through data analysis and based on a crack model;
s3, deducing a multi-component VSP crack parameter inversion algorithm containing frequency-dependent anisotropic characteristics through processing and operation based on a vector convolution model;
s4, inversion derivation is carried out on a direction angle theta of the fast transverse wave and time difference delta tau of the fast transverse wave and the slow transverse wave is calculated;
and S5, analyzing the example effect of the algorithm.
2. The frequency-dependent anisotropic shear wave splitting algorithm of multiple VSP data according to claim 1, wherein: the fracture model in S2 is used to theoretically illustrate that anisotropy is frequency dependent, and has two different sizes, micro and macro: particle size and fracture size, the fluid in the pore space flows according to a jet mechanism during propagation of the elastic wave, and the equivalent elastic stiffness tensor in a fractured porous rock is:
C＝C 0 -Φ p C 1 -ε c C 2 -ε f C 3 (1)
wherein, C 0 Is the elastic constant of isotropic surrounding rock, C 1 、C 2 And C 3 The correction terms are generated when the pores, the microcracks and the cracks act independently; phi p Is the porosity of the surrounding rock, epsilon c Is the density of microcracks,. epsilon. f Is the crack density.
3. The algorithm of frequency dependent anisotropic shear wave splitting of multiple VSP data according to claim 2, characterized by: fluid flow and relaxation time τ of said particle size m In connection with the inflow and outflow of fluid in the fracture with a lower characteristic frequency or a longer relaxation time tau f Correlation, the correlation formula is as follows:
wherein, a f Is the crack radius, a m Is the particle size.
4. The frequency-dependent anisotropic shear wave splitting algorithm of multiple VSP data according to claim 1, wherein: the formula for the far-field vector wavefield in isotropic media generated by a point source in the vector convolution model in S3 is as follows:
d(t)＝Aδ(t-τ)*B(t)*s(t)/r (3)
wherein the seismic source is in a homogeneous medium and the vector s (t) represents the direction of polarization; τ represents the travel time (r/V) of the source wave with origin at the source s (t); r represents the total propagation distance, and V represents the wave velocity; b (t) represents a decay factor dependent on the quality factor Q; a represents a constant amplitude factor, and represents convolution in the time domain;
the polarization vector of the seismic source is:
wherein s is X (t) and s Y (t) projections of the seismic source in the X-axis and Y-axis directions, respectively;
vectors s for equivalent motion produced by the seismic source in the directions of the f-axis and the s-axis, orthogonal according to a coordinate system p (t) represents, i.e.:
s p (t)＝R(θ)s(t)； (5)
where the subscript p indicates that the physical quantity is transformed into a propagation coordinate system consisting of f-axis and s-axis, R (θ) is a rotation matrix, i.e.:
the propagation of the transverse wave along the f-axis and the s-axis is respectively expressed by formula (3), and the displacement vector recorded in the propagation coordinate system is obtained as follows:
d p (t)＝Λ(t)*[R(θ)s(t)] (7)
wherein the operator Λ (t) is:
in the time domain, λ f (t) and lambda s Convolution of (t) and the source waveform s (t) yields the amplitude attenuation and time of the fast and slow transverse waves, respectivelyAnd (3) translating, wherein the expressions are respectively:
for isotropic media, λ f (t)＝λ s (t)＝λ(t)。
5. The frequency-dependent anisotropic shear wave splitting algorithm of multiple VSP data according to claim 4, wherein: said displacement vector d recorded in a propagation coordinate system p (t) transforming back to the coordinate system defined by the detector, resulting in the recording of the detector components in the X-axis and Y-axis directions as:
d(t)＝R T (θ)d p (t)＝[R T (θ)Λ(t)R(θ)]*s(t) (10)
expressed as independent variables:
d(t)＝[R T (θ)Λ(t；t f ,t s ,a f ,a s )R(θ)]*s(t)+n(t) (11)
wherein, the superscript T represents the matrix transposition, and n (T) is an added noise item to represent the influence of random noise or local correlation;
a common factor d is extracted from Λ (t) of equation (11) based on the time independence of R (θ) 0 (t)＝a f b(t)δ(t-t f ) And/r, obtaining:
d(t)＝d 0 (t；t f ,a f )*[R T (θ)Λ(t；Δτ；Δa)R(θ)]*s(t)+n(t) (12)
wherein, Δ τ ═ t s -t f A is the time difference between the fast and slow shear waves, Δ a ═ a s /a f Is the relative amplitude between the fast and slow transverse waves.
6. The frequency-dependent anisotropic shear wave splitting algorithm of multiple VSP data according to claim 1, wherein: the specific steps of inversion derivation of the fast transverse wave direction angle theta and calculation of the time difference delta tau of the fast transverse wave and the slow transverse wave in the step S4 are as follows:
setting the vibration directions of two orthogonal transverse wave seismic sources as the X-axis direction and the Y-axis direction respectively, recording the transverse wave generated by each seismic source by using a detector with two horizontal components to obtain a four-component data matrix, wherein the four-component data matrix comprises: the XX component and the XY component are respectively records which are excited by a seismic source along the X direction and then received along the X direction and the Y direction; the YX and YY components are respectively records which are excited by a seismic source along the Y direction and then received along the X direction and the Y direction, and the data matrix component form is as follows:
ignoring the noise term, the detector recordings for a set of orthogonal sources are:
without loss of generality, writing them in the form of a matrix of dyads, one can obtain:
D(t)＝[d X (t)|d Y (t)]＝[R T (θ)Λ(t)R(θ)]*[s X (t)|s Y (t)] (15)
if the intensity and waveform of the two seismic sources are the same, the following are obtained:
thus, d (t) ═ R T (θ)Λ(t)R(θ)I*s(t) (17)
Where I is the identity matrix and s (t) is a general scalar source time function.
7. The algorithm of claim 6 for frequency dependent anisotropic shear wave splitting of multiple VSP data, characterized by: the R (theta) is an orthogonal matrix, and is multiplied by the R (theta) on the left side and the R (theta) on the right side of the formula (17) T (θ), obtaining:
Λ(t)*s(t)＝R(θ)D(t)R T (θ) (18)
FFT is carried out on two sides of the above formula, and the expression in the frequency domain is as follows:
Λ(ω)*s(ω)＝R(θ)D(ω)R T (θ) (19)
the written component is in the form of:
wherein:
and has s f ＝λ f (ω).s(ω)，s s ＝λ s S (ω) with λ in the frequency domain f ＝A f exp(-iωt f )，λ s ＝A s exp(-iωt s )；
The formula (20) is simplified:
wherein
8. The frequency-dependent anisotropic shear wave splitting algorithm of multiple VSP data according to claim 7, wherein: and obtaining a theta value with minimum energy on the off-diagonal component by using a least square method for the fast transverse wave direction angle theta, wherein an energy function on the off-diagonal component is as follows:
where n is the number of sampling points in the frequency domain, and the extreme point of n satisfies that the first derivative of the function E (θ) is zero, that is:
this gives:
where k is an arbitrary integer, and θ has four solutions spaced by pi/4 in the interval from 0 to pi, corresponding to k being 0,1,2, and 3, respectively, and the minimum point of the function E (θ) satisfies that the second derivative of the function is positive, that is:
the obtained polarization angle is substituted back to the formula (20), and the time difference of the fast and slow transverse waves is calculated as follows:
where real represents the real part of the computation result.
9. The algorithm of claim 1 for frequency dependent anisotropic shear wave splitting of multiple VSP data, characterized by: in the S5, a set of 2 × 2 four-component VSP synthetic records obtained by calculating a multilayer model composed of four anisotropic layers with different forms and parameters by using ANISEIS simulation software is used for the example effect analysis of the algorithm.
10. The frequency-dependent anisotropic shear wave splitting algorithm of multiple VSP data according to claim 9, wherein: the expression defining the anisotropy factor r in the example effect analysis of the algorithm in S5 is:
wherein, V s1 For fast transverse wave velocity, V s2 Is a slow shear wave velocity.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210650897.7A CN115097520A (en) | 2022-06-09 | 2022-06-09 | Frequency-dependent anisotropic shear wave splitting algorithm for multiple VSP data |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210650897.7A CN115097520A (en) | 2022-06-09 | 2022-06-09 | Frequency-dependent anisotropic shear wave splitting algorithm for multiple VSP data |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115097520A true CN115097520A (en) | 2022-09-23 |
Family
ID=83290740
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210650897.7A Pending CN115097520A (en) | 2022-06-09 | 2022-06-09 | Frequency-dependent anisotropic shear wave splitting algorithm for multiple VSP data |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115097520A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117572495A (en) * | 2023-10-24 | 2024-02-20 | 成都理工大学 | Quantitative prediction method for crack scale |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20030195705A1 (en) * | 2002-04-10 | 2003-10-16 | Scott Leaney | Method and apparatus for anisotropic vector plane wave decomposition for 3D vertical seismic profile data |
| WO2015014762A2 (en) * | 2013-07-29 | 2015-02-05 | Cgg Services Sa | Method and device for the generation and application of anisotropic elastic parameters in horizontal transverse isotropic (hti) media |
| CN110687601A (en) * | 2019-11-01 | 2020-01-14 | 中南大学 | Inversion method for fluid factor and fracture parameter of orthotropic medium |
-
2022
- 2022-06-09 CN CN202210650897.7A patent/CN115097520A/en active Pending
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20030195705A1 (en) * | 2002-04-10 | 2003-10-16 | Scott Leaney | Method and apparatus for anisotropic vector plane wave decomposition for 3D vertical seismic profile data |
| WO2015014762A2 (en) * | 2013-07-29 | 2015-02-05 | Cgg Services Sa | Method and device for the generation and application of anisotropic elastic parameters in horizontal transverse isotropic (hti) media |
| CN110687601A (en) * | 2019-11-01 | 2020-01-14 | 中南大学 | Inversion method for fluid factor and fracture parameter of orthotropic medium |
Non-Patent Citations (1)
| Title |
|---|
| 韩开锋: "含分布裂缝岩石中弹性波传播特性研究", 《中国博士学位论文全文数据库-基础科学辑》, no. 5, pages 011 - 26 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117572495A (en) * | 2023-10-24 | 2024-02-20 | 成都理工大学 | Quantitative prediction method for crack scale |
| CN117572495B (en) * | 2023-10-24 | 2024-05-17 | 成都理工大学 | Quantitative prediction method for crack scale |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Carcione et al. | Computational poroelasticity—A review | |
| US20090132169A1 (en) | Methods and systems for evaluating fluid movement related reservoir properties via correlation of low-frequency part of seismic data with borehole measurements | |
| Sun et al. | A sweeping window method for detection of flaws using an explicit dynamic XFEM and absorbing boundary layers | |
| Zerwer et al. | Rayleigh wave propagation for the detection of near surface discontinuities: Finite element modeling | |
| US20130289881A1 (en) | Hydraulic Fracture Characterization Using Borehole Sonic Data | |
| Zheng et al. | Simulation of the borehole quasistatic electric field excited by the acoustic wave during logging while drilling due to electrokinetic effect | |
| Alkhimenkov et al. | Azimuth-, angle-and frequency-dependent seismic velocities of cracked rocks due to squirt flow | |
| CN115097520A (en) | Frequency-dependent anisotropic shear wave splitting algorithm for multiple VSP data | |
| Yang et al. | A strong stability-preserving predictor-corrector method for the simulation of elastic wave propagation in anisotropic media | |
| Barbosa et al. | Impact of fracture clustering on the seismic signatures of porous rocks containing aligned fractures | |
| Yang et al. | Poro-acoustoelasticity finite-difference simulation of elastic wave propagation in prestressed porous media | |
| Lambert et al. | Attenuation and dispersion of P-waves in porous rocks with planar fractures: Comparison of theory and numerical simulations | |
| Lissa et al. | Fluid pressure diffusion in fractured media: The role played by the geometry of real fractures | |
| Nakagawa et al. | Three-dimensional elastic wave scattering by a layer containing vertical periodic fractures | |
| Song et al. | Attenuation and dispersion of P-waves in fluid-saturated porous rocks with a distribution of coplanar cracks—Scattering approach | |
| Tadeu et al. | 2.5 D elastic wave propagation in non-homogeneous media coupling the BEM and MLPG methods | |
| Xu et al. | Effects of single vertical fluid-filled fractures on full waveform dipole sonic logs | |
| Wang et al. | Effects of the outermost boundary on acoustic waves in an artificial cased borehole | |
| Sayers | Increasing contribution of grain boundary compliance to polycrystalline ice elasticity as temperature increases | |
| Rahnema et al. | 3D numerical modeling of multi-channel analysis of surface wave in homogeneous and layered concrete slabs | |
| Bhaumik et al. | Higher-order thin layer method as an efficient forward model for calculating dispersion curves of surface and Lamb waves in layered media | |
| Liu et al. | Frequency-domain FD modeling with an adaptable nearly perfectly matched layer boundary condition for poroviscoelastic waves upscaled from the effective Biot theory | |
| Lu | Analysis and modeling of diffuse ultrasonic signals for structural health monitoring | |
| Verweij et al. | Transmission and reflection of transient elastodynamic waves at a linear slip interface | |
| Bécache et al. | A fictitious domain method with mixed finite elements for elastodynamics |
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 | ||
| CB02 | Change of applicant information |
Address after: 410073 Hunan province Changsha Kaifu District, Deya Road No. 109 Applicant after: National University of Defense Technology Address before: No. 60, Shuanglong Avenue, Jiangning District, Nanjing, Jiangsu 210000 Applicant before: National University of Defense Technology |
|
| CB02 | Change of applicant information |