Seismic Inversion by Multi-dimensional Newtonian Machine Learning

An example of a three-layer autoencoder architecture is shown in Figure 2, which includes an encoder network, latent space, and decoder network. The pink box emphasizes the encoder network that compresses the input data to a lower-dimensional space by using several neural layers. The number of neurons in each layer gradually decreases. The green box represents the latent space which contains the lowest-dimensional representation of the input data. The compressed data is then passed to the decoder network (purple box) to get the decoded waveform that has the same dimension as the input data. The decoder network is often the mirror of the encoder network which is composed of several neural layers with an increasing number of neurons in each layer. The misfit function of the autoencoder $\epsilon=(d_{input}-d_{decode})^{2}$ computes the differences between the input $d_{input}$ and decoded data $d_{decode}$, and the training process stops when the misfit reduces to a certain threshold. For a well-trained autoencoder, it should be able to extract the effective low-dimensional representation of the high-dimensional seismic data.

The reconstruction ability of an autoencoder is affected by several factors, such as the number of neural layers, the number of neurons in each layer, and the LS dimension. Here we only discuss and analyze how the LS dimension affects the reconstruction ability of an autoencoder. We design two autoencoders with the same encoder and decoder network, where each consists of multi-perceptrons. We set one of the autoencoder’s LS dimensions to two (la=2) and another one to four (la=4). Two thousand seismic traces are used to train these two autoencoders, and a subset of the training data is shown in Figure 3a. Figure 3b shows the decoded waveform from the autoencoder with a two-dimensional LS. It shows that the first arrival energy has been well recovered, but some reflection events are missing. Because the two-dimensional LS is too small to preserve the reflection energies. These missing reflection events can be clearly seen in Figure 3c which shows the data difference between the input and decoded data. Figure 3d shows the decoded waveform from the autoencoder with a four-dimensional LS, where both the first-arrival and reflection energies have been well recovered. Its data residual given in Figure 3e shows that the differences of the reflection events are largely minimized compared to Figure 3c. Therefore we conclude that an autoencoder with a larger LS dimension can better preserve and reconstruct input data information than an autoencoder with a smaller LS.

However, the data reconstruction ability of an autoencoder increases slowly when the LS dimension reaches a certain level. Figure 4 shows a curve that describes the data reconstruction error versus the LS dimension changes for this example. It shows that the reconstruction error drops rapidly when the LS dimension goes from one to four. But when the LS dimension is larger than four, the reconstruction error barely changes. This phenomenon suggests that the four-dimensional LS is good enough to preserve and reconstruct the entire information of the input data, and this four-dimensional LS is considered as the optimal LS dimension.

A $nt\times 1$ seismic trace can be viewed as a data point in a $nt$ dimensional space, where $nt$ is the number of time samples of a seismic trace. Therefore the FWI misfit function measures the spatial distance between the observed and predicted data point in an $nt$ dimensional space. Similarly, the MNML method measures the spatial distance of the LS feature of the observed and predicted trace in a $n\times 1$ dimensional space, where $n$ indicates the dimensionality of the LS features which is much smaller than $nt$. Figure 5 illustrates a keystone idea underlying the MNML method. An autoencoder compresses the high-dimensional observed and predicted seismic trace to a two-dimensional LS. The spatial distance between the compressed observed and predict two-dimensional LS features can be used to quantify the accuracy of the inverted model. A large distance means that the inverted model is far away from the true model, but a small distance indicates a strong similarity between the true and inverted model. In general, for the MNML method, the LS dimension could be an arbitrary number that is larger than one and smaller than $nt$. However, we can find the optimal LS dimension $n$ according to the relationship between the LS dimension $n$ and data reconstruction error of the autoencoder, which is demonstrated in Figure 4. We identify the optimal dimension as the ’elbow’ point of the curve which shows a relatively small decrease in data reconstruction error with increasing dimension.

The misfit function of the MNML method can be written as

where $\Delta z$ represents the multi-dimensional LS feature difference of the observed and predicted data. The locations of the source and receiver are represented by $\mathbf{x}_{s}$ and $\mathbf{x}_{r}$, respectively, where $s$ and $r$ indicate the source and receiver index. Here, $k$ indicates the dimension index of the multi-dimensional LS feature. For simplicity, we assume the LS feature dimension equals two. Then equation 1 can be re-written as

The gradient $\gamma(\mathbf{x})$ can be found by taking the derivative of the misfit $\epsilon$ to the velocity $v(\mathbf{x})$ as

Equation 3 shows that the velocity gradient calculation is carried out at each LS dimension.

To compute the gradient $\gamma(\mathbf{x})$, we need to define the relationship between the velocity perturbation and the multi-dimensional LS feature perturbation mathematically. To do so, we build a multi-variable connective function that measures the similarity between the observed and predicted seismic data through cross-correlation:

where $p_{(z_{1},z_{2})}^{pred}(\mathbf{x}_{r},t;\mathbf{x}_{s})$ indicates a predicted seismic trace recorded at the receiver location $x_{r}$ for a source at $\mathbf{x}_{s}$. The subscripts $z_{1}$ and $z_{2}$ represent the LS feature values of this trace for the first and second LS dimensions, respectively. Similarly, $p_{(z_{1}-\tilde{z}_{1},z_{2}-\tilde{z}_{2})}^{obs}(\mathbf{x}_{r},t;\mathbf{x}_{s})$ represents the observed seismic trace at the same source-receiver location with the LS feature value equal to $(z_{1}-\tilde{z}_{1},z_{2}-\tilde{z}_{2})$, where $(\tilde{z}_{1},\tilde{z}_{2})$ is the difference between the observed and predicted multi-dimensional LS features. This difference will be zero ($(\tilde{z}_{1},\tilde{z}_{2})=(0,0)$) if the velocity model we build is accurate enough, otherwise $(\tilde{z}_{1},\tilde{z}_{2})\neq(0,0)$. We need to (1) find the optimal LS feature difference $(\tilde{z}_{1},\tilde{z}_{2})=(\Delta z_{1},\Delta z_{2})$ that maximize equation 4; then (2) recover the accurate velocity model by minimizing the difference. The optimal difference can be found by setting equation $\nabla F(\tilde{z}_{1},\tilde{z}_{2})|_{(\tilde{z}_{1}=\Delta z_{1},\tilde{z}_{2}=\Delta z_{2})}=0$ as

Rewriting equation 5 in matrix form and combining it with the multi-variable implicit function theorem (Krantz and Parks,, 2012; Guo,, 2016) , we can get the Fréchet derivatives

The detailed derivations from equations 5 and 6, and the formulas for the partial derivatives on the right-hand-side of equation 6 can be found in the Appendix B.

Substituting equation 6 into equation 3, we get the gradient formula for the two-dimensional NML

If we assume $z_{1}$ and $z_{2}$ are weakley correlated, then $\frac{\partial^{2}F}{\partial\tilde{z}_{1}\partial\tilde{z}_{2}}\approx 0$. Therefore equation 8 can be re-written as

where $E_{1}=\frac{1}{\frac{\partial^{2}F}{\partial\tilde{z}_{1}^{2}}}$ and $E_{2}=\frac{1}{\frac{\partial^{2}F}{\partial\tilde{z}_{2}^{2}}}$ are the weighting parameters. For a more general $n$-dimensional case, its gradient formula can be found by extending equation 8 to $n$-dimensions

The Fréchet derivative $\frac{\partial p}{\partial v}$ of the first-order acoustic equation can be written as (Chen and Schuster,, 2020)

where $\rho$ and $v$ represent the density and velocity, respectively. The particle velocity is indicated by $\mathbf{v}$, $g_{p}$ is the Green’s function, and $*$ means temporal convolution. Substituting equation 10 into equation 9 allows the MNML gradient to be expressed as

where $\Delta p_{\mathbf{z}}(\mathbf{x},t;\mathbf{x}_{s})=\sum_{k}\frac{\partial p_{(z_{1}-\Delta z_{1},z_{2}-\Delta z_{2})}^{obs}(\mathbf{x}_{r};\mathbf{x}_{s})}{\partial\tilde{z}_{k}}E_{k}\Delta z_{k}$ denotes the virtual source at the receiver location. For each LS dimension, a local virtual source is computed by weighting its LS feature difference $\Delta z_{k}$ on the partial derivative $\frac{\partial p_{(z_{1}-\Delta z_{1},z_{2}-\Delta z_{2})}^{obs}(\mathbf{x}_{r};\mathbf{x}_{s})}{\partial\tilde{z}_{k}}$, and then re-scale its value by dividing the weighting parameter $E_{k}$. The final virtual source $\Delta p_{\mathbf{z}}(\mathbf{x},t;\mathbf{x}_{s})$ can be obtained by summing all the local virtual sources at each LS dimensions together. Once we have the gradient, the velocity model can be updated by using the steepest descent formula

where $k$ represents the iteration index and $\alpha_{k}$ indicates the step length.

The workflow of the MNML inversion is shown in Figure 6, where an autoencoder is trained to generate multi-dimensional LS features of the seismic data. The seismic traces in the observed shot gathers are often used to train and validate the autoencoder network. This autoencoder is only trained once and its parameters are fixed during the MNML inversion. At each iteration of the MNML inversion, we feed the observed and predicted data to the well-trained autoencoder to get their LS features. We then compute the velocity gradient using equation 11 and update the current velocity model. We stop the inversion when the LS feature misfit goes below a certain threshold.

We first test the MNML method on the data generated by checkerboard models with three different acquisition geometries to demonstrate the capability of the MNML method in dealing with different acquisition systems.

A portion of the Pluto model is selected as the true velocity model (shown in Figure 11a) to test the MNML method. A crosswell acquisition system is used to generate the seismic data. The source well is deployed at x=10 m which contains 59 shots at an equal interval of 30 m. Each shot has 177 receivers evenly distributed in the receiver well, which is located at x = 1750 m. A 20-Hz Ricker wavelet is used as the seismic source, and a linearly increasing velocity model is used as the initial velocity model. The minimum and maximum velocity of the initial model are equal to 3100 m/s and 3700 m/s, respectively. The acoustic finite-difference modeling is used to generate the observed data based on the true model.

The field data test uses a crosswell dataset collected by Exxon in Texas (Chen et al.,, 1990). The source and receiver wells are 305 m deep and separated by 183 m from each other. There are 97 shots fired at a shot interval of 3 m from z=9 m to z=305 m using downhole explosive charges. Each shot has 96 receivers distributed in the receiver well (Chen and Schuster,, 2020). The seismic data are recorded for 0.375 s with a time interval of 0.25 ms. The data contains extreme noise and tube waves. We first bandpass the data to 80-400 Hz to remove the extreme linear noise. We then use a nine-point median filter to remove the tube wave. Finally, we separate the up and down-going waves using FK filtering. A similar processing workflow with more details can be found in Cai and Schuster, (1993); Dutta and Schuster, (2014); Chen and Schuster, (2020). We plot one example of the raw and processed data in Figures 14a and 14b, respectively.

We assume the initial model is a homogeneous model with a velocity equals to 1900 m/s. We first use the NML method to recover the background velocity model. We design an autoencoder with four encoder and decoder layers and a single-dimensional LS. The first to the fourth encoder layers have a dimension of $3750\times 1$, $1000\times 1$, $200\times 1$, $20\times 1$, and the decoder network is the mirror of the encoder network. The training dataset contains the envelope of the observed seismic traces. We use the same training strategy described above to train this autoencoder. Figures 16c and 16d show the first iteration gradient and NML tomogram, which mainly contains the low-wavenumber velocity information.

To reconstruct the high-wavenumber velocity details, we apply the MNML method on the NML background tomogram. We build a new autoencoder with a six-dimensional LS (la=6) but keep the encoder and decoder network the same as the previous autoencoder used in the NML. The training dataset for MNML includes the observed seismic traces rather than its envelopes. Figures 17a and 17b show the first iteration gradient and MNML tomogram (la=6), which contains the high-wavenumber velocity details. The MNML tomogram (la=6) shows a noticeable resolution increase compared to the NML tomogram. To validate the robustness of the MNML tomogram, we compare it with the FWI tomogram shown in Figure 17d. Same as the MNML inversion, FWI also uses the NML background tomogram as the starting model. The FWI tomogram is almost identical to the MNML tomogram (la=6), which means the MNML tomogram is correct.