IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, VOL. 1, NO. 3, JULY 2004 147 Automatic P Phase Picking Using Maximum Kurtosis and -Statistics Criteria Christos D. Saragiotis, Student Member, IEEE, Leontios J. Hadjileontiadis, Member, IEEE, Ioannis T. Rekanos, Member, IEEE, and Stavros M. Panas, Member, IEEE Abstract—The identification problem of P seismic phase onset has been addressed, based on the maximum-kurtosis assumption, -statistics, and the Chebyshev inequality. Depending on two sta- tistical decision criteria, the proposed approach provides either a single P onset peak or an interval in which the P onset exists, both with a confidence percentage. Results on real seismic data eval- uated using a performance index justify the contribution of the proposed method toward an accurate and fully automated P onset identification. Index Terms—Higher order statistics (HOS), kurtosis, -statis- tics, P phase identification, seismic signal processing. I. INTRODUCTION T HE DETECTION of the P seismic phase arrival (P onset) has to be as precise as possible in order to facilitate the extraction of useful conclusions on epicentral and hypocentral coordinates of the seismic source, the interior structure of the Earth, seismic hazard assessment, seismological networks oper- ation, early warning systems, etc., and is traditionally performed visually by human experts, the analysts. However, it is usually the case that the seismic signal is embedded in a sufficiently strong noise, such that the visual detection of the P phase ar- rival can be impossible or even, when possible, different ana- lysts will identify different P onsets or the same analyst will identify different P onsets at different times. The deviation in the P onset identification may not be large; however, depending on the distance of the seismic source, even small deviations (of the order of tenths of a second) of the P onset time may produce a deviation of kilometers in the determination of the hypocen- tral or epicentral coordinates. The P onset identification is, thus, a difficult task, due to the variation of the temporal and spectral characteristics of the seismic signal and noise, as well as the ex- istence of later arrivals immediately following the P phase, the huge data volume, and the subjectivity of the analysts. Therefore, a fast and objective detector for the P onset identification is of great significance for the seismologists. Such detectors involve energy ratio criteria, i.e., variations of the ratio of short-term average (STA) to long-term average Manuscript received January 13, 2004; revised March 18, 2004. The work of C. D. Saragiotis was supported by the Hellenic Scholarships State Institute (I.K.Y.) under Grant 1998-2002. C. D. Saragiotis, L. J. Hadjileontiadis, and S. M. Panas are with the Department of Electrical and Computer Engineering, Aristotle University of Thessaloniki, GR 54124 Thessaloniki, Greece (e-mail: csar@auth.gr; leontios@auth.gr; panas@psyche.ee.auth.gr). I. T. Rekanos is with the Department of Informatics and Communications, Technological and Educational Institute of Serres, GR 62124 Serres, Greece (e-mail: rekanos@teiser.gr). Digital Object Identifier 10.1109/LGRS.2004.828915 (LTA) or STA/LTA, the seismic wave polarity assumption, artificial neural networks, maximum-likelihood methods, fuzzy logic theory, the wavelet transform, and higher order statistics (HOS) (a collection of references for those methods can be found in [1]). Among these methods, the most promising are Allen’s energy criterion [2], Earle’s energy criterion [3], and HOS-based methods [1], [4], [5], with the latter looking the most efficient ones, since they are the least affected by noise and are signal-bandwidth independent [1]. In [1], there is also a comparison of a HOS-based algorithm, namely the phase arrival identification-kurtosis (PAI-K) algorithm, with the Allen algorithm [2] in an artificially noisy environment, in which the PAI-K algorithm outperforms the Allen one. The PAI-K algorithm uses a sliding time window whose length has to be predefined. The optimal length of the sliding time window depends on the sampling frequency and the seismic characteristics, while its value is crucial for the calculation of the kurtosis estimate and for the accurate detection of the P onset. In this letter, an alternative approach is proposed, which circumvents the determination of an optimal time-window length. In addition, the proposed method focuses not only on the accurate estimation of the P onset, but on the certainty of this estimation as well. This is achieved by employing -statistics criteria. Furthermore, an index function is proposed that quanti- fies the performance of the estimation procedure. The proposed algorithm is applied to real seismic traces. Its P onset estimates are compared to those derived by the Allen’s algorithm and the P onsets picked up by expert analysts. II. RATIONALE AND MATHEMATICAL BACKGROUND A. HOS Considerations The kurtosis of a fourth-order stationary, zero-mean, sto- chastic process is defined [6] as the normalized fourth-order cumulant for zero time lag or (1) where denotes the expected value operator. Theoretically, the kurtosis of a Gaussian signal is zero; hence, one would only have to compare its kurtosis value to zero, in order to detect the deviation or not from Gaussianity. However, this is not the case at all, as the kurtosis parameter is merely esti- mated using a finite-length time-window, and its value depends on the length of the sample. Thus, the kurtosis estimate is autho- rized to exist in a confidence interval, which is conditioned by the probability properties of the estimator. Asymptotically, this 1545-598X/04$20.00 © 2004 IEEE