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 directionOr
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, tiIs initiated byTime delay, n being 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. The loop calculation is started from sample point A (p, q), where p represents the line number of sample point A and q represents the track number of sample point A, and the radial direction isIntercepting the component data and the transverse component data according to the analysis window range, and respectively storing the intercepted data into an array R ═ R1,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 matrixWherein,
m represents P or Q, N represents P or Q, i.e.
ρ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 setting the fracture strike azimuth angle β to the radial component forward azimuth angle α at an angle θ - α,
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,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.