1 INTRODUCTION

The receiver function(RF)method(Langston, 1979)has become a general technique for imaging the crust and upper mantle structure. Until now, however, the most interpretations of RFs are still based on the laterally stratified earth model(Owens et al., 1984; Duecker and Sheehan, 1997; Zhu and Kanamori, 2000). For example, for the waveform inversion, we have to average the RFs over backazimuths to suppress the scattering induced by the laterally heterogeneous crustal media(Liu et al., 2014). Such an approach will lead to energy loss of the multiple reflections. On the other hand, the importance of seismic anisotropy has been recognized(Tian et al., 2008; Fang and Wu, 2009; Schulte-Pelkum and Mahan, 2014), and yet the scattering brings many troubles for the anisotropic parameters estimation from RFs. The low-pass filtering can suppress the scattering, to some degree, but it would lose a certain amount of sensitivity to the thickness and depth of the anisotropic layer. This is not what we expect.

For the stratified crustal model, the primary noise in the RFs comes from the scattering induced by the lateral heterogeneity, since the ambient noise has been greatly reduced in the process of the RF isolation. In fact, different from the ambient noise, the lateral scattering is attributed to a self-excited secondary wavefield, which has the similar frequency components to the converted waves in the RFs and has random wave paths or orientations. On the contrary, theoretical tests manifest that the converted phases in the RF, whether in isotropic or anisotropic media, hold spatial continuity. This means that the converted phases and their multiple reflections are well traceable, along with the backazimuth or distance(or station offset)(Nagaya et al., 2008; Fang and Wu, 2009). In this paper, we shall introduce a curvelet denoising method to the RF study based on such differences between the converted phases and scattering noise. This will be significant for the practical RF studies.

The curvelet transform is a novel tool of the multiscale geometric analysis for high-dimensional signal processing(Candès and Donoho, 2000, 2004). In comparison with the classical 2D denoising methods, the curvelet-based method does minimum damage to the signal(Zhang, 2005). Also, the curvelet is superior to the wavelet in the directional feature detection, and therefore, particularly suitable for seismic event identification and seismic data denoising. The unique property of the curvelet transform provides us a new way to suppress the scattering induced by lateral heterogeneity. Currently, the curvelet transform has been used for multiple elimination(Herrmann and Verschuur, 2004; Zhang et al., 2006), noise reduction(Neelamani et al., 2008; Peng et al., 2008), seismic migration imaging(Douma and De Hoop, 2007)and ambient noise tomography(Stehly et al., 2011, 2015). However, it is still difficult to directly apply the curvelet transform to the RF data, since it requires the data satisfy uniform sampling in each dimension. In practice, the RF data are usually nonuniformly sampled or sometimes severely missing in some of the spatial orientations. This is because of the global distribution of seismic zones and also constraints of seismic station deployment(especially temporary stations). Therefore, it is necessary to normalize the RF data along the spatial axis before denoising with the curvelet transform. Chai et al.(2015)interpolated the RF wavefield from the transportable array of the western U.S. and adjacent regions using a weighted average approach. Zhang and Zheng(2015)improved the RF imaging with the reconstructed wavefield from sparsely scattered stations with the cubic-spline interpolation. In the following, we shall introduce the compressive sensing theory(Candès, 2006b; Donoho, 2006)into the RF wavefield reconstruction, so that the curvelet transform can be used in the RF data processing. We shall begin our discussion from a brief introduction to some basic concepts of the curvelet transform and compressive sensing.

2 CURVELET THRESHOLD-DENOISING

Curvelets can be considered as a special type of wavelets with directional property. They provide optimally sparse representations of objects with discontinuities along curves(i.e. ‘edges’)(Candès and Donoho, 2004). Different from the wavelets, the curvelets obey a parabolic scaling relation, i.e. the width of a curvelet approximately equals the square of its length in both spatial and frequency domain. The energy of a curvelet displays an oscillatory behavior along its vertical axis. This makes the curvelets look like local plane waves in spatial plane. In frequency plane, however, the curvelets are compactly supported near parabolic wedges away from the origin. For a better understanding, we visualize four different curvelets with the aid of the Curvelab toolbox(Candès et al., 2006a). Fig. 1 shows spatial and frequency viewpoints of the curvelets, respectively. It should be noted that: 1)for an easy view, the curvelets have been translated in spatial plane and only half of the support regions are shown in frequency plane, since the support regions are symmetrical about the origin; 2)the orthogonal relationship of the curvelets in both spatial and frequency domain relies on the time-frequency scaling property of the Fourier transform.

According to Candès et al.(2006a), the curvelet function can be defined as:

(1)

here, the parameters j, l, k correspond to scale, orientation and position, respectively. The scale takes the value as 2^-j and the orientation satisfies with 0 ≤ θ_{j, l} ≤ 2π, where 0 ≤ l ≤ ${{2}^{\left\lfloor j/2 \right\rfloor }}$-1, $\left\lfloor j/2 \right\rfloor $ denotes the integer part of j/2. The position satisfies x_k^{j, l} = R_{θj, l}^-1(k_l · 2^-j, k₂ · 2^-j/2)where R_θ is a rotation matrix, and k₁, k₂ are the translation parameter. Taking φ_j as a ‘mother’ curvelet, all curvelets at scale 2^-j can be obtained through the rotation and translation of φ_j .

Much like in the orthonormal basis, the curvelet transform is a tight frame, so we have the reconstruction formula:

(2)

where is the inner product of the signal f and the curvelet function φ_{j, l, k}. This means that the curvelet coefficients corresponding to the signal are few and high-valued, while those of the random noise are numerous and low-valued. Thus, the overlap between the signal and noise will be minimal in the curvelet domain, so that in terms of the threshold setting or curvelet shrinkage, the signal can be separated from the noise.

Thus, we suppose that the observed RF matrix d satisfies the following linear equation:

(3)

where the row of the matrix corresponds to the time samples and the column to the spatial parameters(such as backazimuth or distance). m is the noise-free RF matrix and n is Gaussian white noise with a mean value of zero and standard deviation σ. The denoising process by the curvelet thresholding can be described by

(4)

where C denotes the forward curvelet transform(i.e. decomposition matrix)and C^† denotes its inverse(i.e. synthesis matrix). S_w is the threshold function. The tight frame guarantees C^† = C^H, H denotes the conjugate transpose. This means that the inverse curvelet transform is the conjugate transpose of its forward transform. We choose the hard threshold function to set the curvelet coefficients less than the threshold to zero. The threshold will be determined by the noise level.

3 COMPRESSIVE SENSING

Traditional signal sampling has to obey the Nyquist-Shannon theorem. This means that the signal can be accurately reconstructed, only when the sampling rate of a band-limited signal is at least twice its bandwidth. The latest theory, however, known as compressive sensing, goes beyond this limitation, and it has received great attention from the seismic exploration community, since it was proposed(Herrmann and Hennenfent, 2008; Naghizadeh and Sacchi, 2010). Also in seismology, the compressive sensing theory has been applied to the source rupture process of major earthquakes(Yao et al., 2011, 2013)and the reverse time migration of earthquake data(Shang, 2014).

According to the compressive sensing theory, a signal can be reconstructed with high probability from less amount of its samples, as long as the signal is sparse or sparse in some transform domain. This involves designing a measurement matrix(i.e. sensing matrix), which is incoherent with the transform basis and capable of projecting the multidimensional signal into the low-dimensional space, as well as solving an optimization problem(Candès and Wakin, 2008). This implies that the signal can be recovered from samples far below the Nyquist-Shannon frequency. In practice, the sparsity means that the signal is redundant, to a certain degree, i.e. the degree of freedom is far less than the length of the signal. If we choose a proper transform basis, the majority of the coefficients will become zero or small enough to be ignored. Such a signal is called compressible or sparse.

Generally, the sampling below the Nyquist frequency causes signal aliasing, and more specifically, the uniform undersampling will lead to the overlap of the spectrum. The random undersampling is equivalent to random noise in frequency domain. The noise introduced by the undersampling can be separated from the original signal in the transform domain, only if the measurement matrix is incoherent with the transform basis. The random matrix is largely incoherent with any fixed basis, so that they are often chosen as measurement matrix. From this point of view, the signal recovery problem in spatial domain can turn into the denoising problem in frequency domain.

It should be noted that the measurement matrix corresponds to the sampling operation. In general, the RF data satisfy the Nyquist sampling theorem along the time axis, but their samples along the spatial axis are usually nonuniformly sampled. Taking d as the observed RF matrix, and its row denotes the time samples, the column represents the spatial parameters, we can have

(5)

where m is the recovered RF matrix, G is the measurement matrix and n is noise.

Following Naghizadeh and Sacchi(2010), we choose the curvelet as the transform basis, then we can have

(6)

where C^H denotes the conjugate transpose of the curvelet transform, c denotes the unknown curvelet coefficients and W is the mask function. If we choose W coinciding with the threshold, the denoising can be incorporated with the wavefield reconstruction, and then the denoised RFs are obtained. Also, note that these two processes are performed simultaneously.

Let A = GC^HW, the unknown coefficients c are obtained by minimizing the object function:

(7)

where λ is the Lagrange multiplier, which leads to the tradeoff between the misfit and sparsity of the solution. Once c is known, the recovered RF matrix m can be calculated from Eq.(6). In this paper, the optimization of the L1 problem(7)is solved with the spectral gradient-projection method provided by the SPGL1 toolbox(Van den Berg and Friedlander, 2008).

4 NUMERICAL TEST

To test our method, firstly, we consider single station case, in which the dataset is composed of RFs from different backazimuths. We calculate the synthetic RFs of the anisotropic crustal model and take them as the numerical test subject. It should be pointed out that the codes we used do not differ from Levin and Park(1997)in principle.

The crustal model for the synthetic RFs is shown in Table 1. For simplicity, we only consider the case when a plane wave incident on a single-layered anisotropic crust. Its horizontal slowness is 0.0618 s·km^-1, which is equivalent to the epicenter distance of 60°. We add the Gaussian white noise to the synthetic RFs to simulate the scattering noise in data. The signal-to-noise ratio of the RFs is 0 dB.

4.1 Curvelet Denoising

Firstly, we consider the uniform sampling case. Fig. 2 shows the transverse component of the synthetic RFs. It can be seen that the converted phases are well traceable in time-backazimuth domain(Fig. 2a), and the corresponding energy is concentrated in a narrow region near the center of the horizontal axis in frequency domain(Fig. 2b). If the data contain the scattering noise induced by the lateral heterogeneity, the denoising can be easily made by means of the curvelet thresholding shrinkage.

Figures 3a and 3b show that the synthetic RFs(Fig. 2a)with Gaussian white noise and the denoised RFs, respectively. Here, we only keep those curvelet coefficients larger than the threshold and inside the region with subvertical orientations(Stehly et al., 2011). Comparing Fig. 3b with Fig. 2a, we can see that the RFs are accurately recovered.

4.2 Wavefield Reconstruction

As mentioned above, different from the seismic exploration, the backazimuth distribution of the RFs relies on global seismicity. The RF wavefield reconstruction, however, needs a spatially complete wavefield. Fig. 4a shows an uneven RF wavefield, which is randomly sampled from Fig. 2a. We shall reconstruct the randomly sampled RFs, which is shown in Fig. 4b, using the compressive sensing technique. Our results show that the signal-to-noise ratio of the reconstructed wavefield is 36.3 dB, and the amplitudes, polarities as well as phases of all the converted waves have been correctly recovered(Herrmann and Hennenfent, 2008).

5 REAL DATA TEST

For single station observations, the reconstructed missing traces mainly rely on the event distribution. Generally, the observation period of a permanent station is longer than temporary one. Therefore, it is easier to collect data with relatively complete coverage. As the first real data test, we choose the data recorded at Station CTAO of the IRIS Global Seismographic Network. This station is located in northeastern Australia(Lat: -20.0882°N, Lon: 146.2545°E, Network code: IU). A total of 151 teleseismic events(Feb. 1999-Jul. 2015, 30° ≤ △ ≤ 90° and M_b > 5.5)are used, and their RFs are isolated by using the iterative time-domain deconvolution method(Ligorría and Ammon, 1999).

Figure 5 shows the 2D Fourier spectrum of such 151 RFs organized by time and backazimuth(from Station CTAO). Fig. 5 indicates that weak noise is widespread across the entire horizontal axis(corresponding to backazimuth). Comparing Fig. 5b with Fig. 2b, we can see that the lateral scattering energy of the transverse component cannot be ignored. Although the noise is present, it is relatively weak in the radial component. Such noise could be caused by the nonuniform backazimuth sampling.

To evaluate the quality of the wavefield reconstruction for the real data, firstly, we stack 151 RFs with interval of 1° into 123 traces. Then, we randomly choose 35 traces from them as the missing data in a manner similar to the jittered undersampling by Hennenfent and Herrmann(2008). Thus, we can simultaneously reconstruct and denoise the remaining 88 traces to obtain the spatial complete wavefield with interval of 1°. For an easy view, Fig. 6 only shows 20 traces of the observed(red color)and reconstructed(black color)RFs among them.

Figure 6 manifests that the missing data can be well recovered, and through wavefield reconstruction, the anomalies of the amplitudes and phases have been corrected. The variations in the amplitudes and phases of the reconstructed data over backazimuths are more continuous and more reasonable. Therefore, in the case of single station observation, our method can not only efficiently suppress the scattering noise, but also correctly recover the missing data, only if the backazimuthal distribution of events used is ideal enough.

Figure 7 shows the observed and denoised RFs with interval of 10° at Station CTAO, respectively, suggesting that: 1)not only the complete coverage of the RFs can be reconstructed, but also clear and better traceable waveforms can be obtained without loss of the original amplitude energy due to that lateral scattering has been well suppressed; 2)the waveform trace by summing all of the denoised RFs over different backazimuths is the same with the original one; 3)the predicted amplitudes and phases of the converted waves related to the missing traces are coordinated with those around each predicted trace. This implies that we can have more data available for the further study.

Different from permanent stations, temporary observations usually have much shorter period, so it is impractical to obtain a complete coverage over backazimuths. For example, in most part of China, the shortage of seismic events from the northwest is quite common, since the seismicity of the Mediterranean seismic belt is relatively low, so that the events with magnitude larger than 5.0 are rare for one or two years’ observation. For this reason, it is necessary to study the situation, when large backazimuthal section is missing.

In fact, for movable observations, the backazimuthal sampling generally does not meet the requirement of randomness. Before dealing with the real data, we need to test the accuracy of wavefield reconstruction of synthetic RFs with actual backazimuthal coverage. We choose a temporary station S05 of the western Sichuan seismic array, which was operated by the Institute of Geology, China Earthquake Administration during Oct. 2006-Nov. 2009. Fig. 8a shows the geographic location of the station and distribution of the teleseismic events we used. The details regarding the western Sichuan seismic array and Station S05 were described in the paper by Liu et al.(2008).

Figure 8a shows that the seismic events recorded by Station S05 are not only nonuniformly sampled by backazimuth, but also absent in the northwest direction(backazimuth: 200°~300°). This blank area may exceed the length of curvelets, so that the seismic recordings from events in this area are unable to be accurately recovered(Naghizadeh and Sacchi, 2010). Fig. 8b shows the resampled synthetic RFs with the backazimuthal coverage from Station S05. To test the limits of compressive sensing in this case, we reconstruct the resampled data as shown in Fig. 8c. Fig. 8d shows the misfits between the wavefield of Fig. 8c and of Fig. 2a. Fig. 8d shows that the data in the relatively dense areas can be correctly reconstructed, but for the blank area, the amplitude misfits between the original and recovered RFs are 6%~8%. In addition, the signal-to-noise ratio of the reconstructed wavefield drops down to 10.4 dB.

To test whether our denoising method can be used for the observed data, a total of 280 events(M_b > 5.0)recorded at Station S05 are used. The RFs are isolated by using the multi-channel maximal likelihood deconvolution method(Liu and Kind, 2004), which is especially suitable for seismic array data. The backazimuth of the events we used is from 30° to 311° and the epicenter distance is 30° ≤ △ < 90°.

Figure 9 shows the group stacking of the observed and denoised RF data at Station S05 with backazimuthal interval of 10°. Our results suggest that: 1)after denoising, the phase axes are much clearer than the original ones; 2)although the RFs from the blank area(backazimuth: 200°~300°)are unable to be accurately recovered, the denoising causes no problems apart from that area; 3)the summation of all RFs stays almost the same before and after denoising; 4)the phases in the recovered RFs from the blank backazimuths are consistent with those of the adjacent traces. This implies the information carried by the predicted phases is valid.

For the earthquake observations with a 2D seismic profile composed of several stations, the curvelet denoising process is something similar to the case of seismic exploration(Herrmann and Hennenfent, 2008; Peng et al., 2008). The only difference between them is that the layout of movable stations along a 2D profile for earthquake observations generally does not satisfy the requirement of equally-spaced distribution. Different from the single station observations mentioned above, the 2D profile data do not require a series of events over a range of orientations. In order to test our denoising method in this case, we consider 19 stations along the profile of 31°N in the Sichuan moveable seismic array. The station interval is about 5~40 km and the profile length is about 420 km(Liu et al., 2009). Their locations are shown in Fig. 8a. Fig. 10 shows the RFs at these stations before and after the denoising process. The event we used is the Kuril Islands earthquake(M_s6.4; 46.26°N, 153.80°E; 16 km; 19:10:23.2 UTC)on December 7th 2006. The RFs at these stations are isolated by using the multi-channel maximal likelihood deconvolution method(Liu and Kind, 2004). Due to that the recordings from this event at Station S04 and S14 are poor, we replace them with the waveform data from another two events located at roughly the same place as that of the event mentioned above. In particular, for Station S04, we use the Kuril Islands earthquake(M_s5.6; 48.38°N, 154.84°E; 44 km; 10:18:2.3 UTC)on April 9th 2007. And for Station S14, we use the east of the Kuril Islands earthquake(Ms5.9; 45.96°N, 153.63°E; 10 km; 13:50:4.0 UTC)on October 25th 2007.

It can be seen from Fig. 10 that 1)after denoising, the signal-to-noise ratio of each trace is enhanced significantly, and the traceability of the converted phases is also improved; 2)different from the case of Stations S05-S13 located in the Songpan-Garzê unit, the RFs at Stations S15-S19 in the Sichuan basin manifest a simple waveform feature, suggesting that the crust beneath the Sichuan basin is relatively homogeneous and weakly anisotropic. In comparison with that beneath the Songpan-Garzê unit, the strong P_s conversions and their multiple reflections in the radial and transverse components of RFs at Stations S05-S13 may be related to the low velocity zone within the middle crust and complex crustal deformation beneath the Songpan-Garzê unit(Liu et al., 2014); 3)comparing the noise level before the onsets with that in Fig. 9, we can see that the denoising of the RFs from the single station observation is better than that from the seismic profile observation.

6 CONCLUSION AND DISCUSSION

We develop a 2D denoising method for RFs by means of the curvelet transform and compressive sensing. Numerical tests and real data experiments show that our method is efficient in suppressing the lateral scattering in the RF data. The denoised data are more suitable for further study based on the stratified crustal model. The method can be used in datasets both from single station and seismic profile with numerous stations. Also, the improvement of the data quality prevents energy loss in multiple reflections in the stacking process and reduces the data quantity for isotropic structure studies, thus is beneficial for shortening the observation period.

For the RF polarization analysis in the crustal anisotropy study(Qi et al., 2009), the interference from scattering noise is substantially suppressed, which enhances the reliability of the polarization analysis. Also, the recovered data fills the gap left by the missing data, which is helpful for crustal anisotropy study by diagnose of the backazimuthal pattern in the RFs(Tian et al., 2008; Schulte-Pelkum and Mahan, 2014).

Our results show that the denoising of single station observations is better than that of seismic profile observations. This is because the scattering wave from different backazimuths may be more randomly oriented than that from a linear profile. Theoretically, the curvelet denoising method tends to protect the coherent signal and suppress the scattering randomly oriented in spatial domain. Also, the spatial sampling of the profile dataset relies on the deployment of stations, while the backazimuthal sampling of seismic events recorded by a single station are much closer to a random distribution. The randomness of the data sampling ensures the quality of wavefield reconstruction.

Subject to the distribution of the global seismicity, for most part of China, the recovered data from the blank area may not be correct in both amplitude and phase. But this does not prevent us from denoising the existing data. How to overcome the difficulties in accurately recovering data from a blank area is still open.

ACKNOWLEDGMENTS

This research is supported by the National Natural Science Foundation of China (41274060). We thank two anonymous reviewers. Their comments substantially improved the paper. The waveform data were accessed via the IRIS Data Management Center and the State Key Laboratory of Earthquake Dynamics, institute of Geology, China Earthquake Administration. Codes and Toolboxes such as Curvelab, Spot and SPGL1 were used to prepare this paper.

Curvelab:http://www.curvelet.org/

Spot:http://www.cs.ubc.ca/labs/scl/spot/

SPGL1:https://www.math.ucdavis.edu/~mpf/spgl1/