1 Introduction

Because both the sea surface and the seabed are strong reflection interfaces, the marine seismic data often contain many strong multiples, which reduce the signal-to-noise ratio (SNR) of seismic data and lead to poor seismic imaging. Therefore, multiple attenuation is always a key processing step in offshore seismic data processing.

Many methods have been developed for multiple elimination, including but not limited to, the filtering methods (Ryu, 1982; Foster and Mosher, 1992; Zhou and Greenhalgh, 1994; Hu and White, 1998; Song, 2019; Geng et al., 2022; Song et al., 2022) and wave-equation-based methods (Wiggins, 1988; Verschuur et al., 1992; Berkhout and Verschuur, 1997; Verschuur and Berkhout, 1997, 2005; Weglein et al., 1997; Weglein, 1999; Hokstad and Sollie, 2006; van Groenestijn and Verschuur, 2009a, 2009b; Sun and Innanen, 2018, 2019; Zhang and Staring, 2018; Zhang and Slob, 2019, 2020a, 2020b; Aaker et al., 2021a, 2021b). Among all these methods, the surface-related multiple elimination (SRME) method (Verschuur et al., 1992; Berkhout and Verschuur, 1997; Lopez and Verschuur, 2015; Zhang and Verschuur, 2021), one of wave-equation-based methods, is most widely used in practical production.

There are typically two steps involved in the multiple suppression based on the SRME method, including multiple prediction and multiple adaptive attenuation. In the paper, we only focus on the latter.

Usually, there are two types of methods used for adaptive attenuation. One is the Wiener filter method, which is performed in time-space domain, while the other is the matching attenuation method performed in other domains. The Wiener filter method was first developed by Verschuur et al. (1992) for attenuating the multiples generated by the sea surface. Based on the least squares criterion, it enables the SRME method for field data processing. Then, Monk (1991) and Wang (2003a) developed the extended multichannel matching (EMCM) filter method, which significantly enhanced the robustness of the adaptive subtraction. In addition, Li et al. (2007), Li et al. (2010), and Li et al. (2011) proposed pseudo-multichannel least-square matching filter methods which further improve the effectiveness of the multiple attenuation.

The matching attenuation method was first proposed by Zhou and Greenhalgh (1996) based on the parabolic Radon transform. This method applies the parabolic Radon transformation to both the original and the multiple records and then the multiples can be attenuated from the original record by performing a mask filter. Following this, Wang (2003b) achieved adaptive subtraction in time domain by setting the data in Radon domain as the initial model of filtering. To select the appropriate parameters for the matching subtraction, Shi and Wang (2012) designed an effective adaptive Butter-worth-type filter. Herrmann et al. (2008) proposed a curvelet transform-based multiple matching attenuation method, which makes a multi-scale and multi-directional transform. It is often better to sparsely represent the edge detail information of an 'image' such as a seismic event, effectively avoiding the damage of the primary wave. Dong et al. (2015) developed the complex curvelet transform-based multiple adaptive subtraction method.

There are other matching attenuation methods with constraints in addition to the two types mentioned above. For example, the matching subtraction method based on event tracing (Tan et al., 2016, 2017), the model-based adaptive subtraction method (Spitz, 1999), principal component analysis (PCA)-based adaptive subtraction method (Lu et al., 2004), the multi-dimensional predicted error filtering-based adaptive subtraction method (Guitton, 2005), the non-stationary regression-based method (Fomels, 2008), the inverse wavelet domain-based method (Jin et al., 2008), the multi-channel convolutional signal blind separation-based method (Li et al., 2012; Li and Lu, 2013), and the multiple adaptive subtraction method based on machine learning (Li, 2019; Jiang et al., 2020).

In general, the Wiener filtering method can effectively attenuate multiples, but it may damage the primaries when the events of the primaries and the multiples intersect or partially coincide. The matching attenuation methods based on mathematics transformation will better preserve the primary information when event intersections occur, but they have higher requirements for the prediction accuracy of the multiples, making it failure when there are errors in the predicted multiples.

To fully leverage the advantages of these two methods, we combine them together by introducing extended filtering into the curvelet domain, and achieve a multiple matching attenuation method based on curvelet domain extended filtering (hereinafter referred to as the 'curvelet domain extended filtering method'). The method takes into account the adaptability of extended filtering to multiple prediction errors and the advantage of curvelet transform in multiple and primary separation, effectively addressing the problem of multiple attenuation in complex sea areas.

2 Principle of Curvelet Domain Extended Filtering

Monk (1993) first identified several common issues in predicted multiple information, such as amplitude differences, phase (constant phase reversal) differences, and time delay (Tan et al., 2017). As a result, the multiples in the original seismic records can be expressed as follows:

$ d\left(t \right) = {f_1}m\left(t \right) + {f_2}m'\left(t \right) + {f_3}{m^{\text{H}}}\left(t \right) + {f_4}{\left[ {{m^{\text{H}}}\left(t \right)} \right]^\prime }, $

(1)

where t is travel-time, d(t) is the multiple signal of original seismic record, m(t) is the predicted multiple signal, m^H(t) is Hilbert transform of m(t), m'(t) and [m^H(t)]' are the time derivatives of m(t) and m^H(t) respectively, which are both given time window length in the single seismic trace record, f₁, f₂, f₃, and f₄, are the weight factors when $d(t)$ is expanded into the sum of m(t), m^H(t), m'(t) and [m^H(t)]', respectively.

Within a small space and time window, the multiple signals of adjacent traces in a shot gather have obvious 'homology' (the same source, originating from the same underground reflection interface or geological body). Therefore, the same multiple reflection signals in adjacent traces also approximately satisfy the relationship between the predicted multiples and the multiples in the original record described in Eq. (1). As a result, the formula can be extended to 2D as follows:

$\begin{align} d\left({x, t} \right) = & {f_1}m\left({x, t} \right) + {f_2}m'\left({x, t} \right) + {f_3}{m^{\text{H}}}\left({x, t} \right) + \\ & \;{f_4}{\left[ {{m^{\text{H}}}\left({x, t} \right)} \right]^\prime }, \end{align}$

(2)

where x is the offset, d(x, t) is the multiple signal in the original seismic record, m(x, t) is the predicted multiple signal, m^H(t) is Hilbert transform of m(x, t), m'(t) and [m^H(t)]' are the time derivatives of m(x, t) and m^H(t) respectively, f₁, f₂, f₃, and f₄ are the weight factors when d(t) is expanded into the sum of m(x, t), m^H(t), m'(t) and [m^H(t)]', respectively.

Composing the curvelet transform (Candès, 2006) on the original seismic record, we can get:

$ c(j, l, k) = \iint {d(x, t)\overline {{\varphi _{j, l, k}}(x, t)} {\text{d}}x{\text{d}}t}, $

(3)

where $\overline {{\varphi _{j, l, k}}(x, t)} $ denotes the curvelet basis function, c(j, l, k) represents the curvelet coefficients, and j, l and k denote the scale, direction and position of the curvelet, respectively.

Substituting Eq. (2) into Eq. (3) yields:

$\begin{align} c(j, l, k) = &\iint {\left\{ {{f_1}m(x, t) + {f_2}m'(x, t) + {f_3}{m^{\text{H}}}(x, t)} \right.} \\&\; \left. { + {f_4}{{\left[ {{m^{\text{H}}}(x, t)} \right]}^\prime }} \right\}\overline {{\varphi _{j, l, k}}(x, t)} {\text{d}}x{\text{d}}t . \end{align}$

(4)

According to the properties of the integral operation, the integral term at the right end of Eq. (4) can be expanded as:

$\begin{align} c\left({j, l, k} \right) =& {f_1}\iint {m\left({x, t} \right)\overline {{\varphi _{j, l, k}}\left({x, t} \right)} {\text{d}}x{\text{d}}t} + {f_2}\iint {m'\left({x, t} \right)\overline {{\varphi _{j, l, k}}\left({x, t} \right)} {\text{d}}x{\text{d}}t} + {f_3}\iint {{m^{\text{H}}}\left({x, t} \right)\overline {{\varphi _{j, l, k}}\left({x, t} \right)} {\text{d}}x{\text{d}}t} \\&\; + {f_4}\iint {{{\left[ {{m^{\text{H}}}\left({x, t} \right)} \right]}^\prime }\overline {{\varphi _{j, l, k}}\left({x, t} \right)} {\text{d}}x{\text{d}}t} . \end{align}$

(5)

Then, we define:

$ \left\{ \begin{array}{l} {c^M}(j, l, k) = \iint {m\left({x, t} \right)\overline {{\varphi _{j, l, k}}(x, t)} {\text{d}}x{\text{d}}t} \hfill \\ {c^{MD}}(j, l, k) = \iint {m'\left({x, t} \right)\overline {{\varphi _{j, l, k}}(x, t)} {\text{d}}x{\text{d}}t} \hfill \\ {c^{MH}}(j, l, k) = \iint {{m^{\text{H}}}\left({x, t} \right)\overline {{\varphi _{j, l, k}}(x, t)} {\text{d}}x{\text{d}}t} \hfill \\ {c^{MHD}}(j, l, k) = \iint {{{\left[ {{m^{\text{H}}}\left({x, t} \right)} \right]}^\prime }\overline {{\varphi _{j, l, k}}(x, t)} {\text{d}}x{\text{d}}t} \hfill \\ \end{array} \right. . $

(6)

Combining Eqs. (5) and (6) as shown in the formula:

$\begin{align} c(j, l, k) =& {f_1}{c^M}\left({j, l, k} \right) + {f_2}{c^{MD}}\left({j, l, k} \right) + \\& {f_3}{c^{MH}}\left({j, l, k} \right) + {f_4}{c^{MHD}}\left({j, l, k} \right), \end{align}$

(7)

where c(j, l, k) denotes the curvelet coefficients of the multiple signals in the original seismic record, c^M(j, l, k), c^MD(j, l, k), c^MH(j, l, k), and c^MHD(j, l, k) are the curvelet coefficients of m(x, t), m^H(t), m'(t) and [m^H(t)]', respectively, f₁, f₂, f₃, and f₄ denote the weight factors when c(j, l, k) is expanded as the sum of c^M(j, l, k), c^MD(j, l, k), c^MH(j, l, k), and c^MHD(j, l, k).

From Eq. (7), if there are any errors in the predicted multiple signals, the original seismic records can be expanded as the sum of the curvelet coefficients of the predicted signals and their respective transforms (Hilbert transform, time derivative and the time derivative of the Hilbert transform result). In other words, if there is an error in the predicted multiples, the transforms of the predicted record can be introduced into the multiple matching attenuation method in the curvelet domain, and each weight factor f₁, f₂, f₃, and f₄ in Eq. (7) is the corresponding filter factor.

3 Multiple Matching Subtraction Based on Curvelet Domain Extended Filtering

In order to effectively match and attenuate the multiples with prediction errors, this paper integrates the advantages of both curvelet transform and extended matching filtering to achieve multiple matching attenuation based on curvelet domain extended filtering. Firstly, the method generates the corresponding Hilbert transform record m^H(t) using the predicted multiple record m(x, t), and then calculates the time derivatives of the two records m'(x, t) and [m^H(t)]', respectively. Since the multiple adaptive subtraction is performed in the curvelet domain, the original seismic record d(x, t), the predicted multiple record m(x, t) and their transformed records (including m'(x, t), m^H(x, t) and [m^H(x, t)]') need to be converted to the curvelet domain to obtain the curvelet coefficients c(j, l, k), c^M(j, l, k), c^MD(j, l, k), c^MH(j, l, k), and c^MHD(j, l, k) of each record.

In order to ensure the stability of the filtering process, for the sample point at position k in the 2D curvelet coefficients matrix of any scale j and direction l, a rectangular window is set centered on it, and the corresponding data block is intercepted to construct the expression for obtaining the filter factor:

$ \left[ {\begin{array}{*{20}{c}} {c_1^M}&{c_1^{MD}}&{c_1^{MH}}&{c_1^{MHD}} \\ {c_2^M}&{c_2^{MD}}&{c_2^{MH}}&{c_2^{MHD}} \\ {...}&{...}&{...}&{...} \\ {c_n^M}&{c_n^{MD}}&{c_n^{MH}}&{c_n^{MHD}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{f_1}} \\ {{f_2}} \\ {{f_3}} \\ {{f_4}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{c_1}} \\ {{c_2}} \\ {...} \\ {{c_n}} \end{array}} \right], $

(8)

where f₁, f₂, f₃, and f₄ are the filter factors to be obtained, c_n, c_n^M, c_n^MD, c_n^MH, and c_n^MHD are the sample points in the intercepted record block, and n is the number of sample points in the intercepted record block (n ≥ 4 is usually required to enable a stable solution to the equation).

Based on the least squares criterion for the filter factors, Eq. (8) can be transformed into the least-squares problem described by the following equation:

$\boldsymbol{f}=\underset{\boldsymbol{f}}{\arg \min }\|C-\boldsymbol{M} \boldsymbol{f}\|_2,$

(9)

where f is the filter factor vector [f₁, f₂, f₃, f₄]^T, M denotes the matrix consisting of the intercepted record blocks c_n^M, c_n^MD, c_n^MH, and c_n^MHD, and C represents the vector consisting of each sample point c_n in the record block.

Eq. (8) is typically an overdetermined equation, so we solve it by the damped least squares method. The expression for calculating the filter factor f can be expressed as:

$\boldsymbol{f}=\left(\boldsymbol{M}^{\mathrm{T}} \boldsymbol{M}+\varepsilon \boldsymbol{I}\right)^{-1}\left(\boldsymbol{M}^{\mathrm{T}} \boldsymbol{C}\right)$

(10)

where M^T is the transpose matrix of M, I denotes the unit matrix, and ε is the damping coefficient (ε∈(0, 1), set to 0.001 times of the average amplitude of the normalized seismic records).

After obtaining the filter factor f, the multiple matching subtraction process can be realized for the target sample points by the original data curvelet coefficients, that is:

$ p(j, l, k) = c(j, l, k) - {f_1}{c^M}(j, l, k) - {f_2}{c^{MD}}(j, l, k) - {f_3}{c^{MH}}(j, l, k) - {f_4}{c^{MHD}}(j, l, k), $

(11)

where p(j, l, k) denotes the adaptive subtraction result for the target sample points in the curvelet domain.

The curvelet coefficients with multiple elimination are obtained by repeating the processing steps of Eqs. (8) to (11) for the curvelet coefficient matrix C in all scales and directions of the original seismic record. The curvelet coefficients are then transformed back into the time-space domain to obtain the seismic record that multiples have been eliminated. The detailed flowchart of the method is shown in Fig.1, and the basic steps include:

1) Generate the Hilbert transform record m^H(x, t) for the predicted multiple record m(x, t), and then calculate the time derivative records m'(x, t) and [m^H(x, t)]';

2) Convert the original seismic record d (x, t), the predicted multiple record m(x, t) and all transform records (including m'(x, t), m^H(x, t) and [m^H(x, t)]') to the curvelet domain;

3) For each sample point in the matrix of the curvelet coefficients, a rectangular block of data centered on it is intercepted, and the multiple matching attenuation results are calculated by Eqs. (8) to (11);

4) Step (3) is repeated until all sample points in the curvelet coefficients are processed and the final result is transformed back into the time-space domain.

4 Method Validation Based on Synthetic Data

We establish a theoretical model (as shown in Fig.2) to test the multiple elimination effect of the method. A Ricker wavelet with a main frequency of 35 Hz is used as the source to generate a set of synthetic data, including primaries and multiples, based on the model. There are 301 shot records in the dataset and each shot gather has 200 traces. The minimum offset is zero and the trace interval is 10 m. The 150th shot record is shown in Fig.3a. There is only a minimal time difference between the multiple events and primary events, making it challenging to attenuate the multiples without damaging the primaries. Based on the simulated original seismic data, the predicted multiple records are obtained by using the conventional surface-related multiple prediction method. The 150th multiple record is shown in Fig.3b. By comparing the multiples in the original record (as shown by the arrows in Fig.3a) with the predicted multiples (as shown by the arrows in Fig.3b), we can find significant differences in both waveform and phase between them.

Based on the original shot gather and the predicted multiple records, we conducted the adaptive subtraction using the time-space domain extended filtering method, the traditional curvelet transform method and the curvelet domain extended filtering method to evaluate the effectiveness of multiple attenuation t. The results are shown in Figs.4a, 4b, and 4c, while the subtracted multiple records are shown in Figs.5a, 5b, and 5c, respectively.

From Figs.4 and 5, it is evident that all three methods effectively attenuate the predicted multiples. However, significant damages to the primary waves occurs at the location where the primary event intersects the multiple event (as marked by the circled ranges in Fig.4 and the arrows in Fig.5) for the time-domain extended filtering method and the traditional curvelet transform method. However, there are no apparent damages to the primary waves in the result of curvelet domain extended filtering method, which demonstrates that the method has higher accuracy in eliminating multiples when multiple events and primary events are very close. In summary, the method proposed in this paper introduces the transform records of the predicted multiples to better match the multiple signals in the curvelet domain, which not only improves the matching attenuation effect of the multiple but also effectively avoids the damage to the primary records at the location of intersection (or partial coincidence) of the primary events and the multiple events.

5 Method Validation Based on the Field Data

A 2D field seismic dataset is used to test the applicability of the method to the field data processing. There are a total of 2362 shots, each with 180 traces. The shot and trace interval are both 26.6 m, and the minimum offset is 218.8 m. Because the seismic data was acquired from a sea area with rigid and rough seabed, there are plenty of strong surface-related multiples and other complex multiples in the data.

Before attenuating surface-related multiples, several data pre-processing steps are performed as follows: 1) the attenuation for the direct and refracted interference waves in the shallow part; 2) high-pass filtering with a high cut-off frequency of 5 Hz; 3) the spherical spreading compensation.

The records before multiple attenuation and the predicted multiple records of the 201st shot are shown in Figs. 6a and 6b. By comparing the rectangular areas, we can find some differences.

Figs.7a–7c display the results of multiple attenuation by using the time-space domain extended filtering method, the traditional curvelet transform method and the curvelet domain extended filtering method, respectively. Although most of the strong multiples are suppressed by the previous two methods, there are still some residual multiples visible in Figs.7a and 7b. However, there are no obvious residual multiples in Fig.7c.

Figs.8a–8c display the multiples removed by three methods mentioned above. Not only are there more removed multiples in Fig.8c compared to Figs.8a and 8b, but there are also fewer primary events remaining mistakenly. This proves that the curvelet domain extended filtering method, which combines the advantages of extended filtering and curvelet transform, is effective in suppressing the multiples in field data while preserving the primary waves.

We performed Kirchhoff prestack time migration based on the data before multiple attenuation and the data after multiple attenuation by using the three methods. The results are shown in Figs.9–12.

Due to the application of the primary migration velocities, even based on the data before multiple attenuation, only a few strong multiples are imaged in the section (as indicated by the arrows in Fig.9). In Figs.10 and 11, there are still some residual multiples imaged, which implies that there must be lots of residual multiples in the data after multiple elimination by the time-space domain extended filtering method or the traditional curvelet transform method. However, no obvious multiple events can be seen in Fig.12, which further proves that the method proposed in the paper is capable of improving the effect of the adaptive subtraction under the premise of preserving primary reflections, even for the areas with rough seabed or complex geological structures.

6 Conclusions

The paper introduces a multiple matching attenuation method based on curvelet domain extended filtering, which combines the extended Wiener filter method and the curvelet transform method, drawing on their respective strengths to effectively suppress the multiples without damaging the primaries. Tests on the model data and field data show that the proposed method has higher accuracy in multiple elimination and is more suitable for the multiple attenuation in areas with complex structures, compared to the time-space domain extended filtering method and the conventional curvelet transform method.

However, there are also some issues with the curvelet domain extended filtering method, notably its low computational efficiency. While the traditional curvelet transform method requires only 2 curvelet transforms for calculation, the curvelet domain extended filtering method necessitates 5 curvelet transform calculations. Therefore, our main focus moving forward is to improve the computational efficiency of the curvelet domain extended filtering method.

Acknowledgements

This research is jointly funded by the Wenhai Program of the ST Fund of Laoshan Laboratory (No. 202204803), the National Natural Science Foundation of China (Nos. 4207 4138, 42206195), the National Key R&D Program of China (No. 2022YFC2803501), and the Research Project of the China National Petroleum Corporation (No. 2021ZG02).