2. Mengcheng National Geophysical Observatory, University of Science and Technology of China, Mengcheng 233500, Anhui, China;
3. Key Laboratory of Seismic Observation and Geophysical Imaging, Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
The Binchuan basin in Yunnan is a typical transtensional structure located in a pullapart basin formed since the Late Cenozoic (Luo Ruijie et al., 2015). This basin is surrounded by high mountains and lakes, with a large elevation difference (about 1000m). Many faults with different trending and movement patterns are distributed in this area, especially the Red River fault (RRF) and the Chenghai fault (CHF), which are the two major faults accommodating significant tectonic deformation and posing serious seismic hazards in this region (Allen C. R. et al., 1984). Due to its special tectonic location and complex geological structure, many earthquakes have occurred in this region and neighboring areas. The largest earthquake occurred on the CHF was the 1515 M7.7 Yongsheng earthquake (Huang Xiaolong et al., 2018). According to the positive correlation between the intensity of seismic and fault activities, former studies considered the Binchuan area as a M7.0 earthquake potential area (Wang Fan et al., 2015). Therefore, the investigation of subsurface structure is necessary to help us better understand the complex fault zone structure and strong motion characteristics.
During the 1980s, numerous seismic experiments have been conducted in Yunnan, such as wideangle seismic reflection and refraction studies (Zhang Zhongjie et al., 2005; Zhang Xi et al., 2009), receiver function analysis (e.g., Hu Jiafu et al., 2005; Wang Chunyong et al., 2010), seismic tomography (e.g., He Zhengqin et al., 2004; Yao Huajian et al., 2008; Huang Zhouchuan et al., 2015), and joint seismic inversion (e.g. Liu Qiyuan et al., 2014; Bao Xuewei et al., 2015). In addition to 3D structures, the China Earthquake Administration (CEA) has built a Fixed Airgun Signal Transmission Station (FASTS) in Binchuan, Yunnan in 2011 (Wang Baoshan et al., 2012; Chen Yong et al., 2017) to better capture temporal velocity and attenuation changes of fault zones due to earthquakes because of high repeatability of waveforms generated by airgun sources. During March 21 to May 30, 2017, a 2D shortperiod dense array (about 400 portable seismic stations) was also performed by the CEA (Fig. 1) (Wang Baoshan et al., 2018; Xu Y. et al., 2018) around the Binchuan basin in Yunnan. The purpose of this experiment was to obtain the highresolution shallow crustal structure of the Binchuan basin and surrounding area as well as to understand complexities of seismic wave propagation and amplification in this region.
Seismic rays would follow the greatcircle path in the horizontally isotropic layered medium. However, in real situations, seismic wave propagation exhibits complex characteristics due to strong heterogeneity of underground media with mountains and basins, for instance, in the Binchuan area. In this study, we intend to focus on investigating the characteristics of seismic wave propagation in this region using airgun signals recorded by dense array stations in this experiment (Wang Baoshan et al., 2018). Two classical methods are used in this study, the first method is the beamforming method (e.g., Rost S. et al., 2002) based on dense array stations. The second is the polarization analysis method (e.g., Jurkevics A., 1988) based on single station waveforms.
1 DATA AND PREPROCESSINGThe seismic data we analyze is part of the dense seismic experiment, carried out at Binchuan in 2017 by CEA. Airgun shots produced by an airgun array (four airguns with a volume of 2000 inch^{3} each) fired from the Binchuan FASTS were used in this study. During the observation of this dense array in the experiment, ~1400 airgun shots were fired. The observation system consists of 400 portable shortperiod seismic stations deployed in a rectangular region (Fig. 1), which spans an area about 25.6°—25.9°N and 100.3°—100.7°E. Considering the limitation of research methods we used, we only chose 67 stations for their dense distribution (see stations in red and blue in Fig. 1).
Because the airgun source is repeatedly fired at the same location, it can generate highly repeatable seismic waves (e.g., Lin Jianmin et al., 2008; Chen Yong et al., 2017). We can therefore enhance the signaltonoise ratio (SNR) of the waveform data and obtain highquality seismic waveforms through stacking. We cut the verticalcomponent data from the excitation time to 20s later, removed the mean and trend of each trace, and linearly stacked all the waveforms recorded at each receiver. Previous studies have shown that the dominant frequency band of the excited airgunshot signal is between 3Hz and 6Hz (e.g., Lin Jianmin et al., 2008; Yang Wei et al., 2013). Therefore, we bandpassfiltered the waveform data in the 0.510.0Hz frequency band (or the 0.12.0s period band) to suppress noise. An example of waveforms from the Binchuan FASTS after the linear stack is shown in Fig. 2(a).
2 METHOD AND RESULTSIn this section, we will use two methods to analyze the propagation direction of the P waves that are excited by the Binchuan FASTS. The first method is the arraybased frequencydomain beamforming analysis (e.g, Rost S. et al., 2002), which will determine the average back azimuth of the incoming (plane) wave to the dense array. The second method is the singlestation based polarization analysis (e.g., Jurkevics A., 1988), which directly determines the back azimuth of the incoming wave to the station. In each subsection, we will first introduce each method and then provide the results.
2.1 Beamforming AnalysisThe Beamforming method (e.g., Rost S. et al., 2002) in the frequency domain was used to calculate the apparent slowness vectors of the first arrival P waveforms generated by the airgun shots after stacking. This method is based on the farfield plane wave propagation hypothesis. We performed a gird search in the apparent slowness space (slowness and back azimuth) and intended to maximize the beamforming energy of the aligned P waveforms.
For each selected slowness S_{i} and back azimuth θ_{j} at grid node (i, j), we can calculate the phase delay Δφ_{k, ij} at each station k for central frequency ω as
$ \Delta \varphi_{k, i j}=\omega \cdot \Delta t_{k, i j} $  (1) 
where Δt_{k, ij} represents the time delay for each station k with respect to the reference location. In order to ensure the integrity of effective signal in the limited time window, we first aligned the waveforms before doing beamforming analysis (Fig. 2(b)). So the time delay Δt_{k, ij} can be divided into two parts: time delay δt_{1k, ij} caused by the location of each station and the time shift δt_{2k} for each trace for aligning the waveforms:
$ \Delta t_{k, i j}=\delta t_{1 k, i j}+\delta t_{2 k} $  (2) 
$ \delta t_{1 k, i j}=l_{x} S_{i, x} \sin \left(\theta_{j}\right)+l_{y} S_{i, y} \cos \left(\theta_{j}\right) $  (3) 
where l_{x} and l_{y} are the x and y component of the distance between station k and the array center, and S_{i, x} and S_{i, y} are, respectively, the x and y component of the slowness for node (i, j).
The position of the node where the maximum beamforming energy P(ω) is reached provides an estimate of the apparent slowness and back azimuth of the incoming wave front.
$ P(\omega)=\left\sum_{k=1}^{N} X_{k}(\omega) \cdot e^{i \Delta \varphi_{k, i j}}\right^{2} $  (4) 
In order to apply this method to array data, we have to select an adequate set of parameters, including the filter, time window and apparent slowness grids. We used a bandpass, zerophase Butterworth filter in the 28Hz frequency band, where most of the energy is concentrated for the airgun signal. We selected a window length of 0.5s for the analysis. This window represents about 2.5 periods at the dominant frequency of 5Hz, which has been suggested as the optimum window length for this type of analysis (Almendros J. et al., 1999). The apparent slowness and azimuth grid interval are chosen to be 0.001s/km and 1°, respectively.
To satisfy the plane wave hypothesis, we only selected twelve subarrays (67 stations) from the dense part of the 2D array of Binchuan to perform beamforming analysis. In this situation, we assume a homogeneous structure under each subdense array. These subdense arrays are distributed from north to south along the Binchuan Basin. The propagation distance is in the range of 719km. Each subarray consists of ~7 stations in a circular region with a radius of about 1.5km. Although we have chosen the denser part of this Binchuan array, the distances between stations are still too large. That will result in a problem of spatial aliasing, that is, the beamforming energy images are not clean enough, typically with multiple peaks. There will be relatively important secondary peaks produced by this spatial aliasing. To ensure the reliable results on spatial aliasing, we did synthetic tests for each subarray using theoretical waveforms. Fig. 3 is an example of synthetic tests. As we see smallaperture arrays show wider peaks (Fig. 3(c)), resulting in lower resolution and larger uncertainties, whereas largeaperture arrays with sparse station distribution display a narrow main peak with some secondary peaks produced by spatial aliasing (Fig. 3(b)). After doing this test, we can find only one peak in northwest represents the true main peak for wave propagation direction and slowness (Fig. 3(e)), and the others are all produced by spatial aliasing.
The procedure described above allows us to calculate the average back azimuth that characterizes the first arrivals of P waveforms generated by the airgun shots to the subseismic arrays. To extract the most significant solutions among this data set, we consider only those stable and reliable solutions. Fig. 4 shows a summary of the results for the twelve selected subarrays. In this figure we place the beamforming result at the right of the corresponding subarray. We also mark the greatcircle path and the real direction of seismic wave propagation to the subarray in order to illustrate a visual representation of the raybending. We observe the existence of back azimuth anomalies from 0° to 30°. The variation angle depends on the position of the subarray. The regions of the 3rd6th and 8th12th subarray show major anomalies of wave propagation direction.
2.2 Polarization AnalysisThe analysis of multicomponent recordings can provide an estimate of polarization properties of seismic arrivals. Knowledge of the polarization state can be used to distinguish between surface wave arrivals and body waves, or more generally for the purpose of wavefield decomposition. The method we used is based on the analysis of the covariance matrix of multicomponent recordings and principal component analysis using singular value decomposition.
The processing is carried out in the time domain. The polarization within a time window is estimated as follows. Let us consider the vector valued signal S(t)=[S_{x}(t), S_{y}(t), S_{z}(t)]^{T}. We assume that this signal represents a threecomponent seismic recording. Here we associate S_{x}(t), S_{y}(t), S_{z}(t) with the usual radial (inline), transverse, and vertical components of a triaxial recording. Polarization analysis of S(t) can be carried out by eigen analysis of the crossenergy matrix M, defined from the data vector as (Flinn E.A., 1965; Kanasewich E.R., 1981; Jurkevics A., 1988; Jackson C. M. et al., 1991; Kulesh M. et al., 2007):
$ {\bf M}(\xi)=\left[\begin{array}{ccc}{I_{x x}(\xi)} & {I_{x y}(\xi)} & {I_{x z}(\xi)} \\ {I_{y x}(\xi)} & {I_{y y}(\xi)} & {I_{y z}(\xi)} \\ {I_{z x}(\xi)} & {I_{y z}(\xi)} & {I_{z z}(\xi)}\end{array}\right] $  (5) 
In the equation above
$ I_{k m}(\xi)=\frac{1}{T} \int\limits_{\xiT / 2}^{\xi+T / 2}\left[S_{k}(\tau)\mu_{k}(\xi)\right]\left[S_{m}(\tau)\mu_{m}(\xi)\right] d \tau $  (6) 
where k, m=(x, y, z).The above covariance matrix is computed within a time window of length T and around its time centre ξ; μ_{k}(ξ) is the mean value of the signal component S_{k} in the window T.
The eigenanalysis performed on M(ξ) yields the principal component decomposition of the energy for the time window T. Such a decomposition produces three eigenvalues λ_{1}(ξ)≥λ_{2}(ξ)≥λ_{3}(ξ) and three corresponding eigenvectors v_{k}(ξ) that fully characterize the magnitudes and directions of the principal components of the ellipsoid that approximates the particle motion in the considered time window. Purely rectilinear ground motion (Pwave) has only one nonzero eigenvalue. The azimuth can be expressed by:
$ \alpha(\xi)=\arctan \left(v_{1, y}(\xi) / v_{1, x}(\xi)\right) $  (7) 
To effectively use the standard covariance method for polarization analysis, one needs to start by selecting the appropriate length of the time window T. Selecting this window may not be easy. If the window length is too short, the stability of the polarization analysis will be low due to a higher sensitivity to local information, whereas we can hardly find the polarity because of the insensitivity to effective signals with too long a window length. The optimal length of the time window must depend on the main frequency of the effective signal in the time window. We selected the window length based on three criteria: (1) only one wave arrives in the time window; (2) the SNR (signal to noise ratio) should be the largest; (3) the effective signal must be completely contained.
Fig. 5 gives an example of the polarization analysis. We plot the seismic records which are rotated to the great circle path and the real path (estimate by polarization analysis), respectively (Fig. 5(a)). The particle motions in time window also shown in Figs. 5(b) and 5(c). We can clearly find that the particle motion of transverse (T) component is corrected to zero in time window after rotated to the real path.
Fig. 6 shows the result of the polarization analysis for each station. We find that the back azimuth anomaly of the Pwave propagation for each station is nearly the same as the beamforming result of the corresponding subarray, which confirms that both methods yield reliable estimates of wave propagation direction. There exists small differences for some stations, which may be caused by the quality of effective signals.
In this study we used 12 subseismic arrays to measure the back azimuth anomalies of the first arrivals generated by the Binchuan FASTS. Different values of back azimuth anomalies have been detected from beamforming and polarization analysis. As we know, in laterally homogeneous media, seismic rays are confined in a vertical plane that includes the source and the receiver (or approximately the subarray centre). Thus the expected back azimuth anomaly is zero everywhere. Any deviation from zero indicates the presence of lateral heterogeneities that affect the ray path.
Apparently, the Binchuan region is not laterally homogeneous. The seismic wave fronts are affected by these heterogeneities. Seismic waves speed up in highvelocity regions and slow down in lowvelocity areas, producing distorted wave fronts and twisted ray paths. For example, if a lowvelocity body was located between the source and the receiver, seismic rays would follow the fastest path, turning around the lowvelocity body instead of across it.
As we have discussed above, those detected azimuth anomalies indicate lowvelocity bodies (sediments in the Binchuan basin) between FASTS and stations or highvelocity bodies (crystalline rocks of the mountains) besides the theoretical path. According to numerical simulation of ray propagation based on a simple 1D velocity model, we find that the depth of Pwave transmissions will not be greater than 1km at the current propagation distance (< 20km). It means that all these azimuth anomalies we observed are attributed to the heterogeneities in the very shallow crust (< 1km).
Apart from the evidence about the highly heterogeneous structure of the Binchuan basin, we have to emphasize that we are using waveform data provided by an active seismic experiment. We are therefore dealing with highfrequency and very shallow seismic source. In this sense, our azimuth anomalies estimated from the data set may constitute a more complicated situation. In general, highfrequency seismic waves are more affected by topography and heterogeneities than lowfrequency waves, due to shorter wavelengths. Moreover, the shallow structure is expected to be more heterogeneous with large velocity variations.
We should point out that there might be some component azimuth errors caused by careless deployments, which is important in single station polarization analysis. In order to validate the reliability of results, we use the teleseismic body waves from an M_{W}6.8 earthquake in the South Pacific region (epicentral distance ~8, 500km) to evaluate this issue. Fig. 7 shows the histograms of azimuth anomalies from the polarization analysis using airgun shots and teleseismic events. Calculated results indicate that more than twothirds of stations (45 stations) show at most 5 degrees, about onefifth of stations (14 stations) show not exceeding 10 degrees and only 8 stations show more than 10 degrees deviation from the larger than 5 degrees but circle path (Fig. 7(b)). Therefore, the azimuth errors caused by careless deployments have little influence on polarization analysis.
In this study, the polarization analysis based on single station and beamforming analysis based on dense array both reveal similar but complex characteristics of seismic wave propagation in the Binchuan basin. This situation is due to the strong structural heterogeneity of underground media. Our research can provide a better understanding of the complex geologic structures in this area and provide guidance for detecting concealed faults and distribution of velocity anomalies. In the future, combining the characteristics of seismic wave propagation with highresolution tomography can provide a reliable model for this complex region. The obtained propagation directions of seismic waves can also be used to access the reliability of future tomographic models in this region.
ACKNOWLEDGEMENTSThe authors thank the editors and two anonymous reviewers for valuable comments that helped improve this article.
Allen C.R., Gillespie A.R., Han Yuan, Sieh K.E., Zhang Buchuan, Zhu Chengnan. Red River and associated faults, Yunnan Province, China:Quaternary geology, slip rates, and seismic hazard[J]. Geological Society of America Bulletin, 1984, 95(6): 686700. DOI:10.1130/00167606(1984)95<686:RRAAFY>2.0.CO;2 
Almendros J., Ibáñez J.M., Alguacil G., Del Pezzo E. Array analysis using circularwavefront geometry:an application to locate the nearby seismovolcanic source[J]. Geophysical Journal International, 1999, 136(1): 159170. DOI:10.1046/j.1365246X.1999.00699.x 
Bao Xuewei, Sun Xiaoxiao, Xu Mingjie, Eaton D.W., Song Xiaodong, Wang Liangshu, Ding Zhifeng, Mi Ning, Li Hua, Yu Dayong, Huang Zhouchuan, Wang Pan. Two crustal lowvelocity channels beneath SE Tibet revealed by joint inversion of Rayleigh wave dispersion and receiver functions[J]. Earth and Planetary Science Letters, 2015, 415: 1624. DOI:10.1016/j.epsl.2015.01.020 
Chen Yong, Wang Baoshan, Yao Huajian. Seismic airgun exploration of continental crust structures[J]. Science China Earth Sciences, 2017, 60(10): 17391751. DOI:10.1007/s1143001690966 
Flinn E.A. Signal analysis using rectilinearity and direction of particle motion[J]. Proceedings of the IEEE, 1965, 53(12): 18741876. DOI:10.1109/PROC.1965.4462 
He Zhengqin, Su Wei, Ye Tailan. Seismic tomography of Yunnan region using shortperiod surface wave phase velocity[J]. Acta Seismologica Sinica, 2004, 17(6): 642650. DOI:10.1007/s1158900400037 
Hu Jiafu, Su Youjin, Zhu Xiongguan, Chen Yun. Swave velocity and Poisson's ratio structure of crust in Yunnan and its implication[J]. Science in China Series D:Earth Sciences, 2005, 48(2): 210218. DOI:10.1360/03yd0062 
Huang Xiaolong, Wu Zhonghai, Wu Kungang. Surface rupture of the 1515 Yongsheng earthquake in Northwest Yunnan, and its Seismogeological implications[J]. Acta Geologica Sinica, 2018, 92(4): 13241333. DOI:10.1111/17556724.13629 
Huang Zhouchuan, Zhao Dapeng, Wang Liangshu. P wave tomography and anisotropy beneath Southeast Asia:Insight into mantle dynamics[J]. Journal of Geophysical Research:Solid Earth, 2015, 120(7): 51545174. DOI:10.1002/2015JB012098 
Jackson G.M., Mason I.M., Greenhalgh S.A. Principal component transforms of triaxial recordings by singular value decomposition[J]. Geophysics, 1991, 56(4): 528533. DOI:10.1190/1.1443068 
Jurkevics A. Polarization analysis of threecomponent array data[J]. Bulletin of the Seismological Society of America, 1988, 78(5): 17251743. 
Time Sequence Analysis in Geophysics[M]. Edmonton, Alta.: University of Alberta, 1981.

Kulesh M., Diallo M.S., Holschneider M., Kurennaya K., Krüger F., Ohrnberger M., Scherbaum F. Polarization analysis in the wavelet domain based on the adaptive covariance method[J]. Geophysical Journal International, 2007, 170(2): 667678. DOI:10.1111/gji.2007.170.issue2 
Lin Jianmin, Wang Baoshan, Ge Hongkui, Tang Jie, Zhang Xiankang, Chen Yong. Study on large volume airgun source characteristics and seismic phase analysis[J]. Chinese Journal of Geophysics, 2008, 51(1): 206212 (in Chinese with English abstract). 
Liu Qiyuan, Van Der Hilst R.D., Li Yu, Yao Huajian, Chen Jiuhui, Guo Biao, Qi Shaohua, Wang Jun, Huang Hui, Li Shuncheng. Eastward expansion of the Tibetan Plateau by crustal flow and strain partitioning across faults[J]. Nature Geoscience, 2014, 7(5): 361365. DOI:10.1038/ngeo2130 
Luo Ruijie, Wu Zhonghai, Huang Xiaolong, Huang Xiaojin, Zhou Chunjing, Tian Tingting. The main active faults and the active tectonic system of Binchuan area, northwestern Yunnan[J]. Geological Bulletin of China, 2015, 34(1): 155170. 
Rost S., Thomas C. Array seismology:methods and applications[J]. Reviews of Geophysics, 2002, 40(3): 21. 
Wang Baoshan, Ge Hongkui, Yang Wei, Wang Weitao, Wang Bin, Wu Guohua, SuYoujin. Transmitting seismic station monitors fault zone at depth[J]. Eos, Transactions American Geophysical Union, 2012, 93(5): 4950. 
Wang Baoshan, Tian Xiaofeng, Zhang Yunpeng, Li Yulan, Yang Wei, Zhang Bo, Wang Weitao, Yang Jun, Li Xiaobin. Seismic signature of an untuned largevolume airgun array fired in a water reservoir[J]. Seismological Research Letters, 2018, 89(3): 983991. DOI:10.1785/0220180007 
Wang Chunyong, Lou Hai, Silver P.G., Zhu Lupei, Chang Lijun. Crustal structure variation along 30°N in the eastern Tibetan Plateau and its tectonic implications[J]. Earth and Planetary Science Letters, 2010, 289(3/4): 367376. 
Wang Fan, Wang Min, Wang Yanzhao, Shen Zhengkang. Earthquake potential of the SichuanYunnan region, western China[J]. Journal of Asian Earth Sciences, 2015, 107: 232243. DOI:10.1016/j.jseaes.2015.04.041 
Xu Y., Wang B., Wang W., Zhang B., Sun T. Multiple seismological evidences of basin effects revealed by Array of Binchuan (ABC), northwest Yunnan, China[J]. Earthquake Science, 2018, 31: 110. DOI:10.29382/eqs201800011 
Yang Wei, Wang Baoshan, Ge Hongkui, Wang Weitao, Chen Yong. The active monitoring system with large volume airgun source and experiment[J]. Earthquake Research in China, 2013, 29(4): 399410 (in Chinese with English abstract). 
Yao Huajian, Beghein C., Van Der Hilst R.D. Surface wave array tomography in SE Tibet from ambient seismic noise and twostation analysisⅡ. Crustal and uppermantle structure[J]. Geophysical Journal International, 2008, 173(1): 205219. DOI:10.1111/gji.2008.173.issue1 
Zhang Xi, Wang Yanghua. Crustal and upper mantle velocity structure in Yunnan, Southwest China[J]. Tectonophysics, 2009, 471(3/4): 171185. 
Zhang Zhongjie, Bai Zhiming, Wang Chunyong, Teng Jiwen, Lü Qingtian, Li Jiliang, Sun Shanxue, Wang Xinzheng. Crustal structure of Gondwana and Yangtzetyped blocks:an example by wideangle seismic profile from Menglian to Malong in western Yunnan[J]. Science in China Series D:Earth Sciences, 2005, 48(11): 1828. DOI:10.1360/03yd0547 