The two-dimensional acoustic equation in the frequency domain has the general form of ∂ 2 ∂ x 2 + ∂ 2 ∂ z 2 + ω 2 v x , z 2 a x , z , ω = s x , z , ω , (1) where ω is the angular frequency, v is the P-wave velocity, a is the pressure field, s represents the seismic source, and x and z represent the horizontal and vertical coordinates, respectively. The matrix equation with seismic source related to Equation 1 is given by F A = S , (2) where F represents the Helmholtz operator matrix, A is the pressure field matrix, and S is the source matrix. The LU factorization techniques are generally used to solve Equation 2 by F = L U , (3) Y = L − 1 S , (4) A = U − 1 Y . (5)

We can simply express Equation 2 as A = F − 1 S . (6)

As illustrated in Figure 1, according to the frequency-domain Kirchhoff integral (Berkhout, 1984), the wavefield at a point (P ₁) situated outside of the closed surface O can be expressed as a P 1 , ω = − 1 4 π ∬ O ∂ a ∂ n + j k ∂ r 1 ∂ n a + 1 r 1 ∂ r 1 ∂ n a G d O , (7) where G is Green’s function related to r ₁ or r ₂, k = ω/v and j is the imaginary unit. In this equation, the variable “a” on the right-hand side represents the wavefield on the closed surface O. We can calculate wavefields at other spatial locations based on the known wavefields on the closed surface O and their directional derivatives. The wavefield at a point (P ₂) within the closed surface O can be expressed as a P 2 , ω = ∭ Ω s ω G d Ω + 1 4 π ∬ O ∂ a ∂ n + j k ∂ r 2 ∂ n a + 1 r 2 ∂ r 2 ∂ n a G d O . (8)

The volume integral term represents the contribution of the source distributed in the volume Ω to the wavefield at P ₂. The first term of the equation can be regarded as the background wavefield, and the second term is the primary or multiple scattering wavefield. If O is infinite, the energy on it has little influence for the wavefield. Therefore, we can ignore the second term of Equation 8 and then a P 2 , ω = ∭ Ω s ω G d Ω . (9)

Comparing Equation 9 with Equation 6, we have Equation 10 as F − 1 = ∭ Ω G d Ω . (10)

When subjected to rigid boundary conditions, the Kirchhoff integral, as illustrated in Equation 7, can be reduced to the first type of RS integral (Berkhout, 1984) a P 1 , ω = − 1 2 π ∬ O 1 ∂ a ∂ n G d O . (11)

This equation demonstrates that the wavefield at P ₁ can be calculated using the wavefields on the plane O ₁ and their directional derivative. Huygens’ principle (Born and Wolf, 1999) states that fluctuations on the surface O ₁ can be treated as sources, thereby allowing ∂ a / ∂ n in Equation 11 to be treated as a source. ∭ Ω G d Ω and ∬ O 1 G d O represent the same wave propagator F − 1 for a given velocity model. In Equation 6, A represents the wavefield for the entire space. In Equation 11, a P 1 , ω is the wavefield at P ₁. Accordingly, Equation 11 can be expressed in the following form: A = β F − 1 E , (12) where A is the reconstructed wavefield, β represents the amplitude, E acts as the reconstructed source, and E = ∂ a ∂ n . (13) Where E represents an element of E. Equation 12 is the frequency-domain formula of wavefield reconstruction. As shown in Figure 2, the red lines represent the positions of the reconstructed sources (E ₁ -E _n) and the blue areas represent where the reconstructed wavefields propagate. According to the first type of RS integral (Berkhout, 1984), the wavefields propagating in the white areas above the reconstructed sources are total reflections of the reconstructed wavefields from the reconstructed sources. Unless otherwise specified, the reconstructed wavefield mentioned in the following section of this paper refers to the wavefield below the reconstructed source. After obtaining the original wavefield, we use the original wavefield on the red line to calculate the reconstructed source (such as E ₁) according to Equation 13. Subsequently, the reconstructed source (such as E ₁) is positioned in the corresponding red line position, thereby yielding the reconstructed wavefield (such as C ₁) in accordance with Equation 12. For simplicity, we set a constant on β to calculate the reconstructed wavefield. The reconstructed source can be placed at any desired location, and different reconstructed wavefields (such as C ₁–C _n) can be generated by loading the appropriate reconstructed sources (such as E ₁–E _n).

For a two-dimensional velocity model as shown in Figure 2, the reconstructed sources (E ₁–E _n) calculated with the original wavefields at different depths are loaded to obtain reconstructed wavefields (C ₁–C _n). We sum a proportion of the reconstructed wavefields to obtain the MRW and enhance reflected waves in the MRW. If we load all the reconstructed sources to a model at the same time, according to Equation 12, we have (Yu and Jia, 2021) C = β F − 1 E 1 + E 2 + · · · + E n = C 1 + C 2 + · · · + C n . (14)

The wavefield calculated by Equation 14 is the MRW. According to Figure 2 and Equation 14, The different depth parts of the MRW are the superposition of different amounts of the reconstructed wavefield. For amplitude preservation, the different depth parts of the MRW should be divided by the number of stacking times. The value of n can be freely set according to the specific needs and the size of the velocity model. The reflected waves in the MRW, which can facilitate the update of the deep velocity (Lian et al., 2018; Dong et al., 2018; Yao et al., 2020), are enhanced. The MRW is used in place of the original wavefield to optimize the gradient of the FWI. For noisy data, dot product of the adjoint wavefield with the MRW can mitigate the effect of noise on FWI gradients.

The FWI problem is to minimize the objective function min ⁡ f q = 1 2 ∑ ω ∑ i R i A i − B i 2 , (15) where q is the squared slowness, B _i denotes the recorded data, R _i represents the operator for extracting simulated data from simulated wavefield A _i , and the subscript i indicates the shot-gather number. The gradient of the objective function (Equation 15) is expressed by the adjoint state method (Plessix, 2006; Métivier et al., 2012) as ∇ f = R e ∑ ω ∑ i ω 2 λ i * ω A i ω , (16) where λ _i ^* is the complex conjugation of the adjoint wavefield, and Re denotes the real part of a complex value. According to Taylor series, the exact model update format is Q l + 1 = Q l − H − 1 l ∇ f l , (17) where l represents the number of updates, Q and H are the parameter matrix and the Hessian matrix, respectively. Using l-BFGS method, we approximate the Hessian matrix from Equation 17 to derive Equation 18, which is Q l + 1 = Q l + α l D l , (18) where α denotes the step length, and D is the search direction. In this paper, the line search method for calculating the step length satisfies the strong Wolfe conditions (Nocedal and Wright, 1999) f Q l + α p D l ≤ f Q l + θ 1 α p ∇ f l T D l , (19)

and ∇ f Q l + α p D l T D l ≤ θ 2 ∇ f l T D l , (20) where 0 < θ ₁ < θ ₂ < 1, T denotes the matrix transpose, and p represents the p-th searching iteration.

We used the l-BFGS to explain the FWI utilizing the MRW (i.e., l-BFGS-MRW). In the process of the l-BFGS-MRW inversion, after the original wavefield (Equation 6) of forward propagation is obtained, the MRW is calculated according to Equation 14 and then applied to the optimized gradient by Equation 16. The forward wavefield A _i (ω) in Equation 16 is replaced by the MRW, and other processes are consistent with the original l-BFGS inversion. The specific details of the workflow of the l-BFGS-MRW inversion are shown in Figure 3.

In order to evaluate the inversion results quantitatively, we introduce the residual sum of squares (RSS), defined by Equation 21 as: RSS = ∑ x ∑ z v x , z − v true x , z 2 , (21) where v x , z and v true x , z represent inversion model and the true model respectively. The RSS represents the discrepancy between the inversion results and the true model. A smaller RSS value indicates a greater similarity between the inversion results and the true model.

Figure 4A shows a two-layer velocity model, where the red hexagon represents the location of the seismic source. The solid magenta line denotes a part of a closed surface O with n as its outer normal vector. Note that the solid magenta line also represents the position of the reconstructed source E. The grid spacing was 15 m. Figure 4B shows the original wavefield for a frequency of 10 Hz. After obtaining the original wavefield, we calculated the reconstructed source and the reconstructed wavefield according to Equations 13 and 12, respectively. Figure 4C shows the reconstructed wavefield for 10 Hz, and the wavefield above the reconstructed source is the totally reflected wavefield. To reconstruct the wavefield accurately, the reconstructed source contained a number of grid points in the n-direction. The residual between the reconstructed wavefield and the original wavefield is displayed in Figure 4D. Although the residual wavefield has different spatial characteristics, the energies of all its components are much smaller than those of the original wavefield (Figure 4B). Figure 4E shows the real parts of the original wavefield and the reconstructed wavefield along the dotted white line in Figure 4A. As shown in Figure 4E, the curves of the original wavefield and the reconstructed wavefield are overlapped and their residual is approximately three orders of magnitude smaller than the original wavefield, which indicates that the reconstructed wavefield is almost equal to the original wavefield when their positions are close to the source. When the position of the reconstruction wavefield is far from the source, as shown in Figure 4F, there is a small amplitude difference between the reconstructed wavefield and the original wavefield. Because we just take attention to a part of the closed surface O and ignore the rest of it, some energy, which exist in the original wavefield, is not involved in the reconstruction wavefield; therefore, the reconstructed wavefield and the original wavefield have little amplitude difference in the area far from the source.

In the example of a complex model, we can get the same conclusion as the former. Figure 5A shows a slice of the Marmousi velocity model, where the solid magenta line represents the location of the reconstructed source. The shot was located at distance of 6.9 km and depth of 0.09 km. The grid spacing was 15 m. Figure 5B displays the original wavefield for 10 Hz, and Figure 5C displays the reconstructed wavefield for 10 Hz. The residual between the reconstructed wavefield and the original wavefield is shown in Figure 5D.

Figure 5E illustrates the real parts of the original wavefield and the reconstructed wavefield on the dotted black line in Figure 5A. The reconstructed wavefield in Figure 5E is close to the source; therefore, the residual of the two wavefields is approximately 2 orders of magnitude smaller than the original wavefield, which suggests that the two wavefields are almost the same. Figure 5F shows the real parts of the original wavefield and the reconstructed wavefield on the dashdotted black line in Figure 5A, which indicates that the two wavefields have amplitude differences in the area far from the source. Nevertheless, as Figure 5F indicates, the reconstructed wavefield still preserves most of the energy of the original wavefield and the phases of these two wavefields are almost identical.

In the reconstructed wavefield, most of the original wavefield are completely preserved. Although the reconstructed wavefield far from the source have small energy differences from the original wavefield, their phases are almost identical. These small differences are not enough to interfere with FWI to obtain accurate results. In addition, according to the propagation path of the deep residual wavefields (the deep parts of Figures 4D, 5D), lots of the residual wavefields propagated outside the study area and could not be detected by the geophones distributed on the surface; the deep energy that cannot be received by the geophones is useless for FWI. The MRW is a superposition of different reconstructed wavefields; in the MRW, the components of the wavefield with the same character, including those far from the source, are enhanced by superposition; among them, the reflected waves are more enhanced because they are reconstructed more accurately compared to the other components. In other words, the proportion of the reflected waves in energy is increased.

Updates to the deep part of the model rely primarily on reflected waves (Lian et al., 2018; Dong et al., 2018; Yao et al., 2020). Noise in the data can lead to unwanted adjoint wavefields; in the deep part of the model, the dot product of these unwanted adjoint wavefields with the deep wavefield components that are useless for the FWI gradient is detrimental. The enhanced reflected waves will be highlighted in the gradient calculation to suppress the detrimental dot product. Therefore, the MRW can be used to mitigate the adverse effects of data noise on the FWI gradient. In addition, the difference between the MRW and the original wavefield will avoid FWI from overfitting the noisy data. For these reasons, when we employ the MRW in FWI, we generally obtain a result with high resolution.

For the FWI experiments below, the source was a Ricker wavelet and the dominant frequency was 7 Hz. 436 geophones were evenly distributed on the surface and 109 seismic sources were sequentially placed on the surface. The grid spacing was 25 m and the time interval was 0.5 ms. The recorded data were generated by forward modeling with uniformly distributed noise added. Some of the data are displayed in Figure 6. We made J to represent the ratio of the mean value of the noise energy to the mean value of the signal energy. We took the model in Figure 7A as the true model and employed the initial model shown in Figure 7B. The initial model is derived by applying a high degree of smoothing to the true model. During the inversion (Figures 2, 3), the reconstructed sources (E ₁−E _n) filled the entire current model space to obtain the MRW. Figures 7C, E, G illustrate the l-BFGS results (1–6 Hz) when J is 8.97%, 53.83%, 269.13%, respectively. Figures 7D, F, H show the l-BFGS-MRW results (1–6 Hz) when J is 8.97%, 53.83%, 269.13%, respectively.

In the experiment on data with a modest noise level (J = 8.97%), the results of l-BFGS-MRW and l-BFGS are not significantly different. The residual sum of squares (RSS) of the l-BFGS-MRW result is 1,689, which is slightly smaller than that (2,159) of the l-BFGS result. For a detailed comparison, Figure 8A exhibits the velocities extracted along a line (Distance = 0–7.335 km, Depth = 2 km) in Figures 7A, C, D. The result of l-BFGS-MRW is closer to the true model than that of l-BFGS and there are some large outliers in the l-BFGS result. When J is 53.83%, the RSS of the l-BFGS-MRW result shown in Figure 7F is 2,187, which is much smaller than that (29,533) of the l-BFGS result shown in Figure 7E. The result of l-BFGS-MRW has fewer artifacts and a clearer structure than that of l-BFGS. Figure 8B displays the velocities extracted along a line (Distance = 0–7.335 km, Depth = 1.85 km) in Figures 7A, E, F. The l-BFGS result has many large outliers and the l-BFGS-MRW result is much better than the l-BFGS result. The l-BFGS-MRW inversion is more robust than the l-BFGS inversion. Despite the high noise level, the l-BFGS-MRW inversion can obtain a relatively accurate initial model for high-frequency inversion. When J is large as 269.13%, the result of l-BFGS (Figure 7G) becomes worse than its former result. Because the data is severely damaged by noise, the l-BFGS-MRW inversion can only use less information from the data and restore the low-wavenumber velocity structure as demonstrated in Figure 7H. The RSS of l-BFGS-MRW result is 2,728, while that of l-BFGS result is 100,370.

The data convergence of the 1–6 Hz inversion corresponding to Figures 7C, D is exhibited in Figure 8C. Figure 8D shows the data convergence of the 1–6 Hz inversion corresponding to Figures 7E, F. The values of the objective function following the initial iteration are distinct, resulting in a misalignment of the curves at the starting point in both figures. As illustrated in Figures 8C, D, despite the disparate descending paths of the objective values for the two methods, all objective values converge to a sufficient degree with the same number of iterations. It is challenging to eliminate the introduction of data noise into the calculation of the objective value. While fitting data noise may reduce the objective value, it also carries the risk of diverging the inversion model from the true model. Because the gradients were optimized by the MRW, the models of the l-BFGS-MRW inversion would resist being updated in the direction of data noise guidance. On the contrary, the l-BFGS inversion would fall into local minima generated by data noise, which can result in inversion models that do not approximate the true model. As a consequence of the aforementioned factors, the objective values of l-BFGS-MRW inversion are larger than those of l-BFGS inversion after numerous iterations. Figure 8E illustrates the same phenomenon as the formers. When J equal to 269.13%, the value of the objective function that utilizes the strong noisy data is meaningless. Our gradient construction can hardly reduce the value of the objective function that utilizes the strong noisy data. Therefore, the objective value of the l-BFGS-MRW was not sufficiently reduced (Figure 8E). The l-BFGS inversion overfitting the recorded data is equivalent to fitting some data noise. Despite the objective values being reduced to relatively low values, the results shown in Figures 7E, G exhibit significant discrepancies from the actual models.

The grid spacing of the low-frequency inversion and the high-frequency inversion is 25 m and 15 m, respectively. RSS is calculated based on the velocity at the grid point and all RSS are kept as integers. The low-frequency inversion has a larger grid spacing and fewer grid points than the high-frequency inversion, so the RSS of the high-frequency inversion result may be larger than that of the low-frequency inversion result for the same inversion method. In the experiment of the 7–10 Hz inversion, Figure 9A demonstrates the l-BFGS result (J = 8.97%) with the initial model shown in Figures 7C, 9B illustrates the l-BFGS-MRW result (J = 8.97%) with the initial model shown in Figure 7D. The RSS of the l-BFGS-MRW result is 3,897, which is smaller than that (5,224) of the l-BFGS result. Although J is small, the l-BFGS result shown in Figure 9A has a low resolution in the deep part of the model. In contrast, the l-BFGS-MRW result has a high resolution in the entire model space. Figure 9C displays the l-BFGS result (J = 53.83%) with the initial model shown in Figures 7E, 9D exhibits the l-BFGS-MRW result (J = 53.83%) with the initial model shown in Figure 7F. Their RSS are 37,822 and 167,925 respectively. Although the noise level has increased, the result of l-BFGS-MRW still has a higher resolution than that of l-BFGS. The detailed comparisons for the results of the two methods are shown in Figures 10A, B. In order to facilitate the display, some extremely large outliers in the results of l-BFGS have been manually suppressed. In the high-frequency inversion with noisy data, the l-BFGS-MRW is still robust. The aforementioned RSS have been synthesized into Table 1 for purposes of clarity and ease of reference.

The data convergence for Figures 9A, B is illustrated in Figures 10C, D shows the data convergence for Figures 9C, D. As shown in Figures 8E, 10C, D, the objective values of l-BFGS-MRW converge quickly, while the objective values of l-BFGS can barely converge even in a long iteration epoch. The RSS of the 1–6 Hz and 7–10 Hz inversion results are plotted as a function of J in Figures 10E, F, respectively. In comparison to the l-BFGS inversion, our FWI is not susceptible to data noise.

FWI is a time-consuming inversion technique, in which LU factorization (Equation 3) of the Helmholtz operator matrix accounts for most of the total time. L ⁻¹ and U ⁻¹ were stored after LU factorization for solving the wavefield using Equations 4, 5. In the FWI with l-BFGS-MRW, we could use L ⁻¹ and U ⁻¹ to solve the MRW without conducting an additional LU factorization. We did not need to calculate a series of reconstructed wavefields (Equation 12) and add them to obtain the MRW. According to Equation 14, we loaded a series of reconstruction sources at the same time to obtain the MRW, which costed the same time as calculating a single reconstructed wavefield. In FWI, a suitable search method can help quickly find the appropriate step length to save calculation time. The line search method for the step length employed in this work satisfied the strong Wolfe conditions as shown in Equations 19, 20 (Nocedal and Wright, 1999). During the inversion process, the step length was dynamically adjusted according to the number of the searches and the objective values. In theory, the total time consumption of our FWI is approximately 1.5 times that of l-BFGS FWI at the same iterations. When the noise is strong in the data, in addition to the time for solving the reconstructed wavefield, the l-BFGS-MRW inversion will cost more calculation time than the l-BFGS inversion in searching for the right step size to reduce the data residuals. The time spent on searching the step lengths accounts for most of the additional time. In some cases, such as that shown in Figures 8E, 10C, D, l-BFGS-MRW inversion requires less computational time than l-BFGS inversion because the objective value of l-BFGS-MRW converges much faster than that of l-BFGS, and we can stop the iteration of l-BFGS-MRW when the convergence occurs.