{'fraction_non_alphanumeric': 0.03481864550864267, 'fraction_numerical': 0.01700198356474922, 'mean_word_length': 4.811650885137917, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Unsupervised Learning of Full-Waveform Inversion: Connecting CNN and Partial Differential Equation in a Loop

This paper investigates unsupervised learning of Full-Waveform Inversion (FWI), which has been widely used in geophysics to estimate subsurface velocity maps from seismic data. This problem is mathematically formulated by a second order partial differential equation (PDE), but is hard to solve. Moreover, acquiring velocity map is extremely expensive, making it impractical to scale up a supervised approach to train the mapping from seismic data to velocity maps with convolutional neural networks (CNN). We address these difficulties by integrating PDE and CNN in a loop, thus shifting the paradigm to unsupervised learning that only requires seismic data. In particular, we use finite difference to approximate the forward modeling of PDE as a differentiable operator (from velocity map to seismic data) and model its inversion by CNN (from seismic data to velocity map). Hence, we transform the supervised inversion task into an unsupervised seismic data reconstruction task. We also introduce a new large-scale dataset OpenFWI, to establish a more challenging benchmark for the community. Experiment results show that our model (using seismic data alone) yields comparable accuracy to the supervised counterpart (using both seismic data and velocity map). Furthermore, it outperforms the supervised model when involving more seismic data.


INTRODUCTION
Geophysical properties (such as velocity, impedance, and density) play an important role in various subsurface applications including subsurface energy exploration, carbon capture and sequestration, (a) (b) Figure 2: An example of (a) a velocity map and (b) seismic measurements (named shot gather in geophysics) and the 1D time-series signal recorded by a receiver.
estimating pathways of subsurface contaminant transport, and earthquake early warning systems to provide critical alerts. These properties can be obtained via seismic surveys, i.e., receiving reflected/refracted seismic waves generated by a controlled source. This paper focuses on reconstructing subsurface velocity maps from seismic measurements. Mathematically, the velocity map and seismic measurements are correlated through an acoustic-wave equation (a second-order partial differential equation) as follows: where p(r, t) denotes the pressure wavefield at spatial location r and time t, v(r) represents the velocity map, and s(r, t) is the source term. Full-Waveform Inversion (FWI) is a methodology that determines high-resolution velocity maps v(r) of subsurface via matching synthetic seismic waveforms to raw recorded seismic data p(r, t), wherer represents the locations of seismic receivers.
A velocity map describes the wave propagation speed in the subsurface region of interest. An example in 2D scenario is shown in Figure 2a. Particularly, the x-axis represents the horizontal offset of a region, and the y-axis stands for the depth. The regions with the same geologic information (velocity) are called a layer in velocity maps. In a sample of seismic measurements (termed a shot gather in geophysics) as depicted in Figure 2b, each grid in the x-axis represents a receiver, and the value in the y-axis is a 1D time-series signal recorded by each receiver.
Existing approaches solve FWI in two directions: physics-driven and data-driven. Physics-driven approaches rely on the forward modeling of Equation 1, which simulates seismic data from velocity map by finite difference. They optimize velocity map per seismic sample, by iteratively updating velocity map from an initial guess such that simulated seismic data (after forward modeling) is close to the input seismic measurements. However, these methods are slow and difficult to scale up as the iterative optimization is required per input sample. Data-driven approaches consider FWI problem as an image-to-image translation task and apply convolution neural networks (CNN) to learn the mapping from seismic data to velocity maps (Wu & Lin, 2019). The limitation of these methods is that they require paired seismic data and velocity maps to train the network. Such ground truth velocity maps are hardly accessible in real-world scenario because generating them is extremely time-consuming even for domain experts.
In this work, we leverage advantages of both directions (physics + data driven) and shift the paradigm to unsupervised learning of FWI by connecting forward modeling and CNN in a loop. Specifically, as shown in Figure 1, a CNN is trained to predict a velocity map from seismic data, which is followed by forward modeling to reconstruct seismic data. The loop is closed by applying reconstruction loss on seismic data to train the CNN. Due to the differentiable forward modeling, the whole loop can be trained end-to-end. Note that the CNN is trained in an unsupervised manner, as the ground truth of velocity map is not needed. We name our unsupervised approach as UPFWI (Unsupervised Physics-informed Full-Waveform Inversion).
Additionally, we find that perceptual loss (Johnson et al., 2016) is crucial to improve the overall quality of predicted velocity maps due to its superior capability in preserving the coherence of the reconstructed waveforms comparing with other losses like Mean Squared Error (MSE) and Mean Absolute Error (MAE).
To encourage fair comparison on a large dataset with more complicate geological structures, we introduce a new synthetic dataset named OpenFWI, which contains 60,000 labeled data (velocity map and seismic data pairs) and 48,000 unlabeled data (seismic data alone). 30,000 of those velocity maps contain curved layers that are more challenge for inversion. We also add geological faults with various shift distances and tilting angles to all velocity maps.
We evaluate our method on this dataset. Experimental results show that for velocity maps with flat layers, our UPFWI trained with 48,000 unlabeled data achieves 1146.09 in MSE, which is 26.77% smaller than that of the supervised baseline H-PGNN+ (Sun et al., 2021), and 0.9895 in Structured Similarity (SSIM), which is 0.0021 higher; for velocity maps with curved layers, our UPFWI achieves 3639.96 in MSE, which is 28.30% smaller, and 0.9756 in SSIM, which is 0.0057 higher.
Our contribution is summarized as follows: • We propose to solve FWI in an unsupervised manner by connecting CNN and forward modeling in a loop, enabling end-to-end learning from seismic data alone. • We find that perceptual loss is helpful to boost the performance comparable to the supervised counterpart. • We introduce a large-scale dataset as benchmark to encourage further research on FWI.

PRELIMINARIES OF FULL-WAVEFORM INVERSION (FWI)
The goal of FWI in geophysics is to invert for a velocity map v ∈ R W ×H from seismic measurements p ∈ R S×T ×R , where W and H denote the horizontal and vertical dimensions of the velocity map, S is the number of sources used to generate waves during data acquisition process, T denotes the number of samples in the wavefields recorded by each receiver, and R represents the total number of receivers.
In conventional physics-driven methods, forward modeling is commonly referred to the process of simulating seismic datap from given estimated velocity mapsv. For simplicity, the forward acoustic-wave operator f can be expressed as Given this forward operator f , the physics-driven FWI can be posed as a minimization problem (Virieux & Operto, 2009) where ||p − f (v)|| 2 2 is the the 2 distance between true seismic measurements p and the corresponding simulated data f (v), λ is a regularization parameter and R(v) is the regularization term which is often the 2 or 1 norm ofv. This requires optimization per sample, which is slow as the optimization involves multiple iterations from an initial guess.
Data-driven methods leverage CNNs to directly learn the inverse mapping as (Adler et al., 2021) where g θ (·) is the approximated inverse operator of f (·) parameterized by θ. In practice, g θ is usually implemented as a CNN (Adler et al., 2021;Wu & Lin, 2019). This requires paired seismic data and velocity maps for supervised learning. However, the acquisition of large volume of velocity maps in field applications can be extremely challenging and computationally prohibitive.

METHOD
In this section, we present our Unsupervised Physics-informed solution (named UPFWI), which connects CNN and forward modeling in a loop. It addresses limitations of both physics-driven and data-driven approaches, as it requires neither optimization at inference (per sample), nor velocity maps as supervision.
3.1 UPFWI: CONNECTING CNN AND FORWARD MODELING As depicted in Figure 1, our UPFWI connects a CNN g θ and a differentiable forward operator f to form a loop. In particular, the CNN takes seismic measurements p as input and generates the corresponding velocity mapv. We then apply forward acoustic-wave operator f (see Equation 2) on the estimated velocity mapv to reconstruct the seismic datap. Typically, the forward modeling employs finite difference (FD) to discretize the wave equation (Equation 1). The details of forward modeling will be discussed Section 3.3. The loop is closed by the reconstruction loss between input seismic data p and reconstructed seismic datap = f (g θ (p)). Notice that the ground truth of the velocity map v is not involved, and the training process is unsupervised. Since the forward operator is differentiable, the reconstruction loss can be backpropagated (via gradient descent) to update the parameters θ in the CNN.

CNN NETWORK ARCHITECTURE
We use an encoder-decoder structured CNN (similar to Wu & Lin (2019) and Zhang & Lin (2020)) to model the mapping from seismic data p ∈ R S×T ×R to velocity map v ∈ R W ×H . The encoder compresses the seismic input and then transforms the latent vector to build the velocity estimation through a decoder. See the implementation details in Appendix A.1.

DIFFERENTIABLE FORWARD MODELING
We apply the standard finite difference (FD) in the space domain and time domain to discretize the original wave equation. Specifically, the second-order central finite difference in time domain ( ∂ 2 p(r,t) ∂t 2 in Equation 1) is approximated as follows: where p t r denotes the pressure wavefields at timestep t, and p t+1 r and p t−1 r are the wavefields at t + ∆t and t − ∆t, respectively. The Laplacian of p(r, t) can be estimated in the similar way on the space domain (see Appendix A.2). Therefore, the wave equation can then be written as where ∇ 2 here denotes the discrete Laplace operator.
The initial wavefield at the initial timestep is set to zero (i.e. p 0 r = 0). Thus, the gradient of loss L with respect to estimated velocity at spatial location r can be computed using the chain rule as where T indicates the total number of timesteps.

LOSS FUNCTION
The reconstruction loss of our UPFWI includes a pixel-wise loss and a perceptual loss as follows: where p andp are input and reconstructed seismic data, respectively. The pixel-wise loss L pixel combines 1 and 2 distance as: where λ 1 and λ 2 are two hyper-parameters to control the relative importance. For the perceptual loss L perceptual , we extract features from conv5 in a VGG-16 network (Simonyan & Zisserman, 2015) pretrained on ImageNet (Krizhevsky et al., 2012) and combine the 1 and 2 distance as: where φ(·) represents the output of conv5 in the VGG-16 network, and λ 3 and λ 4 are two hyperparameters. Compared to the pixel-wise loss, the perceptual loss is better to capture the region-wise structure, which reflects the waveform coherence. This is crucial to boost the overall accuracy of velocity maps (e.g. the quantitative velocity values and the structural information).

OPENFWI DATASET
We introduce a new large-scale geophysics FWI dataset OpenFWI, which consists of 108K seismic data for two types of velocity maps: one with flat layers (named FlatFault) and the other one with curved layers (named CurvedFault). Each type has 54K seismic data, including 30K with paired velocity maps (labeled) and 24K unlabeled. The 30K labeled pairs are splitted as 24K/3K/3K for training, validation and testing respectively. Samples are shown in Appendix A.3.
The shape of curves in our dataset follows a sine function. Velocity maps in CurvedFault are designed to validate the effectiveness of FWI methods on curved topography. Compared to the maps with flat layers, curved velocity maps yield much more irregular geological structures, making inversion more challenging. Both FlatFault and CurvedFault contain 30,000 samples with 2 to 4 layers and their corresponding seismic data. Each velocity map has dimensions of 70×70, and the grid size is 15 meter in both directions. The layer thickness ranges from 15 grids to 35 grids, and the velocity in each layer is randomly sampled from a uniform distribution between 3,000 meter/second and 6,000 meter/second. The velocity is designed to increase with depth to be more physically realistic. We also add geological faults to every velocity map. The faults shift from 10 grids to 20 grids, and the tilting angle ranges from -123 • to 123 • .
To synthesize seismic data, five sources are evenly placed on surface with a 255-meter spacing, and seismic traces are recorded by 70 receivers at each grid with an interval of 15 meter. The source is a Ricker wavelet with a central frequency of 25 Hz (Wang, 2015). Each receiver records time-series data for 1 second, and we use a 1 millisecond sample rate to generate 1,000 timesteps. Therefore, the dimensions of seismic data become 5×1000×70. Compared to existing datasets (Yang & Ma, 2019;Moseley et al., 2020), OpenFWI is significantly larger. It includes more complicated and physically realistic velocity maps. We hope it establishes a more challenging benchmark for the community.

EXPERIMENTS
In this section, we present experimental results of our proposed UPFWI evaluated on the OpenFWI.

IMPLEMENTATION DETAILS
Training Details: The input seismic data are normalized to range [−1, 1]. We employ AdamW (Loshchilov & Hutter, 2018) optimizer with momentum parameters β 1 = 0.9, β 2 = 0.999 and a weight decay of 1 × 10 −4 to update all parameters of the network. The initial learning rate is set to 3.2 × 10 −4 , and we reduce the learning rate by a factor of 10 when validation loss reaches a plateau. The minimum learning rate is set to 3.2 × 10 −6 . The size of a mini-batch is set to 128. All trade-off hyper-parameters λ in our loss function are set to 1. We implement our models in Pytorch and train them on 8 NVIDIA Tesla V100 GPUs. All models are randomly initialized.
Evaluation Metrics: We consider three metrics for evaluating the velocity maps inverted by our method: MAE, MSE and SSIM. Both MAE and MSE have been employed in existing methods (Wu & Lin, 2019;Zhang & Lin, 2020) to measure pixel-wise errors. Considering the layered-structured velocity maps contain highly structured information, degradation or distortion can be easily perceived by a human. To better align with human vision, we employ SSIM to measure perceptual similarity. Note that for MAE and MSE calculation, we denormalize velocity maps to their original scale while we keep them in normalized scale [-1, 1] for SSIM according to the algorithm.
Comparison: We compare our method with three state-of-the-art algorithms: two pure data-driven methods, i.e., InversionNet (Wu & Lin, 2019) and VelocityGAN (Zhang & Lin, 2020), and a physicsinformed method H-PGNN (Sun et al., 2021). We follow the implementation described in these papers and search for the best hyper-parameters for OpenFWI dataset. Note that we improve H-PGNN by replacing the network architecture with the CNN in our UPFWI and adding perceptual loss, resulting in a significant boosted performance. We refer our implementation as H-PGNN+, which is a strong supervised baseline. Our method has two variants (UPFWI-24K and UPFWI-48K), using 24K and 48K unlabeled seismic data respectively.

MAIN RESULTS
Results on FlatFault: Table 1 shows the results of different methods on FlatFault. Compared to InversionNet and VelocityGAN, our UPFWI-24K performs better in MSE and SSIM, but is slightly worse in MAE score. Compared to H-PGNN+, there is a gap between our UPFWI-24K and H-PGNN+ when trained with the same amount of data. However, after we double the size of unlabeled data (from 24K to 48K), a significant improvement is observed in our UPFWI-48K for all three metrics, and it outperforms all three supervised baselines in MSE and SSIM. This demonstrates the potential of our UPFWI for achieving higher performance with more unlabeled data involved.  Figure 4. Consistent with our quantitative analysis, more accurate details are observed in the velocity maps generated by UPFWI-48K. For instance, we find in the visualization results that both InversionNet and Ve-locityGAN generate blurry results in deep region, while H-PGNN+, UPFWI-24K and UPFWI-48K yield much clearer boundaries. We attribute this finding as the impact of seismic loss. We further observe that the slope of the fault in deep region is different from that in the shallow region, yet only UPFWI-48K replicates this result as highlighted by the green square.
Results on CurvedFault: Table 1 shows the results of CurvedFault. Performance degradation is observed for all models, due to the more complicated geological structures in CurvedFault. Although our UPFWI-24K underperforms the three supervised baselines, our UPFWI-48K significantly boosts the performance, outperforming all supervised methods in terms of all three metrics. This demonstrates the power of unsupervised learning in our UPFWI that greatly benefits from more unlabeled data when dealing with more complicated curved structure.
The bottom row of Figure 4 shows the visualized velocity maps in CurvedFault obtained using different methods. Similar to the observation in FlatFault, our UPFWI-48K yields more accurate details compared to the results of supervised methods. For instance, only our UPFWI-24K and UPFWI-48K precisely reconstruct the fault beneath the curve around the top-left corner as highlighted by the yellow square. Although some artifacts are observed in the results of UPFWI-24K around the layer boundary in deep region, they are eliminated in the results of UPFWI-48K. More visualization results are shown in Appendix A.3.  Loss Terms: We study the contribution of each loss term in our loss function: (a) pixel-wise 2 distance (MSE), (b) pixel-wise 1 distance (MAE), and (c) perceptual loss. All experiments are conducted on FlatFault using 24,000 unlabeled data. Figure 5a shows the predicted velocity maps for using three loss combinations (pixel-2 , pixel-1 2 , pixel-1 2 +perceptual) in UPFWI. The ground truth seismic data and velocity map are shown in the left column. For each loss option, we show the difference between the reconstructed and the input seismic data (on the top) and predicted velocity (on the bottom). When using pixel-wise loss in l 2 distance alone, there are some obvious artifacts in both seismic data (around 600 millisecond) and velocity map. These artifacts are mitigated by introducing additional pixel-wise loss in l 1 distance.   With perceptual loss added, more details are correctly retained (e.g. seismic data from 400 millisecond to 600 millisecond, velocity boundary between layers). Figure 5b compares the reconstructed seismic data (in terms of residual to the ground truth) at a slice of 525 meter offset (orange dash line in Figure 5a). Clearly, the combination of pixel-wise and perceptual loss has the smallest residual.

ABLATION STUDY
The quantitative results are shown in Table 2. They are consistent with our observation in qualitative analysis (Figure 5a). In particular, using pixel-wise loss in 2 distance has the worst performance. The involvement of 1 distance mitigates velocity errors but is slightly worse on MSE and SSIM of seismic error. Adding perceptual loss boosts the performance in all metrics by a clear margin. This shows perceptual loss is helpful to retain waveform coherence, which is correlated to velocity boundary, and validates our proposed loss function (combining pixel-wise and perceptual loss).   Feng et al. (2021) and employ the Marmousi velocity map as the style image to construct a low-resolution dataset. Table 3 shows the quantitative results on both datasets. Although our UPFWI achieves good results on Salt dataset with preserved subsurface structures, it has clearly larger errors than the supervised InversionNet. This is due to two reasons: (a) Salt dataset has a small amount of training data (120 samples), which is very challenging for unsupervised methods; (b) the variability between training and testing samples is small, providing a significantly larger favor to supervised methods than the unsupervised counterparts.  Table 6: Quantitative results of our UPFWI tested on seismic inputs with missing traces.
Other Network Architectures: We further conducted experiments by using Vision Transformer (ViT, Dosovitskiy et al., 2020) and MLP-Mixer (Tolstikhin et al., 2021) to replace CNN as the encoder. Table 4 further shows the quantitative results. Solid results are obtained for both network architectures, indicating our proposed method is modelagnostic. Visualization results are shown in Appendix A.4.

Robustness Evaluation:
We validate the robustness of our UPFWI models by two additional tests: (1) testing data contaminated by Gaussian noise and (2) testing data with missing traces. The quantitative results are shown in Table 5 and Table 6, respectively. We observe that in both experiments our model is robust to a certain level of noise and irregular acquisition. Visualization results are shown in Appendix A.4.

DISCUSSION
Our UPFWI has two major limitations. Firstly, it needs further improvement on a small number of challenging velocity maps where adjacent layers have very close velocity values. We find that the lack of supervision is not the cause as our UPFWI yields comparable or even better results compared to its supervised counterparts. Another limitation is the speed and memory consumption for forward modeling, as the gradient of finite difference (see Equation 6) need to be stored for backpropagation. We will explore different loss functions (e.g. adversarial loss) and methods that can balance the requirement of computation resources and the accuracy in the future work. We believe the idea of connecting CNN and PDE to solve full-waveform inversion has potential to be applied to other inverse problems with a governing PDE such as medical imaging and flow estimation.
Data-driven Methods: Recently, some researchers employed neural networks to solve FWI. Those methods can be further divided into supervised and unsupervised methods.
Supervised: One type of supervised methods require labeled samples to directly learn the inverse mapping, and they can be formulated as: where p denotes the seismic measurements, v is the velocity map, θ represents the trainable weights in the inversion network g θ (·), f (·) is the forward modeling, and L(·, ·) is a loss function. One example of supervised methods is the fully connected network proposed by Araya-Polo et al. (2018). Wu & Lin (2019) developed an encoder-decoder structured network to handle more complex velocity maps. Zhang & Lin (2020) adopted GAN and transfer learning to improve generalizability. Li et al. (2020) designed SeisInvNet to solve misaligned issue when dealing sources from different locations. In Yang & Ma (2019), a U-Net architecture was proposed with skip connections. Feng et al. (2021) proposed a multi-scale framework by considering different frequency components. Rojas-Gómez et al. (2020) developed an adaptive data augmentation method to improvegeneralizability. Sun et al. (2021) combined the data-driven and physics-based methods and proposed H-PGNN model.
Another type of supervised methods GANs to learn a distribution from velocity maps in training set as a prior (Richardson, 2018;Mosser et al., 2020). They can be formulated as: v(z * ) = g θ * (z * ) s.t. z * (p) = arg min z L(f (g θ * (z)), p), where Φ v is a training dataset including numerous velocity maps. z and α i are tensors sampled from the normal distribution. The iterative optimization is then performed on z to draw a velocity map sampled from the prior distribution.
Unsupervised: The existing unsupervised methods follow the iterative optimization paradigm and perform FWI per sample. They employ neural networks to reparameterize velocity maps. The networks serve as an implicit regularization and are required to be pretrained on an expert initial guess. Those methods can be formulated as: where a is a random tensor. Different network architectures have been proposed including CNNdomain FWI (Wu & McMechan, 2019) and DNN-FWI (He & Wang, 2021). Zhu et al. (2021) developed NNFWI which does not need pretraining ahead, but the initial guess is still required to be fed into the PDE with estimated velocity maps.

CONCLUSION
In this study, we introduce an unsupervised method named UPFWI to solve FWI by connecting CNN and forward modeling in a loop. Our method can learn the inverse mapping from seismic data alone in an end-to-end manner. We demonstrate through a series of experiments that our UPFWI trained with sufficient amount of unlabeled data outperforms the supervised counterpart on our dataset. The ablation study further substantiates that perceptual loss is a critical component in our loss function and has a great contribution to the performance of our UPFWI.
reflection coefficient chosen to be 10 −7 . With absorbing layers added, Equation 6 can be ultimately written as    Figure 13: Results of adding Gaussian noise to CurvedFault. The model is trained on the clean data (without noise) and tested on different levels (PSNR) of Gaussian noises. Similar to the results of FlatFault, our method is robust to the noise although slight degradation is observed when noise level increases.  Figure 14: Results of randomly missing traces on FlatFault. The model is trained on the clean data (without missing traces) and tested on multiple missing rates from 5% to 25%. Our method is robust to the missing traces. Although the higher missing rate leads to shifts in velocity values, the geological structures are well preserved.  Figure 15: Results of randomly missing traces on CurvedFault. The model is trained on the clean data (without missing traces) and tested on multiple missing rates from 5% to 25%. Similar to the results of FlatFault, our method is robust to the missing traces. Although the higher missing rate leads to shifts in velocity values, the geological structures are well preserved.
Additional experiments to investigate generalization. We conducted two additional experiments: (1) training our model on the CurvedFault dataset and further testing on the FlatFault dataset (visualization results are listed in Figure 16, and quantitative results are shown in Table 7); (2) testing our model on time-lapse imaging problems (visualization results are listed in Figure 17). The results demonstrate that our proposed model yields generalization ability to a certain degree.  Table 7: Quantitative results of our UPFWI models evaluated on FlatFault.
Ground Truth Prediction Figure 16: Results on generalization across datasets. The test is performed on FlatFault by applying a UPFWI model that is trained on CurvedFault dataset. Although the artifact is not negligible, the fault structures and velocity values are well preserved. This demonstrates that our model has generalizability to a certain degree.

Ground Truth
Prediction Ground Truth Prediction t = 0 t = 1 t = 2 t = 3 Figure 17: Results on generalizability over geological anomalies. The test is performed on a dataset where we add additional geological anomalies to simulate time-lapse imaging problems. The velocity maps containing those anomalies are not included during training. However, our model captures the spatial and temporal dynamics of anomalies in prediction. This demonstrates that our model has generalizability to a certain degree.