CN106934183B - Dispersion curve determines that method and apparatus and p-and s-wave velocity determine method and apparatus - Google Patents
Dispersion curve determines that method and apparatus and p-and s-wave velocity determine method and apparatus Download PDFInfo
- Publication number
- CN106934183B CN106934183B CN201710265841.9A CN201710265841A CN106934183B CN 106934183 B CN106934183 B CN 106934183B CN 201710265841 A CN201710265841 A CN 201710265841A CN 106934183 B CN106934183 B CN 106934183B
- Authority
- CN
- China
- Prior art keywords
- frequency
- data
- time
- velocity
- related coefficient
- 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.)
- Expired - Fee Related
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH DRILLING; MINING
- E21B—EARTH DRILLING, e.g. DEEP DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B47/00—Survey of boreholes or wells
- E21B47/12—Means for transmitting measuring-signals or control signals from the well to the surface, or from the surface to the well, e.g. for logging while drilling
- E21B47/14—Means for transmitting measuring-signals or control signals from the well to the surface, or from the surface to the well, e.g. for logging while drilling using acoustic waves
Abstract
The present invention provides dispersion curves to determine that method and device and stratum p-and s-wave velocity determine method and apparatus, is related to geological exploration field.Dispersion curve provided by the invention determines method, it first carries out Short Time Fourier Transform to the sonic data that acoustic logging obtains, to obtain multiple time-frequency numeric field datas, later, Velocity Time related operation is carried out to the multiple time-frequency numeric field data, to obtain the first related data body of characterization frequency-Velocity-time value and the relationship of the first related coefficient, finally, according in the first related data body, the second related coefficient corresponding to each frequency-velocity amplitude, determine the dispersion curve of sound wave corresponding to the sonic data, wherein, second related coefficient is the maximum value of the first related coefficient corresponding to each frequency-velocity amplitude, and then final dispersion curve is determined by way of this mathematical statistics and analysis, improve the accuracy for determining dispersion curve.
Description
Technical field
The present invention relates to geological exploration fields, and method and device, and velocity of wave in length and breadth are determined in particular to dispersion curve
Degree determines method and apparatus.
Background technology
Acoustic logging is a kind of geological exploration field commonly technological means, and sound wave mainly is utilized in different medium
When middle propagation, speed, the acoustic characteristics such as variation of amplitude and frequency differ, and to complete to log well, and then obtain acoustic logging money
Material.
The information such as formation porosity, lithology, formation pore fluid type can be obtained using Sonic Logging Data；It can be with
Foundation is provided for the processing explanation of seismic prospecting data.Certain foundation can also be provided for the explanation of seismic prospecting data, it can
See that the application of acoustic logging is quite extensive.
During acoustic logging there is frequency dispersion in received waveform mostly, when causing to carry out the Wave data of measurement
Between domain analysis when, desired data can not be accurately obtained.
Invention content
The purpose of the present invention is to provide dispersion curves to determine method, obtains the accuracy of dispersion curve.
In a first aspect, an embodiment of the present invention provides dispersion curves to determine method, including：
Multiple initial data are obtained, each initial data is by the sonic data received by different receivers；
Short Time Fourier Transform is carried out to each initial data respectively, to obtain multiple time-frequency numeric field datas；
To the multiple time-frequency numeric field data carry out Velocity Time related operation, with obtain characterization frequency-Velocity-time value with
First related data body of the relationship of the first related coefficient；
According in the first related data body, the second related coefficient corresponding to each frequency-velocity amplitude determines the sound wave
The dispersion curve of sound wave corresponding to data, wherein the second related coefficient is the first phase relation corresponding to each frequency-velocity amplitude
Several maximum values.
With reference to first aspect, an embodiment of the present invention provides the first possible embodiments of first aspect, wherein step
Suddenly Short Time Fourier Transform is carried out to each initial data respectively, includes to obtain multiple time-frequency domain data：
Short Time Fourier Transform is carried out to obtaining the initial data according to following formula,
Wherein, xn(m) it is the received sonic data of n-th of receiver, n=1,2,3 ... number for receiver；m
To represent the variable of time, ω is frequency, and w (m) is preset Hanning window function；W (x) is time-frequency numeric field data.
With reference to first aspect, an embodiment of the present invention provides second of possible embodiments of first aspect, wherein institute
The expression formula for stating preset Hanning window function isWherein, NwFor the Chinese
Peaceful window length.
With reference to first aspect, an embodiment of the present invention provides the third possible embodiments of first aspect, wherein institute
The value for stating Hanning window length is 16.
With reference to first aspect, an embodiment of the present invention provides the 4th kind of possible embodiments of first aspect, wherein step
Suddenly Velocity Time related operation is carried out to the multiple time-frequency numeric field data, to obtain characterization frequency-Velocity-time value and the first phase
The related data body of the relationship of relationship number includes：
The numerical value with the first related coefficient corresponding to each frequency-Velocity-time value is calculated according to following formula,
Wherein, N is the total number of receiver in receiver array, and the position of window when T is, t is time, n=1,2,3 ... N
For receiver serial number, v is speed, and d is receiver spacing, TwFor when window width, ω is frequency；ρ (ω, ν, T）For the first correlation
The numerical value of coefficient, XnFor time-frequency numeric field data.
With reference to first aspect, an embodiment of the present invention provides the 5th kind of possible embodiments of first aspect, wherein step
Suddenly according in the first related data body, the second related coefficient corresponding to each frequency-velocity amplitude determines the sonic data institute
The dispersion curve of corresponding sound wave includes：
Determine the second related data body of the relationship of characterization frequency-velocity amplitude and the second related coefficient；
To in the second related data body, frequency and velocity amplitude corresponding to the second related coefficient maximum value are sought, with
Obtain the dispersion curve.
With reference to first aspect, an embodiment of the present invention provides the 6th kind of possible embodiments of first aspect, wherein step
It is rapid to determine that characterization frequency-velocity amplitude and the second related data body of the relationship of the second related coefficient include：
The second related coefficient is calculated according to following formula,
Wherein, ρ0For the numerical value of the second related coefficient, the position of window when T is, ρ (ω, ν, T）For the number of the first related coefficient
Value.
Second aspect, the embodiment of the present invention additionally provide dispersion curve determining device, including：
Acquisition module, for obtaining multiple initial data, each initial data is received by different receivers
Sonic data；
Conversion module, for carrying out Short Time Fourier Transform to each initial data respectively, to obtain multiple time-frequency domain numbers
According to；
Related operation module, for carrying out Velocity Time related operation to the multiple time-frequency numeric field data, to be characterized
First related data body of frequency-Velocity-time value and the relationship of the first related coefficient；
First determining module, for according in the first related data body, second corresponding to each frequency-velocity amplitude to be related
Coefficient determines the dispersion curve of sound wave corresponding to the sonic data, wherein the second related coefficient is each frequency-velocity amplitude
The maximum value of the first corresponding related coefficient.
The third aspect, the embodiment of the present invention additionally provide stratum p-and s-wave velocity and determine method, include such as first aspect
Dispersion curve determines method, further includes：
It is worth the shear wave velocity that corresponding speed determines stratum according to the low-frequency cutoff of dispersion curve.
Fourth aspect, the embodiment of the present invention additionally provide stratum p-and s-wave velocity determining device, include such as second aspect
Dispersion curve determining device further includes：
Second determining module, for being worth the shear wave speed that corresponding speed determines stratum according to the low-frequency cutoff of dispersion curve
Degree.
Dispersion curve provided in an embodiment of the present invention determines method, and in the prior art when there are multiple for same frequency range
When the mode wave of frequency dispersion, the dispersion curve that can only access a wave is compared, the sound wave number first obtained to acoustic logging
According to Short Time Fourier Transform is carried out, to obtain multiple time-frequency numeric field datas, later, speed is carried out to the multiple time-frequency numeric field data
Time correlation operation, to obtain the first related data body of characterization frequency-Velocity-time value and the relationship of the first related coefficient,
Finally, according in the first related data body, the second related coefficient corresponding to each frequency-velocity amplitude determines the sound wave number
According to the dispersion curve of corresponding sound wave, wherein the second related coefficient is the first related coefficient corresponding to each frequency-velocity amplitude
Maximum value, and then final dispersion curve is determined by way of this mathematical statistics and analysis, improves and determine frequency
The accuracy of non-dramatic song line.
To enable the above objects, features and advantages of the present invention to be clearer and more comprehensible, preferred embodiment cited below particularly, and coordinate
Appended attached drawing, is described in detail below.
Description of the drawings
In order to illustrate the technical solution of the embodiments of the present invention more clearly, below will be to needed in the embodiment attached
Figure is briefly described, it should be understood that the following drawings illustrates only certain embodiments of the present invention, therefore is not construed as pair
The restriction of range for those of ordinary skill in the art without creative efforts, can also be according to this
A little attached drawings obtain other relevant attached drawings.
Fig. 1 shows that the dispersion curve that the embodiment of the present invention is provided determines in method, with brill dipole array cement bond logging
Well theoretical modeling waveform；
Fig. 2-9 respectively illustrates the dispersion curve that the embodiment of the present invention is provided and determines in method, is connect to the first to eight
Receive the oscillogram that the sonic data received by device carries out the time-frequency numeric field data obtained after Short Time Fourier Transform；
Figure 10 shows that the dispersion curve that the embodiment of the present invention is provided determines in method, is carried out to waveform shown in Fig. 1 short
When Fourier transformation after carry out Velocity Time related operation again after, obtained related coefficient 3D data volume (the first related data
Body) paraphrase figure；
Figure 11-15 respectively illustrates the dispersion curve that the embodiment of the present invention is provided and determines in method, and frequency content is
1996.1Hz, 3992.2Hz, 5988.3Hz, 8982.4Hz figure related to time difference time of component when 10979Hz；
Figure 16 shows that the dispersion curve that the embodiment of the present invention is provided determines in method, to waveform shown in Fig. 1 using this
Invention carries out the result map after dispersion analysis；
Figure 17 shows the dispersion curves that the embodiment of the present invention is provided to determine in method, uses and adds to waveform shown in Fig. 1
It weighs frequency spectrum coherent method and carries out the result map after dispersion analysis；
Figure 18 shows that the dispersion curve that the embodiment of the present invention is provided determines in method, and hard formation is with brill dipole acoustic
Well logging oscillogram；
Figure 19 shows that the dispersion curve that the embodiment of the present invention is provided determines in method, to hard formation shown in Figure 17 with brill
Dipole waveform determines that method carries out the result of dispersion analysis using dispersion curve provided by the present invention；
Figure 20 shows that the dispersion curve that the embodiment of the present invention is provided determines in method, to hard formation shown in Figure 17 with brill
Dipole waveform uses the result of the frequency spectrum coherent method progress dispersion analysis in traditional technology；
Figure 21 shows that the dispersion curve that the embodiment of the present invention is provided determines in method, and the soft formation of theoretical modeling is with brill
Dipole acoustic logging oscillogram；
Figure 22 shows that the dispersion curve that the embodiment of the present invention is provided determines in method, to soft formation shown in Figure 20 with brill
Dipole waveform determines that method carries out the result of dispersion analysis using dispersion curve provided by the present invention；
Figure 23 shows that the dispersion curve that the embodiment of the present invention is provided determines in method, to soft formation shown in Figure 20 with brill
Dipole waveform uses the result of the frequency spectrum coherent method progress dispersion analysis in traditional technology；
Figure 24 shows that the dispersion curve that the embodiment of the present invention is provided determines the basic flow chart of method.
Specific implementation mode
Below in conjunction with attached drawing in the embodiment of the present invention, technical solution in the embodiment of the present invention carries out clear, complete
Ground describes, it is clear that described embodiments are only a part of the embodiments of the present invention, instead of all the embodiments.Usually exist
The component of the embodiment of the present invention described and illustrated in attached drawing can be arranged and be designed with a variety of different configurations herein.Cause
This, the detailed description of the embodiment of the present invention to providing in the accompanying drawings is not intended to limit claimed invention below
Range, but it is merely representative of the selected embodiment of the present invention.Based on the embodiment of the present invention, those skilled in the art are not doing
The every other embodiment obtained under the premise of going out creative work, shall fall within the protection scope of the present invention.
In the related technology, occurred much determining waveform using acoustic logging, and then determined acoustic speed information
Method, but lead to the waveform number to measurement the study found that the waveform that acoustic logging receives has frequency dispersion mostly by inventor
When according to the analysis for carrying out time-domain, cannot effectively obtain the velocity information of each mode wave (can do between a mode wave
It disturbs).It is, therefore, necessary to which the velocity profile of waveform can just be fully understood by carrying out dispersion analysis to the waveform received, extraction is accurate
True velocity information.
At present there are many kinds of the methods of sound wave measuring well curve dispersion analysis, it is divided into two class side of parameter Estimation and non-parametric estmation
Method.These two kinds of methods are briefly introduced separately below：
One, the method for parameter Estimation class mainly has Prony methods, pencil of matrix method and the side Prony expanded forward, backward
Method, such method are based on theoretical model, need to match data with theoretical model function；Spatial sampling point is less, processing
As a result alias is easy tod produce；Need to determine the quantity of mode wave in advance；Disadvantages mentioned above makes the stability of such methods poor, is easy
It is interfered by noise.
Two, against the above deficiency, Tang et al. proposes method-Weighted spectral coherent method of non-parametric estmation.This kind of side
Method increases the data sampling points in processing, and mould is determined by finding the maximum value of the coherence factor within the scope of specific slowness
The actual value of formula wave number.The method increase the utilization rates of data, reduce the influence of noise, overcome spatial sampling points not
The problem of foot.This method need not find pattern function, but directly be handled data to obtain dispersion curve, this kind of method
Speed is fast, and anti-noise ability is strong, is the most common method of current sonic data processing.
In practical operation, since Weighted spectral coherent method directly does correlation in frequency domain, in processing acoustic logging while drilling number
According to when, since acoustic logging while drilling same frequency range is there are the mode wave of multiple frequency dispersions, treatment effect is very poor, can only obtain certain
The dispersion curve of one wave is greatly limited in the processing of acoustic logging while drilling data frequency dispersion.
In view of this, this application provides a kind of improved dispersion curves to determine method, as shown in figure 24, this method is main
It is made of following four steps：
S101, obtains multiple initial data, and each initial data is by the sound wave number received by different receivers
According to；
S102 carries out Short Time Fourier Transform, to obtain multiple time-frequency numeric field datas to each initial data respectively；
S103 carries out Velocity Time related operation to multiple time-frequency numeric field datas, to obtain characterization frequency-Velocity-time value
With the first related data body of the relationship of the first related coefficient；
S104, according in the first related data body, the second related coefficient corresponding to each frequency-velocity amplitude determines sound
Wave number according to corresponding sound wave dispersion curve, wherein the second related coefficient is first related corresponding to each frequency-velocity amplitude
The maximum value of coefficient.
In step 1, by receiver come receive the method for initial data (i.e. sonic data) this field widely
It uses, excessive description is no longer done in the application.
Under normal conditions, initial data is shown in the form of specific waveform, and then step S102, firstly for
The receiver waveform (i.e. initial data) being made of N number of receiver that acoustic logging while drilling obtains (can receive it by (1) formula
Device waveform) Short Time Fourier Transform is carried out, to obtain the time-frequency domain waveform (instant frequency domain data) of each receiver waveform, also
It is that the waveform signal of time-domain is transformed to signal about two variables of time and frequency, instant frequency domain data.
Wherein, xn(m) it is the received sonic data of n-th of receiver, n=1,2,3 ... N are that receiver is numbered；m
To represent the variable of time, ω is frequency, and w (x) is preset Hanning window function；Xn(t, ω) is time-frequency numeric field data.
Preferably, the expression formula of preset Hanning window function is
Wherein, NwFor Hanning window length, use Hanning window that can effectively be avoided as window function progress time-frequency conversion herein transformed
Gibbs' effect in journey, can be to avoid the false waveform jittering noise that Gibbs' effect introduces.
It is furthermore preferred that the value of the Hanning window length is 16, which is to pass through a large amount of data processing practice summary
Adaptation acoustic logging while drilling signal frequency length of window.
Step S103, in step S102, signal (time-frequency numeric field data) X after Short Time Fourier Transformn(t,ω)
Velocity Time related operation is carried out by (2) formula, calculates the first related coefficient corresponding with frequency-Velocity-time value, Jin Erke
To obtain a three-dimensional frequency-Velocity-time related data body ρ (ω, v, T) to get to the first related data body.This first
For related data precursor reactant in specified frequency, specified speed, under the specified time, the numerical value of the first dependency number coefficient is big
It is small.
Wherein, N is the total number of receiver in receiver array, and the position of window when T is, t is time, n=1,2,3 ... N
For receiver serial number, v is speed, and d is receiver spacing, TwFor when window width, ω is frequency；ρ (ω, ν, T）For the first correlation
The numerical value of coefficient, XnFor time-frequency numeric field data.
Step S104 exists the three-dimensional frequency obtained in step 2-Velocity-time data volume ρ (ω, v, T) by formula (3)
The maximum value of correlation coefficient ρ (the first related coefficient) is taken at each ω and each v, then it can be by three-dimensional frequency-speed-
The first related data body of time obtains two-dimensional frequency-velocity correlation coefficint data volume ρ0 (ω, ν）(i.e. the second related data body),
The two-dimensional correlation coefficient data body is the correspondence for illustrating frequency-speed of each mode wave in array waveform, that is, is characterized
In the Dispersion Characteristics of specified frequency, specified speed, the numerical values recited of the second dependency number coefficient, that is, each mode wave,
Corresponding frequency at related coefficient local maximum in dispersion relation figure and speed are sought, you can obtain each pattern velocity of wave
Spend information with frequency change, i.e. dispersion curve v (ω).
Wherein, ρ0For the numerical value of the second related coefficient, the position of window when T is, ρ (ω, v, T) is the number of the first related coefficient
Value.
In the following, illustrating that dispersion curve provided herein determines method with a detailed example.
In scheme provided by the present invention, first by way of Short Time Fourier Transform, each receiver institute is utilized
The Wave data (sonic data) received acquires the time-frequency distributions data (time-frequency numeric field data) of array sonic log waveform, i.e.,
The acoustic wave array signal of time-domain is expressed as to the signal about frequency and time.As shown in Figure 1, being certain stratum p-and s-wave velocity
The theoretical waveforms with brill dipole array acoustic logging of respectively 3000m/s and 1400m/s, are denoted as xn(t), wherein t tables
Show time, n=1,2,3,4,5,6,7,8.Data shown in table 1- tables 8 are received by first eight receiver of receiver-the
The sonic data arrived is from left to right the sequence of time.Short Time Fourier Transform is carried out to its (sonic data) using formula (1),
Its time-frequency distributions can be obtained to get then frequency domain data, be denoted as Xn(t, ω), wherein ω indicates frequency.Table 9- tables 16 are respectively
The time-frequency domain number after Short Time Fourier Transform is carried out to the acoustic waveform received by first receiver-the, eight receivers
According in table 9- tables 16, lateral coordinates are frequency, and longitudinal coordinate is the time.Fig. 2-Fig. 9 is respectively first receiver-the 8
Receiver carries out the oscillogram of the time-frequency numeric field data after Short Time Fourier Transform, and Fig. 2-Fig. 9 is corresponding with table 9-16.In Fig. 2-Fig. 9,
Horizontal axis is the time, and the longitudinal axis is frequency.To the time-frequency domain waveform X of this 8 receiversn(t, ω) carries out time speed phase by formula (2)
Close operation, you can obtain a three-dimensional frequency-time difference-Time Correlation Data body ρ (ω, v, T), i.e. the first data volume, such as scheme
Shown in 10, the directions x reference axis indicates that time, the directions y reference axis indicate that speed, the directions z reference axis indicate frequency, color light and shade
Represent the numerical values recited of the first related coefficient.
3D data volume describes ω, v, the functional relation of these three parameters of T, and when expression, it is fixed that should fix an amount
Value since desirable numerical value is excessive, therefore is only illustrated by taking the value of following several ω as an example, take respectively ω be 1996.1Hz,
The data such as table 17- tables of related coefficient data volume ρ (ω, v, T) when 3992.2Hz, 5988.3Hz, 8982.4Hz and 10979Hz
Shown in 21, in table 17- tables 21, abscissa is the time, and ordinate is speed, and data are drawn, and related figure can be obtained as schemed
11, Figure 12, Figure 13, Figure 14, image shown in figure 15, i.e. Figure 11 are frequency-Velocity-time that waveshape obtains shown in Fig. 1
Related data body ρ (ω, v, T), frequency value are the component of 1996.1Hz, i.e. ρ (1996.1Hz, v, T)；Figure 12 is Fig. 1 institutes
The phase when frequency that oscillography shape is calculated-Velocity-time related data body ρ (ω, v, T), i.e. frequency value are 3992.2Hz
Relationship number figure, i.e. ρ (3992.2Hz, v, T)；Figure 13 is frequency-Velocity-time related data that waveshape obtains shown in Fig. 1
Related coefficient figure when body ρ (ω, v, T), i.e. frequency value are 5988.3Hz, i.e. ρ (5988.3Hz, v, T)；Figure 14 is Fig. 1 institutes
The phase when frequency that oscillography shape is calculated-Velocity-time related data body ρ (ω, v, T), i.e. frequency value are 8982.4Hz
Relationship number figure, i.e. ρ (8982.4Hz, v, T)；Figure 15 is frequency-Velocity-time related data that waveshape obtains shown in Fig. 1
Related coefficient figure when body ρ (ω, v, T), i.e. frequency value are 10979Hz, i.e. ρ (10979Hz, v, T).By that can be seen in figure
Go out, the related figure that different frequency contents obtains is different.The time-frequency of array acoustic waveform is obtained by Short Time Fourier Analysis
After distribution, then carries out the correlation analysis of each frequency content and can be obtained the velocity information of each mode wave under different frequency ingredient,
The maximum value (i.e. the maximum value of the first related coefficient) that its related coefficient is sought to the related figure of each frequency content, then can obtain
To two-dimensional frequency-velocity correlation coefficint data volume.The two-dimensional correlation coefficient data body illustrates each mould in array waveform
The correspondence of frequency-speed of formula wave, i.e., the Dispersion Characteristics of each mode wave.To related coefficient data volume ρ shown in Fig. 10
(ω, v, T) seeks the maximum value of its first related coefficient under each frequency and each speed, that is, uses as shown in formula (3)
Formula sought to get to frequency dispersion data, to get to data acquisition system as shown in table 22, in table 22, X direction is frequency
Rate, y direction are speed.Data shown in table 22 are subjected to frequency dispersion figure drafting, obtain image as shown in figure 16, i.e., shown in Fig. 1
The frequency dispersion figure with brill dipole mode wave of waveform.Intermediate dark colour is dispersion curve v (ω), by figure it is found that each
A frequency values can obtain two speed, this example medium velocity is low to be denoted as v1 for stratum dipole wave pattern, fast to be
Instrument dipole wave pattern is denoted as v2, and corresponding data are as shown in table 23.
Processing through the above steps can obtain accurate acoustic logging while drilling dispersion analysis as a result, due to brill sound
Wave log has the mode wave of multiple frequency dispersions in same frequency, can further obtain each mould based on the above method
The dispersion curve of formula wave has significant advantage.
In addition, to illustrate that the technological value of this programme, spy are directed to identical data, dispersion analysis determined using existing scheme
As a result.Below only comparative example is carried out for using frequency spectrum weighted, coherent method.
Frequency dispersion processing is carried out to data shown in table 1 using frequency spectrum weighted, coherent method, obtained Dispersion Characteristics data such as table
Shown in 24, corresponding frequency dispersion figure is as shown in figure 17, is compared it is found that using frequency spectrum with using the result (Figure 16) handled by the present invention
Weighted, coherent method can not obtain instrument dipole wave (in Figure 16, the dark colour region at middle part), can only obtain stratum dipole wave
(in Figure 16, the dark colour region of lower part), i.e., the corresponding related coefficient maximum value of dark colour in frequency dispersion figure.By frequency spectrum weighted, coherent
As shown in Table 25, the only data of stratum dipole wave for the frequency dispersion data that method obtains.It can be seen that being using frequency spectrum weighted, coherent method
It is unable to get accurate dispersion curve.
The 1st receiver waveform x of table 11(t) data
| Time t/ms | Wave data | Time t/ms | Wave data | Time t/ms | Wave data | Time t/ms | Wave data |
| 0 | -0.00058 | 1.232877 | 0.015351 | 2.465753 | 3.46E-05 | 3.70E+00 | -1.74E-05 |
| 0.058708 | 3.63E-05 | 1.29E+00 | 0.072232 | 2.524462 | 0.0002688 | 3.757339 | 8.77E-05 |
| 0.117417 | 0.000313 | 1.350294 | 0.075039 | 2.58317 | 0.0003215 | 3.816047 | 6.57E-05 |
| 0.176125 | -0.0003 | 1.409002 | -0.20121 | 2.641879 | 0.0002805 | 3.874755 | -3.42E-05 |
| 0.234834 | -0.00391 | 1.46771 | -0.3289 | 2.700587 | 0.000199 | 3.933464 | -7.73E-05 |
| 0.293542 | -0.00628 | 1.526419 | 0.21918 | 2.759295 | 8.15E-05 | 3.99E+00 | -1.45E-06 |
| 0.35225 | 0.02433 | 1.585127 | 0.4157 | 2.818004 | 6.29E-05 | 4.05E+00 | 8.50E-05 |
| 0.410959 | 0.040298 | 1.643836 | -0.01981 | 2.876712 | 2.04E-05 | 4.11E+00 | 8.09E-05 |
| 0.469667 | -0.05167 | 1.702544 | -0.20465 | 2.935421 | -0.000104 | 4.168297 | -6.27E-06 |
The 2nd receiver waveform x of table 22(t) data
The 3rd receiver waveform x of table 33(t) data
The 4th receiver waveform x of table 44(t) data
The 5th receiver waveform x of table 55(t) data
The 6th receiver waveform x of table 66(t) data
The 7th receiver waveform x of table 77(t) data
The 8th receiver waveform x of table 88(t) data
The data X of the 1st receiver time-frequency distributions of table 91(t,ω)
The data X of the 2nd receiver time-frequency distributions of table 102(t,ω)
The data X of the 3rd receiver time-frequency distributions of table 113(t,ω)
The data X of the 4th receiver time-frequency distributions of table 124(t,ω)
The data X of the 5th receiver time-frequency distributions of table 135(t,ω)
The data X of the 6th receiver time-frequency distributions of table 146(t,ω)
The data X of the 7th receiver time-frequency distributions of table 157(t,ω)
The data X of the 8th receiver time-frequency distributions of table 168(t,ω)
Related coefficient data volume when 17 frequencies omega of table is 1996.1Hz
Related coefficient data volume when 18 frequencies omega of table is 3992.2Hz
Related coefficient data volume when 19 frequencies omega of table is 5988.3Hz
Related coefficient data volume when 20 frequencies omega of table is 8982.4Hz
Related coefficient data volume when 21 frequencies omega of table is 10979Hz
22 Dispersion Characteristics data of table
Table 23 extracts obtained dispersion curve tables of data
24 Dispersion Characteristics data of table
Table 25
Embodiment one：
In acoustic logging while drilling, the knot of all frequency dispersion mode wave comprehensive effects can only be obtained using Weighted spectral coherent method
Fruit, and the Dispersion Characteristics of each mode wave can all be found out using the obtained handling result of method provided by the present invention
Come, the accurate frequency dispersion information of acoustic logging while drilling can be obtained, therefore this part shows two acoustic logging while drilling data processings
As a result.Figure 17 show the hard formation of theoretical modeling with brill dipole acoustic logging array waveform.Stratum and bore parameters such as table
Shown in 26.
Table 26
It can be obtained by theory analysis, at the common frequency range of acoustic logging while drilling (500-12000Hz frequency bands), with brill cement bond logging
Well dipole measurement pattern is primarily present two patterns, stratum dipole subpattern and stratum second order dipole subpattern.To Figure 18 institutes
Oscillography shape carries out dispersion analysis using the present invention, and it is as shown in figure 19 can to obtain its dispersion curve.The lower pattern of speed is that stratum is even
Pole wavelet, higher speed is stratum second order dipole wave.Curve shown in circle is theoretical stratum dipole wave frequency in figure
Non-dramatic song line, asterisk line are the stratum second order dipole wave curve of theory, as can be seen from Figure, two moulds of the solution of the present invention pair
The dispersion curve extraction of formula wave is accurate, can obtain reliably with brill hard formation dipole dispersion curve.As a comparison, it calculates
Using traditional Weighted spectral coherent method to waveform shown in Figure 18 carry out dispersion analysis as a result, as shown in figure 20, with Figure 19 mono-
Sample, curve shown in circle are theoretical stratum dipole wave dispersion curve, and asterisk line is theoretical stratum second order dipole wave
Curve, as can be seen from Figure, for the dipole wave on stratum, Weighted spectral coherent method can obtain effective dispersion curve knot
Fruit, due to there is multiple frequency dispersion mode waves (stratum second order dipole wave) in identical frequency range, this method is unable to get ground
The frequency dispersion information of layer second order dipole wave when stratum dipole wave signal is weaker, there is formation two only near 3100Hz
The signal of rank dipole wave responds.
Embodiment two：
Figure 21 is soft formation acoustic logging while drilling dipole waveform theoretical modeling curve.Stratum and bore parameters such as 27 institute of table
Show.
Table 27
When soft formation, under dipole measurement pattern, two patterns, stratum one are primarily present in 500-12000Hz frequency bands
The subpattern of rank dipole and instrument dipole subpattern.Since formation velocity is relatively low, instrument speed is relatively high, and when pattern has apparent
Instrument dipole wave be excited out.Figure 22 is (to be schemed with sonic data is bored to soft formation dipole using the solution of the present invention
21) dispersion analysis handling result.The lower pattern of speed is stratum dipole wave, and higher speed is instrument dipole wave.
Curve shown in circle is theoretical stratum dipole wave dispersion curve in figure, and asterisk line is the instrument dipole curve of theory
Line, as can be seen from Figure, the extraction of the dispersion curve of two mode waves of the present invention couple are accurate, can obtain reliably soft with boring
Layer dipole dispersion curve and instrument dipole dispersion curve.As a comparison, it calculates using traditional Weighted spectral coherent method
To waveform shown in Figure 21 carry out dispersion analysis as a result, as shown in figure 23, as can be seen from Figure, for the dipole on stratum
Wave, Weighted spectral coherent method can obtain effective dispersion curve as a result, due to there is multiple frequency dispersions in identical frequency range
Mode wave (instrument dipole wave), the frequency dispersion information that this method is unable to get instrument dipole wave are less than in addition, in low-frequency range
When 2300Hz, dispersion curve is influenced to have apparent shake by instrument dipole wave, is misfitted with theoretical curve.
It, can be with if the function is realized in the form of SFU software functional unit and when sold or used as an independent product
It is stored in a computer read/write memory medium.Based on this understanding, technical scheme of the present invention is substantially in other words
The part of the part that contributes to existing technology or the technical solution can be expressed in the form of software products, the meter
Calculation machine software product is stored in a storage medium, including some instructions are used so that a computer equipment (can be
People's computer, server or network equipment etc.) it performs all or part of the steps of the method described in the various embodiments of the present invention.
And storage medium above-mentioned includes：USB flash disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), arbitrary access are deposited
The various media that can store program code such as reservoir (RAM, Random Access Memory), magnetic disc or CD.
The above description is merely a specific embodiment, but scope of protection of the present invention is not limited thereto, any
Those familiar with the art in the technical scope disclosed by the present invention, can easily think of the change or the replacement, and should all contain
Lid is within protection scope of the present invention.Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims.
Claims (9)
1. dispersion curve determines method, which is characterized in that including：
Multiple initial data are obtained, each initial data is by the sonic data received by different receivers；
Short Time Fourier Transform is carried out to each initial data respectively, to obtain multiple time-frequency numeric field datas；
Velocity Time related operation is carried out to the multiple time-frequency numeric field data, to obtain characterization frequency-Velocity-time value and first
First related data body of the relationship of related coefficient；
According in the first related data body, the second related coefficient corresponding to each frequency-velocity amplitude determines the sonic data
The dispersion curve of corresponding sound wave, wherein the second related coefficient is the first related coefficient corresponding to each frequency-velocity amplitude
Maximum value；
It is described to the multiple time-frequency numeric field data carry out Velocity Time related operation, with obtain characterization frequency-Velocity-time value with
The related data body of the relationship of first related coefficient includes：
The numerical value with the first related coefficient corresponding to each frequency-Velocity-time value is calculated according to following formula,
Wherein, N is the total number of receiver in receiver array, and the position of window when T is, t is the time, and n=1,2,3 ... N are to connect
Device serial number is received, v is speed, and d is receiver spacing, TwFor when window width, ω is frequency；ρ (ω, v, T) is the first related coefficient
Numerical value, XnFor time-frequency numeric field data.
2. according to the method described in claim 1, it is characterized in that, step carries out Fourier in short-term to each initial data respectively
It converts, includes to obtain multiple time-frequency domain data：
Short Time Fourier Transform is carried out to obtaining the initial data according to following formula,
Wherein, xn(m) it is the received sonic data of n-th of receiver, n=1,2,3 ... number for receiver；M is representative
The variable of time, ω are frequency, and w (x) is preset Hanning window function；Xn(t, ω) is time-frequency numeric field data.
3. according to the method described in claim 2, it is characterized in that, the expression formula of the preset Hanning window function isWherein, NwFor Hanning window length.
4. according to the method described in claim 3, it is characterized in that, the value of the Hanning window length is 16.
5. according to the method described in claim 1, it is characterized in that, step according in the first related data body, each frequency-speed
The second related coefficient corresponding to angle value determines that the dispersion curve of sound wave corresponding to the sonic data includes：
Determine the second related data body of the relationship of characterization frequency-velocity amplitude and the second related coefficient；
To in the second related data body, frequency and velocity amplitude corresponding to the second related coefficient maximum value are sought, to obtain
The dispersion curve.
6. according to the method described in claim 5, it is characterized in that, step determines characterization frequency-velocity amplitude and the second phase relation
Second related data body of several relationships includes：
The second related coefficient is calculated according to following formula,
Wherein, ρ0For the numerical value of the second related coefficient, the position of window when T is, ρ (ω, v, T) is the numerical value of the first related coefficient.
7. dispersion curve determining device, which is characterized in that including：
Acquisition module, for obtaining multiple initial data, each initial data is by the sound received by different receivers
Wave number evidence；
Conversion module, for carrying out Short Time Fourier Transform to each initial data respectively, to obtain multiple time-frequency numeric field datas；
Related operation module, for carrying out Velocity Time related operation to the multiple time-frequency numeric field data, to obtain characterization frequency-
First related data body of Velocity-time value and the relationship of the first related coefficient；
First determining module, for according in the first related data body, the second phase relation corresponding to each frequency-velocity amplitude
Number, determines the dispersion curve of sound wave corresponding to the sonic data, wherein the second related coefficient is each frequency-velocity amplitude institute
The maximum value of corresponding first related coefficient；
Related operation module, specifically for according to following formula calculate corresponding to each frequency-Velocity-time value with the first phase
The numerical value of relationship number,
Wherein, N is the total number of receiver in receiver array, and the position of window when T is, t is the time, and n=1,2,3 ... N are to connect
Device serial number is received, v is speed, and d is receiver spacing, TwFor when window width, ω is frequency；ρ (ω, v, T) is the first related coefficient
Numerical value, XnFor time-frequency numeric field data.
Include that the dispersion curve of any one of claim 1-6 such as determines that method, feature exist 8. p-and s-wave velocity determines method
In further including：
It is worth the shear wave velocity that corresponding speed determines stratum according to the low-frequency cutoff of dispersion curve.
9. p-and s-wave velocity determining device, including dispersion curve determining device as claimed in claim 7, which is characterized in that further include：
Second determining module, for being worth the shear wave velocity that corresponding speed determines stratum according to the low-frequency cutoff of dispersion curve.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710265841.9A CN106934183B (en) | 2017-04-21 | 2017-04-21 | Dispersion curve determines that method and apparatus and p-and s-wave velocity determine method and apparatus |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710265841.9A CN106934183B (en) | 2017-04-21 | 2017-04-21 | Dispersion curve determines that method and apparatus and p-and s-wave velocity determine method and apparatus |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106934183A CN106934183A (en) | 2017-07-07 |
| CN106934183B true CN106934183B (en) | 2018-10-23 |
Family
ID=59437745
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710265841.9A Expired - Fee Related CN106934183B (en) | 2017-04-21 | 2017-04-21 | Dispersion curve determines that method and apparatus and p-and s-wave velocity determine method and apparatus |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106934183B (en) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP7030502B2 (en) * | 2017-12-20 | 2022-03-07 | 株式会社安藤・間 | Rock evaluation method |
| CN108345036A (en) * | 2018-01-10 | 2018-07-31 | 长江大学 | A kind of method and system of measurement while drilling crustal stress |
| CN110426740B (en) * | 2019-08-02 | 2022-02-15 | 中铁第四勘察设计院集团有限公司 | Seismic noise imaging exploration method and device and storage medium |
| CN112360447A (en) * | 2020-11-20 | 2021-02-12 | 中国石油天然气集团有限公司 | Method for evaluating reservoir perforation effect |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8392133B1 (en) * | 2010-06-28 | 2013-03-05 | The United States Of America As Represented By The Secretary Of The Navy | Method for measuring shear wavespeed in an isotropic plate |
| CN104678435A (en) * | 2014-10-27 | 2015-06-03 | 李欣欣 | Method for extracting Rayleigh surface wave frequency dispersion curve |
| CN104462667A (en) * | 2014-11-21 | 2015-03-25 | 电子科技大学 | Method for obtaining ultrasonic frequency dispersion curve |
-
2017
- 2017-04-21 CN CN201710265841.9A patent/CN106934183B/en not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| CN106934183A (en) | 2017-07-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106934183B (en) | Dispersion curve determines that method and apparatus and p-and s-wave velocity determine method and apparatus | |
| CN113759425B (en) | Method and system for evaluating filling characteristics of deep paleo-karst reservoir stratum by well-seismic combination | |
| CN108387933B (en) | A kind of method, apparatus and system of definitely interval quality factors | |
| US10145227B2 (en) | Method for estimating permeability of fractured rock formations from induced slow fluid pressure waves | |
| CN109490965B (en) | Method and device for quantitatively evaluating formation heterogeneity | |
| CN105116442B (en) | The reconstructing method of the weak seismic reflection signals of lithologic deposit | |
| CN106855636A (en) | Based on the prototype geological model Seismic forward method that carbonate reservoir is appeared | |
| CN104898161B (en) | Effective sandstone predicting method based on logging response simulator | |
| Spica et al. | Shallow VS imaging of the Groningen area from joint inversion of multimode surface waves and H/V spectral ratios | |
| CN113341455B (en) | Viscous anisotropic medium seismic wave numerical simulation method, device and equipment | |
| CN109991664A (en) | Seismic exploration in desert random noise method for reducing based on noise modeling analysis | |
| US11360232B2 (en) | Mapping wave slowness using multi-mode semblance processing techniques | |
| CN111123354A (en) | Method and equipment for predicting dense gas layer based on frequency-dependent reflection amplitude attenuation | |
| EP3281045A1 (en) | Method and system for the multimodal and multiscale analysis of geophysical data by transformation into musical attributes | |
| CN108594301B (en) | A kind of method and processing terminal of the seismic data fusion with difference characteristic | |
| CN105093300B (en) | A kind of boundary recognition of geological body method and device | |
| CN110320575A (en) | Method and device is determined based on the shale content of organic matter of petrophysical model | |
| CN106125133B (en) | It is a kind of based on gas cloud area constrain under fine velocity modeling method | |
| CN106324663B (en) | A kind of acquisition methods of quality factor | |
| CN110244383B (en) | Geological lithology comprehensive model establishing method based on near-surface data | |
| CN112630840B (en) | Random inversion method based on statistical characteristic parameters and processor | |
| CN109917459A (en) | A kind of method, apparatus and system for suppressing seismic noise | |
| CN114814951A (en) | Method, device, equipment and storage medium for acquiring slowness of fast and slow transverse waves | |
| CN110568491B (en) | Quality factor Q estimation method | |
| CN109885906B (en) | Magnetic resonance sounding signal sparse noise elimination method based on particle swarm optimization |
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 | ||
| TA01 | Transfer of patent application right | ||
| TA01 | Transfer of patent application right |
Effective date of registration: 20180904 Address after: 102249 18 Fu Xue Road, Changping District, Beijing Applicant after: China University of Petroleum (Beijing) Address before: 102249 18 Fu Xue Road, Changping District, Beijing Applicant before: Wang Bing |
|
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20181023 Termination date: 20210421 |