Curvelet domain multiple adaptive attenuation method based on time domain transform channel

Abstract

The invention relates to a curvelet domain multiple adaptive attenuation method based on a time domain transformation channel, which belongs to the field of seismic data processing, and specifically comprises the following steps of: firstly, generating a multi-order time domain integral and time domain differential record by utilizing a multiple record obtained by wave equation prediction, and transforming the record into a curvelet domain through curvelet transformation to obtain a curvelet domain; then applying multi-channel least square filtering in the domain to eliminate multiple components in the original data, and finally inversely transforming a multiple suppression result into a time-space domain. Field actual data processing results show that the method gives consideration to the adaptability of time domain transformation to multiple prediction errors and the separation advantage of curvelet transformation to primary waves and multiple waves, multiple wave interference can be effectively suppressed, and therefore the problem of multiple wave attenuation caused by inaccurate multiple wave prediction in complex sea areas is solved.

Description

Technical Field

[0001] The invention belongs to the field of seismic data processing, and in particular relates to a curvelet domain multiple wave adaptive attenuation method based on time domain transformation traces. Background Art

[0002] Marine seismic exploration is currently the most effective technology for marine oil and gas exploration. However, in marine seismic exploration, since the geological interfaces on the sea surface and on the seabed or below the seabed are all high wave impedance interfaces, strong multiple waves will be generated, which will seriously affect the signal-to-noise ratio of seismic data and bring difficulties to subsequent high-precision seismic imaging.

[0003] The multiple wave suppression of wave equation type multiple wave elimination methods can usually be divided into two main steps: multiple wave prediction and matching attenuation, that is, first, the multiple waves in the seismic data are predicted by wave theory, and then the predicted multiple waves are matched with the original data and attenuated to remove the multiple wave components. In this process, the multiple wave matching attenuation algorithm is a key step in determining the processing effect of the wave equation type multiple wave elimination method.

[0004] When multiple wave prediction is performed using the wave equation method (including the wave field extrapolation method, the free interface multiple wave attenuation method and the inverse scattering series method), the bandwidth of the multiple wave records becomes significantly narrower due to the calculation process of multi-channel weighted summation in the input records, and the main frequency moves from the high frequency end to the low frequency end, which seriously affects the multiple wave elimination effect of the traditional multiple wave matching attenuation method. Summary of the invention

[0005] The technical problem to be solved by the present invention is to provide a method for adaptive attenuation of multiple waves in the curvelet domain based on time domain transformation. First, multiple wave records predicted by the wave equation are used to generate multi-order time domain integral and time domain differential records, and the above records are transformed into the curvelet domain by curvelet transformation. Then, multi-channel least square filtering is applied in the domain to eliminate the multiple wave components in the original data, and finally, the multiple wave suppression results are inversely transformed into the time-space domain. The actual data processing results in the field show that the method of the present invention takes into account the adaptability of time domain transformation to the multiple wave prediction error and the separation advantage of curvelet transformation for primary waves and multiple waves, and can effectively suppress the interference of multiple waves, thereby solving the problem of multiple wave attenuation caused by inaccurate multiple wave prediction in complex sea areas.

[0006] The present invention adopts the following technical solutions:

[0007] A method for adaptive attenuation of multiple waves in the curvelet domain based on time domain transformation, the method specifically comprising the following steps:

[0008] The first step is to give an initial velocity model V(x,z), where x and z represent the spatial position coordinates;

[0009] Given a seismic wavelet W(t), t represents time, based on the initial velocity model, set the detection points and shot points, the shot points excite the seismic wavelet, the detection points receive it, and the original seismic record R(x, t) containing multiple waves is obtained. Based on this seismic record, the multiple wave prediction method of the feedback loop is used to obtain the multiple wave record m(x, t);

[0010] Step 2: Input the multiple wave record m(x, t) in the time-space domain and transform it into the frequency domain using the Fourier transform formula (1). The multiple wave record in the frequency domain is M(x, ω), where ω represents the angular frequency and j is the scale of the curvelet. The formula (1) is:

[0011] Step 3: Calculate the series of time-domain integral and time-domain differential records of the multiple wave prediction record M(x,ω) based on formula (2) Where q is the time domain transform order, q∈[q1,q2], when q<0,

[0012] represents the time domain integral record, its main frequency moves toward the low frequency end, and when q>0, For time domain differential recording, the main frequency moves toward the high frequency end; when q = 0, represents the multiple wave record itself; the formula (2) is:

[0013] Step 4: Based on formula (3), the original seismic record R(x, t) and the series of multiple wave records are combined. All are converted to the curvelet domain to obtain the curvelet coefficient c of the original seismic record and the curvelet coefficient of the multiple wave record In formula (3), d(x,t) represents the earthquake record, is the curvelet basis function, j, l, k represent the scale, direction and position of the curvelet respectively; the formula (3) is:

[0014]

[0015] Step 5: Take each sample point in the curvelet domain as the center and cut a rectangular window with a window size of a×b, where a>1 and b>1. Add the data of each row of the window and divide it by the number of horizontal grid points a in the row window to get the average value of the row. Then the cut rectangular window data becomes column data (c1…c b ), use the intercepted data to construct the least square filter formula (4), and use formula (4) to obtain the filter factors f1, f2, f3, ... f b ,

[0016] To intercept the sample points in the rectangular window record block, a is the number of horizontal grid points of the window, and b is the number of vertical grid points of the window; the formula (4) is:

[0017] Step 6: Use formula (5) to realize the adaptive attenuation of multiple waves in the curvelet domain and obtain the suppression result p of a sample point in the curvelet domain. i , where c i is the column data (i=1...b) after averaging in the fifth step, i is the changing longitudinal grid point; the formula (5) is:

[0018] Step 7: Repeat steps (5) and (6) for each sample point in the curvelet domain to obtain the adaptive subtraction result P of the entire record in the curvelet domain;

[0019] Step 8: Based on formula (3), the adaptive subtraction result P of the overall record in the curvelet domain is subjected to a curvelet domain inverse transform to obtain the final multiple wave suppression result in the time-space domain.

[0020] As a preferred embodiment, a is 10 and b is 6.

[0021] The beneficial effects of the present invention compared with the prior art are as follows:

[0022] First, the multiple wave records predicted by the wave equation are used to generate multi-order time-domain integral and time-domain differential records, and the above records are transformed into the curvelet domain by curvelet transform. Then, multi-channel least square filtering is applied in this domain to eliminate the multiple wave components in the original data, and finally the multiple wave suppression results are inversely transformed into the time-space domain. The experimental results of actual field data processing show that the method of the present invention takes into account the adaptability of time domain transformation to the multiple wave prediction error and the advantage of curvelet transform in separating primary waves and multiple waves, and can effectively suppress the interference of multiple waves, thereby solving the problem of multiple wave attenuation caused by inaccurate multiple wave prediction in complex sea areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0023] Figure 1 This is a flow chart of conventional multiple wave matching attenuation;

[0024] Figure 2 It is a flow chart of a method for adaptive attenuation of multiple waves in the curvelet domain based on time domain transformation channel;

[0025] Figure 3 It is a four-layer horizontal velocity model;

[0026] Figure 4 It is the original artillery collection record;

[0027] Figure 5 Multiple wave records for prediction;

[0028] Figure 6The shot gather record after multiple waves are eliminated by conventional time-space domain expansion filtering method and the subtracted multiple wave interference. (a) is the shot gather record after multiple waves are eliminated by conventional time-space domain expansion filtering method (b) is the subtracted multiple wave interference; Figure 7 The shot gather record after eliminating multiple waves by conventional curvelet transform method and the subtracted multiple wave interference; (a) is the shot gather record after eliminating multiple waves by conventional curvelet transform method (b) is the subtracted multiple wave interference;

[0029] Figure 8 The shot gather record after eliminating multiple waves by the adaptive attenuation method of multiple waves in the curvelet domain based on the time domain transformation channel, and the subtracted multiple wave interference. (a) is the shot gather record after eliminating multiple waves by the adaptive attenuation method of multiple waves in the curvelet domain based on the time domain transformation channel, and (b) is the subtracted multiple wave interference. DETAILED DESCRIPTION

[0030] The present invention takes a four-layer horizontal velocity model as an example (such as Figure 3 The horizontal and vertical lengths of the four-layer horizontal velocity model are 1000m and 400m respectively, and the horizontal and vertical grid steps of the model are both 5m.

[0031] In this model, the earthquake source is placed at (0m, 0m), the detection points are located at the top layer of the model, and the interval between the detection points is 5m.

[0032] The specific implementation mode of the present invention is described in detail below. Figure 2 As shown, the specific steps are as follows:

[0033] (1) A model with a horizontal length of 1000 m and a vertical length of 400 m is gridded, and the intervals in the horizontal and vertical directions of the model are both 5 m. Given a seismic wavelet W(t), where t represents time, the shot point position is (0 m, 0 m) to excite the seismic wavelet, the source is excited, and the receiver receives the seismic record to obtain the original seismic record R(x, t) containing multiple waves. Based on this seismic record, the multiple wave record m(x, t) is obtained using the feedback loop multiple wave prediction method.

[0034] (2) Input the multiple wave record m(x, t) in the time and space domain, and transform it into the frequency domain using the Fourier transform formula (1). The multiple wave record in the frequency domain is M(x, ω), where ω represents the angular frequency. The formula (1) is:

[0035]

[0036] (3) Based on formula (2), the series of time-domain integral and time-domain differential records of the multiple wave prediction record M(x,ω) are calculated. Where q is the time domain transformation order, q∈[-50,50], when q<0, represents the time domain integral record, its main frequency moves toward the low frequency end, and when q>0, For time domain differential recording, the main frequency moves toward the high frequency end; when q = 0, represents the multiple wave record itself; the formula (2) is:

[0037]

[0038] (4) Based on formula (3), the original seismic record R(x, t) is combined with the series of multiple wave records All are converted to the curvelet domain to obtain the curvelet coefficient c of the original seismic record and the curvelet coefficient of the multiple wave record In formula (3), d(x,t) is the earthquake record, is the curvelet basis function, j, l, k represent the scale, direction and position in the curvelet domain respectively; the formula (3) is:

[0039] (5) Taking each sample point in the curvelet domain as the center, cut a rectangular window with a window size of a×b, where a>1 and b>1. Add the data of each row of the window and divide it by the number of horizontal grid points a in the window of that row to obtain the average value of that row. Then the cut rectangular window data becomes column data (c1…c b ), use the intercepted data to construct the least square filter formula (4), and use formula (4) to obtain the filter factors f1, f2, f3, ... f b , To intercept the sample points in the rectangular window record block, a is the number of horizontal grid points of the window, and b is the number of vertical grid points of the window; the formula (4) is:

[0040] (6) Using formula (5) to realize the adaptive attenuation of multiple waves in the curvelet domain, the suppression result p of a sample point in the curvelet domain is obtained. i , where c i is the column data (i=1...b) after averaging in the fifth step, i is the changing longitudinal grid point; the formula (5) is:

[0041] (7) Repeat steps (5) and (6) for each sample point in the curvelet domain to obtain the adaptive subtraction result P of the entire record in the curvelet domain;

[0042] (8) Based on formula (3), the adaptive subtraction result P of the overall record in the curvelet domain is subjected to a curvelet domain inverse transformation to obtain the final multiple wave suppression result in the time-space domain.

[0043] In order to prove the effect of the method described in the present invention, the conventional time-space domain expansion filtering method, the conventional curvelet transform method and the method described in this paper are used to Figure 4 The shot gather records shown are processed and the results are shown in Figure 6 (a) Figure 7 (a) and Figure 8 As shown in (a), the subtracted parts are as follows Figure 6 (b) Figure 7 (b) and Figure 8 As shown in (b) in the figure. By comparing the multiple wave attenuation results of the three methods, it can be seen that at the intersection of the primary wave and the multiple wave events (within the dotted circle), due to the close proximity of the events, the time-space domain expansion filtering will inevitably cause effective wave damage at the intersection; while the conventional curvelet transform can transform the events of different scales and directions to different positions in the curvelet domain, reducing the damage of the effective wave, but due to the difference in waveform phase between the predicted multiple waves and the actual multiple wave records, there are obvious multiple wave remnants at the near-track position pointed by the arrow; and the multiple wave matching attenuation method based on the time domain transform and curvelet domain filtering adopted in this paper can still suppress the multiple wave interference well when the multiple wave prediction is inaccurate, and the effective wave event is more continuous. By comparing the processing results, it can be seen that the method in this paper can greatly retain the information of the primary reflection wave while effectively suppressing the multiple waves, and has higher accuracy in suppressing the multiple waves.

Claims

1. A curvelet domain multiple adaptive attenuation method based on time domain transformation channel, characterized in that: The method specifically comprises the following steps: The first step is to give an initial velocity model V(x,z), where x and z represent spatial position coordinates; given a seismic wavelet W(t), where t represents time, set the detection points and shot points based on the initial velocity model, the shot points excite the seismic wavelets, and the detection points receive them, and obtain the original seismic record R(x,t) containing multiple waves. Based on this seismic record, the multiple wave prediction method of the feedback loop is used to obtain the multiple wave record m(x,t); Step 2: Input the multiple wave record m(x, t) in the time-space domain and transform it into the frequency domain using the Fourier transform formula (1). The multiple wave record in the frequency domain is M(x, ω), where ω represents the angular frequency and j is the scale of the curvelet. The formula (1) is: Step 3: Calculate the series of time-domain integral and time-domain differential records of the multiple wave prediction record M(x,ω) based on formula (2) Where q is the time domain transform order, q∈[q1,q2], when q<0, represents the time domain integral record, its main frequency moves toward the low frequency end, and when q>0, For time domain differential recording, the main frequency moves toward the high frequency end; when q = 0, represents the multiple wave record itself; the formula (2) is: Step 4: Based on formula (3), the original seismic record R(x, t) is combined with the series of multiple wave records All are converted to the curvelet domain to obtain the curvelet coefficient c of the original seismic record and the curvelet coefficient of the multiple wave record In formula (3), d(x,t) represents the earthquake record, is the curvelet basis function, j, l, k represent the scale, direction and position of the curvelet respectively; the formula (3) is: Step 5: Take each sample point in the curvelet domain as the center and cut a rectangular window with a window size of a×b, where a>1 and b>1. Add the data of each row of the window and divide it by the number of horizontal grid points a in the row window to get the average value of the row. Then the cut rectangular window data becomes column data (c1…c b ), use the intercepted data to construct the least squares filter formula (4), and use formula (4) to obtain the filter factors f1, f2, f3, ... f b , To intercept the sample points in the rectangular window record block, a is the number of horizontal grid points of the window, and b is the number of vertical grid points of the window; the formula (4) is: Step 6: Use formula (5) to realize the adaptive attenuation of multiple waves in the curvelet domain and obtain the suppression result p of a sample point in the curvelet domain. i , where c i is the column data (i=1...b) after averaging in the fifth step, i is the changing longitudinal grid point; the formula (5) is: Step 7: Repeat steps (5) and (6) for each sample point in the curvelet domain to obtain the adaptive subtraction result P of the entire record in the curvelet domain; Step 8: Based on formula (3), the adaptive subtraction result P of the overall record in the curvelet domain is subjected to a curvelet domain inverse transform to obtain the final multiple wave suppression result in the time-space domain.

2. The method for adaptive attenuation of multiple waves in the curvelet domain based on time domain transformation according to claim 1, characterized in that: The a is 10 and b is 6.

Images