Earthquake Reaearch in China  2019, Vol. 33 Issue (2): 248-264     DOI: 10.19743/j.cnki.0891-4176.201902005
Improving Airgun Signal Detection with Small-Aperture Seismic Array in Yunnan
WANG Weijun1, ZHOU Qingyun2, KOU Huadong1, LIU Guiping1, ZHU Hongbo2, YAN Kun1     
1. Key Laboratory of Earthquake Prediction, Institute of Earthquake Forecasting, CEA, Beijing 100036, China;
2. Yunnan Earthquake Agency, Kunming 650224, China
Abstract: Repeating airgun sources are eco-friendly sources for monitoring the changes in the physical properties of subsurface mediums, but their signals decay quickly and are buried in the noises soon after traveling short distances. Stacking waveforms from different airgun shots recorded by a single seismic station (shot stacking) is the most popular technique to detect weak signals from noisy backgrounds, and has been widely used to process the data of Fixed Airgun Signal Transmission Stations (FASTS) in China. However, shot stacking sacrifices the time resolution in monitoring to recover a qualified airgun signal by stacking many shots at distance stations, and also suffers from persistent local noises. In this paper, we carried out several small-aperture seismic array experiments around the Binchuan FAST Station (BCFASTS) in Yunnan Province, China, and applied the array technique to improve airgun signal detection. The results show that seismic array processing combining with shot stacking can suppress seismic noises more efficiently, and provide better signal-to-noise ratio (SNR) and coherent airgun signals with less airgun shots. This work suggests that the array technique is a feasible and promising tool in FAST to increase the time resolution and reduce noise interference on routine monitoring.
Key words: Airgun Signal Detection     Small-Aperture Seismic Array     Binchuan Fixed Airgun Signal Transmission Stations(BCFASTS)    

INTRODUCTION

Time-lapse monitoring of subsurface medium with seismic methods is very important for many fields from mineral/oil production to natural hazard mitigation, such as earthquakes, landslides and snow slides. Passive repeating sources, e.g. repeating earthquakes (e.g. Poupinet G. et al., 2008; Uchida N. et al., 2019; and references therein) and virtual sources between two seismic stations constructed by ambient noise cross correlation (e.g. Brenguier F. et al., 2008; Lecocq T. et al., 2014; Mordret A. et al., 2016), are generally utilized to eliminate the source uncertainty and improve the accuracy in measuring physical property changes. Man-made active sources, such as airgun (e.g. Wang Baoshan et al., 2016; Chen Yong et al., 2017) and ACROSS (e.g. Kasahara J. et al., 2017), having the advantages of controllable, highly repeatable and short exciting interval, are becoming promising tools for precise and continuous monitoring. Using the active sources, 4-D seismic survey to evaluate 3-D structure deformation is becoming practical in mining and oil industries nowadays (e.g. Lumley D.E., 2004; Du Wenfeng et al., 2017; Kasahara J. et al., 2017).

In contrast, monitoring active faults for earthquake prediction remains challenging. Because stress accumulation is a long-time period in the earthquake process, the physical change would be too subtle to be detected most times. Moreover, obvious change would happen at the critical stage of earthquake nucleation, but might only last for a short period that cannot be resolved by passive source. Hence, it essentially requires precise and continuous monitoring to capture the earthquake process.

To build practicable platforms for monitoring the earthquake processes, several Fixed Airgun Signal Transmission Stations (FASTS) have been setup and running in continental China during the past decade (Chen Yong et al., 2017). Many experiments on airgun source characteristics, data acquisition and analysis have been conducted (e.g., Wang Baoshan et al., 2016; Chen Huifang et al., 2016; Wang Bin et al., 2016; Su Jinbo et al., 2016; Zhang Yuansheng et al., 2016; Yang Wei et al., 2016; Xia Ji et al., 2016; Zhai Qiushi et al., 2016; Wei Congxin et al., 2018; Wang Baoshan et al., 2018, Yang Wei et al., 2018).

To enlarge the monitoring area with limited FASTS, one key task is to detect weak airgun signals from noisy data at long distances. The main frequency contents of a seismic wave emitted from FAST are at a frequency band above 1Hz, which is related to airgun volume and airgun depth underwater (Wang Baoshan et al., 2016; Wei Congxin et al., 2018). Suffering from high attenuation and anthropogenic noise at high frequencies, airgun signals decay quickly over distance and are usually buried in ambient noises. Generally, single shot signal traveling within ten of kilometers can have qualified SNR for monitoring, and beyond that distance, signal stacking is needed.

Shot stacking, which stacks signals from repeating shots at different times, is the main technique to improve the SNR of airgun signals (e.g. Wang Baoshan et al., 2016). In cases of random noise, the SNR gain from directly/linear stacking is typically proportional to the square root of stacking number. For a good case in Yunnan, the signals can be tracked 300km away from BCFASTS by 100 shots stacking (Wang Baoshan et al., 2016), decreasing the time resolution to 100 times of the shot intervals. For many cases, shot stacking is hard to achieve theoretical gain in complicated noise environments, and often outputs incoherent signals.

To get better noise suppression, phased-weighted stacking (PWS, Schimmel M. et al., 1997) bases on signal phase coherency instead of amplitude, providing a more effective strategy to reduce the noise. Schimmel M. et al., (2007) extended the phase coherency estimation on frequency-dependency time windowing (tf-PWS) with the time-consuming S transform method (Stockwell R.G. et al., 1996). Recently, Ventosa S. et al., (2017) used complex frames wavelet theory to build a time-scale phase-weighted stacking method(ts-PWS). It is more efficient and faster to compute, yet still preserves the performance and flexibility of the tf-PWS.

For signals contaminated by persistent noises, Olivier G. et al., (2018) proposed a machine learning framework to separate and partially remove coherent background noises. In general, noises have complex sources and are still difficult to be fully separated from signals recorded by a single station.

Seismic array processing is another promising method to effectively retrieving airgun signals, yet seldom used in FAST applications. An array is a set of sensors arranged in a regular geometric pattern to receive a wave field of passing signal. The technique has been widely used in many military and civilian applications, such as phased array radar, underwater acoustic array, telescope array, microphone array and so on. Seismic arrays were originally built in order to improve the detection of nuclear tests worldwide since 1960s (e.g., YKA or Yellowknife Seismological Array in Canada, 1962; LASA or Large Aperture Seismic Array in USA, 1965; NORSAR or Norwegian Seismic Array in Norway, 1968). An array with M stations can directly improve the SNR similar to M times shot stacking (Rost S. et al., 2002). Waves from different sources can be separated and tracked according to their directions of arrival (DOAs) and propagating velocities in the array. Utilizing the wave field characters, incoherent noises can be suppressed more effectively, and coherent signals can be enhanced to higher SNR. Nowadays, seismic arrays have been widely used in source processes and local to global scale structure imaging (Rost S. et al., 2002; and references therein). In recent years, microtremor small-aperture seismic array is becoming a popular way to investigate near surface soils and their site response to strong motion, as well as volcano monitoring (e.g., Wathelet M. et al., 2008; Bianco F. et al., 2005).

In this paper, we introduced the seismic array technique to process airgun signals. We setup several small-aperture arrays around BCFASTS to verify the possibility and efficiency of array processing to extract airgun signals in noisy backgrounds. We first described the field works and processing methods, then compared the results between shot stacking and array stacking. Next, we discussed the problems and potential improvements in data processing, and finally gave conclusions.

1 DATA AND METHODS 1.1 Seismic Array Installation

Several small-aperture seismic arrays, as shown in Fig. 1, were setup around the Fixed Airgun Signal Transmission Station in Binchuan, Yunan province, China (BCFASTS). BCFASTS was installed in the Dayindian reservoir (Wang Bin et al., 2016), and generally worked at the time from midnight to early morning with low anthropogenic noises. The emission interval was about 25 minutes. The distances from the arrays to BCFASTS are between 9km and 108.5km. Each array has at least 7 stations equipped with the same kind of instruments.

Fig. 1 Map of airgun source (star) and small-aperture seismic arrays (triangles) Array A-D were discussed in this paper. Red lines indicate active faults

We used three kinds of instruments in these arrays, including domestic EPS all-in-one system, Güralp CMG-6TD and Nanomtrics 120P (120P seismometer with Taurus digitizer). All instruments have a flat frequency response in the frequency band of the airgun signal (1-40Hz). It is noted that EPS has the lowest sensitivity (about 250v/m/s, 4 times more less than others) and the highest self-noise. EPS, with integrated three-component sensors, digitizer and battery into a case except the GPS, is the smallest in terms of size and weight, and is the most portable among these instruments. CMG-6TD integrates the sensor and digitizer, and is the second most portable instrument. Instruments similar to EPS have been frequently used in airgun applications, especially in dense seismic arrays (e.g. Bai Zhiming et al., 2016; Wang Baoshan et al., 2018; Tian Xiaofeng et al., 2016).

The instruments were buried in holes with depths of 30-50cm except the GPS antenna, which was exposed in order to receive the GPS signal. The sample rate was set to 100 or 200Hz, and the clock was corrected continuously with GPS.

Local effects, including shallow soil site responses and instrument-soil coupling, may vary from site to site in the array, and may influence the wave coherency across the array. We applied ambient noise horizontal to vertical spectral ratio (NHVSR) to estimate the local effects (e.g. Wang Weijun et al., 2009). Spectral ratios were calculated from different 100s time windows, and averaged to get the final NHVSR. NHVSR curve will be near 1 at different frequencies when the site is on the rock, but will vary with different soil or soil-instrument coupling conditions. Generally, similar NHVSR curves over the same array will benefit the array processing.

1.2 The Array Response Function

The array response function describes array sensitivity and resolution for seismic signals with different frequency contents and slowness, which is defined as (Capon J., 1969):

$ R\left( {{{\vec k}_m}} \right) = \frac{1}{{{N^2}}}{\left| {\sum\nolimits_{i = 1}^N {{e^{ - j\left( {{{}_i}} \right)}}} } \right|^2} $ (1)

where N is sensor number, $\vec{r}_{i}$ is sensor position vector, $\vec{k}_{m}$is wave number vector. The array performance, such as DOA resolution and SNR gain, depends on its geometry and on the wave eld characteristics. To achieve reasonable results in the microtremor small aperture array, the relationship of the minimum and maximum sensor spacing (Dmin and Dmax) and the minimum and maximum wavelengths (λmin and λmax) are generally considered to (e.g. Wathelet M. et al., 2008):

$ \begin{array}{*{20}{l}} {{\rm{ \mathit{ λ} }}\;\max < 3{\rm{D}}\;\max }\\ {{\rm{ \mathit{ λ} }}\;\min > 2{\rm{D}}\;\min } \end{array} $ (2)

If the array geometry fails the first condition, some array benefits will be lost, e.g. unable to discriminate DOA, but the SNR gain is still available. If the second fails, wave field is spatial aliased and signal would be distorted by array stacking, which has to be avoided. In our field installation with limited sensor numbers, Dmin was the main parameter to considering. Then the geometrics were extended according to local topography and available tracks, resulting diverse array geometries in this paper. Generally, S-waves with velocity greater than 1km/s and frequency lower than 8Hz would not be spatially aliased in this paper.

1.3 Shot Stacking and Array Stacking

For signals from different airgun shots received by a single station, linear shot stacking is defined as:

$ u(t) = \frac{1}{N}\sum\limits_{i = 1}^N {{x_i}} (t) $ (3)

where xi(t) is single shot record starting from the time of shot excitation, and N is the shot number. Phase-weighted stacking (PWS) for shots is defined as:

$ u(t) = \frac{1}{N}\sum\limits_{i = 1}^N {{x_i}} (t){\left| {\frac{1}{N}\sum\nolimits_{k = 1}^N {{e^{i{\emptyset _k}t}}} } \right|^v} $ (4)

Where k is the instantaneous phase, and v is a factor to control the coherency sharpness. The phase coherency can be represented in a time-frequency domain (tf-PWS, Schimmel M. et al., 2007) or time-scale domain (ts-PWS, Ventosa S. et al. 2017). The discrepancy between tf- and ts-PWS is ignorable (Ventosa S. et al. 2017), and here ts-PWS was used.

Assuming a plane wave sweeping across the array with vertical incident angle i and back-azimuth θ, the slowness vector across the array is:

$ \begin{array}{l} \vec S\left( {{s_x},{s_y},{s_z}} \right) = {s_{hor}}\left( {\sin \theta ,\cos \theta ,{{\tan }^{ - 1}}i} \right)\\ \;\; = \frac{1}{v}(\sin i\sin \theta ,\sin i\cos \theta ,\cos i) \end{array} $ (5)

where shor and v are the horizontal slowness and media velocity under the array respectively. The slowness vector characterizes the wavefield in the array, and are the key need to precisely estimate in array processing.

Many methods in regards to array processing has been discussed by Rost S. et al., (2002). When both the back azimuth and horizontal slowness are unknown, frequency-wavenumber(F-K) analysis is applied to estimate the values. In the case of known back azimuth or slowness, we can use slant stacking (Vespa process) to scan the slowness/back azimuth. Finally, the array stacking is beamforming with determined back azimuth and horizontal slowness.

In FAST applications, the coordinates of the source and receivers are known, so the back azimuth can be directly calculated and we estimated the horizontal slowness with linear and phased-weighted slant stacking. Similar to eq.(3), linear slant stacking waveform in array center is defined as:

$ {u_s}(t) = \frac{1}{M}\sum\limits_{i = 1}^M {{x_i}} \left( {t - {t_{s,i}}} \right) = \frac{1}{M}\sum\limits_{i = 1}^M {{y_{i,s}}} (t) $ (6)

where ts, i is the relative time delay between array center and station i for horizontal slowness s, and M is the station number in the array, and yi, s(t) is the waveform in station i delay time corrected with slowness. For unknown slowness, different values of slowness are scanned and the optimal one is defined as the one produces maximum energy contend (amplitude) of slant stacking.

Similarly, the phase-weighted stacking for the array is defined as:

$ {u_s}(t) = \frac{1}{M}\sum\limits_{i = 1}^M {{y_{i,s}}} (t){\left| {\frac{1}{M}\sum\nolimits_{k = 1}^M {{e^{i{\emptyset _{{k^t}}}}}} } \right|^v} $ (7)

where k is the instantaneous phase for yk, s(t), and yk, s(t) is time corrected waveform as eq.(6), and v is a factor to control the coherency sharpness.

Data processing in this paper has 5 steps: (1) cutting 100s waveforms starting from shot time; (2) shot stacking with specific shot number for each component in each station; (3) estimating horizontal slowness with known back-azimuth; (4) array beamforming a shot records with fixed back-azimuth and horizontal slowness; (5) shot stacking over outputs from step 4. The last step combines both the array and shot stacking to achieve airgun signals with higher quality.

Both linear and phase-weighted methods were applied to shot and array stacking, and the phase-weighted parameter v is set to 2 and 1 for shot and array stacking, respectively. In step 3, P-wave and S-wave horizontal slowness was searched in the range of 0-50s/deg with 0.2s/deg step, using time window 2s around the P arrival in the vertical component and 2s-5s around the S arrival in N-S component, separately. In step 4, P-wave and S-wave were beamformed in the vertical component with P slowness and in N-S component with S slowness, respectively. As shown in the following results, amplitude variations caused by site conditions would be obvious across the array. To reduce the influence in array processing, waveforms were normalized with the largest amplitude in band pass filtered S-wave window in each component.

We output time serial waveforms by moving stacking with a specific stacking window and sliding step. E.g., for 26 shot records, moving stacking with 10-shot window and 1-shot step will produce 17 stacked waveforms. The stacking window was defined by trial to get acceptable SNR. For comparison's sake, all waveforms were bandpass filtered to 2-8Hz, which is the frequency band with major energy in airgun signal.

2 RESULTS

We only showed results from Array A-D in Fig. 1 in this paper. The Array A and B are located 9.0km and 82.7km SW away from BCFASTS, and both have 10 stations equipped with Nanometrics 120P. Array C and D are 108.5km and 99.5km north away from BCFASTS. Array C has 9 stations equipped with Gurulap 6TD but only 7 stations had data recorded. Array D has 12 stations equipped with EPS. Array A recorded 5 airgun excited during daytime, and Array B recorded 60 shots over 3 nights. Array C and D recorded 26 shots in the same night.

Four arrays have diverse geometries with different aperture and sensor spacing (top row in Fig. 2). Their corresponding array response functions are shown in the bottom row in Fig. 2. The more energy the response function focuses on the center, the higher the resolution and ability to discriminate and suppress different slowness across the array. According to the response functions, Array B and D have better resolution to distinguish slowness from different directions, while Array A and C have higher sensitivity form waves propagating along the array strike.

Fig. 2 Top row: The geometries of Array A-D Red triangles are seismic stations, and crosses are array centers. Bottom row: Corresponding seismic array transfer functions for Array A-D. Hot to cold colors represent high to low energy

The NHVSR curves for each station in 4 arrays are shown in Fig. 3, obviously not consistent across the array, especially for Array A and B which have large fluctuations. The extremely large and small NHVSRs in Array A and B, respectively, are probably caused by sensor-soil coupling problems. Array C was installed in a mountainous area and has the most uniform NHVSRs across the array. Nevertheless, we still found largely fluctuating waveform amplitudes across the same array for some shots, suggesting the site differences in the array. As completely removing the local effects is still impractical, the normalizing in array processing can partly reduce the contribution from abnormal amplitude and stabilize the beamforming output.

Fig. 3 Curves of NHRSV for each station in Array A-D

Fig. 4 shows an example of the P- and S-wave F-K analysis for one shot recorded in Array A. The back-azimuth angle from F-K deviates the true angle by about 8 degrees, and the horizontal slownesses have several s/deg bias comparing with the results of Fig. 5(a), suggesting low resolution in F-K analysis that are possibly due to limited receiver numbers and array geometry (Fig. 4).

Fig. 4 F-K analysis for P (a) and S waves (b) in Array A Hot colors represent high energy, and black lines show the slowness vectors

Fig. 5 Slowness searched by array slant stacking with fixed back-azimuth for Array A-D The stackings are based on 1-shot, 10-, 20- and 20-shots for Array A-D respectively. P and S in legends represent P and S waves, and dls for linear and pws for phase-weighted stacking separately

Fig. 5 shows the horizontal slowness of the P- and S-wave searched by array slant stacking. Benefit from high SNR of original signals, stable slowness can be estimated from each shot in Array A, and are also comparable with the F-K analysis in Fig. 4. The differences between two stacking methods for S-wave slowness are probably caused by contaminated S-wave by P coda. This would disturb the wave phases which are sensitive for phase-weighted stacking. We selected the linear result as the true slowness.

For Array B-D, slowness cannot be estimated stably with single shot, so we firstly N shot stacking for each station in array separately, then run slant stack to retrieve slowness. For Array B, N=10 is enough to get stable slowness over different stacking windows. For Array C-D, N=20 and larger values were tested and showed similar results. we chose the values appeared most frequently for next processing.

Fig. 6 shows the airgun waveforms in Array A. The waveforms from single shot have high SNR and clear P-wave arrivals, but are different in shape (Fig. 6(a)). It is probably caused by local effects and local noises, as stations are space clustered with almost same ray paths. 5-shot waveforms in single station show high coherency (Fig. 6(b)), indicating the airgun sources are highly repeated. The waveforms stacked from the array preserve the same feature of high coherency over shots (Fig. 6(c) and 6(d)), and are comparable with the original waveforms at the station near the array center (Fig. 6(b)). Waveforms stacked by linear and phase-weighted array methods (Fig. 6(c) and 6(d)) are similar to each other, featured by attenuated P coda and more obvious later phase signals. Phase-weighted array stacking suppresses more unrelated signals than linear method. It should be noted that the S waves in vertical components would be distorted because the component is beamformed with only P-wave slowness, and it is vice versa for P waves in horizontal components.

Fig. 6 Airgun signals in Array A (a) Vertical component waveforms of one shot recorded in 10 stations. (b) Vertical(Z) and north-south(N) component waveforms of 5 shots recorded in station near the array center. (c) Linear and (d) phase-weighted array stacking waveforms for 5 shots in Z and N components. All waveforms are 2-8Hz bandpass filtered and vertical dashed lines mark a time close to P-wave arrivals. Horizontal dash lines separate N and Z components

For Array B, airgun signals were buried in noises and are hard to clearly identify the P- and S-wave from a single shot without stacking. With 5 shots stacking, both waves in a station can be identified more clearly (Fig. 7(a)-(b), only partial results are shown in a-better waveform comparison), but the waveform coherency over time is quite poor. Local noises can be still dominant after stacked (e.g. the top 4 traces in Fig. 7(a)-(b)). Array stacking with over 5 shots produced qualified waveforms for both P and S waves with higher SNR and higher coherency (Fig. 7(c)-(f)). Local noises are dramatically suppressed, especially with phase-weighted array stacking (Fig. 7(e)-(f)). Coherent later phases after P and S arrivals appear to fully explore the array potentials, coherently.

Fig. 7 Airgun signals in Array B Comparisons of linear shot stacking(a-b), linear(c-d) and phase-weighted (PWS, e-f) array stacking. Left panels are for vertical component(Z) and right ones for north-south component (N). 10 shots data were used in each stacking and only 10 of earlier results (skipped every 2 shots) are displayed for clear comparison. Waveforms are 2-8Hz bandpass filtered and vertical dashed lines mark the onset of P-wave arrivals

To further show the power of array processing, waveforms of array stacking with single shot are shown in Fig. 8. Comparing to the results of array stacking with 5 shots (Fig. 7(c)-(f)), these results have relative lower SNR and coherency, especially for P waves. But they may be good enough for time-lapse processing, meaning the time resolution in monitoring can be increased to one shot intervals.

Fig. 8 Results of vertical (left column) and north-earth (right column) components by phase-weighted array stacking with only one shot There are 60 shots totally in Array B. Waveforms are 2-8Hz bandpass filtered and vertical dashed lines mark a time close to P-wave arrivals

Fig. 9 shows stacking results for Array C. 10-shot moving stacking was applied. S-waves from Shot stacking are prominent and coherent in different traces, but their onsets are still blurred by background noises (Fig. 9(b)). Shot stacking P-waves are vague and prone to be misjudged by predominant wave packets before them (Fig. 9(a)). Array stacking with 10 shots apparently suppresses most noises, significantly improves the signal qualities of both P and S waves (Fig. 9(c)-(d)).

Fig. 9 Airgun signals in Array C. Comparisons of phase-weighted (PWS) shot stacking with 10 shots(a)-(b) and phase-weighted (PWS, (c)-(d)) array stacking with 10 shots Left column is vertical component(Z) and right is North-south component (N). Waveforms are 2-8Hz bandpass filtered and vertical dashed lines mark a time close to P-wave arrivals

Stacking the data of Array D is a special case. Array D was setup in weathered red soil similar to Array B, and has a quieter observation environment than Array C which was located less than 1km near a busy highway. The instruments used in this array are EPS. Unlike Array A-C, the P and S waves in Array D are hard to be identified by shot stacking in all stations (e.g. 10 shots stacking in Fig. 10(a)-(b)), or by array stacking with one shot. These explained unstable slowness searched by slant stacking in Fig. 5(d). After array stacking with about 10 shots, the airgun signals emerge with acceptable SNR (Fig. 10(c)-(d)).

Fig. 10 Airgun signals in Array D. Comparisons of phase-weighted (PWS) shot stacking(a)-(b) and phase-weighted (PWS, (c)-(d)) array stacking Left column is for vertical component (Z) and right is for north-south component (N). 5 shots moving stacking were used. Waveforms are 2-8Hz bandpass filtered and vertical dashed lines mark a time close to P-wave arrivals
3 DISCUSSIONS AND CONCLUSIONS

We tested three types of instruments in this paper. Proper instruments are the prerequisite for array application: convenience saves time and money during the deployment, and data quality guarantees reliable signal detection. Our results show that the arrays equipped with these instruments can work and detect airgun signal about 100km, and would have the potential for further distance.

Among these instruments, The EPS system used in Array D has the greatest convenience with its portability in array application, yet we found some potential problems that still need to be further clarified. Firstly, we notice unsuccessful shot stacking in this array with 10 shots (Fig. 10(a)-(b)) or all shots, although this array is closer to airgun source and has a better observation environment in comparison to Array C. We suspect that this may relate to the instrument's high self-noise. When airgun signal attenuates to the level lower or comparable with instrument's self-noise, signal detection would be impossible. High self-noise means the instrument can be used in FAST with limited range close to the source.

Secondly, something abnormal may also relate to the instrument. Short-term P waves delays in the range of about 0.2-0.5s are observed in the middle vertical traces in Array D (Fig. 10(c)), but such changes can be seen neither on the S-waves on the NS component (Fig. 10(d)) nor on the P and S waves in Array C (Fig. 9), excluding possibilities of subsurface media change. Another possibility may be a timing issue that exists within the instrument, which happened occasionally in similar instruments with improper firmware in previous experiences.

Moreover, obvious vP/vS differences are found on similar ray paths from the source to Arrays C and D. The P/S arrivals are about 21s/34s and 20s/35.5s for Array D and C, corresponding to vP/vS about 1.61 and 1.77, respectively (Fig. 9 and Fig. 10). So far, we could not discriminate whether they are caused by abnormal near surface structure under the arrays, or from instrument timing problem in one of the arrays. According to above situations, we suggest that stricter tests should be done before deployment of the EPS and similar instruments in FAST applications.

Nevertheless, we successfully applied the array technique to improve the airgun signal detection for 4 arrays around BCFASTS in the Yunnan area, proved its advantages compared with traditional shot stacking. Apparently, array sacrifices sensors to gain the SNR. For arrays with 7, 10 and 12 sensors in this paper, they can theoretically get about 2.65, 3.26 and 3.46 times SNR gain from a single shot, respectively, hence can reduce shot number and increase the time resolution in monitoring.

But the most important advantage of the array is its ability to suppress noise by wavefield characters. The suppression brings not only the SNR gain, but also the waveform coherency, which should be a key criteria for good waveform in FAST applications. We have shown that nonrandom noises are hard to suppress by shot stacking as shown in array C (Fig. 9(a)-(b)). We also showed that incompletely suppressed noises can lead to low coherency signals although with high SNR (e.g. Fig. 7(a)-(b)), which are also frequently seen in previous studies. Because precise estimation of underground velocity or attenuation changes is highly dependent on tiny disturbance of time or amplitude on repeating waveforms, coherency acts as much important role in monitoring. At this point, array stacking, separating waves from difference sources and enhancing a particular incident wave by utilizing the wavefield characters, can suppress noise more completely and work more effectively for complicated noise sources. Thus, array stacking, especially in combination with shot stacking, could potentially yield higher quality waveforms than purely shot stacking even when they have same theoretical SNR gain.

There are many spaces that can be improved by the array technique. The originations of the phases appeared after the P- or S-phase are still undear, and may be solved by array analysis.Furthermore, airgun signals can be better separated from other signals by fully utillize the vector wavefield with the three components(3C) technique. Overall, seismic array technique is a feasible and promising tool in weak signal detection, and can advance the FAST applications in the future.

ACKNOWLEDGEMENT

We gratefully thank Wang Baoshan, Yang Wei, Yang Jun and Li Xiaobin for their support with the instruments and airgun excitations. We are very thankful to the three anonymous reviewers for their constructive comments and text improvement on the original manuscript.

REFERENCES
Bai Zhiming, Wu Qingju, Xu Tao, Wang Xiao. Basic features of the crustal structure in the lower yangtze and its neighboring area in the Chinese mainland:review of deep seismic sounding research[J]. Earthquake Research in China, 2016, 30(3): 298-315.
Bianco F., Cusano P., Petrosino S., Castellano M., Buonocunto C., Capello M., Del Pezzo E. Small-aperture array for seismic monitoring of Mt. Vesuvius[J]. Seismological Research Letters, 2005, 76(3): 344-355. DOI:10.1785/gssrl.76.3.344
Brenguier F., Campillo M., Hadziioannou C., Shapiro N.M., Nadeau R.M., Larose E. Postseismic relaxation along the San Andreas Fault at Parkfield from continuous seismological observations[J]. Science, 2008, 321(5895): 1478-1481. DOI:10.1126/science.1160943
Brenguier F., Rivet D., Obermann A., Nakata N., Boué P., Lecocq T., Campillo M., Shapiro N. 4-D noise-based seismology at volcanoes:Ongoing efforts and perspectives[J]. Journal of Volcanology and Geothermal Research, 2016, 321: 182-195. DOI:10.1016/j.jvolgeores.2016.04.036
Capon J. High-resolution frequency-wavenumber spectrum analysis[J]. Proceedings of the IEEE, 1969, 57(8): 1408-1418. DOI:10.1109/PROC.1969.7278
Chen Huifang, Lin Binhua, Jin Xing, Wu Lihua, Cai Huiteng. An experimental study on optimization of large-volume airgun source excitation conditions in a reservoir[J]. Earthquake Research in China, 2016, 30(3): 355-363.
Chen Yong, Wang Baoshan, Yao Huajian. Seismic airgun exploration of continental crust structures[J]. Science China Earth Sciences, 2017, 60(10): 1739-1751. DOI:10.1007/s11430-016-9096-6
Coutant O., Doré F., Brenguier F., Fels J.F., Brunel D., Judenherc S., Dietrich M. The high-resolution imaging (HRI) portable array:a seismic (and internet) network dedicated to kilometric-scale seismic imaging[J]. Seismological Research Letters, 2008, 79(1): 47-54. DOI:10.1785/gssrl.79.1.47
Du Wenfeng, Peng Suping. An application of 4D seismic monitoring technique to modern coal mining[J]. Geophysical Prospecting, 2017, 65(3): 823-839. DOI:10.1111/gpr.2017.65.issue-3
Kasahara J., Hasada Y. Time Lapse Approach to Monitoring Oil, Gas, and CO2 Storage by Seismic Methods[M]. Amsterdam: Gulf Professional Publishing, 2017.
Lecocq T., Caudron C., Brenguier F. MSNoise, a python package for monitoring seismic velocity changes using ambient seismic noise[J]. Seismological Research Letters, 2014, 85(3): 715-726. DOI:10.1785/0220130073
Lumley D.E. Business and technology challenges for 4D seismic reservoir monitoring[J]. The Leading Edge, 2004, 23(11): 1166-1168. DOI:10.1190/1.1825942
Mordret A., Mikesell T.D., Harig C., Lipovsky B.P., Prieto G.A. Monitoring southwest Greenland's ice sheet melt with ambient seismic noise[J]. Science Advances, 2016, 2(5): e1501538. DOI:10.1126/sciadv.1501538
Olivier G., Chaput J., Borchers B. Using supervised machine learning to improve active source signal retrieval[J]. Seismological Research Letters, 2018, 89(3): 1023-1029. DOI:10.1785/0220170239
Poupinet G., Got J.L., Brenguier F. Chapter 14 monitoring temporal variations of physical properties in the crust by cross-correlating the waveforms of seismic doublets[J]. Advances in Geophysics, 2008, 50: 373-399. DOI:10.1016/S0065-2687(08)00014-9
Rost S., Thomas C. Array seismology:methods and applications[J]. Reviews of Geophysics, 2002, 40(3): 1008.
Schimmel M., Paulssen H. Noise reduction and detection of weak, coherent signals through phase-weighted stacks[J]. Geophysical Journal International, 1997, 130(2): 497-505. DOI:10.1111/gji.1997.130.issue-2
Schimmel M., Gallart J. Frequency-dependent phase coherence for noise suppression in seismic array data:frequency-dependent phase coherence[J]. Journal of Geophysical Research:Solid Earth, 2007, 112(B4).
She Yuyang, Yao Huajian, Zhai Qiushi, Wang Fuyun, Tian Xiaofeng. Shallow crustal structure of the middle-Lower Yangtze River Region in Eastern China from surface-wave tomography of a large volume airgun-shot experiment[J]. Seismological Research Letters, 2018, 89(3): 1003-1013. DOI:10.1785/0220170232
Stockwell R. G., Mansinha L., Lowe R. P. Localization of the complex spectrum:the S transform[J]. IEEE Transactions on Signal Processing, 1996, 44(4): 998-1001. DOI:10.1109/78.492555
Su Jinbo, Wang Qiong, Wang Haitao, Shi Yongjun, Chen Xiangjun, Feng Lei, Chen Hao. Research on large volume airgun source signal reception in Hutubi, Xinjiang using seismic network data[J]. Earthquake Research in China, 2016, 30(3): 326-332.
Tian Xiaofeng, Wang Fuyun, Liu Baofeng, Yang Zhuoxin, Zheng Chenglong, Gao Zhanyong. Upper crustal velocity structure along the Yangtze River from Ma'anshan to Anqing[J]. Earthquake Research in China, 2016, 30(3): 430-438.
Uchida N., Bürgmann R. Repeating earthquakes[J]. Annual Review of Earth and Planetary Sciences, 2019, 47. DOI:10.1146/annurev-earth-053018-060119
Ventosa S., Schimmel M., Stutzmann E. Extracting surface waves, hum and normal modes:time-scale phase-weighted stack and beyond[J]. Geophysical Journal International, 2017, 211(1): 30-44. DOI:10.1093/gji/ggx284
Wang Baoshan, Ge Hongkui, Wang Bin, Wang Haitao, Zhang Yuansheng, Cai Huiteng, Chen Yong. Practices and advances in exploring the subsurface structure and its temporal evolution with repeatable artificial sources[J]. Earthquake Research in China, 2016, 30(3): 284-297.
Wang Baoshan, Tian Xiaofeng, Zhang Yunpeng, Li Yulan, Yang Wei, Zhang Bo, Wang Weitao, Yang Jun, Li Xiaobin. Seismic signature of an untuned large-Volume airgun array fired in a water reservoir[J]. Seismological Research Letters, 2018, 89(3): 983-991. DOI:10.1785/0220180007
Wang Bin, Li Xiaobin, Liu Zifeng, Yang Jun, Yang Runhai, Wang Baoshan, Yang Wei. The source and observation system of Binchuan Earthquake signal transmitting seismic station and its preliminary observation results[J]. Earthquake Research in China, 2016, 30(3): 316-325.
Wang Weijun, Liu Lanbo, Chen Qifu, Zhang Jie. Applications of microtremor H/V spectral ratio and array techniques in assessing the site effect and near surface velocity structure[J]. Chinese Journal of Geophysics, 2009, 52(6): 1515-1525 (in Chinese with English abstract).
Wathelet M., Jongmans D., Ohrnberger M., Bonnefoy-Claudet S. Array performances for ambient vibrations on a shallow structure and consequences over Vs inversion[J]. Journal of Seismology, 2008, 12(1): 1-19. DOI:10.1007/s10950-007-9067-x
Wei Congxin, Qin Manzhong, Zhang Yuansheng, Zou Rui, Wang Lei, Guo Xiao, Liu Xuzhou, Wang Yahong, Sun Dianfeng. Airgun excitation experiments at different placement depths in the Qilian Mountain of Gansu Province, China[J]. Seismological Research Letters, 2018, 89(3): 974-982. DOI:10.1785/0220170253
Xia Ji, Jin Xing, Cai Huiteng, Xu Jiajun. The time-frequency characteristic of a large volume airgun source wavelet and its influencing factors[J]. Earthquake Research in China, 2016, 30(3): 364-379.
Yang Wei, Wang Baoshan, Liu Zhengyi, Yang Jun, Li Xiaobin, Chen Yong. Study on the source characteristics of downhole airgun with different excitation environments[J]. Earthquake Research in China, 2016, 30(3): 342-354.
Yang Wei, Wang Baoshan, Yuan Songyong, Ge Hongkui. Temporal variation of seismic-wave velocity associated with groundwater level observed by a downhole airgun near the Xiaojiang Fault Zone[J]. Seismological Research Letters, 2018, 89(3): 1014-1022. DOI:10.1785/0220170282
Zhai Qiushi, Yao Huajian, Wang Baoshan. Study on the deconvolution method and processing flow of airgun source data[J]. Earthquake Research in China, 2016, 30(3): 394-404.
Zhang Yuansheng, Guo Xiao, Qin Manzhong, Liu Xuzhou, Wei Congxin, Shen Xuzhang, Yan Wenhua. The construction of active source repeated monitoring in the Qilian Mountain, Gansu Province[J]. Earthquake Research in China, 2016, 30(3): 333-341.