Discussions on the Processing of the Multi-Component Seismic Vector Field

: Multi-component seismic data contain a great deal of vector ﬁeld information that reﬂects the situation of the underground medium. However, the processing methods used for multi-component seismic data are still being developed, and e ﬀ ectively retaining and using this information is the di ﬃ culty and the focus of the task. Currently, the main-stream processing techniques of multi-component seismic data treat the individual components independently as a scalar ﬁeld; in this way, they do not excavate the vector features of the waveﬁeld, thus restricting the potential utilities of the e ﬀ ective information. Research into processing methods that are suitable for use with the vector ﬁeld, which can better retain and use the orientations and the relative amplitude relationship between multi-component seismic data, is urgently needed and represent an important direction for the current development of multi-component seismic data processing techniques. In this paper, we introduce and summarize several existing vector pre-processing techniques, including polarization ﬁltering, de-noising using vector order statistics, group sparse representation, and vector separation of compressional waves and shear waves, to help scholars develop more e ﬀ ective vector ﬁeld processing methods and to promote the development of vector processing techniques for multi-component seismic data.


Introduction
The three-component seismic data in a three-dimensional Descartes system simultaneously record the linear motion of the particle in the vertical and horizontal directions as well as retain a relatively complete record of seismic wave translations. They not only simultaneously provide the travel time (velocity), amplitude, and frequency information of the compressional (P) and shear (S) waves but also provide the relative amplitude relationship of the same wave among different components, as well as the amplitude, frequency, phase, and travel time differences for different waves. The information of these vector wavefields would provide a more detailed and accurate characterization of the structure, lithology, fluid saturation, pore pressure, and fractures of an underground medium [1]; however, because of the limitations of current methods and processing techniques of seismic data, such information cannot be fully obtained. The advantages of the current multi-wave and multi-component seismic techniques are limited to the following [2][3][4][5]: (1) inversions of travelling time and waveform for the P-wave only [6,7]; (2) shear-wave birefringence analysis [8]; (3) compensation of a defect caused by the difficulty of imaging with a PP-wave under a gas cloud, a salt dune, or in a low-acoustic impedance area using PS-wave imaging [9,10]; (4) anisotropy and fracture identification [11]; and (5) prediction of reservoir lithology and fluids [12]. In particular, except for the prediction of the fast and slow

The Polarization Filtering Method
Different types of seismic waves have different polarization manners during propagation. The surface wave exhibits elliptical polarization characteristics, the P-wave exhibits linear polarization characteristics in the direction of the wave propagation, while the S-wave exhibits linear polarization characteristics perpendicular to the direction of the wave propagation, which provides the theoretical basis for identifying and separating various wave types with their differences in polarization characteristics. Polarization filtering is usually used to separate the body and surface waves in the field of earthquake research [36,37] and suppress random noise [38,39], while in exploration seismology, it is mostly used for surface-wave suppression and seldom used for the attenuation of random noise.
The earliest method of polarization analysis occurred in the time domain [40]. Thereafter, Jurkevics [41] and Chen et al. [42] conducted detailed studies. The time-domain polarization analysis constructs a covariance matrix using the multi-component seismic data within a certain time window and computes the polarization property of the wave within this time window using the characteristic value and characteristic vector of the constructed covariance matrix. The difficulty with time-domain polarization filtering is the selection of the time-window length because different lengths result in different polarization parameters. Therefore, Diallo et al. [43] and Ma [44] proposed to adaptively select the length of the time window according to the instantaneous frequency of the signal; Rene et al. [45], Morozov and Smithson [46], Vidale [23], Schimmel and Gallart [47], and Lu et al. [26] proposed the method of Complex Trace Analysis (CTA) based on the Hilbert transform to obtain the instantaneous polarization properties of the multi-component seismic signals.
However, these methods of time-domain polarization analysis cannot differentiate the various wave types that are mutually superposed in time. Park et al. [48] and Reading et al. [27] proposed the frequency-domain polarization filtering method to differentiate the waves that are superposed in time but with different frequency components. However, the body waves and surface waves would be superposed in both time and frequency; therefore, polarization analysis methods based on a time-frequency transform were developed, including the continuous wavelet transform-based method [28,49], the discrete wavelet transform-based method [50], and the S transform-based methods [51,52]. To solve the problem of the inaccurate estimation of the polarization parameter caused by the phase distortion of the conventional wavelet transform coefficient, Galiana-Merino et al. [53] proposed to use the stationary wavelet packet transform to estimate the instantaneous polarization parameters.
Methods of polarization filtering in the time-frequency domain simultaneously use the travel time, frequency, and polarization features of seismic waves, but these methods belong to the single-channel processing mode. Another significant difference between the surface wave and body wave, namely, the apparent velocity, cannot be used. Therefore, Wang et al. [29] proposed to develop a complex seismic data gathering technique by taking the radial component as the real part and the vertical component as the imaginary part, before applying the polarization analysis in the t-f-k-domain for surface-wave attenuation. If the surface wave and body wave are different in either time, frequency, or apparent velocity, the surface waves can be automatically identified and eliminated according to the polarization rate in the t-f-k domain, which fully reflects the advantages of the vector processing method. An example is given here to exhibit the effectiveness of the t-f-k domain polarization filtering method. The synthetic two-component seismic data are shown in Figure 1 where Figure 1a is the vertical component, and Figure 1b is the radial component. The dominant frequency is 30 Hz for P-waves and 22 Hz for PS-waves. Considering that the frequency range of the ground roll is 5-20 Hz, it shows considerable overlap between the frequency of the PS-waves and the ground roll. The most challenging issue is that the apparent velocities of the ground roll and the reflected wave in far offset are very close. A comparison of the t-f-k domain polarization filtering with the f-k filtering and time-frequency domain polarization filtering would convince that the t-f-k domain polarization filtering performs more effectively. As shown in Figure 1c,d, the results of f-k filtering work well but not perfectly. The ground roll is not completely removed from the filtered result. Additionally, some body waves appear in the removed noise profiles, as shown in Figure 1e,f; this is because the apparent velocity of these body waves is close to the apparent velocity of the ground roll. Figure 2a,b show the filtered results of the time-frequency polarization filtering. Because the body wave and ground roll overlap both in the time domain and frequency domain, they cannot be effectively separated in the time-frequency domain. The t-f-k polarization filtering is implemented in the t-f-k domain. Three f-k slices (from time 150, 650, and 950 ms) of the t-f-k transform of the gathered complex seismic data are shown in Figure 3a,c,e. These figures show a meaningful phenomenon; the body waves have equal energy while the ground roll has discrepant energy in the positive and negative frequency planes. The ground roll allocated more energy in the positive frequency plane. The f-k slices taken from different times are different from each other; therefore, the computed filtering windows should change with time, as shown in Figure 3b,d,f. With the automatically identified filtering windows, the filtered results of the t-f-k polarization filtering are shown in Figure 4a,b. Furthermore, the removed ground roll is shown in Figure 4c,d. The test results indicate that the t-f-k polarization filtering method can suppress the ground roll effectively and guarantee the body wave does not drop out.
Appl. Sci. 2019, 9,1770 4 of 17 this is because the apparent velocity of these body waves is close to the apparent velocity of the ground roll. Figure 2a,b show the filtered results of the time-frequency polarization filtering. Because the body wave and ground roll overlap both in the time domain and frequency domain, they cannot be effectively separated in the time-frequency domain. The t-f-k polarization filtering is implemented in the t-f-k domain. Three f-k slices (from time 150, 650, and 950 ms) of the t-f-k transform of the gathered complex seismic data are shown in Figure 3a,c, and e. These figures show a meaningful phenomenon; the body waves have equal energy while the ground roll has discrepant energy in the positive and negative frequency planes. The ground roll allocated more energy in the positive frequency plane. The f-k slices taken from different times are different from each other; therefore, the computed filtering windows should change with time, as shown in Figure 3b,d and f. With the automatically identified filtering windows, the filtered results of the t-f-k polarization filtering are shown in Figure 4a,b. Furthermore, the removed ground roll is shown in Figure 4c,d. The test results indicate that the t-f-k polarization filtering method can suppress the ground roll effectively and guarantee the body wave does not drop out.

The Vector Wavefield Filtering Method Based on Vector Order Statistics
Many studies have revealed that multivariate order statistic filters perform better in color image processing because they treat the color image as vector data whereas other filtering methods treat the color image as scalar data [54]. These multivariate order statistic filters can be classified into several

The Vector Wavefield Filtering Method Based on Vector Order Statistics
Many studies have revealed that multivariate order statistic filters perform better in color image processing because they treat the color image as vector data whereas other filtering methods treat the color image as scalar data [54]. These multivariate order statistic filters can be classified into several categories according to the ordering method for multivariate data, such as the multi-channel α-trimmed mean filter, multichannel Modified-Trimmed Mean (MTM) filter [55], Vector Median Filter (VMF) [56,57], and vector direction filter [58]. A study by Pitas et al. [59] indicated that the VMF is more applicable to suppressing noise that has a long-tail distribution, such as negative index noise and salt-and-pepper noise; the MTM filter is more suitable for suppressing noise that has a short-tail distribution, such as Gaussian white noise.
In recent years, Huo et al. [60] applied the VMF to seismic data processing and separated blended seismic data; Liu [61] used the VMF to eliminate direction anomalies of the local dip vectors of seismic events. In these two applications, however, the filter still processes every single component of the seismic data, and the vector is formed from the mathematical construction of some sampling points of adjacent traces or adjacent time, rather than a vector in physical space. Wang et al. [30] and Xun et al. [31] applied an extension of the vector order statistics-based filter to multi-component seismic data. First, according to the correlation between the vector signals of adjacent traces, the optimum local direction of the event was examined, and then, the vector points along the event were selected for vector order statistics-based filtering. Because every sampling point of a three-component detector records a vector, it is very intuitive, and there are practical advantages to adopting the vector order statistics-based filter to process three-component seismic data. Tests of the model data and real data showed that the noise-suppressing result using vector order statistics-based filtering is better than when using the signal-component approach [30]. The solid lines in Figure 5 show some hodograms of two-component seismic signal segments taken from the de-noised results with different methods, and the dashed lines are the hodogram curves of the uncontaminated seismic data. The first, second, and third columns are the results of low-pass filtering, single-component filtering, and vector order statistics-based filtering, respectively. From top to bottom, there are four seismic signal segments which are taken from different detectors and different events. We can see that vector order statistics-based filtering maintains the vector information of multicomponent seismic data more effectively than the single-component filtering.

Multi-Component Joint De-Noising Based on the Group Sparsity Representation
Expressing the time-domain noise-containing signals in other domains through mathematical transforms is the approach often used in de-noising seismic data and wavefield reconstruction. The mathematical transforms often used include Fourier transform, S transform, wavelet transform, and curvelet transform. These signal representation methods based on specific mathematical transformations can effectively represent the signal of a corresponding structural model according to the structural characteristics of the individual transform dictionaries. However, real seismic signals are very complex; therefore, it is difficult to achieve the sparse representation of seismic signals by only selecting a single transform method. Many overlaps in the representation coefficients of the noise and the effective signal still exist when using a single transform method. We combine the dictionary of multiple transforms to compose an over-complete dictionary or construct a new

Multi-Component Joint De-Noising Based on the Group Sparsity Representation
Expressing the time-domain noise-containing signals in other domains through mathematical transforms is the approach often used in de-noising seismic data and wavefield reconstruction. The mathematical transforms often used include Fourier transform, S transform, wavelet transform, and curvelet transform. These signal representation methods based on specific mathematical transformations can effectively represent the signal of a corresponding structural model according to the structural characteristics of the individual transform dictionaries. However, real seismic signals are very complex; therefore, it is difficult to achieve the sparse representation of seismic signals by only selecting a single transform method. Many overlaps in the representation coefficients of the noise and the effective signal still exist when using a single transform method. We combine the dictionary of multiple transforms to compose an over-complete dictionary or construct a new dictionary according to the signal characteristics. Through searching, we can only use a small number of elements in the over-complete dictionary for the sparse representation of complex signals [62,63]. Because many fast and effective algorithms of sparse representation have been proposed in recent years, such as the matching pursuit algorithm [64], the orthogonal matching pursuit algorithm [65], and the base tracking algorithm [66], sparse representation has been widely applied in seismic data processing. Applications include the reconstruction of the wavefield [67], surface-wave attenuation [68,69], and de-convolution [70]. To retain the relative amplitude relationship among the different components of multi-component seismic data after de-noising, Rodriguez et al. [32] developed and applied the time-frequency transform based on the group sparsity constrained representation to de-noise three-component micro-seismic data. They believe that because the micro-seismic events recorded by the three-component detector have the same arrival time and frequency component for the three components, the coefficient of sparse representation for the effective signal of the three components in the time-frequency domain should have the same distribution mode. In other words, the coefficients of the three components at the same time-frequency point are either all zero or nonzero. With the coefficients of three-component data at the same time-frequency point as one group, we can use the concept of group sparsity proposed by Yuan and Lin [71], Fornasier and Rauhut [72], and Eldar [73] to process the multicomponent seismic data. The three components x, y, and z can form an N × 3 matrix (where N is the number of time sampling points): The group sparsity constrained representation in an over-complete dictionary Φ can be expressed as the following mathematical optimization problem: whereÂ is the coefficient of sparse representation, and each column corresponds to a component. The symbol • p,q represents the mod of the matrix, which is defined as The optimization problem of Equation (2) can be solved by adopting the fast iterative shrinkage-threshold algorithm (FISTA) proposed by Beck and Teboulle [74]. While constructing the dictionary, Rodriguez et al. used the complex Ricker wavelet as the mother wavelet. The complex wavelet dictionary not only improves the sparsity of the representation coefficients but also retains the consistency between the representation coefficients for waves with different phases.
We applied the multi-component joint de-noising method based on the group sparsity constrained representation to synthesized seismic data. Figure 6a shows the three components of the synthesized data, and Figure 6b shows the data contaminated with Gaussian white noise. The first waveform in the y component is completely submerged in the noise. We used the multi-component joint de-noising method based on the group sparsity constrained representation and the single-component de-noising method based on the sparse representation, separately, for the de-noising processing. The results are shown in Figure 6c,d. We can see that the multi-component joint de-noising method can effectively attenuate random noise, and it can preserve weak signals more effectively than the single-component de-noising method. Beyond that, the multi-component joint de-noising method has incomparable advantage in reserving the vector information of the multi-component seismic signal, which is useful for successive processing.

The Vector Separating Method of the P-wave and S-wave
The separation of the P-wave and S-wave cannot only be used as a special de-noising method to eliminate the mutual interference between the P-wave and S-wave, but also as an important tool for the extraction of the pure wave and the reconstruction of the wavefield. For the conventional and scalar processing systems, separation of the P-wave and S-wave is the foundational step, while for the vector processing methods, separation of the P-wave and S-wave can reduce the complexity of the vector field data, which is favorable for the analysis of fast and slow S-waves and can increase the imaging quality. It should be emphasized that the subsequent multi-wave joint application requires the separated P-wave and S-wave, though the data must retain the original vector characteristics.
The wavefield separation method based on the Radon transform independently processes the single-component seismic data, and cannot retain the vector characteristics of multi-component seismic data [75]. Hu et al. [76] separated P-waves and S-waves according to the polarization direction of multi-component seismic data, but the separated P-waves and S-waves retrograded to the scalar field. Yao et al. [77], Sun et al. [78], and Yan and Sava [79] proposed a series of separation methods for P-waves and S-waves based on the divergence and curl of different media with homogeneous, inhomogeneous, and isotropic properties, but the separated P-waves became scalars. Moreover, the divergence and curl operation causes distortions in the phase and amplitude for the wavefield after separation; therefore, the results need to be corrected.
To ensure, after separation, that the P-waves and S-waves still retain the relative vector characteristics as in the original input data, Li [35] proposed a wave equation method for the separation of P-waves and S-waves in the frequency-wavenumber domain, which first separated the P-wavefield and S-wavefield from the total wavefield using the divergence and curl operation according to the Helmholtz decomposition principle. Thereafter, the method projected two wavefields in the vertical and horizontal directions by multiplying them with the normalized wave numbers in the vertical and horizontal directions. This method barely changes the phase and amplitude information of the P-waves and S-waves. The derived P-waves and S-waves are still two components of x and z and retain the vector information of the original wavefield. However, this method is only applicable to the homogeneous isotropic surface medium with known P-velocity and S-velocity. For an inhomogeneous isotropic medium, Wang et al. [20] proposed a vector separation method of P-waves and S-waves based on the wavefield extrapolation. In this method, the P-and S-wave-separated elastic equation was first used to downward extrapolate the elastic wavefield from the surface to a reference depth, and the P-waves and S-waves were separated in the process of downward extrapolation. Only the P-waves were saved at the reference depth; then, the elastic wave equation was used to upward extrapolate the saved P-wave to the surface to obtain the separated P-wave recording. The P-waves were subtracted from the original recording of the total wavefield to obtain the separated S-wave recording. Li et al. [80] proposed a similar vector separation method for P-waves and S-waves, but they first used the elastic wave equation to downward extrapolate the elastic wavefield from the surface to a reference depth. Thereafter, they used the P-and S-wave-separated elastic equation to extrapolate the wavefield upward to the ground surface to obtain the separated P-waves and S-waves. Therefore, the P-waves and S-waves are separated naturally during the upward extrapolation, and their phase and amplitude are the same as before separation. The P-and S-wave-separated elastic equation used by these two vector separation methods is the equivalent form of the conventional coupled elastic wave equation proposed by Ma and Zhu [81], and it is only applicable to an isotropic medium. Using the vector decomposition method for an elastic wavefield developed by Zhang and McMechan [19], the vector separation method for P-waves and S-waves could be extended to an anisotropic medium.
The abovementioned vector wavefield separation methods all imply that the P-waves and S-waves are orthogonally polarized, and this assumption is generally not valid for an anisotropic medium [82]. For the wavefield separation problem of non-orthogonal P-waves and S-waves in an anisotropic medium, Lei [33] and Lu et al. [34] discussed this problem from the perspective of natural earthquakes and exploration seismology and provided a solution. The basic idea of the solution is to apply a rotational transform to transform the vertical and horizontal components based on the Cartesian coordinate system to the affine coordinate system composed of the wave vector directions of realistic P-waves and S-waves, thereby achieving separation of the P-wave and S-wave. Although the separated P-waves and S-waves given in their paper are scalar fields, the wave vector directions of the different waves in the Cartesian coordinate system have been calculated in the separation process, and according to the directions of the wave vectors, the scalar field could be restored to the vector field. Because this method can calculate the actual wave vector direction at every sample point, it provides a new approach for the de-noising technique based on the characteristics of wave vectors.

Suggestion
Great achievements have been obtained in vector field processing with the unremitting efforts of many scientific researchers [30,32,35,83]. However, generally speaking, the existing methods of multi-component seismic data vector processing are based on a relatively simple theoretical foundation, and so they can only process simple seismic signals. The processing capability is not sufficient for complex seismic signals modulated by complex inhomogeneous anisotropic media.
The quaternion theory [84] and multilinear algebra theory [85] developed in recent years are very suitable for processing high-dimensional vector signals, and they have been successfully applied in the fields of electromagnetic physics [86], multi-channel digital communication [87], and color image processing [88,89]. Therefore, it will be of significant interest to study more sophisticated and systematic methods of multi-component seismic vector field processing based on the quaternion or multilinear algebra theory, thereby exploiting the potential capability of multi-component seismic data in solving the problem of complex anisotropic media.

Conclusions and Perspectives
The multi-component seismic wavefield is a physical vector field. To fully demonstrate the advantages of multi-component seismic data, we first have to retain the seismic vector field information during pre-processing, at least without interrupting the relative spatial direction and amplitude relationship between the multiple components. Secondly, we need an advanced technique to convert the kinematic and dynamic information hidden in the vector wavefield into the elastic parameters of the underground medium, such as velocity anisotropy, stress distribution, pore pressure, fluid saturation, and fracture size, which directly relates to oil and gas reservoirs. In this paper, we introduced several pre-processing methods for a multi-component seismic data vector field. The examples of real multi-component seismic data indicated that the vector processing technique could better avoid the distortion of vector field information and give more accurate predictions than the scalar processing technique.
Simply because studies on the extraction and utilization of seismic vector information are still in their infancy, vector processing techniques are currently outside of mainstream industrial interests. However, we believe this field will continue to develop, taking a key role in the future of fossil energy exploration, earthquake research, and micro-seismic monitoring.