2. Shandong Earthquake Agency, Jinan 250102, China;
3. School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China
Since February 14, 2017, seismic activities have occurred continually in Chang Island area, Shandong Province. Until August, 2017, more than 2, 300 earthquakes have been recorded in the local catalog, including 43 M_{L}≥3.0 events. The largest earthquake was the M_{L}4.5 earthquake occurred on March 3, 2017 (Fig. 1). The frequency and intensity of the Chang Island earthquake swarm exceed those before the 1976 Tangshan M_{S}7.8 earthquake in the same area, which is unwonted in the history. Chang Island area has always been considered as an information window of the seismicity investigation in North China (Su Luansheng, 1993). Therefore, detailed analysis of the characteristics of Chang Island earthquake swarm is conducive to understand the seismogenic mechanism of earthquakesand further assist estimate the seismic risk in this area.
The Chang Island swarm is located at the eastern end of the NWWoriented ZhangjiakouPenglai fault zone. The faults in this area are complex and mainly consist of the NWWoriented PenglaiWeihai fault and several NEtrending small faults (Fig. 1). Due to the difficulty of tectonic stress measurement as well as the complexity of the seismic triggering process, it is difficult to establish a physical model for accurate description. However, the essential rules can be explored from complex phenomena using statistical methods. Since Ogata (Ogata Y., 1988) proposed the Epidemic Type Aftershock Sequence (ETAS) model based on the "OmoriUtsu" law (Omori F., 1894; Utsu T., 1961), the model has received extensive attention and gained considerable development. With the support of ETAS model, the parameters of seismic sequence can be accurately estimated, the variation law of seismic activity can be effectively found (Ogata Y., 1992, 2001; Jiang Haikun et al., 2007; Jiang Changsheng et al., 2013(a), 2013(c), 2015, 2017), the relative quietness of seismic activity can be monitored (Zhuang Jiancang, 2000; Jiang Changsheng et al., 2013(b)), slight stress change can be detected (Helmstetter A. et al., 2003; Lei Xinglin et al., 2008(a), 2017; Jia Ke et al., 2014, 2018), and induced earthquakes related to fluid action can be identified (Hainzl S. et al., 2005; Lei Xinglin et al., 2008(b); Long Feng et al., 2010; Jiang Haikun et al., 2011, 2012; Peng Yajun et al., 2012)
The statistical seismic modelling method depends largely on the earthquake catalogue. The seismic sequence activity characteristics and the aftershock probability prediction of large earthquakes are usually affected by the magnitude of completeness (M_{c}) (Jiang Changsheng et al., 2013(c); Zhuang Jiancang et al., 2017). Also, the artificial setting of the magnitude of completeness may affect the estimation of ETAS model parameters (Jiang Haikun et al., 2012; Li Zhichao et al., 2014), which mainly depends on the monitoring ability of regional seismic network. In areas like Chang Island, the signaltonoise ratios of seismic waveforms are usually low, and seismic stations are sparse or the distribution is not ideal. Thus, the clustering occurrence of earthquake sequences in a short time usually results in the omission of earthquake in the catalog. However, the reliability of analysis can be improved by using a more complete earthquake catalogue (Schaff D. P., 2008).
The template matching filtering technique based on waveform crosscorrelation is suitable for microseismic detection (Gibbons S. J. et al., 2006), which has developed rapidly and has been widely used recently. It was applied for detecting tremors (Shelly D. R. et al., 2007), studying the aftershock sequence changes (Peng Zhigang et al., 2009; Hou Jinxin et al., 2017; Wang Peng et al., 2017), determining the seismogenic structure (Yang Hongfeng et al., 2009; Tan Yipei et al., 2016), analyzing the stress changes in the source region (Meng Xiaofeng et al., 2012), identifying repeating earthquakes (Ma Tengfei, 2015), and monitoring nuclear explosions (Zhang Miao et al., 2015).
In this paper, a GPUbased template matching method is used to detect the missing earthquakes of the Chang Island earthquake swarm. The statistical parameters of the ETAS model of the earthquake swarm are analyzed based on the completed earthquake catalogue, and the possible seismogenic mechanism of the Chang Island earthquake swarm is then discussed.
1 GEOLOGICAL BACKGROUND AND SEISMIC ACTIVITY IN THE STUDY AREAThe Chang Island earthquake swarm is located on the PenglaiWeihai fault zone, which is the southeast segment of the NWWtrending ZhangjiakouBohai fault zone (Xu Jie et al., 1998). It is composed of a series of NWtrending faults developed between Penglai and Weihai sea area in the northern Shandong Peninsula. The major part is located between the Chang Island and the Dazhu Island; most faults are normal and strikeslip faults, some of which are thrusting (Zheng Jianchang et al., 2018). The PenglaiWeihai fault zone has hosted many historical earthquakes. Among them, the Chang Island M_{S}7.0 earthquake in 1948 occurs in the west of the fault zone and the Weihai M_{S}6.0 in 1948 occurs in the east of the fault zone (Wang Zhicai et al., 2006).
From February 14 to August 30, 2017, an earthquake swarm with 2, 377 earthquakes (from the Shandong Seismic Network) has occurred in the Chang Island, including 1, 459 M_{L}1.01.9 earthquakes, 271 M_{L}2.02.9 earthquakes, 39 M_{L}3.03.9 earthquakes and 43 M_{L}≥4.0 earthquakes. Among them, the largest M_{L}4.5 earthquake occurred on March 3. There are 16 seismic stations within 150km from the swarm, all are threecomponent broadband seismographs with a sampling rate of 100Hz. They are mainly distributed in the land area to south of the swarm. In the following, we select Chang Island Station (CHD), Beihuangcheng Station (BHC), Longkou Station (LOK), and Yantai Station (YTA), which are close to and well distributed to the earthquake swarm, to detect missing earthquakes by using the threecomponent continuous waveforms.
2 METHODS 2.1 Template Matching MethodThe procedure of the template matching method used in this paper follows the work by Peng and Zhao (Peng Zhigang and Zhao Peng, 2009), and uses the similar processing steps as Hou and Wang (Hou Jinxin and Wang Baoshan, 2017). The procedure is shown in Fig. 2, which is roughly divided into four steps and is briefly described as follows:
(1) Data preprocessing and template selection. Using the threecomponent records of the four stations mentioned above (blue triangles in Fig. 1), the mean values of the continuous waveform data from the period February 9 to August 30, 2017 are removed. Then, the data are detrended and bandpass filtered to 210Hz. In total, there are 302 earthquakes with M_{L}≥2.0 and SNR>5 from Chang Island swarm selected as templates, and the P and Swave arrival times of the template earthquakes are obtained from the local catalog.
(2) The acquisition of crosscorrelation coefficients. The template waveforms are then crosscorrelated to the corresponding continuous waveforms. The correlating time window is 1s before and 5s after the P and Swave arrival times of the template. The crosscorrelation coefficient of the four stations are then shifted and averaged.
(3) The detection of missing earthquakes.We calculate the Median Absolute Deviation (MAD) of averagedcrosscorrelation coefficients, and select 12 times of the MAD as threshold. When the crosscorrelation coefficient is higher than the threshold, it is marked as a tentative event; the maximum value of the average correlation coefficient within 5s is set as a detected event.
(4) The determination of earthquake time and magnitude. Based on the assumption that the missing event has similar travel times as the corresponding template event, the origin time of the missing event is obtained. The ratio of maximum amplitude of the horizontal component of the missing event to the corresponding template event is used to estimate the magnitude of the missing earthquake.
2.2 ETAS ModelThe conditional intensity function λ(t) of the ETAS model is (Ogata Y., 1988; Jiang Haikun et al., 2012):
$ \lambda \left( t \right) = \mu \left( t \right) + K\sum\limits_{{t_i} < t} {\frac{{{e^{\alpha \left( {{M_i}  {M_z}} \right)}}}}{{{{\left( {t  {t_i} + c} \right)}^p}}}} , $  (1) 
The first item on the right μ(t) is the theoretical number of earthquakes per unit time after eliminating the clustering of earthquakes. It also indicates that the intensity of microseismic activity triggered by external force (fluid), which has no correlation with the selfexcitation process between earthquakes (Hainzl S. et al., 2005). The second item describes the contribution of the ith earthquake to the occurrence of following earthquakes. i traverse all earthquakes in the earthquake sequence. M_{z} is the reference magnitude no less than the magnitude of completeness M_{c}. p represents the attenuation factor of seismic sequence: large p means fast attenuation, while small pmeans slow attenuation. In the mainshockaftershock type earthquake sequence, when M_{z} takes the magnitude of the mainshock, the smaller α value has higher secondary aftershock excitation ability (Song Jin et al., 2009); when M_{z} takes the complete magnitude of M_{c}, α indicates the ability of triggering secondary aftershocks, and greater α indicates the stronger ability (Jiang Haikun et al., 2012). The initial physical implication of c is a very small positive number, ensuring that the denominator of the second term of equation (1) is not zero, the exact value of c is determined by the completeness of the data in a short time after the mainshock (Jiang Haikun et al., 2007) and the stress state of the source region (Narteau C. et al, 2009). k is the level of seismic activity associated with the magnitude of the mainshock and the b, p, α, c, μ parameters, where b is the scale factor in the GR law (Song Jin et al., 2009).
One observed seismic sequence can be fitted with multiple statistical models mathematically, andmultiple sets of qualified parameters can be obtained. Ogata (Ogata Y., 1988) uses the maximum likelihood method to estimate the parameters of the ETAS model and uses the Akaike information criterion (AIC) to minimize the optimal model parameters (Akaike H., 1974). The AIC criterion not only considers the pros and cons of the model fitting of the observed data, but also punishes the behavior of increasing the model parameters in an unrestricted manner for improving the fitting degree. The log likelihood of the ETAS model are:
$ \lg L = \sum\limits_{i = 1}^{{N_{{\rm{AIC}}}}} {\lg \lambda \left( t \right)  \int\limits_{{t_{\rm{b}}}}^{{t_{\rm{e}}}} {\lambda \left( t \right){\rm{d}}t} } , $  (2) 
$ {\rm AIC} =  \lg L + 2k, $  (3) 
In equation (2), t_{b} and t_{e} are the beginning and end time used for AIC calculation, and N_{AIC} is the number of earthquakes in the calculation time period. It is necessary to ensure that the earthquakes in the calculation period are basically complete. In formula (3), k is the number of parameters to be fitted. In this paper, k=5, and the model with minimum AIC is taken as the final model.
It can be seen from equation (1) that the seismic rate represented by ETAS model consists of two parts, the first term is independent of the aftershock excitation triggered by the external forces, while the second item is the seismic cluster caused by the "aftershock excitation aftershock". Therefore, through the separation of the two parts of the ETAS model, the ratio of the externalforcetriggered earthquake R_{b} can be expressed by formula (4), and the proportion of selfexcited earthquake is R_{t}=1R_{b}, The relationship between the triggering of external factors (fluid action) and seismic activity can be investigated according to their respective proportions. (Jiang Haikun et al., 2011).
$ {{R}_b} = \frac{{\int_{{t_{\rm{b}}}}^{{t_{\rm{e}}}} {\mu \left( t \right){\rm{d}}t} }}{{\int_{{t_{\rm{b}}}}^{{t_{\rm{e}}}} {\lambda \left( t \right){\rm{d}}t} }}, $  (4) 
Using the template matching technique, we detected 15, 286 earthquakes from February 9 to August 20, 2017, including 302 selfdetected templates and 14, 984 missing events. The seismic catalogue obtained by template matching method is hereinafter referred to as detection catalog. The event number in the detection catalog is more than 6 times compared with the number in Shandong Seismological Network Catalog. Fig. 3 is an example showing the detection of a missing M_{L}2.6 earthquake. The seismic rates of the detection and network catalogues are shown in Fig. 4, from which we suggest that there are few earthquakes when magnitude M_{L} > 2.5. The difference between the detection and network catalogues is related to the monitoring ability of the network in this area. The monitoring ability of one network can be described with the magnitude of completeness. To evaluate the magnitude of completeness, we used the maximum curvature method, and the corresponding magnitude with the maximum slope in the cumulative frequencymagnitude curve is usually selected as the magnitude of completeness(M_{c})(Wiemer S. et al., 2000). In practical applications, this magnitude usually corresponds to the magnitude of the noncumulative frequencymagnitude distribution with the largest seismic frequency (Huang Yilei et al., 2016). The estimated completeness magnitudes of the detection catalogue and the network catalogue are 0.5 and 1.0, respectively. The cumulative frequencymagnitude curve is then fitted using the maximum likelihood method, and the bestfitted b values for detection catalogue and the network catalogue are 1.04±0.01 and 0.79±0.02, respectively.
Based on the detection catalogue, we analyze the temporal evolution of the seismicity of Chang Island earthquake swarm. Judging from the variation of cumulative frequency and magnitude with time (Fig. 5), we divide the Chang Island earthquake swarm into four periods: in the first stage, the seismic activity with low magnitudes started to increase from February 13, 2017, the activity continued until the occurrence of the largest M_{L}4.5 earthquake (March 3, 2017). During this time, the network catalog revealed a relative quiet scenario (black point in Fig. 5); the second stage is mainly concentrated between March and April, 2017. After the M_{L}4.5 earthquake on March 3, the seismic activity increased obviously with higher intensity and frequency. However, the seismic activity gradually weakened after reaching the highest value in April, 2017; the third stage occurred in May, 2017, and a new round of enhancement occurred in the earthquake activity from May 2 to 4, and then gradually weakened. The fourth stage is from June to August, 2017. Similar to the third phase, the seismic activity appeared again in the process of enhancementattenuation, which slowly attenuated after a short time increase in early June. The detection catalogue and the network catalogue have the same seismic activity period, although the number of earthquakes in the lowmagnitude section is different. According to the development process of the earthquake swarm, each stage can be further divided into several substages (Fig. 5). In the following, the ETAS model will be used to model these different stages separately; we will obtain the variations of each model parameter and discuss the influence of external force (fluid) action and selfexcitation on the development of the earthquake swarm.
The cutoff magnitude has great influence in ETAS modelling (Jiang Haikun et al., 2012; Jiang Changsheng et al., 2013c) and is usually set greater than the magnitude of completeness. To model the seismic activities of all stages, we need to estimate the M_{c} of the corresponding stage. We use the maximum curvaturemethod and 'goodness of fit' test method to estimate the change of complete magnitude (M_{c}) with time (Fig. 6). From the results of the optimal solution (the black solid line in Fig. 6), it can be seen that the magnitude of completeness (M_{c}) is basically constant in the process of swarm development, all are no greater than 0.5. Therefore, in the subsequent parameter inversion of ETAS model, we can set the cutoff magnitude as 0.5.
According to the stages defined in Fig. 5, the earthquake catalogues of each stage M_{L}≥0.5 are taken, and the respective ETAS model parameters are estimated (Table 1). We select the first and second stages as examples for detailed analysis (Fig. 7).
For the first stage, the ratio of random components is 31.9%, indicating that 31.9% of the seismic activity is triggered by external forces, and 68.1% of the earthquakes are generated by selfexcitation of aftershocks. The avalue is 0.608, which is consistent with the activity characteristics of the earthquake cluster (0.350.85) (Ogata Y., 1992) and lower than the general tectonic earthquake (2.0). A small avalue indicates that the aftershock excitation ability is not strong. Meanwhile, low αvalue means low dependence on earthquake magnitude, therefore, the stress variation caused by the earthquake has only a slight change (Lei Xinglin et al., 2008a). The pvalue is larger than the typical tectonic earthquakes (around 1.0), indicating that the sequence decays rapidly. The fitting between theoretical and observed values of the model is poor, and the observed earthquake occurrence rate (blue solid line in Fig. 7(a)) is significantly higher than that of the theoretical curve (grey point line in Fig. 7(a)). This indicates that the selftriggering of Omori's law has only a slight promoting effect, and the proportion of seismic activity caused by external force may be underestimated (Lei Xinglin et al., 2008a).
The second stage corresponds to the period when the seismic activity rate gradually increases and then decreases. According to the activity rate, we divided stage Ⅱ into five subphases (Ⅱ_{1}Ⅱ_{5}, Fig. 5), of which Ⅱ_{2} has the highest seismic activity rate. The observations and theoretical values of these five subphases are well fitted, and the proportion of external force (fluid) action becomes significantly larger and gradually increases, from 46.9% in Ⅱ_{1} to 63.5% in Ⅱ_{3}, indicating the presence of active external force (fluid). The number of earthquakes triggered by fluids increases, and even more than half of the earthquakes are caused by external forces, instead of being influenced by "aftershocks triggering aftershocks" scenario in Omori's law. At the same time, the corresponding α values are higher than that in the first stage, suggesting that the excitation ability of aftershocks is enhanced. The pvalue remains high from stage Ⅰ to Ⅱ_{3}, but decreases significantly in stage Ⅱ_{4}and Ⅱ_{5}, indicating that the decay rate of the sequence decreases.
We calculate the bvalue of each stage using maximum likelihood method, and analyze the evolution of the bvalue of Chang Island earthquake swarm as well as R_{b}, α, p and μ from ETAS model in the whole process of swarm development (Fig. 8). The results show that, except for the individual activity stages, the R_{b}, μ and α are stable and positively correlated. All the values increase at first and then decrease (Fig. 8(a)(c)), the trend of bvalue and pvalue are not obviously synchronized, but their value are relative high at early stages and low in later stages. (Fig. 8(d)(e)). The μvalue is the background seismicity rate. Larger μ value represents stronger external force, and R_{b} is defined by the ratio of "triggered" earthquakes and "natural" earthquakes unrelated to the external force, which also indicates the triggering intensity (Formula 4). It can be seen from R_{b} values of different stages that the fluid plays a role in the development of the earthquake swarm to some extent, the effect of fluid increase at first and then weakened, indicating that the extent of external triggering intensity varies at different times.
We use the ratio of earthquakes induced by external force (fluid) and selfexcited earthquakes R_{b} (Formula (4)) to investigate the influence of external force on the swarm activity. Firstly (stage Ⅰ), the proportion of earthquake triggered by external force is not high (31.9%), then the triggering ratio increases gradually until stage Ⅱ_{3}, reaching the maximum (63.5%). Subsequently, the triggering effect of fluids i weakened obvioously.
The bvalue usually reflects the stress level of the study area: a low bvalue reflects high stress state, while a high bvalue corresponds to a low stress state. The high bvalue in the early stage of the earthquake swarm indicates that the stress at the beginning of the earthquake swarm is not high (Fig. 8(d)), but the largest M_{L}4.5 earthquake in the swarm occurred, which could be attributed to fluid triggering. R_{b} indicates the external force (fluid) effect and is generally thought to be related to the stress state of the source region (Peng Yajun et al., 2012; Jia Ke et al., 2014, 2018). In the initial stage of the earthquake swarm, the bvalue is high, while the R_{b} value is relatively low, and the actual earthquake occurrence rate at this stage may be significantly higher than the fitted theoretical curve (Fig. 7(a)) and the external force may be underestimated (Lei Xinglin et al., 2008a). The bvalue in other stages is also consistent with the process of increasing at first and then decreasing. Besides, the change of external force is also consistent. When the proportion of fluid action effect is large, the bvalue stays in a high level, indicating the fluid reduces the effective normal stress of the fault. Fluids greatly reduce the strength imposed by the fault hence even relatively low local stress can cause earthquakes.
The process of increase and decrease of the earthquake proportion of external force (fluid) triggering may be related to the fluid infiltration process. With the penetration of the fluid, the pore pressure gradually increases, the strength of the fault (fracture) decreases, leading to a larger possibility of earthquake occurrence. Fluid effect at the beginning of the second stage is obviously stronger than that in the firststage; however, after a certain period of time, the fluids tends to be saturated in the pore, and the change of pore pressure decreases gradually, which may weaken the triggering effect of external factors of fluid, and the ratio coefficient R_{b} decreases from 63.5% to less than 30%. Seismic activities in the third and fourth stages are mainly caused by the selfexcitation of earthquakes.
In each stage of swarm development, R_{b} is positively correlated with μ value, indicating the strength of external force. Larger value corresponds to stronger fluid effect. This is also consistent with the trend of earthquake selfexcitation (α). The parameters vary in different stages, however, the fluid interaction may generally have a greater impact on the early stage of the swarm, and the selfexcitation of the earthquake plays a more important role in the whole process of swarm development. In the initial stage of the swarm, the p value has been at a high level, result that indicates the swarm attenuates rapidly, suggesting that it may be related to the larger proportion of microearthquakes in the swarm. Alternatively, the p value is positively correlated with the focal depth and the temperature at the focal depth (Creamer F. H. et al., 1993; Nyffenegger P. et al., 2000; Song Jin et al., 2009). Seismic sequences with larger source depths attenuate faster, stress relaxation in regions with higher crustal temperature is faster, and aftershock activity attenuates relatively faster. The focal depth of Chang Island swarm is about 11km, which is larger than that of earthquakes induced by reservoirs and water injection(Jiang Haikun et al., 2012; Lei Xinglin et al., 2008b); Meanwhile, whether the source of fluids is related to highpressure crustal melting materials in the study area (Liu Lihua et al., 2015) remains questionable. Both of these factors may lead to the accelerated attenuation of earthquake swarms. Detailed changes of background earthquake occurrence rate and pvalue can be observed, and subsequent discussions can be completed in the way of iteration and continuous sliding window proposed by Peng Yajun et al. (2012).
5 CONCLUSIONSIn this paper, the GPUaccelerationbased template matching method is used to detect microseismic event from the continuous seismic waveforms of the Chang Island earthquake swarm from February 9 to August 20, 2017. A total of 15, 286 earthquake events are identified, of which 302 are selfdetection of template events and 14, 984 are missing earthquake events. Events in detected catalogue is more than 6 times of that of the network catalogue, and the magnitude of completeness of the earthquake swarm reduces from 1.0 to 0.5.
Based on the detected earthquake catalogue, the ETAS model is used to analyze the characteristics of the Chang Island earthquake swarm. It is preliminarily inferred that fluid triggering plays a role in the development of the Chang Island swarm to some extent, and has a greater impact in the early stage of the swarm, while the selfexcitation of the swarm also plays an important role in the whole development of the swarm. In this paper, the seismogenic mechanism of Chang Island swarm is explored only from the statistical point of view. As for the detailed mechanism and further analyzation of fluid triggering and the change of driving forces in the physical process of swarm development, more methods are needed.
ACKNOWLEDGEMENTThanks to the two reviewers for their valuable comments. In this paper, the template matching algorithm is developed by Professor Peng Zhigang, and the ETAS model is calculated by the GeoTaos software developed by Professor Lei Xinglin.
Akaike H.Akaike H. A new look at the statistical model identification[J]. IEEE Transactions on Automatic Control, 1974, 19(6): 716723. DOI:10.1109/TAC.1974.1100705 
Creamer F.H., Kisslinger C.Creamer F.H., Kisslinger C. The relation between temperature and the Omori decay parameter for aftershock sequences nearJapan[J]. EOS, 1993, 74(S43): 417. 
Gibbons S.J., Ringdal F.Gibbons S.J., Ringdal F. The detection of low magnitude seismic events using arraybased waveform correlation[J]. Geophysical Journal International, 2006, 165(1): 149166. DOI:10.1111/j.1365246X.2006.02865.x 
Hainzl S., Ogata Y.Hainzl S., Ogata Y. Detecting fluid signals in seismicity data through statistical earthquake modeling[J]. Journal of Geophysical Research:Solid Earth, 2005, 110(B5): B05S07. DOI:10.1029/2004JB003247 
Helmstetter A., Sornette D., Grasso J.R.Helmstetter A., Sornette D., Grasso J.R. Mainshocks are aftershocks of conditional foreshocks:How do foreshock statistical properties emerge from aftershock laws[J]. Journal of Geophysical Research:Solid Earth, 2003, 108(B1): 2046. DOI:10.1029/2002JB001991 
Hou Jinxin, Wang BaoshanHou Jinxin, Wang Baoshan. Temporal evolution of seismicity before and after the 2014 Ludian M_{s}6.5 earthquake[J]. Chinese Journal of Geophysics, 2017, 60(4): 14461456 (in Chinese with English abstract). 
Huang Yilei, Zhou Shiyong, Zhuang JiancangHuang Yilei, Zhou Shiyong, Zhuang Jiancang. Numerical tests on catalogbased methods to estimate magnitude of completeness[J]. Chinese Journal of Geophysics, 2016, 59(4): 13501358 (in Chinese with English abstract). 
Jia Ke, Zhou Shiyong, Zhuang Jiancang, Jiang Changsheng, Guo Yicun, Gao Zhaohui, Gao SheshengJia Ke, Zhou Shiyong, Zhuang Jiancang, Jiang Changsheng, Guo Yicun, Gao Zhaohui, Gao Shesheng. Did the 2008 M_{W}7.9 Wenchuan earthquake trigger the occurrence of the 2017 M_{W}6.5 Jiuzhaigou earthquake in Sichuan, China?[J]. Journal of Geophysical Research:Solid Earth, 2018, 123(4): 29652983. DOI:10.1002/2017JB015165 
Jia Ke, Zhou Shiyong, Zhuang Jiancang, Jiang ChangshengJia Ke, Zhou Shiyong, Zhuang Jiancang, Jiang Changsheng. Possibility of the Independence between the 2013 Lushan Earthquake and the 2008 Wenchuan Earthquake on Longmen Shan Fault, Sichuan, China[J]. Seismological Research Letters, 2014, 85(1): 6067. DOI:10.1785/0220130115 
Jiang Changsheng, Zhuang Jiancang, Long Feng, Han Libo, Guo LujieJiang Changsheng, Zhuang Jiancang, Long Feng, Han Libo, Guo Lujie. Statistical analysis of ETAS parameters in the early stage of the 2013 Lushan M_{S}7.0 earthquake sequence[J]. Acta Seismologica Sinica, 2013a, 35(5): 661669 (in Chinese with English abstract). 
Jiang Changsheng, Wu Zhongliang, Zhuang JiancangJiang Changsheng, Wu Zhongliang, Zhuang Jiancang. ETAS model applied to the EarthquakeSequence Association (ESA) problem:the Tangshan sequence[J]. Chinese Journal of Geophysics, 2013b, 56(9): 29712981 (in Chinese with English abstract). 
Jiang Changsheng, Wu Zhongliang, Han Libo, Guo LujieJiang Changsheng, Wu Zhongliang, Han Libo, Guo Lujie. Effect of cutoff magnitude M_{c} of earthquake catalogues on the early estimation of earthquake sequence parameters with implication for the probabilistic forecast of aftershocks:the 2013 MinxianZhangxian, Gansu, M_{S}6.6 earthquake sequence[J]. Chinese Journal of Geophysics, 2013c, 56(12): 40484057 (in Chinese with English abstract). 
Jiang Changsheng, Wu Zhongliang, Yin Fengling, Guo Lujie, Bi Jinmeng, Wang YawenJiang Changsheng, Wu Zhongliang, Yin Fengling, Guo Lujie, Bi Jinmeng, Wang Yawen. Stability of earlyestimation sequence parameters for continuous forecast of the aftershock rate:a case study of the 2014 Ludian, Yunnan M_{S}6.5 earthquake[J]. Chinese Journal of Geophysics, 2015, 58(11): 41634173 (in Chinese with English abstract). 
Jiang Changsheng, Zhuang Jiancang, Wu Zhongliang, Bi JinmengJiang Changsheng, Zhuang Jiancang, Wu Zhongliang, Bi Jinmeng. Application and comparison of two shortterm probabilistic forecasting models for the 2017 Jiuzhaigou, Sichuan, M_{S}7.0 earthquake[J]. Chinese Journal of Geophysics, 2017, 60(10): 41324144 (in Chinese with English abstract). 
Jiang Haikun, Zheng Jianchang, Wu Qiong, Qu Yanjun, Li YongliJiang Haikun, Zheng Jianchang, Wu Qiong, Qu Yanjun, Li Yongli. Earlier statistical features of ETAS model parameters and their seismological meanings[J]. Chinese Journal of Geophysics, 2007, 50(6): 17781786 (in Chinese with English abstract). 
Jiang Haikun, Yang Maling, Sun Xuejun, Lü Jian, Yan Chunheng, Wu Qiong, Song Jin, Zhao Yong, Huang Guohua, Zhang Hua, Yao Hong, Mu Jianying, Li Jin, Qu JunhaoJiang Haikun, Yang Maling, Sun Xuejun, Lü Jian, Yan Chunheng, Wu Qiong, Song Jin, Zhao Yong, Huang Guohua, Zhang Hua, Yao Hong, Mu Jianying, Li Jin, Qu Junhao. A typical example of locally triggered seismicity in the boundary area of Lingyun and Fengshan following the large rainfall event of June 2010[J]. Chinese Journal of Geophysics, 2011, 54(10): 26062619 (in Chinese with English abstract). 
Jiang Haikun, Song Jin, Wu Qiong, Li Jin, Qu JunhaoJiang Haikun, Song Jin, Wu Qiong, Li Jin, Qu Junhao. Quantitative investigation of fluid triggering on seismicity in the Three Gorge Reservoir area based on ETAS model[J]. Chinese Journal of Geophysics, 2012, 55(7): 23412352 (in Chinese with English abstract). 
Lei Xinglin, Ma Shengli, Wen Xueze, Su Jinrong, Du FangLei Xinglin, Ma Shengli, Wen Xueze, Su Jinrong, Du Fang. Integrated analysis of stress and regional seismicity by surface loadinga case study of Zipingpu reservoir[J]. Seismology and Geology, 2008a, 30(4): 10461064 (in Chinese with English abstract). 
Lei Xinglin, Yu Guozheng, Ma Shengli, Wen Xueze, Wang QiangLei Xinglin, Yu Guozheng, Ma Shengli, Wen Xueze, Wang Qiang. Earthquakes induced by water injection at ~3km depth within the Rongchang gas field, Chongqing, China[J]. Journal of Geophysical Research:Solid Earth, 2008b, 113(B10): B10310. DOI:10.1029/2008JB005604 
Lei Xinglin, Huang Dongjian, Su Jinrong, Jiang Guomao, Wang Xiaolong, Wang Hui, Guo Xin, Fu HongLei Xinglin, Huang Dongjian, Su Jinrong, Jiang Guomao, Wang Xiaolong, Wang Hui, Guo Xin, Fu Hong. Fault reactivation and earthquakes with magnitudes of up to Mw4.7 induced by shalegas hydraulic fracturing in Sichuan Basin, China[J]. Scientific Reports, 2017, 7(1): 7971. DOI:10.1038/s4159801708557y 
Li Zhichao, Huang QinghuaLi Zhichao, Huang Qinghua. Assessment of detectability of the Capitalcircle Seismic Network by using the probabilitybased magnitude of completeness (PMC) method[J]. Chinese Journal of Geophysics, 2014, 57(8): 25842593 (in Chinese with English abstract). 
Liu Lihua, Hao Tianyao, Lü Chuanchuan, You Qingyu, Pan Jun, Wang Fuyun, Xu Ya, Zhao Chunlei, Zhang JianshiLiu Lihua, Hao Tianyao, Lü Chuanchuan, You Qingyu, Pan Jun, Wang Fuyun, Xu Ya, Zhao Chunlei, Zhang Jianshi. Crustal structure of Bohai Sea and adjacent area (North China) from two onshoreoffshore wideangle seismic survey lines[J]. Journal of Asian Earth Sciences, 2015, 98: 457469. DOI:10.1016/j.jseaes.2014.11.034 
Long Feng, Du Fang, Ruan Xiang, Deng Yuanqing, Zhang TiebaoLong Feng, Du Fang, Ruan Xiang, Deng Yuanqing, Zhang Tiebao. Water injection triggered earthquakes in the Zigong mineral wells in ETAS model[J]. Earthquake Research in China, 2010, 26(2): 164171 (in Chinese with English abstract). 
Ma Tengfei. "Repeating Earthquakes" Extracted from Continuous Waveforms: Analysis of the Wenchuan Sequence in Connection to the WIFSD Drilling Site[D]. Beijing: Institute of Geophysics, China Earthquake Administration, 2015 (in Chinese).

Meng Xiaofeng, Yu Xiao, Peng Zhigang, Hong BoMeng Xiaofeng, Yu Xiao, Peng Zhigang, Hong Bo. Detecting earthquakes around Salton Sea following the 2010 M_{W}7.2 El MayorCucapah earthquake using GPU parallel computing[J]. Procedia Computer Science, 2012, 9: 937946. DOI:10.1016/j.procs.2012.04.100 
Narteau C., Byrdina S., Shebalin P., Schorlemmer D.Narteau C., Byrdina S., Shebalin P., Schorlemmer D. Common dependence on stress for the two fundamental laws of statistical seismology[J]. Nature, 2009, 462(7273): 642646. DOI:10.1038/nature08553 
Nyffenegger P., Frohlich C.Nyffenegger P., Frohlich C. Aftershock occurrence rate decay properties for intermediate and deep earthquake sequences[J]. Geophysical Research Letters, 2000, 27(8): 12151218. DOI:10.1029/1999GL010371 
Ogata Y.Ogata Y. Statistical models for earthquake occurrences and residual analysis for point processes[J]. Journal of the American Statistical Association, 1988, 83(401): 927. DOI:10.1080/01621459.1988.10478560 
Ogata Y.Ogata Y. Detection of precursory relative quiescence before great earthquakes through a statistical model[J]. Journal of Geophysical Research:Solid Earth, 1992, 97(B13): 1984519871. DOI:10.1029/92JB00708 
Ogata Y.Ogata Y. Increased probability of large earthquakes near aftershock regions with relative quiescence[J]. Journal of Geophysical Research:Solid Earth, 2001, 106(B5): 87298744. DOI:10.1029/2000JB900400 
Omori F.Omori F. On the aftershocks of earthquake[J]. Journal of the College of Science, Imperial University of Tokyo, 1894, 7: 11200. 
Peng Yajun, Zhou Shiyong, Zhuang Jiancang, Shi JiaPeng Yajun, Zhou Shiyong, Zhuang Jiancang, Shi Jia. An approach to detect the abnormal seismicity increase in Southwestern China triggered coseismically by 2004 Sumatra M_{W}9.2 earthquake[J]. Geophysical Journal International, 2012, 189(3): 17341740. DOI:10.1111/j.1365246X.2012.05456.x 
Peng Zhigang, Zhao PengPeng Zhigang, Zhao Peng. Migration of early aftershocks following the 2004 Parkfield earthquake[J]. Nature Geoscience, 2009, 2(12): 877881. DOI:10.1038/ngeo697 
Schaff D.P.Schaff D.P. Semiempirical statistics of correlationdetector performance[J]. Bulletin of the Seismological Society of America, 2008, 98(3): 14951507. DOI:10.1785/0120060263 
Shelly D.R., Beroza G.C., Ide S.Shelly D.R., Beroza G.C., Ide S. Nonvolcanic tremor and lowfrequency earthquake swarms[J]. Nature, 2007, 446(7133): 305307. DOI:10.1038/nature05666 
Song Jin, Jiang HaikunSong Jin, Jiang Haikun. A review on decay and generation of aftershock activity[J]. Seismology and Geology, 2009, 31(3): 559571 (in Chinese with English abstract). 
Su LuanshengSu Luansheng. Relationship between the fractal structure variation in Changdao area and the arealsei smicity and shortterm and impending precursor feature[J]. Journal of Seismological Research, 1993, 16(1): 3340 (in Chinese with English abstract). 
Tan Yipei, Deng Li, Cao Jingquan, Shan LianjunTan Yipei, Deng Li, Cao Jingquan, Shan Lianjun. Seismological mechanism analysis of 2015 Luanxian swarm, Hebei province[J]. Chinese Journal of Geophysics, 2016, 59(11): 41134125 (in Chinese with English abstract). 
Utsu T.Utsu T. A statistical study on the occurrence of aftershocks[J]. Geophysical Magazine, 1961, 30: 526605. 
Wang Peng, Hou Jinxin, Wu PengWang Peng, Hou Jinxin, Wu Peng. Temporal evolution of the seismicity of the 2017 Jiuzhaigou M_{S}7.0 earthquake sequence[J]. Earthquake Research in China, 2017, 33(4): 453462 (in Chinese with English abstract). 
Wang Zhicai, Deng Qidong, Chao Hongtai, Du Xiansong, Shi Ronghui, Sun Zhaomin, Xiao Lanxi, Min Wei, Ling HongWang Zhicai, Deng Qidong, Chao Hongtai, Du Xiansong, Shi Ronghui, Sun Zhaomin, Xiao Lanxi, Min Wei, Ling Hong. Shallowdepth sonic reflection profiling studies on the active PenglaiWeihai fault zone offshore of the northern Shandong peninsula[J]. Chinese Journal of Geophysics, 2006, 49(4): 10921101 (in Chinese with English abstract). 
Wiemer S., Wyss M.Wiemer S., Wyss M. Minimum magnitude of completeness in earthquake catalogs:Examples from Alaska, the Western United States, and Japan[J]. Bulletin of the Seismological Society of America, 2000, 90(4): 859869. DOI:10.1785/0119990114 
Xu Jie, Song Changqing, Chu QuanzhiXu Jie, Song Changqing, Chu Quanzhi. Preliminary study on the Seismotectonic characters of the Zhangjiakou Penglai fault zone[J]. Seismology and Geology, 1998, 20(2): 146154 (in Chinese with English abstract). 
Yang Hongfeng, Zhu Lupei, Chu RishengYang Hongfeng, Zhu Lupei, Chu Risheng. Faultplane determination of the 18 April 2008 Mount Carmel, Illinois, earthquake by detecting and relocating aftershocks[J]. Bulletin of the Seismological Society of America, 2009, 99(6): 34133420. DOI:10.1785/0120090038 
Zhang Miao, Wen LianxingZhang Miao, Wen Lianxing. Seismological evidence for a lowyield nuclear test on 12 May 2010 in North Korea[J]. Seismological Research Letters, 2015, 86(1): 138145. DOI:10.1785/02201401170 
Zheng Jianchang, Xu Changping, Li Dongmei, Zhou CuiyingZheng Jianchang, Xu Changping, Li Dongmei, Zhou Cuiying. Recent change of stress field in North China derived from station composite focal mechanisms[J]. Recent Developments in World Seismology, 2018(4): 4857 (in Chinese with English abstract). 
Zhuang JiancangZhuang Jiancang. Statistical modeling of seismicity patterns before and after the 1990 Oct 5 Cape Palliser earthquake, New Zealand[J]. New Zealand Journal of Geology and Geophysics, 2000, 43(3): 447460. DOI:10.1080/00288306.2000.9514901 
Zhuang Jiancang, Ogata Y., Wang TingZhuang Jiancang, Ogata Y., Wang Ting. Data completeness of the Kumamoto earthquake sequence in the JMA catalog and its influence on the estimation of the ETAS parameters[J]. Earth, Planets and Space, 2017, 69(1): 36. 