Pseudo-3D Receiver Deghosting of Seismic Streamer Data Based on l 1 Norm Sparse Inversion

Featured Application: Authors are encouraged to provide a concise description of the speciﬁc application or a potential application of the work. This section is not mandatory. Abstract: The ghost effect in marine seismic data causes low-frequency suppression and frequency notch, resulting in incomplete frequency information for seismic records, which poses challenges for high-resolution imaging. The deghosting effect depends on the approximation of the ghost delay operator. Due to the strict requirements of dense sampling, the 2D deghosting method for a densely sampled inline dataset is still the mainstream method for marine data processing. As the trade-off, inversion-based methods are widely used in the industry to reduce the inﬂuence of the inaccurate ghost delay operator. In order to overcome the sampling limits and improve the 2D deghosting effect, we propose a pseudo-3D deghosting framework based on an l 1 norm sparse inversion. In the framework, the inversion problem is divided into two subproblems, i.e., pseudo-3D operator building and optimization inversion. Considering the data coherence along the shot direction, we derive a pseudo-3D ghost delay operator to deghost simultaneously for multi-shot gathers. We then introduce a sparse inversion method in the pseudo-3D radon domain (multi-shot gathers) to further improve the inversion accuracy. The proposed framework is easy to implement, is not sensitive to noise, and shows superior performance in synthetic and ﬁeld examples.


Introduction
For signal-to-noise ratio considerations, the receivers are sunk to a certain depth below the air-water interface during marine seismic acquisition [1]. When the upgoing wavefield propagates to the air-water interface, it reflects and propagates downward. The downgoing wavefield is called ghost [2], and it interferes with the upgoing wavefield, resulting in frequency notch and reducing the vertical resolution of the seismic dataset [3].
In marine acquisition, the sampling interval in the crossline direction is usually much larger than that in the inline direction (usually between 4:1 and 16:1) [4], and the sparse sampling poses a large challenge for deghosting [5]. Therefore, the 2D deghosting method for a densely sampled inline dataset is still the mainstream method for marine data processing. The current deghosting methods are mainly implemented in the frequency domain, and the idea of plane wave decomposition is directly or indirectly applied [6]. The deghosting methods based on plane wave decomposition can be roughly divided into two categories: the deconvolution-based method and the inversion-based method.
The basic idea of the deconvolution-based method is to transform the seismic dataset containing ghost into other data domains (time domain [7], frequency domain [8], f -k domain [9], etc.), construct a ghost delay operator, use the deconvolution algorithm to achieve deghosting in this domain, and transform it into the space-time domain. However,  , ) The conventional 2D deghosting method indicates a shot-by-shot deghosting manner. Because every shot gather is processed independently, the inversion accuracy of the 2D deghosting algorithm is easily limited by the field noise and accuracy of the 2D operator, which is necessary to introduce a priori information to constrain the inversion. During the horizontal streamer acquisition, the location of the sources and receivers is relatively stable, and so the offset parameters of the different shots are consistent. According to this property, we can rearrange the dataset into the data cube shown in Figure 2. In this data cube, the seismic records are arranged in a source-offset-time manner. Then, we can adopt a Fourier transform to the data cube, and the frequency slice pse-3D ( ) can be expressed as: The delay time t ip in the radon domain is calculated by To simplify Equation (1), we redefine the 2D deghosting operator G in the radon domain, which can be regarded as a combination of the radon transform operator and the ghosting operator as follows: Equation (1) can then be rewritten as a form of matrix multiplication: The conventional 2D deghosting method indicates a shot-by-shot deghosting manner. Because every shot gather is processed independently, the inversion accuracy of the 2D deghosting algorithm is easily limited by the field noise and accuracy of the 2D operator, which is necessary to introduce a priori information to constrain the inversion. During the horizontal streamer acquisition, the location of the sources and receivers is relatively stable, and so the offset parameters of the different shots are consistent. According to this property, we can rearrange the dataset into the data cube shown in Figure 2. In this data cube, the seismic records are arranged in a source-offset-time manner. Then, we can adopt a Fourier transform to the data cube, and the frequency slice D pse−3D (ω k ) can be expressed as: · · · · · · · · · · · · · · · D ns (x 1 , ω k ) · · · D ns (x j , ω k ) · · · D ns (x ntr , ω k ) where ns represents the number of common shot gathers (CSGs) and ntr represents the number of common offset gathers (COGs). Since the offsets and receiver depths of the  (1) is applicable to the data slice described in Figure 2: where s i and x j represent the location of the source and offset, respectively. As shown in Equation (7), the deghosting case for multi-shot gathers is integrated into a linear system of equations, which can be simplified as: where G represents the ghost delay operator. In this way, we can simultaneously deghost for the data cube. In the source direction of the data cube, each slice represents a COG. The characteristics of events in COGs are similar to those of the underground interface, which means the COGs are composed as several smooth and approximately linear events [26]. The event coherence in COGs provides the basis for the introduction of sparse representations, and we can utilize the data coherence along the shot direction as an inversion constraint. In the common offset gathers, more information that can be used for coherence analysis is obtained [27], which effectively improves the inversion accuracy.
where ns represents the number of common shot gathers (CSGs) and ntr represent the number of common offset gathers (COGs). Since the offsets and receiver depths of th different shot gathers data are consistent, the ghost delay operator in Equation (1) is ap plicable to the data slice described in Figure 2: where i s and j x represent the location of the source and offset, respectively. As show in Equation (7), the deghosting case for multi-shot gathers is integrated into a linear sys tem of equations, which can be simplified as: where G represents the ghost delay operator. In this way, we can simultaneousl deghost for the data cube. In the source direction of the data cube, each slice represents COG. The characteristics of events in COGs are similar to those of the underground inter face, which means the COGs are composed as several smooth and approximately linea events [26]. The event coherence in COGs provides the basis for the introduction of spars representations, and we can utilize the data coherence along the shot direction as an in version constraint. In the common offset gathers, more information that can be used fo coherence analysis is obtained [27], which effectively improves the inversion accuracy.

Compressed Sensing-Based Sparse Inversion Method
In recent years, compressed sensing-based inversion methods, which improve inver sion accuracy by sparse representation, have been emerging. The basic idea is to discretiz the seismic signal by a sampling matrix, then remove redundant information and retai fewer elements to represent the main characteristics of the signal. Because the signal be comes sparse, it does not easily fall into the local minimum during inversion. In this way the signal becomes sparse and redundant information is not included during inversion and it does not easily fall into the non-information local minimum during inversion [28] Sparse inversion requires the sparse representation ( Figure 3) of the seismic signal m y R ∈ using the sampling matrix . The sparse representation can be describe as a linear equation:

Compressed Sensing-Based Sparse Inversion Method
In recent years, compressed sensing-based inversion methods, which improve inversion accuracy by sparse representation, have been emerging. The basic idea is to discretize the seismic signal by a sampling matrix, then remove redundant information and retain fewer elements to represent the main characteristics of the signal. Because the signal becomes sparse, it does not easily fall into the local minimum during inversion. In this way, the signal becomes sparse and redundant information is not included during inversion, and it does not easily fall into the non-information local minimum during inversion [28].
Sparse inversion requires the sparse representation ( Figure 3) of the seismic signals y ∈ R m using the sampling matrix A ∈ R n×m . The sparse representation can be described as a linear equation: where the sampling matrix A can be expressed as where x ∈ R m represents the coefficients in the transformed domain and in x 0 m, x 0 represents the l 0 norm of x, which measures the number of non-zero elements and represents the degree of matrix sparsity, and Φ y represents the vector of the sampling matrix A. Because the l 0 norm is difficult to solve optimally (non-deterministic polynomial-hard problem, NP-hard), an l 1 norm is usually used to replace the l 0 norm as the constraint condition to simplify the optimization problem, as follows: where m x R ∈ represents the coefficients in the transformed domain and in 0 x m  , 0 x represents the 0 l norm of x , which measures the number of non-zero elements and represents the degree of matrix sparsity, and y Φ represents the vector of the sampling matrix A . Because the l0 norm is difficult to solve optimally (non-deterministic polynomial-hard problem, NP-hard), an l1 norm is usually used to replace the l0 norm as the constraint condition to simplify the optimization problem, as follows: To demonstrate the sparsity change before and after the sparse representation of the seismic signal, we adopt a curvelet transform [29] to noisy synthetic data for sparse representation, where the noise distribution is . We select 5% as the curvelet and time-space coefficients, in order of amplitude, to reconstruct the data. From Figure 4, the sparsity of the curvelet coefficients is significantly better than that of the time-space coefficients. In the optimization inversion process, the objective function is generally minimized in the sense of least squares. However, the seismic signals not only contain useful signals, but also noise. The sparse representation of seismic signals regards the noise as redundant information, which means inversion in the transform domain can obtain more stable inversion results. Therefore, the sparse inversion methods are suitable for solving such problems. To demonstrate the sparsity change before and after the sparse representation of the seismic signal, we adopt a curvelet transform [29] to noisy synthetic data for sparse representation, where the noise distribution is N(µ, σ) =N(0, 0.0003) . We select 5% as the curvelet and time-space coefficients, in order of amplitude, to reconstruct the data. From Figure 4, the sparsity of the curvelet coefficients is significantly better than that of the time-space coefficients.
where m x R ∈ represents the coefficients in the transformed domain and in 0 x m  , 0 x represents the 0 l norm of x , which measures the number of non-zero elements and represents the degree of matrix sparsity, and y Φ represents the vector of the sampling matrix A . Because the l0 norm is difficult to solve optimally (non-deterministic polynomial-hard problem, NP-hard), an l1 norm is usually used to replace the l0 norm as the constraint condition to simplify the optimization problem, as follows: To demonstrate the sparsity change before and after the sparse representation of the seismic signal, we adopt a curvelet transform [29] to noisy synthetic data for sparse representation, where the noise distribution is . We select 5% as the curvelet and time-space coefficients, in order of amplitude, to reconstruct the data. From Figure 4, the sparsity of the curvelet coefficients is significantly better than that of the time-space coefficients. In the optimization inversion process, the objective function is generally minimized in the sense of least squares. However, the seismic signals not only contain useful signals, but also noise. The sparse representation of seismic signals regards the noise as redundant information, which means inversion in the transform domain can obtain more stable inversion results. Therefore, the sparse inversion methods are suitable for solving such problems. In the optimization inversion process, the objective function is generally minimized in the sense of least squares. However, the seismic signals not only contain useful signals, but also noise. The sparse representation of seismic signals regards the noise as redundant information, which means inversion in the transform domain can obtain more stable inversion results. Therefore, the sparse inversion methods are suitable for solving such problems.

Pseudo-3D Deghosting Based on Sparse Inversion
To highlight the data coherence, we introduce the compressed sensing-based sparse inversion methods into the deghosting method. In the pseudo-3D deghosting framework, the linear system of Equation (6) can be rewritten as where m represents the vector of the deghosting result M pse−3D in the radon domain, d represents the vector of the input dataset D pse−3D in the time-space domain, A represents the identity matrix, ⊗ represents the Kronecker product of the matrix, f and f H t represent the Fourier transform pair in the time direction, and L represents the combination of the ghost delay operator and the Fourier transform operator.
In actual seismic acquisition, field noise cannot be avoided. Assuming that the input datasets d are a collection of the noise n and useful signals s = Lm, then: Based on the compressed sensing theory, the seismic signals are over-complete and satisfy the requirement of sparse representation. Therefore, the useful signal can be reconstructed by the sampling matrix A: where x 0 represents the coefficients in the transformed domain. In general, A can be constructed by a sparse transform (wavelet transform [30], curvelet transform [29], etc.). Without considering the noise, the main purpose of this problem is to solve the sparse solution of LAx = d, where LA represents a matrix of the size m × n and d represents a vector of the size n × 1. When m n, this problem is under-determined, which can be summarized as a basic pursuit problem (BP). To solve this BP problem [31], we propose the objective function: min Considering the existence of noise, the BP problem can be transformed into a BPDN problem (basis pursuit denoising) [32] to avoid the noise interference during inversion: where σ represents the l 2 norm residual of the objective function, which is generally set as 5-10% of the l 2 norm of d. The proposed method uses the l 2 norm to approximate the reconstructed data to the original data while the l 1 norm constrains the noise. In this way, the noise n in the data can be suppressed. To solve the above problem, we adopt the spectral gradient projection method (SPGL1) [33] min wherem andx represent the inversion results in the data domain and the sparse domain, respectively, σ represents the noise level, which controls the sparsity, and SLGL1 is a solver for the l 1 -regularized least-squares inversion problem, which is suitable for largescale datasets. In order to improve the inversion accuracy, we adopt a sparse transform to achieve the sparse representation of the curved events with multi-scale and multi-directional characteristics. It can be seen from Section 2.2 that sparse representation can increase the sparseness of the signal. Considering representation accuracy and computational efficiency, we combine a 2D curvelet transform and a 1D wavelet transform as the sampling matrix. Therefore, the pseudo-3D deghosting based on a sparse inversion can be expressed as where x represents the 3D coefficients in the sparse domain, S represents the combined operators of the 2D curvelet transform and the 1D wavelet transform, and S = C 2 ⊗ W 1 and S H are the inverse transforms of S. The inversion result is the coefficients in the radon domain, and then we obtain the result in the time-space domain by adopting a conventional inverse radon transform.

Results
In order to verify the effectiveness of the proposed method for deghosting, we first apply the pseudo-3D algorithm to the synthetic data obtained from the forward modeling of the saltdome model ( Figure 5). The saltdome model is a classic acoustic model for testing multiple attenuation methods for marine data [34], and it has a width of 5400 m, a depth of 1250 m, and a grid spacing of 5 m. The model contains high-speed layers and faults, and the corresponding records are complex. As for the observation system, in order to highlight the influence of ghost, the sinking depth of the sources is 14 m and the sinking depth of the receivers is 25 m. The maximum offset is 2550 m and the trace interval is 10 m. To further verify the effectiveness of the algorithm, we also adopt the algorithm to the field streamer data. For data preparation, we choose the open-source software Madagascar, which is applicable to digital image and data processing in geophysics and related fields, as a tool for forward modeling and data preprocessing.
 where x represents the 3D coefficients in the sparse domain, S represents the combined operators of the 2D curvelet transform and the 1D wavelet transform, and 2 1 = ⊗ S C W and H S are the inverse transforms of S . The inversion result is the coefficients in the radon domain, and then we obtain the result in the time-space domain by adopting a conventional inverse radon transform.

Results
In order to verify the effectiveness of the proposed method for deghosting, we first apply the pseudo-3D algorithm to the synthetic data obtained from the forward modeling of the saltdome model ( Figure 5). The saltdome model is a classic acoustic model for testing multiple attenuation methods for marine data [34], and it has a width of 5400 m, a depth of 1250 m, and a grid spacing of 5 m. The model contains high-speed layers and faults, and the corresponding records are complex. As for the observation system, in order to highlight the influence of ghost, the sinking depth of the sources is 14 m and the sinking depth of the receivers is 25 m. The maximum offset is 2550 m and the trace interval is 10 m. To further verify the effectiveness of the algorithm, we also adopt the algorithm to the field streamer data. For data preparation, we choose the open-source software Madagascar, which is applicable to digital image and data processing in geophysics and related fields, as a tool for forward modeling and data preprocessing.

Synthetic Example
To ensure the completeness of the synthetic data, we choose the wave equation-based simulation method, i.e., the fine difference algorithm, to conduct forward modeling. The synthetic multi-shot synthetic records include 256 CSGs, 256 COGs, 512 samplings, and a sampling interval of 0.004 s. Then, we apply the conventional 2D method and the sparse inversion method to obtain the deghosting results, and we compare the results of the two methods from the CSGs, COGs, f-k spectrum, etc. Figure 6 shows the original data before deghosting. Figure 6a-c shows the data cube, the zoom section of the COG, and the corresponding f-k spectrum of the CSG. From Figure  6a, we can observe that when the sinking depth is large, ghost is easy to recognize. The separation of the primaries and the ghosts can be clearly observed in the near-offset section. Figure 6b shows that the frequency bandwidth of the deghost-interfered datasets is

Synthetic Example
To ensure the completeness of the synthetic data, we choose the wave equation-based simulation method, i.e., the fine difference algorithm, to conduct forward modeling. The synthetic multi-shot synthetic records include 256 CSGs, 256 COGs, 512 samplings, and a sampling interval of 0.004 s. Then, we apply the conventional 2D method and the sparse inversion method to obtain the deghosting results, and we compare the results of the two methods from the CSGs, COGs, f -k spectrum, etc. Figure 6 shows the original data before deghosting. Figure 6a-c shows the data cube, the zoom section of the COG, and the corresponding f -k spectrum of the CSG. From Figure 6a, we can observe that when the sinking depth is large, ghost is easy to recognize. The separation of the primaries and the ghosts can be clearly observed in the near-offset section. Figure 6b shows that the frequency bandwidth of the deghost-interfered datasets is narrowed, with a lack of low frequency (5-10 Hz) and a main frequency notch . Therefore, in order to obtain a frequency-rich seismic record, it is necessary to deghost. Figures 7 and 8 respectively demonstrate the deghosting results of the conventional 2D algorithm and the proposed method. Comparing Figures 7a and 8a, it can be seen that both methods have suppressing effects on deghosting. However, the conventional algorithm has some artifacts in the deep position (blue arrow). The reason is that the complexity of the underground structure increases and the inversion process is unstable. These artifacts in the sparse inversion result are effectively suppressed, which proves that the introduction of the proposed method improves the deghosting effect. From Figures 7b and 8b, we see that the notch frequency is mainly 25-50 Hz. Both results compensate for the notch effect, but some energy leakage can still be observed in the conventional results. The sparse inversion results have the complete frequency band and more balanced high-frequency energy, which proves the effectiveness of the proposed method. Figures 7c and 8c show the zoomed section of the COGs. The difference in the deep events (1.5 s, 150-250th shots) is obvious. The conventional method has artifacts but the sparse inversion results do not. narrowed, with a lack of low frequency (5-10 Hz) and a main frequency notch . Therefore, in order to obtain a frequency-rich seismic record, it is necessary to deghost.  Comparing Figures 7a and 8a, it can be seen that both methods have suppressing effects on deghosting. However, the conventional algorithm has some artifacts in the deep position (blue arrow). The reason is that the complexity of the underground structure increases and the inversion process is unstable. These artifacts in the sparse inversion result are effectively suppressed, which proves that the introduction of the proposed method improves the deghosting effect. From Figures 7b  and 8b, we see that the notch frequency is mainly 25-50 Hz. Both results compensate for the notch effect, but some energy leakage can still be observed in the conventional results. The sparse inversion results have the complete frequency band and more balanced highfrequency energy, which proves the effectiveness of the proposed method. Figures 7c and  8c show the zoomed section of the COGs. The difference in the deep events (1.5 s, 150-250th shots) is obvious. The conventional method has artifacts but the sparse inversion results do not.
To further demonstrate the effectiveness of the proposed method, we calculate the frequency spectrum in the time window (1.25-1.5 s) to quantitatively compare the frequency compensation ability of the conventional deghosting method and the pseudo-3D deghosting method (Figure 9). Figure 9b clearly shows the energy comparison of the three curves. Compared with the f-k spectrum shown above, we can observe a good correspondence between them. It is shown in Figure 6b that the frequency bandwidth of deghostinterfered datasets is narrowed at the main frequency notch , and in Figure 9b, the original data has the lowest energy (near 0.25) at approximately 40 Hz. Both the red curve and the yellow curve in Figure 9b have higher energy than the original data in the main frequency notch , which is consistent with the phenomenon seen in Figures 7b and 8b, where both methods compensate for the notch effect. Similarly, the energy     To further demonstrate the effectiveness of the proposed method, we calculate the frequency spectrum in the time window (1.25-1.5 s) to quantitatively compare the frequency compensation ability of the conventional deghosting method and the pseudo-3D deghosting method ( Figure 9). Figure 9b clearly shows the energy comparison of the three curves.
Compared with the f -k spectrum shown above, we can observe a good correspondence between them. It is shown in Figure 6b that the frequency bandwidth of deghost-interfered datasets is narrowed at the main frequency notch , and in Figure 9b, the original data has the lowest energy (near 0.25) at approximately 40 Hz. Both the red curve and the yellow curve in Figure 9b have higher energy than the original data in the main frequency notch , which is consistent with the phenomenon seen in Figures 7b and 8b, where both methods compensate for the notch effect. Similarly, the energy leakage caused by traditional methods also leads to a large difference between the 2D method and the pseudo-3D method in a quantitative comparison of the main frequency energy. The results of the quantitative comparison also prove the effectiveness of the proposed method.

Field Example
In this section, we adopt a field dataset to test the deghosting effect. The dataset contains a total of 193 shot gathers, 128 offset gathers, 1000 sampling points, and a sampling interval of 4 ms. The sink depth of the geophone is 7 m. The maximum spacing is 3175 m, and the track spacing is 25 m. Figures 10-12 show the original data, the conventional 2D

Field Example
In this section, we adopt a field dataset to test the deghosting effect. The dataset contains a total of 193 shot gathers, 128 offset gathers, 1000 sampling points, and a sampling interval of 4 ms. The sink depth of the geophone is 7 m. The maximum spacing is 3175 m, and the track spacing is 25 m. Figures 10-12 show the original data, the conventional 2D method, and the deghosting result of the algorithm, respectively.   As a result, there are three main differences: frequency extension, event continuity, and signal-to-noise ratio (SNR). In terms of frequency extension, since the source depth of the data is 7 m, the first non-zero notch frequency is relatively large. The influence of ghost on the seismic records is mainly manifested in the suppression of the low frequency (2-10 Hz) and high frequency (60-80 Hz). The two methods have similar effects on lowfrequency energy compensation, but the proposed method has better effects on highfrequency energy compensation. The continuity of the events (Figures 11c and 12c) also proves the effectiveness of the proposed method. Because noise is unavoidable in the field datasets, the continuity of the event becomes worse. The proposed method uses the sparse representation of the seismic signal and retains the main components for reconstruction to improve the event continuity. The SNR difference can be clearly observed in the f -k spectrum, which means the sparse inversion method has a stronger anti-noise performance. The noise reduces the image quality of the f -k spectrum, and the f -k spectrum of the sparse inversion result is smoother, indicating that the redundant information (noise) is attenuated during the inversion process. In summary, the deghosting method in this paper has more advantages in frequency extension. In addition, this method is more suitable for field data processing, which can simultaneously improve the SNR in the process of deghosting.

Discussion
In this paper, we introduce the l1 norm sparse inversion method to improve the version accuracy. The curvelet and wavelet transforms are excellent tools for the spa representation of seismic signals, which is fully data-driven. The core assumption of sparse representation is that the noise energy is smaller than the signal energy in the tra formed domain, otherwise signal damage may occur during the inversion. For comp geological conditions, the wavefield contains some weak signals generated by the sma

Discussion
In this paper, we introduce the l 1 norm sparse inversion method to improve the inversion accuracy. The curvelet and wavelet transforms are excellent tools for the sparse representation of seismic signals, which is fully data-driven. The core assumption of the sparse representation is that the noise energy is smaller than the signal energy in the transformed domain, otherwise signal damage may occur during the inversion. For complex geological conditions, the wavefield contains some weak signals generated by the small-scale structure, which more easily covered by noise. In this case, the parameter setting criteria of the proposed method is a trade-off between being signal-preserving and noise-free.
In addition, although sparse representation improves the sparsity of seismic signals, the redundancy is also improved, which means an increase in data size. Therefore, the calculation efficiency of the proposed deghosting method is limited. To further improve the proposed method, we recommend two future research directions: (1) Based on symmetry arguments, one would expect that they could also be applied for source-side deghosting in the common-receiver domain. (2) A combination of this method with a multiple suppression method to propose a set of pre-processing flow data suitable for offshore exploration, including degusting and demultiple.

Conclusions
In this paper, we propose a novel pseudo-3D deghosting framework based on 3D sparse inversion to improve the effect of the conventional 2D deghosting method. The pseudo-3D data cube required by the framework is easily implemented from conventional 2D streamer data. To improve the inversion accuracy, the data coherence along the shot direction and the compressed sensing-based inversion method are introduced.
The synthetic and field examples show that the main advantages of the proposed method are the frequency compensation and noise stability. Field noise cannot be avoided in marine seismic acquisition, and the processing process requires the amplitude-preservation of useful events. The algorithms with low noise sensitivity are more suitable for deghosting in industrial data processing.
Author Contributions: D.W. and R.W. conceived and designed the study; W.Z. programmed the code of the deghosting algorithm; Y.L. evaluated and interpreted the results; B.H. provided the idea; L.W. processed the synthetic and field datasets; and R.W. wrote the paper. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The seismic data used in this paper are the property of Jilin University.