CN113552631A - Time-frequency double-domain regularization sparse deconvolution method and device for narrow-band signals - Google Patents

Time-frequency double-domain regularization sparse deconvolution method and device for narrow-band signals Download PDF

Info

Publication number
CN113552631A
CN113552631A CN202110945575.0A CN202110945575A CN113552631A CN 113552631 A CN113552631 A CN 113552631A CN 202110945575 A CN202110945575 A CN 202110945575A CN 113552631 A CN113552631 A CN 113552631A
Authority
CN
China
Prior art keywords
wavelet
gradient
reflection coefficient
sparse
domain
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
Application number
CN202110945575.0A
Other languages
Chinese (zh)
Other versions
CN113552631B (en
Inventor
金丹
王保利
王云宏
豆旭谦
覃思
代晨昱
裴跟弟
张庆庆
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Research Institute Co Ltd of CCTEG
Original Assignee
Xian Research Institute Co Ltd of CCTEG
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Xian Research Institute Co Ltd of CCTEG filed Critical Xian Research Institute Co Ltd of CCTEG
Priority to CN202110945575.0A priority Critical patent/CN113552631B/en
Publication of CN113552631A publication Critical patent/CN113552631A/en
Application granted granted Critical
Publication of CN113552631B publication Critical patent/CN113552631B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/63—Seismic attributes, e.g. amplitude, polarity, instant phase
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
    • Y02D30/00—Reducing energy consumption in communication networks
    • Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks

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)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

The invention relates to a time-frequency double-domain regularization sparse deconvolution method and a time-frequency double-domain regularization sparse deconvolution device for narrow-band signals. The method and the device convert the conventional reflection coefficient sparse regularization into the reflection coefficient gradient sparse regularization, the stagnation point of the gradient is output as sparse regularization during calculation, the stagnation point in the reflection coefficient curve is reserved through a threshold function, then the stagnation point value larger than a given threshold value is selected through a composite threshold function, and the given threshold value is calculated by the gradient self-adaption of input data. The method gradually iterates from a larger gradient to a smaller gradient of the weak signal by iterative reservation of the gradient of the reflection coefficient, avoids loss of the weak signal, and has better adaptivity and noise immunity. In addition, aiming at the characteristics of narrow-band signals, the wavelets are considered to be narrow-band, so that the seismic wavelets have sparse characteristics in a frequency domain, and therefore, the wavelet gradient is sparsely constrained in the frequency domain, and the method can be suitable for the narrow-band signals.

Description

Time-frequency double-domain regularization sparse deconvolution method and device for narrow-band signals
Technical Field
The invention relates to a signal processing method and a signal processing device, belongs to the technical field of geological detection, and particularly relates to a time-frequency double-domain regularization sparse deconvolution method and a time-frequency double-domain regularization sparse deconvolution device for narrow-band signals.
Background
In the seismic prospecting process, the seismic signals are generally required to be filtered to select effective signals, so that the time resolution is reduced by filtering, and the spatial resolution after offset is also reduced. In order to solve the problem of resolution reduction, the wavelet is compressed mainly by deconvolution in seismic data processing, so that the resolution is improved. Limited by an effective frequency band, the resolution of a filtered narrow-band signal is difficult to be greatly improved by a traditional deconvolution method, meanwhile, deconvolution is a ill-conditioned anti-problem due to the existence of noise in actual recording, and a pre-whitening method is usually adopted to add a small disturbance amount to a wavelet matrix to improve the stability of an algorithm, but the obtained resolution is not high enough. In addition, the conventional deconvolution is implemented based on the assumption that wavelets are minimum phase and reflection coefficients are white noise sequences, but the actual data does not conform to such characteristics, so that the deconvolution effect is not ideal.
Many researches have been made to improve the deconvolution effect by using other statistical characteristics to relax the constraints on the reflection coefficients, for example, in chinese patent 201610190426.7, the reflection coefficient estimation method and apparatus based on the bayesian inversion frame are used to perform sparse blind deconvolution under the bayesian inversion frame, and combine with the geostatistics, thereby reducing the unreasonable constraints on wavelets and reflection coefficients and improving the deconvolution accuracy. However, the reflection coefficient obtained in this patent does not satisfy the sparsity. Chinese patent 201710572136.3 discloses a simultaneous inversion method for a multi-channel seismic record reflection coefficient sequence, which fully considers the sparse characteristic of the reflection coefficient and improves the lateral continuity by utilizing multi-channel joint inversion. However, the method disclosed by the patent is complex, the related zero norm approximation parameter is difficult to control, and the deconvolution effect is influenced by too large and too small parameters, so that self-adaptation cannot be achieved.
The sparse deconvolution assumes that the reflection coefficient is a series of sparse pulses, extracts the amplitude and the position of the reflection coefficient from the noisy record through error least square fitting, suppresses noise while solving, and improves the resolution of the result compared with the traditional deconvolution. Based on the concept of sparse reflection coefficients, many different algorithms appear, wherein sparse deconvolution based on L0 norm is a more ideal sparse optimization algorithm. Direct solution based on the L0 norm minimization algorithm is difficult, and an approximate solution method is generally adopted. Among various approximate solving methods, the iterative hard threshold algorithm with a simpler solving form is more suitable for large-scale seismic data processing, so that the method is widely applied. The iteration hard threshold algorithm needs to select a proper initial threshold in the calculation process, the initial threshold is adjusted when the data are different, and the parameter adjustment needs to be carried out by investing more workload when a large amount of seismic data are faced, so that the self-adaptability is poor.
In summary, the conventional sparse deconvolution technology mainly has the following problems:
1. when the reflection coefficient sparse feature adopts a non-zero norm or an approximate zero norm, the sparse feature cannot be satisfied;
2. the zero norm sparse deconvolution technology can fully meet the coefficient property, but threshold parameters and the like introduced in the solving process cannot be adaptively adjusted according to data, so that the deconvolution convergence speed and precision are influenced;
3. the characteristics of the narrow-band signal are not considered, so that the sparse deconvolution effect of the narrow-band signal is poor, and the stability problem exists.
Disclosure of Invention
The following presents a simplified summary of one or more aspects in order to provide a basic understanding of such aspects. This summary is not an extensive overview of all contemplated aspects, and is intended to neither identify key or critical elements of all aspects nor delineate the scope of any or all aspects. Its sole purpose is to present some concepts of one or more aspects in a simplified form as a prelude to the more detailed description that is presented later.
The present invention is to solve the above problems in the prior art, and provide a deconvolution method and apparatus based on gradient sparse regularization for narrowband signals. The method and the device convert the conventional reflection coefficient sparse regularization into the reflection coefficient gradient sparse regularization, the stagnation point of the gradient is output as sparse regularization during calculation, the stagnation point in the reflection coefficient curve is reserved through a threshold function, then the stagnation point value larger than a given threshold value is selected through a composite threshold function, and the given threshold value is calculated by the gradient self-adaption of input data. The method gradually iterates from a larger gradient to a smaller gradient of the weak signal by iterative reservation of the gradient of the reflection coefficient, avoids loss of the weak signal, and has better adaptivity and noise immunity. In addition, aiming at the characteristics of narrow-band signals, the wavelets are considered to be narrow-band, so that the seismic wavelets have sparse characteristics in a frequency domain, and therefore, the wavelet gradient is sparsely constrained in the frequency domain, and the method can be suitable for the narrow-band signals.
In order to solve the problems, the scheme of the invention is as follows:
a time-frequency double-domain regularization sparse deconvolution method for narrow-band signals, comprising:
calculating the reflection coefficient gradient of the input seismic data;
calculating to obtain a conjugate gradient;
sparse regularization is carried out on the reflection coefficient conjugate gradient in a time domain;
iteratively updating the reflection coefficient based on the step length of the conjugate gradient of the reflection coefficient obtained by linear search;
the wavelets of the seismic data are updated based on the updated reflection coefficients.
Preferably, in the above time-frequency two-domain regularization sparse deconvolution method for narrow-band signals, the reflection coefficient gradient Δ r is calculated based on the following formula:
Δr＝ωT(y-ωr)
in the formula, r is a reflection coefficient, omega is a seismic wavelet, y is observed seismic data, and T is a transpose of a wavelet matrix.
Preferably, in the above time-frequency two-domain regularization sparse deconvolution method for narrow-band signals, the nonlinear conjugate gradient method uses a linear combination of a currently calculated negative gradient direction and a previous conjugate gradient direction as a search direction of a next iteration, that is:
Figure BDA0003213311170000041
in the formula, Δ γkRepresenting the current conjugate gradient, Δ γk-1Denotes the last conjugate gradient, Δ rkRepresenting the current k-th reflection coefficient gradient, betakExpressing the linear weighting coefficient, and calculating the formula as follows:
Figure BDA0003213311170000042
preferably, the aforementioned time-frequency two-domain regularization sparse deconvolution method for narrow-band signals, where the sparse regularization of the reflection coefficient conjugate gradient in the time domain includes:
Δγ＝Hλ(Δγ)
Hλfor an improved threshold function, the values are:
Hλ＝H1(H2)
Figure BDA0003213311170000043
Figure BDA0003213311170000044
in the formula, theta is a function independent variable, and when the function independent variable is substituted into the reflection coefficient conjugate gradient, theta is the reflection coefficient conjugate gradient. β ═ α κ + (1- α) γ, κ is a nonzero mean of | θ | γ, γ is a maximum of | θ |, α gives a relative threshold coefficient, α ∈ [0,1 ∈ γ]，
Figure BDA0003213311170000045
Is the derivative of theta with respect to time t.
Preferably, the above time-frequency two-domain regularization sparse deconvolution method for narrow-band signals calculates the reflection coefficient iteration step size μ based on the following formula:
Figure BDA0003213311170000051
in the formula, ω is the seismic wavelet and Δ γ is the conjugate gradient of the reflection coefficient.
Preferably, the above time-frequency two-domain regularization sparse deconvolution method for narrow-band signals updates the reflection coefficient r based on the following formula:
r＝r+μΔγ
in the formula, r is a reflection coefficient, mu is an iteration step length, and delta gamma is a conjugate gradient of the reflection coefficient.
Preferably, the above-mentioned time-frequency two-domain regularization sparse deconvolution method for narrow-band signals, the updating wavelet initial values of seismic data based on the updated reflection coefficients includes:
calculating the gradient of the wavelet;
calculating the conjugate gradient of the wavelet;
sparse regularization is carried out on the wavelet conjugate gradient in a frequency domain, and then the wavelet conjugate gradient after frequency domain sparse is obtained by utilizing inverse Fourier transform;
calculating a wavelet iteration step length based on the thinned wavelet conjugate gradient;
updating the wavelet initial value based on the wavelet iteration step size.
Preferably, in the foregoing time-frequency two-domain regularization sparse deconvolution method for narrow-band signals, the wavelet gradient is calculated based on the following formula:
Δω＝(y-ωr)rT
in the formula, reflection coefficients r and y are observed seismic data, omega is a seismic wavelet, r is a reflection coefficient, and T is a transpose of a wavelet matrix.
Preferably, in the foregoing time-frequency two-domain regularization sparse deconvolution method for narrow-band signals, the computation of the wavelet conjugate gradient is based on the following formula:
Figure BDA0003213311170000061
in the formula (I), the compound is shown in the specification,
Figure BDA0003213311170000062
which is indicative of the current conjugate gradient,
Figure BDA0003213311170000063
represents the last conjugate gradient, Δ ωkRepresenting the current k-th sub-wave gradient, σkExpressing the linear weighting coefficient, and calculating the formula as follows:
Figure BDA0003213311170000064
preferably, the above time-frequency two-domain regularization sparse deconvolution method for narrow-band signals, where calculating a wavelet iteration step based on a sparse wavelet gradient includes:
Figure BDA0003213311170000065
wherein r is a reflection coefficient,
Figure BDA0003213311170000066
is the conjugate gradient of the seismic wavelet;
and/or:
updating the seismic wavelet based on:
Figure BDA0003213311170000067
where, ω is the seismic wavelet,
Figure BDA0003213311170000068
is a seismic waveletλ is the iteration step of the wavelet.
Therefore, the beneficial effects of the invention are as follows: in reserving L0While the norm-constrained sparse deconvolution has the advantage, sparse regularization of the reflection coefficients is changed to sparse regularization of the reflection coefficient gradients. By iterative preservation of the reflection coefficient gradient, the weak signal is gradually iterated from a larger gradient to a smaller gradient of the weak signal, and the loss of the weak signal is avoided. More importantly, the method avoids manual selection of threshold parameters, has great adaptivity, improves convergence rate, increases calculation efficiency, and has better practicability compared with other sparse deconvolution methods.
Drawings
The accompanying drawings, which are incorporated herein and form a part of the specification, illustrate embodiments of the present invention and, together with the description, further serve to explain the principles of the invention and to enable a person skilled in the pertinent art to make and use the disclosure.
FIG. 1 is a technical flow chart of the present invention;
FIG. 2 is a time-frequency two-domain regularization sparse deconvolution processed result (lower) of a synthetic seismic record (middle) convolved with a reflection coefficient sequence (upper);
FIG. 3 shows the result of a time-frequency two-domain regularization sparse deconvolution process on a synthetic seismic record (top) with noise added (bottom);
FIG. 4 is a comparison of the effect of the deconvolution method on actual data.
FIG. 5 is a comparison of time-frequency two-domain regularization sparse deconvolution method and L0 norm-constrained fast iteration hard threshold deconvolution method iteration convergence curves
Embodiments of the present invention will be described with reference to the accompanying drawings.
Detailed Description
Examples
As shown in fig. 1, a time-frequency two-domain regularization sparse deconvolution method for narrowband signals proposed in this embodiment.
Considering that the linear mathematical operation of the sparse signal is still the sparse signal, the gradient of the sparse reflection coefficient has sparsity under the condition that the initial value is sparse. For this purpose, the present embodiment converts the conventional sparse regularization of the reflection coefficient into sparse regularization of the reflection coefficient gradient. During calculation, the stagnation point of the gradient is output as sparsification, meanwhile, a threshold function is improved, the stagnation point in the reflection coefficient curve is reserved through the threshold function, then the stagnation point value larger than a given threshold value is selected through the composite threshold function, and the given threshold value is calculated in a self-adaptive mode through the gradient of input data. By iterative reservation of the gradient of the reflection coefficient, the gradient is gradually iterated from a larger gradient to a smaller gradient of the weak signal, so that the loss of the weak signal is avoided, and the method has better self-adaptability and noise immunity. In addition, the wavelet of the narrow-band signal is also narrow-band, so that the seismic wavelet has sparse characteristics in a frequency domain, and therefore sparse constraint is carried out on the wavelet gradient in the frequency domain, and the embodiment can adapt to the narrow-band signal.
The time-frequency double-domain regularization sparse deconvolution method of the narrow-band signal of the embodiment specifically comprises the following steps:
calculating the reflection coefficient gradient of the input seismic data, and performing sparse regularization on the reflection coefficient gradient in a time domain;
iteratively updating the reflection coefficient based on the reflection coefficient gradient step length obtained by linear search;
and updating wavelet initial values of the seismic data based on the updated reflection coefficients.
The steps are further described in detail below with reference to FIGS. 1-3.
Step 1, reading in seismic data;
step 2, inputting a path of seismic data, giving an initial value of a reflection coefficient r as a zero vector, and giving an initial value of a seismic wavelet omega as a pulse vector;
step 3, calculating the gradient of the reflection coefficient and the conjugate gradient thereof according to the following formula:
Δr＝ωT(y-ωr)
in the formula, Δ r is a reflection coefficient gradient, r is a reflection coefficient, ω is a seismic wavelet, y is observed seismic data, and T is a transpose of a wavelet matrix.
The conjugate gradient was calculated from the reflection coefficient gradient:
Figure BDA0003213311170000081
in the formula, Δ γkRepresenting the current conjugate gradient, Δ γk-1Denotes the last conjugate gradient, Δ rkRepresenting the current k-th reflection coefficient gradient, betakExpressing the linear weighting coefficient, and calculating the formula as follows:
Figure BDA0003213311170000091
step 4, sparse regularization is carried out on the reflection coefficient conjugate gradient in the time domain:
Δγ＝Hλ(Δγ)
wherein HλFor improved threshold function, first pass threshold function H2Retaining the stagnation point in the reflection coefficient curve, i.e. the position of wave crest and wave trough, and then using H1Selecting greater than Δ r ═ Hλ(Δ r) the stagnation value of the threshold β, whereas the given threshold is calculated adaptively from the gradient of the input data:
Hλ＝H1(H2)
Figure BDA0003213311170000092
Figure BDA0003213311170000093
in the formula, θ is a function independent variable, and when we perform sparse regularization on the reflection coefficient conjugate gradient, θ is the reflection coefficient conjugate gradient. β ═ α κ + (1- α) γ, κ is a nonzero mean of | θ | γ, γ is a maximum of | θ |, α gives a relative threshold coefficient, α ∈ [0,1 ∈ γ]，
Figure BDA0003213311170000094
Is the derivative of theta with respect to time t.
And 5, calculating the iterative step size mu of the reflection coefficient by using a linear search method according to the following formula:
Figure BDA0003213311170000095
wherein | | | purple hair2Representing a two-norm.
And 6, iteratively updating the reflection coefficient r (when the reflection coefficient r is iterated for the first time r ═ r0)：
r＝r+μΔγ
That is, the sparse regularization of the reflection coefficient is changed into the sparse regularization of the reflection coefficient gradient, so that the absolute filtering of the reflection coefficient by the threshold function is converted into relative filtering. By iterative preservation of the reflection coefficient gradient, the weak signal is gradually iterated from a larger gradient to a smaller gradient of the weak signal, and the loss of the weak signal is avoided.
The sparse deconvolution comprises two steps of reflection coefficient inversion and wavelet inversion, wherein the step 1-6 is the iterative inversion of the reflection coefficient, and the wavelet progressive iterative inversion is introduced by combining the step 7-12;
step 7, calculating the wavelet gradient by using the following formula:
Δω＝(y-ωr)rT
step 8, calculating the wavelet conjugate gradient by using the following formula:
Figure BDA0003213311170000101
in the formula (I), the compound is shown in the specification,
Figure BDA0003213311170000102
which is indicative of the current conjugate gradient,
Figure BDA0003213311170000103
represents the last conjugate gradient, Δ ωkRepresenting the current k-th sub-wave gradient, σkExpressing the linear weighting coefficient, and calculating the formula as follows:
Figure BDA0003213311170000104
step 9, after the reflection coefficient updating is completed, utilizing H in the frequency domain1The wavelet conjugate gradient is subjected to sparse regularization by a domain value function, and then inverse Fourier transform is carried out to obtain the wavelet conjugate gradient after frequency domain sparseness:
Figure BDA0003213311170000105
in the formula, FFT is Fourier transform, IFFT is inverse Fourier transform.
Step 10, calculating wavelet iteration step size lambda according to the following formula
Figure BDA0003213311170000106
Step 11, updating the seismic wavelet with the following formula (ω ═ ω at the first iteration)0)：
Figure BDA0003213311170000107
Step 12, repeating the steps 3-11 until a more ideal deconvolution result is obtained through iteration;
and step 13, returning to the step 2, calculating the next channel until all channels are calculated, and outputting the calculated reflection coefficient and each channel of wavelet.
As shown in fig. 2 and 3, a sparsely distributed reflection coefficient sequence is constructed, and a rake wavelet with a main frequency of 100Hz is selected, and convolution is performed to obtain a corresponding synthetic seismic record (in fig. 2), wherein the sampling interval is 2 ms. The seismic record with random noise is shown in fig. 3 (middle), and the above deconvolution method based on gradient sparse regularization is applied to perform calculation, and the obtained results are shown in fig. 2 (lower) and fig. 3 (lower). It can be seen that the synthetic seismic record without noise can completely restore the initial reflection coefficient sequence after deconvolution computation. Noise-added synthetic seismic records, in which weak signals are mostly annihilated in the noise due to weak energy. Through deconvolution calculation, the original reflection coefficient sequence is effectively recovered, the weak reflection coefficient is well protected, and stronger anti-noise performance is embodied.
FIG. 4(a) is a migration profile of a measured three-dimensional seismic survey of an oil field, where the survey line 1.5s-2.0s is the survey target section, the dominant frequency of the data is 38Hz, the CMP spacing is 40m, and the selected survey line length is 20 km. The wave group characteristics on the imaging section are rich, and the integral inclination of the stratum reaches 10 degrees. FIG. 4(b) and FIG. 4(c) are each L0Fast iterative hard threshold deconvolution of norm constraints and the results of processing it by the methods herein. It can be seen that the two methods obviously improve the overall resolution of the section, enhance the energy focusing performance of the in-phase axis and make the detail description more clear. Comparing the arrow positions in FIGS. 4(b) and 4(c), the method compares L0The fast iteration hard threshold deconvolution method based on norm constraint has a better recovery effect on weak reflection signals. Meanwhile, the calculation efficiency of the gradient sparse regularization deconvolution method is greatly improved based on the self-adaptability of the method, so that the method has better applicability.
To further compare the performance of the two deconvolution algorithms on computational efficiency, an iterative convergence curve of the two is given herein (fig. 5). Wherein the horizontal axis is iteration times, the vertical axis is normalized error, i.e. | | | y- ωnrn||2/||y-ω1r1||2. It can be seen that the deconvolution method based on gradient sparse regularization has a faster convergence rate. The same normalization error is achieved, the iteration times of the algorithm can be reduced to be within ten times from hundreds of times of fast iteration hard threshold algorithm, the iteration stop condition is met earlier, smaller errors can be achieved finally, and the calculation precision is higher.
In this embodiment, while, for purposes of simplicity of explanation, the methodologies are shown and described as a series of acts, it is to be understood and appreciated that the methodologies are not limited by the order of acts, as some acts may, in accordance with one or more embodiments, occur in different orders and/or concurrently with other acts from that shown and described herein or not shown and described herein, as may be understood by those of ordinary skill in the art.
It is noted that references in the specification to "one embodiment," "an example embodiment," "some embodiments," etc., indicate that the embodiment described may include a particular feature, structure, or characteristic, but every embodiment may not necessarily include the particular feature, structure, or characteristic. Moreover, such phrases are not necessarily referring to the same embodiment. Further, when a particular feature, structure, or characteristic is described in connection with an embodiment, it is submitted that it is within the knowledge of one skilled in the art to effect such feature, structure, or characteristic in connection with other embodiments whether or not explicitly described.
The previous description of the disclosure is provided to enable any person skilled in the art to make or use the disclosure. Various modifications to the disclosure will be readily apparent to those skilled in the art, and the generic principles defined herein may be applied to other variations without departing from the spirit or scope of the disclosure. Thus, the disclosure is not intended to be limited to the examples and designs described herein but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims (10)

1. A time-frequency double-domain regularization sparse deconvolution method for narrow-band signals, comprising:
calculating the reflection coefficient gradient of the input seismic data;
calculating to obtain a conjugate gradient;
sparse regularization is carried out on the reflection coefficient conjugate gradient in a time domain;
iteratively updating the reflection coefficient based on the step length of the conjugate gradient of the reflection coefficient obtained by linear search;
the wavelets of the seismic data are updated based on the updated reflection coefficients.
2. A time-frequency two-domain regularization sparse deconvolution method for a narrowband signal as claimed in claim 1 wherein said reflection coefficient gradient ar is calculated based on the following equation:
Δr＝ωT(y-ωr)
in the formula, r is a reflection coefficient, omega is a seismic wavelet, y is observed seismic data, and T is a transpose of a wavelet matrix.
3. The time-frequency double-domain regularization sparse deconvolution method of claim 1, wherein said nonlinear conjugate gradient method utilizes a linear combination of a currently computed negative gradient direction and a previous conjugate gradient direction as a search direction for a next iteration, namely:
Figure FDA0003213311160000011
in the formula, Δ γkRepresenting the current conjugate gradient, Δ γk-1Denotes the last conjugate gradient, Δ rkRepresenting the current k-th reflection coefficient gradient, betakExpressing the linear weighting coefficient, and calculating the formula as follows:
Figure FDA0003213311160000012
4. a time-frequency two-domain regularization sparse deconvolution method for a narrowband signal as claimed in claim 1 wherein sparsely regularizing the reflection coefficient conjugate gradient in the time domain comprises:
Δγ＝Hλ(Δγ)
Hλfor an improved threshold function, the values are:
Hλ＝H1(H2)
Figure FDA0003213311160000021
Figure FDA0003213311160000022
where β ═ α κ + (1- α) γ, κ is a nonzero mean value of | θ | γ, γ is a maximum value of | θ |, α gives a relative threshold coefficient, α ∈ [0,1 ∈ γ]，
Figure FDA0003213311160000023
Is the derivative of theta with respect to time t.
5. The time-frequency two-domain regularization sparse deconvolution method for narrowband signals according to claim 1, characterized in that the reflection coefficient iteration step μ is calculated based on the following equation:
Figure FDA0003213311160000024
in the formula, ω is the seismic wavelet and Δ γ is the conjugate gradient of the reflection coefficient.
6. The time-frequency two-domain regularization sparse deconvolution method for narrowband signals according to claim 1, characterized in that the reflection coefficient r is updated based on the following equation:
r＝r+μΔγ
in the formula, r is a reflection coefficient, mu is an iteration step length, and delta gamma is a conjugate gradient of the reflection coefficient.
7. The time-frequency two-domain regularization sparse deconvolution method for narrow band signals of claim 1, wherein updating wavelet initial values for seismic data based on updated reflection coefficients comprises:
calculating the gradient of the wavelet;
calculating the conjugate gradient of the wavelet;
sparse regularization is carried out on the wavelet conjugate gradient in a frequency domain, and then the wavelet conjugate gradient after frequency domain sparse is obtained by utilizing inverse Fourier transform;
calculating a wavelet iteration step length based on the thinned wavelet conjugate gradient;
updating the wavelet initial value based on the wavelet iteration step size.
8. A time-frequency two-domain regularization sparse deconvolution method for narrow-band signals as claimed in claim 6 wherein said wavelet gradients are calculated based on the following equation:
Δω＝(y-ωr)rT
in the formula, reflection coefficients r and y are observed seismic data, omega is a seismic wavelet, r is a reflection coefficient, and T is a transpose of a wavelet matrix.
9. A time-frequency two-domain regularization sparse deconvolution method for narrow-band signals as claimed in claim 6 wherein said wavelet conjugate gradient is calculated based on the following equation:
Figure FDA0003213311160000031
in the formula (I), the compound is shown in the specification,
Figure FDA0003213311160000032
which is indicative of the current conjugate gradient,
Figure FDA0003213311160000033
represents the last conjugate gradient, Δ ωkRepresenting the current k-th sub-wave gradient, σkExpressing the linear weighting coefficient, and calculating the formula as follows:
Figure FDA0003213311160000034
10. the time-frequency two-domain regularization sparse deconvolution method for narrow-band signals of claim 6, wherein computing a wavelet iteration step based on a sparsified wavelet gradient comprises:
Figure FDA0003213311160000035
wherein r is a reflection coefficient,
Figure FDA0003213311160000036
is the conjugate gradient of the seismic wavelet;
updating the seismic wavelet based on:
Figure FDA0003213311160000041
where, ω is the seismic wavelet,
Figure FDA0003213311160000042
is the conjugate gradient of the seismic wavelet and λ is the iteration step of the wavelet.
CN202110945575.0A 2021-08-16 2021-08-16 Time-frequency double-domain regularized sparse deconvolution method and device for narrowband signals Active CN113552631B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110945575.0A CN113552631B (en) 2021-08-16 2021-08-16 Time-frequency double-domain regularized sparse deconvolution method and device for narrowband signals

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110945575.0A CN113552631B (en) 2021-08-16 2021-08-16 Time-frequency double-domain regularized sparse deconvolution method and device for narrowband signals

Publications (2)

Publication Number Publication Date
CN113552631A true CN113552631A (en) 2021-10-26
CN113552631B CN113552631B (en) 2023-11-03

Family

ID=78133959

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110945575.0A Active CN113552631B (en) 2021-08-16 2021-08-16 Time-frequency double-domain regularized sparse deconvolution method and device for narrowband signals

Country Status (1)

Country Link
CN (1) CN113552631B (en)

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150198729A1 (en) * 2014-01-13 2015-07-16 Cgg Services Sa Regularization of spatially aliased seismic data
US20160116620A1 (en) * 2014-10-24 2016-04-28 Ion Geophysical Corporation Methods and systems for seismic inversion and related seismic data processing
CN107238862A (en) * 2016-03-29 2017-10-10 中国石油化工股份有限公司 Reflectance factor method of estimation and device based on Bayes's inverting framework
CN108333628A (en) * 2018-01-17 2018-07-27 中国石油大学(华东) Elastic wave least square reverse-time migration method based on regularization constraint
CN110865409A (en) * 2019-12-02 2020-03-06 怀化学院 Seismic wave impedance inversion method based on wave impedance low-rank regularization
CN111368247A (en) * 2020-03-12 2020-07-03 电子科技大学 Sparse representation regularization prestack AVO inversion method based on fast orthogonal dictionary
CN113077386A (en) * 2021-04-06 2021-07-06 电子科技大学 Seismic data high-resolution processing method based on dictionary learning and sparse representation

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150198729A1 (en) * 2014-01-13 2015-07-16 Cgg Services Sa Regularization of spatially aliased seismic data
US20160116620A1 (en) * 2014-10-24 2016-04-28 Ion Geophysical Corporation Methods and systems for seismic inversion and related seismic data processing
CN107238862A (en) * 2016-03-29 2017-10-10 中国石油化工股份有限公司 Reflectance factor method of estimation and device based on Bayes's inverting framework
CN108333628A (en) * 2018-01-17 2018-07-27 中国石油大学(华东) Elastic wave least square reverse-time migration method based on regularization constraint
CN110865409A (en) * 2019-12-02 2020-03-06 怀化学院 Seismic wave impedance inversion method based on wave impedance low-rank regularization
CN111368247A (en) * 2020-03-12 2020-07-03 电子科技大学 Sparse representation regularization prestack AVO inversion method based on fast orthogonal dictionary
CN113077386A (en) * 2021-04-06 2021-07-06 电子科技大学 Seismic data high-resolution processing method based on dictionary learning and sparse representation

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
DAN JIN等: "The Effectiveness and Application of the Gradient Sparse Regularization-Based Deconvolution Method", IEEE GEOSCIENCE AND REMOTE SENSING LETTERS *
XIN XU等: "Multichannel Reflectivity Inversion With Sparse Group Regularization Based on HPPSG Algorithm", IEEE GEOSCIENCE AND REMOTE SENSING LETTERS *
梁东辉，陈生昌: "基于L0范数稀疏约束的地震数据反褶积", 石油物探 *

Also Published As

Publication number Publication date
CN113552631B (en) 2023-11-03

Similar Documents

Publication Publication Date Title
Chen Iterative deblending with multiple constraints based on shaping regularization
US6636810B1 (en) High-resolution Radon transform for processing seismic data
CN108919347A (en) Seismic signal stochastic noise suppression method based on vmd
Zhang et al. A patch based denoising method using deep convolutional neural network for seismic image
CN111368247B (en) Sparse representation regularization prestack AVO inversion method based on fast orthogonal dictionary
CN110208862B (en) Seismic inversion method based on mixed high-order fractional-order ATpV sparse regularization
CN110542923A (en) Rapid high-precision post-stack seismic impedance inversion method
CN107179550B (en) A kind of seismic signal zero phase deconvolution method of data-driven
Chen et al. Iterative deblending using shaping regularization with a combined PNMO-MF-FK coherency filter
CN110244353B (en) Seismic data regularization method based on sparse norm optimization algorithm
CN114418886A (en) Robustness denoising method based on deep convolution self-encoder
CN113077386A (en) Seismic data high-resolution processing method based on dictionary learning and sparse representation
CN113866826B (en) Hybrid domain seismic migration hessian matrix estimation method
Nose-Filho et al. Algorithms for sparse multichannel blind deconvolution
CN113552631A (en) Time-frequency double-domain regularization sparse deconvolution method and device for narrow-band signals
CN115327624B (en) Inversion method and inversion system for seismic wavelets and reflection coefficients
CN108226996B (en) Self-adaptive anisotropic frequency division partition filtering method based on energy frequency band distribution
CN111856568A (en) MWV model-based frequency domain multi-channel reflection coefficient joint inversion method and system
CN111856559B (en) Multi-channel seismic spectrum inversion method and system based on sparse Bayes learning theory
Chen* Deblending by iterative orthogonalization and seislet thresholding
CN113156514B (en) Seismic data denoising method and system based on dominant frequency wavenumber domain mean value filtering
CN110596756B (en) Desert seismic exploration noise suppression method based on self-adaptive mixed complex diffusion model
CN109991657B (en) Seismic data high-resolution processing method based on inverse binary recursive singular value decomposition
CN109856673B (en) High-resolution Radon transformation data separation technology based on dominant frequency iterative weighting
CN112363217A (en) Random noise suppression method and system for seismic data

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