CN108710851B - Seismic signal random noise attenuation method, terminal device and storage medium - Google Patents

Seismic signal random noise attenuation method, terminal device and storage medium Download PDF

Info

Publication number
CN108710851B
CN108710851B CN201810487212.5A CN201810487212A CN108710851B CN 108710851 B CN108710851 B CN 108710851B CN 201810487212 A CN201810487212 A CN 201810487212A CN 108710851 B CN108710851 B CN 108710851B
Authority
CN
China
Prior art keywords
seismic signal
random noise
noise attenuation
denoising
tgv
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.)
Expired - Fee Related
Application number
CN201810487212.5A
Other languages
Chinese (zh)
Other versions
CN108710851A (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.)
Minnan Normal University
Original Assignee
Minnan Normal University
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 Minnan Normal University filed Critical Minnan Normal University
Priority to CN201810487212.5A priority Critical patent/CN108710851B/en
Publication of CN108710851A publication Critical patent/CN108710851A/en
Application granted granted Critical
Publication of CN108710851B publication Critical patent/CN108710851B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing
    • G06F2218/04—Denoising
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/20—Image preprocessing
    • G06V10/30—Noise filtering

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Signal Processing (AREA)
  • General Engineering & Computer Science (AREA)
  • Multimedia (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Image Processing (AREA)

Abstract

The invention relates to a seismic signal random noise attenuation method based on overlapping group sparse generalized total variation, which comprises the following steps: s1, inputting a seismic signal graph G to be denoised; s2, initializing parameters in the solved and output image model; s3, F is obtained through calculation(k)，F(k)Obtaining a seismic signal diagram for k iterations, and judging whether to satisfy | | | F(k+1)‑F(k)||2/||F(k)||2>tol, if yes, turning to the step S11, otherwise, turning to the step S4; updating F, Vx，Vy(ii) a S5, judging whether the requirement is met
Figure DDA0002931501030000011
If yes, go to step S6, otherwise go to step S8; s6, utility model
Figure DDA0002931501030000012
Updating
Figure DDA0002931501030000013
S7, where n is n +1, and the process returns to step S5; s8, calculating
Figure DDA0002931501030000014
S9, using formula
Figure DDA0002931501030000015
Updating
Figure DDA0002931501030000016
S10, k equals k +1, and the process returns to step S3; and S11, outputting the denoised seismic signal diagram F. The method of the invention fully utilizes the neighborhood similarity of the first-order and second-order gradients of the image, and improves the difference between a smooth region and a boundary region, thereby improving the robustness of the denoising algorithm and obtaining better denoising performance compared with the classical TGV.

Description

Seismic signal random noise attenuation method, terminal device and storage medium
Technical Field
The invention relates to the field of image processing, in particular to a seismic signal random noise attenuation method based on overlapping group sparse generalized total variation, terminal equipment and a storage medium.
Background
There is a lot of random noise, such as gaussian noise, gamma noise, etc., in the actual seismic signal. Total variation regularization (TV) proved to be a regularization term effective in removing random noise, and has attracted a great deal of attention from students since it was proposed by Rudin et al. The total variation regularization term fully excavates the transverse and longitudinal gradient information of the two-dimensional image, well fits the prior knowledge of local smoothness, gradient sparsity and the like of the natural image, and is widely applied to the fields of seismic image denoising, seismic wave impedance inversion, weak and small target detection, super-resolution analysis and the like. Since the TV model assumes that the image is a piecewise smooth constant, the model has a more serious step effect. In order to suppress the step effect and improve the denoising effect of the TV model, Bredies, Kunisch and pack propose a generalized Total variant (TGV) model. Among the methods for removing the TV step effect, the generalized full variation model is proved to be a more effective method by theory and practice. The method simultaneously constrains the first-order gradient and the second-order gradient of the image, thereby effectively relieving the step effect of the total variation model. The model is a popularization of a total variation model, has a plurality of excellent mathematical properties such as convexity, lower semi-continuity, rotation invariance and the like, can approach any polynomial, arouses wide attention of scholars, and applies a TGV regular term to a plurality of fields.
In recent years, Selesnick and Chen propose overlapping sets of sparse regularization terms. The regular term is a non-separation regular term and can better keep the sparsity of the target function. The overlapping sparse regular term not only considers the sparsity of image difference domains, but also excavates neighborhood difference information of each point, thereby excavating the structural sparse characteristic of image gradient. The difference between the smooth region and the boundary region can be improved by overlapping the combined gradient, thereby suppressing the step effect of the TV model. Liu and the like use the work of Selesnick and Chen for reference, a one-dimensional overlapping group sparse regular term is popularized to a two-dimensional overlapping group sparse regular term, and the two-dimensional overlapping group sparse regular term is introduced into an anisotropic total variation model and is used for denoising and deconvolution of salt-and-pepper noise. Liu et al use an overlapping group sparsity regularization term for speckle noise removal. All the work is based on the popularization of the anisotropic total variation model, however, the anisotropic total variation model has no constraint of second-order difference, and the inhibition capability of the introduction of the overlapping group sparse gradient on the step effect is limited. It is noted that the TGV and the overlap group sparse shrinkage have different suppression mechanisms for the step effect, the former utilizes the first-order and second-order differential constraints of the image to alleviate the step effect, and the latter suppresses the step effect through the structural characteristics of the image gradient. At present, the cross-over study of the TGV model and the group sparse shrinkage is still in the starting stage, and it is noted that the classical TGV model does not consider the neighborhood structural characteristics of the image gradient.
Disclosure of Invention
The invention aims to provide a seismic signal random noise attenuation method based on overlapping group sparse generalized total variation, so as to solve the problems of the existing seismic signal denoising method. Therefore, the invention adopts the following specific technical scheme:
the seismic signal random noise attenuation method based on the overlapping group sparse generalized total variation comprises the following steps:
s1, inputting a seismic signal graph G to be denoised;
s2 solving output image model
Figure GDA0001873387840000021
Wherein F represents a denoised seismic signal diagram to be output,
Figure GDA0001873387840000022
a two-dimensional inverse fast fourier transform operator is represented,
Figure GDA0001873387840000023
Figure GDA0001873387840000024
Figure GDA0001873387840000025
Figure GDA0001873387840000031
Figure GDA0001873387840000032
Z1＝Kh*F-Vx，Z2＝Kv*F-Vy，Z3＝Kh*Vx，Z4＝Kv*Vy，Z5＝Kv*Vx+Kh*Vy，μ0＝μα0，μ1＝μα1，Λi(i ═ 1,2, 5) is the shrunken Lagrange multiplier, Kh＝[1,-1]，
Figure GDA0001873387840000033
μ denotes the regular coefficient, Vx,VyRespectively represents the intermediate variable of the overlapping group sparse regular term in the horizontal and vertical directions, alpha0＝1,α1＝0.5,β0＝1,β1＝1，
Figure GDA0001873387840000034
The fourier transform of the expression variable, the notation represents the point multiplier, and the parameters specifically needed to be initialized are as follows: k, Zi,Λi,F,μ0,μ1γ, tol, where K is the number of overlapping combinations, and tol is the threshold for the end of the iteration of the algorithm;
s3, according to
Figure GDA0001873387840000035
Is calculated to obtain F(k)，F(k)Obtaining a seismic signal diagram for k iterations, and judging whether to satisfy | | | F(k+1)-F(k)||2/||F(k)||2Tol, if yes, go to step S11, otherwise go to step S4, whichIn | · | non-conducting filament2Represents the matrix element L2A norm;
s4, using formula
Figure GDA0001873387840000036
Updating F, Vx,Vy；
S5, judging whether the requirement is met
Figure GDA0001873387840000037
If yes, go to step S6, otherwise go to step S8, wherein,
Figure GDA0001873387840000038
wherein
Figure GDA0001873387840000039
Representing an identity matrix, mat represents a vector matrixing operator,
Figure GDA00018733878400000310
represents Z after the k-th outer loop and the n-th inner loop (group sparse contraction iteration loop) are updated iterativelyiAnd is and
Figure GDA00018733878400000311
zirepresents ZiColumn vector form of (1);
s6, utility model
Figure GDA00018733878400000312
Updating
Figure GDA00018733878400000313
S7, where n is n +1, and the process returns to step S5;
s8, calculating
Figure GDA0001873387840000041
S9, using formula
Figure GDA0001873387840000042
Updating
Figure GDA0001873387840000043
S10, k equals k +1, and the process returns to step S3;
and S11, outputting the denoised seismic signal diagram F.
Further, K ═ 3, μ0＝1,μ1＝0.5,γ＝0.1,tol＝10-4。
Furthermore, the invention also provides a terminal device for seismic signal random noise attenuation, comprising a memory, a processor and a computer program stored in the memory and operable on the processor, wherein the processor implements the steps of the seismic signal random noise attenuation method as described above when executing the computer program.
Furthermore, the invention also provides a computer readable storage medium, in which a computer program is stored, wherein the computer program, when being executed by a processor, realizes the steps of the seismic signal random noise attenuation method as described above.
By adopting the technical scheme, the invention has the beneficial effects that: the method of the invention fully utilizes the neighborhood similarity of the first-order and second-order gradients of the image, and improves the difference between a smooth region and a boundary region, thereby improving the robustness of the denoising algorithm, obtaining better denoising performance compared with the classical TGV, and being particularly suitable for the denoising processing of the seismic signal image.
Drawings
FIG. 1 is a flow chart of a seismic signal random noise attenuation method based on overlapping group sparse generalized total variation according to the present invention;
FIG. 2 is a two-dimensional overlapping group sparse schematic in which (a) a smooth region contaminated with noise; (b) a boundary region free of noise contamination;
FIG. 3 is a schematic diagram of a synthetic seismic signal;
FIG. 4 is a graph of the denoising effect of four denoising algorithms, wherein (a) the original image; (b) a noisy map; (c) an ATV denoising result graph; (d) an ATV-OGS denoising result graph; e) a TGV-L1 denoising result graph; (f) a TGV-OGV denoising result graph;
FIG. 5 is a one-pass denoising contrast map of four denoising algorithms, wherein, (a) the original image; (b) a noisy map; (c) an ATV denoising result graph; (d) an ATV-OGS denoising result graph; e) a TGV-L1 denoising result graph; (f) a TGV-OGV denoising result graph;
FIG. 6 is a graph comparing PSNR with SSIM, where (a) PSNR for different Ks; (b) SSIM of different K;
FIG. 7 is a comparison graph of the four denoising algorithms denoised, wherein (a) is a two-dimensional seismic signal graph; (b) a noisy map; (c) an ATV denoising result graph; (d) an ATV-OGS denoising result graph; (e) a TGV-L1 denoising result graph; (f) and (5) a TGV-OGV denoising result graph.
Detailed Description
To further illustrate the various embodiments, the invention provides the accompanying drawings. The accompanying drawings, which are incorporated in and constitute a part of this disclosure, illustrate embodiments of the invention and, together with the description, serve to explain the principles of the embodiments. Those skilled in the art will appreciate still other possible embodiments and advantages of the present invention with reference to these figures.
1. Preparing knowledge:
1.1 second order TGV model
The standard second order TGV model is defined as equation (1). For ease of discussion, the processed image is assumed to be a square matrix.
Figure GDA0001873387840000051
Wherein
Figure GDA0001873387840000052
Representing the de-noised image to be solved,
Figure GDA0001873387840000053
in order to view the image of the seismic signal,
Figure GDA0001873387840000054
is the intermediate variable of the TGV regularization term.
Figure GDA0001873387840000061
In order to be a gradient operator, the method comprises the following steps,
Figure GDA0001873387840000062
the lateral gradient operator is represented by the lateral gradient operator,
Figure GDA0001873387840000063
representing a longitudinal gradient operator.
Figure GDA0001873387840000064
Representing fidelity terms, TGV2(f) And the second-order generalized total variation regular term is represented, and mu represents a regular coefficient and is used for balancing the fidelity term and the regular term.
Figure GDA0001873387840000065
To increase the algorithm speed using fourier transforms, we assume that the image satisfies periodic boundary conditions. To facilitate the subsequent discussion, we write the second order TGV model in matrix form and represent the differential form as a convolution form. That is to say that the first and second electrodes,
Figure GDA0001873387840000066
wherein F, G, Vx,Vy∈RN×N，Kh＝[1,-1]，
Figure GDA0001873387840000067
Symbol represents two-dimensional matrix convolution operation and satisfies
Figure GDA0001873387840000068
vec is the vectorization operator. (ii) (| · non-conducting phosphor) in formula2Represents the matrix element L2And (4) norm.
1.2 two-dimensional matrix overlapping group sparse neighbor operator
Suppose V0For the matrix to be shrunk, its overlapped set of sparse neighboring operators is recorded as,
Figure GDA0001873387840000069
wherein
Figure GDA00018733878400000610
Representing overlapping groups of sparse regularization terms.
Figure GDA00018733878400000611
An overlapping group sparse matrix of scale K × K is represented, defined as follows,
Figure GDA00018733878400000612
wherein.
Figure GDA00018733878400000613
Representing the lower integer operator.
As can be seen from equation (4), the overlapping group sparse regularization term takes into account the neighborhood K of the image gradient information2And points are adopted, so that the neighborhood structure similarity of the image gradient is fully mined. To demonstrate the denoising mechanism of the overlapped set of sparse regularization terms, we present a set of schematic diagrams as shown in fig. 2, assuming that K is 3 in fig. 2. Assume that fig. 2(a) represents the longitudinal gradient of a smooth region of an image, and that (i) and (ii) are heavily contaminated with noise. Fig. 2(b) shows the longitudinal gradient of the edge region in the image without noise pollution, and the two pixel points (c) and (d) show the pixel gradient of the edge region. In fig. 2, the light color circles represent pixel points with lower gray values, and the dark color circles represent pixel points with higher gray values, and for convenience of discussion, it is assumed that the gray value of the light color area is 0, and the gray value of the dark color area is 1. In fig. 2(a), if the conventional TV denoising method is adopted, since the gray values of the first and second pixels are equivalent to the gradient values of the third and fourth pixels at the boundary, the two pixels are very easily determined as the edge of the image by mistake and are retained, and the noise cannot be effectively removed. In the conventional TV model, the 4 points in question have exactly the same lateral and vertical differences. No matter how large the threshold is, the 4 point gradients will have consistent change regularity.
It is noted that in a smooth region, the probability that consecutive points are simultaneously contaminated by high noise is extremely small, and the noise can be removed by utilizing the structural similarity. Obviously, the combined gradient of the pixels (i) and (ii) obtained from the formula (4) is
Figure GDA0001873387840000071
The gradients of the two pixels are
Figure GDA0001873387840000072
At this time, as long as we set the threshold to be
Figure GDA0001873387840000073
The noise of the pixels can be well removed, and the gradient of the pixels is shrunk to
Figure GDA0001873387840000074
The image boundaries are well preserved. In conclusion, when the gradient of the pixel points is calculated in the overlapping combination mode, the difference between the high-noise pollution points in the smooth region and the pixel points in the boundary region can be highlighted, so that the denoising is performed more robustly.
Equation (3) can be efficiently calculated using an optimization minimization (MM) algorithm. According to the MM algorithm, to minimize p (V), a function is first found, Q (V, U) ≧ p (V) for all V, U, and the equality holds if and only if U ═ V. Accordingly, for optimization with the minimum value of Q (V, U) for each calculation being P (V), the calculation of equation (3) can be converted into an iteration,
Figure GDA0001873387840000081
considering the particularity of the group sparse regularization term, and noting the existence of the following inequality,
Figure GDA0001873387840000082
wherein an equal sign holds if and only if U ═ V.
By observing the formula (3) in combination with the formula (6), the compounds can be obtained
Figure GDA0001873387840000083
The optimization term of (a) is shown as the following formula,
Figure GDA0001873387840000084
we rewrite S (V, U) to,
Figure GDA0001873387840000085
where V is the vector form of matrix V, and c (u) is independent of V and can be considered a constant term with respect to V.
Figure GDA0001873387840000086
Is a diagonal matrix, whose diagonal elements are defined as follows,
Figure GDA0001873387840000087
combining equation (5) and equation (7), equation (3) can be transformed into the following iterative optimization problem,
Figure GDA0001873387840000088
the above equation is a quadratic programming problem, the iterative optimal solution of which is shown below,
V(k+1)＝mat{(I+γD2(V(k)))-1v0} (11)
wherein
Figure GDA0001873387840000089
Representing the identity matrix, v0Is V0In the form of vectors, mat denotes a vector matrixing operator。
2. Proposed model and solving method thereof
2.1 improved TGV model based on overlapping group sparsity
In this section, we improve the L1 constraint term in TGV to the overlap group sparse constraint term, so as to better mine the difference information of the first order gradient and the second order gradient of the image. We modify the second order TGV modeling as follows,
Figure GDA0001873387840000091
wherein TGV2-OGS(F) Representing a second-order generalized total variation regularization term based on an overlapping set sparse shrinkage operator,
Figure GDA0001873387840000092
representing overlapping groups of sparse regularization terms.
In contrast to the classical TGV model, the model proposed herein combines the first-order gradient with each pixel neighborhood K of the second-order gradient matrix2The information points are overlapped and combined, so that the structural characteristics of the first-order gradient matrix and the second-order gradient matrix are more fully mined. Obviously, the neighborhood gradient of each pixel point has certain similarity with the second-order neighborhood gradient, and the sparse priori knowledge of the structure is better characterized by the overlapped sparse model.
2.2 ADMM solution of proposed model
To solve the improved TGV model defined by equation (12), we solve the model using the ADMM framework, which transforms the complex problem into several simple sub-problems to solve by introducing decoupled split variables. The split variable is defined as: z1＝Kh*F-Vx，Z2＝Kv*F-Vy，Z3＝Kh*Vx，Z4＝Kv*Vy，Z5＝Kv*Vx+Kh*VyThe original problem translates into a constraint problem, namely:
Figure GDA0001873387840000093
according to the ADMM principle, the constraint problem described by equation (13) can be converted into an unconstrained augmented lagrange function, and the objective function becomes the following expression:
Figure GDA0001873387840000101
wherein mu0＝μα0，μ1＝μα1In general, alpha0＝1,α1＝0.5，Λi(i ═ 1,2, 5) is a shrunken Lagrangian multiplier,<X,Y>representing the inner product of two matrices X, Y.
2.2.1 F,Vx,VySub-problems
Under ADMM framework, separate between variables and their dual variables from triples F, Vx,VyAre decoupled. These three variables are coupled to each other, so it is necessary to establish a system of equations of three equations in a linear three-way manner, in which, for the sub-problem of F, its sub-objective function is degenerated as:
Figure GDA0001873387840000102
it is assumed herein that the process data satisfies periodic boundary conditions, and thus convolution calculations can be performed using a fast fourier transform. According to the convolution law, two matrixes are convoluted in a space domain and correspond to a frequency domain, the frequency spectrum of the convolution result is the dot product of the frequency spectrums of the two matrixes, and two sides of the formula (15) are subjected to Fourier transform:
Figure GDA0001873387840000103
wherein
Figure GDA0001873387840000104
Represents the fourier transform of the variable. The symbol represents a point multiplier. Order to
Figure GDA0001873387840000105
About
Figure GDA0001873387840000106
The derivative is zero to obtain
Figure GDA0001873387840000107
Order:
Figure GDA0001873387840000108
Figure GDA0001873387840000109
the mixture is obtained by finishing the raw materials,
Figure GDA00018733878400001010
for VxThe sub-problem, sub-objective function is:
Figure GDA0001873387840000111
similarly, Fourier transform is performed on two sides of the formula (19),
Figure GDA0001873387840000112
order to
Figure GDA0001873387840000113
About variables
Figure GDA0001873387840000114
The partial derivatives are calculated and made to be zero.
Figure GDA0001873387840000115
Order to
Figure GDA0001873387840000116
Figure GDA0001873387840000117
Is finished to obtain
Figure GDA0001873387840000118
For VyThe sub-problem, sub-objective function is:
Figure GDA0001873387840000119
fourier transformation is carried out on two sides of the formula (23),
Figure GDA00018733878400001110
order to
Figure GDA00018733878400001111
About variables
Figure GDA00018733878400001112
And (3) solving the partial derivatives, and setting zero to obtain:
Figure GDA00018733878400001113
order to
Figure GDA00018733878400001114
Figure GDA00018733878400001115
Is finished to obtain
Figure GDA0001873387840000121
By combining (18), (22), (26) to obtain information about F, Vx,VyThe system of equations for the three variables, that is,
Figure GDA0001873387840000122
F,Vx,Vythe solution can be solved using the crime law and the fast inverse fourier transform, i.e.:
Figure GDA0001873387840000123
the division in the above formula is a division by element points.
Figure GDA0001873387840000124
Representing a two-dimensional inverse fast fourier transform operator. Wherein
Figure GDA0001873387840000125
2.2.2 Zi(i-1, 2, 5) subproblem solving
For Z1A sub-problem, whose objective sub-function is,
Figure GDA0001873387840000126
according to formula (11) having1The update formula of (a) is as follows,
Figure GDA0001873387840000127
wherein the content of the first and second substances,
Figure GDA0001873387840000128
represents Z after the k-th outer loop and the n-th inner loop (group sparse contraction iteration loop) are updated iteratively1And is and
Figure GDA0001873387840000129
z1represents Z1In the form of a column vector. mat denotes the operator that converts the column vector into a matrix.
In the same way, Z2，Z3，Z4，Z5The update formula of (a) is as follows,
Figure GDA0001873387840000131
wherein the content of the first and second substances,
Figure GDA0001873387840000132
Figure GDA0001873387840000133
2.2.3 Λi(i-1, 2, 5) subproblem solving
For Z1Lagrange variable Λ of1The objective sub-function of which is,
Figure GDA0001873387840000134
the updated formula can be obtained by using the gradient ascent method,
Figure GDA0001873387840000135
similarly, the Lagrangian variable Λ2,Λ3,Λ4,Λ5Is shown as equation (34),
Figure GDA0001873387840000136
where γ is the learning rate.
All the sub-problems of the proposed model are solved so far. We summarize the algorithm as algorithm 1,
Figure GDA0001873387840000137
Figure GDA0001873387840000141
in summary, the steps of the seismic signal random noise attenuation method based on the overlapping group sparse generalized total variation of the present invention are shown in fig. 1, and are specifically described as follows:
s1, inputting a seismic signal graph G to be denoised;
s2, initializing and solving parameters of the output image model, wherein the specific parameters are as follows:
K,Zi＝0,Λi＝0(i＝1,2,,5),F＝0,μ0,β0,μ1,β1γ, tol, E ═ 1, K ═ 1, and n ═ 0 where tol is the threshold at which the algorithm iteration ends, and in one specific example K ═ 3, μ0＝1,β0＝1,μ1＝0.5,β1＝1,γ＝0.1,tol＝10-4,k＝1,n＝0；
S3, judging whether | | | F is satisfied(k+1)-F(k)||2/||F(k)||2If yes, go to step S11, otherwise go to step S4;
s4, updating F, V using equation (28)x,Vy；
S5, judging whether the requirement is met
Figure GDA0001873387840000142
If yes, go to step S6, otherwise go to step S8;
s6, updating by using equations (30) - (31)
Figure GDA0001873387840000143
S7, where n is n +1, and the process returns to step S5;
s8, calculating
Figure GDA0001873387840000144
S9, updating by equations (33) - (34)
Figure GDA0001873387840000145
S10, k equals k +1, and the process returns to step S3;
and S11, outputting the denoised seismic signal diagram F.
In addition, the invention also discloses a terminal device for seismic signal random noise attenuation, which comprises a memory, a processor and a computer program stored in the memory and capable of running on the processor, wherein the processor realizes the steps of the seismic signal random noise attenuation method when executing the computer program.
Further, the terminal device may be a desktop computer, a notebook, a palm computer, a cloud server, or other computing devices. The terminal device may include, but is not limited to, a processor, a memory. It will be understood by those skilled in the art that the above-mentioned structure of the terminal device is only an example of the terminal device for random noise attenuation of seismic signals, and does not constitute a limitation to the terminal device for random noise attenuation of seismic signals, and may include more or less components than the above, or combine some components, or different components, for example, the terminal device for random noise attenuation of seismic signals may further include an input-output device, a network access device, a bus, etc., which is not limited in this embodiment of the present invention.
Further, the Processor may be a Central Processing Unit (CPU), other general purpose Processor, a Digital Signal Processor (DSP), an Application Specific Integrated Circuit (ASIC), an off-the-shelf Programmable Gate Array (FPGA) or other Programmable logic device, discrete Gate or transistor logic device, discrete hardware component, or the like. The general purpose processor may be a microprocessor or the processor may be any conventional processor or the like, the processor being the control center for the terminal equipment for seismic signal random noise attenuation, various interfaces and lines connecting the various parts of the overall terminal equipment for seismic signal random noise attenuation.
The memory may be used to store the computer programs and/or modules, and the processor may implement the various functions of the terminal device for seismic signal random noise attenuation by executing or executing the computer programs and/or modules stored in the memory and invoking data stored in the memory. The memory may mainly include a program storage area and a data storage area, wherein the program storage area may store an operating system, an application program required for at least one function, and the like. In addition, the memory may include high speed random access memory, and may also include non-volatile memory, such as a hard disk, a memory, a plug-in hard disk, a Smart Media Card (SMC), a Secure Digital (SD) Card, a Flash memory Card (Flash Card), at least one magnetic disk storage device, a Flash memory device, or other volatile solid state storage device.
The invention also provides a computer readable storage medium having a computer program stored thereon, wherein the computer program when executed by a processor implements the steps of the method for random noise attenuation of seismic signals as described above.
The terminal device integrated modules/units for random noise attenuation of seismic signals may be stored in a computer readable storage medium if implemented in the form of software functional units and sold or used as a stand-alone product. Based on such understanding, all or part of the flow of the method according to the embodiments of the present invention may also be implemented by a computer program, which may be stored in a computer-readable storage medium, and when the computer program is executed by a processor, the steps of the method embodiments may be implemented. Wherein the computer program comprises computer program code, which may be in the form of source code, object code, an executable file or some intermediate form, etc. The computer-readable medium may include: any entity or device capable of carrying the computer program code, recording medium, usb disk, removable hard disk, magnetic disk, optical disk, computer Memory, Read-Only Memory (ROM), Random Access Memory (RAM), electrical carrier wave signals, telecommunications signals, software distribution medium, and the like. It should be noted that the computer readable medium may contain content that is subject to appropriate increase or decrease as required by legislation and patent practice in jurisdictions, for example, in some jurisdictions, computer readable media does not include electrical carrier signals and telecommunications signals as is required by legislation and patent practice.
3. Numerical experiment
In this section, the proposed algorithm is applied to seismic signals, the denoising performance of the proposed algorithm is objectively evaluated through a peak signal-to-noise ratio, a structural similarity indication coefficient and running time, and comprehensive comparison is made with ATV (active wavelet transform), an anisotropic overlapping group sparse total variation denoising method and other various TGV denoising methods.
3.1 Experimental platform and evaluation index
The most common evaluation indexes in the denoising field are peak signal-to-noise ratio, structural similarity indexes, relative errors and operation time. Wherein, the peak signal-to-noise ratio and the structural similarity index are shown as formula (35) and formula (36).
Figure GDA0001873387840000171
Where X is the original and Y is the reconstructed image.
Figure GDA0001873387840000172
Wherein u isXIs the average value of X, uYIs the average value of the values of Y,
Figure GDA0001873387840000173
is the variance of X and is,
Figure GDA0001873387840000174
is the variance of Y, σXYIs the covariance of X and Y, c1＝(k1L)2，c2＝(k2L)2Is a constant used to maintain the denominator non-zero, L512, k1＝0.05，k20.05. In the experiment, L is 255.
In the experiment of the section, four algorithms, namely ATV, TGV constrained by L1 norm, ATV-OGS and TGV-OGS, need to be compared, in order to ensure the objectivity and fairness of evaluation, the iteration stopping condition of the algorithm is defined as,
||F(k+1)-F(k)||2/||F(k)||2＜10-4 (37)
the hardware environment is 6G memory, and the main parameters of the CPU are main frequency 2.5G. The software platform is matlab 2016.
3.2 synthetic seismic signal denoising performance contrast test
This section presents synthetic seismic records obtained by convolving 20Hz rake wavelets with artificially synthesized reflection coefficients. We cut out the first 50 tracks of the composite record to observe the denoising effect, as shown in fig. 3.
3.2.1 overlap group sparse regularization term effect test
This section of experiment is intended to verify the validity of the overlap group sparse regularization term. We added 0 mean, standard deviation σ, white gaussian noise to the synthetic recordings. In this section, we respectively test the following algorithms, namely an anisotropic ATV model, an anisotropic variation model based on an overlapped group sparse shrinkage operator (ATV-OGS for short), an TGV model based on ADMM and L1 norm constraints (TGV-L1 for short), and the enhanced TGV denoising method based on overlapped group sparse shrinkage proposed herein. The test results are shown in table 1, and the indexes with the best indexes are marked by black bold. In this section of the experiment, we set K to 3.
TABLE 1 Algorithm Performance test Table
Figure GDA0001873387840000181
As can be seen from Table 1, the overlap group sparse regularization term not only improves the performance of the classical anisotropic TV model algorithm, but also improves the TGV-L1 denoising method. To reflect the effect intuitively, we show the denoising effect of the above four algorithms when σ is 30 in fig. 4.
By carefully observing the area outlined by the rectangular frame in fig. 4, it is obvious that, in the original image, the area is a smooth area, the ATV denoising method has a relatively obvious blocking effect in the area, while the blocking effect is well suppressed by using the overlapping group sparse regular term, and it is worth pointing out that the anisotropic overlapping group sparse method does not completely suppress the blocking effect of the ATV. And observing the denoising result of the TGV-L1, the TGV-L1 denoising method has a good inhibition effect on the blocking effect of the ATV. The method provided by the invention combines the advantages of sparse overlapping groups and TGV to obtain better denoising effect than ATV-OGS and TGV-L1, and as can be seen from the observation of FIG. 4(f), by carefully comparing the area outlined by the rectangular frame, it can be found that the denoising result obtained by the algorithm provided by the invention is closer to the original image in the area.
To further observe the denoising effect, we extract the first pass of the six subgraphs in fig. 4 to observe, as shown in fig. 5. It can be seen from fig. 5 that there is a significant step effect in ATV, and the group sparse regularization term better mitigates the step effect of ATV, but there are still more spikes as in fig. 5. Similar problems exist with the conventional TGV method. The method provided by the invention combines a group sparse shrinkage method on a TGV theoretical framework, so that the structural information of the data neighborhood gradient is mined, the advantages of the two are fully combined, and the reconstruction quality of the signal is greatly improved.
Through the experiment, the conclusion is drawn that the anisotropic overlapped sparse denoising method can effectively utilize the structural characteristics of the first-order gradient of the image, the first-order and second-order image gradients exist in the generalized fully-variable model, the gradients also contain neighborhood structural similarity, the overlapped sparse regular term is fully combined with the TGV model, and the denoising capability of the TGV model is further improved.
3.2.2 analysis of parameter sensitivity
The section tests and compares important parameters of the proposed algorithm and the overlapping combination number K to evaluate the overall effect of the algorithm. In this section, PSNR and SSIM are used to objectively evaluate the algorithm, and we continuously change K from 1 to 13 for different noise levels, and record PSNR and SSIM values, as shown in fig. 6.
As can be seen from fig. 6, as K increases, K has different roles in PSNR and SSIM under different noises, and it is obvious that at low noise (for example, σ ═ 10), K reaches 11, PSNR and SSIM both reach a peak, and PSNR begins to fall back if K increases. It can be seen that when the noise is low, the neighborhood information of the image plays a positive role in the performance of the algorithm, however, K should not be too large, otherwise, a region with severe image change in the neighborhood may be obtained, and in this case, PSNR and SSIM are reduced. When the noise is large, the structural characteristics of the neighborhood gradient are damaged more seriously, in this case, the performance of the algorithm is reduced due to a slightly large K, for example, when σ is 40, K is 7, and the highest values of PSNR and SSIM are obtained.
Notably, the TGV-OGS (K ═ 3) algorithm we show in section 3.2.2 has not yet reached its highest performance. The PSNR of the proposed algorithm reaches a maximum of 37.7846dB, 1.905dB above the TGV-L1 algorithm, when σ is 10.
3.3 actual seismic Signal denoising Performance test
In this subsection, we use two-dimensional seismic signals as test signals, which are shown in FIG. 7 (a). The seismic signal is a post-stack seismic signal that has been denoised. In this section, we add white gaussian noise with σ of 30 to the seismic signal, compare the noise removal performance indicators of the following algorithms ATV, ATV-OGS, TGV-L1 and TGV-OGS, and each indicator is recorded in table 2. Table 2 shows that the denoising effect of the algorithm proposed herein is best for the actual seismic signal.
TABLE 2 Algorithm Performance test Table
Figure GDA0001873387840000201
Fig. 7 shows the two-dimensional image denoising effect of various algorithms, and it can be seen that the algorithm proposed herein effectively removes artificially increased gaussian noise. Compared with the original seismic signal, the method provided by the invention can remove high-frequency interference in the seismic signal better, the blocking effect problem of the ATV method is avoided, and the in-phase axis has better transverse continuity.
4. Conclusion
An improved generalized total variation denoising algorithm is provided under an ADMM framework by starting from a group of sparse regular terms and combining the definition of generalized total variation, and is applied to seismic signal denoising. The algorithm fully utilizes the neighborhood similarity of the first-order and second-order gradients of the image, and improves the difference between a smooth region and a boundary region, so that the robustness of the denoising algorithm is improved, and better denoising performance compared with the classical TGV is obtained.
While the invention has been particularly shown and described with reference to a preferred embodiment, it will be understood by those skilled in the art that various changes in form and detail may be made therein without departing from the spirit and scope of the invention as defined by the appended claims.

Claims (4)

1. The seismic signal random noise attenuation method based on the overlapping group sparse generalized total variation comprises the following steps:
s1, inputting a seismic signal graph G to be denoised;
s2 solving output image model
Figure FDA0002948022850000011
Wherein F represents a denoised seismic signal diagram to be output,
Figure FDA0002948022850000012
a two-dimensional inverse fast fourier transform operator is represented,
Figure FDA0002948022850000013
Figure FDA0002948022850000014
Figure FDA0002948022850000015
Figure FDA0002948022850000016
Figure FDA0002948022850000017
Figure FDA0002948022850000018
Z1＝Kh*F-Vx，Z2＝Kv*F-Vy，Z3＝Kh*Vx，Z4＝Kv*Vy，Z5＝Kv*Vx+Kh*Vy，μ0＝μα0，μ1＝μα1，Λi(i ═ 1,2, …,5) is a shrunken Lagrangian multiplier, Kh＝[1,-1]，
Figure FDA0002948022850000019
μ denotes the regular coefficient, Vx,VyRespectively represents the intermediate variable of the overlapping group sparse regular term in the horizontal and vertical directions, alpha0＝1,α1＝0.5,β0＝1,β1＝1，
Figure FDA00029480228500000110
Fourier transform, sign of expression variable
Figure FDA00029480228500000111
The parameters representing the point multiplier, which need to be initialized specifically, are as follows: k, Zi,Λi,F,μ0,μ1γ, tol, where K is the number of overlapping combinations, and tol is the threshold for the end of the iteration of the algorithm;
s3, according to
Figure FDA00029480228500000112
Is calculated to obtain F(k)，F(k)Obtaining a seismic signal diagram for k iterations, and judging whether to satisfy | | | F(k+1)-F(k)||2/||F(k)||2If so, go to step S11, otherwise go to step S4, where | · | | survival |2Represents the matrix element L2A norm;
s4, using formula
Figure FDA0002948022850000021
Updating F, Vx,Vy；
S5, judging whether the requirement is met
Figure FDA0002948022850000022
If yes, go to step S6, otherwise go to step S8, wherein,
Figure FDA0002948022850000023
wherein
Figure FDA0002948022850000024
Representing an identity matrix, mat represents a vector matrixing operator,
Figure FDA0002948022850000025
z after the k-th outer loop and the n-th inner loop are updated iterativelyiAnd is and
Figure FDA0002948022850000026
zirepresents ZiColumn vector form of (1);
s6, utility model
Figure FDA0002948022850000027
Updating
Figure FDA0002948022850000028
S7, where n is n +1, and the process returns to step S5;
s8, calculating
Figure FDA0002948022850000029
S9, using formula
Figure FDA00029480228500000210
Updating
Figure FDA00029480228500000211
S10, k equals k +1, and the process returns to step S3;
and S11, outputting the denoised seismic signal diagram F.
2. The overlapping group sparse generalized total variation-based seismic signal random noise attenuation method of claim 1, wherein K is 3, μ0＝1,μ1＝0.5,γ＝0.1,tol＝10-4。
3. A terminal device for random noise attenuation of seismic signals, comprising a memory, a processor and a computer program stored in the memory and executable on the processor, characterized in that the processor, when executing the computer program, implements the steps of the method of random noise attenuation of seismic signals according to any of claims 1 to 2.
4. A computer-readable storage medium, in which a computer program is stored which, when being executed by a processor, carries out the steps of the method for random noise attenuation of seismic signals according to any one of claims 1 to 2.
CN201810487212.5A 2018-05-21 2018-05-21 Seismic signal random noise attenuation method, terminal device and storage medium Expired - Fee Related CN108710851B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810487212.5A CN108710851B (en) 2018-05-21 2018-05-21 Seismic signal random noise attenuation method, terminal device and storage medium

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810487212.5A CN108710851B (en) 2018-05-21 2018-05-21 Seismic signal random noise attenuation method, terminal device and storage medium

Publications (2)

Publication Number Publication Date
CN108710851A CN108710851A (en) 2018-10-26
CN108710851B true CN108710851B (en) 2021-04-16

Family

ID=63868325

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810487212.5A Expired - Fee Related CN108710851B (en) 2018-05-21 2018-05-21 Seismic signal random noise attenuation method, terminal device and storage medium

Country Status (1)

Country Link
CN (1) CN108710851B (en)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109146797B (en) * 2018-06-15 2019-10-25 闽南师范大学 A kind of impulsive noise ancient book image inpainting method sparse based on Lp pseudonorm and overlapping group
CN109934815B (en) * 2019-03-18 2023-04-14 电子科技大学 Tensor recovery infrared small target detection method combined with ATV constraint
CN110084756B (en) * 2019-04-15 2020-12-29 闽南师范大学 Image denoising method based on high-order overlapping group sparse total variation
CN110208862B (en) * 2019-07-04 2021-01-29 电子科技大学 Seismic inversion method based on mixed high-order fractional-order ATpV sparse regularization

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5282424A (en) * 1991-11-18 1994-02-01 Neill Gerard K O High speed transport system
CN104395779A (en) * 2012-07-10 2015-03-04 雪佛龙美国公司 System and method for estimating and attenuating noise in seismic data
CN104880730A (en) * 2015-03-27 2015-09-02 西安交通大学 Seismic data time-frequency analysis and attenuation estimation method based on Synchrosqueezing transform
CN107390267A (en) * 2017-07-27 2017-11-24 西安交通大学 A kind of seismic data attenuation compensation method of synchronous extruding transform domain

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5282424A (en) * 1991-11-18 1994-02-01 Neill Gerard K O High speed transport system
CN104395779A (en) * 2012-07-10 2015-03-04 雪佛龙美国公司 System and method for estimating and attenuating noise in seismic data
CN104880730A (en) * 2015-03-27 2015-09-02 西安交通大学 Seismic data time-frequency analysis and attenuation estimation method based on Synchrosqueezing transform
CN107390267A (en) * 2017-07-27 2017-11-24 西安交通大学 A kind of seismic data attenuation compensation method of synchronous extruding transform domain

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Four-directional fractional order total variation regularization for image denoising;Wu L 等;《Journal of Electronic imaging》;20170531;第1-2页 *
基于稀疏表示及正则约束的图像去噪方法综述;陈颖频 等;《数据采集与处理》;20180131;第33卷;第1-2页 *

Also Published As

Publication number Publication date
CN108710851A (en) 2018-10-26

Similar Documents

Publication Publication Date Title
CN108710851B (en) Seismic signal random noise attenuation method, terminal device and storage medium
Liu et al. Image restoration using total variation with overlapping group sparsity
Jin et al. Sparse and low-rank decomposition of a Hankel structured matrix for impulse noise removal
Zhang et al. 100+ times faster weighted median filter (WMF)
Ding et al. Total variation with overlapping group sparsity for deblurring images under Cauchy noise
US9443286B2 (en) Gray image processing method and apparatus based on wavelet transformation
Shukla et al. Generalized fractional filter-based algorithm for image denoising
CN102945548A (en) Directional pyramid filtering-based image processing method and device
Chen et al. Image denoising via local and nonlocal circulant similarity
Ullah et al. A new variational approach for multiplicative noise and blur removal
Raj et al. Medical image denoising using multi-resolution transforms
Hidane et al. Nonlinear multilayered representation of graph-signals
Li et al. A novel weighted anisotropic total variational model for image applications
Zhang et al. High-order total bounded variation model and its fast algorithm for Poissonian image restoration
Steidl Combined first and second order variational approaches for image processing
Jalab et al. Fractional conway polynomials for image denoising with regularized fractional power parameters
Ji et al. Image recovery via geometrically structured approximation
Zhao et al. Constant time texture filtering
Xu Sparse regularization with the ℓ 0 norm
Wang et al. A unidirectional total variation and second-order total variation model for destriping of remote sensing images
Chen et al. SAR image despeckling by using nonlocal sparse coding model
CN113139920B (en) Ancient book image restoration method, terminal equipment and storage medium
Li et al. Video denoising using shape-adaptive sparse representation over similar spatio-temporal patches
Satti et al. DIBS: Distance-and intensity-based separation filter for high-density impulse noise removal
CN113627357B (en) High-spatial-high-spectral-resolution intrinsic decomposition method and system for remote sensing image

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20210416