Improving full-waveform inversion based on sparse regularization for geophysical data (2024)

Abstract

1. Introduction

Full-waveform inversion (FWI) is a powerful technique for constructing high-resolution models of subsurface seismic velocity structures, which is extensively employed in geophysical exploration, seismic monitoring, and oil and gas reservoir characterization (Virieux and Operto 2009). In the frequency-domain FWI, the seismic modelling technique simulates the forward and backward wavefield solution by solving the wave equation (Wang and Rao 2009) . Subsequently, simulated and observed seismic records are matched iteratively to solve the inverse problem and obtain a velocity model (Warner and Guasch 2016).

FWI encountered several challenges, including the tendency to settle into local minima, constraints imposed by azimuth and acquisition systems, and the emergence of imaging artefacts. The inversion outcomes are notably compromised by numerical inaccuracies in forward data and the inefficiencies of low-density acquisition systems (Zhang 2010). The frequency-domain forward modelling generates forward wavefields, describing the aspects of seismic wavefield propagation relationship between the source and receiver at a specific frequency (Wu and McMechan 2021). The three-dimensional wavefield stack mainly contains multi-scale information, ranging from low to high frequencies, which can represent the subsurface's structures and properties at varying scales (Ovcharenko et al. 2018). However, in real cases, the forward and backward wavefield is often influenced by numerical error, leading to inaccurate forward wavefields and further resulting in degraded quality of imaging and inversion outcomes (Alkhalifah et al. 2021). Although many scholars have tried to develop new schemes for accurate seismic modelling in the frequency domain (Jo et al. 1996, Chen 2012, Takekawa and Mikada 2018), incrementing the computational costs is generally inevitable, especially in field studies. It makes it challenging to apply FWI in practice. More importantly, in situations where high-density acquisition is challenging to attain, the effect of spatial aliasing must be considered. The spatial aliasing effect arises from the unsatisfied Nyquist–Shannon sampling theorem during the sampling process, which leads to the possibility of high-frequency components being misinterpreted as low-frequency components. It causes low-frequency information to be lost, consequently affecting the accurate recovery of the subsurface velocity structure (Shaiban et al. 2022). The most intuitive effect of this phenomenon is that there will be high-frequency artefacts in the high-frequency band. The traditional approach of increasing the sampling rate will augment the amount of data but may also elevate the computational cost (Liu and Fomel 2011).

Furthermore, FWI encounters numerous issues in cases of low-density sampling (which implies that the seismic data are sparsely sampled). First, it is well known that low-frequency information is essential for FWI, as it helps to overcome the local minima problem. The abundance of low-frequency information facilitates the convergence of the inversion to the global optimal solution rather than a locally optimal solution (Kan et al. 2023). Additionally, low-density sampling may contribute to the limited offset distance issue, that is, insufficient receiver coverage and a decline in the inversion resolution of the subsurface structure (Operto et al. 2009). These limitations can render the FWI update process more reliant on the initial model, diminish the accuracy of the inversion results, and affect the identification and interpretation of subsurface structures by FWI (Operto and Miniussi 2018).

On the basis of these problems, we propose to use the sparse relaxation regularized regression (SR3) algorithm in FWI to optimize the forward and backward wavefield obtained from the frequency-domain numerical modelling to improve the accuracy of FWI and achieve high-resolution inversion (Zheng et al. 2019). As a more abundant and comprehensive sparse regression algorithm, the SR3 algorithm is designed to optimize the forward wavefield in FWI, including denoising and interpolation functions. Like the quadratic penalty method, this improved algorithm provides multiple specification options for different conditions and models with broader applicability and superior convergence speed. The optimized forward wavefield shows better continuity, thus providing a solid foundation for improved performance of FWI. The SR3 algorithm, as an optimization algorithm, enhances the fundamental regression problem into a new regression problem with additional constraints through the assistance of auxiliary matrices. The main difference between the SR3 algorithm and conventional regression problems is that it tightens the function space of the error function, which can compress the nonlinear wavefield loss space, thereby shrinking the range of singular values. As a result, the iteration speed of regression based on the SR3 algorithm is faster than that of conventional regression algorithms, and it also encourages sparse solutions (Champion et al. 2020). Specifically, a more appropriate image processing technique for the spatial aliasing effect that may occur in high-frequency data is to use matrix reconstruction and interpolation techniques (Liu and Fomel 2011). The SR3, as a regularization algorithm, has a specific interpolation effect by introducing regularization terms and sparse constraints. The SR3 can reduce the effect of aliasing, while sparse constraints help to reduce the complexity of the solution space, thus improving the accuracy of matrix reconstruction and interpolation (Brunton et al. 2016). In addition, processing low-frequency data completion and interpolating down-sampled data are also crucial for FWI. The SR3 algorithm adopts multi-regularization algorithms through regression and iterative processes to regularize the sparse basis of the simulated data and, through fast iterations, the data are recovered and completed to solve the optimization problem (Erichson et al. 2020). Moreover, this algorithm can find a balance between the sparse constraints and the smoothness of the solution by choosing an appropriate regularization parameter, thereby providing a more stable denoising and data completion effect (Champion et al. 2020).

In this paper, based on the frequency-domain FWI, the SR3 algorithm is innovatively introduced to achieve data complementation and reconstruction of FWI for low-density acquisition data, anti-spatial blending effect for high-frequency data, and improved resolution for high background noise datasets. Two hom*ogeneous medium models and two synthesis models (2004 BP Model and SEG/EAGE overthrust model) are used to verify the superiority of our approach by comparing Tikhonov regularized FWI with SR3-based FWI.

2. Sparse relaxation regularized regression

In FWI technology, the frequency-domain seismic modelling technique has garnered considerable attention and research due to its outstanding precision, efficient computational speed, and independence from overly idealized time sampling rates (Jiang et al. 2021). The frequency-domain seismic modelling technique does not rely on an excessively idealized time sampling rate, thus offering greater flexibility in practical applications (Ravaut et al. 2004). A crucial part of this algorithm typically involves solving the Helmholtz equation (Appelö et al. 2020). This foundational equation |${{\bf U}}({{\bf X}},{{\bf Z}};{{\bf W}})$|provides an essential mathematical framework for our understanding of geophysical phenomena:

where |${{\bf W}}$| is the angular frequency, |${{\bf M}}$| is the squared slowness, |${\mathop{\rm{F}}\nolimits} ({{\bf W}})$| is the source signature, |$\Delta $| is the Laplacian operator, |${{\bf U}}$| denotes the wavefield, |${\mathop{\rm{G}}\nolimits} ({{\bf X}} - {{{{\bf X}}}_{\rm{S}}})$|⁠, and |${\mathop{\rm{g}}\nolimits} ({{\bf Z}} - {{{{\bf Z}}}_{\rm{S}}})$| denote the source at points |${{{{\bf X}}}_{\rm{S}}}$| and |${{{{\bf Z}}}_{\rm{S}}}$|⁠.

In the subsequent steps, the conventional algorithm ingeniously maps the model space to the data space, enabling a better understanding and handling of the model. Within the data space, FWI performs a least squares calculation on the discrepancy between observed data and calculated data to assess their agreement (Alkhalifah et al. 2021). The resultant least squares value reflects the model increment in the model space, representing the extent of improvement from the model's current state to a more optimal state. FWI then uses this model increment to compute the update gradient, which expresses the direction of model change and reflects the rate of model change, providing a guide for optimizing the model (Choi and Alkhalifah 2012). Equipped with this gradient, we are well-prepared for the next iteration: exploring a more optimal model based on the current model and gradient, enhancing the agreement between observed data and calculated data, and thereby attaining a more precise and reliable model:

where |${\mathop{\rm{F}}\nolimits} ({{\bf M}})$| is the calculated data and |${{{{\bf D}}}_{{\rm{obs}}}}$| is the observed data.

However, in the minimization process, the effectiveness and accuracy of the inversion results are directly affected if the quality of the calculated data matrix we are using is low or whether data loss has occurred due to low-density sampling, even when more advanced inversion algorithms are applied (Operto and Miniussi 2018). The integrity and quality of data play a crucial role in this process, with data loss and noise causing the inversion results to deviate from actual values (Liu and Fomel 2011). Therefore, we employ an optimization algorithm to complete data and denoise the calculated data matrix before the minimization operation, thereby enhancing the quality and completeness of the data and consequently improving the accuracy of inversion results.

In industrial and scientific domains, regularization regression problems are widespread, with sparse regression occupying a crucial position in applications like compressive sensing and data completion. Sparse regression is particularly advantageous in selecting the most significant features from extensive datasets, often containing numerous unimportant or redundant features, which is precious in handling high-dimensional data scenarios. In data completion, sparse regression can help us predict possible values of missing data based on existing partial data, thereby completing the dataset; thus, sparse regression plays a highly significant role in these domains (Erichson et al. 2020). We incorporate this approach into the optimization of the forward modelling results in FWI:

where |${{\bf A}} \in {{\mathbb{R}}^{m \times d}}$| is the data-generating model for the actual data |${{\bf D}}$|⁠, |${{\bf U}} \in {{\mathbb{R}}^d}$| is the recovery data, |$\lambda $| is the regularization parameter, |${\mathop{\rm{R}}\nolimits} ({g} )$| is any regularization form, and |${{\bf S}} \in {{\mathbb{R}}^{n \times d}}$| is a linear map.

In modern optimization, non-smooth regularizations are widely used in various applications. We innovatively apply it to seismic surveys to replace intensive acquisition with optimization while denoising the signal to achieve seismic data processing and waveform inversion at a lower cost (Gholami et al. 2021). The SR3 algorithm demonstrates more pronounced advantages than other regression algorithms (Zheng et al. 2019). One of its primary benefits is its capacity to effectively identify sparse signals, a significant effect brought about by introducing an auxiliary matrix. Applying such an auxiliary matrix becomes particularly crucial when the data contains noise, and the model parameter matrix is ill-conditioned. The presence of this auxiliary matrix aims to aid us in more efficiently recovering the data matrix, reducing the artefacts induced by noise and data loss, thereby improving and optimizing the final inversion results (Aravkin et al. 2014).

Additionally, our approach offers flexibility, allowing us to adjust the regularization function to accommodate various scenarios. We can employ plural forms of regularization functions depending on the sparsity characteristics of different models. The flexibility in selecting regularization functions reflects our approach's strong adaptability and flexibility, enabling us to achieve optimal results under many different settings and conditions of the constraint function and the fidelity of the relaxation problem (Esser et al. 2018). First, we need to construct an auxiliary variable |${{\bf W}}$| to relax Equation(3):

where |${{{{\bf W}}}^k} = {\rm{pro}}{{{\rm{x}}}_{\lambda k{\mathop{\rm{R}}\nolimits} }}({{{{\bf U}}}^k})$| is the auxiliary variable, which gradually approaches |${{\bf U}}$|⁠; |${\rm{pro}}{{{\rm{x}}}_{\lambda k{\mathop{\rm{R}}\nolimits} }}$| is the proximity operator (prox) for R. λ is the penalty parameter, and k is the relaxation parameter, where k controls the degree of relaxation. The scalar |${{\| \centerdot \|}_{\mathop{\rm{P}}\nolimits} }$| represents a different form of the regularization used in the optimization, which can be flexibly used as different regularization functions depending on the sparsity of the data matrix, such as L1 and L2 norms and even a nuclear norm.

Associated withEquation(4) is a form of the value function, which allows us to depict this relaxed framework precisely (Champion et al. 2020). A value function offers a quantitative representation of a problem, capable of revealing the nature and characteristics of the problem, aiding in our understanding of its essence, and providing direction for problem-solving. Specifically, this value function is obtained by minimizingEquation(4) in recovery data. This minimum corresponds to the value function we are discussing, which to some extent reflects the properties of the optimal solution to the problem:

We assume that |${{{{\bf H}}}_k} = {{{{\bf A}}}^{\rm{T}}}{{\bf A}} + k{{{{\bf S}}}^{\rm{T}}}{{\bf S}}$| is invertible, so that |${{\bf U}}({{\bf W}}) = {{\bf H}}_k^{ - 1}({{{{\bf A}}}^{\rm{T}}}{{\bf D}} + k{{{{\bf S}}}^{\rm{T}}}{{\bf W}})$| is unique, and

which provides a closed-form expression for

Equation(4) then reduces to

To elucidate the superiority of the SR3 algorithm more accurately, we demonstrate that partial minimization improves the condition of the problem inFig.1. InFig.1a, the coloured ellipses depict the contours of |$\| {{{\bf AU}} - {{\bf D}}} \|_{\mathop{\rm{F}}\nolimits} ^2$|⁠, while inFig.1b, the contours of |$\| {{{{{\bf F}}}_k}{{\bf W}} - {{{{\bf g}}}_k}} \|_{\mathop{\rm{F}}\nolimits} ^2$| are vividly portrayed as a circle (Zheng et al. 2019). InFig.1a, we exhibit the contour lines of the quadratic part similar to the LASSO problem (coloured ellipses) and the approximate solution paths (red solid line) in horizontal projection. InFig.1b, the contour lines of the quadratic part of the SR3 loss function (coloured circle) and the corresponding approximate paths (red solid line) of the SR3 solution in the relaxed coordinates |${{\bf W}}$| are shown. Additionally, the grey diamonds indicate the contour lines of the L1-norm of the LASSO problem in each coordinate group.

Figure1 shows that for the widely applied class of LASSO-like problems, the properties of |${{\bf F}}$| are generally superior to |${{\bf A}}$|⁠, particularly in terms of the condition number, where |${{\bf F}}$| is typically smaller than |${{\bf A}}$|⁠. The ratio of the maximum to minimum singular values of |${{\bf F}}$| is smaller, thereby compressing the contour lines into a shape closer to a sphere, which accelerates convergence and enhances performance. Moreover, executing proximal gradient descent solely in |${{\bf W}}$| naturally resolves these types of problems. The formulas for |${{\bf F}}$| and |${{\bf G}}$| can also be applied to acceleration methods, such as the algorithm fast iterative shrinkage-thresholding (FISTA). Overall, the SR3 algorithm reduces the singular values of F relative to |${{\bf A}}$| and has a weaker impact on small singular values. This effect ‘squeezes’ the contour lines into a near-spherical shape, accelerating convergence and enhancing performance.

Finally, the algorithm requires two parameters to be specified simultaneously; the parameter λ determines the strength of the regularization, while k determines the degree of relaxation. We set the coefficient threshold to |$\mathcal{J}$|⁠, then,

WithEquation(11), we can change the two-parameter selection problem to a single-parameter selection and significantly improve the algorithm. In addition, we suggest a cross-validation approach to achieve parameter tuning of the automatic strategy.Equation(11) is an empirical formula regarding parameter selection that we derived based on our experience in practical computations. During cross-validation, we used an intermediate parameter to represent the ratio of two unknown parameters, λ and k. We first set a range for one of the parameters and then input it into the algorithm. The algorithm iteratively calculates and optimizes the error between the computed and observed wavefield. We then select the value of |$\mathcal{J}$| that corresponds to the minimum error, which gives us the optimal ratio of λ to k. After fixing one of these parameters, then the value of the other parameter can be determined by these steps.

3. Numerical examples

We used two types of hom*ogeneous medium and two distinct artificially synthesized velocity models to evaluate the performance of FWI optimized on the basis of the SR3 algorithm. The specific models include the single-layer hom*ogeneous medium and double-layer hom*ogeneous medium mentioned in this section and the 2004 BP and SEG/EAGE Overthrust models. We employed our experiments' frequency-domain forward modelling algorithm and introduced the perfectly matched layer as the absorption boundary condition (Pratt 1999).

Furthermore, FWI is notably susceptible to the interference of noise. To articulate the impact of noise on inversion outcomes quantitatively and to juxtapose the noise-reducing capabilities of distinct algorithms, we define the intensity of noise present in a signal using the formula for signal-to-noise ratio (SNR) (de Ridder and Dellinger 2011):

where the |${{\bf E}}$| is the noise. To further illustrate the outlier suppression capabilities of the SR3-based FWI, we employ a signal with a low SNR of 10 dB across different models.

To evaluate the simulation accuracy of different algorithms for the model, we employ the model error:

where the |${{{{\bf M}}}_{{\rm{true}}}}$| and |${{{{\bf M}}}_{{\rm{inv}}}}$| represent the actual and inversion model velocity, respectively (Warner and Guasch 2016).

When implementing regression algorithms, we have several options to optimize our model. Specifically, we can consider using different regularizations or a combination of various regularization methods, which can help us achieve more effective optimization inEquation(5).Figure1 demonstrates a three-dimensional visualization of two common regularizations, L1 regularization inFig.2a and L2 regularization inFig.2b. Both regularizations are widely applied in geophysics and engineering, and each possesses unique advantages.

For the L1 regularization, we can define this as the |${\rm{P}}$| regularization inEquation(5). The L1 regularization emphasizes the sparsity of regularization, a feature that can significantly enhance the denoising capability of the algorithm, helping the model to focus on the most essential features and, thus, improve the interpretability of the model. On the other hand, we could also select the L2-norm as the |${\rm{P}}$| regularization. The L2 regularization emphasizes the generalization capability of the algorithm. It can improve the smoothness and stability of the model and effectively prevent overfitting.

To demonstrate the core theory of the SR3 algorithm we adopted, inFig.1 we present an illustrative diagram of the gradient iteration process using the LASSO problem as an example.

Figure1a represents a schematic diagram of the update path of the proxy gradient under the action of the L1 regularization for an elliptical contour loss function.Figure1b highlights the core idea of the SR3 algorithm, which ‘tightens’ the elliptical contour of the loss function to an approximate circle, accelerating the convergence speed and performance of regression computations. This accelerated convergence enables us to embed this algorithm into conventional FWI without the pressure of computational memory and time consumption.

3.1. Single-layer hom*ogeneous medium

First, the effectiveness of the SR3 algorithm in the application of moderately high-frequency spatial mixing phenomena has been thoroughly tested in the context of a single-layer hom*ogeneous medium.Figure3a shows the basic structure of a single-layer hom*ogeneous medium.Figure3b provides a one-dimensional velocity model of a single-layer hom*ogeneous medium that achieves a velocity of 2 km s⁻¹. Based on this,Fig.3 parts c–e further show the wavefield for the analytical solution of the monolayer hom*ogeneous medium at frequencies of 10, 13, and 15 Hz.

Figure4 parts a–d show the results of the conventional algorithm for the monolayer hom*ogeneous medium at 10 Hz. The discretized colour map is intended to improve recognition performance.Figure4 parts e–h show the results after processing by the SR3 algorithm for the monolayer hom*ogeneous medium at 10 Hz.

To better compare the forward results at different frequencies in the case of the multi-scale algorithm and the processing effects of the SR3 algorithm, we have also compared the wavefield results and processing results at 13 and 15 Hz, respectively. As shown inFig.5a–h, the wavefield of the conventional algorithm at 13 Hz is compared with the wavefield processed by the SR3 algorithm using the same methods and order of comparisons as inFig.4a–h.

In addition,Fig.6a–h compare the conventional wavefield results and the SR3 processed wavefield for the 15 Hz case, with the same methods and order of comparisons as inFig.4a–h. The discretized colour map is intended to improve recognition performance.

3.2. 2004 BP model

Further, we test our proposed algorithm in synthetic data to verify its performance in near-real situations.Figure7a shows the actual velocity model for the 2004 BP model, andFig.7b shows the initial velocity model used for the FWI with a linear velocity increase ranging from 1.5 to 5 km s⁻¹.

The 2004 BP model, as a highly representative salt body model, is extensively used in FWI. The choice of the salty model is motivated by the presence of a high-velocity salt anomaly in its central part, which creates a significant velocity contrast with its surroundings. This contrast leads to considerable difficulties in accurately delineating the salt body contours and surrounding velocities under conditions of extremely low SNR. High noise levels exacerbate these challenges. Consequently, the salty model is an ideal test bed in this experiment for assessing the precision and effectiveness of both conventional and proposed algorithms in high-accuracy inversion under complex conditions. To test the effect of noise, we add random background noise to the recorded data set, as shown inFig.8. The middle part of the 2004 BP model describes the geological characteristics of the Eastern/Central Gulf of Mexico and offshore Angola. Because the middle part of the model is composed of a high-velocity salt body, the division of the salt body is one of the difficulties in the inversion of the model. In addition, the channels are located near the salt body, which will further affect the velocity inversion and increase the difficulty of inversion (Billette and Brandsberg-Dahl 2005).

The source–receiver wavefield of the 2004 BP model is demonstrated for the 3 Hz scenario, as illustrated inFig.8a. Using this as a basis, we integrated 10 dB of random noise to mimic the interference that would typically arise in actual production. The noise equation is provided inEquation(12).Figure8b depicts the 10 dB random noise, andFig.8c represents the wavefield after the clean wavefield has been linearly combined with the random noise. Consequently,Fig.8d presents the subsampled wavefield, which we simulate as a low-quality wavefield influenced by downsampling. This low-quality wavefield is evenly distributed to incorporate the missing trace, mixed with random noise, and overlaid linearly onto the clean wavefield, resulting in the wavefield matrix containing 10 dB noise and the missing trace, as shown inFig.8e. The subsequent FWI procedures hinge on the inversion iterations of the low-quality forward wavefield established through this process.

We compute the source–receiver wavefield of the 2004 BP model within the frequency domain, as illustrated inFig.9.Figure9a presents the clean wavefield, whileFig.9b offers a side view, adjusting the perspective from (a) to a view from (105, 1). In view (a, b), the ‘a’ represents the horizontal rotation angle used to rotate the view around the z-axis, while ‘b’ is the vertical rotation angle used for tilting the view in a direction perpendicular to the horizontal-vertical plane. We follow the previously mentioned simulation steps and linearly overlay the clean wavefield with random noise and downsampling simulations, resulting in a low-quality wavefield depicted inFig.9c.Figure9d provides a side view of this low-quality wavefield inFig.9c. In the following step, we employ the SR3 optimization algorithm, as proposed in our method, to enhance the quality of the suboptimal wavefield, culminating in a superior version shown inFig.9e. As can be seen, the continuity of the wavefield improves progressively from low to high frequency, effectively denoising and supplementing the matrix for reconstruction.Figure9f offers a side view of this improved wavefield inFig.9e.

Particularly noteworthy inFig.9 parts c and9d are the low-quality seismic wavefield data contaminated with noise. In these two images, the continuity of the wavefield information is significantly disrupted, mainly due to the sparsification, which has a severe impact on the lateral continuity of the wavefield data. These low-quality issues undoubtedly present considerable challenges for the seismic data inversion process. In contrast,Fig.9 parts e and f show the wavefield data optimized using the SR3 algorithm. These images demonstrate the remarkable effectiveness of the proposed algorithm in interpolating and denoising low-frequency wavefield data, particularly in compensating for the gaps and discontinuities caused by data sparsification. The results show significant enhancement in the lateral continuity of the data, as is particularly evident inFig.9f. In this study, subsequent inversion processes will be based on the untreated wavefield data and the wavefield data optimized with the SR3 algorithm. The optimized wavefield data will provide more abundant and precise information for the subsequent processing steps, laying a solid foundation for high-precision imaging inversion.

In addition, we also give the side view of the simulated subsampled wavefield as inFig.10, whereFig.10a is the clean wavefield at 1 Hz,Fig.10b is the low-mass wavefield obtained by simulated downsampling, andFig.10c is the wavefield optimized by SR3 for the low-mass wavefield, and our proposed algorithm significantly improves the continuity and uniformity of the wavefield.Figure10 parts d–f are the case at 2 Hz, and the conclusion is consistent.

In applying FWI, low and ultra-low-frequency data play a crucial role. The significance lies in that FWI relies on accurate and comprehensive frequency content to reconstruct the subsurface velocity structure. High-quality, low-frequency data enhance the precision of the inversion process and are crucial in addressing the challenges of low-frequency paucity and over-dependence on the initial model in FWI.Figure10 illustrates low-frequency seismic wavefield data impacted by noise and sampling issues, where the characteristics of data sparsification and fragmentation significantly reduce its overall quality. The decline in wavefield quality affects the data reliability and negatively affects the subsequent inversion process and interpretation. In this context, the proposed algorithm effectively enhances the quality of low-frequency wavefield data by efficiently interpolating and reconstructing missing or damaged information. Improving the SR3 strengthens the continuity and integrity of the data and compensates for lost information caused by data sparsification and noise introduction. Therefore, the application of the SR3 algorithm has the potential to enhance the quality of inversion results in FWI, leading to more accurate velocity models and providing a more reliable foundation for the detailed imaging and interpretation of complex geological structures.

Additionally, eigenvalues and eigenvectors are crucial for the dimensionality reduction of the wavefield matrix, as they can capture the most significant directions of variation in the matrix, making them extremely useful for interpolation and reconstruction algorithms. We present the results of singular value decomposition applied to the source–receiver wavefield at different frequencies inFig.11.

As can be seen, there are no extreme size differences in the eigenvalues of the wavefield matrix, but the eigenvalues decrease more rapidly at higher frequencies. Consequently, the reconstruction and interpolation of low-frequency matrices become particularly important. In other words, optimizing low-frequency wavefield matrices can significantly enhance the accuracy and precision of the inversion.

Following these, we present the inversion results of the multi-scale frequency-domain FWI, as shown in Fig.12.

The frequency-domain multi-scale inversion strategy in FWI is widely used in geophysical exploration. The key to this strategy is progressively using data across various frequencies to construct an underground model. Multi-scale inversion typically starts with low-frequency data, establishing a rough model of subsurface velocities to capture large-scale geological structures. Subsequently, the frequency of the data is gradually increased, building on the initial low-frequency model. As the frequency increases, the inversion provides finer resolution details. This multi-scale inversion process is iterative, with each step further refining the model based on the preceding one.Figure12 parts a1–a8 represent the broadband inversion results based on the conventional Tikhonov regularization FWI at frequencies of 1.20, 2.99, 3.58, 5.16, 7.43, 10.70, 12.84, and 15.41 Hz, respectively. The forward modelling data were linearly added with 10 dB of random background noise and missing traces, with the processing flow and the results shown inFigs.8 and9. In addition,Fig.12 parts b1–b8 represent the FWI inversion results obtained from the wavefield processed by SR3, with the same frequency range. To showcase the difference in performance, we calculated the differences between the velocity model obtained by the conventional algorithm and the actual velocity model (Fig.12c1–c8), as well as the differences between the FWI inversion results based on SR3 and the actual velocity model (Fig.12d1–d8), both using the same colour map. The 2004 BP model is a sharp contrast velocity model, where the velocity of the salt body part is significantly different from the surrounding velocity, making delineating the salt body contour a challenging aspect of inversion. Additionally, the structure of the water channels on both sides also poses a challenge for inversion. FromFig.12, we can see that the conventional FWI algorithm, affected by background noise and downsampling, has indistinct inversion results for the deep salt column structure and exhibits inversion anomalies for the structures of the water channels on both sides, especially an error in the structure on the left side. By contrast, the FWI optimized by the SR3 algorithm performs excellently; the structure of the water channel on the left is correct and precise, the continuity of the two salt columns is better, and its velocity difference model is cleaner compared with the true velocity. Meanwhile, the velocity difference model of the conventional algorithm, compared with the actual velocity, shows more distinct salt body contours and water channel structures, indicating it is not as clean in its processing.Figure13 shows the convergence rates of the misfit error and model error for both methods, indicating that the convergence rate of SR3 (represented by the red solid line) is faster than that of the conventional FWI (represented by the blue solid line).

Finally, we compare the one-dimensional velocity models of the two methods across six data sets, both horizontally and vertically.Figure14 presents the comparison of velocities at six different x locations. Anomalies and missing traces can lead to apparent errors in the inversion of the conventional FWI, such as a velocity anomaly at x=17.88 km. However, the FWI based on the SR3 algorithm can effectively optimize and reconstruct the forward-modelled data, thereby correcting these velocity anomalies, which display a standard velocity curve unimpacted by missing traces. Therefore, this demonstrates the superiority of our proposed algorithm over the conventional one.

The comparison of lateral velocities is one of the more challenging aspects, as shown inFig.15.

Owing to the propagation characteristics of seismic waves, the lateral velocity comparison of conventional FWI is often of inferior quality. However, the lateral velocity comparison of FWI, based on SR3 preprocessing, closely matches the actual velocity model. This match reflects sufficient velocity compensation at the shallow part of the model and better control of anomalies at deeper parts. By contrast, the lateral velocity of the conventional algorithm, especially at the deeper part, still exhibits noticeable velocity anomalies due to the influence of noise and missing traces.

To evaluate the performance of our proposed algorithm under lower SNR conditions, we also conducted a set of comparative tests in the presence of 5 dB of random noise. Figure16 illustrates these test results.Figure16a shows the clean wavefield information at 2 Hz;Fig.16b depicts the information of 5 dB random noise;Fig.16c represents the wavefield after adding 5 dB random noise;Fig.16d demonstrates the simulated missing-trace matrix, leading toFig.16e, which is the subsampled wavefield matrix with 5 dB noise and missing traces. It is evident that 5 dB of noise significantly disrupts the wavefield, making high-resolution inversion exceedingly challenging under these conditions.

Subsequently, we present the test results as illustrated inFig.17.Figure17 parts a1–a4 depict the inversion results at frequencies 2.99, 3.58, 5.16, and 7.43 Hz using conventional algorithms, whileFig.17 parts b1–b4 show the inversion results at the same frequencies using our proposed algorithm. These results indicate that noise and missing traces significantly affect the inversion outcomes, especially the artefacts resulting from extremely low SNR (5 dB of random noise). Moreover, the structures of the bottom salt bodies appear blurred in the conventional algorithm results. By contrast, the results obtained with our proposed algorithm are notably superior, particularly evident in the higher resolution and more precise structure of the two salt columns in the deeper part of the model. The results highlight the effectiveness of our algorithm in wavefield interpolation reconstruction and noise reduction, efficiently compensating for the inversion anomalies caused by missing information.

Finally, we conducted a one-dimensional velocity analysis and comparison of the inversion results, as shown inFig.18. The comparison reveals that the red solid line, representing our proposed algorithm, aligns more closely with the black line, indicating the actual velocity, especially in deep ranges with significant velocity differences. Our proposed algorithm compensates for the velocity discrepancies evident in conventional methods, suggesting that our improved algorithm more accurately depicts the critical details of the model.

3.3. SEG/EAGE overthrust model

We proceeded to evaluate the performance of our refined approach using the SEG/EAGE overthrust model. This model is a widely adopted geological construct in seismic exploration, predominantly employed to characterize geological formations with overthrust structures. The primary rationale for selecting this model lies in its distinct representativeness as a geological model. The complexity of the inversion process for the overthrust model is primarily concentrated on accurately inverting the velocity of high-velocity layers obscured by multiple overlying strata. The inversion difficulty for the overthrust model involves precisely delineating the velocities of various strata within the covering layers and characterizing the features of the high-velocity basem*nt layer. The P-wave velocity of the SEG/EAGE overthrust model, as depicted inFig.19a, is widely used for testing and validating various geophysical algorithms (Yuan et al. 2015). The overthrust model exhibits numerous complex structures, including many thrust faults and fluvial sediments, which facilitates our exploration of the improved FWI's performance and representation under intricate geological conditions.Figure19b indicates the initial velocity model.

We showcase the SEG/EAGE overthrust model's clean wavefield at 3 Hz, as illustrated inFig.20a. Accompanying this,Fig.20b displays the same wavefield with an addition of 10 dB random noise, providing a comparative analysis of the wavefield under different conditions. The noise-included wavefield,Fig.20c, is obtained by linearly superimposing the random noise onto the clean wavefield. Subsequently, we simulate a subsampled wavefield by uniformly downsampling the matrix, resulting inFig.20d. Finally, we superimpose this onto the clean wavefield to produce a subsampled wavefield matrix with 10 dB of noise and a missing trace, as represented inFig.20e. The inversion for the SEG/EAGE overthrust model is based precisely on this, and the SR3 algorithm is also used to optimize the low-quality wavefield shown inFig.20e.

In the context of FWI, when data are missing, the quantity and information content of the seismic wavefield data decreases. This reduction mainly affects certain aspects of the wavefield, notably attenuating and eliminating some low-frequency information. Low-frequency information often plays a pivotal role in achieving high-precision imaging. The primary purpose of this treatment is to simulate real scenarios where the quality of seismic data collection is compromised due to limitations in the data acquisition system or environmental constraints. Consequently, this approach facilitates the evaluation of both conventional and proposed algorithms in terms of their effectiveness and impact on the inversion results, mainly when dealing with low-quality, low-frequency data.

Similarly, we present a 3D stack of the source–receiver wavefield, as depicted inFig.21a, representing the clean wavefield from 0.5 to 5 Hz.Figure21b offers a clearer side view, obtained by adjusting the perspective from (155, 20) to (105, 5) inFig.21a.Figure21c demonstrates the subsampled wavefield matrix stack with 10 dB of noise and the missing trace, obtained by following the two-step processing procedure fromFig.20. We employ the proposed SR3 to optimize it in response, resulting in a superior wavefield stack, as shown inFig.21e. Our proposed algorithm yields commendable matrix interpolation and reconstruction results in both the 2004 BP model and the overthrust model, exhibiting improved continuity in the wavefield stack, superior noise suppression, and more apparent wavefield details.

Then, we present the inversion results of the multi-scale FWI.Figure22 shows images (Fig.22a1–a8) depicting the inversion results obtained using conventional Tikhonov regularization FWI, whereasFig.22b1–b8 exhibit the inversion results procured via FWI based on SR3 optimization. The results correspond to the inversions at frequencies of 2.15, 3.09, 5.35, 7.70, 9.24, 11.09, 13.31, and 15.97 Hz, respectively. From the results, it is evident that the enhanced algorithm distinctly outperforms in handling complex structures. It not only delineates shallow layered structures more clearly, enhancing their contours, but also sketches and outlines velocity disparities more effectively at depth, avoiding the blurry deep structures often obtained from conventional algorithms. To emphasize the differences, we calculated the discrepancies between the two algorithms and the actual velocity model, as shown inFig.22c1–c8 and d1–d8.

Under the same colourmap conditions, the discrepancy between the inversion results of the improved algorithm and the actual velocity model is visibly superior, exhibiting indistinct layering of the discrepancy, uniform velocity differences between layers, and no notable velocity anomalies at the bottom. Moreover, the improved algorithm ensures the stability of FWI in the high-frequency range without inversion anomalies caused by the spatial aliasing effect and always maintains stability and accuracy. These indeed validate the advantages and efficacy of the improved algorithm. We compare the misfit and model errors of the two algorithms inFig.23, in which the model error is given byEquation(13).

Last, we compare the one-dimensional velocity models derived from the two algorithms at various x- and y-positions, representing the vertical and horizontal dimensions. InFig.24, velocity model discrepancies at six different x-positions are presented, at 3.13, 4.17, 7.77, 10.97, 17.67, and 18.33 km precisely. From the curves, it is apparent that the velocity model generated by the improved algorithm conforms more closely to the actual velocity at many locations with substantial velocity disparities. By contrast, the conventional algorithm often performs poorly in in-depth fitting, sometimes resulting in inversion errors.

Figure25 presents the more challenging cross-sectional velocity comparisons. Owing to the propagation nature of seismic waves and the constraints of data acquisition, seismic waves emanate spherically from the source, i.e. they cover a greater vertical distance compared to the horizontal distance. Moreover, in practical seismic exploration, we typically cannot place a receiver at every point on the ground due to economic and implementation restrictions, resulting in a non-uniform data sampling on the surface, affecting the accuracy of lateral velocity determination. Despite these challenges, the velocity model obtained from the improved algorithm inFig.25 still fits the actual velocity model better. We present data for six groups, respectively: (i) y=0.40 km; (ii) y=0.81 km; (iii) y=1.66 km; (iv) y=3.07 km; (v) y=3.44 km; and (vi) y=3.68 km. In every group, the FWI based on SR3 optimization demonstrates superior performance. The conventional method is affected by the inferior quality of the wavefield data, which includes noise and missing traces, and fails to meet the challenges under extreme conditions. In contrast, our proposed new algorithm successfully addresses these challenges, thus demonstrating its effectiveness, robustness, and practicality.

4. Discussions

FWI, as a research direction that has received much attention, many scholars have made outstanding contributions to this field. However, the truth is that even though it has been proposed for several decades, it has not been widely applied due to certain shortcomings and defects. As an algorithm composed of multiple steps, in FWI, good forward modelling does not necessarily depend on good inversion, but good inversion certainly depends on good forward modelling (Warner and Guasch 2016). Therefore, the results of forward modelling are crucial for FWI. In practical production, it is well known that external factors can create many troubles and difficulties for data collection and quality. Interference from noise, limitations on the azimuth angle, and other factors can all gradually reduce the data quality, reducing the accuracy of forward modelling and increasing its difficulty (Yang et al. 2010). On the other hand, due to the large scale of seismic exploration, FWI, which relies on grid methods, highly demands computer memory (Brossier 2011). Under the current situation where the theory of grid methods and computer memory are close to their limit, improving the quality of the forward-modelled wavefield and laying a better foundation for inversion are urgent problems to be solved in FWI (Takekawa and Mikada 2018).

Regularized regression problems have widespread applications in statistical modelling, signal processing, and machine learning, particularly sparse regression, which plays a significant role in scientific model discovery, including applications in compressed sensing, variable selection, and high-dimensional analysis (Zheng et al. 2019). This paper proposes a preprocessing step before optimization, the SR3 algorithm. The SR3 algorithm can adequately complete and denoise seismic wavefields without affecting the overall computational efficiency and memory consumption through the empowerment of an auxiliary matrix and the idea of compressing singular value space to speed up the convergence of the regression algorithm. We conducted tests in hom*ogeneous media and two types of artificially synthesized data, and the results aligned with our expectations.

The reason for proposing such an approach is precisely to solve the numerous problems encountered in the practical application of FWI. For example, FWI is often affected by data noise, so denoising effectively becomes a problem. The proposed algorithm can effectively denoise from low to high frequencies, allowing FWI to achieve denoising effects without relying on regularization algorithms. Furthermore, the concept of full azimuth is often mentioned in FWI. Although this is an excellent system to improve FWI resolution, we must consider the issues of cost-effectiveness and practicality. Full azimuth necessarily brings higher time and financial costs, contrary to the original intention of proposing FWI. The original reason for proposing FWI was to reduce time and financial costs through more innovative algorithms. Therefore, updating algorithms within FWI may be a more cost-effective way. Thus, completing or reconstructing low-quality data through regression algorithm-like methods has an excellent research background and significance. The SR3 algorithm proposed in this paper is a higher-level sparse regression algorithm, mainly reflected in three aspects. (i) For the denoising problem of FWI, the SR3 algorithm has already been completed in preprocessing, so it is not necessary to carry it out during the minimization process, which avoids the computational pressure brought about by the need to introduce the model increment twice during the minimization process. (iii) The SR3 algorithm is also essentially a multi-constraint optimization problem, so a composite regularization can be used, making the algorithm very flexible and robust. (iii) Compared with other regression algorithms, the SR3 algorithm tightens the singular value space, as shown inFig.2, significantly speeding up computational efficiency without bringing computational pressure to the overall FWI.

However, although the SR3 algorithm is more advanced than other methods, selecting parameters is as complex as the other nonlinear algorithms. We do not object to using heuristic parameter selection methods, but a more logical and reasonable hyperparameter selection method should be considered. Although cross-validation and the L-curve method, among others, are both excellent solutions, choosing the appropriate hyperparameters without increasing computational time and without the need for repeated computations is one of the aspects worthy of more profound future research.

5. Conclusions

Our research introduced a more comprehensive approach to broadband FWI, integrating a preprocessing procedure that handled frequency-domain source–receiver wavefields. This process was achieved by applying the SR3 algorithm. We proposed a method to confront several of the most familiar challenges in seismic data processing: denoising, data reconstruction, and resistance to spatial aliasing effects. In the numerical part, we initiated our analysis by comparing the treatment effects of two hom*ogeneous media cases: single-layer and double-layer media. We specifically highlighted the variance between the wavefield before and after processing and the numerical solution in the high-frequency part. Subsequently, our proposed algorithm was applied to two sets of credible benchmarks widely used for further testing. The results from the forward modelling and inversion processes compellingly demonstrated the efficacy of our algorithm, with prominent features including significant optimization of noisy data and more comprehensive inversion of model details. In practical operations, the quality of seismic data encountered often parallels the quality of the simulation forward modelling data used in our research, which is typically low. Handling such information becomes extremely difficult in these highly challenging scenarios, especially when dealing with low-quality seismic data. Consequently, the innovative optimization algorithm introduced in our study emerges as a promising and more logical methodology for the expanded industrial application of FWI. The proposed approach opens up a promising avenue for significantly improved seismic inversion results, especially in environments where the raw data quality does not meet our requirements.

Conflict of interest statement. The authors declare that they have no conflict of interest.

Funding

This study is supported by Kansai Research Foundation for Technology Promotion (grant no. 2024P001).

Data availability

The model used in this article is available in Zenodo at [https://zenodo.org/]. The code underlying this paper will be shared on reasonable request to the corresponding author.

References

© The Author(s) 2024. Published by Oxford University Press on behalf of the SINOPEC Geophysical Research Institute Co., Ltd.