Optimal Inversion Method Based on Joint Waveform Inversion and Least Squares Reverse Time Migration

Joint Waveform inversion (JWI) uses the results of reflection waveform inversion (RWI) as the initial model for full waveform inversion (FWI). Compared with the FWI, the JWI method can obtain more information about the structure of the subsurface medium. The reason is that the reflected waveform inversion can invert the long wavelength component in the middle and deep areas. In JWI, reflected waveform inversion is used to obtain the reflected wave information in the simulation record by demigration, which is computationally more expensive than FWI; the least squares reverse time migration (LSRTM) also obtains the reflected wave information in the simulated record by demigration. In order to effectively use the reflected wave information brought by the high computational amount of reverse migration in JWI, this paper proposes a simultaneous inversion method of JWI and LSRTM (JWI-LSRTM). This method can simultaneously perform an iterative update of the subsurface medium velocity of JWI and the migration imaging of LSRTM, which improves the calculation data utilization rate of each forward and inversion process. In the model test, the effectiveness of the method is proved.Reflected waveform inversion (RWI) is a technique that obtains reflected wave information from simulated records through reverse migration. It can invert long wavelength components in the middle and deep layers [1-3], but its computational cost is higher than full waveform inversion (FWI). FWI is a high-precision inversion method [4-6] with the potential to provide accurate models of subsurface media parameters. Joint waveform inversion (JWI) [7] leverages the results of reflection waveform inversion (RWI) as an initial model for FWI. Compared to FWI, JWI can retrieve more information about the subsurface media structure [8]. This is because RWI can obtain long wavelength components in the middle and deep layers [9, 10]. Ren (2019) proposed an adaptive JWI method that automatically switches between RWI and FWI by adjusting the weight factor with the number of iterations and allowable errors, without manually pausing the switch [11]. This approach addresses the limitations of traditional waveform inversion methods and improves the efficiency and accuracy of subsurface media modeling.The LSRTM [12, 13] is based on the Born approximation, and the reflection coefficient is solved by many iterations with the known background velocity. LSRTM is also used to obtain the reflected wave information from simulated records.In the above three methods, the simulated data and observed data are inverted using the generalized least squares method to obtain the corresponding gradient. However, the computational cost of forward modeling the wave equation and reverse migration is substantial, accounting for at least 90% of the total computation time in these three methods. As a result, the inversion cycle is often prolonged [14, 15]. Now there are a variety of programming techniques (MPI, Openmp, Openacc, GPU) or optimization methods (L-BFGS method, pseudo-Newton direction correction mixed conjugate gradient method, etc.) can shorten this cycle, and there are also ways to improve the utilization of calculated data (FWI imaging, RWI, and LSRTM joint inversion) [16, 17].In order to make more effective use of the reflected wave information and wave equation forward simulation data brought by the high-computational demigration in JWI, a JWI and LSRTM synchronous inversion method is proposed in this paper. This method can simultaneously update the velocity of subsurface media in the JWI and the migration profile of the LSRTM, which improves the utilization rate of the calculated data in the forward and inversion processes. The effectiveness of the proposed method is proved by model test and noisy data.This paper takes acoustic, which is easier to realize, as an example. Here are the wave equation and demigration equation of sound wave:where x and z are the space coordinates, t is the time, p0 is the background field, p is the disturbance field, and vp is the p-wave velocity. In equation (2), dcalxr,t;xs is the simulated data, dcalrefxr,t;xs is the scattered data, dobsxr,t;xs is the observed data, and m is the reflection coefficient, whose expression is: mimage(x)=Δv(x)vp(x) , the scattered data can be calculated through equation (2).The FWI obtains the subsurface model parameters by minimizing the error of the observed data and the simulated data. Its adjoint equation and gradient formula can be obtained by the adjoint method. The following objective function can be obtained:where, T is the maximum recording time dcal and m is the subsurface medium parameter, which are respectively the simulated data (obtained by equation dobs(1)) and the observed data. Replace the source term of the adjoint equation with the data residuals to obtain the corresponding back propagation data (the acoustic adjoint equation is consistent with equation (1) apart from the initial conditions). The gradient formula in FWI can be obtained by the adjoint method:where dcal* is the value of the backpropagation wavefield. When m is the velocity, the v iterative formula is as follows [18-20]:Among them, αk is the first step and dk is the search direction of step 1. To sum up, the process of FWI [21-23] can be summarized as follows: 1. Obtain simulated data through equation (1); 2. Obtain the error between the observed data and the simulated data obtained by Formula (1), and obtain the sum of the residuals of the two data by formula (3); 3. Replace the source term of the adjoint equation with the data residuals to obtain the corresponding backpropagation data, and obtain the dcal* corresponding gradient according to formula (4); 4. Update the vp according to formula (5).The LSRTM can obtain the subsurface imaging result by minimizing the reflected wave information error of observed data and simulated data. Its minimization target functional is as follows:where dcalref is reflected wave information of simulated data (obtained by equation (1) and (2)) and reflected wave information of dobsref is observed data (obtained by windowing observed data or reverse migration of model data) respectively. Replace the source term of the adjoint equation with the data residual to obtain the corresponding backpropagation data dcalref* (the acoustic adjoint equation is consistent with equations (1) and (2)) [23]. The gradient formula in mimage LSRTM can be obtained by the adjoint method. The gradient formula of LSRTM under the velocity stress equation of time order is as follows:When m is velocity, the mimage iterative formula is as follows:To sum up, the process of LSRTM can be summarized as follows: 1. Obtain the reflected wave information of simulated data through equations (1) and (2); 2. Obtain the error between the reflected wave information of observed data and the reflected wave information of simulated data obtained by formulas (1) and (2), and obtain the sum of the residual difference of the two data by formula (6); 3. Replace the source term of the adjoint equation with the data residuals to obtain the corresponding backpropagation data, and obtain the dcalref* corresponding gradient according to formula (7); 4. Update the mimage according to formula (8).The process of reflection wave inversion is similar to LSRTM, except that RWI is updating the velocity parameter (that is, the physical meaning of the gradient is different, and its iterative formula is the same as that of FWI). The gradient formula in RWI can be obtained by the adjoint method [24]. The gradient formula of RWI under the velocity stress equation of time order is as follows:The process of reflection waveform inversion can be summarized as follows: 1. Obtain the reflected wave information of simulated data through equations (1) and (2); 2. Obtain the error between the reflected wave information of observed data and the reflected wave information of simulated data calculated by formulas (1) and (2), and obtain the sum of the residual difference of the two data by formula (6); 3. Replace the source term of the adjoint equation with the data residuals to obtain the corresponding backpropagation data (⁠dcalref* , dcal*⁠), and obtain the corresponding gradient according to formula (9); 4. Update the vp according to formula (5).The objective function of the adaptive joint inversion can be expressed as follows [9]:where witer is the weight operator, which is related to setting the number of iterations niter and setting the allowable residuals.When the number of iterations is less than the set number of iterations or the objective function value is greater than the set allowable error, the witer is 1, namely RWI. Otherwise, the witer is 0, or FWI. This paper gives the following new formula of gradient:As can be seen from the above, the key of the above methods lies in different wavefields obtained from the forward and reverse migration when the gradient is obtained. Among them, the wavefield used in the gradient inversion of reflected waveform cover all the wavefields used in the gradient formulas of FWI and LSRTM. Therefore, when RWI is carried out, LSRTM operation is not carried out at the same time, which will lead to the waste of computing resources in forward and demigration (forward and demigration calculation amount is very large). Similarly, when LSRTM is carried out, FWI is not carried out at the same time, which will lead to the waste of computing resources.In order to make more effective use of the reflected wave information and wave equation forward simulation data brought by the high-computation demigration in JWI, a JWI-LSRTM synchronous inversion method is proposed in this paper. This method can simultaneously update the subsurface media velocity of the JWI and the migration profile of the LSRTM, so as to improve the utilization rate of the calculated data in the forward and inversion processes. The specific process of this method is described as follows: in the first stage, synchronous inversion of LSRTM and RWI is carried out according to formulas (1), (10), (11), (13), and (14); in the second stage, joint inversion with FWI is carried out according to formulas (1), (10), (11), (13), and (14). The specific objective function strategy is the same as formula (10). Since the iteration efficiency of LSRTM and RWI is higher than that of FWI, only the FWI inversion is carried out when the synchronous inversion of LSRTM and RWI reaches the residual threshold or the number of iterations. In this case, there is no need for demigration in the second stage, which can reduce the computation cost. LSRTM and RWI can achieve the corresponding effect. Therefore, the forward equation of the wave equation needs to be used all the time, while the strategy of the demigration formula is as follows:That is, demigration is carried out in the first stage and no reverse migration is carried out in the second stage. The corresponding gradient formula is as follows:According to the number of iterations or the allowable error, the velocity update of RWI and FWI is switched. The gradient of LSRTM also determines whether to update the imaging according to the number of iterations or the allowable error. To sum up, the specific process is shown in the Figure 1 below:Next, we test the correctness and feasibility of our method in this paper. The Canada Overthrust model (Figure 2) has a size of 8200 × 3500 m and a spatial interval of 10 m. The source is a Ricker wavelet with a dominant frequency of 20 Hz, and the recorded traces are sampled in time intervals of 1 ms for a total time of 2.20 s. There was a total of 82 shots, with 100 m interval; the longitudinal coordinate of the shot is 10 m. The total number of receivers is 820, separated by 10 m, and the longitudinal coordinate of the receivers is 0 m. We performed 40 radius smoothing on Figure 2(a) to obtain Figure 2(b); we use Figure 2(b) as the initial model. We add Gaussian noise (S/N=20, S/N is signal-to-noise ratio) to the seismic records obtained from the model of ( Figures 2 and 3).In Figure 4, the number of iterations is 20. The iterative results of the LSRTM+JWI method are the best, but there are errors due to noise interference in the results. We compare the 20th iteration of the single-channel Vp at x=4000m in Figure 5. From Figures, the result of LSRTM+JWI method has better Vp accuracy. In summary, it is proved that the proposed method has relatively higher velocity inversion accuracy in addition to higher data utilization rate.Figure 6 shows the LSRTM results after Laplacian filtering. All three have relatively poor performance due to the same Gaussian noise shot data. However, the interface of the LSRTM+JWI method is more pronounced and continuous compared to the other two methods. We extract a single line of image results in Figure 6 for comparison. As can be seen from Figure 7, compared with conventional LSRTM, the proposed method has a more accurate velocity reduction effect. The accuracy of the results of the LSRTM+RWI method is basically the same compared to that of the LSRTM+JWI method.This paper also uses the field data equivalent model for testing. The Field data equivalent model has a certain complexity. The size of the Field data equivalent model (Figure 8(a)) is 8000 × 4000m, and the spatial sampling interval is 10 m. The source is a Riker wavelet with a main frequency of 25 Hz, with a time sampling interval of 1 ms and a total sampling time of 2.0 s. A total of 80 shots were fired, with an interval of 100 m; A total of 800 geophones, spaced 10 m apart. Figure 8 (b) is obtained by applying 30 layers of Gaussian smoothing to Figure 8 (a), which is used as the initial model in this paper. Figure 8 (c)is the initial image obtained through RTM (i.e., the initial mimage⁠.The gradient optimization method used in the test in this paper is the hybrid conjugate gradient method with pseudo-Newton direction correction, which has higher iteration stability and iteration efficiency than L-BFGS-C. Because of the high iteration efficiency of this method, and for the convenience of comparison, the inversion strategy is as follows: the nonjoint synchronous class method has 20 iterations in total; The first stage of the method in this paper (RWI+LSRTM synchronous inversion) was iterated for 6 times (according to experience, the residual percentage of the sixth iteration generally decreased to 1%) or less than 1%, and the second stage was iterated for at least 14 times, that is, the 20th iteration was stopped. Such a strategy can ensure that the number of iterations of the two is consistent, so as to facilitate the comparison of the final effect. In this paper, the two-point interpolation method is used to calculate the step length, which has higher convergence stability and less computation than the three-point interpolation method.By comparing Figure 9(c), Figure 10(b), and Figure 11(b), it can be seen that compared with conventional FWI, the proposed method has more accurate velocity inversion results in a deeper layer below 2 km, and the velocity boundary at the velocity reversal layer and strong reflection layer is more obvious. The reason is that the synchronous inversion of RWI+LSRTM in the first stage provides an initial velocity model containing reflected wave information of subsurface interface, which improves the utilization rate of reflected wave information in FWI. The effect of the RWI+FWI method differs less compared to that of the LSRTM+JWI method, we extract a single line of velocity for comparison. As can be seen from Figure 12, compared with conventional FWI and RWI+FWI, the proposed method has a more accurate velocity effect. In summary, it is proved that the proposed method has relatively higher velocity inversion accuracy in addition to higher data utilization rate.By comparing Figure 9(b), Figure 10(a), and Figure 11(a), it can be seen that compared with conventional LSRTM, the proposed method has higher resolution, more uniform longitudinal energy at the velocity reversal layer and strong reflection layer, and at the right high-velocity body, and better axial continuity at the strong reflection layer. The effect of the LSRTM+RWI method differs less compared to that of the LSRTM+JWI method, we extract a single line of mimage for comparison. As can be seen from Figure 13, compared with conventional LSRTM, the proposed method has a more accurate velocity reduction effect. The accuracy of the results of the LSRTM+RWI method is basically the same compared to that of the LSRTM+JWI method. The reason for this is that the residuals are difficult to decrease after reaching a certain threshold. The effects of these two methods are basically the same.In Figures 14 and 15, although the conventional method has a very high residual convergence rate, the proposed method has a higher residual convergence rate. The reason is that the RWI can invert the long wavelength components in the middle and deep layers, and provide more information about subsurface reflected waves for FWI and LSRTM. In the 8th iteration, based on the judgment criteria, the RWI method was switched to FWI for subsequent iterations, resulting in a turning point in the convergence curve. In conclusion, in addition to higher data utilization rate, the proposed method has relatively higher residual convergence rate and inversion accuracy.In order to make effective use of the reflected wave information and wave equation forward simulation data caused by the high-computational demigration in JWI, a method of JWI-LSRTM synchronous inversion is proposed in this paper. This method can simultaneously update the velocity of subsurface media and the migration imaging of LSRTM, and improve the utilization rate of calculated data in the process of forward and inversion. In the model test, the effectiveness of this method is proved. At the same time, the proposed method can obtain more information on subsurface media structure, with better inversion accuracy and higher residual convergence rate. The reason is that reflected waveform inversion can invert long wavelength components in the middle and deep layers, and provide more reflected wave information for FWI and LSRTM.The first test data is the Canadian overthrust fault model, which is publicly available and the shot data wasdata were obtained through our own forward modeling. And the second model is an equivalent model of real data, which is also self simulated, and there is no problem of referencing data from others.No potential competing interests were reported by the authors.The authors express a sincere thanks to the SWPI Research Group of China University of Petroleum (East China) and National Key Laboratory of Deep Oil and Gas for their support. This work was supported by the Natural science foundation of Shandong Province (grant no. ZR2022MD037).