CN107894616B - Multi-component converted wave crack prediction method - Google Patents
Multi-component converted wave crack prediction method Download PDFInfo
- Publication number
- CN107894616B CN107894616B CN201711122137.4A CN201711122137A CN107894616B CN 107894616 B CN107894616 B CN 107894616B CN 201711122137 A CN201711122137 A CN 201711122137A CN 107894616 B CN107894616 B CN 107894616B
- Authority
- CN
- China
- Prior art keywords
- component
- transverse
- seismic data
- theta
- wave
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 46
- 238000011161 development Methods 0.000 claims abstract description 37
- 238000005070 sampling Methods 0.000 claims abstract description 31
- 238000004458 analytical method Methods 0.000 claims abstract description 16
- 238000004364 calculation method Methods 0.000 claims abstract description 12
- 238000006243 chemical reaction Methods 0.000 claims abstract description 10
- 206010017076 Fracture Diseases 0.000 claims description 36
- 208000010392 Bone Fractures Diseases 0.000 claims description 35
- 239000011159 matrix material Substances 0.000 claims description 9
- 239000000126 substance Substances 0.000 claims description 5
- 238000011426 transformation method Methods 0.000 claims description 4
- 238000013507 mapping Methods 0.000 claims 1
- 238000004422 calculation algorithm Methods 0.000 abstract description 2
- 230000008569 process Effects 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 4
- 239000004215 Carbon black (E152) Substances 0.000 description 3
- 229930195733 hydrocarbon Natural products 0.000 description 3
- 150000002430 hydrocarbons Chemical class 0.000 description 3
- 230000010287 polarization Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 239000013598 vector Substances 0.000 description 3
- 238000003491 array Methods 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 239000003208 petroleum Substances 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 230000001902 propagating effect Effects 0.000 description 2
- 238000007619 statistical method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 230000006872 improvement Effects 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
- 238000005215 recombination Methods 0.000 description 1
- 230000006798 recombination Effects 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 230000003068 static effect 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
- 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 the field of petroleum exploration, relates to a fracture system detection method for a fractured oil and gas reservoir, and particularly relates to a multi-component converted wave fracture prediction method.
Background
From the perspective of field statistics developed in the world, fractured reservoirs account for approximately half of the total world volume in terms of hydrocarbon resources and production capacity. Therefore, fractured hydrocarbon reservoirs are the focus of the geophysical exploration and research of petroleum nowadays. How to accurately judge the existence state of the fracture in the underground has a crucial role in reservoir evaluation, block evaluation and oil and gas recovery improvement. Most of the current common fracture prediction methods are methods based on post-stack seismic attributes, such as coherence attributes, curvature attributes, texture attributes and the like, but the methods based on the post-stack seismic attributes are mostly suitable for predicting large-scale fracture zones, pre-stack information and anisotropic information in the stratum cannot be fully utilized, and the prediction result is not fine enough. Therefore, a more elaborate crack prediction technique needs to be explored.
With the continuous and deep multi-wave research, the multi-wave exploration technology is increasingly applied to oil and gas exploration and development. Shear wave splitting is a powerful technique in the field of natural seismic and oil exploration and is of great interest to geologists and geophysicists. Transverse wave splitting, also called transverse wave birefringence, is that when a transverse wave passes through an azimuthal anisotropic medium, the transverse wave is split into a fast transverse wave propagating in a direction parallel to a fracture and a slow transverse wave propagating in a direction perpendicular to the fracture, and the propagation speeds of the fast transverse wave and the slow transverse wave are different, so that a fast transverse wave time difference, an amplitude attenuation and the like are formed in the anisotropic medium, and the characteristics are related to anisotropic properties. The azimuth and the density of the fracture system can be researched reversely according to the time difference, the waveform, the amplitude attenuation and the like of the fast transverse wave and the slow transverse wave. Predicting reservoir fractures using shear wave splitting becomes a straightforward and reliable method. However, shear wave sources are expensive, and large-scale application of shear wave splitting exploration by using pure shear wave sources is difficult. The converted wave exploration overcomes the defects of difficult excitation, high cost, large static correction value and the like of pure transverse wave exploration, has the advantages of rich information, combination of the advantages of longitudinal waves and transverse waves, high information signal-to-noise ratio and the like, and enables the multicomponent converted wave exploration to become a powerful tool for oil and gas reservoir exploration.
In the converted wave exploration fracture prediction, the accurate calculation of the fracture development direction and the time difference has important significance for accurately depicting the underground distribution form of a fracture system. Currently, the commonly used methods are a cross correlation method and an energy ratio method. In the correlation or covariance matrix calculated in the traditional method, a linear correlation mode is adopted, the similarity between waveforms is nonlinear, the data sample capacity calculated by parameters is limited, the characteristics of normal distribution are not met, and the sample data per se can not meet the requirements in the mathematical sense.
Disclosure of Invention
In view of the deficiencies in the prior art, it is an object of the present invention to address one or more of the problems in the prior art as set forth above.
In order to achieve the above object, the present invention provides a multi-component converted wave fracture prediction method, which may include the following steps: selecting an analysis window containing a target depth range, and intercepting the radial component seismic data and the transverse component seismic data of the sampling point to be processed according to the analysis window; performing angle rotation on the intercepted radial component seismic data and transverse component seismic data, and performing time delay conversion on the transverse component seismic data 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 from slow transverse waves of the sampling points to be processed, and calculating the development strength of the crack according to the travel time difference of the fast transverse waves and the slow transverse waves; repeating the step of selecting the analysis window to the step of calculating the crack development strength, and circulating to the next sampling point to be processed until the seismic data is calculated, wherein the step of constructing the target function and calculating the minimum value of the target function comprises the following steps:
recording the component of the intercepted radial component seismic data after angle rotation as P (theta)i+nT), the component of the time delay conversion of the transverse component seismic data after the angle rotation is Q (theta)i+n,t+ti+n) (ii) a Respectively to P (theta)i+n,t)、Q(θi+n,t+ti+n) Performing a ranking operation, which is recorded as rgP and rgQ, wherein thetai+nDenotes the angle of rotation, t denotes the time, ti+nRepresenting the time delay, n representing an integer; constructing matrices
Wherein the content of the first and second substances,
m represents P or Q, N represents P or Q, PrgM,rgNRepresenting Pearson correlation coefficients, cov (rgM, rgN) representing the covariance function, σrgPDenotes the standard deviation, σ, of rgPrgQRepresents a standard deviation of rgQ; computing matrix U (θ)i+n,t+ti+n) Characteristic value eig (U (theta))i+n,t+ti+n) ); wherein, thetai+n＝θi+nΔθ，ti+n＝ti+nΔt，θmin≤θi+n≤θmaxAnd t ismin≤ti+n≤tmax，θiTo start the rotation angle, tiFor initial delay, Δ θ is the rotation angle cycle step, Δ t is the delay cycle step, [ θ [ O ] ]min,θmax]Is a threshold range of rotation angle, [ tmin,tmax]Is a delay threshold range; construction ofThe objective function F (θ, t) is min [ eig (U (θ)i+n,t+ti+n))]N starts to take a value from 0, n increases by 1 each time to obtain the minimum value of the objective function, wherein min [ eig (U (theta) ]i+n,t+ti+n))]The expression finds the eigenvalue eig (U (theta)i+n,t+ti+n) ) of the minimum value.
Compared with the prior art, the method introduces a non-parameter statistical theory into the multi-wave fracture prediction process, adopts a level correlation mode to replace a linear correlation method in the traditional algorithm, is more sensitive to the measure of waveform similarity, avoids the requirement of linear correlation on sample data, improves the calculation accuracy and provides a basis for fine description of the fractured oil-gas reservoir.
Drawings
The above and other objects and features of the present invention will become more apparent from the following description taken in conjunction with the accompanying drawings, in which:
FIG. 1 shows a flow diagram of a multi-component converted wave fracture prediction method according to an exemplary embodiment of the invention.
FIG. 2 shows a schematic diagram of shear wave splitting vs. fracture development orientation in accordance with an exemplary embodiment of the present invention.
Detailed Description
Hereinafter, a multicomponent converted wave fracture prediction method according to the present invention will be described in detail with reference to the accompanying drawings and exemplary embodiments.
Specifically, the multi-component converted wave can generate a splitting phenomenon in an EDA medium, and fracture development zone prediction can be carried out by using the converted shear wave, so that a basis is provided for fine description of a fractured hydrocarbon reservoir. The converted transverse wave is excited by longitudinal wave in field exploration, and seismic wave is obliquely incident to the elastic interface to generate reflected transverse wave, wherein the uplink transverse wave is converted transverse wave. The converted transverse wave has all properties of the transverse wave, and when the upward converted transverse wave is obliquely crossed with the trend of the crack, the converted transverse wave is also split into a fast transverse wave and a slow transverse wave. And the split fast and slow waves are transmitted to the earth surface, decomposed and recombined according to an earth surface observation system coordinate system, and received by the multi-component detector. Therefore, when the measuring line is not parallel to the crack direction, the ground horizontal component records contain fast and slow wave information, then fast and slow wave separation is carried out by utilizing multi-component records, the information such as the direction, the density, the time difference and the like of the crack system is researched reversely according to the time difference, the waveform, the amplitude attenuation and the like of the fast wave and the slow wave, and the crack development is predicted according to the information.
The method is characterized in that a non-parameter statistical method is introduced into a crack prediction method, rank operation is carried out on converted data when an objective function is constructed, then a covariance matrix based on the rank is constructed, linear correlation is changed into nonlinear correlation, the calculation accuracy is improved, and a more reliable basis is provided for underground crack identification and prediction.
FIG. 1 shows a flow diagram of a multi-component converted wave fracture prediction method according to an exemplary embodiment of the invention. FIG. 2 shows a schematic diagram of shear wave splitting vs. fracture development orientation in accordance with an exemplary embodiment of the present invention.
The invention provides a multi-component converted wave fracture prediction method, as shown in fig. 1, in an exemplary embodiment of the invention, the method may include:
step S01, rotating the horizontal components X and Y of the converted wave three-dimensional three-component seismic data to obtain the converted wave radial component seismic data R (namely radial component R) and transverse component seismic data T (namely transverse component T) after rotation.
In the above, three-component receivers, respectively the X, Y horizontal component receiver and the Z vertical component receiver, may be used in three-dimensional three-component seismic surveys. Because the Inline and CrossLine directions are not consistent with the polarization direction of the transverse wave, the transverse wave recorded in the X, Y direction has both P-SV waves and P-SH waves, and has no clear polarization meaning, and wave fields in all polarization directions are mixed with each other, which is not beneficial to processing of converted wave data and research and extraction of transverse wave splitting information. Therefore, in order to obtain a consistent converted wave field, when three-dimensional three-component data is processed, the energy of two components in the horizontal direction needs to be redistributed, and the converted wave coordinate is rotated to realize the rotation from a field acquisition coordinate system (X-Y coordinate system) to a radial-transverse coordinate system (R-T coordinate system). After coordinate rotation processing, the converted wave has a main energy distribution on a radial component, and a transverse component mainly represents the influence of anisotropy. The coordinate rotation type may be:
wherein, the included angle between the radial component R and the X component direction
Or
The rotated converted wave may be subjected to an anisotropic treatment operation, which is a conventional method in the art.
And step S02, determining the range of an analysis window, and intercepting the radial component seismic data R and the transverse component seismic data T of the sampling point to be processed according to the analysis window.
In the present exemplary embodiment, the selection of the analysis window needs to include a target depth range, that is, a fractured anisotropic formation, and the target depth range may be controlled by a horizon control or a time period. And after the analysis window is determined from the sampling points to be processed, intercepting the radial component seismic data and the transverse component seismic data of all the sampling points respectively according to the determined analysis window range, and respectively storing the intercepted data into different arrays so as to perform subsequent calculation.
In the above, in the same working area, the horizon ranges selected by the analysis window ranges of all sampling points are consistent, for example, through horizon control, the time window ranges of all sampling points all include horizon 1 and horizon 2, but the range value corresponding to each sampling point is different in horizon 1 and horizon 2, for example, the range of sampling point C may take 1200ms (horizon 1) to 1300ms (horizon 2), and the range of sampling point D may take 1300ms (horizon 1) to 1400ms (horizon 2).
And step S03, performing rotation angle and time delay conversion on the radial component seismic data and the transverse component seismic data to obtain converted components.
In the present exemplary embodiment, since the data recorded in the radial component and the lateral component is a recombination of the fast shear wave and the slow shear wave information at the time of surface reception, that is, a superposition of the information of both the fast shear wave and the slow shear wave in the radial component and the lateral component. Therefore, it is necessary to separate the fast and slow shear waves by means of rotation. Meanwhile, the fast transverse wave and the slow transverse wave are obtained by splitting the same converted transverse wave, the waveforms of the converted transverse wave and the converted transverse wave are similar, the propagation speed of the slow transverse wave is lower than that of the fast transverse wave, and the fast transverse wave and the slow transverse wave have time difference. Therefore, delay compensation is required for slow shear waves based on the angular rotation.
Specifically, firstly, the angular rotation is performed on the intercepted radial component R and the transverse component T at the sampling point to be processed, and the rotation type may be:
the component P (theta) is obtainedi+nT) and Q (θ)i+nT) and then time-delaying the rotated transverse component data, Q (θ)i+n,t+ti+n) Is Q (theta)i+nEach value in t) is time-shifted up by ti+nPoint of/S, S is the sampling rate, thetai+nIs a rotation angle, ti+nFor time delay, θiTo start the rotation angle, tiFor initial delay, n denotes an integer, θi+n＝θi+nΔθ，ti+n＝ti+nΔt，θmin≤θi+n≤θmaxAnd t ismin≤ti+n≤tmaxWhere Δ θ is the rotation angle circulation step, Δ t is the delay circulation step, [ θ [ [ theta ]min,θmax]Is a threshold range of rotation angle, [ tmin,tmax]Is a delay threshold range.
Above, Q (θ)i+n,t+ti+n) Is Q (theta)i+nEach value in t) is time-shifted up by ti+nthe/S point refers to, for example, a time shift ti+nTaking 1 point of/S, Q (theta)i+nData set Q in t)1,Q2,Q3,Q4.. time-shift by one point, Q (θ)i+n,t+ti+n) The corresponding array becomes from Q2Is started, namely Q2,Q3,Q4...。
Starting from an initial angle and an initial time delay, namely an initial value, and obtaining a component P (theta) after the angular rotation of the intercepted radial component seismic data through an angular rotation type and a time delay calculation formulai+nT) and the component after time delay conversion of the transverse component seismic data after the angular rotation is Q (theta)i+n,t+ti+n)。
Here, the time delay calculation is mainly for the rotated lateral component, i.e., the diagonally-rotated lateral component seismic data P (θ)i+nT) the component after time delay conversion is Q (theta)i+n,t+ti+n). Angle of rotation thetai+nAnd a time delay ti+nValues are required to be taken within a threshold range, which can be a set value or an empirical value. For example, the threshold range for the rotation angle may typically be [ - π, π]Or [0,2 π]。
For example, assuming the sampling rate of the processed multi-component data is S, the analysis window is controlled by horizon hor1At the top of the window, hor2Is the bottom of the timing window. And circularly calculating from a sampling point A (p, q), wherein p represents the line number of the sampling point A, q represents the track number of the sampling point A, intercepting the radial component data and the transverse component data according to the analysis window range, and respectively storing the intercepted data into an array R [ R ] respectively1,r2,r3,......,rn]And array T ═ T1,t2,t3,......,tn]To perform the calculation.
Then, θ is performed on the arrays R and Ti+nThe rotation of the angle may be set to a threshold range of 0 < θ, for examplei+nLess than pi; then, time delay compensation is carried out on the rotated transverse component data, and the threshold range of the time delay can be set to 0, T]Obtaining the angular rotated component P (theta) of the intercepted radial component seismic datai+nT), the component of the transverse component seismic data after the angle rotation is subjected to time delay conversion is Q (theta)i+n,t+ti+n)。
In step S04, an objective function is constructed, and the minimum value of the objective function is obtained.
In the present exemplary embodiment, the construction of the objective function may include:
in step S04.1, first, component P (θ) is subjected toi+nT) and Q (θ)i+n,t+ti+n) Performing a ranking operation, i.e. separately obtaining the component P (theta)i+nT) and Q (θ)i+n,t+ti+n) Is recorded as rgP and rgQ, where θ isi+nDenotes the angle of rotation, t denotes the time, ti+nRepresenting the time delay and n representing an integer.
Step S04.2, constructing a matrix
Wherein the content of the first and second substances,
ρrgM,rgNRepresenting Pearson correlation coefficients, cov (rgM, rgN) representing the covariance function, σrgPDenotes the standard deviation, σ, of rgPrgQThe standard deviation of rgQ is shown.
Step S04.3, calculating a matrix U (theta)i+n,t+ti+n) Characteristic value eig (U (theta))i+n,t+ti+n)). Wherein, thetai+n＝θi+nΔθ，ti+n＝ti+nΔt，θmin≤θi+n≤θmaxAnd t ismin≤ti+n≤tmax，θiTo start the rotation angle, tiFor initial delay, Δ θ is the rotation angle cycle step, Δ t is the delay cycle step, [ θ [ O ] ]min,θmax]Is a threshold range of rotation angle, [ tmin,tmax]Is a delay threshold range.
Here, it is necessary to determine the rotation angle θi+nAnd a time delay ti+nIf the value is within the threshold range, if the threshold condition is satisfied, the new rotation angle and the new time delay value can be continuously and circularly taken, and the step S03 is carried out, so that different components P (theta) are obtainedi+nT) and Q (θ)i+n,t+ti+n) And then different characteristic values are obtained.
In step S04.4, an objective function F (θ, t) is constructed as min [ eig (U (θ)i+n,t+ti+n))]N starts to take a value from 0, n increases by 1 each time to obtain the minimum value of the objective function, wherein min [ eig (U (theta) ]i+n,t+ti+n))]The expression finds the eigenvalue eig (U (theta)i+n,t+ti+n) ) of the minimum value.
In step S05, the rotation angle and the time shift corresponding to the minimum value of the objective function F (θ, t) are recorded.
Here, when the objective function takes the minimum value, the rotation angle at this time is the azimuth angle of crack development, and the time delay is the time difference of fast and slow waves. When the conversion wave is split, the larger the difference between the fast wave and the slow wave is, the more the crack develops, therefore, the difference between the fast wave and the slow wave can be used as an important index of the crack development degree.
And step S06, separating fast transverse waves and slow transverse waves of the sampling points to be processed, and calculating the development strength of the crack according to the travel time difference of the fast transverse waves and the slow transverse waves.
In the present exemplary embodiment, the radial component seismic data and the transverse component seismic data each contain fast and slow shear wave information, and therefore, it is necessary to separate fast and slow shear waves from the radial component seismic data and the transverse component seismic data.
When the crack development at the point to be sampled is confirmed, a trigonometric transformation method or an Alford rotation method can be adopted. The triangular transformation method and the Alford rotation method are used for separating the transverse component seismic data and the radial component seismic data at the crack in a coordinate rotation mode to obtain fast transverse waves and slow transverse waves.
The above-mentioned trigonometric transformation method is adoptedObtaining a fast transverse wave and a slow transverse wave, wherein S1Representing fast transverse waves, S2Representing slow shear waves.
The basic idea of the Alford rotation method may include: an included angle between the fracture strike azimuth angle beta and the radial component forward azimuth angle alpha is set as beta-alpha,
defining the orthogonal rotation matrix G as:
constructing an acceptance signal matrix:
wherein R is1、T1And R2、T2The two azimuth angles are respectively a radial vector and a transverse vector of two different azimuth angles, and the two azimuth angles meet the orthogonal relation.
According to the Alford theory of rotation: V-RSRT，
Wherein the content of the first and second substances,
S1representing fast transverse waves, S2Which represents a slow transverse wave, is shown,
pair formula V ═ RSRTAfter orthogonal rotation, S ═ R is obtainedTVR。
The method for calculating the fracture development strength comprises the following steps: crack development Strength KC＝ΔT/TS1Wherein Δ T ═ TS1-TS2，TS1The time of existence of reflected wave in the separated fast transverse wave section, TS2The existing time of the reflected wave of the separated slow shear wave section is shown.
And step S07, repeating the steps S01 to S06, and circulating to the next sample point to be processed until the seismic data are calculated.
Here, when the next sample point to be processed is processed, the line number and track number of the last sample point to be processed need to be determined, for example, assuming that a sample point a (p, q) is satisfied, when start _ inline < p < end _ inline, and start _ crossline < q < end _ crossline, then the method may cycle to the next sample point to be processed until the seismic data is completely calculated, where p is the line number of the sample point to be processed, q is the track number of the sample point to be processed, start _ inline is the work area start line number, end _ inline is the work area end line number, start _ crossline is the work area start track number, and end _ crossline is the work area end track number.
It should be noted that, when sampling points are processed in a cycle, the cycle sequence of the sampling points may be cycled according to the line number and the track number. For example, during the circulation, a certain line number p is kept unchanged, the circulation is performed in sequence according to the track number in the track number range of start _ crossline < q < end _ crossline, after the track number circulation is completed, the line number is changed into p +1, and then all the track numbers are circulated in sequence.
And step S08, drawing a comprehensive fracture development distribution map.
After all sampling points are calculated, the data result of each sampling point is stored to form two data, namely a crack development azimuth angle and crack development strength, and then a crack development comprehensive distribution map is drawn in the comprehensive crack development azimuth and crack development strength. The method for drawing the comprehensive fracture development distribution map comprises the following steps: the crack development of each sampling point is represented by short lines, the azimuth of the short lines represents the azimuth of the crack development, and the length of the short lines represents the crack development strength.
The relationship between shear wave splitting and crack development orientation is shown in FIG. 2. After splitting, the fast transverse wave is along the crack development direction, the slow transverse wave is perpendicular to the crack development direction, and according to the wave superposition principle, the radial component and the transverse component data both contain fast and slow wave components.
In conclusion, the invention introduces the non-parametric statistical method into the multi-component converted wave fracture prediction, avoids the requirement of the traditional linear correlation method on the seismic data distribution, replaces the grade correlation calculation result when carrying out covariance matrix calculation, can well describe the correlation between two vectors, improves the calculation accuracy and ensures that the prediction result is more accurate and reliable.
Although the present invention has been described above in connection with exemplary embodiments, it will be apparent to those skilled in the art that various modifications and changes may be made to the exemplary embodiments of the present invention without departing from the spirit and scope of the invention as defined in the appended claims.
Claims (8)
1. A multi-component converted wave fracture prediction method is characterized by comprising the following steps:
selecting an analysis window containing a target depth range, and intercepting the radial component seismic data and the transverse component seismic data of the sampling point to be processed according to the analysis window;
performing angle rotation on the intercepted radial component seismic data and transverse component seismic data, and performing time delay conversion on the transverse component seismic data 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 from slow transverse waves of the sampling points to be processed, and calculating the development strength of the crack according to the travel time difference of the fast transverse waves and the slow transverse waves;
repeating the steps from the step of selecting the analysis window to the step of calculating the crack development strength, circulating to the next sampling point to be processed until the seismic data is calculated, wherein,
the step of constructing the objective function and calculating the minimum value of the objective function comprises the following steps:
recording the component of the intercepted radial component seismic data after angle rotation as P (theta)i+nT), after angular rotationThe component of the transverse component seismic data subjected to time delay conversion is Q (theta)i+n,t+ti+n)；
Respectively to P (theta)i+n,t)、Q(θi+n,t+ti+n) Performing a ranking operation, which is recorded as rgP and rgQ, wherein thetai+nDenotes the angle of rotation, t denotes the time, ti+nRepresenting the time delay, n representing an integer;
constructing matrices
Wherein the content of the first and second substances,
computing matrix U (θ)i+n,t+ti+n) Characteristic value eig (U (theta))i+n,t+ti+n))；
Wherein, thetai+n＝θi+nΔθ，ti+n＝ti+nΔt，θmin≤θi+n≤θmaxAnd t ismin≤ti+n≤tmax，θiTo start the rotation angle, tiFor initial delay, Δ θ is the rotation angle cycle step, Δ t is the delay cycle step, [ θ [ O ] ]min,θmax]Is a threshold range of rotation angle, [ tmin,tmax]Is a delay threshold range;
constructing an objective function F (theta, t) ═ min [ eig (U (theta) ]i+n,t+ti+n))]N starts to take a value from 0, n increases by 1 each time to obtain the minimum value of the objective function, wherein min [ eig (U (theta) ]i+n,t+ti+n))]The expression finds the eigenvalue eig (U (theta)i+n,t+ti+n) ) of the minimum value.
2. The multi-component converted wave fracture prediction method of claim 1, wherein the radial component seismic data and the transverse component seismic data are obtained by rotating horizontal components X and Y of three-dimensional three-component seismic data, and the coordinate rotation type of the rotation of the horizontal components X and Y is as follows:
wherein the content of the first and second substances,
is the angle of the radial component to the horizontal component X,
or
R represents radial component seismic data and T represents transverse component seismic data.
3. The multi-component converted wave fracture prediction method of claim 1, wherein the component P (θ)i+nT) and Q (θ)i+n,t+ti+n) The calculation formula of (A) is as follows:
Q(θi+n,t+ti+n) Is Q (theta)i+nTime-shifting each value in t) by ti+nAnd S is the sampling rate.
4. The multi-component converted wave fracture prediction method of claim 1, wherein the method of calculating fracture development strength comprises:
crack development Strength KC＝ΔT/TS1Wherein Δ T ═ TS1-TS2，TS1The time of existence of reflected wave in the separated fast transverse wave section, TS2The existing time of the reflected wave of the separated slow shear wave section is shown.
5. The multi-component converted wave fracture prediction method of claim 1, wherein the method for separating fast shear waves from slow shear waves of the sample points to be processed comprises an Alford rotation method or a triangular transformation method.
6. The multi-component converted wave fracture prediction method of claim 1, wherein the rotation angle threshold range is [ -pi, pi ] or [0,2 pi ], and the time delay threshold range is an empirical value or a given value.
7. The method of claim 1, further comprising the step of generating a synthetic profile of fracture development after the step of computing the seismic data.
8. The multi-component converted wave fracture prediction method of claim 7, wherein the method of mapping the fracture development comprehensive profile comprises: the crack development of each sampling point is represented by short lines, the azimuth of the short lines represents the azimuth of the crack development, and the length of the short lines represents the crack development strength.
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 CN107894616A (en) | 2018-04-10 |
| CN107894616B true 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) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112083485B (en) * | 2019-06-14 | 2024-03-01 | 中国石油天然气集团有限公司 | Oil gas distribution detection method and device |
| CN113009572B (en) * | 2021-02-23 | 2022-07-01 | 成都理工大学 | Method for predicting fracture azimuth angle based on transverse wave polarization analysis |
| 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 |
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旋转的转换波各向异性校正技术;刘红爱 等;《天然气工业》;20111231;第31卷(第5期);第41-43页 * |
| 转换横波分裂分析及校正技术研究;刘红爱;《中国优秀硕士学位论文全文数据库 基础科学辑》;20120415(第4期);第011-510页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107894616A (en) | 2018-04-10 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107894616B (en) | Multi-component converted wave crack prediction method | |
| CN103076623B (en) | Crack detection method based on prestack coherence | |
| NO315767B1 (en) | Procedure for Processing Interface-Reflected Cormorants from an Azimuthal-Anisotropic Geological Formation | |
| CN101551466A (en) | Method for improving prediction precision of oil and gas reservoir by using seismic attribute related to offset distance | |
| CN105629303A (en) | Prestack crack quantitative forecast method and system based on rock physics | |
| CN105425292A (en) | Oil and gas prediction method and oil and gas prediction device | |
| CN102967883A (en) | Method for predicting rock brittleness probability through pre-stack elasticity parameter inversion of shale gas | |
| CN103645503A (en) | Three-dimensional time domain illumination analysis and amplitude compensation method | |
| CN103869362A (en) | Method and equipment for obtaining body curvature | |
| CN104237936B (en) | A kind of frequency of oil and gas detection becomes inversion method | |
| Diaz-Acosta et al. | Investigation of fractured carbonate reservoirs by applying shear-wave splitting concept | |
| CN102053262A (en) | Method for acquiring azimuth velocity of seismic converted wave and method for processing seismic data | |
| Grossman et al. | Integration of multicomponent time-lapse processing and interpretation: Focus on shear-wave splitting analysis | |
| Yang et al. | Fracture prediction based on walkaround 3D3C vertical seismic profiling data: A case study from the Tarim Basin in China | |
| CN103869366A (en) | Method and device for determining fracture strike | |
| Zhenwu et al. | Practices and expectation of high-density seismic exploration technology in CNPC | |
| Banerjee et al. | Anisotropy and fracture analysis for coalbed methane reservoir development in Bokaro coalfield, India | |
| Yang et al. | A novel method for determining geophone orientations from zero-offset VSP data constrained by scalar field | |
| CN112684498A (en) | Reservoir fracture prediction method and system based on wide-azimuth seismic data | |
| Nadri et al. | Estimation of stress-dependent anisotropy from P-wave measurements on a spherical sample | |
| Sun et al. | The application of amplitude-preserved processing and migration for carbonate reservoir prediction in the Tarim Basin, China | |
| Bignardi et al. | Free and improved computer codes for HVSR processing and inversions | |
| Lu et al. | Inversion of Rayleigh waves using a genetic algorithm in the presence of a low-velocity layer | |
| Luo et al. | High-dimensional co-occurrence matrix: A new tool for 3D seismic target visualization and interpretation | |
| Pei et al. | Research and application of 5D seismic prediction technology |
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: GEOPHYSICAL EXPLORATION COMPANY OF CNPC CHUANQING DRILLING ENGINEERING Co.,Ltd. Patentee before: China National Petroleum Corp. Patentee before: CNPC Chuanqing Drilling Engineering Co.,Ltd. |