CN110646841A - Time-varying sparse deconvolution method and system - Google Patents
Time-varying sparse deconvolution method and system Download PDFInfo
- Publication number
- CN110646841A CN110646841A CN201810679707.8A CN201810679707A CN110646841A CN 110646841 A CN110646841 A CN 110646841A CN 201810679707 A CN201810679707 A CN 201810679707A CN 110646841 A CN110646841 A CN 110646841A
- Authority
- CN
- China
- Prior art keywords
- seismic
- reflection coefficient
- wavelet
- optimized
- optimization model
- 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 52
- 238000005457 optimization Methods 0.000 claims abstract description 74
- 239000011159 matrix material Substances 0.000 claims abstract description 54
- 238000004590 computer program Methods 0.000 claims description 3
- 238000010521 absorption reaction Methods 0.000 abstract description 5
- 230000036039 immunity Effects 0.000 abstract description 5
- 230000002452 interceptive effect Effects 0.000 abstract description 5
- 238000001228 spectrum Methods 0.000 description 6
- 238000001914 filtration Methods 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- 230000002238 attenuated effect Effects 0.000 description 3
- 230000008901 benefit Effects 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000014759 maintenance of location Effects 0.000 description 2
- 230000008859 change Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000010363 phase shift Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
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
A time-varying sparse deconvolution method and system are disclosed. The method can comprise the following steps: step 1: inputting a seismic record, initial seismic wavelets and sparsity; step 2: estimating a Q value vector, and constructing an attenuation matrix; and step 3: according to the attenuation matrix, establishing a reflection coefficient optimization model and a seismic wavelet optimization model; and 4, step 4: obtaining an optimized reflection coefficient according to the initial seismic wavelet, the sparsity and the reflection coefficient optimization model; and 5: obtaining optimized seismic wavelets according to the optimized reflection coefficient and the seismic wavelet optimization model; step 6: and repeating the step 4-5 until the reflection coefficient difference is smaller than the reflection coefficient difference threshold or the seismic wavelet difference is smaller than the seismic wavelet difference threshold. The method carries out absorption attenuation estimation, carries out interactive iteration on reflection coefficient optimization and seismic source wavelet optimization to obtain a more accurate result, improves the resolution ratio of seismic data, well keeps the relative amplitude of the reflection coefficient and the phase deviation of the reflection coefficient, and has good fault tolerance and noise immunity.
Description
Technical Field
The invention relates to the technical field of oil-gas geophysical, in particular to a time-varying sparse deconvolution method and a time-varying sparse deconvolution system.
Background
In seismic exploration data processing, the resolution of seismic data is very important, people are exploring a link-attack target by taking high signal-to-noise ratio, high resolution and high fidelity, and the resolution of the seismic data directly influences the results of subsequent inversion imaging and the like. With the deep oil-gas exploration, people have higher and higher requirements on the resolution of seismic data, pay more attention to amplitude, frequency, phase and other important information, and meanwhile, higher requirements are provided for accurately extracting the position of a reflection coefficient, the extracted reflection coefficient has good relative amplitude retention and the like, and the real lithologic information and the thin-layer structure before and after processing are expected to be reflected.
Deconvolution is one of effective means for improving the resolution of seismic data, and conventional classical methods are steady-state deconvolution methods, i.e. assuming that the earth layers are completely elastic, seismic wavelets are time-invariant in the propagation process, but the method is far from the real situation. Due to the effects of earth filtering and the like, the seismic wavelets are not invariable in the propagation process, but have amplitude attenuation, frequency change and phase distortion all the time, so that some traditional classical methods are not very suitable for practical industry. The results obtained after conventional deconvolution are also not convincing, because the resulting deconvolution results destroy the relative amplitude retention of the reflection coefficient, while greatly affecting the position of the reflection coefficient. At present, many techniques estimate the time-varying wavelet in the frequency domain or estimate the amplitude spectrum and phase spectrum of the instantaneous wavelet by time-frequency analysis to obtain the time-varying wavelet, and then obtain the reflection coefficient. Therefore, there is a need to develop a time-varying sparse deconvolution method and system based on compressed sensing.
The information disclosed in this background section is only for enhancement of understanding of the general background of the invention and should not be taken as an acknowledgement or any form of suggestion that this information forms the prior art already known to a person skilled in the art.
Disclosure of Invention
The invention provides a time-varying sparse deconvolution method and a time-varying sparse deconvolution system, which can decompose an unsteady state into an attenuation estimation and a steady state, perform absorption attenuation estimation, and perform interactive iteration on reflection coefficient optimization and seismic source wavelet optimization, so as to obtain a relatively accurate result.
According to an aspect of the present invention, a time-varying sparse deconvolution method is presented. The method may include: step 1: inputting a seismic record, initial seismic wavelets and sparsity; step 2: estimating a Q value vector, and constructing an attenuation matrix; and step 3: establishing a reflection coefficient optimization model and a seismic wavelet optimization model according to the attenuation matrix; and 4, step 4: obtaining an optimized reflection coefficient according to the initial seismic wavelet, the sparsity and the reflection coefficient optimization model; and 5: obtaining an optimized seismic wavelet according to the optimized reflection coefficient and the seismic wavelet optimization model; step 6: and repeating the step 4-5 until the reflection coefficient difference is smaller than the reflection coefficient difference threshold or the seismic wavelet difference is smaller than the seismic wavelet difference threshold.
Preferably, the reflection coefficient optimization model is:
where r represents the reflection coefficient, y represents the seismic record, A (Q) is the attenuation matrix for the Q value, Q is the quality factor, and W is the seismic wavelet matrix.
Preferably, the seismic wavelet optimization model is:
wherein R (w) is seismic wavelet constraint, w is seismic wavelet, y represents seismic record, Q is quality factor,i is the matrix with the ith diagonal element being 1 and the remainder being 0.
Preferably, the step 4 obtains the optimized reflection coefficient through an OMP algorithm, and includes: and calculating an optimized seismic record according to the seismic record, the initial seismic wavelet and the sparsity, seeking an index with the largest atomic inner product in a matrix of the residual error and the initial seismic wavelet, updating an index set, and solving the optimized reflection coefficient of the reflection coefficient optimization model under the condition of least square.
Preferably, the step 5 comprises: and calculating the difference between the re-optimized seismic record and the actual seismic record according to the optimized reflection coefficient, and solving the optimized seismic wavelet by using a least square method on the premise that the 2 norm in the seismic wavelet optimization model is the minimum.
According to another aspect of the invention, a time-varying sparse deconvolution system is proposed, having stored thereon a computer program which when executed by a processor performs the steps of: step 1: inputting a seismic record, initial seismic wavelets and sparsity; step 2: estimating a Q value vector, and constructing an attenuation matrix; and step 3: establishing a reflection coefficient optimization model and a seismic wavelet optimization model according to the attenuation matrix; and 4, step 4: obtaining an optimized reflection coefficient according to the initial seismic wavelet, the sparsity and the reflection coefficient optimization model; and 5: obtaining an optimized seismic wavelet according to the optimized reflection coefficient and the seismic wavelet optimization model; step 6: and repeating the step 4-5 until the reflection coefficient difference is smaller than the reflection coefficient difference threshold or the seismic wavelet difference is smaller than the seismic wavelet difference threshold.
Preferably, the reflection coefficient optimization model is:
where r represents the reflection coefficient, y represents the seismic record, A (Q) is the attenuation matrix for the Q value, Q is the quality factor, and W is the seismic wavelet matrix.
Preferably, the seismic wavelet optimization model is:
wherein R (w) is seismic wavelet constraint, w is seismic wavelet, y represents seismic record, Q is quality factor,
i is the matrix with the ith diagonal element being 1 and the remainder being 0.
Preferably, the step 4 obtains the optimized reflection coefficient through an OMP algorithm, and includes: and calculating an optimized seismic record according to the seismic record, the initial seismic wavelet and the sparsity, seeking an index with the largest atomic inner product in a matrix of the residual error and the initial seismic wavelet, updating an index set, and solving the optimized reflection coefficient of the reflection coefficient optimization model under the condition of least square.
Preferably, the step 5 comprises: and calculating the difference between the re-optimized seismic record and the actual seismic record according to the optimized reflection coefficient, and solving the optimized seismic wavelet by using a least square method on the premise that the 2 norm in the seismic wavelet optimization model is the minimum.
The present invention has other features and advantages which will be apparent from or are set forth in detail in the accompanying drawings and the following detailed description, which are incorporated herein, and which together serve to explain certain principles of the invention.
Drawings
The above and other objects, features and advantages of the present invention will become more apparent by describing in more detail exemplary embodiments thereof with reference to the attached drawings, in which like reference numerals generally represent like parts.
FIG. 1 shows a flow chart of the steps of a time-varying sparse deconvolution method according to the present invention.
FIGS. 2a, 2b, 2c, 2d show schematic diagrams of a 30Hz Rake wavelet, reflection coefficients, a stationary unattenuated seismic recording, and an actual attenuated seismic recording, respectively, according to one embodiment of the invention.
FIGS. 3a and 3b show a comparison of a wavelet extracted according to the prior art and a wavelet extracted according to the present method and a real wavelet, respectively, according to an embodiment of the present invention.
Fig. 4a, 4b show graphs comparing deconvolution results obtained by the prior art and the present method with true reflection coefficients, respectively, according to an embodiment of the present invention.
FIGS. 5a, 5b show graphs comparing a real wavelet to an estimated wavelet of 110ms and 700ms, respectively, according to one embodiment of the present invention.
FIG. 6 shows a graph comparing a real wavelet to an estimated wavelet at Q49 according to one embodiment of the present invention.
Fig. 7 shows a graph comparing the true reflection coefficient to the estimated reflection coefficient at Q49 according to one embodiment of the present invention.
FIG. 8 shows a graph comparing a real wavelet to an estimated wavelet at a Q of 55 according to one embodiment of the present invention.
Fig. 9 shows a graph comparing the true reflection coefficient to the estimated reflection coefficient at Q49 according to one embodiment of the present invention.
FIG. 10 shows a comparison of a real wavelet to an estimated wavelet at 5% noise according to one embodiment of the present invention.
FIG. 11 shows a graph comparing the true reflection coefficient to the estimated reflection coefficient at 5% noise according to one embodiment of the invention.
Detailed Description
The invention will be described in more detail below with reference to the accompanying drawings. While the preferred embodiments of the present invention are shown in the drawings, it should be understood that the present invention may be embodied in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the invention to those skilled in the art.
FIG. 1 shows a flow chart of the steps of a time-varying sparse deconvolution method according to the present invention.
In this embodiment, the time-varying sparse deconvolution method according to the present invention may comprise: step 1: inputting a seismic record, initial seismic wavelets and sparsity; step 2: estimating a Q value vector, and constructing an attenuation matrix; and step 3: according to the attenuation matrix, establishing a reflection coefficient optimization model and a seismic wavelet optimization model; and 4, step 4: obtaining an optimized reflection coefficient according to the initial seismic wavelet, the sparsity and the reflection coefficient optimization model; and 5: obtaining optimized seismic wavelets according to the optimized reflection coefficient and the seismic wavelet optimization model; step 6: and repeating the step 4-5 until the reflection coefficient difference is smaller than the reflection coefficient difference threshold or the seismic wavelet difference is smaller than the seismic wavelet difference threshold.
In one example, the reflection coefficient optimization model is:
where r represents the reflection coefficient, y represents the seismic record, A (Q) is the attenuation matrix for the Q value, Q is the quality factor, and W is the seismic wavelet matrix.
In one example, the seismic wavelet optimization model is:
wherein R (w) is seismic wavelet constraint, w is seismic wavelet, y represents seismic record, Q is quality factor,i is the matrix with the ith diagonal element being 1 and the remainder being 0.
In one example, step 4 obtains the optimized reflection coefficient through an OMP algorithm, including: and calculating and optimizing the seismic record according to the seismic record, the initial seismic wavelet and the sparsity, seeking an index with the largest atomic inner product in a matrix of the residual error and the initial seismic wavelet, updating an index set, and solving the optimized reflection coefficient of the reflection coefficient optimization model under the condition of least square.
In one example, step 5 comprises: and calculating the difference between the re-optimized seismic record and the actual seismic record according to the optimized reflection coefficient, and solving the optimized seismic wavelet by using a least square method on the premise of the minimum 2 norm in the seismic wavelet optimization model.
Specifically, the classical convolution model is proposed by Robinson, that is, the seismic record is obtained by convolution of seismic wavelets and reflection coefficients, where y (t) is the seismic record, w (t) is the seismic wavelets, and r (t) is the reflection coefficients, and thus formula (3):
y(t)＝w(t)*r(t)+n(t) (3)
where n (t) is seismic noise and denotes convolution operators. In conventional processing, many methods use this model, assuming that the seismic wavelet is time-invariant, but in the actual propagation process of the seismic wavelet, the seismic wavelet is time-variant due to earth filtering and the like, which causes the certainty of its high-frequency components and phase distortion. Therefore, Margrave proposes a non-steady-state convolution model, and introduces Q theory into the steady-state convolution model, and the attenuation function is formula (4):
|α(t,f)|＝e-π|f|t/Q (4)
where t is time, f is frequency, and Q is quality factor. In fact, the amplitude spectrum and the phase spectrum of the decay function depend only on the Q value, which is in the form of equation (5) in the time domain:
a(t,τ)＝∫α(t,f)e2πifτdf (5)
it is assumed here that the amplitude spectrum and the phase spectrum of the decay function satisfy the assumption of minimum phase. Since the seismic wavelets vary during propagation, the entire convolution model becomes more appropriate for equation (6):
y(t)＝w(t,Q)*r(t)+n(t) (6)
considering the splitting of the seismic wavelet, the seismic wavelet at each time is obtained by attenuating the source wavelet, and therefore, the model is transformed into formula (7):
y(t)＝w(t)*a(t,τ)·r(t)+n(t) (7)
wherein the operation of ". about." is still convolution, which can be understood as filtering of the source wavelet and the attenuation function at each moment, the use of ". about." is only beneficial for understanding, and is not an actual operation, because the seismic wavelet is attenuated during the convolution process with the reflection coefficient, and therefore, the convolution operation is not in the conventional sense, so that we convert it into a matrix form and solve it more easily, namely, formula (8):
y＝WA(Q)r (8)
where y is the seismic record, W is the Toeplitz matrix formed by the source wavelets, A (Q) is the attenuating filter matrix for the Q value, and each column of A (Q) is effectively the inverse Fourier transform of α (t, f).
Therefore, the overall model is formula (9):
the Q value in each trace set is estimated by using the theory in the existing inverse Q filtering, and then the solution of the model is completed, so that the model (9) has small error disturbance on the Q value and can still obtain an accurate result.
Due to sparse deconvolution, the reflection coefficients in the model are constrained by L0, and R (w) is constrained by seismic wavelets.
The model (9) is a blind problem and is not beneficial to direct solution, the model is divided into reflection coefficients and seismic source wavelets for interactive optimization solution, in the theory of compressed sensing, original data are recovered from signal observation, W can be used as an observation matrix for reflection coefficient optimization to solve, and because W is not full rank, although not in accordance with RIP conditions, the generalized RIP conditions are met, and the solution can still be realized by using a compressed sensing method.
Thus, a time-varying sparse deconvolution method according to the present invention may comprise:
step 1: and inputting the seismic record, the initial seismic wavelet and the sparsity.
Step 2: and estimating a Q value vector, and constructing an attenuation matrix of A (Q).
And step 3: according to the attenuation matrix, a reflection coefficient optimization model is established as formula (1), and due to the assumption of sparse reflection coefficients, in the theory of compressed sensing, the L0 norm and the L1 norm can be equivalent, so that solving the formula (1) can use | | | r | | survival1Since W is a Toeplitz matrix, the compressive sensing method for solving the observation matrix into the Toeplitz matrix is generally a fast iteration threshold (FISTA, L1) and an orthogonal matching pursuit algorithm (OMP, L0), wherein the FISTA method needs to use a backsracking format, and the conventional format cannot be used for solving. Whether FISTA&backsracking or OMP methods, all of which require known sparsity, and the empirical relationship between the number of reflection coefficients and the seismic record peaks under sparse deconvolution is as follows:
K＝cK0 (10)
wherein K is the number of reflection coefficients, K0And c is the empirical range of 1.5-3 for the number of seismic recording wave crests, and sparsity can be obtained according to the number of reflection coefficients.
And (3) after obtaining the optimized reflection coefficient according to the formula (1), optimizing according to the optimized reflection coefficient to obtain the optimized seismic wavelet.
Since W is a Toeplitz matrix, let the source wavelet be W ═ ω (ω)-(l-1),…,ω0,…,ωl-1) The thus constructed Toeplitz matrix shape is:
wherein L < < n, after the whole wavelet is expanded to 2n-1 length by 0, the whole source wavelet shows sparse characteristic, and meanwhile, Ricker wavelet in theory is very smooth, so TV norm and L2 norm are introduced, and the constraint of the whole source wavelet is as follows:
R(w)＝β||w||1+β1∑(ωi-ωi-1)+β2||w||2 (12)。
since the Toeplitz matrix has the same diagonal elements in each diagonal, it can be split, so
The initial model for optimizing the source wavelet is therefore:
and (3) carrying out deformation simplification on the model to obtain a seismic source wavelet optimization model as a formula (2).
And 4, step 4: obtaining an optimized reflection coefficient through an OMP algorithm according to the initial seismic wavelet, the sparsity and the reflection coefficient optimization model: and calculating and optimizing the seismic record according to the seismic record, the initial seismic wavelet and the sparsity, seeking an index with the largest atomic inner product in a matrix of the residual error and the initial seismic wavelet, updating an index set, and solving the optimized reflection coefficient of the reflection coefficient optimization model under the condition of least square.
And 5: and calculating the difference between the re-optimized seismic record and the actual seismic record according to the optimized reflection coefficient, and solving the optimized seismic wavelet by using a least square method on the premise of the minimum 2 norm in the seismic wavelet optimization model.
Step 6: and repeating the step 4-5 until the reflection coefficient difference is smaller than the reflection coefficient difference threshold or the seismic wavelet difference is smaller than the seismic wavelet difference threshold.
The invention decomposes the unsteady state into attenuation estimation and a steady state condition, carries out absorption attenuation estimation, and carries out interactive iteration on reflection coefficient optimization and seismic source wavelet optimization, thereby obtaining more accurate results, improving seismic data resolution, better keeping the relative amplitude of the reflection coefficient and the phase deviation of the reflection coefficient, and having good fault tolerance and noise immunity.
Application example
To facilitate understanding of the solution of the embodiments of the present invention and the effects thereof, a specific application example is given below. It will be understood by those skilled in the art that this example is merely for the purpose of facilitating an understanding of the present invention and that any specific details thereof are not intended to limit the invention in any way.
FIGS. 2a, 2b, 2c, 2d show schematic diagrams of a 30Hz Rake wavelet, reflection coefficients, a stationary unattenuated seismic recording, and an actual attenuated seismic recording, respectively, according to one embodiment of the invention.
The time-varying sparse deconvolution method according to the present invention comprises:
step 1: inputting seismic records, initial seismic wavelets and sparsity, as shown in FIGS. 2a-2 d;
step 2: estimating a Q value vector, and constructing an attenuation matrix;
and step 3: according to the attenuation matrix, establishing a reflection coefficient optimization model as a formula (1) and a seismic wavelet optimization model as a formula (2);
and 4, step 4: obtaining an optimized reflection coefficient through an OMP algorithm according to the initial seismic wavelet, sparsity and reflection coefficient optimization model, wherein the method comprises the following steps: calculating and optimizing the seismic record according to the seismic record, the initial seismic wavelet and the sparsity, seeking an index with the largest atomic inner product in a matrix of the residual error and the initial seismic wavelet, updating an index set, and solving an optimized reflection coefficient of a reflection coefficient optimization model under the condition of least square;
and 5: according to the difference between the re-optimized seismic record obtained by calculating the optimized reflection coefficient and the actual seismic record, on the premise that the 2 norm in the seismic wavelet optimization model is minimum, the optimized seismic wavelet is obtained by using a least square method;
step 6: and repeating the step 4-5 until the reflection coefficient difference is smaller than the reflection coefficient difference threshold or the seismic wavelet difference is smaller than the seismic wavelet difference threshold.
FIGS. 3a and 3b show a comparison of a wavelet extracted according to the prior art and a wavelet extracted according to the present method and a real wavelet, respectively, which are obtained by the present method with a smaller difference.
Fig. 4a and 4b are graphs showing the comparison of the deconvolution result obtained by the prior art and the present method with the true reflection coefficient, respectively, according to an embodiment of the present invention, the deconvolution result obtained by the present method being nearly coincident with the true reflection coefficient.
FIGS. 5a and 5b show graphs comparing a real wavelet to an estimated wavelet for 110ms and 700ms, respectively, for different moments in time, the extracted wavelet and the real wavelet are nearly coincident.
FIG. 6 shows a graph comparing a real wavelet to an estimated wavelet at Q49 according to one embodiment of the present invention.
Fig. 7 shows a graph comparing the true reflection coefficient to the estimated reflection coefficient at Q49 according to one embodiment of the present invention.
FIG. 8 shows a graph comparing a real wavelet to an estimated wavelet at a Q of 55 according to one embodiment of the present invention.
Fig. 9 shows a graph comparing the true reflection coefficient to the estimated reflection coefficient at a Q of 55 according to one embodiment of the present invention.
Fig. 6-9 show that if the Q value estimation has a small error, the method can also obtain good results, i.e. the method has a certain fault tolerance. Since the true Q value is 50, the results of the separate tests using Q49 and Q55 match well with both the true wavelet and the reflection coefficient.
FIG. 10 shows a comparison of a real wavelet to an estimated wavelet at 5% noise according to one embodiment of the present invention.
FIG. 11 shows a graph comparing the true reflection coefficient to the estimated reflection coefficient at 5% noise according to one embodiment of the invention.
10-11 show the noise immunity of the method, with the addition of 5% white Gaussian noise, and the estimated source wavelet with little phase shift.
In summary, the invention decomposes the unsteady state into the attenuation estimation and a steady state, performs the absorption attenuation estimation, and performs the interaction iteration of the reflection coefficient optimization and the source wavelet optimization, thereby obtaining a more accurate result, improving the seismic data resolution, better maintaining the relative amplitude of the reflection coefficient and the reflection coefficient phase deviation, and having good fault tolerance and noise immunity in the seismic data processing.
It will be appreciated by persons skilled in the art that the above description of embodiments of the invention is intended only to illustrate the benefits of embodiments of the invention and is not intended to limit embodiments of the invention to any examples given.
A time-varying sparse deconvolution system according to the present invention, having stored thereon a computer program which when executed by a processor performs the steps of: step 1: inputting a seismic record, initial seismic wavelets and sparsity; step 2: estimating a Q value vector, and constructing an attenuation matrix; and step 3: according to the attenuation matrix, establishing a reflection coefficient optimization model and a seismic wavelet optimization model; and 4, step 4: obtaining an optimized reflection coefficient according to the initial seismic wavelet, the sparsity and the reflection coefficient optimization model; and 5: obtaining optimized seismic wavelets according to the optimized reflection coefficient and the seismic wavelet optimization model; step 6: and repeating the step 4-5 until the reflection coefficient difference is smaller than the reflection coefficient difference threshold or the seismic wavelet difference is smaller than the seismic wavelet difference threshold.
In one example, the reflection coefficient optimization model is:
where r represents the reflection coefficient, y represents the seismic record, A (Q) is the attenuation matrix for the Q value, Q is the quality factor, and W is the seismic wavelet matrix.
In one example, the seismic wavelet optimization model is:
wherein R (w) is seismic wavelet constraint, w is seismic wavelet, y represents seismic record, Q is quality factor,
i is the matrix with the ith diagonal element being 1 and the remainder being 0.
In one example, step 4 obtains the optimized reflection coefficient through an OMP algorithm, including: and calculating and optimizing the seismic record according to the seismic record, the initial seismic wavelet and the sparsity, seeking an index with the largest atomic inner product in a matrix of the residual error and the initial seismic wavelet, updating an index set, and solving the optimized reflection coefficient of the reflection coefficient optimization model under the condition of least square.
In one example, step 5 comprises: and calculating the difference between the re-optimized seismic record and the actual seismic record according to the optimized reflection coefficient, and solving the optimized seismic wavelet by using a least square method on the premise of the minimum 2 norm in the seismic wavelet optimization model.
The system decomposes the unsteady state into attenuation estimation and a steady state, carries out absorption attenuation estimation, and carries out interactive iteration on reflection coefficient optimization and seismic source wavelet optimization, thereby obtaining a more accurate result, improving seismic data resolution, better keeping the relative amplitude of the reflection coefficient and the phase deviation of the reflection coefficient, and having good fault tolerance and noise immunity.
Having described embodiments of the present invention, the foregoing description is intended to be exemplary, not exhaustive, and not limited to the embodiments disclosed. Many modifications and variations will be apparent to those of ordinary skill in the art without departing from the scope and spirit of the described embodiments.
Claims (10)
1. A time-varying sparse deconvolution method, comprising:
step 1: inputting a seismic record, initial seismic wavelets and sparsity;
step 2: estimating a Q value vector, and constructing an attenuation matrix;
and step 3: establishing a reflection coefficient optimization model and a seismic wavelet optimization model according to the attenuation matrix;
and 4, step 4: obtaining an optimized reflection coefficient according to the initial seismic wavelet, the sparsity and the reflection coefficient optimization model;
and 5: obtaining an optimized seismic wavelet according to the optimized reflection coefficient and the seismic wavelet optimization model;
step 6: and repeating the step 4-5 until the reflection coefficient difference is smaller than the reflection coefficient difference threshold or the seismic wavelet difference is smaller than the seismic wavelet difference threshold.
2. The time-varying sparse deconvolution method of claim 1, wherein the reflection coefficient optimization model is:
where r represents the reflection coefficient, y represents the seismic record, A (Q) is the attenuation matrix for the Q value, Q is the quality factor, and W is the seismic wavelet matrix.
3. The time-varying sparse deconvolution method of claim 1, wherein the seismic wavelet optimization model is:
wherein R (w) is seismic wavelet constraint, w is seismic wavelet, y represents seismic record, Q is quality factor,
i is the matrix with the ith diagonal element being 1 and the remainder being 0.
4. The time-varying sparse deconvolution method of claim 1, wherein said step 4 of obtaining an optimized reflection coefficient by an OMP algorithm comprises:
and calculating an optimized seismic record according to the seismic record, the initial seismic wavelet and the sparsity, seeking an index with the largest atomic inner product in a matrix of the residual error and the initial seismic wavelet, updating an index set, and solving the optimized reflection coefficient of the reflection coefficient optimization model under the condition of least square.
5. The time-varying sparse deconvolution method of claim 3, wherein said step 5 comprises:
and calculating the difference between the re-optimized seismic record and the actual seismic record according to the optimized reflection coefficient, and solving the optimized seismic wavelet by using a least square method on the premise that the 2 norm in the seismic wavelet optimization model is the minimum.
6. A time-varying sparse deconvolution system having a computer program stored thereon, wherein the program when executed by a processor performs the steps of:
step 1: inputting a seismic record, initial seismic wavelets and sparsity;
step 2: estimating a Q value vector, and constructing an attenuation matrix;
and step 3: establishing a reflection coefficient optimization model and a seismic wavelet optimization model according to the attenuation matrix;
and 4, step 4: obtaining an optimized reflection coefficient according to the initial seismic wavelet, the sparsity and the reflection coefficient optimization model;
and 5: obtaining an optimized seismic wavelet according to the optimized reflection coefficient and the seismic wavelet optimization model;
step 6: and repeating the step 4-5 until the reflection coefficient difference is smaller than the reflection coefficient difference threshold or the seismic wavelet difference is smaller than the seismic wavelet difference threshold.
7. The time-varying sparse deconvolution system of claim 6, wherein the reflection coefficient optimization model is:
where r represents the reflection coefficient, y represents the seismic record, A (Q) is the attenuation matrix for the Q value, Q is the quality factor, and W is the seismic wavelet matrix.
8. The time-varying sparse deconvolution system of claim 6, wherein the seismic wavelet optimization model is:
wherein R (w) is seismic wavelet constraint, w is seismic wavelet, y represents seismic record, Q is quality factor,
i is the matrix with the ith diagonal element being 1 and the remainder being 0.
9. The time-varying sparse deconvolution system of claim 6, wherein said step 4 obtains optimized reflection coefficients by an OMP algorithm comprising:
and calculating an optimized seismic record according to the seismic record, the initial seismic wavelet and the sparsity, seeking an index with the largest atomic inner product in a matrix of the residual error and the initial seismic wavelet, updating an index set, and solving the optimized reflection coefficient of the reflection coefficient optimization model under the condition of least square.
10. The time-varying sparse deconvolution system of claim 8, wherein the step 5 comprises:
and calculating the difference between the re-optimized seismic record and the actual seismic record according to the optimized reflection coefficient, and solving the optimized seismic wavelet by using a least square method on the premise that the 2 norm in the seismic wavelet optimization model is the minimum.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810679707.8A CN110646841B (en) | 2018-06-27 | 2018-06-27 | Time-varying sparse deconvolution method and system |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810679707.8A CN110646841B (en) | 2018-06-27 | 2018-06-27 | Time-varying sparse deconvolution method and system |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110646841A true CN110646841A (en) | 2020-01-03 |
| CN110646841B CN110646841B (en) | 2021-05-25 |
Family
ID=69009119
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810679707.8A Active CN110646841B (en) | 2018-06-27 | 2018-06-27 | Time-varying sparse deconvolution method and system |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110646841B (en) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111880218A (en) * | 2020-07-13 | 2020-11-03 | 西南石油大学 | Construction Method of Inversion Wavelet Dictionary Based on Quality Factor |
| CN112305616A (en) * | 2020-09-23 | 2021-02-02 | 中国石油天然气集团有限公司 | Method and device for acquiring seismic data profile in optical fiber well |
| CN113219530A (en) * | 2020-02-05 | 2021-08-06 | 中国石油天然气股份有限公司 | Unsteady-state blind deconvolution method and device |
| CN113341463A (en) * | 2021-06-10 | 2021-09-03 | 中国石油大学(北京) | Pre-stack seismic data non-stationary blind deconvolution method and related components |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20100097888A1 (en) * | 2007-04-10 | 2010-04-22 | Exxonmobil Upstream Research Company | Separation and Noise Removal for Multiple Vibratory Source Seismic Data |
| US20140204711A1 (en) * | 2011-03-25 | 2014-07-24 | Saudi Arabian Oil Company | Simultaneous wavelet extraction and deconvolution processing in the time domain |
| CN105334532A (en) * | 2015-10-21 | 2016-02-17 | 中国石油化工股份有限公司 | Seismic wavelet estimation method |
| CN105467442A (en) * | 2015-12-09 | 2016-04-06 | 中国石油大学(北京) | A globally optimized time-varying sparse deconvolution method and apparatus |
| US20160116620A1 (en) * | 2014-10-24 | 2016-04-28 | Ion Geophysical Corporation | Methods and systems for seismic inversion and related seismic data processing |
-
2018
- 2018-06-27 CN CN201810679707.8A patent/CN110646841B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20100097888A1 (en) * | 2007-04-10 | 2010-04-22 | Exxonmobil Upstream Research Company | Separation and Noise Removal for Multiple Vibratory Source Seismic Data |
| US20140204711A1 (en) * | 2011-03-25 | 2014-07-24 | Saudi Arabian Oil Company | Simultaneous wavelet extraction and deconvolution processing in the time domain |
| US20160116620A1 (en) * | 2014-10-24 | 2016-04-28 | Ion Geophysical Corporation | Methods and systems for seismic inversion and related seismic data processing |
| CN105334532A (en) * | 2015-10-21 | 2016-02-17 | 中国石油化工股份有限公司 | Seismic wavelet estimation method |
| CN105467442A (en) * | 2015-12-09 | 2016-04-06 | 中国石油大学(北京) | A globally optimized time-varying sparse deconvolution method and apparatus |
Non-Patent Citations (4)
| Title |
|---|
| N. KAZEMI: "Blind Deconvolution with Toeplitz-structured Sparse Total Least Squares Algorithm", 《80TH EAGE CONFERENCE & EXHIBITION 2018》 * |
| 刘喜武，等: "基于带状混合矩阵ICA实现地震盲反褶积", 《地球物理学进展》 * |
| 李艳琴: "Bayesian框架下的地震盲反褶积算法研究", 《中国博士学位论文全文数据库》 * |
| 陈祖庆，等: "基于压缩感知的稀疏脉冲反射系数谱反演方法研究", 《石油物探》 * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113219530A (en) * | 2020-02-05 | 2021-08-06 | 中国石油天然气股份有限公司 | Unsteady-state blind deconvolution method and device |
| CN113219530B (en) * | 2020-02-05 | 2023-09-26 | 中国石油天然气股份有限公司 | Unsteady state blind deconvolution method and device |
| CN111880218A (en) * | 2020-07-13 | 2020-11-03 | 西南石油大学 | Construction Method of Inversion Wavelet Dictionary Based on Quality Factor |
| CN112305616A (en) * | 2020-09-23 | 2021-02-02 | 中国石油天然气集团有限公司 | Method and device for acquiring seismic data profile in optical fiber well |
| CN112305616B (en) * | 2020-09-23 | 2024-03-01 | 中国石油天然气集团有限公司 | Method and device for acquiring seismic data section in optical fiber well |
| CN113341463A (en) * | 2021-06-10 | 2021-09-03 | 中国石油大学(北京) | Pre-stack seismic data non-stationary blind deconvolution method and related components |
| CN113341463B (en) * | 2021-06-10 | 2023-05-26 | 中国石油大学(北京) | Non-stationary blind deconvolution method for pre-stack seismic data and related components |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110646841B (en) | 2021-05-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110646841B (en) | Time-varying sparse deconvolution method and system | |
| US9429668B2 (en) | Iterative dip-steering median filter for seismic data processing | |
| Huang et al. | Erratic noise suppression using iterative structure‐oriented space‐varying median filtering with sparsity constraint | |
| CN112882099B (en) | Earthquake frequency band widening method and device, medium and electronic equipment | |
| CN108828670B (en) | A kind of seismic data noise-reduction method | |
| CN108693555A (en) | Intelligent time-varying blind deconvolution wideband processing method and processing device | |
| CN113077386A (en) | Seismic data high-resolution processing method based on dictionary learning and sparse representation | |
| Li et al. | Seismic noise suppression using weighted nuclear norm minimization method | |
| Kuruguntla et al. | Erratic noise attenuation using double sparsity dictionary learning method | |
| CN111505709A (en) | Attenuation qualitative analysis method based on sparse spectral decomposition | |
| CN113341463B (en) | Non-stationary blind deconvolution method for pre-stack seismic data and related components | |
| CN113109873B (en) | Desert seismic signal noise suppression method based on rank residual error constraint | |
| CN107390261B (en) | Surface multiple and wavelet estimation method and system based on linear Bregman algorithm | |
| CN116859461B (en) | Multiple imaging method and system | |
| Hao et al. | Denoising Method Based on Spectral Subtraction in Time‐Frequency Domain | |
| Guo et al. | Seismic Noise Suppression Method Based on Wave-Unet and Attention Mechanism | |
| CN111190226B (en) | Three-dimensional seismic data surface wave noise suppression method | |
| CN111766624B (en) | A method, device, storage medium and electronic device for frequency extension processing of seismic data | |
| Zhou et al. | Flattening the seismic data for optimal noise attenuation | |
| CN114662045A (en) | Multidimensional seismic data denoising method based on p-order tensor deep learning of frame set | |
| CN109459788B (en) | Stratum quality factor calculation method and system | |
| Sang et al. | Noise attenuation of seismic data via deep multiscale fusion network | |
| Aghayan et al. | Noise suppression using a near-source wavelet | |
| CN118915142B (en) | A spatiotemporal joint sparse deconvolution method and system for enhancing geological body recognition accuracy | |
| CN109991657A (en) | High resolution processing method of seismic data based on inverse bisection recursive singular value decomposition |
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 |