{'fraction_non_alphanumeric': 0.0365400687157582, 'fraction_numerical': 0.015600334292877705, 'mean_word_length': 4.577162092180218, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, '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}}
CNN-tuned spatial filters for P- and S-wave decomposition and applications in elastic imaging

P- and S-wave decomposition is essential for imaging multi-component seismic data in elastic media. A data-driven workflow is proposed to obtain a set of spatial filters that are highly accurate and artifact-free in decomposing the P- and S-waves in 2D isotropic elastic wavefields. The filters are formulated initially by inverse Fourier transforms of the wavenumber-domain operators, and then are tuned in a convolutional neural network to improve accuracy using synthetic snapshots. Spatial filters are flexible for decomposing P-and S-waves at any time step without performing Fourier transforms over the entire wavefield snapshots, and thus are suitable for target-oriented imaging. Snapshots from synthetic data show that the network-tuned spatial filters can decompose P- and S-waves with improved decomposition accuracy compared with other space-domain PS decomposition methods. Elastic reverse-time migration using P- and S-waves decomposed from the proposed algorithm shows reduced artifacts in the presence of a high velocity contrast.


I. INTRODUCTION
P -and S-waves coexist in seismic waves propagating through the solid earth. The causes, polarization directions, and propagation velocities, of P-and S-waves are different, and thus they carry different information from the subsurface at a variety of scales. For example, in seismic imaging for oil and gas exploration, P-and S-wave reflections can be imaged independently. In the presence of gas clouds and fractures, the P-waves are strongly attenuated, while Swaves are less affected, and S-wave images provide valuable structural information from reservoirs [1], [2]. In computational seismology, most imaging methods for earthquake and microseismic events involve picking the P-and S-wave arrival times. Artman et al. [3] and Yang and Zhu [4] propose time-reversal imaging of earthquake ruptures and hydraulic fractures by separating and propagating the recorded P-and S-waves back into the earth, and applying the elastic imaging conditions. All of these applications require separating the Pand S-waves accurately and efficiently.
One approach to separate and utilize both P-and S-wave data is the ray-based method. Examples include Kirchhoff migration [5], [6], [7], in which the PS separation is implicit by using P-and S-velocity models separately for ray-tracing to compute the image times. Ray-based methods are efficient, but their high frequency assumption makes ray theory unable to fully describe wave phenomena, especially for complicated models [8]. Another category is wave-equation based solutions [9], [10], [11], which reconstruct the full elastic receiver wavefields from boundary conditions [12], and separates the P-and S-waves before, or as a part of, application of the image condition [13].
The divergence and curl operators, which are based on Helmholtz decomposition [14], are widely applied in PS separations. However, the phases and amplitudes of the separated P-and S-waves are distorted by the divergence and curl operators, causing problems for true amplitude imaging and subsequent interpretations [15]. To solve this problem, PS wavefield decomposition is proposed, which preserves the vector components of the input elastic wavefield [16], [17], [18].
There are two groups of algorithms for PS decomposition depending on the domain of calculation. In the wavenumber domain, elastic wavefields are projected to the P-and Spolarization directions [17]; Zhu [19] decomposes the wavefields by solving a Poisson's equation. Both of these methods involve Fourier transforms. However, it is expensive, or even prohibitive, to perform Fourier transforms on large (3D) datasets. Xiao and Leaney [20] introduce an auxiliary stress wavefield into the elastodynamic wave equations to obtain decoupled wavefields, but the decomposition is embedded in the elastic extrapolations (decoupled propagation), and artifacts are likely to be generated at large velocity contrasts in the model [18]. Yan and Sava [13] propose spatial filters for separating P-and S-waves in transversely isotropic media with a vertical axis (VTI). The spatial filters don't require Fourier transforms, but the accuracy is limited by the truncation of the spatial filters.
PS separation and decomposition can also be achieved with deep learning algorithms. Wei et al. [21] reconstruct P-waves in VSP data using conditional generative adversarial network. Wei et al. [22] extend the algorithm to separate P-and Swaves in multi-component VSP data using multi-task learning. However, these methods need a lot of training samples from a variety of velocity models to compensate the lack of physical information. The spatial filters are convolutional operators, which are widely used in the computer vision community [i.e. convolutional neural networks (CNNs)] [23]. Thus, CNNs are suited for separating P-and S-wave snapshots. Kaur et al. [24] apply a generative adversarial network (GAN) to match the predicted P-waves with the decomposed ones from lowrank approximation in anisotropic media; the network needs to be re-trained to apply to different models. We propose to formulate the PS decomposition as a single-layer CNN with the coupled wavefield snapshots as inputs and the decomposed P-wave snapshots as outputs. The spatial filters are tuned to optimize PS decomposition accuracy by training in a training set. Different from the existing methods, our network only needs to be trained once on a complex model and can be applied to any models, because the physics in PS decomposition are embedded in the network structure. The trained network can decompose elastic wavefields in different models without re-training.
This paper is organized as follows, we first illustrate the wavenumber domain operators for P-and S-wave decomposition, and their relation to spatial filters. Then we describe a CNN structure to tune the spatial filters. Wavefield snapshot decomposition and elastic reverse time migration are performed to document the improvements in accuracy by comparing the decomposition results with those from two other space-domain PS decomposition methods. To limit the scope, we focus on 2D isotropic media, although PS decompositions in anisotropic and anelastic media are also possible [25], [26], [27].

II. PS DECOMPOSITION
Different from the commonly used divergence and curl operators [14] which distort the amplitudes and phases during wavefield separation, PS decomposition preserves all the components of the elastic vector wavefield. In the wavenumber domain, the explicit equations to obtain 2D decomposed Pwaves are [17] where K x and K z are the components of the normalized wavenumber vector K, which represents the P-wave polarization direction in isotropic media. U x and U z are the space-domain elastic wavefield components with coupled Pand S-waves, and a tilde over the wavefield indicates the wavenumber-domain representation. The S-waves can be decomposed in a similar way, or by subtracting the P-waves from the coupled wavefields, component-by-component-bỹ in either the space or the wavenumber domain. Thus only the P-wave decomposition is discussed and demonstrated in the following parts of this paper.

III. SPATIAL FILTERS FOR P-AND S-WAVE DECOMPOSITION
In the 2D equations 1 and 2, three wavenumber operators K 2 x , K 2 z and K x K z are involved; the K 2 z is the transpose of K 2 x when the spatial intervals along x and z directions are the same. K 2 x and K x K z are plotted in Figures 1a-1b, respectively. The wavenumber domain operators can be transformed into the space-domain by inverse Fourier transforms to formulate the spatial filters [13]. The equivalent expressions of equations 1 and 2 in the space-domain are and where L x , L z and L xz are filters that are designed to decompose the elastic wavefields. The square brackets represent spatial filtering of the wavefields. Different from the anisotropic media, the polarizations of P-and S-waves in isotropic wavefield have much simpler relations with the wavenumbers, they do not depend on the model parameters (V P , V S or ρ). Thus the spatial filters are stationary and universal for different models. The spatial filters are trained only once using synthetic labelled data pairs, and the trained filters can be applied to decompose P-and S-waves in any isotropic wavefields. Full representations of the wavenumber domain operators are infinitely large 2D patches in the space-domain which make them infeasible to be applied. Yan and Sava [13] propose to truncate the filters to balance the computation efficiency and the separation accuracy. A set of truncated L x and L xz filters of size 15 × 15 are shown in Figures 2a-2b (L z is the transpose of L x , and thus is not plotted).

IV. TUNING THE SPATIAL FILTERS WITH CNN
Truncations of the filters reduce the accuracy and completeness of the decomposition. To improve the accuracy without increasing the size of the filters, we propose to tune the spatial filters L x and L xz by training in a CNN. Following equations 5 and 6, a CNN is constructed to obtain the P-waves U P from the coupled wavefield snapshot U of the form The two components (U x and U z ) of the elastic wavefield snapshot U are input to the network as two channels, and the output is expected to be the decomposed 2-component P-wave snapshot U P . The CNN for PS decomposition is shown in Figure 3. Two spatial filters (L x and L xz ) are to be trained. Most neural networks in deep learning are black-boxes and cannot be explained. However, the CNN which we use contains only one CNN layer, which can be simply interpreted as spatial filters. The physical meaning of our CNN is a spatial representation of a Fourier domain operator, which is similar to the work of Yan and Sava [13]. Thus our network does not contain pooling or non-linear activation functions because of the well understood physics that is involved. Another difference is that, instead of initializing the spatial filters with random numbers, they are obtained by transforming the wavenumber domain operators in Figure 2a and 2b, and are set as initial values before training.
The network is trained by minimizing the loss function where U P andÛ P are the predicted and ground truth Pwaves, respectively. | · | calculates the modulus arguments of vector difference. N is the number of grid points in a wavefield snapshot. The filters are tuned in a CNN training process. We expect the training to reduce the errors that are associated with the filter truncations. We generate a training set with 800 samples, which contains coupled and decomposed P-wave snapshot pairs. The coupled wavefield snapshots are generated by simulating wave propagation in a portion of the Sigsbee model [28] (Figure 4) using an eighth-order in space, second-order in time, stress-particlevelocity, staggered-grid, finite-difference solution [29]. An explosive source with a 10 Hz Ricker wavelet is placed at (x, z) = (1.28, 0.02) km. The spatial intervals in both x and z directions are 10 m, and the time interval is 1 ms. The 800 2-component snapshots used are captured at random time steps from T = 0.0 s to T = 1.4 s. Although reflections from model edges can also be separated into P-and S-waves, we use Fourier transforms to generate the ground truth (decomposed P-and S-wave snapthots), so they may have wrap-around errors at the model edges. Thus, convolutional perfectly matched layer (CPML) absorbing boundary conditions [30] are used on all four boundaries to reduce unwanted reflections.
The training doesn't need a massive training set, because (1) approximate spatial filters (Figures 2a and 2b) are set as initial filters which are close to the optimized filters; (2) the CNN structure in Figure 3 imposes a physical contraint on the training; (3) only 2 filters (L x and L xz ) need to be tuned. After training, the CNN-tuned spatial filters can be applied to any wavefield snapshot that has the same spatial intervals as the training set.
The ground truth P-waves are the decomposition results  using the wavenumber decomposition algorithm (equations 1 and 2), which is, by far, the most accurate algorithm for Pand S-wave decomposition [17]. Samples of the training set are plotted in Figure 5. The training set is sufficiently diverse, as the Sigsbee model contains a wide range of structures and velocities.
The size of the spatial filter should be at least one wavelength to maximize the resolution, but large spatial filters increase the computation cost. Three sets of spatial filters are trained with dimension 9 × 9, 15 × 15 and 21 × 21. The networks are trained using an Adam optimizer [31] in the Pytorch (www.pytorch.org) platform. Figure 6 shows that losses decrease as the training proceeds, and that the decomposition accuracy of all these sizes of filters is improved during training. As expected, spatial filters with larger sizes can reach lower losses. Big improvements can be made when the filters expand from 9 × 9 to 15 × 15, but improvements are less significant for filter size 21 × 21 ( Figure 6). In the following tests, only the updated L x and L xz filters of the dimension 15 × 15 are used, which are plotted in Figures 7a  and 7b, respectively. Figures 7c and d show the differences between the filters before and after the tuning.

V. SYNTHETIC TESTS AND COMPARISONS
We perform tests on synthetic data using the CNN-tuned spatial filters (CNN-SF), and the results are compared with the decomposed P-waves from untuned spatial filters (SF), and those from decoupled propagation (DP) [16]. DP involves solving an auxiliary equation along with the elastodynamic wave equations, and thus has to be applied at each time step from the beginning of the elastic wavefield extrapolation. Note that CNN-SF is tuned only from the Sigsbee model in the previous section, and can be directly applied to the following (different) models without re-training.

A. Two-layer example
The first test is performed on an isotropic two-layer model with V P = 3 km/s, V S =2.1 km/s, and ρ = 2.  g/cm 3 for the lower layer. The model has 10 m grid spacing in both the x-and z-directions. An explosive source with a 10 Hz Ricker wavelet is placed at (x, z) = (1.28, 0.9) km. A snapshot of the PS coupled x and z particle velocity components at t = 0.42 s is shown in Figure 8a, and the ground truth P-wave using equations 1 and 2 are shown in Figure 8b as benchmarks.
The decomposed P-waves using DP, SF and CNN-SF methods are shown in Figures 9a-9c, respectively. To measure the decomposition accuracy, we subtract the decomposition results (U P i ) in Figures 9a-9c from the ground truth (Û P i ) P-wave (Figure 8b), and the residuals are shown in Figures 9d-9f. DP has high accuracy at homogeneous parts of the medium, but it generates serious artifacts along the interface between layers, which are caused by the added auxiliary equation [18], indicating that the DP method can be applied only to smooth models. The spatial filter algorithms SF and CNN-SF don't have this problem. The residuals in Figure 9f are smaller than in Figure 9e, indicating that the CNN-SF achieves higher accuracy than the SF. We quantify the accuracy by The calculated accuracies of DP, SF and CNN-SF in this example are 92.2%, 90.1%, and 98.6%, respectively. The computation times of DP, SF, and CNN-SF are 0.01, 0.25 and 0.25 seconds, respectively, for each decomposition of the above model. DP has the best efficiency, but it involves solving an auxiliary equation along with solving the wave equations, thus decomposition with DP is not independent from wavefield extrapolation. CNN-SF and SF algorithms, on the other hand, don't require extrapolations, and they have the same filter size, and thus share the same computation time.

B. Marmousi-2 model example
The second test is performed on a portion of the Marmousi-2 model [32] (Figure 10). The model has 10 m grid spacing in both x-and z-directions. An explosive source with a 15 Hz Ricker wavelet is placed at (x, z) = (1.28, 0.04) km. A snapshot of the x and z particle velocity components at t = 0.21 s is shown in Figure 11.
We decompose the representative coupled wavefields in Figure 11 using the DP, SF and CNN-SF methods, respectively; the decomposed P-waves are shown in Figures 12a-12c, and their residuals are in Figures 12d-f, respectively. Results of DP contain artifacts along the velocity interfaces. The decomposition results from SF and CNN-SF methods don't have such artifacts, and the decomposition accuracy of CNN-SF is better than SF.
To have a comprehensive analysis of the decomposition accuracy, 575 snapshots of coupled wavefields at different time steps are captured from the wavefield simulation in the Marmousi-2 model, and we use DP, SF and CNN-SF to generate the decomposed P-waves separately from all the snapshots. The accuracy of the three algorithms with increasing time steps are plotted in Figure 13. As the wave propagates through the model, more converted S-waves are generated, and the accuracies of all three algorithms decrease. Because the Marmousi-2 model has more velocity interfaces than the two-layer model, the accuracy of DP decreases significantly as the wave propagates deeper into the model. CNN-SF has the best accuracy among the three algorithms at all time steps.

C. Application to elastic reverse-time migration
In elastic reverse-time migration [10], the source wavefield and receiver wavefield are extrapolated forward and backward in time, respectively. Both the source and receiver wavefields are decomposed into P-and S-wave components, and various vector-based imaging conditions [33], [34], [35] can be applied to generate images. We apply the 2D dot-product imaging condition for the PP image and for the PS image, where T is the number of discrete time steps during wavefield extrapolation; U P src , U P rec and U S rec are the decomposed P-waves from the source wavefield, and decomposed P-and S-waves from the receiver wavefield, respectively. All the wavefields in equations 10 and 11 are a) b) c) The Pluto model [36] is used in this test ( Figure 14). The original Pluto model has only the P-wave velocity defined. We approximate the S-wave velocity model by multiplying the Pwave values by 0.6 at each grid point, and the density has a constant value of 2.1 g/cm 3 . The grid increments are h = 10 m in both spatial dimensions and time increment dt = 1 ms. One hundred explosive sources with 15 Hz dominant frequency are initiated sequentially from (x, z) = (0.5, 0.0) km to (13.5, 0.0) km with 0.13 km separation in the x direction. 1361 receivers are placed from (x, z) = (0.0, 0.0) to (13.6, 0.0) km with 0.01 km separation. The stacked PP and PS images using DP for PS decomposition are shown in Figure 15a and 15b, respectively. The PS separation artifacts are visible along the high velocity contrasts in the zoomed-in sections of Figure 15c and d. Compare with Figures 16a-16d, which are the PP, PS and zoomed-in images using CNN-SF for PS decomposition; the image boundaries are more coherent than in Figure 15 as the result of applying CNN-SF for PS decomposition.

VI. DISCUSSION
In the above tests, the dimensions of the CNN-tuned filters is fixed to be 15 × 15. The accuracy of decomposition can be further improved by using larger filters at the cost of computation efficiency. The ground truth snapshots for training the CNN are the P-waves decomposed in the wavenumber domain, thus the accuracy of CNN-tuned filters can never surpass that of the wavenumber domain algorithm.
The computational complexity for PS decomposition using the 15 × 15 spatial filter set is 225N , while the cost for the wavenumber domain algorithm is 4N log(N ), thus the latter is also much more efficient than the former for small models. However, wavenumber domain algorithms require Fourier transforms which consume large RAM and have difficulty  being applied locally in the model, while the spatial filters are HPC-friendly and can be easily applied in parallel using multiple cores. Spatial filters are suitable for processing targetoriented data sets, where PS decomposition can be performed locally on a portion of the wavefield and thus further reducing the cost.
The idea of using CNN-tuned spatial filters for PS decomposition can be extended to anisotropic media. In that case, the spatial filters are non-stationary across different positions within the model, because in anisotropic wavefields, the polarization directions of P-and S-waves depend on the model parameters [13], and several sets of filters can be tuned separately by training with parameters from representative points in the model. The proposed CNN-tuned spatial filters can also be extended to 3D, for which the number of filters increases to 6 (L x , L y , L z , L xy , L xz and L yz ) which need to be tuned as 3D convolutional operators.

VII. CONCLUSIONS
We propose network-trained spatial filters for P-and S-wave decomposition in isotropic elastic wavefields. The accuracy of decomposition using spatial filters is improved after CNNtuning with no extra computation costs. Migration tests with synthetic data show that the proposed spatial filters do not generate artifacts along high velocity contrasts as in the decoupled propagation method. CNN-tuned spatial filters are suitable for parallel computation and target-oriented imaging.   Figure 15, the artifacts are much reduced.