CN105334537A - Primary wave and multiple wave separation method based on alternative splitting Bregman iterative algorithm - Google Patents

Abstract

The invention belongs to the field of seismic signal processing in seismic exploration technologies, and specifically discloses a primary wave and multiple wave separation method based on an alternative splitting Bregman iterative algorithm. With regard to a multiple wave self-adaptive subtraction method based on a 3D matched filter, the primary wave and multiple wave separation method utilizes the alternative splitting Bregman iterative algorithm to solve an optimization problem of applying sparsity constraint on primary waves, achieves the estimation of the 3D matched filter, and utilizes the estimated 3D matched filter to separate the primary waves and the multiple waves in a 3D data window in a self-adaptive manner. Compared with the traditional iterative reweighted least squares algorithm, the alternative splitting Bregman iterative algorithm adopted by the primary wave and multiple wave separation method only needs to calculate matrix-matrix multiplication and matrix inversion once when estimating the 3D matched filter at each 3D data window, can effectively reduce the calculation complexity of solving the optimization problem, and improves the calculating efficiency of primary wave and multiple wave self-adaptive separation.

Description

Based on the primary reflection and the multiple reflection separation method that alternately divide Bregman iterative algorithm

Technical field

The invention belongs to the seismic data processing field in seismic exploration technique, being specifically related to a kind of primary reflection and multiple reflection separation method based on alternately dividing Bregman iterative algorithm.

Background technology

SRME (SurfaceRelatedMultipleElimination) is the multiple reflection drawing method extensively adopted in offshore shooting.The multiple reflection of SRME prediction and true multiple reflection life period and spatial diversity, can utilize matched filter prediction multiple reflection self-adaptation from raw data to be deducted.Usually, along with the increase (from 1D, 2D to 3D matched filter) of matched filter dimension, repeatedly wave energy better from raw data self-adaptation deduct, but need more computing time.Mathematical model based on the multiple reflection self-adaptation subtractive method of 3D matched filter is (Li, Z., andW.Lu, 2013, Adaptivemultiplesubtractionbasedon3Dblindseparationofcon volvedmixtures:Geophysics, 78, V251-V266):

v＝d-Hx(1)

Wherein, v represents estimation primary reflection, and d represents raw data, and x represents 3D matched filter, and H represents the convolution matrix of prediction multiple reflection.

Based on the multiple reflection self-adaptation subtractive method of 3D matched filter in overlapped 3D data window, by estimating that 3D matched filter carrys out self-adaptation and is separated primary reflection and multiple reflection.For estimating 3D matched filter, traditional multiple reflection self-adaptation subtractive method applies energy minimization constraint to estimation primary reflection.In addition, the instability estimated for avoiding wave filter, filter coefficient is also supposed to meet energy minimization constraint.Corresponding optimization problem is:

$\arg \underset{x}{min}\[\left|\right|d-Hx|{|}_2^2+\mathrm{\μ}|\left|x\right|{|}_2^2\]---\left(2\right)$

Wherein, μ is regularization parameter.3D matched filter in equation (2) can adopt least-squares algorithm to solve:

x＝(H ^TH+μI) ^-1H ^Td(3)

Wherein, I is unit matrix.

Least-squares algorithm needs the orthogonality of primary reflection and multiple reflection to suppose.When primary reflection with multiple reflection is overlapped or when having strong primary reflection lineups to exist, least-squares algorithm can produce remaining multiple reflection or cause the damage of primary reflection.For overcoming the shortcoming of orthogonality hypothesis, sparse constraint being applied to primary reflection and has been incorporated in multiple reflection self-adaptation subtractive method.In addition, suppose that 3D matched filter coefficient meets energy minimization constraint and guarantees the stability that 3D matched filter is estimated, corresponding optimization problem is:

$\arg \underset{x}{min}\[\left|\right|d-Hx|{|}_1+\mathrm{\λ}|\left|x\right|{|}_2^2\],---\left(4\right)$

Wherein, λ is regularization parameter.Iteration heavy weighted least square algorithm (Guitton, A., andD.J.Verschuur, 2004, AdaptivesubtractionofmultiplesusingtheL can be adopted ₁-norm:GeophysicalProspecting, 52,27-38) carry out in estimator (4) 3D matched filter.But the heavy weighted least square algorithm of iteration all needs calculating matrix-matrix multiple and matrix inversion in every single-step iteration, and computation complexity is higher.

Summary of the invention

For the above-mentioned technical matters existed in prior art, the present invention proposes a kind of primary reflection and multiple reflection separation method based on alternately dividing Bregman iteration, the method effectively can reduce the computation complexity of optimization problem primary reflection being applied to sparse constraint, improves the counting yield that primary reflection is separated with multiple reflection.

To achieve these goals, the present invention adopts following technical scheme:

Based on the primary reflection and the multiple reflection separation method that alternately divide Bregman iterative algorithm, comprise the steps:

A arranges variable initial value, needs the variable arranging initial value to comprise the time span T of 3D data window ₀, the space length X of 3D data window ₀, the road collection number Y of 3D data window ₀, the space length R of the time span K of 3D matched filter, 3D matched filter, the road collection number G of 3D matched filter, primary reflection threshold value s _β, damping factor β and iterations

B inputs raw data road and concentrates data d in a pending 3D data window, then utilizes the parameter T of prediction multiple reflection data and 3D data window ₀, X ₀and Y ₀, 3D matched filter parameter K, R and G construct convolution matrix H, and adopt Cleskey decomposition computation inverse matrix

C, based on alternately dividing Bregman iterative algorithm, utilizes the inverse matrix that step b obtains the data d of all 3D data windows of seismic channel set is processed one by one;

D judges that raw data road concentrates the data d in all 3D data windows whether to be all disposed; If not, step c is returned; If be all disposed, then 3D Hanning window is first adopted to be weighted by the primary reflection estimated in each 3D data window, and a data volume that permeates then a data volume that permeated by 3D Hanning window is in the same way adopted final primary reflection estimated result is expressed as: wherein ,/represent the operation of being divided by of element one by one.

Preferably, described step c specifically comprises:

C1 arranges number of iterations m=1, utilizes the inverse matrix that step b obtains ask for the initial estimate v of primary reflection ⁽⁰⁾: $v^{\left(0\right)}=d-H\left(\stackrel{\‾}{H}\right(H^{T}d\left)\right),$ And make b ₁ ⁽⁰⁾=0, u ₁ ⁽⁰⁾=0;

C2 compute vector y ^(m)=v ^(m-1)+ b ₁ ^(m-1);

C3 is to y ^(m)utilize following distance operator prox compute vector u ₁ ^(m):

${u_1}^{\left(m\right)}={\mathrm{prox}}_s^1\left(y^{\left(m\right)}\right)=\left\{{\left(\right|y_{i,j,n}^{\left(m\right)}|-s_{\mathrm{\β}}{C}^{\left(m\right)})}_{+}\mathrm{sgn}\left({y_{i,j,n}}^{\left(m\right)}\right)\right\},$

Wherein, i=1,2 ..., T ₀, j=1,2 ..., X ₀, n=1,2 ..., Y ₀, s=s _βc ^(m), 0 < s _β< 1; y ^(m)={ y _{i, j, n} ^(m), y _{i, j, n} ^(m)represent vectorial y ^(m)in under be designated as the element of (i, j, n), C ^(m)=max (| y _{i, j, n} ^(m)|), ${\left(b\right)}_{+}=\left\{\begin{array}{cc}b,& b\&GreaterEqual;0\\ 0,& b<0\end{array}\right.,$ $\mathrm{sgn}\left(b\right)=\left\{\begin{array}{cc}1,& b>0\\ 0,& b=0\\ -1,& b<0\end{array}\right.;$

C4 compute vector b ₁ ^(m)=b ₁ ^(m-1)-[u ₁ ^(m)-v ^(m-1)];

C5 calculates 3D matched filter $X^{\left(m\right)}=\stackrel{\&OverBar;}{H}\left(H^T\right(d+{b_1}^{\left(m\right)}-{u_1}^{\left(m\right)}\left)\right);$

C6 estimates primary reflection v ^(m)=d-Hx ^(m);

C7 makes m=m+1, if turn back to step c2; If export the primary reflection result that current 3D data window is estimated.

Tool of the present invention has the following advantages:

For the multiple reflection self-adaptation subtractive method based on 3D matched filter, the present invention utilizes alternately division Bregman iterative algorithm to solve optimization problem primary reflection being applied to sparse constraint, realize the estimation of 3D matched filter, and utilize the primary reflection in the 3D matched filter self-adaptation separation 3D data window estimated and multiple reflection.Compared to the heavy weighted least square algorithm of traditional iteration, alternately division Bregman iterative algorithm in the present invention is when 3D matched filter estimated by each 3D data window, only need calculate a matrix-matrix multiple and matrix inversion, effectively can reduce the computation complexity of optimization problem, improve the counting yield that primary reflection is separated with multiple reflection self-adaptation.

Accompanying drawing explanation

Fig. 1 is based on alternately dividing the primary reflection of Bregman iterative algorithm and the process flow diagram of multiple reflection separation method in the present invention;

Fig. 2 a is the common offset road collection figure of raw data;

Fig. 2 b is the common offset road collection figure of prediction multiple reflection;

Fig. 3 a is the common offset road collection figure of the primary reflection estimated result based on least-squares algorithm;

Fig. 3 b is the common offset road collection figure of the primary reflection estimated result based on the heavy weighted least square algorithm of iteration;

Fig. 3 c is the common offset road collection figure based on the primary reflection estimated result alternately dividing Bregman iterative algorithm;

Fig. 4 a for based on least-squares algorithm remove the common offset road collection figure of multiple reflection;

Fig. 4 b for based on the heavy weighted least square algorithm of iteration remove the common offset road collection figure of multiple reflection;

Fig. 4 c for based on alternately division Bregman iterative algorithm remove the common offset road collection figure of multiple reflection;

Fig. 5 a is amplification display result figure (corresponding to black box in Fig. 2 a) of raw data;

Fig. 5 b is the amplification display result figure of prediction multiple reflection;

Fig. 5 c for estimated by least-squares algorithm primary reflection amplification display result figure;

Fig. 5 d is for weighing the amplification display result figure of primary reflection estimated by weighted least square algorithm based on iteration;

Fig. 5 e is for showing result figure based on the amplification of alternately dividing primary reflection estimated by Bregman iterative algorithm.

Embodiment

Basic thought of the present invention is: 3D data window ground is separated primary reflection and multiple reflection one by one, builds optimization problem primary reflection being applied to sparse constraint: wherein, d is raw data, and x represents 3D matched filter, and H represents the convolution matrix of prediction multiple reflection, and λ is regularization factors.Solve optimization problem in above formula to estimate 3D matched filter, and adopt the 3D matched filter estimated to be separated with multiple reflection primary reflection in 3D data window, finally the primary reflection estimated result in all 3D data windows is merged, obtain final primary reflection estimated result.

Below in conjunction with accompanying drawing and embodiment, the present invention is described in further detail:

As shown in Figure 1, based on the primary reflection and the multiple reflection separation method that alternately divide Bregman iterative algorithm, step is comprised:

A arranges variable initial value, needs the variable arranging initial value to comprise the time span T of 3D data window ₀, the space length X of 3D data window ₀, the road collection number Y of 3D data window ₀, the space length R of the time span K of 3D matched filter, 3D matched filter, the road collection number G of 3D matched filter, primary reflection threshold value s _β, damping factor β and iterations

B inputs raw data road and concentrates data d in a pending 3D data window, then utilizes the parameter T of prediction multiple reflection data and 3D data window ₀, X ₀and Y ₀, 3D matched filter parameter K, R and G construct convolution matrix H, the line number of matrix H is T ₀x ₀y ₀, columns is KRG, and adopts Cleskey decomposition computation inverse matrix

The inverse matrix H that c utilizes step b to obtain is to the data d process of the 3D data window of seismic channel set;

The present invention adopts alternately division Bregman iterative algorithm to solve 3D matched filter, realize primary reflection to be separated with the self-adaptation of multiple reflection, namely in every single-step iteration, alternately division Bregman iterative algorithm alternately solves and applies the optimization problem of sparse constraint and the optimization problem to the constraint of 3D matched filter applying energy minimization to primary reflection; Specifically, in every single-step iteration, alternately division Bregman iterative algorithm adopts distance operator to solve optimization problem primary reflection being applied to sparse constraint, adopts least-squares algorithm to solve optimization problem 3D matched filter being applied to energy minimization constraint; Suppose that the multiple proportion between regularization parameter simplifies iterative step simultaneously; The 3D matched filter self-adaptation of estimation is finally utilized to be separated primary reflection and multiple reflection.

Its concrete processing procedure is:

C1 arranges number of iterations m=1, utilizes the inverse matrix that step b obtains ask for the initial estimate v of primary reflection ⁽⁰⁾: $v^{\left(0\right)}=d-H\left(\stackrel{\&OverBar;}{H}\right(H^{T}d\left)\right),$ And make b ₁ ⁽⁰⁾=0, u ₁ ⁽⁰⁾=0;

C2 compute vector y ^(m)=v ^(m-1)+ b ₁ ^(m-1);

C3 is to y ^(m)utilize following distance operator prox compute vector u ₁ ^(m):

${u_1}^{\left(m\right)}={\mathrm{prox}}_s^1\left(y^{\left(m\right)}\right)=\left\{{\left(\right|{y_{i,j,n}}^{\left(m\right)}|-s_{\mathrm{\&beta;}}{C}^{\left(m\right)})}_{+}\mathrm{sgn}{{(}_{i,j,n}}^{\left(m\right)}\right)\},$

Wherein, i=1,2 ..., T ₀, j=1,2 ..., X ₀, n=1,2 ..., Y ₀, s=s _βc ^(m), 0 < s _β< 1; y ^(m)={ y _{i, j, n} ^(m), y _{i, j, n} ^(m)represent vectorial y ^(m)in under be designated as the element of (i, j, n), C ^(m)=max (| y _{i, j, n} ^(m)|), ${\left(b\right)}_{+}=\left\{\begin{array}{cc}b,& b\&GreaterEqual;0\\ 0,& b<0\end{array}\right.,$ $\mathrm{sgn}\left(b\right)=\left\{\begin{array}{cc}1,& b>0\\ 0,& b=0\\ -1,& b<0\end{array}\right.;$

C4 compute vector b ₁ ^(m)=b ₁ ^(m-1)-[u ₁ ^(m)-v ^(m-1)];

C5 calculates 3D matched filter $X^{\left(m\right)}=\stackrel{\&OverBar;}{H}\left(H^T\right(d+{b_1}^{\left(m\right)}-{u_1}^{\left(m\right)}\left)\right);$

C6 estimates primary reflection v ^(m)=d-Hx ^(m);

C7 makes m=m+1, if turn back to step c2; If export the primary reflection result of current data window estimation.

D judges that raw data road concentrates the data d in all data windows whether to be all disposed; If not, step c is returned; If be all disposed, then 3D Hanning window is first adopted to be weighted by the primary reflection estimated in each 3D data window, and a data volume that permeates then a data volume that permeated by 3D Hanning window is in the same way adopted final primary reflection estimated result is expressed as: wherein ,/represent the operation of being divided by of element one by one.

Herein, same mode refers to and the data of same time, space and Dao Ji position is added.

2D real data is utilized to verify the validity of the inventive method below:

Fig. 2 a is the common offset road collection figure of raw data.Fig. 2 b is the common offset road collection figure of prediction multiple reflection.Fig. 3 a is the common offset road collection figure of the primary reflection estimated result based on least-squares algorithm.Fig. 3 b is the common offset road collection figure of the primary reflection estimated result based on the heavy weighted least square algorithm of iteration.Fig. 3 c is the common offset road collection figure based on the primary reflection estimated result alternately dividing Bregman iterative algorithm.Fig. 4 a for based on least-squares algorithm remove the common offset road collection figure of multiple reflection.Fig. 4 b for based on the heavy weighted least square algorithm of iteration remove the common offset road collection figure of multiple reflection.Fig. 4 c for based on alternately division Bregman iterative algorithm remove the common offset road collection figure of multiple reflection.Obtain close primary reflection estimated result in Fig. 3 a, Fig. 3 b and Fig. 3 c, be respectively 8.88 seconds, 37.10 seconds and 9.63 seconds computing time.Compared to the heavy weighted least square algorithm of iteration, the alternately division Bregman iterative algorithm in the present invention, while maintenance primary reflection estimated accuracy, can improve counting yield effectively.In addition, Fig. 5 a is the amplification display result figure of raw data, corresponding to the black box in Fig. 2 a.Fig. 5 b is the amplification display result figure of prediction multiple reflection.Fig. 5 c for estimated by least-squares algorithm primary reflection amplification display result figure.Fig. 5 d is for weighing the amplification display result figure of primary reflection estimated by weighted least square algorithm based on iteration.Fig. 5 e is for showing result figure based on the amplification of alternately dividing primary reflection estimated by Bregman iterative algorithm.Wherein, the white arrow in Fig. 5 a to Fig. 5 e shows, relative to traditional least-squares algorithm, the alternately division Bregman iterative algorithm in the present invention can remove remaining multiple reflection better in primary reflection estimated result.Wherein, in Fig. 2 a to 2b, Fig. 3 a to 3c, Fig. 4 a to 4c and Fig. 5 a to 5e, horizontal ordinate TraceNumber represents Taoist monastic name, and ordinate Time represents the time, and unit is millisecond (ms).

Certainly; more than illustrate and be only preferred embodiment of the present invention; the present invention is not limited to enumerate above-described embodiment; should be noted that; any those of ordinary skill in the art are under the instruction of this instructions; made all equivalently to substitute, obvious form of distortion, within the essential scope all dropping on this instructions, protection of the present invention ought to be subject to.

Claims (2)

1., based on the primary reflection and the multiple reflection separation method that alternately divide Bregman iterative algorithm, it is characterized in that, comprise the steps:

A arranges variable initial value, needs the variable arranging initial value to comprise the time span T of 3D data window ₀, the space length X of 3D data window ₀, the road collection number Y of 3D data window ₀, the space length R of the time span K of 3D matched filter, 3D matched filter, the road collection number G of 3D matched filter, primary reflection threshold value s _β, damping factor β and iterations

B inputs raw data road and concentrates data d in a pending 3D data window, then utilizes the parameter T of prediction multiple reflection data and 3D data window ₀, X ₀and Y ₀, 3D matched filter parameter K, R and G construct convolution matrix H, and adopt Cleskey decomposition computation inverse matrix

C, based on alternately dividing Bregman iterative algorithm, utilizes the inverse matrix that step b obtains the data d of all 3D data windows of seismic channel set is processed one by one;

D judges that raw data road concentrates the data d in all 3D data windows whether to be all disposed; If not, step c is returned; If be all disposed, then 3D Hanning window is first adopted to be weighted by the primary reflection estimated in each 3D data window, and a data volume that permeates then a data volume that permeated by 3D Hanning window is in the same way adopted final primary reflection estimated result is expressed as: wherein ,/represent the operation of being divided by of element one by one.

2. primary reflection and the multiple reflection separation method based on alternately dividing Bregman iterative algorithm according to claim 1, it is characterized in that, described step c specifically comprises:

C1 arranges number of iterations m=1, utilizes the inverse matrix that step b obtains ask for the initial estimate v of primary reflection ⁽⁰⁾: $v^{\left(0\right)}=d-H\left(\stackrel{\&OverBar;}{H}\left(H^{T}d\right)\right),$ And make b ₁ ⁽⁰⁾=0, u ₁ ⁽⁰⁾=0;

C2 compute vector y ^(m)=v ^(m-1)+ b ₁ ^(m-1);

C3 is to y ^(m)utilize following distance operator prox compute vector u ₁ ^(m):

${u_1}^{\left(m\right)}={\mathrm{prox}}_s^1\left(y^{\left(m\right)}\right)=\left\{{\left(\left|{y_{i,j,n}}^m\right|-s_{\mathrm{\&beta;}}{C}^{\left(m\right)}\right)}_{+}\mathrm{sgn}\left({y_{i,j,n}}^{\left(m\right)}\right)\right\},$

Wherein, i=1,2 ..., T ₀, j=1,2 ..., X ₀, n=1,2 ..., Y ₀, s=s _βc ^(m), 0 < s _β< 1; y ^(m)={ y _{i, j, n} ^(m), y _{i, j, n} ^(m)represent vectorial y ^(m)in under be designated as the element of (i, j, n), C ^(m)=max (| y _{i, j, n} ^(m)|), ${\left(b\right)}_{+}=\left\{\begin{array}{cc}b,& b\&GreaterEqual;0\\ 0,& b<0\end{array}\right.,$ $\mathrm{sgn}\left(b\right)=\left\{\begin{array}{cc}1,& b>0\\ 0,& b=0\\ -1,& b<0\end{array}\right.;$

C4 compute vector b ₁ ^(m)=b ₁ ^(m-1)-[u ₁ ^(m)-v ^(m-1)];

C5 calculates 3D matched filter $x^{\left(m\right)}=\stackrel{\&OverBar;}{H}\left(H^T\left(d+{b_1}^{\left(m\right)}-{u_1}^{\left(m\right)}\right)\right);$

C6 estimates primary reflection v ^(m)=d-Hx ^(m);

C7 makes m=m+1, if turn back to step c2; If export the primary reflection result that current 3D data window is estimated.