**1. Introduction**

Multiple reflections are usually regarded as coherent noise in reflection seismic data especially in marine exploration. The existence of multiple reflections brings difficulty to seismic data processing [1]. It makes the structural illusion appear in the migrated-stacked section, which can affect the accuracy of subsequent data interpretation [2]. The Radon transform is a common method for multiple attenuation because of its high efficiency and due that it is easily implemented [3,4]. Due to the kinematic differences between primary and multiple reflections, most multiples are considered as coherent noise with low velocity [5]. Therefore, the primary reflections are upturned or flat, and multiple reflections are parabolic after normal movement correction with certain velocities in the CMP gather, which provides a theoretical basis for multiple attenuation using Radon transform. In 1986, Hampson [6] used the parabolic Radon transform to suppress multiples. The parabolic Radon transform makes their seismic events of different curvatures projected in different regions of the Radon domain, so as to separate primary and multiple reflections. However, finite seismic acquisition aperture of seismic data causes energy diffusion in the Radon domain, and it causes the smearing shown by the red arrows in Figure 1b. The energy diffusion limits the effect of multiple attenuation. To improve the above problem, the conception of optimization inversion is introduced into Radon transform [7,8].

When solving the inverse problems of Radon transform, the solution of inversion Radon transform is not unique. That is to say, there will be multiple sets of different solutions in the Radon domain, and these solutions do not describe the model function well.

Hampson [6] proposed the least squares inversion method, and the original data are reconstructed by optimizing the data in the Radon domain. This method reduces the smearing in the radon domain, but there are still some artefacts and data outside the

**Citation:** Wu, Q.; Hu, B.; Liu, C.; Zhang, J. Sparse Parabolic Radon Transform with Nonconvex Mixed Regularization for Multiple Attenuation. *Appl. Sci.* **2023**, *13*, 2550. https://doi.org/10.3390/ app13042550

Academic Editor: Roberto Zivieri

Received: 12 January 2023 Revised: 10 February 2023 Accepted: 14 February 2023 Published: 16 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

offset that are not restricted. After Radon inverse transform, the data within a restricted offset still consist of the original data, but they are different when the data are outside the restricted range. In addition, the solution of the Radon domain is not optimal. Therefore, to solve this problem, it is hoped that the solution of the Radon domain will be sparser [2]. In 1985, the Radon transform was considered as a sparse inversion problem by Thorson et al. [9]. However, a lot of computation is disadvantageous. Sacchi et al. [10] adopted the method of Thorson to implement the sparse inversion Radon transform, and it is in the frequency domain. A high-resolution Radon transform was proposed by Herrmann to distinguish primary and multiple reflections, and it tackled the aliasing and resolution issues of the Radon transform [11]. Trad et al. [12] discussed fast implementations for the Radon transform in the time and frequency domains. With this sparse constraint, the seismic events have better localization characteristics in the Radon domain and the separation effect of primaries and multiples is guaranteed to a certain extent. Lu [13] proposed iterative 2D model shrinkage to solve the sparse inverse problem of RT and obtained the analogously sparse Radon model. Xiong et al. [14] established a mathematical model to combine *L*2-norm and the adaptive multiple subtraction method on *L*1-norm by weighted combination. This method can suppress the energy of multiples relatively better than non-combination.

**Figure 1.** (**a**) CMP gather after NMO; (**b**) Radon domain data after parabolic Radon transform. The red arrows indicate the smearing.

In order to further improve the separation effect, researchers use the regularization methods of L-norm as the sparsity constraint condition that fully reflects the basic characteristics of the seismic data [15–17]. We can obtain the high-resolution Radon transform, and it is constrained sparsely, such as in *L*<sup>0</sup> regularization, but it is generally not used because it is hard to solve. Tisbshirani [18] proposed *L*<sup>1</sup> regularization, and it provides an alternative. *L*<sup>1</sup> regularization [19] is a convex optimization problem. This means that *L*1-norm regularization is easy to solve. Donoho et al. [20] proved that sometimes the solution of the *L*<sup>0</sup> regularization is equivalent to that of the *L*<sup>1</sup> regularization for the sparsity problem. *L*<sup>2</sup> regularization can avoid the overfitting of the model, but its solutions do not have the sparse property. Moreover, some numerical experiments showed that sparse signals are recovered from fewer linear measurements by *Lq* (0 < *q* < 1) regularization [21–24]. *L*1/2 regularization can be taken as a representative of the *Lq* (0 < *q* < 1) regularization, and it has been proved that the solution of *L*1/2 regularization is sparser and more stable than the

solution of *L*<sup>1</sup> regularization [25]. However, the *L*1/2 regularization leads to a nonconvex, nonsmooth and non-Lipschitz optimization problem that is difficult to solve quickly and efficiently. Some scholars proposed corresponding solving methods. The reweighted iteration algorithm [25], the iterative reweighted least squares method [26], the iterative half thresholding algorithm [27] and the generalized iterated shrinkage algorithm [28,29] are efficient methods for solving the *L*1/2 regularization.

Generally, sparse regularization constrains the sparsity of the whole seismic wavefield, which means a unique regularization parameter for different wavefields. However, the different wavefields have different amplitudes, and the same regularization parameter setting easily causes signal damage or redundant residual information. Sparse regularization for multi-tasking exploits differences in features of the data to demix the two distinct components for constraining the sparsity. For the multiple attenuation case in the Radon domain, primary and multiple reflections are distributed in different zones of the Radon domain and have different amplitudes. The unique regularization parameter may limit the performance of multiple attenuation.

In this paper, we propose a sparse parabolic Radon transform in the frequency domain with the nonconvex *Lq*<sup>1</sup> − *Lq*<sup>2</sup> (0 < *q*1, *q*<sup>2</sup> < 1) mixed regularization (SPRT*Lq*<sup>1</sup> − *Lq*<sup>2</sup> ) that constrains the sparsity of primary and multiple reflections, respectively. In addition, we use the alternating direction method of multipliers algorithm (ADMM) [30–32] to approximately solve the multi-task regularization problem. Furthermore, conditions for convergence [33–36] are indicated. The rationality and effectiveness of the proposed method in multiple attenuation are verified by synthetic and real data. The method in this paper is compared with least squares parabolic Radon transform (LSPRT) and sparse parabolic Radon transform based on *L*<sup>1</sup> regularization (SPRT*L*1) for multiple attenuation. This method improves the precision and focusing ability of Radon transform, obtains a better result of multiple suppression, reduces the loss of primary reflections, and improves the interpretability of seismic data.
