Source-independent efficient wavefield inversion

Full-waveform inversion (FWI) is an effective tool to retrieve a high-resolution subsurface velocity model. The source wavelet accuracy plays an important role in reaching that goal. So we often need to estimate the source function before or within the inversion process. Source estimation requires additional computational cost, and an inaccurate source estimation can hamper the convergence of FWI. We develop a sourceindependent waveform inversion utilizing a recently introduced wavefield reconstruction based method we refer to as efficient wavefield inversion (EWI). In EWI, we essentially reconstruct the wavefield by fitting it to the observed data as well as a wave equation based on iterative Born scattering. However, a wrong source wavelet will induce errors in the reconstructed wavefield, which may lead to a divergence of this optimization problem. We use a convolution-based source-independent misfit function to replace the conventional data fitting term in EWI to formulate a source-independent EWI (SIEWI) objective function. By convolving the observed data with a reference trace from the predicted data and convolving the predicted data with a reference trace from the observed data, the influence of the source wavelet on the optimization is mitigated. In SIEWI, this new formulation is able to mitigate the cycle-skipping issue and the source wavelet uncertainty, simultaneously. We demonstrate those features on the Overthrust model and a modified Marmousi model. Application on a 2D real dataset also shows the effectiveness of the proposed method.


INTRODUCTION
Full-waveform inversion (FWI) is a highly non-linear problem which aims at minimizing the misfit between predicted and observed data to retrieve a high-resolution velocity model (Tarantola 1984;Pratt et al. 1998;Virieux & Operto 2009). The knowledge of the true source function is an important prerequisite to achieve a successful wave equation based inversion. Thus, the source wavelet estimation becomes an essential step for practical FWI applications. It is often, however, a difficult task when dealing with noisy data, imperfect acquisition, or unpreserved amplitude during the data preprocessing (Chi et al. 2018). To get a good estimation of the source wavelet in the field data, we often utilize the direct arrivals in streamer data to extract the source signature. However, in shallow water, ocean bottom cable survey or land exploration, the direct arrivals are usually contaminated with reflections and near surface noise (Skopintseva et al. 2016), and the source wavelet extraction is not reliable. We can also estimate the source signature during the model update process (Pratt 1999), and this approach is widely used in the FWI applications (Kamei et al. 2013(Kamei et al. , 2014. However, a simultaneous inversion of the velocity and the source wavelet may diverge when the initial velocity is poor and low frequency components in the data are missing. To mitigate the the source estimation error influence, deconvolution-based source-independent objective functions were developed (Lee & Kim 2003;Zhou & Greenhalgh 2003;Choi et al. 2005;Xu et al. 2006). For the crosshole seismic data, they mitigate the source wavelet effect by normalizing the data amplitude with a reference wavefield. However, the reference trace needs to be chosen wisely. If we implement FWI in the time-domain, it requires an additional Fourier transform. A convolution-based approach also requires to choose a reference trace from the observed and predicted data (Choi & Alkhalifah 2011). In their implementation, they convolve the observed data with the selected predicted data reference trace, and vice versa, they convolve the predicted data with the selected observed reference trace to formulate the objective function. In this source-independent objective function, the source function of both observed and predicted data are equally distributed in both terms, and thus, the effect of the source wavelet is mitigated. This approach is also successfully applied in the passive seismic inversion, in which the source function is unknown (Wang & Alkhalifah 2018). As the source origin time is unknown in passive seismic scenario, this source-independent objective function can be used to update the velocity model. Besides the deconvolution-and convolution-based methods, a new source-independent FWI using an amplitude-semblance objective function was proposed. In this approach, no reference trace selection is required (Chi et al. 2018). All of these approaches were focused on optimization problems with respect to the velocity model, often referred to as FWI.
Besides the uncertainty of the source wavelet, there is always a greater challenge from cycle-skipping issue in FWI, a direct result of the high non-linearity of the conventional l2 norm objective function with respect to the velocity model. Many methods developed recently aim to solve the cycle-skipping problem and retrieve the low-wavenumber model update (Shin & Cha 2008;Ma & Hale 2013;Wu et al. 2014;Biondi & Almomin 2014;Warner & Guasch 2016;Alkhalifah & Wu 2016). Among these methods, wavefield reconstruction inversion (WRI) was proposed to reduce the cycle-skipping issue by relaxing the wave equation constraint transforming it to a regularization term . To mitigate the source wavelet error influence in WRI, a source estimation method is specially designed for WRI (Fang et al. 2018). In WRI, the wavefield is still highly dependent on the velocity, and this dependency can be mitigated by replacing the velocity with a modified source function, which includes the veloctiy perturbation, we refer to as efficient wavefield inversion (EWI) . EWI is able to calculate the medium parameter perturbations directly, and it also proves its effectiveness in complex media (Song et al. 2019a;Song & Alkhalifah 2020a).
In this paper, we propose to combine the convolution-based source-independent approach and EWI to remove the source function effect in inverting the wavefield. This source-independent EWI (SIEWI) formulation replaces the data fitting term by a convolution-based source-independent misfit function represented by multiplications in the frequency domain. This new formulation not only maintains the advantage of efficiency and cycle-skipping immune feature in EWI, but also mitigates the source wavelet effect. We test the proposed approach on the Overthrust and a modified Marmousi model, and the results demonstrate the accuracy of the proposed method. An application on a 2D real dataset also demonstrates the effectiveness of SIEWI.

Full-waveform inversion (FWI)
Conventional FWI measures the misfit between observed and predicted data using the wave equation as a constraint on the wavefield, as follows (Tarantola 1984): where i and j denote the source and receiver indexes in which we sum over. C acts as the mapping operator projecting the background wavefield ui vector to the receiver positions. L(V ) = (V0 + V )ω 2 + ∇ 2 = L0 + ω 2 V is the modelling operator, and L0 = V0ω 2 + ∇ 2 corresponds to the background modelling operator. V denotes the squared slowness perturbation, and we set it to zero at the beginning of the inversion process. According to the adjoint-state method, the squared slowness gradient is evaluated by (R.E. Plessix 2006): where ∆di,j represents the data residual, and (L −1 ) T (C T ∆d * i,j ) represents the back-propagated wavefield vector. T denotes ordinary matrix transpose, the superscript * represents the complex conjugate, and means taking the real part of a complex value. Equation 2 shows that the gradient of FWI is obtained by calculating the product of the back-propagated wavefield and background wavefield (Choi et al. 2008). FWI requires a large amount of computational cost and disk memory in storing the background wavefields. The backgournd modelling operator used to calculate the background wavefield is dependent on the velocity. This dependency leads to the high non-linearity of the l2 norm objective function.

Wavefield reconstruction inversion (WRI)
WRI mitigates these issues by transforming the wave equation to an additional regularization term , as follows: By controlling the positive weighting factor 2 , the accuracy of the wave equation is relaxed to enforce more data fitting even with an inaccurate velocity model. However, in this formulation, the wavefield still has the dependency on the modelling operator L(V ). To mitigate this dependency, we linearise the objective function by replacing the true source function with a modified one(Alkhalifah & Song 2019).

Efficient wavefield inversion (EWI)
EWI presents an optimization for the wavefield independent of the velocity perturbation in the wave equation operator, given by : where, fei represents the modified source function. In equation 4, there are two independent variables: the wavefield ui and the modified source function fei. We calculate the wavefield using a linear equation given by an augmented wave equation: The least-squares solution of equation 5 is to solve (C T C + 2 L * 0 L0)ui = C T di,j + 2 L * 0 fei. L * 0 means the conjugate transpose form of L0. In each frequency, the modelling operator L0 in equation 5 remains stationary. If the wavefield u is accurately reconstructed, the modified source fe can be derived using: As the objective function E is a quadratic function of fei in equation 4, by setting ∂E ∂f ei = 0, fei can be evaluated using: Initially, the wavefield only has single scattering components from equation 5. As medium parameter perturbations are stored in fe, inner iterations between equations 5 and 7 are used to add the multiscattering components to the wavefield corresponding to secondary sources. This is discussed in details in . The squared slowness perturbation is not treated as an unknown parameter in EWI. After optimizing the two independent variables u and fe using inner iterations, V can be directly calculated using a direct division process according to equation 6: The squared slowness perturbation V can be calculated by each individual source. However, a single source has a limited illumination to the subsurface, so we calculate V by stacking all the sources, given by: where λ is a small positive real number to avoid dividing by zero. Usually, we use one percentage of the maximum value of the wavefield amplitude u * i ui. Multiscattering information provides more illumination to the model, so V benefits from more low wavenumber updates (Alkhalifah & Wu 2016;Song et al. 2019b). However, the calculated squared slowness perturbation is only accurate when the wavefield reconstruction result from equation 5 is reliable.

Source-independent efficient wavefield inversion (SIEWI)
The accuracy of the wavefield reconstruction relies on the initial velocity and the source function. If a wrong source wavelet is used in equation 5, the reconstructed wavefield will have a significant error compared to the true wavefield. In this case, the error in u will propagate to fe and V through equations 5 to 8, and hamper the convergence of the inversion process. This is a great challenge in real data applications, because the true wavelet is difficult to retrieve. To solve this problem, we propose to use a source-independent EWI (SIEWI). In time-domain, a source-independent based objective function is stated as (Choi & Alkhalifah 2011): Downloaded from https://academic.oup.com/gji/advance-article-abstract/doi/10.1093/gji/ggaa196/5829858 by King Abdullah University of Science and Technology user on 13 May 2020 Algorithm 1 Source-independent efficient wavefield inversion Input: Observed data do; Initial velocity model V0 and the source function f ; Number of outer iteration number N iter and inner iteration number N ; Selected frequency band; weighting factor 2 ; Output: Inverted velocity model.
for ifre=1:N f re; (Selected frequencies loop) do Initialize the source function fei = fi, the initial operator L0(V0); for i=1:N ;(Inner iteration loop) do Solve the source-independent augmented wave equation to obtain the wavefield ui; Invert for the modified source function fei = L0ui; end for end for di,j and di,j represent the predicted and the observed data, and d i,k and d i,k are defined as the reference trace at the k th receiver position from the predicted and observed data, respectively. The symbol represents the convolution operation in time. For each shot gather, the predicted data are convolved with a selected observed data reference trace, and vice versa. By doing so, the source function is represented equally in both terms, and its effect is removed from the objective function (Choi & Alkhalifah 2011). However, with this objective function we may still encounter the cycle-skipping problem caused by a poor initial velocity. We formulate an SIEWI objective function by replacing the l2 norm data fitting term with a sourceindependent objective function. This new formulation aims to avoid the source estimation and mitigate the cycle-skipping problem, simultaneously. It can be formulated in the frequency domain using the same reference trace concept as: In equation 11, Cui is equivalent to di,j, which denotes the predicted data. C k is a mapping operator which just takes the value from the wavefield ui at the reference trace location, k. So we have d i,k = C k ui. In frequency domain, the convolution process becomes an element by element multiplication, and we use the symbol to denote this. The wavefield u and modified source fe are still two independent variables in SIEWI. As a result, we can calculate the wavefield by solving a source-independent augmented wave equation: We solve (P * k P k + 2 L * 0 L0)ui = 2 L * 0 fei to get the least-squares solution of equation 11. Here, P k = C d i,k − di,j C k This source-independent augmented wave equation is able to reconstruct the wavefield with limited source function influence. Then we can calculate the modified source function and squared slowness perturbation using equations 7 and 9. In equation 9, we calculate the squared slowness perturbation V pixel by pixel, while the source function fi is zero everywhere except at the source location in the space domain. As a result, the accuracy of fi barely affects the calculation of V as long as the wavefield ui can be accurately reconstructed by SIEWI. Finally, the background model is updated:

IMPLEMENTATION
In this section, we explain the concept and functionality of the inner and outer iterations used on EWI and elaborate on the algorithm within the context of the source independent implementation.

The inner iterations
The objective function of either EWI or the SIEWI (in equations 4 and 10) has two terms. The first one aims to extract the wavefield information from the data. The second one is devoted to making sure this wavefield fits the wave equation corresponding to the current model. To extract an accurate wavefield that fits the data, and yet satisfies the wave equation for an inaccurate background velocity, the EWI objective function uses a modified source function fe. Specifically, fe holds the the velocity perturbations information necessary to produce a wavefield that fits the data for an inaccurate background velocity   model. In the first inner iteration, we start with fe = f (the true source). In this case, the calculated wavefield u is a compromise between fitting the wave equation and data by controlling 2 . For a smooth background model, the resulting wavefield and modified source function reflect energy corresponding to single scattering of what is required to fit the data. If only a single inner iteration is used, the original EWI reduces to WRI in the wavefield inversion part of the implementation, with respect to inversion for the source. Within each selected frequency, any update in the velocity requires a new matrix inversion. EWI introduces an inner iteration to update u and fe, and this operation is equivalent to adding scattering components to the wavefield. As fe is updated with the perturbations of the model, these perturbations act as secondary sources in fe weighted by 2 . The larger the 2 , the more secondary scattering we add to the wavefield. However, the larger the 2 , the smaller the  perturbations we extract from fitting the data. How much of the scattering we extract from the data to fe and how much of these secondary sources are used to develop the wavefield need to be considerably balanced. The wavefield initially will contain only single scattering components between the receivers and secondary sources in fe. Another iteration of the linear inversion for u, with the updated fe will add second-order scattering energy to the wavefield (Alkhalifah 2019). In practice, two or three inner iterations is enough, as scattering energy higher than the second-order is reasonably weak to effectively contribute the velocity update. Thus, in all the examples shown in this paper, we use two inner iterations.

The outer iterations
Although inner iterations in EWI can include multi-scattered energy in the reconstruted wavefield u, EWI is not completely immune to the cycle-skipping problem when the velocity perturbation is very large. In this case, we use outer iterations to update the velocity model throughout the whole selected frequency band to gradually improve the inversion result (Song & Figure 6. The inverted velocity using EWI after one sweep over the all frequencies with one iteration per frequency with the true wavelet.  Alkhalifah 2020b). The utilization of outer iterations in frequency domain waveform inversion implementations is sometimes needed to further cement the data fitting for low frequencies after including updates from the high frequencies (Pratt 1999). The outer iterations involve repeating EWI from low to high frequencies once again starting from the velocity model obtained in the previous EWI sweep. Since, in our EWI implementation, we update the velocity once per frequency, the outer iterations are expected to add additional information as the velocity model matures. The outer iteration also can be used to include more data in the inversion as we can include frequencies in the second sweep not used in the first sweep. The inner and outer loops are demonstrated in Algorithm 1.

EXAMPLES
In this section, we show results from synthetic data generated from the Overthrust model and a modified Marmousi model. Finally, we demonstrate the effectiveness of the proposed approach on field data from Western Australia.

Overthrust model
We apply the proposed method on the Overthrust model. The true velocity is shown in Fig. 1a, and the initial velocity we use is a highly smoothed version of the true velocity shown in Fig. 1b. The true time-domain wavelet is shown in Fig. 2a. The peak frequency of the Ricker wavelet is 10 Hz, and its frequency spectrum is shown in Fig. 2b. We transform the time-domain wavelet using fast Fourier transform (FFT) to generate synthetic frequency-domain data. In practical FWI implementations, the true wavelet information is usually unavailable, so we might end up using a wrong wavelet such as the one in Fig. 3a to perform the wavefield reconstruction and inversion. The frequency spectrum of the wrong wavelet is shown in Fig. 3b, and we can easily notice the phase distortion and amplitude change compared to the frequency spectrum of the original wavelet.
The key point of wavefield reconstruction-based methods is to invert for a good wavefield at the beginning. Figs. 4a and 4b Figure 10. The inverted velocity using SIEWI after one sweep over the all frequencies and one iteration per frequency with the wrong wavelet and reference traces at a fixed receiver location. are true wavefields calculated using the wave equation with the true and initial velocities, respectively, at 4 Hz. If we use the augmented wave equation 5 to invert the wavefield using the wrong wavelet and the initial velocity, the inverted wavefield is shown in Fig. 4c. As the black arrow points out, we see that there is obvious polarity change near the source location. By comparison, the wavefield calculated by the source-independent augmented wave equation using equation 11 is shown in Fig. 4d, and it is obvious that no polarity change exists near the source location.
Downloaded from https://academic.oup.com/gji/advance-article-abstract/doi/10.1093/gji/ggaa196/5829858 by King Abdullah University of Science and Technology user on 13 May 2020        We use 21 sources uniformly distributed on the surface, and all the grid points on the surface act as receivers. Receivers' positions are fixed like ocean botton cable (OBC) acquisition. We transform the true wavelet shown in Fig. 2a into the frequency domain, and generate the true data ranging from 3 Hz to 15 Hz with a sampling interval of 0.5 Hz. The 2 in this example is set 10 7 . A detailed study on the proper epsilon is given in . We first perform a standard WRI and FWI using the accurate wavelet. After one iteration through the whole frequency band, the inverted results are shown in Figs. 5a and 5b for WRI and FWI, respectively. Fig. 6 shows the inversion result using EWI. As the initial model is reasonably good, all the three methods reconstruct the main structures in the true model. If we use the wrong wavelet to repeat the inversion, the results using WRI and FWI are shown in Figs. 7a and 7b, respectively. We observe that both methods fail to get a reasonable model. Using the wrong wavelet, the inverted velocity using EWI after one sweep over the all frequencies with one iteration per frequency is shown in Fig. 8a. There are strong source-signature artifacts in the shallow part of the inverted model. Using SIEWI, the inversion result is shown in Fig. 8b. We can see that the proposed method recovers detailed structure of the true model even with one iteration over all frequencies. After five outer iterations, the inverted velocity using SIEWI is shown in Fig. 9, and all structures in the Overthrust model are reasonably inverted with high resolution.
In SIEWI, we need to choose reference traces from both observed and predicted data.
In the examples shown above, we choose the reference trace closest to the source location for both observed and predicted data, and this means we use reference traces at different receiver locations for different shot gathers. If we choose a reference trace at a fixed receiver location (receiver at the middle of the model) for all the shot gathers, the inversion result is shown in Fig. 10. Although there are no strong artifacts, no obvious update exists in the model.

Marmousi model
We further apply SIEWI on a modified Marmousi model, which is shown in Fig. 11a. The initial velocity we use is linearly increasing with depth, which is shown in Fig. 11b. There are 36 sources with 250 m spacing interval and 370 receivers fixed evenly distributed on the surface. We still use the wavelet in Fig. 2a to generate the observed data. The frequency band we use in this example is from 4 Hz to 10 Hz, and the sampling interval is still 0.5 Hz. The 2 we use here is also 10 7 . If we perform the WRI and FWI after one sweep over the all frequencies with one iteration per frequency, the inverted velocity models are shown in Figs. 12a and 12b, respectively. The structures of the true model are recovered more or less by both two methods, but they clearly suffer from the cycle skipping issue. Using the same inversion setup and strategy, EWI can recover more details of the true model due to utilizing the multiscattering components in the reconstructed wavefields , which is shown in Fig. 13. If we use the wrong wavelet (Fig. 3a) to repeat the inversion process using WRI and FWI, we end up with the inverted models shown in Fig. 14a and 14b, respectively. It is obvious that using the wrong wavelet, we cannot get reasonable inversion results. Fig. 15a shows the inversion result using EWI with the wrong wavelet. We observe that a wrong wavelet will be a more severe problem for wavefield reconstruction based methods, including WRI and EWI. For EWI, there are strong artifacts near the source location. In SIEWI implementation, we still first use both reference traces at the source location. By comparison, SIEWI is able to invert for a reasonable velocity model after one sweep over the all frequencies with one iteration per frequency, which is shown in Fig. 15b. After another four outer iterations, the final inverted velocity model is shown in Fig. 16. It is obvious that the proposed method can recover a high-resolution velocity model regardless of the wavelet accuracy. Besides the source-independent based objective functions, we can also perform a source estimation step within the inversion process. With a source estimation step in WRI and FWI (Fang & Herrmann 2015;Pratt 1999), the inversion results after five outer iterations are shown in Figs. 17a and 17b, respectively. Even with a source estimation step, FWI still fails to recover a reasonable velocity model. Though WRI benefits a lot from the source estimation, it is still not comparable to SIEWI. To demonstrate the importance of the reference trace selection, we show the inversion result with a reference trace at a fixed receiver location (receiver at the middle of the model) after one sweep over the all frequencies with one iteration per frequency in Fig. 18. There is barely any useful information getting recovered in the model compared to the result from a dynamic reference trace near the source location.
In order to avoid committing the inverse crime and demonstrate the validity of the proposed method in more realistic cases, we generate the observed data using a time-domain acoustic wave equation and elastic wave equation solver with an explosive source. We use a time-domain Ricker wavelet to generate the data. One shot gather from the acoustic data and the z-component elastic data at the distance of 2.75 km are shown in Figs. 19a,and 19b,respectively. We also show one shot gather from the acoustic data and the z-component elastic data at the same location corresponding to the initial velocity in Figs. 20a, and 20b, respectively. The absence of reflection waves indicates that the initial velocity model is smooth and generally crude. Clearly, we see that the z-component of the elastic data contains surface waves and S-wave artifacts. Thus, we use this z-component of the elastic data in our acoustic inversion. We build the S-wave velocity the same way as previous by dividing it by a scalar, and use a constant density. We use the same inversion setup as what we used before in this example. The initial source function used to reconstruct the wavefield in all frequencies is Delta function. Figs. 21a and 21b show the inversion results after one and five outer iterations (sweep from low to high frequency) using SIEWI. We observe that the main structures of the Marmousi model are reasonably recovered, especially in the shallow part. In this example, we find that wavefield reconstruction based methods (WRI and EWI) are more vulnerable to the wrong source function, so we repeat the no inverse crime inversion using FWI with the same inversion setup. The inversion results after one and five outer iterations are shown in Figs. 22a and 22b, respectively. It is obvious that FWI fails to recover a reasonable velocity model, and strong artifacts are shown in the shallow part. Figure 32. Filtered observed data (left, right) at 4.6 km, and corresponding predicted data (middle) using the inverted velocity.

CGG Field data
The 2D real data we use is acquired by CGG using a Broadseis acquisition system from North-Western Australia Continental Shelf (Soubaras & Dowle 2010). Frequencies below 2.5 Hz has been filtered out because of the low signal-to-noise ratio. Fig. 23 shows one shot gather of the real data at the distance of 4.6 km. The original data set contains 1824 shot gathers with an interval of 0.01875 km, and we choose 100 shot gathers to cover the 2D survey area. The target survey area we focus is 12.5 km long and 3.75 km in depth, and the initial velocity model is borrowed from Kalita & Alkhalifah (2017, as shown in Fig. 24. We resample the data and use 324 receivers from each streamer with an interval of 25 m. This real data is collected using a marine acquisition system, and the largest offset is 8.1 km. Without the source wavelet estimation, we perform SIEWI from 3 Hz to 12 Hz with a dense sampling of 0.15 Hz. The 2 used here is set to 10 8 based on our tests. The reference trace location is at the smallest offset for every shot gather. The final inversion result is shown in Fig. 25. It is obvious that the low velocity layers between 2.0 km to 2.5 km depth are well inverted. The check shot velocity from a well at 10.5 km is shown in Fig. 26, and the velocity between 1 km and 2 km is obtained by interpolation. The SIEWI inverted velocity model fits the well-log velocity well below 2 km, which demonstrates the effectiveness of SIEWI on the real data. We perform reverse time migration (RTM) using the initial and final velocities, and the images are shown in Figs. 27 and 28, respectively. As the red arrows point out, we observe multiple improvements in the RTM image using the inverted velocity. We also show the angle domain common-image gathers (ADCIGs) at locations x = {3, 4, 5, 6, 7, 8, 9} km ranging between 0 and 40 • corresponding to the initial and inverted velocity models in Figs. 29 and 30, respectively. The ADCIGs from the inverted velocity are flattened in some areas (as the red arrows point out) compared to those from the initial velocity, and these features indicate the improvement of the inverted velocity. Finally, we compare one filtered shot gather with the predicted data generated from the initial and inverted velocity models. We invert for the source wavelet in the time domain by back propagating the observed data to the source location using the inverted velocity model. From Fig. 31, we see that the predicted data using the initial velocity lacks most of the reflections apparent in observed data, and their diving waves at far offset don't match. Using the inverted model, the predicted data recovers the reflection events seen in the observed data, as shown in Fig. 32. As the kinematic features of the inverted velocity improves, the diving waves start to fit better for the observed and predicted data.

DISCUSSION
In the wavefield reconstruction based methods including WRI and EWI, 2 balances the weight between the data and the wave equation fitting. If the initial velocity is poor, a small 2 will reduce the cycle-skipping issue. Numerical experiments show that a single fixed value of 2 is sufficient (Leeuwen & Herrmann 2015). We use trial and error to find an appropriate 2 in different scenarios. One quick way to find a proper 2 is to use a scaling the highest eigenvalue of the auxiliary linear system Leeuwen & Herrmann (2015). This weighting factor 2 can also be automatically determined using a multiplicative cost function da Silva & Yao (2017). In EWI, for the inner iterations for updating u and fe, we start with fe = f . In the first inner iteration, the resulting wavefield contains single scattering information to fit the data. As fe is updated with parameter perturbations in the model, these perturbations act as secondary sources in the following inner iterations. As a result, EWI can reconstruct a wavefield including multiscattering information by using cheap inner iterations , only a few outer iterations will achieve convergence.
The source-independent objective function in equation 10 should theoretically work in combination with the original WRI. The source independence we refer to in this paper pertains to the source function (or wavelet). As EWI admits advantages in efficiency and accuracy over WRI, we directly use a source-independent EWI (SIEWI) objective function. The proposed method requires two reference traces from the observed and predicted data. From what we experienced in the examples, we find that it is better to choose the reference trace near the source location. This is because the reference trace at the near offsets contain most of the source function information, and this reference trace selection strategy is also suggested by Wang & Alkhalifah (2018). We cannot use a referenced trace at a fixed receiver position for all the shot gathers. Choosing a reference trace at a fixed receiver location will allow us to extract less useful information of the subsurface model, and slow down the convergence rate of the optimization.

CONCLUSIONS
We proposed a source-independent efficient wavefield inversion (SIEWI) to mitigate the source wavelet effect by incorporating source-independent objective into the original EWI. If we start with a wrong source wavelet in the original EWI, the source wavelet error will propagate into the reconstructed wavefield. In this case, the theory of EWI breaks down and the optimization problem diverges from the desired solution. The proposed method removes the source wavelet accuracy dependency while maintaining the EWI features. Applications to synthetic data generated from the Overthrust model and a modified Marmousi model show that SIEWI is capable of obtaining reasonable inversion results regardless of the source wavelet accuracy. The test on a 2D field data also demonstrates the effectiveness of the proposed approach.

ACKNOWLEDGEMENT
We thank KAUST for its support and the SWAG group for the collaborative environment. This work utilized the resources of the Supercomputing Laboratory at King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi Arabia, and we are grateful for that. We also thank CGG for providing the field data set and Geoscience Australia for providing the well-log information. We thank the assistant editor, the editor, Herve Chauris, and the reviewers for their critical and helpful review of the manuscript.