An efficient plug-and-play regularization method for full waveform inversion

Abstract

Nonlinear inverse problems arise in various fields ranging from scientific computation to engineering technology. Inverse problems are intrinsically ill-posed, and effective regularization techniques are necessary. The core of a suitable regularization method is to introduce the prior information of the model via an explicit or implicit regularization function. Plug-and-play regularization is a flexible framework that integrates the most effective denoising priors into an iterative algorithm, and it has recently shown great potential in the solution of linear ill-posed problems. Unlike traditional regularization methods, plug-and-play regularization does not require an explicit regularization function to represent the prior information of the model. In this work, by using total variation, block-matching and three-dimensional filtering, and fast and flexible denoising convolutional neural network denoisers, we propose a novel iterative regularization algorithm based on the alternating direction method of multipliers method. The combination of total variation and block-matching three-dimensional filtering regularizers can take advantage of the sparsity and nonlocal similarity in the solution of inverse problems. When combined with traditional and novel regularizers, deep neural networks have been shown to be an effective regularization approach, which can achieve state-of-the-art performance. Finally, we apply the proposed algorithm to the full waveform inversion problem to show the effectiveness of our method. Numerical results demonstrate that the proposed algorithm outperforms existing inversion methods in terms of quantitative measures and visual perceptual quality.

1. Introduction

A nonlinear inverse problem is generally expressed as the following least-squares optimization problem

where the model |$\bf {m}$| needs to be inverted from the measured data |$\bf {d}$|⁠. The forward modeling operator |$\mathcal {F(\bf {m})}$| is parametrized by |$\bf {m}$|⁠, which is governed by appropriate partial differential equations. Mathematically speaking, the problem (1) is usually nonconvex and severely ill-posed. Regularization methods are necessary to overcome the ill-posedness of the problem (1) and steer the overall procedure for optimization into an ideal direction (Tikhonov & Arsenin 1977).

After introducing regularization, problem (1) becomes

where the regularizer or regularization function |${\mathcal {R}(\bf {m})}$| includes prior knowledge about the expected solution |$\bf {m}$|⁠, and |$\lambda \gt 0$| is the regularization parameter (Wang 2016).

A popular choice for the regularization term |$\mathcal {R}(\bf {m})$| is the total variation (TV) function |$\Vert \nabla \bf {m}\Vert _{1}$|⁠, i.e., |$\mathcal {R}(\bf {m})$| is the |$l^{1}$|-norm of the model gradient (Rudin et al.1992). TV regularization assumes that the desired model is sparse in the derivative domain, and thus promotes piece-wise smooth solutions while preserving sharp discontinuances (Li & Han 2016). It is also very effective in full waveform inversion (FWI) problems, especially for large-velocity contrast structures, such as salt bodies embedded in sedimentary layers (Esser et al.2016; Lin & Huang 2014; Peng et al.2018). The TV norm is nondifferentiable and some proximal gradient-type algorithms should be used to solve the nonsmooth optimization problem (2) with |$\mathcal {R}(\bf {m})=\Vert \nabla \bf {m}\Vert _{1}$|⁠, where the proximal operator of |$\mathcal {R}$| is defined as

Note that the above proximal operator can be viewed as a Gaussian denoiser, which means it is possible to solve (1) with powerful regularizers embedded in excellent denoising methods such as block-matching three-dimensional (BM3D; Burger et al.2012; Dabov et al.2007) or deep learning denoisers (Sun et al.2019; Wu & Lin 2020; Zhang & Alkhalifah 2019). The core of the BM3D is to simultaneously exploit both the local sparsity of wavelet coefficients and the inherent nonlocal self-similarity properties of grouped image blocks. These two priors are very useful for image denoising, and the BM3D is often more effective compared with other existing denoisers.

The newest trend is deep learning, which has arisen as a promising mechanism and provides new ideas for solving ill-posed inverse problems. Deep-learning-based image denoisers have been developed rapidly, and it may be even possible to provide more performance gain than classical state-of-the-art denoisers such as TV or BM3D. For instance, Zhang et al. (2018) proposed a fast and flexible denoising convolutional neural network (FFDNet) by introducing a tunable noise level map as input. As a result, FFDNet exhibits excellent performance on both synthetic noisy images and real-world noisy images, and has flexibility in handling various types of noises.

With the mathematical equivalence between the proximal operator and denoiser, Venkatakrishnan et al. (2013) proposed plug-and-play (PnP) priors, and have received great success for regularized linear ill-posed problems (Cascarano et al.2022; Chan et al.2017; Rasti-Meymandi et al.2022). Plug-and-play is a highly flexible and effective framework for regularizing ill-posed problems by using more advanced denoisers within the alternating direction method of multipliers (ADMM; Boyd et al.2010) or other proximal-type algorithms. In this paper, we exploit the fact that different pieces of prior information could promote each other, and we further explore the effectiveness of PnP for solving the nonlinear ill-posed problem. More specifically, three popular denoisers (i.e., TV, BM3D, and FFDNet) are integrated into the ADMM algorithm by means of the PnP prior framework. Apparently, the proposed method takes full advantage of the sparsity (both global and local) and nonlocal self-similarity priors of the desired model, as well as the deep priors characterized by deep denoiser networks. The global sparsity is enforced by the TV norm for constructing the entire consistency and preserving edges, while the local sparsity is captured by the BM3D denoiser to preserve more features by exploiting model self-similarity prior. Such advantages suggest that there is the possibility of applying a powerful CNN denoiser to improve the ill-posedness of the inverse problem.

The paper is organized as follows. Section 2 depicts the proposed PnP regularization, which combines TV, BM3D, and FFDNet denoisers. Section 3 focuses on the applications to the FWI problem, and the experimental results and parameter selection are given. Section 4 presents the conclusion.

2. The proposed PnP regularization method

With the sparsity, nonlocal self-similarity, and deep priors all incorporated, problem (2) is formulated as a composite regularization problem

where |$\alpha$|⁠, |$\beta$|⁠, and |$\gamma$| are three regularization parameters, |$\Vert \cdot \Vert _{TV}=\Vert \nabla \bf {m}\Vert _{1}$| is the TV norm, and the last two items |$\Phi (\bf {m})$| and |$\mathcal {D}(\bf {m})$| are implicit regularizers that are deduced from the BM3D and FFDNet denoisers, respectively.

We set |$\mathcal {R}(\bf {m})=\alpha \Vert \bf {m}\Vert _{TV}+ \beta \Phi (\bf {m})+\gamma \mathcal {D}(\bf {m})$| in problem (3), by introducing an auxiliary variable |$\bf {v}$|⁠. Problem (3) can be converted into the following constrained optimization problem

Consider its augmented Lagrangian function

where |$\bf {u}$| is the Lagrange multiplier, and |$\rho$| is a penalty parameter.

To solve problem (4), we adopt the ADMM framework, which decomposes this optimization problem into simpler denoising subproblems, and then we solve the subproblems one by one. Algorithm 1 presents the pseudo code of our method, called plug-and-play regularization, or the PPR method. In step (a), the data fidelity term can be computed with any gradient-based optimization method including the steepest-descent method (Choi et al.2005) or quasi-Newton method, to obtain an initial value as a noisy solution. In steps (b)–(d), the variable |$\bf {v}$| can be split into three variables, |$\bf {v}_{1}$|⁠, |$\bf {v}_{2}$|⁠, and |$\bf {v}_{3}$|⁠, and the corresponding denoising operation is performed over each of them. In step (e), the Lagrange multiplier |$\bf {u}$| is updated. The convergence of the PPR method is guaranteed by the ADMM framework.

In fact, in the PPR method, each regularization subproblem corresponds to a denoising problem (Brifman et al.2016). Generally speaking, the proximal map

is equivalent to a denoising operation, where |$\bf {m}$| is a noisy intermediate solution, and |$\mathcal {R}$| denotes a regularizer corresponding to a denoiser with a denoising threshold |$\sigma$|⁠. In this paper, the denoising thresholds of denoising steps (b)–(d) are represented by |$\sigma _{i}(i=1,2,3)$| respectively, which are decided by the regularization parameter and penalty parameter |$\rho$| through the following equations:

In particular, in the implementation of the above denoising operations, it can be observed that the current denoising results are integrated into the noisy image of the next denoising step in Algorithm 1 to improve the inversion efficiency. Therefore, inspired by the spirit of the proposed PPR method, the solutions of optimization problems (b)–(d) of denoising stage can be denoted by

3. Applications in full waveform inversion

3.1. Acoustic full waveform inversion

Mathematically, FWI is a typical nonlinear ill-posed problem, which seeks to obtain a physical model that agrees with observed seismic data as well as possible.

Succinctly, frequency-domain FWI can be written as the following nonlinear least-squares optimization problem

where |$\bf {m}$| denotes the velocity model to be inverted, |$\omega$| is the angular frequency, |$A_{\omega }(\bf {m})=({\omega ^{2}}/{\bf {m}^{2}}) +\nabla ^{2}$| is the Helmholtz operator with a perfectly matched layer (PML) boundary condition, |$\bf {d}$| is the observed seismic data, |$\bf {Q}$| is the source matrix, and |$\bf {P}$| denotes a measurement operator that projects the solution of the Helmholtz equation on to the receiver locations. In this paper, we ignore the density parameter and only invert for acoustic velocity. In the discrete case, the unknown velocity model |$\bf {m}$| is a matrix of dimension |$N_{z}\times N_{x}$|⁠, where |$N_{z}$| and |$N_{x}$| are the numbers of grid points along the z and x directions, respectively. The limited memory quasi-Newton L-BFGS method (Zhu et al.1997) is utilized to solve step (a) in Algorithm 1 when the gradient of |$\mathcal {F}(\bf {m})={\bf {P}}A_{\omega }(\bf {m})^{-1}\bf {Q}$| is computed via the adjoint-state method (Plessix 2006).

The developed PPR algorithm for solving the acoustic FWI problem is illustrated in the flowchart in figure 1. Algorithm 2 (called the PPR–FWI algorithm) summarizes the PPR regularization procedure that allows a computer to solve FWI problems (5).

3.2. Numerical experiments

The quality of the inverted velocity models is measured using the peak signal-to-noise ratio (PSNR), the structural similarity index measure (SSIM), and the root mean square error (RMSE). We compared the proposed PPR–FWI method with three other algorithms: (1) FWI with the TV denoiser (TV–FWI); (2) FWI with the BM3D denoiser (BM3D–FWI); (3) FWI with FFDNet networks (FFDNet–FWI). In the experiments, we set the outer loop as |$L=4$| and employ 10 iterations in the L-BFGS iteration stage, which is consistent among four FWI methods to ensure a fair and valid comparison. From a practical perspective, noise is inevitable in the collected seismic data. In all experiments, we also consider a noise version of the seismic data |$\bf {d}$|⁠, to which 5% white Gaussian noise is added. In the frequency domain FWI, we parallelize over frequencies and invert all frequency data components simultaneously. In the simulation, the FFDNet model and parameters from https://github.com/cszn/FFDNet are directly used.

In the PPR–FWI method, to satisfy the condition required in the convergence theorem stated in Cascarano et al. (2022), we choose |$\rho _{l+1}=(l+1)(1+\epsilon )^{l+1}$| with |$\epsilon \gt 0$|⁠; pertinent details can be found in the appendix section of Cascarano et al. (2022). For the parameter |$\epsilon$|⁠, we find that |$\epsilon =10^{-3}$| will lead to a good numerical performance. The three regularization parameters for the PPR–FWI algorithm are set as |$\alpha =1\times 10^{-3}$|⁠, |$\beta =0.02$|⁠, and |$\gamma =10^{-4}$|⁠, respectively. In addition, the regularization parameters of comparative inversion methods for different numerical examples can vary, which will be described in the following subsections.

3.2.1. Marmousi model

We perform the modified Marmousi synthetic benchmark model as the first numerical example. The dimensions of the model are |$3060 \times 9600$| m|$^2$|⁠, with grid spacings of 4 m, as shown in figure 2(a). The frequency bandwidth is [2.9,24] Hz. The initial velocity model (figure 2(b)) is estimated by smoothing the original model. We place |$N_{s}=75$| sources and |$N_{r}=225$| receivers located uniformly on the model surface. For the implementation of the TV–FWI, BM3D–FWI, and FFDNet–FWI methods, the regularization parameters are set to 0.8, 0.2, and 0.15, respectively.

The inversion results with four different regularization schemes are shown in figures 3 and 4. Specifically, figure 3 shows inversion results without Gaussian noise in the data, and results with noise are given in figure 4. As shown in figures 3 and 4, the inversion results of TV–FWI, BM3D–FWI, and FFDNet–FWI produce a reasonable velocity model, while the features in the deeper part of the model are poorly restored. By contrast, FWI with PPR regularization improves the reconstruction quality for seismic imaging to some extent. The inverted model obtained by the PPR–FWI method (depicted in figures 3(d) and 4(d)) contains more structural details and fewer inversion artifacts than other methods, which further demonstrates that our PPR–FWI algorithm obtains the best performance. A comparison of the velocity profiles shown in figures 5 and 6 further demonstrates the superiority of our PPR–FWI method.

To further illustrate the performance of the proposed algorithm, the quantitative results of the Marmousi model are reported in Table 1. We can observe that the proposed algorithm realizes a relatively highly competitive performance. It basically outperforms the other three algorithms according to PSNR, SSIM, and RMSE values.

3.2.2. Modified BP 2004 model

As a second example, a numerical experiment based on a modified BP 2004 Model (see figure 7(a)) is presented to test the performance of the PPR method. The model is discretized to |$N_{z}\times N_{x}=128\times 512$| with a grid spacing of 8 m. The frequency bandwidth is [2,26] Hz. The initial model (see figure 7(b)) is obtained by applying Gaussian smoothing to the true model. We evenly distribute 43 sources and 129 receivers in the horizontal direction. The manually selected regularization parameters for the TV–FWI, BM3D–FWI, and FFDNet–FWI methods are set as |$10^{-4}$|⁠, 0.5, and |$6\times 10^{-3}$|⁠, respectively.

Figures 8 and 9 show the inversion velocity model with four different regularization methods under noiseless and noise conditions. If we compare the inversion results carefully, in the areas shown by the rectangular boxes, the PPR–FWI method (depicted in figures 8(d) and 9(d)) can better reconstruct the detailed features and reveal the “tooth” part of the velocity model. Figures 10 and 11 show the velocity profiles of the inversion results at |$z=960$| m under noise-free and noisy cases, respectively. Compared with the TV–FWI, BM3D–FWI, and FFDNet–FWI methods, the inversion result obtained by the proposed PPR–FWI algorithm is closer to the true velocity of the model. Furthermore, Table 2 shows the quantitative inversion results. As can be seen from these metrics, the proposed PPR–FWI method achieves a satisfactory improvement, which further indicates that prior information utilized in the PPR–FWI algorithm can obtain the highest reconstruction quality.

3.2.3. SEG/EAGE salt model

We finally experiment by using our method on a two-dimensional SEG/EAGE salt model. The true model is depicted in figure 12(a), which is discretized to |$N_{z}\times N_{x}=128\times 256$| with a spacing of 5 m. The frequency bandwidth is [2,26] Hz. The synthetic data are generated by 31 sources and 93 receivers near the surface. Figure 12(b) shows the initial velocity model. Note that, to mitigate the reliance on the initial value, there is a total lack of information about lateral variation in the initial velocity model. For the implementation of TV–FWI, BM3D–FWI, and FFDNet–FWI methods, the corresponding regularization parameters are set as |$2\times 10^{-4}$|⁠, 0.5, and |$5\times 10^{-5}$|⁠, respectively.

Similarly, this experiment is also studied in the noise-free and noisy cases (see figures 13 and 14). We make the following conclusions in terms of visual effects and performance metrics. In the noiseless situation, as shown in figures 13(a) and (c), the inversion of TV–FWI and FFDNet–FWI leads to lower visual clarity while strong artifacts occur in the background part. By examining the inversion result obtained by the BM3D–FWI method (figure 13(b)), the salt structure is better maintained while the background part of the inverted model is poor. The inversion results for the PPR–FWI method are shown in figure 13(d). Notably, we find that our proposed method results in a more accurate reconstructed image, and the structural details and boundaries of the model are well preserved. In the noisy case, the inversion results obtained by the proposed PPR–FWI algorithm (see figure 14(d)) generate the slightest deviation with a more precise background and more obvious observation of the inverted model.

Figures 15 and 16 show the velocity profiles of the inversion results obtained by four inversion algorithms at |$x=900$| m. Although the inversion velocities of the four algorithms are all close to the true velocity, the PPR–FWI algorithm obtains better fitting effects between the inversion velocity and the true velocity. Table 3 shows the quantitative evaluation of the reconstructed results. As can be seen from Table 3, our PPR–FWI method achieves the best performance for the quantitative metrics measured by PSNR, SSIM, and RMSE.

4. Conclusions

We design a flexible PnP regularization algorithm by incorporating the sparsity and nonlocal self-similarity priors, as well as the deep priors learned by deep denoiser networks. What is more, three powerful denoisers (TV, BM3D, and FFDNet) are plugged to solve regularization subproblems. The application of the PnP regularization algorithm in the frequency-domain FWI problem leads to promising performance. In the experiments based on the Marmousi model, SEG/EAGE salt model, and modified BP 2004 model, our algorithm shows a superior ability in reducing artifacts and preserving more details to a great extent. Furthermore, because performing L-BFGS iterations in the FWI process is the most time-consuming stage, the computation time for all inversion methods is almost the same. However, the proposed method still has room for improvement in several aspects. The convergence analysis has drawn more attention and it is challenging to verify whether general denoisers meet the required conditions. In addition, the adaptive selection of regularization parameters and the improvement of computational efficiency still need further exploration in practical applications.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (grant no. 42274166)

Data Availability

Data will be made available on request.

Conflict of interest statement. None.

References