Quick Phase Identification for Dense Seismic Array with Aidfrom Long Term Phase Records of Co-located Sparse Permanent Stations1

CHEN Yini ,LI Jun* ,WANG Zhenjie ,ZHAO Mengqi ,YU Tiehong ,LI Danning ,YU Junyi and WANG Weitao

1) Zhejiang Earthquake Agency,Hangzhou 310013,China

2) Key Laboratory of Earthquake Source Physics,Institute of Geophysics,China Earthquake Administration,Beijing 100086,China

3) Yunnan Earthquake Agency,Kunming 650224,China

The phase identification and travel time picking are critical for seismic tomography,yet it will be challenging when the numbers of stations and earthquakes are huge.We here present a method to quickly obtain P and S travel times of pre-determined earthquakes from mobile dense array with the aid from long term phase records from co-located permanent stations.The records for 1 768 M≥2.0 events from 2011 to 2013 recorded by 350 ChinArray stations deployed in Yunnan Province are processed with an improved AR-AIC method utilizing cumulative envelope and rectilinearity.The reference arrivals are predicted based on phase records from 88 permanent stations with similar spatial coverage,which are further refined with AR-AIC.Totally,718 573 P picks and 512 035 S picks are obtained from mobile stations,which are 28 and 22 times of those from permanent stations,respectively.By comparing the automatic picks with manual picks from 88 permanent stations,for M≥3.0 events,81.5%of the P-pick errors are smaller than 0.5 second and 70.5%of S-pick errors are smaller than 1 second.For events with a lower magnitude,76.5%P-pick errors fall into 0.5 second and 69.5%S-pick errors are smaller than 1 second.Moreover,the Pn and Sn phases are easily discriminated from directly P/S,indicating the necessity of combining traditional auto picking and integrating machine learning method.

Key words:Phase picking;Travel time;Dense array;Spatial overlap

INTRODUCTION

Body wave tomography based on seismic wave travel time is an important tool for studying the interior structure of the earth.Teleseismic tomography can provide deep information within a depth of several hundred kilometers,while imaging studies based on regional and local seismic data can study the fine scale of the crust and mantle (Zhao Dapeng,2015).

The density and spatial distribution of body wave rays from an earthquake epicenter to a station are the key factors that determine the quality of body wave tomography.In order to obtain more detailed underground structure,researchers often choose long-term records to increase the number of earthquakes,as well as installing dense observation stations so that the ray coverage will be improved for better imaging (Huang Jinli et al.,2006;Zigone D.et al.,2019).In recent years,with the deployment of dense arrays such as USArray,ChinArray and other large-scale broadband seismic arrays,the observation basis of body wave imaging has been greatly improved.Furthermore,some small-scale and ultra-dense observations also provide excellent data for shallow structure imaging (Li Yonghua et al.,2014;Lin Fanchi et al.,2013;Liu Zhen et al.,2017).

The increasing number of array observations also poses challenges to the traditional manual travel time picking work.The program ChinArray Phase I has been conducted in Yunnan Province and adjacent regions by China Earthquake Administration since the September 2010.From 2011 to 2013,350 broadband portable stations have been deployed in the southern section of the North-South Seismic Belt (Ding Zhifeng et al.,2013).Numerous earthquakes have occurred in Yunnan,yet the identification of seismic phases on portable stations could not be conducted due to the limitation of manual picking efficiency with only conventional manual quick reports and cataloging processes in the permanent network of Yunnan.In the previous studies of body wave tomography based on ChinArray,teleseismic events are used to image the deep structure of the earth,because of the limited quantity of teleseismic body wave phases and the low workload of manual picking (Zhang Fengxue et al.,2018;Huang Zhouchuan et al.,2015).As for regional studies on ChinArray,major events with large magnitudes are analyzed.For example,Gao Jiayi et al.(2016) manually pick the travel times of 45M≥3.0 earthquakes in the region,and construct a one-dimensional velocity model in Yunnan.For a large number of smallmagnitude earthquakes,limited by the efficiency of manual recognition,the automated processing methods still need to be developed to achieve rapid phase recognition and travel time extraction.

Automatic identification of seismic phases has always been an important part of seismology (Allen R.,1982).As a result,the method of automatic identification has developed rapidly with the increase of seismic data volume.The classical STA/LTA (Short Term Averaging/Long Term Averaging) algorithm for automatic detection of microseismic signals based on waveform features(Allen R.,1982;Baer M.et al.,1987;Saragiotis C.D.et al.,2002;Baillard C.et al.,2014),AR-AIC (Leonard M.et al.,1999;Sleeman R.et al.,1999;Zhang Haijiang et al.,2003) have been widely used.In detail,the STA/LTA model is sensitive to variations in amplitude.AR-AIC is not only based on changes in amplitude,but also related to changes in frequency and phase.Normally we use STA/LTA to detect events,which can remove interference from pulse signals,and apply AR-AIC to mark the arrivals of specific phases (Akazawa T.,2004).However,when it comes to dealing with complex observation waveforms,the ability to pick up S-wave phases needs to be improved.By polarizing the three-component data and suppressing P code wave interference before S-wave arrives,the success rate of S detection can be effectively improved (Jurkevics A.,1988;Ross Z.E.et al.,2014).In recent years,deep learning techniques have also been applied to seismic phase recognition research to identify P and S phases automatically,including support vector machines,convolutional neural networks,etc (Yu Ziye et al.,2018;Ross Z.E.et al.,2018;Zhu Lijun et al.,2019).

Seismic phase recognition is mainly used in two aspects.Firstly,it is effective to identify seismic events from continuous data,distinguish seismic phases and pick up travel times,and then realize phase association and location monitoring through joint analysis of multiple stations (Zhang Miao et al.,2019).Secondly,it is used to pick up the travel time of different phases in the seismic record with known original time and location of the seismic event,including the fast picking of the travel time from active source signals (Xu Zhen et al.,2019),and the travel time picking of seismic events at different stations.Li Zefeng and Peng Zhigang (2016) use the existing catalog and initial one-dimensional velocity model to estimate the arrival time of events from each station,and then use a signal-to-noise-ratio-based characteristic function to pick up,updating the velocity model.After several times of iteration,more accurate arrival time and 1-D model can be obtained at the same time.

There are 127 seismic stations in Yunnan and its adjacent regions where the ChinArray stations are deployed (Fig.1).These stations have been running steadily since 2008,and seismic phase reports and earthquake catalogs are regularly produced.Due to the overlap of the space area between the permanent and ChinArray stations,the underground space sample is relatively consistent,as well as the travel time of seismic events.Therefore,using the existing observation reports and earthquake catalogs of fixed stations,the corresponding seismic phases on portable stations of ChinArray can be quickly identified.Based on this,an improved AR-AIC method is applied to quickly analyze the travel time of the 2006M≥2.0 earthquakes of on ChinArray from 2011 to 2013.

1 DATA AND METHOD

The target area,Yunnan and its adjacent regions (96°-108°E,20.5°-30°N),includes 127 permanent seismic stations.From March 2011 to November 2013,the China Earthquake Administration has deployed 350 mobile broadband stations in the southern section of the North-South Seismic Belt as China Earthquake Scientific Exploration Array(ChinArray) Phase I,with an aperture of about 1 000 km and an average station distance of 35 km (Ding Zhifeng et al.,2013;Wang Weitao et al.,2018).The distribution of stations is shown in Fig.1,in which all the stations are recorded with three components with a sampling rate of 100 Hz.

In the area of 20.5°-30.5°N and 96°-108°E,we select 1 796 events ofML≥2.0,among which 326 events areML≥3.0.The distribution of seismic magnitude against time is shown in Fig.2.The catalog information is provided by China Earthquake Networks Center(CENC).

Figure 1 Station configuration in Yunnan area Solid triangles represent ChinArray distribution,hollow triangles indicate permanent stations;the circles mark the local seismic events occurred from 2011 to 2013

Figure 2 Time-Magnitude distribution between March 2011 and December 2013

The theoretical travel times of the corresponding seismic phases are calculated using TauP (Crotwell HP et al.,1999) with the 1-D velocity model from Yunnan Station Network (Fig.3).The results are roughly consistent with the manually recognized travel times.The crustal thickness in Yunnan is relatively thick with large lateral variation.We select the phase of the events with an epicenter larger than 200 km to analyze the Pn phases.

Figure 3 (a) 1-D velocity model from bulletin in Yunnan;(b) theoretical travel times corresponding to the 1-D model Black dots indicate the travel times from manual picking of events above ML2.0 between the year 2011 and 2013

According to the theoretical calculation,it is possible to estimate the travel time of various seismic phases of the local earthquakes,and conduct seismic phase identification.We firstly take event locations from CENC bulletin and a preliminary velocity model of Yunnan and its adjacent regions as input.Given the location of an event and the preliminary velocity model,the phase arrivals are predicted at several stations.We then apply an improved AR-AIC detector based on the characteristic function to search the actual arrival within a window around the predicted arrival.To determine the window for searching the real phase arrival,we assume that the actual velocity deviates within a percentage compared with the velocity model

in whichvis a given velocity model and Δvmaxis the maximal allowed difference between the true velocity structure and the given velocity model.By propagating of the velocity perturbation to arrival-time perturbation,the time window becomes

P phase is picked on the vertical component.1) We first apply bandpass filter from 1 to 10 Hz in the time window;2) and then calculate the cumulative envelope function of the filtered record (Akazawa T.,2004):

in whichav1(t) is the filtered record from step 1),and av2(t) is the achieved cumulative envelope function.3) Finally,we apply improved AR-AIC detection to the cumulative envelope (Fig.4).

After P phase arrival is obtained,the corresponding S phase arrival can be estimated based on the theoretical P-S time difference.And then the picking time window for S phase can be determined by Equa.(2) as well.For S phase identification,an S-wave polarization filter is applied to the horizontal components to suppress the energy of P coda as much as possible.After that,cumulative envelope function calculation and AR-AIC picker is applied as the same procedure in P phase identification (Fig.4).The design of the S filter requires polarization analysis of the three-component waveform to calculate rectilinearityrand apparent vertical incidence angleφ.The details are introduced as follows (Jurkevics A.,1988;Ross Z.E.et al.,2014):

Firstly,take a sliding window of 3s for the three-component waveform,and calculate the covariance matrixσof the waveform in each window

and then calculate the eigenvaluesλ1,λ2,λ3(λ1≥λ2≥λ3) and eigenvectorsu=(u1,u2,u3) of this covariance matrix.Then,the linearity r and the apparent normal incidence angle φ are calculated as

The S-wave polarization filter is defined as

Figure 4 Picking procedure

2 RESULTS

2.1 The Number of Picked Phases

The statistical result shows that based on the improved AR-AIC method,P-wave and S-wave can be well identified.Fig.5 shows the ray coverage distributions without and with ChinArray stations of aML3.3 earthquake at 19 ∶54 on November 9th,2012,and its waveform plotting with automatic P and S pickers marked.It can be seen from Fig.5 that the automatically picked seismic phases are in good agreement with the actual seismic phases.

Figure 5 The picking result and ray path coverage map of ML3.3 event The black triangles mark permanent stations,and the blue triangles mark the mobile stations

The density of ChinArray stations is much higher than that of the permanent stations.Through automatic phase picking,the number of regional phases and seismic rays can be greatly increased.Fig.6 is aML2.1 event compared with the event showed in Fig.5.It also shows the comparison between the number of phases picked up at fixed stations and the number of phases picked up by the ChinArray stations and their ray path coverages respectively.The number of phases corresponding to the two earthquakes is shown in Table 1.The original seismic stations are relatively sparse,thus the number of phases P recorded for a single event is between 10 and 40,and the number of S phases is between 10 and 30.After adding dense array data,the number of P phases that can be effectively picked up is between 200 and 400,and the number of S phases is between 200 and 300,which is approximately 10 times more than that with only fixed stations.

Table 1 Comparison between the number of phases at permanent stations and the number of phases adding ChinArray

Figure 6 Ray path coverage distribution of a ML2.1 event

The integration of mobile array data can also expand the recording scale of vibrations triggered by seismic events.For theML3.3 earthquake,we compare the distribution of seismic phases with the epicentral distance.Fig.7 shows that when only permanent stations are used,the results can identify the distribution of seismic phases as far as about 290 km.After combining the array data,the phases picked can reach as far as 600 km,and it increases greatly as the epicentral distance getting further.

2.2 The Accuracy of Phase Picking

Manual phase picking is relatively precise,but it requires massive labor cost.On average,the time cost for one phase is approximately 10-30 seconds,including identification,confirmation,labeling,review,and archiving.By contrast,such complex process only takes no longer than 1 second in automatic picking.Phase picking not only needs to be efficient,but also requires high accuracy for subsequent applications.Therefore,P-wave phase picking error is required to be within 0.5 second and the S phase picking error is required to be within 1 second.

Since the manual picking results are not recorded in the mobile stations,we also apply automatic picking on the permanent local network records with the same parameters and procedures as those for mobile network,and then compare the results with the bulletin from the local network.Theoretically,the accuracy of automatic picking of mobile network and local network should be consistent.It is indicated that for the events withML≥3.0,81.5%deviation time of P-wave arrival is less than 0.5 s,and for distance within 200 km,the value is 83.8%;70.5% deviation time of S-wave arrival is less than 1s,for distance within 200 km,the value reaches 80.7% (Fig.8(a)).For the events withML≥2.0 andML＜3.0,76.5% deviation time of P-wave arrival is less than 0.5 s,and with distances within 200 km the percentage is 76.9%;69.5% deviation time of S-wave arrival is less than 1s,and for distance within 200 km the value is 79.9% (Fig.8(b)).Obviously,the picking error for S-waves significantly increases with distances greater than 200 km,while P-wave picking do not show much difference.It's mainly because that after 200 km,the S-wave phase is relatively blurring and with lower signal-to-noise ratio,which is also a limitation for manual picking.Computer picking is often at the place where the amplitude changes the most,as performed in characteristic curve like AIC,STA/LTAs.In these cases,the smoothness of S-wave is better than the P-wave's.Thus,the manually picked points are usually prior to the computer-recognized time.In actual work,seismic phase identification is often corrected based on manual experience according to the seismic phase cycle changes.

Figure 7 Distance distribution as a function of phase time Black traces are recorded by permanent stations,while blue traces are from dense array stations

In seismic phase recognition the signal to noise ratio is an important factor that determines the accuracy of travel time picking no matter the phase is picked manually or automatically.Using fixed stations to automatically identify the signal-to-noise ratio of seismic phases,the relationship between the pickup error and the signal-to-noise ratio is analyzed.In the signal-to-noise ratio calculation,we select the signal waveform within T seconds after the arrival time of the signal to calculate the sum of squared amplitudes,while the noise is selected as T seconds before the signal arrives.The waveform is filtered through the frequency of 1-20 Hz,and the signal-to-noise ratio is calculated as

in whicht0is the signal arrival time,x(t) is the bandpass filtered waveform.The window time T is 0.5 second for P-wave,and 1 second for S-wave.

Figure 8 (a) Arrival time deviation of events with ML≥3.0;(b) Arrival time deviation of events with 2.0≤ML＜3.0

Fig.9 shows the relationship between the SNR and the picking error of 8 819 P waves generated by 267 earthquakes withML＞3.0.It can be seen that when the signal-to-noise ratio is relatively low,the pick-up error is relatively high.As the SNR increases,the pickup error gradually decreases.Overall,most errors are distributed within 0.5s.

Figure 9 Relationship between picking error and signal to noise ratio (SNR) of events with ML≥3.0

Figure 10 Plots showing T travel time against distance distribution

We extract the travel times of 2 006 earthquakes withML＞2.0 recorded by ChinArray and local network.Finally,718 573 P-wave phases and 512 035 S-wave phases are extracted from the dense array,which are 28 times and 22 times of those in the bulletin reported by local network respectively.

Fig.10 shows the travel time distribution of various phases within the epicentral distance of 1 200 km.It can be seen that compared with the number of phases at fixed stations,the number of various phases increases greatly after the automatic picking of phases combined with the array,and the corresponding epicentral distance also expands.Among various types of seismic phases,the number of Sn and S phases increase most,but the travel time deviation of the S phase is higher than that of the P-wave phase,which is associated with the complex shape of S-wave and the difficulty of identification.For different distances,the automatic recognition of seismic phases and the artificial recognition of seismic phases in the range of 0-400 km are relatively consistent.In the range of 600-700 km,the recognition of Pg and Sg phases becomes inferior,while the number of P and S phases increases correspondingly.The first arrivals of Pn and Sn are difficult to be recognized in both manual picking and machine learning results.By the assistance of the predicted travel time,in this paper,we can also effectively identify the first wave signals,especially the Pn wave.The number of Pn phases greatly increases in the range of 150-500 km,providing excellent data for subsequent imaging studies on the top of the Moho (Li Zhiwei et al.,2012).

Figure 11 Results of phase picking at several stations

The identification of seismic phase is closely related to the signal-to-noise ratio and the complexity of the waveform.Although some seismic events are difficult to be identified manually,the automatic identification method provides the possibility to identify Pg,Sg,Pn and Sn phases (Fig.11(a),(b)).Meanwhile,there are some relatively clear phases on the waveform that are not automatically recognized (Fig.11(c),(d)).In the next work,it is hopeful to combine the traditional automatic recognition method with the newly developed method based on machine learning so that the recognition rate may be further improved.

3 CONCLUSION

In this paper,the improved AR-AIC method is used to automatically pick the phase arrival time fromML＞2.0 earthquakes recorded by dense arrays deployed in the Yunnan area.Credit to the spatial coincidence of permanent network and mobile stations,existing fixed stations can be used to obtain earthquake prediction travel time and rapid phase recognition.A variety of seismic phases such as Pg,Sg,P,S,Pn,Sn in the mobile array can be obtained by using related methods.The automatic picking accuracy for most of the P-wave phases is within 0.5 s,and the S phase picking accuracy is within 1 s.The number of P-wave and S-wave phases recognized on the obtained mobile dense array stations is 28 times and 22 times of that of the original fixed station,respectively.

The method proposed in this paper has fast processing speed,and much higher efficiency than manual picking,providing foundation for the following body wave tomography.In order to achieve fast recognition,this paper uses a one-dimensional velocity model to predict arrival times.Due to the fact that Yunnan area has a complex velocity structure,three-dimensional velocity structure information should be added in the future.It is expected to further improve the recognition accuracy and perform detailed analysis on the travel time of events with smaller magnitudes.