CN110187388A - A kind of stabilization earthquake quality factor q estimation method based on variation mode decomposition - Google Patents
A kind of stabilization earthquake quality factor q estimation method based on variation mode decomposition Download PDFInfo
- Publication number
- CN110187388A CN110187388A CN201910494524.3A CN201910494524A CN110187388A CN 110187388 A CN110187388 A CN 110187388A CN 201910494524 A CN201910494524 A CN 201910494524A CN 110187388 A CN110187388 A CN 110187388A
- Authority
- CN
- China
- Prior art keywords
- frequency
- function
- quality factor
- intrinsic mode
- seismic
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 59
- 238000000354 decomposition reaction Methods 0.000 title claims abstract description 43
- 230000006641 stabilisation Effects 0.000 title claims abstract description 14
- 238000011105 stabilization Methods 0.000 title claims abstract description 14
- 238000004458 analytical method Methods 0.000 claims abstract description 25
- 230000003595 spectral effect Effects 0.000 claims abstract description 15
- 238000005314 correlation function Methods 0.000 claims abstract description 9
- 238000001228 spectrum Methods 0.000 claims description 59
- 230000009466 transformation Effects 0.000 claims description 7
- 230000008569 process Effects 0.000 claims description 6
- 238000009499 grossing Methods 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 3
- 238000007405 data analysis Methods 0.000 claims 1
- 230000001419 dependent effect Effects 0.000 claims 1
- 238000003672 processing method Methods 0.000 abstract description 2
- 238000010521 absorption reaction Methods 0.000 description 3
- 230000008901 benefit Effects 0.000 description 3
- 238000009434 installation Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 239000004215 Carbon black (E152) Substances 0.000 description 1
- 101100340434 Saccharomyces cerevisiae (strain ATCC 204508 / S288c) IFM1 gene Proteins 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 229930195733 hydrocarbon Natural products 0.000 description 1
- 150000002430 hydrocarbons Chemical class 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000035699 permeability Effects 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
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 stabilization earthquake quality factor q estimation method based on variation mode decomposition that the present invention provides a kind of, belong to oil-gas exploration geophysics processing method field, by using variation mode decomposition, seismic data is decomposed into a series of intrinsic mode functions first, then preferably go out the main contributions intrinsic mode function component of the seismic signal using correlation function algorithm, preferably go out high frequency intrinsic mode function in conjunction with frequency analysis, utilize the Time-Frequency Analysis Method based on variation mode decomposition that 5 slip methods is combined preferably to go out the log-magnitude spectral coverage for estimating the decaying of seismic data adjacent layer bit data point the high frequency intrinsic mode function again, the good ground quality factor q estimated value of the seismic channel is provided in conjunction with least square method.Present invention reduces time-frequency combination variable estimation quality factors to need propagation time long disadvantage, avoids influencing each other for different frequency range, while having evaded the influence of noise, improves the estimated accuracy and stability of quality factor q.
Description
Technical field
It is specifically a kind of based on the steady of variation mode decomposition the present invention relates to oil-gas exploration geophysics processing method field
Determine earthquake quality factor q estimation method.
Background technique
The attenuation by absorption on stratum will lead to the reflected energy that ground receiver arrives and reduce, and dispersion occurs for speed.Quality factor q
This intrinsic medium character commonly used to the decaying for describing that seismic wavelet distortion and amplitude is caused to weaken, the inside of it and stratum
Structure feature, it is closely related containing fluid properties such as permeability, saturation degree, porosity etc., in the side such as layer description and hydrocarbon indication
Face has great importance.
Conventional time-domain earthquake quality factor q estimation method, such as amplitude damped method, analytic signal method, Rise time
Etc., true amplitude information is generally required, and true amplitude information is due to by wavefront expansion, transmission consumption and in addition to by inhaling
The influence of other factors except the decaying of frequency dependence caused by receiving and scattering is difficult to obtain from real data.Frequency domain
Quality factor q estimation method is shaken, such as it is suitable to generally require selection than method, peak (frequency) deviation method, centroid frequency deflection method for spectrum
Frequency range, and to noise-sensitive, in addition, frequency domain estimation Q also depends on the acquisition pattern of frequency spectrum.Utilize one-dimensional Fourier transform
In conjunction with window adding technology spectrum analysis often due to fixed window function length and cause Q estimated accuracy to decline.Utilize two-dimentional time-frequency
The frequency spectrum that spectrum obtains seismic signal has more advantage compared with one-dimensional spectrum.With it is variable when window two-dimentional time-frequency spectrum such as wavelet transformation, S become
It changes etc. in Q estimation to show and relatively consolidates timing window two dimension time-frequency spectrum such as Short Time Fourier Transform, Gabor transformation etc. can be preferably
Estimate attenuation characteristic.But these methods do not consider different frequency range seismic signal and decay different situations, do not avoid difference
The case where frequency range seismic signal influences each other causes quality factor q estimated accuracy low, and stability is not high.
Summary of the invention
It is an object of the invention to overcome the above-mentioned interactional feelings of different frequency range seismic signal existing in the prior art
Condition provides a kind of stabilization earthquake quality factor q estimation method based on variation mode decomposition from time-frequency domain, and one kind being based on variation mould
The stabilization earthquake quality factor q estimation method that state is decomposed, which is characterized in that steps are as follows,
(1) comprehensive utilization well logging, synthetic seismogram and geological information accurate calibration seismic channel target zone, determine earthquake number
According to analysis time range；
(2) within the analysis time of seismic data, using variation mode decomposition by the seismic data in step (1) by
Road is decomposed into a series of intrinsic mode functions；
(3) the combination frequency analysis of correlation function analysis method is utilized, high frequency sheet is preferably gone out from a series of intrinsic mode functions
Levy mode function；
(4) five are combined using the Time-Frequency Analysis Method based on variation mode decomposition to the high frequency intrinsic mode function preferably gone out
Point slip method preferably goes out the log-magnitude spectral coverage for estimating the decaying of seismic data adjacent layer bit data point；
(5) it using the preferred log-magnitude spectral coverage of least square method fitting high frequency intrinsic mode function, obtains in step (1)
The good ground quality factor q estimated value of the seismic channel target zone.
The technical solution of the application, according to comprehensive utilization well logging, synthetic seismogram and geological information accurate calibration earthquake
Road target zone determines the analysis time range of seismic data, within the analysis time of seismic data, by using variation mould
State is decomposed, and seismic data is decomposed into a series of intrinsic mode functions first, then preferably goes out the earthquake using correlation function algorithm
The main contributions intrinsic mode function component of signal preferably goes out high frequency intrinsic mode function in conjunction with frequency analysis, then to the high frequency
Intrinsic mode function utilizes the Time-Frequency Analysis Method based on variation mode decomposition that 5 slip methods is combined preferably to go out for estimating
The log-magnitude spectral coverage for counting the decaying of seismic data adjacent layer bit data point, the stabilization of the seismic channel is provided in conjunction with least square method
Stratum quality factor q estimated value.This method is that one kind combines decaying reason on the basis of variation mode decomposition adaptive decomposition characteristic
By the quality factor extracting method of development.
Using variation mode decomposition generate preferred intrinsic mode function carry out decay behavior, in conjunction with least square method when
Frequency domain estimates quality factor q, the seismic wave high band attenuation by absorption characteristic fast compared with low-frequency range is taken full advantage of, by preferred
High frequency intrinsic mode function carries out decay behavior and then estimates Q, shortens time-frequency combination variable estimation quality factor and needs to propagate
The disadvantage of time length, avoids influencing each other for different frequency range, while having evaded the influence of noise, improves quality factor q
Estimated accuracy and stability.
Key problem is the decomposition number for selecting suitable variation mode decomposition, in variation mode decomposition, best mode letter
Number decomposes the plyability of the installation warrants mode function frequency spectrum of number or orthogonality carries out judgement acquisition, utilizes correlation function analysis
Method combination frequency analysis preferably goes out to reflect the intrinsic mode function of very fast dampening information, then in time-frequency domain from the intrinsic mode
Dampening information is accurately estimated in function, and the quality factor q of seismic data adjacent layer position is estimated in conjunction with least square method, reaches and mentions
The purpose of high Q estimated value precision and stability.
Preferably, decomposable process is specific as follows in step (2):
1. variation mode decomposition is the constraint variation problem being expressed by the following equation
It is limited by
Wherein, f is original earthquake data, ukIt is k-th of mode, ωkIt is centre frequency, and ukAround ωkTight branch is presented
Support, δ (t) are impulse function, and j is imaginary unit；The bandwidth of each mode is believed by the Hilbert complementation analysis that its base band shifts
Number the only square H1 norm with positive frequency determine；
Unconstrained problem is converted by the constraint variation problem in formula (1) 2. being introduced into secondary punishment and Lagrange's multiplier:
Wherein, α is the balance parameters of data fidelity constraint；
3. using multiplier alternating direction algorithm seek formula (2) without constraint variation problem, and in each screening process
Different decomposition mode and centre frequency are generated, is expressed as from spectrum domain by each mode that solution obtains:
4. for each mode function u (t) of acquisition, earthquake instantaneous attribute by IMF component u (t) and IMF Martin Hilb
Spy's transformation y (t) is calculated by calculating following analytic signal:
Z (t)=u (t)+iy (t)=A (t) exp [i φ (t)], (4)
Wherein,P indicates Cauchy's principal value, and A (t) is instantaneous amplitude, and φ (t) is instantaneous phase；
Instantaneous amplitude A (t) and instantaneous phase φ (t) and instantaneous frequency ω (t) can be calculated by following equations:
It enables
H (ω, t)=Re { A (t) ei∫ω(t)dt, (6)
Wherein, Re indicates to take the real part of result；
A three-dimensional space [t, ω (t), A (t)], then, three-dimensional space are defined using time, instantaneous amplitude and frequency
Function H (ω, t) by being transformed into function [t, ω, H (ω, t)] Lai Shixian of three variables by [t, ω (t), A (t)], wherein
A (t)=H [ω (t), t)], thus, obtain the joint time-frequency distribution of each mode function of the seismic channel.
Preferably, step (3) specifically:
I seismic data indicates after variation mode decomposition are as follows:
Wherein, n is the number of intrinsic mode function；
If the centre frequency of each intrinsic mode function obtained by Fourier transformation is [ω1,ω2,L,ωn], pass through
The related coefficient of each intrinsic mode function and original earthquake data that correlation function obtains is [R1,R2,L,Rn], then it is logical first
It crosses related coefficient and preferably goes out the intrinsic mode function component that intrinsic mode function meets R >=0.3, then to the eigen mode preferably gone out
Wherein the maximum intrinsic mode function of centre frequency ω is estimated for subsequent quality factor for state function selection；
II in time-frequency domain, and the amplitude spectrum S (ω, τ) of seismic wave when traveling to current time τ from τ=0 is expressed as following formula:
Wherein, S0Amplitude spectrum when (ω) is τ=0, Q travel to being averaged between current time τ from τ=0 for seismic wave
Quality factor；
χ=ω τ is enabled, then the amplitude along variable χ follows following decaying exponential function relationship
Wherein, S (χa) it is in χ=χaWhen amplitude spectrum,
Logarithm is taken to formula (9) both sides, is arranged:
Wherein,
In the intrinsic mode function logarithmic spectrum preferably gone outIn be at the time of take when logarithmic spectrum amplitude maximum
χa, one section of frequency spectrum is taken since the logarithmic spectrum amplitude maximum moment, smoothly, preferably go out using 5 points of average gliding smoothing methods
For estimating the log-magnitude spectral coverage of the decaying of seismic data adjacent layer bit data point, the length selection of the wavelength coverage is depended on
The signal-to-noise ratio of actual seismic data.
Preferably, it is fitted preferred log-magnitude spectral coverage using least square method in step (5), the slope of acquisition is k, then institute
State the stabilization average quality factor that seismic channel is traveled to from τ=0 between current time τ are as follows:
In the technical solution of the application:
Compared with the prior art, the beneficial effects of the present invention are:
(1) decay behavior is carried out using the preferred intrinsic mode function that variation mode decomposition generates, in conjunction with least square method
Quality factor q is estimated in time-frequency domain, the seismic wave high band attenuation by absorption characteristic fast compared with low-frequency range is taken full advantage of, by excellent
The high frequency intrinsic mode function of choosing carries out decay behavior and then estimates Q, shortens time-frequency combination variable estimation quality factor needs
The disadvantage of propagation time length, avoids influencing each other for different frequency range, while having evaded the influence of noise, improves quality factor
The estimated accuracy and stability of Q；
(2) the decomposition number of suitable variation mode decomposition is selected, in variation mode decomposition, best mode function is decomposed
The plyability or orthogonality of several installation warrants mode function frequency spectrums carry out judgement acquisition, combine frequency using correlation function analysis method
Rate analysis preferably goes out to reflect the intrinsic mode function of very fast dampening information, then quasi- from the intrinsic mode function in time-frequency domain
Really estimation dampening information, the quality factor q of seismic data adjacent layer position is estimated in conjunction with least square method, is reached and is improved Q estimation
It is worth the purpose of precision and stability；
(3) the variation mode decomposition strong with high time frequency resolution and energy accumulating, noise robustness has been used, more
It is suitble to the processing of nonlinear and nonstationary seismic signal, ensure that the more accurate of calculated result；
(4) the algorithm speed of service is fast, is suitble to high-volume seismic data processing.
Detailed description of the invention
Fig. 1 is flow chart of the invention；
Fig. 2 is the synthetic seismic record sought for ground interval quality factors and its sheet generated after variation mode decomposition
Levy mode function；
Fig. 3 is the frequency spectrum of each intrinsic mode function；
The amplitude spectrum changed with variable χ and preferred amplitude section that Fig. 4 is preferred second intrinsic mode function IMF2
And the equivalent quality factor q estimated value of slope and stratum using the preferred amplitude section of logarithm that least square fitting method provides；
Fig. 5 is to estimate Q to the earthquake record than method using conventional spectrum；
Fig. 6 is the noisy synthetic seismic record sought for ground interval quality factors and corresponding after variation mode decomposition
Three intrinsic mode functions generated；
Fig. 7 is the frequency spectrum of each intrinsic mode function of noisy earthquake record；
Fig. 8 be the preferred second intrinsic mode function IMF2 of noisy earthquake record with variable χ variation amplitude spectrum and
The equivalent quality of slope and stratum of preferred amplitude spectral coverage and the preferred amplitude section of logarithm provided using least square fitting method because
Sub- Q estimated value；
Fig. 9 is to estimate Q to noisy earthquake record than method using conventional spectrum.
Specific embodiment
It is right combined with specific embodiments below in order to make those skilled in the art more fully understand technical solution of the present invention
The present invention is described in further detail.
Embodiment 1
Referring to Fig.1, a kind of stabilization earthquake quality factor q estimation method based on variation mode decomposition, steps are as follows,
(1) comprehensive utilization well logging, synthetic seismogram and geological information accurate calibration seismic channel target zone, determine earthquake number
According to analysis time range；
(2) within the analysis time of seismic data, using variation mode decomposition by the seismic data in step (1) by
Road is decomposed into a series of intrinsic mode functions；
(3) the combination frequency analysis of correlation function analysis method is utilized, high frequency sheet is preferably gone out from a series of intrinsic mode functions
Levy mode function；
(4) five are combined using the Time-Frequency Analysis Method based on variation mode decomposition to the high frequency intrinsic mode function preferably gone out
Point slip method preferably goes out the log-magnitude spectral coverage for estimating the decaying of seismic data adjacent layer bit data point；
(5) it using the preferred log-magnitude spectral coverage of least square method fitting high frequency intrinsic mode function, obtains in step (1)
The good ground quality factor q estimated value of the seismic channel target zone；
Decomposable process is specific as follows in step (2):
1. variation mode decomposition is the constraint variation problem being expressed by the following equation
It is limited by
Wherein, f is original earthquake data, ukIt is k-th of mode, ωkIt is centre frequency, and ukAround ωkTight branch is presented
Support, δ (t) are impulse function, and j is imaginary unit；
Unconstrained problem is converted by the constraint variation problem in formula (1) 2. being introduced into secondary punishment and Lagrange's multiplier:
Wherein, α is the balance parameters of data fidelity constraint；
3. using multiplier alternating direction algorithm seek formula (2) without constraint variation problem, and in each screening process
Different decomposition mode and centre frequency are generated, is expressed as from spectrum domain by each mode that solution obtains:
4. for each mode function u (t) of acquisition, earthquake instantaneous attribute by IMF component u (t) and IMF Martin Hilb
Spy's transformation y (t) is calculated by calculating following analytic signal:
Z (t)=u (t)+iy (t)=A (t) exp [i φ (t)], (4)
Wherein,P indicates Cauchy's principal value, and A (t) is instantaneous amplitude, and φ (t) is instantaneous phase；
Instantaneous amplitude A (t) and instantaneous phase φ (t) and instantaneous frequency ω (t) can be calculated by following equations:
It enables
Wherein, Re indicates to take the real part of result；
A three-dimensional space [t, ω (t), A (t)], then, three-dimensional space are defined using time, instantaneous amplitude and frequency
Function H (ω, t) by being transformed into function [t, ω, H (ω, t)] Lai Shixian of three variables by [t, ω (t), A (t)], wherein
A (t)=H [ω (t), t)], thus, obtain the joint time-frequency distribution of each mode function of the seismic channel；
Step (3) specifically:
I seismic data indicates after variation mode decomposition are as follows:
Wherein, n is the number of intrinsic mode function；
If the centre frequency of each intrinsic mode function obtained by Fourier transformation is [ω1,ω2,L,ωn], pass through
The related coefficient of each intrinsic mode function and original earthquake data that correlation function obtains is [R1,R2,L,Rn], then it is logical first
It crosses related coefficient and preferably goes out the intrinsic mode function component that intrinsic mode function meets R >=0.3, then to the eigen mode preferably gone out
Wherein the maximum intrinsic mode function of centre frequency ω is estimated for subsequent quality factor for state function selection；II in time-frequency domain, from τ
The amplitude spectrum S (ω, τ) of=0 seismic wave when traveling to current time τ is expressed as following formula:
Wherein, S0Amplitude spectrum when (ω) is τ=0, Q travel to being averaged between current time τ from τ=0 for seismic wave
Quality factor；
χ=ω τ is enabled, then the amplitude along variable χ follows following decaying exponential function relationship
Wherein, S (χa) it is in χ=χaWhen amplitude spectrum,
Logarithm is taken to formula (9) both sides, is arranged:
Wherein,
In the intrinsic mode function logarithmic spectrum preferably gone outIn be at the time of take when logarithmic spectrum amplitude maximum
χa, one section of frequency spectrum is taken since the logarithmic spectrum amplitude maximum moment, smoothly, preferably go out using 5 points of average gliding smoothing methods
For estimating the log-magnitude spectral coverage of the decaying of seismic data adjacent layer bit data point, the length selection of the wavelength coverage is depended on
The signal-to-noise ratio of actual seismic data；
It is fitted preferred log-magnitude spectral coverage using least square method in step (5), the slope of acquisition is k, the then earthquake
Road travels to the stabilization average quality factor between current time τ from τ=0 are as follows:
Embodiment 2
Referring to Fig. 2, the seismic data sought for ground interval quality factors and its sheet generated after variation mode decomposition
Mode function is levied, first eigenfunction IFM1, second eigenfunction IFM2 and the third generated using variation mode decomposition
A eigenfunction IFM3, respectively as shown in (b), (c) and (d) in Fig. 2, (a) is that the earthquake generated using minimum phase wavelet is closed
At record, Noise, theoretical qualities factor Q are not 50 to the composite traces, sample frequency 500HZ.
Embodiment 3
The correlativity of each eigenfunction and original synthetic seismic record known by table 1, first intrinsic mode function
(IMF1) meet the condition that related coefficient is greater than 0.3 with second intrinsic mode function (IMF2).
The correlativity of table 1 each eigenfunction and original synthetic seismic record
Embodiment 4
Referring to Fig. 3, the frequency spectrum of each intrinsic mode function, in Fig. 3, a, b and c are respectively the frequency of the frequency spectrum of IMF1, MF2
The frequency spectrum of spectrum and IMF3, from the figure, it can be seen that the centre frequency of IMF1, MF2 and IMF3 are respectively 27.7361Hz,
49.2254Hz and 31.2344Hz, in conjunction with IMF1 the and IMF2 component that table 1 preferably goes out, since the centre frequency of IMF2 is greater than IMF1
Centre frequency, so finally be preferably used in estimation ground interval quality factors intrinsic mode function be IMF2.
Embodiment 5
Referring to Fig. 4, preferred second intrinsic mode function IMF2 with the amplitude distribution figure of variable χ=ft and from variable
χaThe smooth amplitude section a of 5 moving average methods is started with, and quasi- using least square method to the smooth preferred amplitude section
The result figure b of conjunction, the slope of the fitting a straight line of acquisition are Kf=-0.0627, used herein is frequency f rather than angular frequency
ω, so needing to convert lower coefficient according to the π of ω=2 f, formula 11 becomes after conversion when calculating using formula 11Q value so as to estimate is 50.06, and estimated value meets theoretical value.
Embodiment 6
Referring to Fig. 5, Q, a, which are the time-frequency spectrum of earthquake record, to be estimated to noisy earthquake record than method using conventional spectrum, b is logarithm
Spectrum and least square fitting straight line extract the frequency spectrum of 0.2s and 2.2s from time-frequency spectrum, using spectrum than method to preferred frequency band
It is fitted, the slope of the fitting a straight line of acquisition is 0.1246, and the Q estimated value being calculated is 50.44, is given with Fig. 4 this technology
Q estimated value out is compared, it can be seen that in the case where not Noise, conventional method and this technology can be accurate
Estimate Q, and the Q estimate error that this technology provides is smaller.
Embodiment 7
Referring to Fig. 6, the noisy synthetic seismic record sought for ground interval quality factors and corresponding variation mode of passing through are divided
The white Gaussian noise that signal-to-noise ratio is 15dB is added to the synthetic seismogram in Fig. 2 in three intrinsic mode functions generated after solution
Generate noisy earthquake record.
Embodiment 8
The related coefficient of each intrinsic mode function and noisy synthetic seismic record, it is known from Table 2 that first intrinsic mode letter
Number (IMF1) and second intrinsic mode function (IMF2) meet the condition that related coefficient is greater than 0.3.
The related coefficient of table 2 each intrinsic mode function and noisy synthetic seismic record
Embodiment 9
It is the frequency spectrum of each intrinsic mode function of noisy earthquake record referring to Fig. 7, Fig. 7, a is the frequency spectrum of IMF1；B is
The frequency spectrum of IMF2；C is the frequency spectrum of IMF3；As can see from Figure 7, the centre frequency of each intrinsic mode function is respectively
27.7361Hz, 49.4753Hz and 159.1704Hz, in conjunction with IMF1 the and IMF2 component preferably gone out in table 2, due in IMF2
Frequency of heart is greater than the centre frequency of IMF1, so the intrinsic mode function for being finally preferably used in estimation ground interval quality factors is
IMF2。
Embodiment 10
Referring to Fig. 8, amplitude distribution figure with variable χ=ft that Fig. 8 is preferred second intrinsic mode function IMF2 and
From variable χaThe smooth amplitude section a of 5 moving average methods is started with, and minimum two is utilized to the smooth preferred amplitude section
The result figure b of multiplication fitting, the slope obtained after being fitted to the smooth preferred amplitude section using least square method are Kf=-
0.0611, the Q value so as to estimate is 51.41, it is known that, the quality factor q that the technology of the present invention is extracted is low in signal-to-noise ratio
In the case of the estimation Q value that provides still very close to theoretical true value.
Embodiment 11
Referring to Fig. 9, Fig. 9 is to estimate Q to noisy earthquake record than method using conventional spectrum, the time-frequency spectrum of the noisy earthquake record of a,
B logarithmic spectrum and least square fitting straight line；The frequency spectrum that 0.2s and 2.2s is extracted from time-frequency spectrum, using spectrum than method to preferred frequency
Rate section is fitted, and the slope of the fitting a straight line of acquisition is 0.0927, and the Q estimated value being calculated is 67.75, with this skill of Fig. 8
The Q estimated value that art provides is compared, it can be seen that the quality factor q that the technology of the present invention is extracted in the case where signal-to-noise ratio is low to
Q estimated value out is stablized, and precision is high.
The specific embodiment of the application above described embodiment only expresses, the description thereof is more specific and detailed, but simultaneously
The limitation to the application protection scope therefore cannot be interpreted as.It should be pointed out that for those of ordinary skill in the art
For, under the premise of not departing from technical scheme design, various modifications and improvements can be made, these belong to this
The protection scope of application.
Claims (4)
1. a kind of stabilization earthquake quality factor q estimation method based on variation mode decomposition, which is characterized in that steps are as follows,
(1) comprehensive utilization well logging, synthetic seismogram and geological information accurate calibration seismic channel target zone, determine seismic data
Analysis time range；
(2) within the analysis time of seismic data, using variation mode decomposition by the seismic data in step (1) by road point
Solution is a series of intrinsic mode functions；
(3) the combination frequency analysis of correlation function analysis method is utilized, high frequency eigen mode is preferably gone out from a series of intrinsic mode functions
State function；
(4) combine put down at 5 points using the Time-Frequency Analysis Method based on variation mode decomposition the high frequency intrinsic mode function preferably gone out
Equal sliding scale preferably goes out the log-magnitude spectral coverage for estimating the decaying of seismic data adjacent layer bit data point；
(5) it using the preferred log-magnitude spectral coverage of least square method fitting high frequency intrinsic mode function, obtains described in step (1)
The good ground quality factor q estimated value of seismic channel target zone.
2. a kind of stabilization earthquake quality factor q estimation method based on variation mode decomposition according to claim 1, special
Sign is that decomposable process is specific as follows in step (2):
1. variation mode decomposition is the constraint variation problem being expressed by the following equation
It is limited by
Wherein, f is original earthquake data, ukIt is k-th of mode, ωkIt is centre frequency, and ukAround ωkCompact schemes, δ (t) is presented
For impulse function, j is imaginary unit；
Unconstrained problem is converted by the constraint variation problem in formula (1) 2. being introduced into secondary punishment and Lagrange's multiplier:
Wherein, α is the balance parameters of data fidelity constraint；
3. seeking generating without constraint variation problem, and in each screening process for formula (2) using multiplier alternating direction algorithm
Different decomposition mode and centre frequency are expressed as from spectrum domain by each mode that solution obtains:
4. earthquake instantaneous attribute is become by the Hilbert of IMF component u (t) and IMF for each mode function u (t) of acquisition
It changes y (t) and is calculated by calculating following analytic signal:
Z (t)=u (t)+iy (t)=A (t) exp [i φ (t)], (4)
Wherein,P indicates Cauchy's principal value, and A (t) is instantaneous amplitude, and φ (t) is instantaneous phase；Instantaneous vibration
Width A (t) and instantaneous phase φ (t) and instantaneous frequency ω (t) can be calculated by following equations:
It enables
H (ω, t)=Re { A (t) ei∫ω(Ddt, (6)
Wherein, Re indicates to take the real part of result；
A three-dimensional space [t, ω (t), A (t)], then, three-dimensional space [t, ω are defined using time, instantaneous amplitude and frequency
(t), A (t)] by the way that function H (ω, t) to be transformed into function [t, ω, H (ω, t)] Lai Shixian of three variables, wherein A (t)
=H [ω (t), t)], thus, obtain the joint time-frequency distribution of each mode function of the seismic channel.
3. a kind of stabilization earthquake quality factor q estimation method based on variation mode decomposition according to claim 1, special
Sign is, step (3) specifically:
I seismic data indicates after variation mode decomposition are as follows:
Wherein, n is the number of intrinsic mode function；
If the centre frequency of each intrinsic mode function obtained by Fourier transformation is [ω1, ω2, L, ωn], pass through correlation
The related coefficient of each intrinsic mode function and original earthquake data that function obtains is [R1, R2, L, Rn], then pass through phase first
Relationship number preferably goes out the intrinsic mode function component that intrinsic mode function meets R >=0.3, then to the intrinsic mode letter preferably gone out
Wherein the maximum intrinsic mode function of centre frequency ω is estimated for subsequent quality factor for number selection；
For II in time-frequency domain, the amplitude spectrum S (ω, τ) of seismic wave when traveling to current time τ from τ=0 is expressed as following formula:
Wherein, S0Amplitude spectrum when (ω) is τ=0, Q be seismic wave from τ=0 travel to the average quality between current time τ because
Son；
χ=ω τ is enabled, then the amplitude along variable χ follows following decaying exponential function relationship
Wherein, S (χa) it is in χ=χaWhen amplitude spectrum,
Logarithm is taken to formula (9) both sides, is arranged:
Wherein,
In the intrinsic mode function logarithmic spectrum preferably gone outIn when taking logarithmic spectrum amplitude maximum at the time of for χa, from
The logarithmic spectrum amplitude maximum moment starts to take one section of frequency spectrum, carries out smoothly, preferably going out to be used for using 5 points of average gliding smoothing methods
Estimate the log-magnitude spectral coverage of the decaying of seismic data adjacent layer bit data point, the length selection of the wavelength coverage is dependent on practical
The signal-to-noise ratio of seismic data.
4. a kind of stabilization earthquake quality factor q estimation method based on variation mode decomposition according to claim 1, special
Sign is, is fitted preferred log-magnitude spectral coverage using least square method in step (5), and the slope of acquisition is k, then the seismic channel
The stabilization average quality factor between current time τ is traveled to from τ=0 are as follows:
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910494524.3A CN110187388B (en) | 2019-06-06 | 2019-06-06 | Stable seismic quality factor Q estimation method based on variational modal decomposition |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910494524.3A CN110187388B (en) | 2019-06-06 | 2019-06-06 | Stable seismic quality factor Q estimation method based on variational modal decomposition |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110187388A true CN110187388A (en) | 2019-08-30 |
| CN110187388B CN110187388B (en) | 2021-01-05 |
Family
ID=67720836
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910494524.3A Expired - Fee Related CN110187388B (en) | 2019-06-06 | 2019-06-06 | Stable seismic quality factor Q estimation method based on variational modal decomposition |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110187388B (en) |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110515127A (en) * | 2019-09-26 | 2019-11-29 | 中国石油大学(北京) | A kind of earthquake quality factor determines method, apparatus, equipment, medium |
| CN112711070A (en) * | 2019-10-24 | 2021-04-27 | 中国石油化工股份有限公司 | Oil-gas detection method and device based on seismic signal decomposition |
| CN112712047A (en) * | 2021-01-08 | 2021-04-27 | 自然资源部第一海洋研究所 | Marine mammal echo positioning signal detection method based on image processing |
| CN114089416A (en) * | 2021-11-17 | 2022-02-25 | 成都理工大学 | Method for estimating seismic wave attenuation gradient by utilizing Schrodinger equation |
| CN114114406A (en) * | 2020-08-26 | 2022-03-01 | 中国石油化工股份有限公司 | Reservoir permeability estimation method and device |
| CN117370737A (en) * | 2023-12-08 | 2024-01-09 | 成都信息工程大学 | Unsteady state non-Gaussian noise removing method based on self-adaptive Gaussian filter |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105044777A (en) * | 2015-07-01 | 2015-11-11 | 中国石油天然气股份有限公司 | Method for detecting strong reflection amplitude elimination of seismic marker layer based on empirical mode decomposition |
| CN105093294A (en) * | 2015-06-04 | 2015-11-25 | 成都信息工程大学 | Method for estimating attenuation gradient of seismic waves based on variable mode decomposition |
| CN108646289A (en) * | 2018-03-19 | 2018-10-12 | 中国海洋石油集团有限公司 | A method of estimation earthquake quality factor |
| CN108845357A (en) * | 2018-06-13 | 2018-11-20 | 成都信息工程大学 | A method of the equivalent quality factor in stratum is estimated based on the synchronous wavelet transformation that squeezes |
-
2019
- 2019-06-06 CN CN201910494524.3A patent/CN110187388B/en not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105093294A (en) * | 2015-06-04 | 2015-11-25 | 成都信息工程大学 | Method for estimating attenuation gradient of seismic waves based on variable mode decomposition |
| CN105044777A (en) * | 2015-07-01 | 2015-11-11 | 中国石油天然气股份有限公司 | Method for detecting strong reflection amplitude elimination of seismic marker layer based on empirical mode decomposition |
| CN108646289A (en) * | 2018-03-19 | 2018-10-12 | 中国海洋石油集团有限公司 | A method of estimation earthquake quality factor |
| CN108845357A (en) * | 2018-06-13 | 2018-11-20 | 成都信息工程大学 | A method of the equivalent quality factor in stratum is estimated based on the synchronous wavelet transformation that squeezes |
Non-Patent Citations (2)
| Title |
|---|
| KONSTANTIN DRAGOMIRETSKIY ET AL.: "Variational Mode Decompositiom", 《IEEE TRANSACTIONS ON SIGNAL PROCESSING》 * |
| 何元 等: "变分模态分解及其在地震去噪中的应用", 《中国地球科学联合学术年会2014》 * |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110515127A (en) * | 2019-09-26 | 2019-11-29 | 中国石油大学(北京) | A kind of earthquake quality factor determines method, apparatus, equipment, medium |
| CN112711070A (en) * | 2019-10-24 | 2021-04-27 | 中国石油化工股份有限公司 | Oil-gas detection method and device based on seismic signal decomposition |
| CN112711070B (en) * | 2019-10-24 | 2024-02-20 | 中国石油化工股份有限公司 | Oil gas detection method and device based on seismic signal decomposition |
| CN114114406A (en) * | 2020-08-26 | 2022-03-01 | 中国石油化工股份有限公司 | Reservoir permeability estimation method and device |
| CN114114406B (en) * | 2020-08-26 | 2023-09-01 | 中国石油化工股份有限公司 | Reservoir permeability estimation method and device |
| CN112712047A (en) * | 2021-01-08 | 2021-04-27 | 自然资源部第一海洋研究所 | Marine mammal echo positioning signal detection method based on image processing |
| CN112712047B (en) * | 2021-01-08 | 2022-09-16 | 自然资源部第一海洋研究所 | Marine mammal echo positioning signal detection method based on image processing |
| CN114089416A (en) * | 2021-11-17 | 2022-02-25 | 成都理工大学 | Method for estimating seismic wave attenuation gradient by utilizing Schrodinger equation |
| CN114089416B (en) * | 2021-11-17 | 2023-02-21 | 成都理工大学 | Method for estimating attenuation gradient of seismic waves by utilizing Schrodinger equation |
| CN117370737A (en) * | 2023-12-08 | 2024-01-09 | 成都信息工程大学 | Unsteady state non-Gaussian noise removing method based on self-adaptive Gaussian filter |
| CN117370737B (en) * | 2023-12-08 | 2024-02-06 | 成都信息工程大学 | Unsteady state non-Gaussian noise removing method based on self-adaptive Gaussian filter |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110187388B (en) | 2021-01-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110187388A (en) | A kind of stabilization earthquake quality factor q estimation method based on variation mode decomposition | |
| RU2579164C1 (en) | Handling method for determining quality of geologic environment | |
| US20080270033A1 (en) | Methods of hydrocarbon detection using spectral energy analysis | |
| CN102305941B (en) | Method for determining stratum stack quality factor by direct scanning of prestack time migration | |
| Wang | The W transform | |
| Niu et al. | Block sparse Bayesian learning for broadband mode extraction in shallow water from a vertical array | |
| CN107132577B (en) | A kind of seismic attenuation method of estimation based on area under spectrum variation | |
| CN104820242B (en) | A kind of road collection amplitude towards prestack inversion divides compensation method | |
| Zhang | Time‐phase amplitude spectra based on a modified short‐time Fourier transform | |
| CN108845357A (en) | A method of the equivalent quality factor in stratum is estimated based on the synchronous wavelet transformation that squeezes | |
| CN102692647A (en) | Stratum oil-gas possibility prediction method with high time resolution | |
| CN102169188A (en) | Method for surveying oil and gas based on Morlet spectrum | |
| Li et al. | Wigner-Ville distribution and its application in seismic attenuation estimation | |
| CN108594301B (en) | A kind of method and processing terminal of the seismic data fusion with difference characteristic | |
| CN108957540B (en) | Method for efficiently extracting attenuation quality factors in complex reservoir | |
| Yang et al. | Q-factor estimation using bisection algorithm with power spectrum | |
| CN112711070B (en) | Oil gas detection method and device based on seismic signal decomposition | |
| CN117452491A (en) | Combined exploration method for identifying characteristics of gas reservoirs of coal series under complicated mountain land surface conditions | |
| CN108957529B (en) | Attribute-based wellless wavelet estimation method | |
| Cheng et al. | Q estimation based on the logarithmic spectral area double difference | |
| CN107272060B (en) | The extracting method and system of a kind of ground interval quality factors | |
| CN110568491B (en) | Quality factor Q estimation method | |
| CN110515127B (en) | Method, device, equipment and medium for determining seismic quality factor | |
| CN112684501A (en) | Q value estimation method based on spectral ratio area and application | |
| CN105572742A (en) | Method and device for determining depth of seawater |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20210105 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |