CN107894616A - Multi-component converted wave crack prediction method - Google Patents
Multi-component converted wave crack prediction method Download PDFInfo
- Publication number
- CN107894616A CN107894616A CN201711122137.4A CN201711122137A CN107894616A CN 107894616 A CN107894616 A CN 107894616A CN 201711122137 A CN201711122137 A CN 201711122137A CN 107894616 A CN107894616 A CN 107894616A
- Authority
- CN
- China
- Prior art keywords
- component
- wave
- mrow
- mtd
- msub
- 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 40
- 238000011161 development Methods 0.000 claims abstract description 35
- 238000004458 analytical method Methods 0.000 claims abstract description 16
- 238000006243 chemical reaction Methods 0.000 claims abstract description 6
- 238000005070 sampling Methods 0.000 claims abstract description 5
- 206010017076 Fracture Diseases 0.000 claims description 42
- 208000010392 Bone Fractures Diseases 0.000 claims description 41
- 239000011159 matrix material Substances 0.000 claims description 11
- 238000009826 distribution Methods 0.000 claims description 7
- 238000000926 separation method Methods 0.000 claims description 7
- 230000015572 biosynthetic process Effects 0.000 claims description 4
- 230000004927 fusion Effects 0.000 claims description 4
- 238000009987 spinning Methods 0.000 claims description 4
- 238000003786 synthesis reaction Methods 0.000 claims description 4
- 238000011426 transformation method Methods 0.000 claims description 4
- 238000004364 calculation method Methods 0.000 abstract description 4
- 238000004422 calculation algorithm Methods 0.000 abstract description 2
- 238000012545 processing Methods 0.000 description 6
- 238000007796 conventional method Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 239000003208 petroleum Substances 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 239000013598 vector Substances 0.000 description 3
- 238000001514 detection method Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000010287 polarization Effects 0.000 description 2
- 230000000644 propagated effect Effects 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 125000004122 cyclic group Chemical group 0.000 description 1
- 230000001351 cycling effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
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/301—Analysis for determining seismic cross-sections or geostructures
-
- 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
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/64—Geostructures, e.g. in 3D data cubes
- G01V2210/646—Fractures
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 provides a multi-component converted wave crack prediction method, which comprises the following steps: selecting an analysis window containing a target depth range, and intercepting a radial component and a transverse component of a sampling point to be processed according to the analysis window; carrying out angle rotation on the intercepted radial component and transverse component, and carrying out time delay conversion on the transverse component after the angle rotation; constructing an objective function and calculating the minimum value of the objective function; recording a rotation angle and time delay corresponding to the minimum value of the target function; separating fast transverse waves and slow transverse waves of the sampling points to be processed, and calculating the development strength of the crack; and repeating the steps, and circulating to the next sampling point to be processed until the seismic data is calculated. The invention introduces a non-parameter statistical theory into the multi-wave fracture prediction process, adopts a grade correlation mode to replace a linear correlation method in the traditional algorithm, improves the calculation accuracy and provides a basis for the fine description of the fractured oil and gas reservoir.
Description
Technical field
The invention belongs to petroleum exploration field, is related to a kind of Fracture System detection method for Fractured oil and gas reservoir,
More particularly, it is related to a kind of multi-component converted wave crack prediction method.
Background technology
From from the point of view of development of fields statistics, fractured reservoir is in petroleum resources and production capacity side in the world
Face, account for the half of world's total amount.Therefore, Fractured oil and gas reservoir turns into the emphasis of current geophysical prospecting for oil research.
How accurately to judge that existence of the crack in underground has for reservoir evaluation, the assessment of block and the raising rate of oil and gas recovery
Vital effect.The currently conventional crack prediction method mostly method based on poststack seismic properties, such as coherence properties,
Curvature attributes, texture properties etc., but the method based on poststack seismic properties is applied to prediction large scale fracture belt more, can not be abundant
It is not fine enough using the anisotropy information in prestack information and stratum, prediction result.Therefore, it is necessary to explore more fine split
Stitch Predicting Technique.
With deepening continuously for more ripples research, multiband fusion technology is applied in oil-gas exploration and development more and more.It is horizontal
Wavefront splitting has the technology of powerful potentiality by geologist and geophysics as earthquake, one of petroleum exploration field
Family's extensive concern.Shear wave splitting, also known as shear-wave birefringence, i.e., it can split into edge when shear wave passes through direction anisotropy medium
The fast transverse wave propagated parallel to fissure direction and the slow shear-wave propagated along vertical fissure direction, both spread speeds are different, therefore
The Shear Waves Splitting time difference, amplitude decay etc. are formed in anisotropic medium, these characteristics are all relevant with anisotropic properties.Can root
Study the orientation and density of Fracture System in turn according to fast transverse wave, the time difference of slow shear-wave, waveform, amplitude decay etc..Utilize shear wave
Dividing Predicating Reservoir Fractures turns into a kind of directly reliable method.But shear wave source is costly, using pure shear wave source come
Carry out shear wave splitting exploration and be also difficult to large-scale application.Converted P-SV-Waves Exploration overcome pure shear wave exploration excite difficult, cost is high,
The defects of static correction value is big, while there is abundant information, have compressional wave concurrently and the advantage of shear wave, data signal to noise ratio are higher,
So that multi-component converted wave exploration turns into the powerful of oil and gas reservoir detection.
In Converted P-SV-Waves Exploration FRACTURE PREDICTION, fracture development orientation and the time difference are calculated exactly for accurately portraying crack
Distributional pattern of the system in underground is significant.At present, conventional method has cross-correlation method and energy ratio function.Conventional method
When the correlation or covariance matrix of middle calculating, by the way of linearly related, and the similitude between waveform is non-linear
, and the data sample finite capacity that parameter calculates, the feature of normal distribution is unsatisfactory for, from mathematical meaning, sample data
Inherently it is unsatisfactory for requiring.
The content of the invention
For the deficiencies in the prior art, an object of the present invention is solve present in above-mentioned prior art
One or more problems.
To achieve these goals, it is of the invention to provide a kind of multi-component converted wave crack prediction method, the prediction
Method may comprise steps of：An analysis window for including target depth range is chosen, by the radial direction of pending sampled point
Component earthquake data and cross stream component geological data are intercepted according to analysis window；To the radial component geological data of interception and
Cross stream component geological data carries out angle rotation, and carries out time delay conversion to the postrotational cross stream component geological data of angle；Structure
Build object function, calculating target function minimum value；Record the anglec of rotation and time delay corresponding to object function minimum value；To pending
Sampled point carries out fast transverse wave and slow shear-wave separation, and fracture development intensity is calculated according to the travel-time difference of fast transverse wave and slow shear-wave；Weight
The step of the step of multiple Analysis on Selecting window to calculating fracture development intensity, it is recycled to next pending sampled point
Finished until geological data calculates, wherein, the structure object function, include the step of calculating target function minimum value：
It is P (θ that the radial component geological data of note interception, which carries out the postrotational component of angle,i+n, t), the postrotational horizontal stroke of angle
The component changed to component earthquake data progress time delay is Q (θi+n,t+ti+n)；Respectively to P (θi+n,t)、Q(θi+n,t+ti+n) enter
Row compiles order operation, is designated as rgP and rgQ, wherein, θi+nThe anglec of rotation is represented, t represents time, ti+nTime delay is represented, n represents integer；
Build matrixWherein,M represents P or Q,
N represents P or Q, ρrgM,rgNPearson correlation coefficients are represented, cov (rgM, rgN) represents to ask covariance function, σrgPRepresent rgP's
Standard deviation, σrgQRepresent rgQ standard deviation；Calculating matrix U (θi+n,t+ti+n) characteristic value eig (U (θi+n,t+ti+n))；Wherein,
θi+n=θi+ n Δs θ, ti+n=ti+ n Δs t, θmin≤θi+n≤θmaxAnd tmin≤ti+n≤tmax, θiTo originate the anglec of rotation, tiTo rise
Beginning time delay, Δ θ are that the anglec of rotation circulates step-length, and Δ t is that time delay circulates step-length, [θmin,θmax] it is anglec of rotation threshold range,
[tmin,tmax] it is delay threshold scope；Build object function F (θ, t)=min [eig (U (θi+n,t+ti+n))], n takes since 0
Value, n increase by 1, obtain object function minimum value every time, wherein, min [eig (U (θi+n,t+ti+n))] represent to seek characteristic value eig (U
(θi+n,t+ti+n)) minimum value.
Compared with prior art, during Non-parametric approach is introduced more ripple FRACTURE PREDICTIONs by the present invention, using grade
Related mode substitutes the linear correlative method in traditional algorithm, and rank correlation mode is quicker for estimating for waveform similarity
Sense, while evaded the linearly related requirement to sample data, the degree of accuracy of calculating is improved, is the essence of Fractured oil and gas reservoir
Thin description provides foundation.
Brief description of the drawings
By the description carried out below in conjunction with the accompanying drawings, above and other purpose of the invention and feature will become more clear
Chu, wherein：
Fig. 1 shows multi-component converted wave crack prediction method schematic flow sheet according to an exemplary embodiment of the present invention.
Fig. 2 shows shear wave splitting according to an exemplary embodiment of the present invention and fracture development position relation schematic diagram.
Embodiment
Hereinafter, it will describe in detail with reference to accompanying drawing and exemplary embodiment and changed according to a kind of multi -components of the present invention
Ripple crack prediction method.
Specifically, multi-component converted wave can produce separating phenomenon in EDA media, can be carried out using converted shear wave
Fractured Zone predicts that the fine description for Fractured oil and gas reservoir provides foundation.Converted shear wave refers to use during ground observation
Compressional wave excites, and inclined seismic wave is converted shear wave to that can produce reflection wave, up shear wave therein during elastic interface.
Converted shear wave possesses whole properties of shear wave, and when up converted shear wave and fracture strike oblique, separating phenomenon can equally occur,
It is split into fast transverse wave and slow shear-wave.Fast, slow wave after division travel to earth's surface, decompose and lay equal stress on according to earth's surface observation system coordinate system
New synthesis, is received by multi -components wave detector.Therefore, when survey line and fractuer direction are not parallel, in ground level component recording all
Containing fast, slow wave information, then carry out fast, slow wave using multi -components record and separate, according to the time difference of fast wave and slow wave, waveform,
The information such as orientation, density and the time difference of Fracture System are studied in amplitude decay etc. in turn, and predict crack hair according to information above
Educate.
The innovation of the present invention is that the method for nonparametric statistics is introduced crack prediction method by method of the invention
In, the data after conversion are carried out when building object function to compile order operation, the covariance matrix based on order are then built, by line
Property correlation be changed into nonlinear correlation, improve the degree of accuracy of calculating, for subterranean fracture identify and predict more reliable foundation is provided.
Fig. 1 shows multi-component converted wave crack prediction method schematic flow sheet according to an exemplary embodiment of the present invention.
Fig. 2 shows shear wave splitting according to an exemplary embodiment of the present invention and fracture development position relation schematic diagram.
The invention provides a kind of multi-component converted wave crack prediction method, as shown in figure 1, in the example of the present invention
In property embodiment, methods described can include：
Step S01, the horizontal component X and Y of converted wave Three-dimendimal fusion geological data are rotated, obtained postrotational
Converted wave radial component geological data R (i.e. radial component R) and cross stream component geological data T (i.e. cross stream component T).
More than, three-component receiver can be used in D3C seismic exploration, three-component receiver is respectively water
Flat X, Y-component receiver, and vertical Z component receiver.Due to Inline and CrossLine directions and shear wave polarization side
To inconsistent, cause X, the existing P-SV ripples of shear wave recorded in Y-direction also there are P-SH ripples, do not polarize implication clearly, it is each
Wave field on polarization direction mutually mixes, and is unfavorable for the processing of converted wave data and the research and extraction of shear wave splitting information.Cause
This, in order to obtain consistent conversion wave field, when handling Three-dimendimal fusion data, it is necessary to two component energies by horizontal direction
Redistribute, converted wave coordinate is rotated, realize from field acquisition coordinate system (X-Y coordinate) to radial direction-lateral coordinates system (R-
T coordinate systems) rotation.After coordinate rotation processing, converted wave main energetic is distributed on the radial component, cross stream component main generation
The anisotropic influence of table.Coordinate can be rotatably：
Wherein, radial component R and the angle in X-component directionOr
Anisotropy processing operation can be carried out to postrotational converted wave, the anisotropy processing operation belongs to ability
The conventional method in domain.
Step S02, analysis window scope is determined, the radial component earthquake number of pending sampled point is intercepted according to analysis window
According to R and cross stream component geological data T.
In the present example embodiment, the selection of analysis window need to include target depth range, i.e., comprising Fractured respectively to
Different in nature stratum, can be controlled or period control targe depth bounds by layer position.Since pending sampled point, work as analysis
After window determines, the radial component geological data and cross stream component geological data to all sampled points are respectively according to the analysis of determination
Window ranges are intercepted, and the data of interception are stored in different arrays to carry out follow-up calculating respectively.
More than, in same work area, the layer position scope that the analysis window scope of all sampled points is chosen is consistent, example
Such as, controlled by layer position, the window scope of all sampled points is taken comprising layer position 1 and layer position 2, but model corresponding to each sampled point
It is different in layer position 1 and layer position 2 to enclose value, for example sampled point C scope can take 1200ms (layer position 1)~1300ms
(layer position 2), sampled point D scope can take 1300ms (layer position 1)~1400ms (layer position 2).
Step S03, the anglec of rotation is carried out to radial direction component earthquake data and cross stream component geological data and time delay is changed, is obtained
Component after to conversion.
In the present example embodiment, because the data recorded in radial component and cross stream component are fast transverse wave and slow shear-wave
Information reconfiguring when earth's surface receives, i.e., be in radial component and cross stream component fast transverse wave and slow shear-wave information it is folded
Add.Therefore, it is necessary to be separated fast transverse wave and slow shear-wave by way of rotation.Fast transverse wave and slow shear-wave are by same simultaneously
One splitting of converted shear and obtain, their waveform is similar, and the spread speed of slow shear-wave is slower than fast transverse wave spread speed, two
There is the time difference between person.Therefore, it is necessary to carry out delay compensation to slow shear-wave on the basis of angle rotation.
Specifically, the radial component R and cross stream component T of the interception to pending sample point respectively first carry out angle rotation
Turn, can be rotatably：
Obtain component P (θi+n, t) and Q (θi+n, t), time delay, Q (θ then are carried out to postrotational cross stream component datai+n,t
+ti+n) it is Q (θi+n, t) in the upward time shift t of each valuei+n/ S point, S are sample rate, θi+nFor the anglec of rotation, ti+nFor time delay,
θiTo originate the anglec of rotation, tiFor initial delay, n represents integer, θi+n=θi+ n Δs θ, ti+n=ti+ n Δs t, θmin≤θi+n≤
θmaxAnd tmin≤ti+n≤tmax, Δ θ is that the anglec of rotation circulates step-length, and Δ t is that time delay circulates step-length, [θmin,θmax] it is the anglec of rotation
Spend threshold range, [tmin,tmax] it is delay threshold scope.
More than, Q (θi+n,t+ti+n) it is Q (θi+n, t) in the upward time shift t of each valuei+n/ S point refer to for example, when
Move ti+n/ S takes 1 point, then Q (θi+n, t) in data group Q1,Q2,Q3,Q4... one point of time shift, Q (θi+n,t+ti+n) corresponding to
Array is changed into from Q2Start, as Q2,Q3,Q4...。
Since start angle and initial delay, i.e., initial value starts, and by angle rotary and time-delay calculation formula, obtains
The postrotational component of radial component geological data angle of interception is P (θi+n, t), and to the postrotational cross stream component earthquake of angle
Data carry out the component after time delay is changed as Q (θi+n,t+ti+n)。
Here, time-delay calculation is mainly for postrotational cross stream component, i.e., to the postrotational cross stream component earthquake number of angle
According to P (θi+n, t) and the component after time delay is changed is carried out as Q (θi+n,t+ti+n).Anglec of rotation θi+nWith time delay ti+nIt is required in threshold value
In the range of value, threshold range can be setting value or empirical value.For example, the threshold range of the anglec of rotation typically can be [- π,
π] or [0,2 π].
For example, set processing after multi-component data sample rate as S, analysis window is controlled by layer position, layer position hor1For
When window top, hor2For when window bottom.Calculating is started the cycle over from sampled point A (p, q), wherein, p represents sampled point A wire size, and q is represented
Sampled point A Taoist monastic name, by radial component data and cross stream component data, window ranges are intercepted according to the above analysis, interception
Data are stored in array R=[r respectively1,r2,r3,......,rn] and array T=[t1,t2,t3,......,tn] in calculated.
Then, θ is carried out to array R and Ti+nThe rotation of angle, for example, the threshold range of the anglec of rotation can be set as 0 ＜
θi+n＜ π；Then delay compensation is carried out to postrotational cross stream component data, the threshold range of time delay can be set to [0, T], obtain
The postrotational component of radial component geological data angle to interception is P (θi+n, t), the postrotational cross stream component earthquake number of angle
Component after being changed according to progress time delay is Q (θi+n,t+ti+n)。
Step S04, object function is built, asks for object function minimum value.
In the present example embodiment, the structure of object function can include：
Step S04.1, first, to component P (θi+n, t) and Q (θi+n,t+ti+n) carry out compiling order operation, i.e., ask for respectively point
Measure P (θi+n, t) and Q (θi+n,t+ti+n) order, be designated as rgP and rgQ, wherein, θi+nThe anglec of rotation is represented, t represents time, ti+n
Time delay is represented, n represents integer.
Step S04.2, build matrixWherein,
M represents that P or Q, N represent P or Q, i.e.,
ρrgM,rgNPearson correlation coefficients are represented, cov (rgM, rgN) represents to ask covariance function, σrgPRepresent rgP mark
Poor, the σ of standardrgQRepresent rgQ standard deviation.
Step S04.3, calculating matrix U (θi+n,t+ti+n) characteristic value eig (U (θi+n,t+ti+n)).Wherein, θi+n=θi+
N Δs θ, ti+n=ti+ n Δs t, θmin≤θi+n≤θmaxAnd tmin≤ti+n≤tmax, θiTo originate the anglec of rotation, tiFor initial delay,
Δ θ is that the anglec of rotation circulates step-length, and Δ t is that time delay circulates step-length, [θmin,θmax] it is anglec of rotation threshold range, [tmin,tmax]
For delay threshold scope.
Here, it is necessary to judge anglec of rotation θi+nWith time delay ti+nWhether in threshold range, if meeting threshold condition,
It can continue cycling through and take the new anglec of rotation and time delay value, bring step S03 into, obtain different component P (θi+n, t) and Q (θi+n,
t+ti+n), and then obtain different characteristic values.
Step S04.4, structure object function F (θ, t)=min [eig (U (θi+n,t+ti+n))], n values since 0, n is every
Secondary increase by 1, obtains object function minimum value, wherein, min [eig (U (θi+n,t+ti+n))] represent to seek characteristic value eig (U (θi+n,t
+ti+n)) minimum value.
Step S05, the corresponding anglec of rotation and time shift when record object function F (θ, t) takes minimum value.
Here, when object function takes minimum value, the anglec of rotation now is the azimuth of fracture development, and time delay is
The time difference of Concerning With Fast-slow Waves.When converting wavefront splitting, the Concerning With Fast-slow Waves time difference is bigger, shows that crack is more developed, and therefore, time differences of Concerning With Fast-slow Waves can be with
A kind of important indicator as development degree of micro cracks in oil.
Step S06, fast transverse wave is carried out to pending sampled point and slow shear-wave separates, according to the travelling of fast transverse wave and slow shear-wave
The time difference calculates fracture development intensity.
In the present example embodiment, radial component geological data and cross stream component geological data include fast, slow horizontal stroke
Ripple information, therefore, it is necessary to parallel-vertical shear wave is separated from radial component geological data and cross stream component geological data.
When confirming to have fracture development to be sampled, triangular transformation method or Alford spinning solutions can be used.Three
Angular transformation method and Alford spinning solutions be cross stream component geological data by way of coordinate rotation at fracture with
And radial component geological data carries out isolated fast transverse wave and slow shear-wave.
More than, triangular transformation method is using formulaFast transverse wave and slow shear-wave are obtained,
Wherein, S1Represent fast transverse wave, S2Represent slow shear-wave.
The basic ideas of Alford spinning solutions can include：If fracture strike azimuthal angle beta and radial component forward direction orientation
Angle α angle is θ=β-α,
Defining orthogonal spin matrix G is：
Construction receives signal matrix：
Wherein, R1、T1And R2、T2The radius vector and transversal vector of respectively two different orientations are full between two azimuths
Sufficient orthogonality relation.
According to Alford spin theories：V=RSRT,
Wherein,S1Represent fast transverse wave, S2Represent slow shear-wave,
To Formula V=RSRTAfter carrying out orthogonal rotation, S=R is obtainedTVR。
The method for calculating fracture development intensity includes：Fracture development intensity KC=Δ T/TS1, wherein, Δ T=TS1-
TS2, TS1For the fast transverse wave section back wave existence time after separation, TS2In the presence of the slow shear-wave section back wave after separation
Between.
Step S07, repeat step S01 are recycled to next pending sampled point until geological data calculates to step S06
Finish.
Here, under treatment one to be sampled when, it is also necessary to upper one processing sampled point wire size and Taoist monastic name enter
Row judges, for example, it is assumed that sampled point A (p, q), when meeting start_inline ＜ p ＜ end_inline, start_
During crossline ＜ q ＜ end_crossline, then next pending sampled point is may loop to, until geological data calculates
Finish, wherein, p is the wire size of pending sampled point, and q is the Taoist monastic name of pending sampled point, and start_inline is work area starting
Wire size, end_inline are that work area terminates wire size, and start_crossline is work area starting Taoist monastic name, and end_crossline is work
Area terminates Taoist monastic name.
It should be noted that in circular treatment sampled point, sampling the order of dot cycle can enter according to wire size and Taoist monastic name
Row circulation.For example, in cyclic process, first keep certain wire size p constant, in Taoist monastic name scope start_crossline ＜ q ＜ end_
Circulated successively according to Taoist monastic name in crossline, after the completion of Taoist monastic name circulation, wire size is changed into p+1, then circulates all Taoist monastic names successively.
Step S08, draw fracture development synthesis distribution map.
Here, after all sampled points calculate, the data result of each sampled point is preserved, forms two numbers
According to i.e. fracture development azimuth and fracture development intensity, then drawing and split in comprehensive fracture development orientation and fracture development intensity
The comprehensive distribution map of seam development.The method for drawing fracture development synthesis distribution map includes：The fracture development of each sampled point is used
Short-term represents that the azimuth of the orientation references fracture development of short-term, the length of short-term represents fracture development intensity.
Shear wave splitting and the graph of a relation in fracture development orientation are as shown in Figure 2.Fast transverse wave is along fracture development side after shear wave splitting
Position, slow shear-wave and fracture development oriented perpendicular, according to the principle of stacking of ripple, include in radial component and cross stream component data
Speed wave component.
In summary, the present invention is evaded by the way that nonparametric statistical method is incorporated into multi-component converted wave FRACTURE PREDICTION
Requirement of the conventional linear correlation technique to earthquake data distribution, when carrying out covariance matrix calculating, instead of rank correlation
Result of calculation, the correlation between two vectors can be described well, improve accuracy in computation, make prediction result it is more accurate,
Reliably.
Although above by describing the present invention with reference to exemplary embodiment, those skilled in the art should be clear
Chu, in the case where not departing from the spirit and scope that claim is limited, the exemplary embodiment of the present invention can be carried out each
Kind modifications and changes.
Claims (8)
1. a kind of multi-component converted wave crack prediction method, it is characterised in that the Forecasting Methodology comprises the following steps：
An analysis window for including target depth range is chosen, by the radial component geological data and transverse direction of pending sampled point
Component earthquake data is intercepted according to analysis window；
Radial component geological data and cross stream component geological data to interception carry out angle rotation, and to the postrotational horizontal stroke of angle
Time delay conversion is carried out to component earthquake data；
Build object function, calculating target function minimum value；
Record the anglec of rotation and time delay corresponding to object function minimum value；
Fast transverse wave and slow shear-wave separation are carried out to pending sampled point, crack is calculated according to the travel-time difference of fast transverse wave and slow shear-wave
Growth strength；
The step of the step of repeating the Analysis on Selecting window to calculating fracture development intensity, it is recycled to next pending
Sampled point finishes until geological data calculates, wherein,
The structure object function, include the step of calculating target function minimum value：
It is P (θ that the radial component geological data of note interception, which carries out the postrotational component of angle,i+n, t), angle postrotational transverse direction point
The component that amount geological data progress time delay is changed is Q (θi+n,t+ti+n)；
Respectively to P (θi+n,t)、Q(θi+n,t+ti+n) carry out compiling order operation, rgP and rgQ are designated as, wherein, θi+nRepresent the anglec of rotation
Degree, t represent time, ti+nTime delay is represented, n represents integer；
Build matrixWherein,
M represents that P or Q, N represent P or Q, ρrgM,rgNRepresent Pearson correlation coefficients,
Cov (rgM, rgN) represents to ask covariance function, σrgPRepresent rgP standard deviation, σrgQRepresent rgQ standard deviation；
Calculating matrix U (θi+n,t+ti+n) characteristic value eig (U (θi+n,t+ti+n))；
Wherein, θi+n=θi+ n Δs θ, ti+n=ti+ n Δs t, θmin≤θi+n≤θmaxAnd tmin≤ti+n≤tmax, θiTo originate the anglec of rotation
Degree, tiFor initial delay, Δ θ is that the anglec of rotation circulates step-length, and Δ t is that time delay circulates step-length, [θmin,θmax] it is anglec of rotation threshold
It is worth scope, [tmin,tmax] it is delay threshold scope；
Build object function F (θ, t)=min [eig (U (θi+n,t+ti+n))], n values since 0, n increases by 1 every time, obtains mesh
Scalar functions minimum value, wherein, min [eig (U (θi+n,t+ti+n))] represent to seek characteristic value eig (U (θi+n,t+ti+n)) minimum
Value.
2. multi-component converted wave crack prediction method according to claim 1, it is characterised in that the radial component earthquake
Data and cross stream component geological data are that the horizontal component X and Y of Three-dimendimal fusion geological data are rotated and obtained, described
The coordinate that horizontal component X and Y are rotated is rotatably：
Wherein,For radial component and horizontal component X angle,OrR represents radially to divide
Geological data is measured, T represents cross stream component geological data.
3. multi-component converted wave crack prediction method according to claim 1, it is characterised in that the component P (θi+n,t)
With Q (θi+n,t+ti+n) calculating formula be：
<mrow>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mi>P</mi>
<mo>(</mo>
<msub>
<mi>&theta;</mi>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mi>n</mi>
</mrow>
</msub>
<mo>,</mo>
<mi>t</mi>
<mo>)</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>Q</mi>
<mo>(</mo>
<msub>
<mi>&theta;</mi>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mi>n</mi>
</mrow>
</msub>
<mo>,</mo>
<mi>t</mi>
<mo>)</mo>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>cos&theta;</mi>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mi>n</mi>
</mrow>
</msub>
</mrow>
</mtd>
<mtd>
<mrow>
<msub>
<mi>sin&theta;</mi>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mi>n</mi>
</mrow>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>-</mo>
<msub>
<mi>sin&theta;</mi>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mi>n</mi>
</mrow>
</msub>
</mrow>
</mtd>
<mtd>
<msub>
<mi>cos</mi>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mi>n</mi>
</mrow>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mi>R</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>T</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
</mrow>
Q(θi+n,t+ti+n) it is Q (θi+n, t) in each value time shift ti+n/ S, S are sample rate.
4. multi-component converted wave crack prediction method according to claim 1, it is characterised in that the calculating fracture development
The method of intensity includes：
Fracture development intensity KC=Δ T/TS1, wherein, Δ T=TS1-TS2, TS1Exist for the fast transverse wave section back wave after separation
Time, TS2For the slow shear-wave section back wave existence time after separation.
5. multi-component converted wave crack prediction method according to claim 1, it is characterised in that described to pending sampling
Point, which carries out fast transverse wave and the method for slow shear-wave separation, includes Alford spinning solutions or triangular transformation method.
6. multi-component converted wave crack prediction method according to claim 1, it is characterised in that the anglec of rotation threshold value
Scope is [- π, π] or [0,2 π], and the delay threshold scope is empirical value or set-point.
7. multi -components according to claim 1 convert ripple crack prediction method, it is characterised in that the Forecasting Methodology is also wrapped
Include and fracture development synthesis distribution map is drawn after the geological data calculates the step of finishing.
8. multi-component converted wave crack prediction method according to claim 7, it is characterised in that the drafting fracture development
The method of comprehensive distribution map includes：The fracture development of each sampled point represents with short-term, the orientation references fracture development of short-term
Azimuth, the length of short-term represent fracture development intensity.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711122137.4A CN107894616B (en) | 2017-11-14 | 2017-11-14 | Multi-component converted wave crack prediction method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711122137.4A CN107894616B (en) | 2017-11-14 | 2017-11-14 | Multi-component converted wave crack prediction method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107894616A true CN107894616A (en) | 2018-04-10 |
| CN107894616B CN107894616B (en) | 2020-01-10 |
Family
ID=61804402
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201711122137.4A Active CN107894616B (en) | 2017-11-14 | 2017-11-14 | Multi-component converted wave crack prediction method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107894616B (en) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112083485A (en) * | 2019-06-14 | 2020-12-15 | 中国石油天然气集团有限公司 | Oil gas distribution detection method and device |
| CN113009572A (en) * | 2021-02-23 | 2021-06-22 | 成都理工大学 | Method for predicting fracture azimuth angle based on transverse wave polarization analysis |
| CN115903024A (en) * | 2022-12-26 | 2023-04-04 | 成都理工大学 | Shear wave splitting analysis method based on gradient descent method |
| CN115980852B (en) * | 2022-12-29 | 2023-08-15 | 成都理工大学 | Efficient calculation method and system for fracture parameters based on transverse wave splitting |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103424776A (en) * | 2013-08-16 | 2013-12-04 | 中国石油大学(华东) | Carbonatite oil and gas reservoir crack earthquake detection method |
| CN105116448A (en) * | 2015-08-11 | 2015-12-02 | 中国石油天然气集团公司 | Conversion wave azimuthal anisotropy correction method and device thereof |
| CN106291673A (en) * | 2015-05-18 | 2017-01-04 | 中国石油化工股份有限公司 | Crack based on shear-wave birefringence attribute factor extracting method and device |
| CN106443775A (en) * | 2016-05-25 | 2017-02-22 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | High resolution converted wave crack prediction method |
-
2017
- 2017-11-14 CN CN201711122137.4A patent/CN107894616B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103424776A (en) * | 2013-08-16 | 2013-12-04 | 中国石油大学(华东) | Carbonatite oil and gas reservoir crack earthquake detection method |
| CN106291673A (en) * | 2015-05-18 | 2017-01-04 | 中国石油化工股份有限公司 | Crack based on shear-wave birefringence attribute factor extracting method and device |
| CN105116448A (en) * | 2015-08-11 | 2015-12-02 | 中国石油天然气集团公司 | Conversion wave azimuthal anisotropy correction method and device thereof |
| CN106443775A (en) * | 2016-05-25 | 2017-02-22 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | High resolution converted wave crack prediction method |
Non-Patent Citations (2)
| Title |
|---|
| 刘红爱 等: "基于Alford旋转的转换波各向异性校正技术", 《天然气工业》 * |
| 刘红爱: "转换横波***分析及校正技术研究", 《中国优秀硕士学位论文全文数据库 基础科学辑》 * |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112083485A (en) * | 2019-06-14 | 2020-12-15 | 中国石油天然气集团有限公司 | Oil gas distribution detection method and device |
| CN112083485B (en) * | 2019-06-14 | 2024-03-01 | 中国石油天然气集团有限公司 | Oil gas distribution detection method and device |
| CN113009572A (en) * | 2021-02-23 | 2021-06-22 | 成都理工大学 | Method for predicting fracture azimuth angle based on transverse wave polarization analysis |
| CN115903024A (en) * | 2022-12-26 | 2023-04-04 | 成都理工大学 | Shear wave splitting analysis method based on gradient descent method |
| CN115903024B (en) * | 2022-12-26 | 2023-08-15 | 成都理工大学 | Transverse wave splitting analysis method based on gradient descent method |
| CN115980852B (en) * | 2022-12-29 | 2023-08-15 | 成都理工大学 | Efficient calculation method and system for fracture parameters based on transverse wave splitting |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107894616B (en) | 2020-01-10 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Eisner et al. | Uncertainties in passive seismic monitoring | |
| CN106094029B (en) | Utilize the method for offset distance vector piece geological data Predicating Reservoir Fractures | |
| CN103076623B (en) | Crack detection method based on prestack coherence | |
| CN101551466B (en) | Method for improving prediction precision of oil and gas reservoir by using seismic attribute related to offset distance | |
| CN107678063B (en) | A kind of multi-component converted wave crack prediction method based on Rank correlation | |
| CN107894616A (en) | Multi-component converted wave crack prediction method | |
| Shih et al. | Observation of shear wave splitting from natural events: South Moat of Long Valley Caldera, California, June 29 to August 12, 1982 | |
| She et al. | Shallow crustal structure of the middle‐lower Yangtze River region in eastern China from surface‐wave tomography of a large volume airgun‐shot experiment | |
| CN104853822A (en) | Method for evaluating shale gas reservoir and searching sweet spot region | |
| CN104977618A (en) | Method for evaluating shale gas reservoir and finding dessert area | |
| CN105629303A (en) | Prestack crack quantitative forecast method and system based on rock physics | |
| CN107045143A (en) | Method and device for predicting crack development | |
| CN104678434A (en) | Method for predicting storage layer crack development parameters | |
| CN102721979B (en) | Seismic data-based thin layer automatic interpretation and thickness prediction method and device | |
| CN106556861A (en) | A kind of azimuthal AVO inversion method based on Omnibearing earthquake auto data | |
| CN103149592A (en) | Method for separating variable offset vertical seismic profile (VSP) wave fields | |
| CN103257363A (en) | Method for detecting inclination angle of fissure in underground fissure type reservoir stratum | |
| CN104570086B (en) | Method for predicting pre-stack cracks in common offset and common azimuth angle domain | |
| CN103869362A (en) | Method and equipment for obtaining body curvature | |
| Zhang et al. | Shallow structure of the Longmen Shan fault zone from a high‐density, short‐period seismic array | |
| CN105445787A (en) | Crack prediction method for preferred orientation daughter coherence | |
| CN103869366A (en) | Method and device for determining fracture strike | |
| Xu et al. | Joint microseismic moment-tensor inversion and location using P-and S-wave diffraction stacking | |
| Yang et al. | Analysis of quality factors for Rayleigh channel waves | |
| ZHAO | Application of multi-component seismic exploration in the exploration and production of lithologic gas reservoirs |
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 | ||
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20201123 Address after: 100007 Beijing, Dongzhimen, North Street, No. 9, No. Patentee after: CHINA NATIONAL PETROLEUM Corp. Patentee after: BGP Inc., China National Petroleum Corp. Patentee after: CNPC Chuanqing Drilling Engineering Co.,Ltd. Address before: 610213, No. 1-4, Huayang Road, Huayang Town, Chengdu, Sichuan, Shuangliu County Patentee before: CNPC Chuanqing Drilling Engineering Co.,Ltd. Patentee before: China National Petroleum Corp. Patentee before: CNPC Chuanqing Drilling Engineering Co.,Ltd. |