CN113219530B - Unsteady state blind deconvolution method and device - Google Patents
Unsteady state blind deconvolution method and device Download PDFInfo
- Publication number
- CN113219530B CN113219530B CN202010080341.XA CN202010080341A CN113219530B CN 113219530 B CN113219530 B CN 113219530B CN 202010080341 A CN202010080341 A CN 202010080341A CN 113219530 B CN113219530 B CN 113219530B
- Authority
- CN
- China
- Prior art keywords
- reflection coefficient
- derivative
- obtaining
- curve
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 66
- 230000002238 attenuated effect Effects 0.000 claims abstract description 37
- 230000006870 function Effects 0.000 claims description 56
- 239000011159 matrix material Substances 0.000 claims description 20
- 238000004590 computer program Methods 0.000 claims description 16
- 238000005070 sampling Methods 0.000 claims description 11
- 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
- 238000004422 calculation algorithm Methods 0.000 claims description 8
- 238000004458 analytical method Methods 0.000 claims description 5
- 230000009466 transformation Effects 0.000 claims description 5
- 238000007781 pre-processing Methods 0.000 claims description 4
- 238000012545 processing Methods 0.000 description 20
- 238000010586 diagram Methods 0.000 description 15
- 230000000694 effects Effects 0.000 description 12
- 230000008859 change Effects 0.000 description 7
- 238000004891 communication Methods 0.000 description 7
- 238000003705 background correction Methods 0.000 description 6
- 230000015572 biosynthetic process Effects 0.000 description 6
- 238000001914 filtration Methods 0.000 description 5
- 230000008569 process Effects 0.000 description 5
- 238000012360 testing method Methods 0.000 description 5
- 239000000872 buffer Substances 0.000 description 4
- 238000009826 distribution 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
- 239000004973 liquid crystal related substance Substances 0.000 description 2
- 230000001902 propagating effect 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
- 238000007796 conventional method 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
- 230000001419 dependent effect Effects 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
- 238000003384 imaging method Methods 0.000 description 1
- 238000011065 in-situ storage Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000000463 material Substances 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
- 230000001568 sexual effect Effects 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/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 application provides an unsteady state 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 attenuated sub-wave group through a preset quality factor range and the initial seismic wavelet, and inverting the reflection coefficient of the attenuated sub-wave group to obtain an inverted reflection coefficient; substituting the reflection coefficient into L 1‑2 Obtaining a derivative curve by obtaining a second derivative after the norm function; minimizing L in the derivative curve 1‑2 And obtaining a quality factor value and a corresponding reflection coefficient inversion result by using the norm high-order derivative.
Description
Technical Field
The application relates to the field of geophysical exploration, in particular to a method based on L 1-2 Provided are a non-steady state blind deconvolution method and device for a norm high-order derivative criterion.
Background
Due to the incompletely elastic nature of the stratum, the seismic waves can generate phenomena of amplitude energy attenuation, phase distortion and band narrowing when propagating in the stratum. The quality factor Q is an important parameter describing the elastic properties of the formation and is used to quantitatively evaluate the degree of attenuation of seismic waves as they propagate underground. The Q value may reflect not only subsurface fracture information, but may also be used to predict hydrocarbon resources. The AVO effect can be better explained by applying the Q value, and the method can also be used for detecting the saturation degree of gas or fluid and the geometric distribution of phases. The closure of the fracture under different stresses significantly alters the Q of the rock, so that attenuating observations of core samples from a well can provide information about the fracture state of the reservoirThe material may alternatively provide a description of in situ stresses. Information such as lithology, pressure, permeability, etc. of the formation can also be predicted using the Q value. In contrast, the presence of formation attenuation makes subsurface imaging and interpretation more difficult, while reliable estimates of attenuation help to improve the resolution of seismic data. Reflection coefficient inversion, inverse Q filtering and the like are effective attenuation compensation methods and are widely developed and applied, but the processing effect of the attenuation compensation methods is always dependent on the accuracy of Q value estimation. In addition, amplitude and phase changes caused by the intrinsic attenuation of the stratum are embodied in wavelets, and the reflection coefficient is an inherent property of the stratum and does not change with the change of the Q value. According to the wave theory, when the wavelets are convolved with the corresponding reflection coefficient sets of the thin interbed, interference phenomenon occurs. The interference effect also causes variations in the amplitude and phase of the seismic waveform, resulting in an estimated attenuation being the Q value of the coupling of the intrinsic attenuation to the interference effect, which in turn can affect the subsequent inversion of the unsteady reflection coefficient. The Q value may be indicative of fluid information and the reflection coefficient may in turn be indicative of formation information. For unsteady data, the two are mutually influenced. In order to peel off the mutual influence effect of the two, the application is based on the sparse formation assumption and utilizes L 1-2 An nth derivative of the norm (e.g., n=2) quantifies the relationship between sparsity of the reflection coefficient inversion result and the Q value. The simultaneous processing of Q value estimation and deconvolution is performed by the strategy of scanning Q.
Traditional steady state reflection coefficient inversion methods, such as sparse layer inversion, base tracking inversion and the like, require that the input seismic data is steady, and in the actual field data acquisition process, the obtained seismic data is unsteady. Thus, if still based on steady state wavelet assumptions, the resulting seismic data processing results are inaccurate, which can affect subsequent seismic interpretation efforts. For unsteady seismic data, the conventional method is to perform attenuation compensation first and then perform steady-state reflection coefficient inversion, but the attenuation compensation needs to know the Q value, and the conventional Q value estimation method can be mainly divided into the following three categories: the first type of method mainly includes a time domain method and a frequency domain method. The second class of methods can be summarized as estimating the Q value using pre-stack CMP data. The third class of methods are some time-frequency domain methods. Q value estimation has been a hotspot studied by geophysicists, and an attenuation model thereof quantifies the process of energy attenuation and phase change when seismic waves propagate in a stratum, but further increases the discomfort of solving the inverse problem. When the seismic wave propagates underground, the precision of the quality factor obtained is reduced due to the factors such as scattering, spherical diffusion and the like. Taner and Treitel have listed 17 methods for calculating Q values, most of the methods consider attenuation in seismic waves, but the calculation accuracy and stability are limited by the pre-processing effect, the original quality, the rock properties and the like of seismic data, and the methods have respective application ranges and have no universality. In practice, the interference of the stratum also causes waveform changes of the seismic waves, and also affects the accuracy of the Q value estimation. The estimated Q value is accurate or not, so that the attenuation degree of the wavelet at different time positions is different from the actual attenuation degree, the accuracy of the time-varying wavelet matrix is affected, and finally, the reflection coefficient obtained by the unsteady deconvolution method is deviated.
Disclosure of Invention
The application aims to provide a method based on L 1-2 The unsteady state blind deconvolution method and device of the norm high-order derivative criterion are used for estimating the Q value of the quality factor through the deformation of the criterion function from the angle of reflection coefficient inversion and for unsteady state reflection coefficient inversion, and the resolution processing of post-stack seismic data is realized.
In order to achieve the above object, the unsteady state blind deconvolution method provided by the present application specifically comprises: acquiring seismic data to be analyzed, and presetting an initial seismic wavelet in the seismic data; constructing an attenuated sub-wave group through a preset quality factor range and the initial seismic wavelet, and inverting the reflection coefficient of the attenuated sub-wave group to obtain an inverted reflection coefficient; substituting the reflection coefficient into L 1-2 Obtaining a derivative curve by obtaining a second derivative after the norm function; minimizing the corresponding L in the derivative curve 1-2 And obtaining the quality factor value and the corresponding reflection coefficient inversion result by the second derivative of the norm function.
In the above unsteady state blind deconvolution method, preferably, the L 1-2 The norm function comprises:
in the above formula, 0 < sigma < 1 is a weight factor, N is a time sampling point number, r is an inverted reflection coefficient sequence, r i For the i-th element in r, Q is the quality factor, J (Q) is the sparsity criterion L 1-2 And a norm function value.
In the above-mentioned unsteady state blind deconvolution method, it is preferable that the reflection coefficient after inversion is substituted into L 1-2 Obtaining the derivative curve by obtaining the second derivative after the norm function comprises the following steps: substituting the inverted reflection coefficient into L 1-2 Obtaining a scanning curve by a norm function; obtaining a second derivative of the scanning curve to obtain a second-order curve; and respectively carrying out smoothing treatment on the scanning curve and the second-order curve, and then carrying out superposition treatment to obtain a derivative curve.
In the above unsteady state blind deconvolution method, preferably, inversion of the reflection coefficient is performed on the attenuator wave group, and obtaining the inverted reflection coefficient includes: and carrying out reflection coefficient inversion on the attenuated sub-wave groups by using a sparse Bayesian learning algorithm to obtain inverted reflection coefficients.
In the above-mentioned unsteady state blind deconvolution method, preferably, the attenuated wavelet group includes:
in the above formula, M is the number of frequency sampling points, F H For the anti-Fourier operator, W is Toeplitz matrix formed by the initial wavelet after Fourier transformation, H is conjugate transpose, s represents unsteady seismic record, A is time-frequency domain attenuation matrix related to quality factor, r is reflection coefficient sequence, and G (Q) is time-varying wavelet matrix.
The application also provides an unsteady state blind deconvolution device, which comprises an acquisition module and deconvolutionThe system comprises a modeling module, a computing module and an analyzing module; the acquisition module is used for acquiring the seismic data to be analyzed, and an initial seismic wavelet is preset in the seismic data; the inversion module is used for constructing an attenuated sub-wave group through a preset quality factor range and the initial seismic wavelet, and inverting the reflection coefficient of the attenuated sub-wave group to obtain an inverted reflection coefficient; the calculation module is used for substituting the reflection coefficient into L 1-2 Obtaining a derivative curve by obtaining a second derivative after the norm function; the analysis module is used for minimizing the corresponding L in the derivative curve 1-2 And obtaining the quality factor value and the corresponding reflection coefficient inversion result by the second derivative of the norm function.
In the above unsteady state blind deconvolution device, preferably, the calculation module further includes a preprocessing unit for substituting the inverted reflection coefficient into L 1-2 Obtaining a scanning curve by a norm function; obtaining a second derivative of the scanning curve to obtain a second-order curve; and respectively carrying out smoothing treatment on the scanning curve and the second-order curve, and then carrying out superposition treatment to obtain a derivative curve.
In the above-mentioned unsteady state blind deconvolution device, preferably, the inversion module includes: and carrying out reflection coefficient inversion on the attenuated sub-wave groups by using a sparse Bayesian learning algorithm to obtain inverted reflection coefficients.
The application 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 application also provides a computer readable storage medium storing a computer program for executing the above method.
The beneficial technical effects of the application are as follows: under the framework of a Bayesian learning method, a time-varying wavelet group is constructed by L 1-2 Quantizing the internal relation between the sparsity of the reflection coefficient inversion result and the Q value by using a norm high-order derivative criterion function, estimating the Q value or time-varying wavelet, and simultaneously completing the unsteady reflection coefficient inversion to obtainThe purpose of improving the resolution of the earthquake is to improve the resolution of the earthquake. 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 application and are incorporated in and constitute a part of this specification, illustrate and together with the description serve to explain the application. In the drawings:
FIG. 1 is a schematic flow chart of an unsteady state blind deconvolution method according to an embodiment of the present application;
FIG. 2 is a flowchart of a derivative curve calculation process according to an embodiment of the present application;
FIG. 3 is a schematic diagram illustrating an application flow of the unsteady-state blind deconvolution method according to an embodiment of the present application in actual operation;
FIG. 4 is a schematic flow chart of an unsteady-state blind deconvolution method according to an embodiment of the present application;
FIGS. 5A-5F are schematic diagrams illustrating a statistical numerical model of an unsteady state blind deconvolution method according to an embodiment of the present application;
fig. 6A to fig. 6D are schematic diagrams illustrating a test of an unsteady state blind deconvolution method applied to three-dimensional physical model data according to an embodiment of the present application;
fig. 7A to fig. 7D are schematic diagrams illustrating a test of an unsteady state blind deconvolution method applied to three-dimensional actual data according to an embodiment of the present application;
FIG. 8 is a schematic diagram of an unsteady-state blind deconvolution device according to an embodiment of the present application;
fig. 9 is a schematic structural diagram of an electronic device according to an embodiment of the application.
Detailed Description
The following will describe embodiments of the present application in detail with reference to the drawings and examples, thereby solving the technical problems by applying technical means to the present application, and realizing the technical effects can be fully understood and implemented accordingly. It should be noted that, as long as no conflict is formed, each embodiment of the present application and each feature of each embodiment may be combined with each other, and the formed technical solutions are all within the protection scope of the present application.
Additionally, the steps illustrated in the flowcharts 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 flowcharts, in some cases the steps illustrated or described may be performed in an order other than that herein.
Referring to fig. 1, the unsteady blind deconvolution method provided by the present application specifically includes: s101, acquiring seismic data to be analyzed, and presetting an initial seismic wavelet in the seismic data; s102, constructing an attenuated sub-wave group through a preset quality factor range and the initial seismic wavelet, and inverting the reflection coefficient of the attenuated sub-wave group to obtain an inverted reflection coefficient; s103 substituting the reflection coefficient into L 1-2 Obtaining a derivative curve by obtaining a second derivative after the norm function; s104, minimizing the corresponding L in the derivative curve 1-2 And obtaining the quality factor value and the corresponding reflection coefficient inversion result by the second derivative of the norm function. Referring to FIG. 2, in the step S103, the inverted reflection coefficient is substituted into L 1-2 Obtaining the derivative curve by obtaining the second derivative after the norm function further comprises: s201 substituting the inverted reflection coefficient into L 1-2 Obtaining a scanning curve by a norm function; s202, obtaining a second derivative of the scanning curve to obtain a second-order curve; s203, respectively carrying out smoothing treatment on the scanning curve and the second-order curve, and then carrying out superposition treatment to obtain a derivative curve. In step S102, inversion of the reflection coefficient of the attenuated sub-wave group is performed, and obtaining the inverted reflection coefficient may include: and carrying out reflection coefficient inversion on the attenuated sub-wave groups by using a sparse Bayesian learning algorithm to obtain inverted reflection coefficients.
In the above embodiment, the L 1-2 The norm function may be:
in the above formula, 0 < sigma < 1, N is the number of time sampling points, and r is the inverted reflectionCoefficient sequence, r i For the i-th element in r, Q is the quality factor, J (Q) is the sparsity criterion L 1-2 And a norm function value.
The attenuated wavelet group comprises:
in the above formula, M is the number of frequency sampling points, F H For the inverse Fourier operator, W is a diagonal matrix formed by the initial wavelets after Fourier transformation, H is a conjugate transpose, s represents an attenuated 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 attenuated wavelet matrix.
Specifically, the unsteady state blind deconvolution method provided by the application mainly comprises two parts in principle: unsteady seismic data forward modeling based on wave equation and L-based method 1-2 Q value estimation of a norm higher derivative criterion;
the forward modeling of unsteady seismic data based on wave equation mainly comprises the following steps:
for a plane wave propagating along the z-axis direction, the analytical expression of the one-dimensional wave equation is:
where x is the displacement function, ω is the angular frequency, t is time, i is the imaginary unit, and k is the wave number. When the seismic wave propagates in the underground medium, viscoelastic absorption attenuation and dispersion can occur; one way to describe this phenomenon is to express the wavenumber k as an actual wavenumber k r (ω) 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 conventional convolution model does not consider the attenuation characteristic of the stratum, so in one embodiment of the application, the final unsteady seismic data simulation formula can be deduced mainly according to the acoustic wave equation theory and the deformation of the formula (1) and the formula (2):
where s represents the decaying seismic record, α (ω, τ) z ,Q z ) An exponential function representing the variation with frequency and time, r (z) being the reflection coefficient of zero angle at depth z; further, the unsteady convolution model of the time domain may be described as:
wherein: m is the number of frequency samples to be taken,
is an inverse Fourier operator, where f m = (M-1) Δf (m=1, 2, … M) represents the sampling frequency of the mth point, M represents the frequency sampling point number; h represents a conjugate transpose;representing a diagonal matrix formed by the initial wavelets after Fourier transformation, wherein H is a conjugate transpose; />For a time-frequency domain attenuation matrix related to Q, the attenuation equation is:
when the Q value is infinite, alpha approaches to 1, namely the traditional steady-state deconvolution formula; g (Q) is an attenuated wavelet matrix. From the formula(4) As can be seen, the unsteady seismic record can be seen as a weighted linear superposition of N time-varying wavelets, each weighted with a weight r n . Because of the band-limited nature of the seismic data and the introduction of Q values, it is not straightforward to solve equation (4); therefore, the application is mainly based on stratum sparse hypothesis, adopts a Sparse Bayesian Learning (SBL) method, and obtains the reflection coefficient sequence through inversion.
With respect to L-based 1-2 The Q value estimation of the norm higher derivative criterion mainly comprises:
amplitude and phase changes caused by the intrinsic attenuation of the formation are manifested in the attenuated seismic wavelets, while the reflection coefficient does not change with the change in Q. In practice, when the wavelet is convolved with the corresponding reflection coefficient sequence of the thin interbed, an interference effect is generated, and the interference effect also causes the amplitude and phase of the seismic waveform to change, so that the estimated Q value is a pseudo value generated by coupling effect of the intrinsic attenuation of the stratum and the interference effect, which affects the subsequent inversion of the unsteady reflection coefficient. Meanwhile, on the basis of sparse assumption, correlation exists between the Q value and the sparsity of the reflection coefficient inversion result, when the Q value is not right, the inversion result and the real model are different, and the correlation is mainly represented on the pulse amplitude and the number of non-zero pulses; when inversion is carried out by using an attenuation wavelet matrix constructed by smaller Q values, deep wavelets are excessively attenuated. In order to match waveforms, the pulse amplitude of the deep reflection coefficient is correspondingly increased, and the number of pulses is increased; when inversion is carried out by using an attenuation wavelet matrix constructed by a larger Q value, wavelet attenuation is reduced, so that a plurality of homopolar small pulses appear; for this reason, for the above phenomena, the present application quantifies by a sparse criterion function, in particular:
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, use L 1-2 N-th derivative of norm (n is taken to be 2) to quantify the relationship between Q value and inverted reflection coefficient, L 1-2 The norm may be expressed as:
wherein: 0 < sigma < 1, N represents the number of time sampling points, r represents the inverted reflection coefficient sequence, r i Is the i-th element in r. However, L 1-2 The norms do not characterize the correct Q value. The application discovers that the convex points of the second derivative can correctly indicate the position and the size of the Q value. The expressions of the first derivative and the second derivative are respectively:
L 1-2 the second derivative of the norm has a good law of convex function with the change of the Q value. Therefore, the Q value is estimated by minimizing the formula (7), and the optimal reflection coefficient is inverted at the same time, so that the aim of decoupling inherent attenuation and interference attenuation is fulfilled, and finally, a high-precision Q field and an inversion result with high resolution are obtained.
In combination with the embodiment, the application provides an L-based 1-2 The application flow of the unsteady state blind deconvolution method of the norm high-order derivative criterion in actual work is specifically shown in fig. 3:
s301, inputting attenuated seismic data;
s302, estimating an initial seismic wavelet from the seismic data;
s303, constructing a time-varying sub-wave group by adopting a Q value scanning strategy. Inversion of reflection coefficients is carried out by using a sparse Bayesian learning method, and time-varying wavelets are eliminated;
s304 substituting the reflection coefficient obtained by inversion into L 1-2 And (5) a norm function to obtain a scanning curve. Obtaining 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, the scanning curve and the curve after the second derivative is obtained are respectively subjected to smoothing and superposition treatment;
s305 minimizing L 1-2 And obtaining an estimated Q value and a corresponding reflection coefficient inversion result by using the norm high-order derivative.
Referring to fig. 4, in the application process of the unsteady state blind deconvolution method provided by the application, initial wavelets are estimated from input attenuated seismic data, a range of scan Q is provided, and an attenuated sub-wave group is constructed by combining the initial wavelets and the attenuated sub-wave group for unsteady state reflection coefficient inversion; secondly, substituting the reflection coefficient obtained by inversion into a criterion function, and solving a second derivative; and then, minimizing the second derivative to obtain the required Q value, and finally obtaining a high-precision Q field and a high-resolution inversion result. Based on the flow, the application tests by respectively using a statistical numerical model, three-dimensional physical model data and three-dimensional actual data; specifically, through 1000 geological models designed in advance, each model has 30 nonzero random reflection coefficient pulses, the maximum interval between nonzero pulses is not more than 30ms, for practical reasons, the nonzero pulse amplitude is set between-0.3 and 0.3, and the total length of the model is 1s. The method comprises the steps of selecting a Rake wavelet with a main frequency of 30Hz as an initial wavelet, enabling a Q value to be 50, constructing an attenuated seismic wavelet matrix, and respectively carrying out convolution on the attenuated seismic wavelet matrix and 1000 reflection coefficient models to obtain 1000 unsteady seismic data. Referring to fig. 5A to 5F, fig. 5A is a graph comparing the raw data and the inversion result, and fig. 5B is a graph comparing the raw reflectance and the inversion reflectance. As can be seen from fig. 5A, the seismic data synthesized by the inversion result matches very well with the raw data. On the premise of data matching, the reflection coefficient obtained by inversion is basically identical with the original reflection coefficient, as shown in FIG. 5B; FIG. 5C is a block diagram of the substitution of the inverted reflectance results for each corresponding Q value into L 1-2 The 1000 scan curves obtained in the norm criterion function are shown in fig. 5D, which is a superposition of the 1000 scan curves. As can be seen from fig. 5D, the superimposed curve does not have significant convexity nor indicates the location of the correct Q value; fig. 5E is a graph obtained by performing second derivative on the scan curve of fig. 5C, and fig. 5F is a graph obtained by superimposing fig. 5E. As is apparent from FIG. 5F, by finding L 1-2 The second derivative of the norm criterion function has obvious convexity, and the inflection point position exactly corresponds to the position of the true Q value. By statisticsThe test of the sexual value shows that the Q value can be accurately estimated, and the Q value is not easily influenced by the reflection coefficient model, so that the method has high applicability and accuracy.
Referring to fig. 6A to 6D, the unsteady state blind deconvolution method is applied to the three-dimensional physical model data for testing. FIG. 6A is a graph of the result of the inverse reflection coefficient for each corresponding Q value substituted into L 1-2 In the norm criterion function, a superimposed curve; FIG. 6B is a graph of the second derivative of FIG. 6A; FIG. 6C is raw three-dimensional physical model data; fig. 6D shows the result of the high resolution processing according to the present application. Comparing fig. 6C and 6D, the processed results have a clearer layer interface and show more construction details, and are more 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 three-dimensional actual data, and fig. 7A is a graph of substituting the inversion reflection coefficient result of each corresponding Q value into L 1-2 In the norm criterion function, a superimposed curve; FIG. 7B is a graph of the second derivative of FIG. 7A; FIG. 7C is a slice of raw real data along layer amplitude; fig. 7D shows the high resolution processing result, and by comparing the original seismic record with the high resolution processing result, the processed river edge and river tip are clearer, and the seismic resolution is obviously improved.
Referring to fig. 8, the present application 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 the seismic data to be analyzed, and an initial seismic wavelet is preset in the seismic data; the inversion module is used for constructing an attenuated sub-wave group through a preset quality factor range and the initial seismic wavelet, and inverting the reflection coefficient of the attenuated sub-wave group to obtain an inverted reflection coefficient; the calculation module is used for substituting the reflection coefficient into L 1-2 Obtaining a derivative curve by obtaining a second derivative after the norm function; the analysis module is used for minimizing the corresponding L in the derivative curve 1-2 Obtaining quality factor taking by second derivative of norm functionAnd inverting the value and the corresponding reflection coefficient.
In the above unsteady state blind deconvolution device, the calculation module may further include a preprocessing unit for substituting the inverted reflection coefficient into L 1-2 Obtaining a scanning curve by a norm function; obtaining a second derivative of the scanning curve to obtain a second-order curve; and respectively carrying out smoothing treatment on the scanning curve and the second-order curve, and then carrying out superposition treatment to obtain a derivative curve. The inversion module may include: and carrying out reflection coefficient inversion on the attenuated sub-wave groups by using a sparse Bayesian learning algorithm to obtain inverted reflection coefficients.
The application 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 application also provides a computer readable storage medium storing a computer program for executing the above method.
The beneficial technical effects of the application are as follows: under the framework of a Bayesian learning method, a time-varying wavelet group is constructed by L 1-2 The norm high-order derivative criterion function quantifies the internal relation between the sparsity of the reflection coefficient inversion result and the Q value, estimates the Q value or time-varying wavelet, and simultaneously completes the inversion of the unsteady reflection coefficient, 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 device 600 according to an embodiment of the present application. As shown in fig. 9, the electronic device 600 may include a central processor 100 and a memory 140; memory 140 is coupled to central processor 100. Notably, the diagram is exemplary; other types of structures may also be used in addition to or in place of the structures to implement telecommunications functions or other functions.
In one embodiment, the CT image shading correction function may be integrated into the central processor 100. Wherein the central processor 100 may be configured to control as follows: performing image texture removal operation on the original CT reconstructed image to obtain a smooth image; dividing the structural components of the original CT reconstructed image according to human tissues to construct a template image; and performing shading correction according to the smooth image and the template image.
Performing image texture removal operation on the original CT reconstructed image to obtain a smooth image, wherein the image texture removal operation 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 for constructing the template image comprises the following steps of: dividing the original CT reconstruction image into a plurality of human tissue areas by adopting an image dividing method; and filling CT values of corresponding tissues under the voltages of the X-ray tube in different human tissue areas respectively to obtain the template image.
Wherein performing shading correction according to the smoothed image and the template image includes: the smooth image and the template image are subjected to difference to obtain a residual image, wherein the residual image comprises image shadows and tissue structure errors; performing low-pass filtering treatment on the tissue structure errors 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.
The low-pass filtering processing is performed on the tissue structure error of the residual image, and the low-pass filtering processing comprises the following steps: and carrying out low-pass filtering processing on the tissue structure errors of the residual image by using 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:
the texture-free smooth image is calculated by an objective function, which is as follows:
wherein, the liquid crystal display device comprises a liquid crystal display device,
S p index of the p-th pixel of the non-texture smooth image S, I p Index the p-th pixel of the original CT reconstructed image I,the number of pixel indices p, λ is a smoothing factor, ">Is->The partial derivatives in both the x and y directions are indexed for the pixel.
In another embodiment, the CT image shading correction device may be configured separately from the central processor 100, for example, the CT image shading correction device may be configured as a chip connected to the central processor 100, and the CT image shading correction function is implemented by control of the central processor.
As shown in fig. 9, the electronic device 600 may further include: a communication module 110, an input unit 120, an audio processing unit 130, a display 160, a power supply 170. It is noted that the electronic device 600 need not include all of the components shown in fig. 9; in addition, the electronic device 600 may further include components not shown in fig. 9, to which reference is made to the related art.
As shown in fig. 9, the central processor 100, sometimes also referred to as a controller or operational control, may include a microprocessor or other processor device and/or logic device, which central processor 100 receives inputs and controls 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 about failure may be stored, and a program for executing the information may be stored. And the central processor 100 can execute the program stored in the memory 140 to realize information storage or processing, etc.
The input unit 120 provides an input to the central processor 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 for displaying display objects such as images and characters. The display may be, for example, but not limited to, an LCD display.
The memory 140 may be a solid state memory such as Read Only Memory (ROM), random Access Memory (RAM), SIM card, or the like. But also a memory which holds information even when powered down, can be selectively erased and provided with further data, an example of which is sometimes referred to as EPROM or the like. Memory 140 may also be some other type of device. Memory 140 includes a buffer memory 141 (sometimes referred to as a buffer). The memory 140 may include an application/function storage 142, the application/function storage 142 for storing application programs and function programs or a flow for executing operations of the electronic device 600 by the central processor 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 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 applications, address book applications, etc.).
The communication module 110 is a transmitter/receiver 110 that transmits and receives signals via an antenna 111. A 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, etc., 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 to receive audio input from the microphone 132 to implement usual telecommunication functions. The audio processor 130 may include any suitable buffers, decoders, amplifiers and so forth. In addition, the audio processor 130 is also coupled to the central processor 100 so that sound can be recorded locally through the microphone 132 and so that sound stored locally can be played through the speaker 131.
It will be appreciated by those skilled in the art that embodiments of the present application may be provided as a method, system, or computer program product. Accordingly, the present application may take the form of an entirely hardware embodiment, an entirely software embodiment or an embodiment combining software and hardware aspects. Furthermore, the present application 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 application is described with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems) and computer program products according to embodiments of the application. It will be understood that each flow and/or block of the flowchart illustrations and/or block diagrams, and combinations of flows and/or blocks in the flowchart illustrations 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 foregoing description of the embodiments has been provided for the purpose of illustrating the general principles of the application, and is not meant to limit the scope of the application, but to limit the application to the particular embodiments, and any modifications, equivalents, improvements, etc. that fall within the spirit and principles of the application are intended to be included within the scope of the application.
Claims (4)
1. A method of unsteady state blind deconvolution, said method comprising:
acquiring seismic data to be analyzed, and presetting an initial seismic wavelet in the seismic data;
constructing an attenuated sub-wave group through a preset quality factor range and the initial seismic wavelet, and inverting the reflection coefficient of the attenuated sub-wave group to obtain an inverted reflection coefficient;
substituting the reflection coefficient into L 1-2 Obtaining a derivative curve by obtaining a second derivative after the norm function;
minimizing L corresponding to the derivative curve 1-2 Obtaining the quality factor value and the corresponding reflection coefficient inversion result by the second derivative of the norm function;
the L is 1-2 The norm function comprises:
；
in the above formula, 0 < sigma < 1 is a weight factor, N is a time sampling point number, r is an inverted reflection coefficient sequence, r i For the i-th element in r, Q is the quality factor, J (Q) is the sparsity criterion L 1-2 A norm function value;
substituting the inverted reflection coefficient into L 1-2 Obtaining the derivative curve by obtaining the second derivative after the norm function comprises the following steps:
substituting the inverted reflection coefficient into L 1-2 Obtaining a scanning curve by a norm function;
obtaining a second derivative of the scanning curve to obtain a corresponding second-order curve;
respectively carrying out smoothing treatment on the scanning curve and the second-order curve, and then carrying out superposition treatment to obtain a derivative curve;
inversion of reflection coefficients is carried out on the attenuated sub-wave groups, and the obtained inversion of the reflection coefficients comprises the following steps: performing reflection coefficient inversion on the attenuated sub-wave groups by using a sparse Bayesian learning algorithm to obtain inverted reflection coefficients;
the attenuated wavelet group comprises:
；
in the above formula, M is the number of frequency sampling points, F H For the inverse Fourier operator, W is Toeplitz matrix formed by the initial wavelet after Fourier transformation, H is conjugate transpose, s represents unsteady seismic record, A is time-frequency domain attenuation matrix related to quality factor, r is inverse reflection coefficient sequence, and G (Q) is time-varying wavelet matrix.
2. The unsteady state 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 the seismic data to be analyzed, and an initial seismic wavelet is preset in the seismic data;
the inversion module is used for constructing an attenuated sub-wave group through a preset quality factor range and the initial seismic wavelet, and inverting the reflection coefficient of the attenuated sub-wave group to obtain an inverted reflection coefficient;
the calculation module is used for substituting the reflection coefficient into L 1-2 Obtaining a derivative curve by obtaining a second derivative after the norm function;
the analysis module is used for minimizing the corresponding L in the derivative curve 1-2 Obtaining the quality factor value and the corresponding reflection coefficient inversion result by the second derivative of the norm function;
the L is 1-2 The norm function comprises:
；
in the above formula, 0 < sigma < 1 is a weight factor, N is a time sampling point number, r is an inverted reflection coefficient sequence, r i For the i-th element in r, Q is the quality factor, J (Q) is the sparsity criterion L 1-2 A norm function value;
the calculation module also comprises a preprocessing unit for substituting the inverted reflection coefficient into L 1-2 Obtaining a scanning curve by a norm function; obtaining a second derivative of the scanning curve to obtain a second-order curve; respectively carrying out smoothing treatment on the scanning curve and the second-order curve, and then carrying out superposition treatment to obtain a derivative curve;
the inversion module includes: performing reflection coefficient inversion on the attenuated sub-wave groups by using a sparse Bayesian learning algorithm to obtain inverted reflection coefficients;
the attenuated wavelet group comprises:
；
in the above formula, M is the number of frequency sampling points, F H For the inverse Fourier operator, W is Toeplitz matrix formed by the initial wavelet after Fourier transformation, H is conjugate transpose, s represents unsteady seismic record, A is time-frequency domain attenuation matrix related to quality factor, r is inverse reflection coefficient sequence, and G (Q) is time-varying wavelet matrix.
3. 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 unsteady state blind deconvolution method of claim 1 when the computer program is executed by the processor.
4. A computer readable storage medium, wherein the computer readable storage medium stores a computer program for performing the unsteady state blind deconvolution method of claim 1.
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 CN113219530A (en) | 2021-08-06 |
| CN113219530B true 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范数极小化的稀疏信号重建条件;周珺 等;《合肥工业大学学报（自然科学版）》;第43卷(第1期);第137-140页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113219530A (en) | 2021-08-06 |
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 | |
| Komatitsch et al. | Anelastic sensitivity kernels with parsimonious storage for adjoint tomography and full waveform inversion | |
| RU2422856C2 (en) | Recording method of seismic pulse stretching in seismic data | |
| US7082368B2 (en) | Seismic event correlation and Vp-Vs estimation | |
| CN103293552B (en) | A kind of inversion method of Prestack seismic data and system | |
| Kamath et al. | Elastic full-waveform inversion for VTI media: A synthetic parameterization study | |
| Shragge et al. | Full-wavefield modeling and reverse time migration of laser ultrasound data: A feasibility study | |
| Lefeuve-Mesgouez et al. | Retrieval of the physical properties of an anelastic solid half space from seismic data | |
| CN103018773A (en) | Method for processing coloured compensation of seismic record driven by logging signal | |
| Tejero et al. | Appraisal of instantaneous phase-based functions in adjoint waveform inversion | |
| CN113219530B (en) | Unsteady state blind deconvolution method and device | |
| Nadri et al. | Estimation of stress-dependent anisotropy from P-wave measurements on a spherical sample | |
| CN113671572B (en) | Seismic data imaging method and device based on indoor sand box | |
| Wei et al. | A high-resolution microseismic source location method based on contrast source algorithm | |
| Chen et al. | Near-Borehole Formation Acoustic Logging Imaging: A Full Waveform Inversion Algorithm in Cylindrical Coordinates | |
| Zyserman et al. | Analysis of the numerical dispersion of waves in saturated poroelastic media | |
| CN112363244B (en) | Wave impedance inversion method and carbonate heterogeneous reservoir prediction method and system | |
| CN113009579B (en) | Seismic data inversion method and device | |
| Li et al. | Influence of upscaling on identification of reservoir fluid properties using seismic-scale elastic constants | |
| CN111736220B (en) | Reverse time migration imaging method and device | |
| Mo et al. | The development and testing of a 2D laboratory seismic modelling system for heterogeneous structure investigations | |
| US20220390631A1 (en) | Laplace-fourier 1.5d forward modeling using an adaptive sampling technique | |
| Cui | Marchenko focusing for target-oriented data processing and full-waveform inversion | |
| Han et al. | Fractional-order velocity-dilatation-rotation viscoelastic wave equation and numerical solution based on a constant-Q model |
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 |