Seismic shot-encoding schemes for waveform inversion

A shot-encoding technique can be used in seismic waveform inversion to significantly reduce the computational cost by reducing the number of seismic simulations in the inversion procedure. Here we developed two alternative shot-encoding schemes to perform simultaneous-sources waveform inversion. The first scheme (I) encodes shot gathers with random-phase rotations applied to seismic traces. The second scheme (II) encodes shot gathers with random static time shifts. The well-known polarity encoding scheme (III) is just a special case of the random-phase rotation scheme. The second scheme is a variation of the conventional static shift encoding (IV), but the static time shifts in the second scheme are limited to one period of the dominant frequency. All encoded shot gathers are added up into a single super-shot gather for seismic waveform inversion. We perform the time-domain waveform inversion, using these shot-encoding schemes in conjunction with a restarted L-BFGS algorithm in the iterative inversion. The effectiveness and efficiency analyses demonstrate that the two shot-encoding schemes (I and II) proposed in this paper may improve the convergence of the iterative inversion, reduce the crosstalk effect among shots and consequently produce a subsurface velocity model with a high resolution.


Introduction
Seismic waveform tomography, also known as full waveform inversion (FWI), is an inversion technique for reconstructing wave velocities and other physical properties based on the field seismic waveforms (Tarantola 2004;Virieux & Operto 2009). This inversion technique is formulated as an optimisation problem, which minimises the data misfit between the simulated seismic waveforms and the field seismic waveforms. Most of the optimisation methods are iterative methods based on the gradient of the data misfit (Wang 2016). The most-used technique for gradient calculation is the adjoint-state method, in which the simulated forward wavefield is correlated with the backpropagation wavefield of the data misfit (Ren et al. 2013). The operation for the gradient calculation is performed for each shot gather independently, and thus the amount of computational time spent on the gradient calculation is linearly proportional to the number of shots in the seismic survey.
The gradient calculation for a single shot gather requires two wavefield simulations. One simulation is the forward wave propagation started from the source point, and another simulation is for the backward propagation of the data misfit. For a seismic survey with N S shots, it is necessary to perform wave simulations 2N S times to build the gradient vector. Therefore, the standard waveform inversion is a computationally expensive problem, mainly because it requires a lot of computational resources and time for the gradient calculation in a realistic sized application (Krebs et al. 2009;Wang & Rao 2009;Ben-Hadj-Ali et al. 2011;da Costa et al. 2018).
One strategy to accelerate the gradient calculation is to use a simultaneous-shot scheme in which all shots are added up to form a super shot and this super shot is used to perform forward and backward propagation (Krebs et al. 2009;Rao & Wang 2017). Using the simultaneous-shot scheme, the number of wavefield simulations required at each iteration of the inversion is reduced to just two. In reality, often only one super shot is enough to calculate the gradient (Castellanos et al. 2014) if the group of shots shares the same set of receivers such as OBC/OBN acquisitions. Nevertheless, there are strategies to deal with shots with different sets of receivers (Choi & Alkhalifah 2012;Son et al. 2014).
However, a drawback of the simultaneous-shot scheme for the gradient calculation is the crosstalk effect that arises from the wavefield correlation between the forward and backward propagations coming from different shots. Romero et al. (2000) and Krebs et al. (2009) suggest use a shot-encoding scheme to mitigate the crosstalk effect in migration and inversion, respectively. In the shot-encoded waveform inversion, it applies a mathematical transformation to the amplitudes of each shot before summing them up.
One of the shot-encoding schemes in waveform inversion is the random-polarity encoding (Krebs et al. 2009), which multiplies a shot gather by a random number either +1 or −1. Another shot-encoding scheme is applying a static shift to each shot chosen randomly (Schuster et al. 2011). In addition, strategies exist to apply the codes for each shot; for example, using a multiscale strategy (Wang 2016), creating several super shots instead of one (Castellanos et al. 2014) or combining encoding schemes with deblurring filters in order to mitigate the crosstalk effect in waveform inversion.
In this paper, we propose two alternative shot-encoding schemes. The first scheme (I) applies random-phase rotations to seismic traces and their Hilbert transforms. The second scheme (II) applies random static time shifts to seismic traces. The conventional random-polarity encoding scheme (III) is a special case of the first scheme. The conventional static shift scheme (IV) (Schuster et al. 2011) is a general case of the second scheme. But the second scheme (II) limits the static values to one period of the dominant frequency of the observed seismic data. We will demonstrate that in this way we can improve the result of seismic waveform inversion. We use the random-polarity encoding scheme (III) as a benchmark to compare to other schemes (I, II and IV).
In order to further mitigate the crosstalk effect through the iterations, we combine the shot-encoding strategies with a restarted L-BFGS algorithm proposed by Rao & Wang (2017). If comparing this restarted algorithm to the previous works such as Castellanos et al. (2014), Rao & Wang's restarted algorithm calculates the gradient difference in L-BFGS based on the secant equation, which represents the gradient differential rather than the gradient difference. Using the gradient differential will improve the accuracy and, in turn, the stability of convergence of the iterative inversion. The central idea of a restarted algorithm is to effectively increase the randomness in an accumulative way by re-setting the codes every few iterations.

Two shot-encoding schemes
We develop two shot-encoding schemes. One is the randomphase rotation scheme (I), and the other is the random static time-shifting limited to one period of the dominant frequency (II). We compare these two schemes to the randompolarity encoding scheme (III) and the conventional static shift scheme (IV).
The random-phase rotation scheme (I) is applied for the first time to seismic waveform inversion. To illustrate the effect of shot-encoding schemes, we show in figure 1a a super shot gather formed from a series of shot gathers, each of which only collects 20 near-offset (the source-receiver distance) traces. Figures 1b-1d show the super shot gather after random-polarity change, random-phase rotation and random static time shift. The waveform change of a shot gather is equivalent to the change of the associated source signature.
In the random-phase rotation scheme (II), all seismic traces within each shot gather are rotated with a random phase in the interval 0 ≤ ≤ (figure 1c). The phase rotation of a time series x(t) is performed using a time-domain equation: where H[x(t)] is the Hilbert transform of x(t) (Barnes 2007). This random-phase rotation scheme with the extreme values of cos = ±1 corresponds to the random-polarity encoding scheme (III). The value of phase is constant for all trace within a single shot gather but is randomly different for any other shot. After phase rotation, all traces at the same position but from different shot gathers are stacked so as to create a super-shot gather of the observed seismic data. In the wavefield simulations that are needed in the waveform inversion, the source signatures are phase rotated by the chosen phase angle and then deployed at the shot positions. In the static time-shift scheme (II), a randomly assigned time shift is applied to all traces of an individual shot gather (figure 1d). The term static means that the time shift is applied to an entire seismic trace. The value of static time is a constant for all traces within a shot gather but is randomly different for each shot. In the wavefield simulation during the waveform inversion, the same static time shift is applied to the source signature at the corresponding shot position. This encoding strategy is used on the simultaneous-source acquisition (Stefani et al. 2007).
In the conventional static shift scheme (IV), the static Δs is a function the record length of data (Schuster et al. 2011). In this paper, the second scheme (II) limits the static Δs to one period of the dominant frequency of observed data. That is where f 0 is the dominant frequency of the observed seismic data. We have found from our experiments that this constraint to the static values leads to a better result of seismic waveform inversion. We note that both the phase-rotation encoding scheme (I) and the static time-shift encoding (II) are time-domain versions of the frequency-domain phase encoding scheme (Ben-Hadj-Ali et al. 2011;Castellanos et al. 2014), where the amplitudes of seismic traces in the frequency domain are multiplied by a pure phase term e i with a randomly chosen phase . To show this assertion, we noticed that a time-shifted function satisfies the following Fourier transform property (Bracewell 1999) where F[x(t)]( ) is the Fourier transform of x(t). In the same way the phase rotation of equation (1) also has a frequency-domain version as (Purves 2014) Therefore, both the random-phase rotation scheme (I) and the random static time-shift scheme (II) are a multiplication of seismic amplitudes by an imaginary exponential with random phase and, consequently, are the time-domain versions of the frequency-domain phase encoding scheme. Figure 2a displays unmigrated constant-offset profile, in which each trace is collected from a shot gather and is annotated by its mid-point (of a source-receiver pair) coordinate. Figures 2b-2d show the profiles after random-polarity change, random-phase rotation and random static time shift. Each trace in these constant-offset profiles represents the waveform changes of all seismic traces within a shot gather after any shot encoding.

Waveform inversion with the restarted L-BFGS algorithm
The shot-encoding schemes can help to mitigate the crosstalk effect in a simultaneous-source waveform inversion to some extent. But it is necessary to change the shot codification throughout the inversion iterations (Krebs et al. 2009). Rao & Wang (2017) proposed a restarted L-BFGS algorithm for an iterative updating in the inversion, which can increase the effective randomness of a shot-encoding scheme and, thus, eventually achieve a complete crosstalk mitigation. This algorithm is different from the strategy of Krebs et al. (2009) (2017), Rao et al. (2019). In this paper, we choose Rao-Wang's restarted L-BFGS algorithm, because this algorithm outperforms other similar schemes, as shown by Rao et al. (2019).
In the shot-encoded waveform inversion, the objective function is defined by the data misfit between simulated waveforms and observed waveforms as where P is the simulated super shot, ∑ N S i=1 p obs (i) is the super gather of encoded field seismic data and ‖ ⋅ ‖ 2 is the L 2 norm. The simulated super shot P is the solution of the following wave equation: where ⃗ x is the spatial coordinate, and e i ⊗ s i (t) is the shotencoding operation.
The gradient vector of the objective function (5) is calculated as where R is the back-propagated wavefield and is calculated by solving equation (8), The T − t in the argument on the source term means that the wavefield is time reversed. Taking linearity of the wave equation, the gradient vector of equation (7) can be expressed as The first term is the correlations between the forward and backward propagation of the same shot i, whereas the second term is the crosstalk effect between shots i and j. For any shotencoding strategy to work properly, the first term should be exactly the gradient vector in the conventional waveform inversion as where P i (⃗ x, t) and R i (⃗ x, t) are the forward and backward wavefields of shot i.
To show that the shot-encoding schemes satisfy this requirement, we can apply the Parseval relation to the first term of equation (9) and obtain the relation whereê i ( ) is the frequency-domain encoding operator, andê * i ( ) is the frequency-domain reverse-time encoding operator. All three shot-encoding schemes satisfy that e * i ( )ê i ( ) = 1. In the same way, the second term of equation (9) satisfies the requirementê * i ( )ê j ( ) = 0. The crosstalk term approaches zero because of the variation of shot-encoding throughout the iterations. We will prove in the next section that the shot-encoded waveform inversion can produce the same result as that obtained from the conventional waveform inversion, but with a remarkable high efficiency.

The effectiveness of shot-encoded waveform inversion
For demonstrating the effectiveness, seismic waveform inversion is performed in the time domain and full waveform data are used as the input to the inversion, rather than inverting frequency-banded data sequentially for field seismic data (Wang 2011).
The SEG/EAGE Overthrust velocity model (figure 3a) is used to demonstrate the effectiveness of waveform inversion with different shot-encoding schemes. The 2D model consists of 401 × 98 cells and the cell size is 50 m. Therefore, the model dimension was 20 × 4.85 km 2 . The first five rows have a constant-velocity value of 2.1 km s −1 . This constantvelocity layer is also set as a known layer in the initial model, so that the gradient is zero in the area next to the source and receiver positions.
We use a finite-difference method to solve the acoustic wave equation (equations (6) and (8)) and apply an absorbing boundary layer to all sides of the model. We use 30-cell CPML method for the boundary absorbing layer (Liu & Sen 2018;He et al. 2020). There are 201 shots disposed at the surface. The source signature is a Ricker wavelet with the peak frequency of 6 Hz (Wang 2015a(Wang , 2015b. One receiver for each column of the model is positioned at the surface. Each trace has 1500 samples in time with a sample interval of 4 ms. A synthetic data set generated from the SEG/EAGE Overthrust velocity model (figure 3a) is used as the observed seismic data for waveform inversion. We implement the inversion iteratively and set a smooth version of the true model to be the initial model (figure 3b). The restarted L-BFGS algorithm has two parameters: , the number of iterations of each segment before re-starting the calculation and m, the number 910 of iterations with re-encoding within a segment. For all inversions presented in this paper, we use = 5 and m = 2, as the same as those used in the original paper (Rao & Wang 2017). Figure 3c is the result of the conventional waveform inversion without shot encoding. The inversion procedure is performed on 201 shots sequentially and implemented for 150 iterations for each shot gather.
The result obtained from this conventional waveform inversion can be used as a benchmark. Figure 4a is the result of waveform inversion with the random-polarity scheme (III), and figures 4b-4d show the results of waveform inversion with the random-phase-rotation scheme (I), the random static time-shift scheme (II) and random static timeshift scheme with ten periods of dominant frequency (IV), respectively. The stopping criterion is set as where J k is the misfit value where the iteration stop, J 0 is the misfit for the initial velocity, 2 is a threshold value to stop the iteration and = 0.03 is used in the numerical experiments presented in this paper. As observed in figure 4, all encoding strategies show a similar quality to the conventional FWI. The inversion with the random-polarity encoding (III) takes 775 iterations to reach the stopping criterion and is the slowest among the encoding strategies. The inversions with either the proposed random-phase-rotation encoding (I) or the random shift encoding (III) converge faster than the inversion with the polarity encoding, and need less than a half of the iterations than the inversion with the polarity encoding. The inversion with scheme IV spends around two times the number of iterations in the inversion with scheme II. This fact supports the idea of limiting the range of the static to one period of dominant frequency. Figure 5a shows the data misfit versus the iteration number for all different schemes (I, II, III and IV). All shotencoded inversions have a slow decay of the data misfit than the conventional waveform inversion, if considering only the number of iterations. However, figure 5b indicates that the conventional waveform inversion without shot encoding needs number of wavefield simulations two magnitude orders more than any shot-encoded waveform inversions need. Figure 5 clearly indicates that the inversion using the random-polarity scheme (III) and the inversion using the conventional static shift scheme with ten periods of dominant period (IV) have similar convergence rates. However, the inversions with two schemes proposed in this paper (I and II) show even faster convergence rates than those the two schemes mentioned above.
In these inversions, the step length for updating is fixed to unity for each iteration. An optimum step-length calculation needs a proper line search, and the linear search is affected by the crosstalk effect present on the gradient. However, as explained in Nocedal & Wright (2006), the secant equation for the L-BFGS method scales well enough the gradient and thus the calculation of step length can be avoided in most cases.
The correlation coefficient calculated position by position is a quantitative comparison between a waveform inversion result and the true velocity model (Gao & Wang 2020). The localised correlation coefficient at position (i, j) is calculated by  there is no any remarkable difference between the model from the conventional waveform inversion without shot encoding (figure 6a) and the model from any shot-encoded waveform inversions (figure 6b), and (ii) there is no any remarkable difference among all three shot-encoded waveform inversions. Therefore, figure 6 displays only one of these correlation-coefficient images (figure 6b). In fact, figure 6 clearly demonstrates the effectiveness of seismic waveform inversion when the inversion em-ploys seismic reflection data. It is referred to as the resolving power of the waveform inversion (Rao et al. 2006). One of the difficulties in seismic reflection waveform inversion is the reconstruction of the deep part of the subsurface model. We adopt the depth-dependent weighting scheme for the model update, closely following the strategy proposed in Wang & Rao (2009), to improve the reconstruction of deep structures. Figure 6 indicates an excellent recovery of the deep part velocity model from the 912 seismic reflection waveform inversion with or without shot encoding.

Conclusions
We have proposed the random-phase-rotation scheme (I) and the random static time-shift scheme constrained to one period of dominant frequency (II), for seismic waveform inversion, and compared them with the inversions using the random-polarity scheme (III) and the conventional static shift scheme defined as a function the record length of data (IV). We have demonstrated that all shot-encoding schemes, in conjunction with the restarted L-BFGS algorithm, can mitigate the crosstalk effect in the waveform inversion. Waveform inversion with any of these schemes produces a similar result to that obtained from the conventional waveform inversion without shot encoding. But the waveform inversion with any shot-encoding scheme has a remarkable efficiency in computation.
Among those four shot-encoding schemes, waveform inversion with two shot-encoding schemes proposed in this paper (I and II) have a better convergence, if compared to the inversion with other two shot-encoding schemes (III and IV). However, the inversion experiments have indicated that waveform inversion method with any shot-encoding scheme can make an excellent recovery of the deep part velocity model with the same effectiveness as the conventional waveform inversion without shot encoding.