CN113219530A - Unsteady-state blind deconvolution method and device - Google Patents
Unsteady-state blind deconvolution method and device Download PDFInfo
- Publication number
- CN113219530A CN113219530A CN202010080341.XA CN202010080341A CN113219530A CN 113219530 A CN113219530 A CN 113219530A CN 202010080341 A CN202010080341 A CN 202010080341A CN 113219530 A CN113219530 A CN 113219530A
- Authority
- CN
- China
- Prior art keywords
- reflection coefficient
- derivative
- curve
- unsteady
- inversion
- 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 75
- 230000006870 function Effects 0.000 claims description 56
- 238000012545 processing Methods 0.000 claims description 35
- 239000011159 matrix material Substances 0.000 claims description 17
- 238000004590 computer program Methods 0.000 claims description 16
- 238000003860 storage Methods 0.000 claims description 11
- 238000004364 calculation method Methods 0.000 claims description 10
- 238000009499 grossing Methods 0.000 claims description 10
- 238000005070 sampling Methods 0.000 claims description 9
- 238000004422 calculation algorithm Methods 0.000 claims description 8
- 238000004458 analytical method Methods 0.000 claims description 6
- 238000007781 pre-processing Methods 0.000 claims description 4
- 238000009795 derivation Methods 0.000 claims description 3
- 238000010586 diagram Methods 0.000 description 16
- 230000000694 effects Effects 0.000 description 12
- 230000008859 change Effects 0.000 description 10
- 238000004891 communication Methods 0.000 description 7
- 238000003705 background correction Methods 0.000 description 6
- 230000008569 process Effects 0.000 description 6
- 230000015572 biosynthetic process Effects 0.000 description 5
- 238000012360 testing method Methods 0.000 description 5
- 239000000872 buffer Substances 0.000 description 4
- 238000001914 filtration Methods 0.000 description 4
- 230000002238 attenuated effect Effects 0.000 description 3
- 238000009826 distribution Methods 0.000 description 3
- 239000010410 layer Substances 0.000 description 3
- 230000011218 segmentation Effects 0.000 description 3
- 238000006467 substitution reaction Methods 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 239000012530 fluid Substances 0.000 description 2
- 230000014509 gene expression Effects 0.000 description 2
- 239000011435 rock Substances 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 239000004215 Carbon black (E152) Substances 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000001413 cellular effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 230000001808 coupling effect Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 229930195733 hydrocarbon Natural products 0.000 description 1
- 150000002430 hydrocarbons Chemical class 0.000 description 1
- 238000003709 image segmentation Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000011065 in-situ storage Methods 0.000 description 1
- 239000011229 interlayer Substances 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000010295 mobile communication Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000035699 permeability Effects 0.000 description 1
- 230000000644 propagated effect Effects 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 238000013139 quantization Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
-
- 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
-
- 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention provides an unsteady blind deconvolution method and a device, wherein the method comprises the following steps: acquiring seismic data to be analyzed, and presetting an initial seismic wavelet in the seismic data; constructing an attenuation sub-wave group through a preset quality factor range and the initial seismic wavelets, and performing reflection coefficient inversion on the attenuation sub-wave group to obtain an inverted reflection coefficient; substituting the reflection coefficient into L1‑2After the norm function, calculating a second derivative to obtain a derivative curve; minimizing L in the derivative curve1‑2And obtaining a high-order derivative of the norm to obtain a quality factor value and a corresponding reflection coefficient inversion result.
Description
Technical Field
The invention relates to the field of geophysical exploration, in particular to a method based on L1-2An unsteady state blind deconvolution method and apparatus for norm high order derivative criterion.
Background
Due to the non-perfect elasticity of the stratum, the phenomena of amplitude energy attenuation, phase distortion and frequency band narrowing are generated when seismic waves propagate in the stratum. The quality factor Q is an important parameter describing the elastic properties of the stratum and is used for quantitatively evaluating the attenuation degree of seismic waves during underground propagation. The Q value may not only reflect subsurface fracture information, but may also be used to predict hydrocarbon resources. The use of Q values may better explain the AVO effect and may also be used to detect the saturation of a gas or fluid, the geometric distribution of the phases. The closure of fractures under different stresses produces a significant change in the Q of the rock, so that attenuating observations of a core sample from a well can provide information about the fracture state of the reservoir or provide a description of the in situ stress. Information about the lithology, pressure, permeability, etc. of the formation may also be predicted using the Q value. In contrast, the presence of formation attenuation makes subsurface imaging and interpretation more difficult, and a reliable estimate of attenuation helps to improve the resolution of the seismic data. Reflection coefficient inversion, inverse Q filtering, etc. are effective attenuation compensation methods, and have been widely developed and applied, but the processing effect of these attenuation compensation methods always depends on the accuracy of Q value estimation. In addition, the amplitude and phase changes caused by the intrinsic attenuation of the stratum are reflected in the wavelets, and the reflection coefficient is the inherent property of the stratum and does not change along with the change of the Q value. According to the wave theory, when the wavelets are folded with the reflection coefficient sets corresponding to the thin inter-layer, an interference phenomenon occurs. The interference effect can also cause the amplitude and phase of the seismic waveform to change, so that the estimated attenuation is the Q value of the coupling of the intrinsic attenuation and the interference effect, and further the subsequent unsteady state reflection coefficient inversion can be influenced. The Q value may indicate fluid information and the reflection coefficient may indicate formation information. For non-steady state data, the two are interacting. In order to strip the effect of mutual influence of the two, the invention is based on stratum sparsity hypothesis and utilizes L1-2Derivative of norm on order n (for example, n-2) to quantify the sparsity and Q of the reflection coefficient inversion resultThe relationship between them. And performing simultaneous processing of Q value estimation and deconvolution through a Q scanning strategy.
The traditional inversion method of the steady-state reflection coefficient, such as sparse layer inversion, basis pursuit inversion and the like, requires that the input seismic data are steady, and the obtained seismic data are all unsteady in the actual field data acquisition process. Therefore, if the stationary wavelet assumption is still used, the obtained seismic data processing result is inaccurate, and the subsequent seismic interpretation work is influenced. For unsteady seismic data, the traditional method is to perform attenuation compensation first and then perform steady-state reflection coefficient inversion, but the attenuation compensation needs a known Q value, and the traditional Q value estimation method can be mainly classified into the following three categories: the first category of methods mainly includes time domain methods and frequency domain methods. The second category of methods can be summarized as using pre-stack CMP data to estimate the Q value. A third class of methods are some time-frequency domain methods. The Q value estimation is always a hot point of research of geophysicists, and an attenuation model of the Q value estimation quantifies the process of energy attenuation and phase change when seismic waves propagate in a stratum, but further increases the inadequacy of solving an inverse problem. When seismic waves are propagated underground, factors such as scattering and spherical diffusion can cause the precision of the solved quality factor to be reduced. Taner and Treitel have listed and compared 17 methods for calculating Q value, most of these methods consider the attenuation inside the seismic wave, but the calculation accuracy and stability are limited by the early processing effect of the seismic data, the original quality, the rock attribute and the like, and all have respective application range and do not have universality. In fact, the interference of the stratum can also cause the waveform of the seismic wave to change, and the accuracy of the Q value estimation can also be influenced. Whether the estimated Q value is accurate or not can cause the attenuation degree of the wavelet at different time positions to be different from the real attenuation degree, so that the accuracy of the time-varying wavelet matrix is influenced, and finally the reflection coefficient obtained by the unsteady state deconvolution method is deviated.
Disclosure of Invention
The invention aims to provide a method based on L1-2The unsteady state blind deconvolution method and device of norm high-order derivative criterion are based on the angle of reflection coefficient inversion,and estimating the Q value of the quality factor through the deformation of the criterion function, and using the Q value for inversion of the unsteady state reflection coefficient to realize resolution processing of the post-stack seismic data.
To achieve the above object, the unsteady-state blind deconvolution method provided by the present invention specifically comprises: acquiring seismic data to be analyzed, and presetting an initial seismic wavelet in the seismic data; constructing an attenuation sub-wave group through a preset quality factor range and the initial seismic wavelets, and performing reflection coefficient inversion on the attenuation sub-wave group to obtain an inverted reflection coefficient; substituting the reflection coefficient into L1-2After the norm function, calculating a second derivative to obtain a derivative curve; minimizing a corresponding L in the derivative curve1-2And obtaining a quality factor value and a corresponding reflection coefficient inversion result by using the second derivative of the norm function.
In the above-mentioned unsteady-state blind deconvolution method, preferably, L is the above-mentioned1-2The norm function includes:
in the above formula, 0 & ltsigma & lt 1 is a weight factor, N is the number of time sampling points, r is an inverted reflection coefficient sequence, and riIs the ith element in r, Q is the quality factor, J (Q) is the sparsity criterion L1-2A norm function value.
In the above-mentioned unsteady-state blind deconvolution method, it is preferable that the inverted reflection coefficient is substituted into L1-2After the norm function, the second derivative is solved, and the derivation curve is obtained by: substituting the inverted reflection coefficient into L1-2Obtaining a scanning curve by a norm function; solving a second derivative of the scanning curve to obtain a second-order curve; and respectively smoothing the scanning curve and the second-order curve and then performing superposition processing to obtain a derivative curve.
In the above unsteady blind deconvolution method, preferably, the inverting the reflection coefficient of the attenuation sub-wave group to obtain an inverted reflection coefficient includes: and performing reflection coefficient inversion on the attenuation sub-wave group by using a sparse Bayesian learning algorithm to obtain an inverted reflection coefficient.
In the above-mentioned unsteady-state blind deconvolution method, preferably, the attenuation wavelet group includes:
in the above formula, M is the number of frequency sampling points, FHThe method is characterized in that the method is an inverse Fourier operator, W is a Toeplitz matrix formed by Fourier transform of initial wavelets, H is a conjugate transpose, s represents unsteady seismic record, A is a time-frequency domain attenuation matrix related to quality factors, r is a reflection coefficient sequence, and G (Q) is a time-varying wavelet matrix.
The invention also provides an unsteady blind deconvolution device, which comprises an acquisition module, an inversion module, a calculation module and an analysis module; the acquisition module is used for acquiring seismic data to be analyzed, and presetting an initial seismic wavelet in the seismic data; the inversion module is used for constructing an attenuation sub-wave group through a preset quality factor range and the initial seismic wavelets, and inverting the reflection coefficient of the attenuation sub-wave group to obtain an inverted reflection coefficient; the calculation module is used for substituting the reflection coefficient into L1-2After the norm function, calculating a second derivative to obtain a derivative curve; the analysis module is configured to minimize a corresponding L in the derivative curve1-2And obtaining a quality factor value and a corresponding reflection coefficient inversion result by using the second derivative of the norm function.
In the above unsteady blind deconvolution device, preferably, the calculation module further includes a preprocessing unit, and the preprocessing unit is configured to substitute the inverted reflection coefficient into L1-2Obtaining a scanning curve by a norm function; solving a second derivative of the scanning curve to obtain a second-order curve; and respectively smoothing the scanning curve and the second-order curve and then performing superposition processing to obtain a derivative curve.
In the above unsteady blind deconvolution device, preferably, the inversion module includes: and performing reflection coefficient inversion on the attenuation sub-wave group by using a sparse Bayesian learning algorithm to obtain an inverted reflection coefficient.
The invention also provides an electronic device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, the processor implementing the above method when executing the computer program.
The present invention also provides a computer-readable storage medium storing a computer program for executing the above method.
The invention has the beneficial technical effects that: under the framework of a Bayes learning method, through constructing a time-varying sub-wave group and L1-2And quantizing the internal relation between the sparsity and the Q value of the reflection coefficient inversion result by the norm high-order derivative criterion function, estimating the Q value or the time-varying wavelet, and completing the inversion of the unsteady reflection coefficient so as to further achieve the purpose of improving the seismic resolution. The method can be widely applied to seismic data processing, inversion and interpretation, and provides technical support for improving the resolution of non-stationary seismic data.
Drawings
The accompanying drawings, which are included to provide a further understanding of the invention and are incorporated in and constitute a part of this application, illustrate embodiment(s) of the invention and together with the description serve to explain the principles of the invention. In the drawings:
FIG. 1 is a schematic flow chart of an unsteady-state blind deconvolution method according to an embodiment of the present invention;
fig. 2 is a schematic diagram of a derivation curve calculation process according to an embodiment of the present invention;
fig. 3 is a schematic flow chart of an application of the unsteady-state blind deconvolution method in practical work according to an embodiment of the present invention;
FIG. 4 is a schematic flow chart illustrating a schematic flow chart of an unsteady-state blind deconvolution method according to an embodiment of the present invention;
fig. 5A to 5F are schematic diagrams illustrating a statistical numerical model testing the unsteady-state blind deconvolution method according to an embodiment of the present invention;
fig. 6A to fig. 6D are schematic diagrams illustrating a test of applying an unsteady-state blind deconvolution method to three-dimensional physical model data according to an embodiment of the present invention;
fig. 7A to 7D are schematic diagrams illustrating an application of the unsteady-state blind deconvolution method to three-dimensional actual data according to an embodiment of the present invention;
fig. 8 is a schematic structural diagram of an unsteady-state blind deconvolution device according to an embodiment of the present invention;
fig. 9 is a schematic structural diagram of an electronic device according to an embodiment of the present invention.
Detailed Description
The following detailed description of the embodiments of the present invention will be provided with reference to the drawings and examples, so that how to apply the technical means to solve the technical problems and achieve the technical effects can be fully understood and implemented. It should be noted that, unless otherwise specified, the embodiments and features of the embodiments of the present invention may be combined with each other, and the technical solutions formed are within the scope of the present invention.
Additionally, the steps illustrated in the flow charts of the figures may be performed in a computer system such as a set of computer-executable instructions and, although a logical order is illustrated in the flow charts, in some cases, the steps illustrated or described may be performed in an order different than here.
Referring to fig. 1, the unsteady-state blind deconvolution method provided by the present invention specifically includes: s101, acquiring seismic data to be analyzed, and presetting an initial seismic wavelet in the seismic data; s102, constructing an attenuation sub-wave group through a preset quality factor range and the initial seismic wavelets, and performing reflection coefficient inversion on the attenuation sub-wave group to obtain an inverted reflection coefficient; s103, substituting the reflection coefficient into L1-2After the norm function, calculating a second derivative to obtain a derivative curve; s104 minimizing the corresponding L in the derivative curve1-2And obtaining a quality factor value and a corresponding reflection coefficient inversion result by using the second derivative of the norm function. Referring to fig. 2, in step S103, the inverted reflection coefficient is substituted into L1-2After the norm function, calculating a second derivative, and obtaining a derivative curve further comprises: s201 substituting the inverted reflection coefficient into L1-2Obtaining a scanning curve by a norm function; s202 pair the scanningThe second derivative is obtained by the curve to obtain a second-order curve; s203, smoothing the scanning curve and the second-order curve respectively, and then performing superposition processing to obtain a derivative curve. In step S102, performing reflection coefficient inversion on the attenuation sub-wave group to obtain an inverted reflection coefficient may include: and performing reflection coefficient inversion on the attenuation sub-wave group by using a sparse Bayesian learning algorithm to obtain an inverted reflection coefficient.
In the above embodiment, L1-2The norm function may be:
in the above formula, 0 & ltsigma & lt 1, N is the number of time sampling points, r is the inverted reflection coefficient sequence, riIs the ith element in r, Q is the quality factor, J (Q) is the sparsity criterion L1-2A norm function value.
The attenuating wavelet group comprises:
in the above formula, M is the number of frequency sampling points, FHThe method is characterized in that the method is an inverse Fourier operator, W is a diagonal matrix formed by Fourier transform of initial wavelets, H is a conjugate transpose, s represents an attenuation seismic record, A is a time-frequency domain attenuation matrix related to a quality factor, r is an inverted reflection coefficient sequence, and G (Q) is an attenuation wavelet matrix.
Specifically, the unsteady blind deconvolution method provided by the present invention mainly comprises two parts in principle: unsteady seismic data forward modeling based on wave equation and based on L1-2Estimating the Q value of norm high-order derivative criterion;
the forward modeling of the unsteady seismic data based on the wave equation mainly comprises the following steps:
for a plane wave propagating along the z-axis, the analytical expression of the one-dimensional wave equation is as follows:
where x is a displacement function, ω is an angular frequency, t is time, i is an imaginary unit, and k is a wave number. When seismic waves propagate in an underground medium, viscoelastic absorption attenuation and dispersion can occur; one way to describe this is to express the wave number k as a real wave number kr(ω) and a product related to the quality factor Q, expressed as:
in the field of seismic exploration, convolution models are widely used to describe the propagation of seismic waves in the subsurface. The attenuation characteristic of the stratum is not considered in the conventional convolution model, so that in one embodiment of the invention, a final unsteady seismic data simulation formula can be derived mainly according to the wave equation theory of the sound wave and the deformation of the formula (1) and the formula (2):
where s represents the decaying seismic record, α (ω, τ)z,Qz) Representing an exponential function as a function of frequency and time, r (z) being the reflection coefficient at zero angle at depth z; further, the unsteady convolution model in the time domain can be described as:
wherein: m is the number of frequency sampling points,
when the Q value is infinite, alpha approaches to 1, namely the traditional steady state deconvolution formula is obtained; g (Q) is an attenuating wavelet matrix. As can be seen from equation (4), the unsteady seismic record can be regarded as a weighted linear superposition of N time-varying wavelets with weights rn. Due to the band-limited nature of seismic data, and the introduction of Q values, direct solution of equation (4) is ill-posed; therefore, the invention is mainly based on stratum sparse assumption, adopts a Sparse Bayesian Learning (SBL) method, and obtains the reflection coefficient sequence through inversion.
With respect to based on L1-2The Q value estimation of the norm high order derivative criterion mainly comprises:
amplitude and phase changes caused by the intrinsic attenuation of the formation are reflected in the attenuated seismic wavelets, while the reflection coefficient does not change with the change in the Q value. In fact, when the wavelet is convoluted with the reflection coefficient sequence corresponding to the thin interbed, an interference effect is generated, and the interference effect also causes the amplitude and the phase of the seismic waveform to change, so that the estimated Q value is a pseudo value generated by the coupling effect of the formation intrinsic attenuation and the interference effect, and the subsequent unsteady state reflection coefficient inversion can be influenced. Meanwhile, based on the sparse hypothesis, correlation exists between the Q value and the sparsity of the reflection coefficient inversion result, and when the Q value is not accurate, the inversion result is different from a real model, which is mainly expressed in pulse amplitude and non-zero pulse number; when the attenuation wavelet matrix constructed by smaller Q value is used for inversion, the attenuation of the deep wavelet is overlarge. In order to match the waveforms, the amplitude of the deep reflection coefficient pulse is correspondingly increased, and the number of pulses is increased; when the attenuation wavelet matrix constructed by a larger Q value is used for inversion, the wavelet attenuation is reduced, so that a plurality of homopolar small pulses appear; for this reason, for the above phenomena, the present invention performs quantization by a sparseness criterion function, specifically:
whether the Q value estimation is accurate or not directly influences the accuracy and precision of reflection coefficient inversion, and further influences subsequent high-resolution processing and interpretation. Thus, with L1-2The nth derivative of the norm (where n is 2 in the present invention) quantifies the relationship between the Q value and the inverted reflection coefficient, L1-2The norm may be expressed as:
wherein: 0 < sigma < 1, N represents the number of time sampling points, r represents the inverted sequence of reflection coefficients, riIs the ith element in r. However, L1-2The norm does not characterize the correct Q value. The invention finds that the salient points of the second derivative can correctly indicate the position and the size of the Q value. The expressions for its first and second derivatives are:
L1-2the second derivative of the norm has a good convex function rule along with the change of the Q value. Therefore, the Q value is estimated by solving the minimum of the formula (7), the optimal reflection coefficient is inverted at the same time, the purpose of decoupling the inherent attenuation and the interference attenuation is achieved, and finally, a high-precision Q field and a high-resolution inversion result are obtained at the same time.
In combination with the above embodiments, the present invention provides an L-based solution1-2Unsteady-state blind deconvolution method of norm high-order derivative criterionThe application flow in actual work is specifically shown in fig. 3:
s301, inputting attenuation seismic data;
s302, estimating an initial seismic wavelet from seismic data;
s303, constructing the time-varying sub-wave group by adopting a Q value scanning strategy. Performing reflection coefficient inversion by using a sparse Bayesian learning method, and eliminating time-varying wavelets;
s304, substituting the reflection coefficient obtained by inversion into L1-2And obtaining a scanning curve by using the norm function. Solving a second derivative of the scanning curve to obtain a corresponding curve; in order to reduce the influence of singular values and overcome the instability of single-channel inversion, smoothing and overlaying are respectively carried out on a scanning curve and a curve obtained after a second derivative is solved;
s305 minimization of L1-2And obtaining the high-order derivative of the norm to obtain an estimated Q value and a corresponding reflection coefficient inversion result.
Referring to fig. 4, in the application process of the unsteady-state blind deconvolution method provided by the present invention, first, initial wavelets are estimated from the input attenuated seismic data, a range of the sweep Q is given, and the two are combined to construct an attenuated wavelet group for unsteady-state reflection coefficient inversion; secondly, substituting the reflection coefficient obtained by inversion into a criterion function to obtain a second derivative; and then, minimizing the second derivative, namely obtaining the required Q value, and finally obtaining a high-precision Q field and a high-resolution inversion result at the same time. Based on the flow, the method uses a statistical numerical model, three-dimensional physical model data and three-dimensional actual data to carry out testing respectively; specifically, through 1000 pre-designed geological models, each model has 30 non-zero random reflection coefficient pulses, the maximum interval between the non-zero pulses does not exceed 30ms, the amplitude of the non-zero pulses is set between-0.3 and 0.3 for practical purposes, and the total length of the model is 1 s. Selecting a Rake wavelet with a main frequency of 30Hz as an initial wavelet, setting the Q value to be 50, constructing an attenuation seismic wavelet matrix, and respectively convolving the attenuation seismic wavelet matrix with 1000 reflection coefficient models to obtain 1000 unstable seismic data. Referring to fig. 5A to 5F, fig. 5A is a comparison graph of the original data and the synthesized record of the inversion result, and fig. 5B is an original inverse recordAnd (3) a comparison graph of the reflection coefficient and the inversion reflection coefficient. As can be seen from fig. 5A, the inversion results synthesized seismic data that closely matched the original data. On the premise of data matching, the reflection coefficient obtained by inversion is basically consistent with the original reflection coefficient, as shown in fig. 5B; FIG. 5C shows the substitution of the inversion reflection coefficient results for each corresponding Q value into L1-2The norm criterion function obtains 1000 scan curves, and fig. 5D is a superposition result of the 1000 scan curves. As can be seen from fig. 5D, the superimposed curve has no obvious convexity and cannot indicate the location of the correct Q value; fig. 5E is a second derivative of the scan curve of fig. 5C, and fig. 5F is a superimposed curve of fig. 5E. As is apparent from FIG. 5F, by finding L1-2The second derivative of the norm criterion function is obviously convex, and the position of the inflection point just corresponds to the position of the real Q value. As can be seen from a statistical numerical test, the method can accurately estimate the Q value, is not easily influenced by a reflection coefficient model, and has high applicability and accuracy.
Referring to fig. 6A to 6D, the unsteady-state blind deconvolution method is applied to three-dimensional physical model data for testing. FIG. 6A shows the substitution of the inversion reflection coefficient results for each corresponding Q value into L1-2In the norm criterion function, a curve after superposition; FIG. 6B is a graph of FIG. 6A after the second derivative is taken; FIG. 6C is the original three-dimensional physical model data; figure 6D shows the high resolution processing results obtained with the inventive process. Comparing fig. 6C and fig. 6D, the processed results have clearer layer interfaces and show more construction details, and are better consistent with the real physical simulation.
Next, referring to fig. 7A to 7D, the effectiveness of the steady-state blind deconvolution method is tested by using three-dimensional actual data, and fig. 7A is a graph obtained by substituting the inversion reflection coefficient result corresponding to each Q value into L1-2In the norm criterion function, a curve after superposition; FIG. 7B is a graph of FIG. 7A after the second derivative is taken; FIG. 7C is an amplitude slice along the layer of the raw actual data; FIG. 7D shows the result of the high resolution processing, and comparing the original seismic record with the result of the high resolution processing shows that the processed edge and end of the river are more visibleClearly, the seismic resolution is improved significantly.
Referring to fig. 8, the present invention further provides an unsteady-state blind deconvolution device, which includes an acquisition module, an inversion module, a calculation module, and an analysis module; the acquisition module is used for acquiring seismic data to be analyzed, and presetting an initial seismic wavelet in the seismic data; the inversion module is used for constructing an attenuation sub-wave group through a preset quality factor range and the initial seismic wavelets, and inverting the reflection coefficient of the attenuation sub-wave group to obtain an inverted reflection coefficient; the calculation module is used for substituting the reflection coefficient into L1-2After the norm function, calculating a second derivative to obtain a derivative curve; the analysis module is configured to minimize a corresponding L in the derivative curve1-2And obtaining a quality factor value and a corresponding reflection coefficient inversion result by using the second derivative of the norm function.
In the above-mentioned unsteady blind deconvolution device, the calculation module may further include a preprocessing unit for substituting the inverted reflection coefficient into L1-2Obtaining a scanning curve by a norm function; solving a second derivative of the scanning curve to obtain a second-order curve; and respectively smoothing the scanning curve and the second-order curve and then performing superposition processing to obtain a derivative curve. The inversion module may include: and performing reflection coefficient inversion on the attenuation sub-wave group by using a sparse Bayesian learning algorithm to obtain an inverted reflection coefficient.
The invention also provides an electronic device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, the processor implementing the above method when executing the computer program.
The present invention also provides a computer-readable storage medium storing a computer program for executing the above method.
The invention has the beneficial technical effects that: under the framework of a Bayes learning method, through constructing a time-varying sub-wave group and L1-2Norm high-order derivative criterion function for quantifying sparsity and Q value of reflection coefficient inversion resultAnd estimating a Q value or a time-varying wavelet, and completing inversion of an unsteady state reflection coefficient at the same time, thereby achieving the purpose of improving the seismic resolution. The method can be widely applied to seismic data processing, inversion and interpretation, and provides technical support for improving the resolution of non-stationary seismic data.
Fig. 9 is a schematic block diagram of a system configuration of an electronic apparatus 600 according to an embodiment of the present invention. As shown in fig. 9, the electronic device 600 may include a central processor 100 and a memory 140; the memory 140 is coupled to the central processor 100. Notably, this diagram is exemplary; other types of structures may also be used in addition to or in place of the structure to implement telecommunications or other functions.
In one embodiment, the CT image shadow correction function may be integrated into the central processor 100. The central processor 100 may be configured to control as follows: carrying out image texture removal operation on the original CT reconstructed image to obtain a smooth image; carrying out segmentation processing on the structural components of the original CT reconstructed image according to human tissues to construct a template image; and carrying out shading correction according to the smooth image and the template image.
The method for removing image texture from an original CT reconstructed image to obtain a smooth image comprises the following steps: and performing edge protection and image texture removal on the original CT reconstructed image by using an L0 norm smoothing algorithm to obtain the smoothed image.
The method comprises the following steps of carrying out segmentation processing on structural components of the original CT reconstructed image according to human tissues to construct a template image, wherein the segmentation processing comprises the following steps: segmenting the original CT reconstructed image into a plurality of human tissue regions by adopting an image segmentation method; and filling CT values of corresponding tissues under the voltage of the X-ray bulb tube in different human tissue areas respectively to obtain the template image.
Wherein, the shading correction is carried out according to the smooth image and the template image, and the shading correction comprises the following steps: subtracting the smooth image from the template image to obtain a residual image, wherein the residual image comprises image shadows and organizational structure errors; carrying out low-pass filtering processing on the organizational structure error of the residual image to obtain CT image shadow distribution; and performing compensation processing on the original CT reconstructed image by utilizing the CT image shadow distribution to obtain a corrected CT image.
Wherein, the low-pass filtering processing is carried out on the organizational structure error of the residual image, and comprises the following steps: and carrying out low-pass filtering processing on the tissue structure error of the residual image by utilizing a Savitzky-Golay local low-pass filter.
The image texture removal of the original CT reconstructed image by using an L0 norm smoothing algorithm comprises the following steps:
calculating a texture-free smoothed image by an objective function, the objective function being as follows:
wherein the content of the first and second substances,
Spfor the p-th pixel index, I, of the texture-free smoothed image SpReconstruct the p-th pixel index of image I for the original CT,
the number of pixel indices p, λ is the smoothing factor,
and
the partial derivatives in both the x and y directions are indexed for the pixel.
In another embodiment, the CT image shading correction apparatus may be configured separately from the central processing unit 100, for example, the CT image shading correction apparatus may be configured as a chip connected to the central processing unit 100, and the CT image shading correction function is realized by the control of the central processing unit.
As shown in fig. 9, the electronic device 600 may further include: communication module 110, input unit 120, audio processing unit 130, display 160, power supply 170. It is noted that the electronic device 600 does not necessarily include all of the components shown in FIG. 9; furthermore, the electronic device 600 may also comprise components not shown in fig. 9, which may be referred to in the prior art.
As shown in fig. 9, the central processor 100, sometimes referred to as a controller or operational control, may include a microprocessor or other processor device and/or logic device, the central processor 100 receiving input and controlling the operation of the various components of the electronic device 600.
The memory 140 may be, for example, one or more of a buffer, a flash memory, a hard drive, a removable media, a volatile memory, a non-volatile memory, or other suitable device. The information relating to the failure may be stored, and a program for executing the information may be stored. And the central processing unit 100 may execute the program stored in the memory 140 to realize information storage or processing, etc.
The input unit 120 provides input to the cpu 100. The input unit 120 is, for example, a key or a touch input device. The power supply 170 is used to provide power to the electronic device 600. The display 160 is used to display an object to be displayed, such as an image or a character. The display may be, for example, an LCD display, but is not limited thereto.
The memory 140 may be a solid state memory such as Read Only Memory (ROM), Random Access Memory (RAM), a SIM card, or the like. There may also be a memory that holds information even when power is off, can be selectively erased, and is provided with more data, an example of which is sometimes called an EPROM or the like. The memory 140 may also be some other type of device. Memory 140 includes buffer memory 141 (sometimes referred to as a buffer). The memory 140 may include an application/function storage section 142, and the application/function storage section 142 is used to store application programs and function programs or a flow for executing the operation of the electronic device 600 by the central processing unit 100.
The memory 140 may also include a data store 143, the data store 143 for storing data, such as contacts, digital data, pictures, sounds, and/or any other data used by the electronic device. The driver storage portion 144 of the memory 140 may include various drivers of the electronic device for communication functions and/or for performing other functions of the electronic device (e.g., messaging application, address book application, etc.).
The communication module 110 is a transmitter/receiver 110 that transmits and receives signals via an antenna 111. The communication module (transmitter/receiver) 110 is coupled to the central processor 100 to provide an input signal and receive an output signal, which may be the same as in the case of a conventional mobile communication terminal.
Based on different communication technologies, a plurality of communication modules 110, such as a cellular network module, a bluetooth module, and/or a wireless local area network module, may be provided in the same electronic device. The communication module (transmitter/receiver) 110 is also coupled to a speaker 131 and a microphone 132 via an audio processor 130 to provide audio output via the speaker 131 and receive audio input from the microphone 132 to implement general telecommunications functions. Audio processor 130 may include any suitable buffers, decoders, amplifiers and so forth. In addition, an audio processor 130 is also coupled to the central processor 100, so that recording on the local can be enabled through a microphone 132, and so that sound stored on the local can be played through a speaker 131.
As will be appreciated by one skilled in the art, embodiments of the present invention may be provided as a method, system, or computer program product. Accordingly, the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment or an embodiment combining software and hardware aspects. Furthermore, the present invention may take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, and the like) having computer-usable program code embodied therein.
The present invention is described with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each flow and/or block of the flow diagrams and/or block diagrams, and combinations of flows and/or blocks in the flow diagrams and/or block diagrams, can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.
These computer program instructions may also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instruction means which implement the function specified in the flowchart flow or flows and/or block diagram block or blocks.
These computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.
The above-mentioned embodiments are intended to illustrate the objects, technical solutions and advantages of the present invention in further detail, and it should be understood that the above-mentioned embodiments are only exemplary embodiments of the present invention, and are not intended to limit the scope of the present invention, and any modifications, equivalent substitutions, improvements and the like made within the spirit and principle of the present invention should be included in the scope of the present invention.
Claims (10)
1. An unsteady-state blind deconvolution method, comprising:
acquiring seismic data to be analyzed, and presetting an initial seismic wavelet in the seismic data;
constructing an attenuation sub-wave group through a preset quality factor range and the initial seismic wavelets, and performing reflection coefficient inversion on the attenuation sub-wave group to obtain an inverted reflection coefficient;
substituting the reflection coefficient into L1-2After the norm function, calculating a second derivative to obtain a derivative curve;
minimizing L corresponding to the derivative curve1-2And obtaining a quality factor value and a corresponding reflection coefficient inversion result by using the second derivative of the norm function.
2. The unsteady-state blind deconvolution method of claim 1, wherein L is1-2The norm function includes:
in the above formula, 0 & ltsigma & lt 1 is a weight factor, N is the number of time sampling points, r is an inverted reflection coefficient sequence, and riIs the ith element in r, Q is the quality factor, J (Q) is the sparsity criterion L1-2A norm function value.
3. The unsteady-state blind deconvolution method of claim 1, characterized in that the inverted reflection coefficients are substituted into L1-2After the norm function, the second derivative is solved, and the derivation curve is obtained by:
substituting the inverted reflection coefficient into L1-2Obtaining a scanning curve by a norm function;
solving a second derivative of the scanning curve to obtain a corresponding second-order curve;
and respectively smoothing the scanning curve and the second-order curve and then performing superposition processing to obtain a derivative curve.
4. The unsteady-state blind deconvolution method of claim 1, wherein performing reflection coefficient inversion on the attenuation sub-wave packet to obtain an inverted reflection coefficient comprises: and performing reflection coefficient inversion on the attenuation sub-wave group by using a sparse Bayesian learning algorithm to obtain an inverted reflection coefficient.
5. The unsteady-state blind deconvolution method of claim 4, wherein the attenuating wavelet group comprises:
in the above formula, M is the number of frequency sampling points, FHThe method is characterized in that the method is an inverse Fourier operator, W is a Toeplitz matrix formed by Fourier transform of initial wavelets, H is a conjugate transpose, s represents unsteady seismic record, A is a time-frequency domain attenuation matrix related to quality factors, r is a reflection coefficient sequence, and G (Q) is a time-varying wavelet matrix.
6. An unsteady blind deconvolution device is characterized by comprising an acquisition module, an inversion module, a calculation module and an analysis module;
the acquisition module is used for acquiring seismic data to be analyzed, and presetting an initial seismic wavelet in the seismic data;
the inversion module is used for constructing an attenuation sub-wave group through a preset quality factor range and the initial seismic wavelets, and inverting the reflection coefficient of the attenuation sub-wave group to obtain an inverted reflection coefficient;
the calculation module is used for substituting the reflection coefficient into L1-2After the norm function, calculating a second derivative to obtain a derivative curve;
the analysis module is configured to minimize a corresponding L in the derivative curve1-2And obtaining a quality factor value and a corresponding reflection coefficient inversion result by using the second derivative of the norm function.
7. The unsteady-state blind deconvolution device of claim 6, wherein the computation module further comprises a preprocessing unit for substituting the inverted reflection coefficients into L1-2Obtaining a scanning curve by a norm function; obtaining a second derivative of the scan curveObtaining a second-order curve; and respectively smoothing the scanning curve and the second-order curve and then performing superposition processing to obtain a derivative curve.
8. The unsteady blind deconvolution device of claim 6, wherein the inversion module comprises: and performing reflection coefficient inversion on the attenuation sub-wave group by using a sparse Bayesian learning algorithm to obtain an inverted reflection coefficient.
9. An electronic device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, wherein the processor implements the method of any of claims 1 to 5 when executing the computer program.
10. A computer-readable storage medium, characterized in that the computer-readable storage medium stores a computer program for executing the method of any one of claims 1 to 5.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010080341.XA CN113219530B (en) | 2020-02-05 | 2020-02-05 | Unsteady state blind deconvolution method and device |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010080341.XA CN113219530B (en) | 2020-02-05 | 2020-02-05 | Unsteady state blind deconvolution method and device |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113219530A true CN113219530A (en) | 2021-08-06 |
| CN113219530B CN113219530B (en) | 2023-09-26 |
Family
ID=77085568
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010080341.XA Active CN113219530B (en) | 2020-02-05 | 2020-02-05 | Unsteady state blind deconvolution method and device |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113219530B (en) |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108693555A (en) * | 2018-05-16 | 2018-10-23 | 中国石油大学(北京) | Intelligent time-varying blind deconvolution wideband processing method and processing device |
| CN110471113A (en) * | 2019-08-01 | 2019-11-19 | 中国石油大学(北京) | Bearing calibration, device and storage medium are moved in inverting based on unstable state seismic data |
| CN110542924A (en) * | 2019-09-02 | 2019-12-06 | 成都理工大学 | High-precision longitudinal and transverse wave impedance inversion method |
| CN110542923A (en) * | 2019-09-02 | 2019-12-06 | 成都理工大学 | Rapid high-precision post-stack seismic impedance inversion method |
| CN110646841A (en) * | 2018-06-27 | 2020-01-03 | 中国石油化工股份有限公司 | Time-varying sparse deconvolution method and system |
-
2020
- 2020-02-05 CN CN202010080341.XA patent/CN113219530B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108693555A (en) * | 2018-05-16 | 2018-10-23 | 中国石油大学(北京) | Intelligent time-varying blind deconvolution wideband processing method and processing device |
| CN110646841A (en) * | 2018-06-27 | 2020-01-03 | 中国石油化工股份有限公司 | Time-varying sparse deconvolution method and system |
| CN110471113A (en) * | 2019-08-01 | 2019-11-19 | 中国石油大学(北京) | Bearing calibration, device and storage medium are moved in inverting based on unstable state seismic data |
| CN110542924A (en) * | 2019-09-02 | 2019-12-06 | 成都理工大学 | High-precision longitudinal and transverse wave impedance inversion method |
| CN110542923A (en) * | 2019-09-02 | 2019-12-06 | 成都理工大学 | Rapid high-precision post-stack seismic impedance inversion method |
Non-Patent Citations (1)
| Title |
|---|
| 周珺 等: "基于l1-l2范数极小化的稀疏信号重建条件", 《合肥工业大学学报（自然科学版）》, vol. 43, no. 1, pages 137 - 140 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113219530B (en) | 2023-09-26 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US9702993B2 (en) | Multi-parameter inversion through offset dependent elastic FWI | |
| Groos et al. | The role of attenuation in 2D full-waveform inversion of shallow-seismic body and Rayleigh waves | |
| RU2422856C2 (en) | Recording method of seismic pulse stretching in seismic data | |
| EP3028071B1 (en) | Method and device for the generation and application of anisotropic elastic parameters in horizontal transverse isotropic (hti) media | |
| CN106405651B (en) | Full waveform inversion initial velocity model construction method based on logging matching | |
| US20050273266A1 (en) | Seismic event correlation and Vp-Vs estimation | |
| WO2010075412A2 (en) | Automatic dispersion extraction of multiple time overlapped acoustic signals | |
| CN105467442B (en) | The time-varying sparse deconvolution method and device of global optimization | |
| RU2570827C2 (en) | Hybrid method for full-waveform inversion using simultaneous and sequential source method | |
| Liu et al. | Effects of conjugate gradient methods and step-length formulas on the multiscale full waveform inversion in time domain: Numerical experiments | |
| CN110687597B (en) | Wave impedance inversion method based on joint dictionary | |
| Djebbi et al. | Frequency domain multiparameter acoustic inversion for transversely isotropic media with a vertical axis of symmetry | |
| Bindi et al. | Seismic input motion determined from a surface–downhole pair of sensors: A constrained deconvolution approach | |
| Yang et al. | Denoising low SNR percussion acoustic signal in the marine environment based on the LMS algorithm | |
| Tejero et al. | Appraisal of instantaneous phase-based functions in adjoint waveform inversion | |
| CN113219530B (en) | Unsteady state blind deconvolution method and device | |
| CN114740528A (en) | Pre-stack multi-wave joint inversion method based on ultramicro Laplace block constraint | |
| CN111025393B (en) | Reservoir prediction method, device, equipment and medium for stratum containing thin coal seam | |
| CN113009579B (en) | Seismic data inversion method and device | |
| CN112649848A (en) | Method and apparatus for solving seismic wave impedance using wave equation | |
| CN112363244B (en) | Wave impedance inversion method and carbonate heterogeneous reservoir prediction method and system | |
| CN111077573A (en) | Method, device and system for determining stratum elastic parameters | |
| Li et al. | Robust Q-compensated multidimensional impedance inversion using seislet-domain shaping regularization | |
| CN117434592B (en) | Seismic data processing method and device and electronic equipment | |
| CN114114406B (en) | Reservoir permeability estimation method and device |
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 |