Anton Zaicenco
Weir-Jones Engineering Consultants, Ltd.
2040 West 10th Avenue, Vancouver, British Columbia
Tel: 604-732-8821, Fax: 604-732-4801

Sharlie Huffman
Bridge seismic engineer, Ministry of
Transportation, 940 Blanshard St., Victoria, British Columbia, V8W 3E6, Canada
Tel: 250-356-0761, Fax: 250-387-7735

Iain Weir-Jones
President, Weir-Jones Engineering Consultants, Ltd. 2040 West 10th Avenue, Vancouver, British Columbia
Tel: 604-732-8821, Fax: 604-732-4801

Design of an on-site seismic early warning system for critical facilities, based on P-wave detection, requires knowledge of local seismicity, analysis of strong-motion records, and local geological settings. The scope of current paper is to present an efficient and reliable methodology for P-wave detection and discrimination. It is tolerant to the environmental, traffic and other types of noise at the location of the critical facility.

The basic idea that is exploited in the proposed algorithm is the degree of P-wave polarization, which is obtained with Karhunen- Loève transform of the orthogonal triaxial records. A system consisting of borehole sensor array that has been designed based on these ideas is described and some preliminary results are provided. This system has been already implemented to provide earthquake early warning for a tunnel on a major highway in British Columbia, Canada.

Due to the increased interest in the earthquake early warning systems (EEWS) 14th World Conference on Earthquake Engineering in Beijing China had a special session on the topic. 11 papers were presented at this session, mainly by authors from Japan. In what follows, a review of several papers will be presented to provide a snapshot of major topics.

Kanda et al. (2008) were interested in accurate prediction of the seismic intensity using different source locations and site- specific empirical method. On-site warning system has been developed for sites with small epicentral distance in addition to JMA system to pick P-waves and predict seismic intensity. The prototype of the system has been validated using vibration tests, numerical simulations, and observations. Authors provide empirical relationships between RMS value of P-wave velocity VSD, maximum velocity Vmax of 1-3 sec time series after P-wave and measured intensity:

I = 2 log(VSD) + α (1)

I = 2 log(Vmax) + β (2)

where α=6.12 and β=5.14 for 2 seconds of P-wave. Yet, standard deviation was found relatively large.

Nakamura (2008) gives an overview of the P-wave early warning systems. The author introduces the notion of “On-Site Alarm”, which is the alarm based on the observation at the side of the objects to be warned, and emphasizes that it is more important that the network alarm. Such capabilities of EEWS as estimation of azimuth, predominant frequency for magnitude evaluation, and vertical to horizontal ratio for discrimination between P- and S-waves are discussed.

UrEDAS system described by the author can detect earthquakes in P-wave triggering and then estimates magnitude, epicentral and hypocentral distance, depth and back azimuth. Tolerance of EEWS to electrical thunder noise and small pre-shocks has been emphasized. Author concludes that “on-site alarm ... is useful for anywhere, even epicentral area”.

Saita et al. (2008) discuss the applications of EEWS. On-site alarm is viewed as an extension of network alarm system important to protect rescue teams. FREQL, Fast Response Equipment against Quake Load, is a new generation of EEWS, which is an all-in-one seismometer with the power unit. It is capable of detecting P-wave from a single station and issues an alarm in 0.2 seconds. EEWS by JMA that includes a large number of seismometers nationwide is pseudo real-time system because the system stores the waveform of the first 2 seconds and calculates parameters. If the system does not converge to the acceptable results, it uses 3, 5, 7, 10 seconds repeatedly. Therefore, the authors recommend that it is better to have an on-site EEWS at each facility, while the network information must be used as backup.

Motosaka (2008) discusses application of EEWS based on the nationwide earthquake observation networks in Japan for disaster prevention in schools. Availability of several EEWS in Japan, Mexico, Taiwan, South California and Istanbul is mentioned. The author demonstrates the application of EEWS to schools in Miyagi prefecture, Japan.

Shanyou and Jindong (2008) present a methodology for predicting magnitude from the first few seconds of the P- wave, which is important in EEWS applications. They propose to estimate the value of magnitude from a single station as
function of amplitude Pmax, predominant period τc and epicentral distance based on a linear regression model:

M = a log Pmax + b log τc + c log + d (3)

The waveforms used in the study were collected from 17 earthquakes occurring between 1999 and 2006 from NSMP by USGS. Automatic P-wave picker has been employed (Allen (1978)) to detect the P-wave arrivals. The proposed technique has less uncertainty in comparison with current methods based
on τc or Pmax only.

Kuyuk and Motosaka (2008) identify another important application of EEWS -- its ability to provide information to reduce seismic response. Forecasting Fourier Amplitude Spectrum is important for structural control systems. The authors developed a regional warning system in Miyagi Prefecture against Miyagi-ken offshore earthquakes and integrated it with JMA national Japan EEWS.

Fujinawa et al. (2008) describe the EEWS built around JMA network in cooperation with NIED. Real-time data collection allows distributing warning to various customers. Since October 1, 2007 the information has been put to general use via broadcasting companies (TV, radio, mobile phones), and contracted use by private companies using special application instruments. Meguro (2008) discusses the problems of EEWS based on JMA network related to proper dissemination of the information to the public. In view of the author, the EEWS should finally become a driving force to promote disaster mitigation countermeasures. Regarding the technical aspects of the EEWS, it is emphasized that accurate determination of the magnitude from initial phase of the rupture is not possible. There is a trade-off between available time and accuracy.

Alarm network is described by Masaki et al. (2008). It is based on JMA EEWS and targets central area of Japan -- Mikawa area. The system is useful for earthquake disaster reduction by performing evacuation of workers and stopping machines before strong shaking arrival.


The equation of the seismic wave propagation in the elastic medium under assumption that Lame parameters λ and µ are
constant is:

ρü = (λ+µ)(∇⋅u) + µ∇×(∇×u) (4)

where ρ is the density of material and u is the displacement vector. This is a 3D homogeneous vector wave equation for a uniform, isotropic linear medium.

There are two propagation modes: P- (associated with div: ∇⋅u and governed by a scalar wave equation) and S-waves (associated with curl: ∇×u, and governed by vector wave equation), which are pressure and shear waves, Fig.1.

Fig.1 Example of seismic wave propagation in the inhomogeneous elastic medium from a point source: numerical solution using spectral method, Zaicenco and Alkaz (2008)
Fig.1 Example of seismic wave propagation in the inhomogeneous elastic medium from a point source: numerical solution using spectral method, Zaicenco and Alkaz (2008)

P-waves are associated with the propagation of deformations in terms of dilatation:

∇⋅u = u1/x1 + u2/x2 + u3/x3 (5)

This is a scalar value and is related to trace of strain tensor
tr(ε), which is invariant under coordinate transformation. By
taking the divergence of Eq.(4) a wave equation for P waves is
obtained. The propagation speed of P wave is:

vP2 = (λ+2µ)/ρ (6) S-waves are associated with the propagation of deformations
related to change of shape/rotation. They are described by off- diagonal terms of antisymmetric rotation tensor:

ξij = 1/2 (ui/xj - uj/xi) (7)

By taking the curl of Eq.(4) a wave equation for S-waves is
obtained. The propagation speed of S-wave is:

2 = µ/ρ (8)

The difference in propagation speeds given by Eqs.(6) and (7)
is usually exploited in early warning systems.

Along with speed of propagation, polarization is another key factor allowing one to distinguish between P- and S-waves. Since P-wave is related to change of volume, it propagates perpendicular to the wave front, and is linearly polarized for a homogeneous isotropic linear medium. S-wave, being related to rotational deformations, is polarized in the plane parallel to the wave front for a homogeneous isotropic linear medium. An example of orthogonality in polarization of P- and S-waves is demonstrated for 1989 Loma Prieta earthquake record in Fig.2. Axes of polarization of 3D ground motion for both waves are visualized with the help of ellipsoids.


b,c, Fig.2 (a) Three components of 1989 Loma Prieta(b),(c)

Fig.2 (a) Three components of 1989 Loma Prieta
earthquake record measured by the instrument aligned parallel to the fault, (b) time-frequency distributions, and (b) scaled polarization ellipses and their axes for P- and S-waves


The method based on polarization analysis has been used for source location by Christoffersson et al. (1988), Ruud et al. (1988). They applied probability estimator to the energy of a temporal window corresponding to P-wave onset. Jackson et al. (1991), based on several previous papers by Flinn (1965), Montalbetti and Kanasewich (1970), Esmersoy (1984), have achieved a polarization analysis using singular value decomposition for a single-station triaxial recording.

When the location of the source is known, a simple rotation of the coordinate system can be used to identify SH-, SV-wave axes. This method has been used by Cassidy and Rogers (1999) to study the site response spectra, which were computed using a 15 s data window, beginning about 2 s before S-wave arrival. The authors analyzed data from 1996 Duvall and 1997 Georgia Strait earthquakes and concluded that strong amplification, up to 11 times, relative to the bedrock is present at frequencies 1-3 Hz.

When the source direction is not known, polarization analysis of the seismic wave can be done with the Karhunen-Loéve transform (KLT). KLT is treating stochastic process as a combination of the orthogonal functions, and the expansion basis depends on the process itself. Polarization of the recorded 3 orthogonal components of ground motion is computed with discrete KLT, which is based on the eigenvalue decomposition of the covariance matrix C:

C = TΛT-1 (9)

Eigenvector matrix T allows to uncorrelate initial orthogonal components. Coordinate rotation aligns transformed axis with the directions of maximum variance. Since components with larger variance correspond to interesting dynamics and lower ones represent noise, signal-to-noise ratio (SNR) is critical in such type of analysis.


The EEWS system includes 2 borehole strings with 4 triaxial sensors each. Total 8 sensors provide 24 channels of real-time data. The basic steps of the algorithm are as follows:
1. Band pass filter is applied. Non-recursive FIR band pass filter is applied to each channel yi:

vij = Σk ckyi-k,j (10)

Coefficients ck are computed outside of the real-time computations.
2. Cramer-Leadbetter (1967) envelope function qij is computed, Fig.3.
3. Peak values of qij are compared with threshold Q.
4. STA/LTA (Allen (1978)) on the envelope function produces
rij . If rij > R, event is detected at time ti* for each channel of 8 sensors, Fig.4.

5. Extract portion of signal Pk×3 that includes the first-motion

P-wave from each triaxial record, which corresponds to time frame tp = [ti* - t1; ti* + t2].

6. Apply discrete KLT to covariance matrix of Pk×3:

Σ = Cov(Pi) = QΛQT (11)

Compute normalized variances:
wi = Var(Yi)= max(Var(Yi)) (5)

and decide if polarization ellipsoid is strongly stretched in one dimension (linear polarization).
7. If q out of 8 sensors show linear polarization and peak levels on q sensors out of 8 are greater then triggering level, detected event is confirmed as P-wave.



Fig.3 (a) Raw and filtered data, (b) Cramer-Leadbetter
enveloped function for filtered data

Fig.4 Event location with STA/LTA ratio

The records of ambient noise obtained by the sensors of the EEWS, Fig.5, indicate the presence of two dominant peaks around 3 Hz and 12 Hz. They seem to originate from the highway R/C structure excited by the traffic, and response of thick alluvial soft soil deposits, which are characteristic for the location of the early warning system.

Fig.5 Normalized PSD of the records of ambient noise obtained by EEWS sensors: (a) vertical components, (b) horizontal components
Fig.5 Normalized PSD of the records of ambient noise obtained by EEWS sensors: (a) vertical components, (b) horizontal components

Horizontal axis calibration and system tests were performed by applying a vertical short-duration force on the free surface, Fig.6. The 40 J pulse was generated by hitting a metallic plate with a 20 kg rod at a distance of 6 m from the borehole collar, which corresponds to magnitude ML=-2. In spite of the fact that the site is contaminated by traffic noise from a highway, the signal was picked by the bottom geophones at depth 40m. The soil profile at shallow depths is represented by 20 m of
loose sand at the top followed by sandy silt, in both cases saturated. The signals with frequency around 50 Hz recorded by the borehole sensors are provided in Fig.7.

Fig.6 In situ free surface P-wave pulse test
Fig.6 In situ free surface P-wave pulse test

Fig.7 P-wave recorded by 3C sensors of the EEWS from the pulse test



Projections of the polarized P-wave component recorded by four 3C sensors on the horizontal plane are shown in Fig.8. The signal picked by the surface sensor exhibits higher elliptical polarization due to the strong contribution of the surface waves. The orientation of the sensors in horizontal plane is computed from polarized components by fitting a straight line with the equation y=kx using least-squares method. Coefficient k for each sensor is given in Fig.8.

Fig.8 Horizontal projections of polarized P-wave component recorded by four 3C sensors of the EEWS from the pulse test
Fig.8 Horizontal projections of polarized P-wave component recorded by four 3C sensors of the EEWS from the pulse test

The ratios of vector lengths of peak amplitudes ||max (ax,y,z)|| for sensors at depths 4 m, 20 m and 40 m are, correspondingly, 1:6:9, while distance to the source ratios are ≈1:2:4. Higher than 1/R attenuation is partly explained by reflection of body wave from the interface between sandy clay and clay at depth 15m, and presence of water table close to this interface.

Knowledge of the sensor horizontal axis orientation allows to compute the direction of the incoming seismic wave, while its amplitude can serve as an estimate of the magnitude of the event and its destruction capacity.

Several multi-channel 3C records available from other locations were used to test the algorithm and confirm the applied methodology. Records 1 through 5 contain P-waves obtained from various sources. Results are summarized in Table 1, and Figure 9.

Table 1 Test of EEWS algorithm using 3C records

table earth quake warning system algorithm

Fig.9 Test multi-channel records: (a,b,c) underground explosion, (d) calibration pulse test, (e) heavy weight drop test(a), (b), (c), (d), (e)
Fig.9 Test multi-channel records: (a,b,c) underground explosion, (d) calibration pulse test, (e) heavy weight drop test

Several records of the Nisqually earthquake (Feb 2001) were kindly provided by GSC Pacific (Onur et al., 2001). Fig.10 shows time-histories of the ground motion record from the station KID, and STFT analysis of these components. These examples demonstrate presence of higher frequency components in the P phase in comparison with the S phase. Polarization analysis of P- and S-phase of the record showed orthogonality of these components and in-plane polarization of the S-wave. The algorithm designed for the EEWS was able to identify successfully the P-wave from this record.


Fig.10 3C record of Nisqually earthquake. Site KID, southern shore of North Fraser Arm: (a) longitudinal, transversal and vertical components, (b) STFT analysis of these components (b)
Fig.10 3C record of Nisqually earthquake. Site KID, southern shore of North Fraser Arm: (a) longitudinal, transversal and vertical components, (b) STFT analysis of these components


The following describes the possible expansion of the on-site EEWS to include earthquake parameters estimation and its ability to predict the severity of the ground motion. Least-squares method applied to the arrival time data is a standard technique for a hypocentral location problem. Adding additional information -- the azimuth and emersion angle of the seismic ray coupled with ray-tracing technique -- allows to calculate more accurately coordinates of the hypocenter, Alessandrini et al. 1994. These additional parameters are computed with the polarization analysis, which is part of the proposed EEWS.

The maximum predominant period of the first 1-4 s of the P- wave τc is widely used as a key parameter to estimate the
magnitude of the seismic event (Allen and Kanamori (2003), Wu and Kanamori (2005)):

τc = 2π/r1/2, r = u2(t)dt / u2(t)dt (6)

where u(t) is the ground-motion displacement. Expression of τc seems to be an obvious extend of the Parseval’s theorem, frequency-domain derivation u=ω2u and an assumption that P-wave component is a narrow-band signal.

It has been demonstrated that the peak initial displacement amplitude, Pd, correlates well with PGV, and its attenuation functions are provided for southern California by Wu and Zhao (2006). Such relationship can be used to estimate the earthquake magnitude.

In order to expand the capabilities of the proposed EEWS, authors consider expanding the on-site system to a regional
level and incorporating parameters τc and Pd into the existing algorithm.

The on-site EEWS for a tunnel on a major highway in British Columbia, Canada has been presented in the current paper. The algorithm for detection and discrimination of P-wave is based on the polarization analysis of band-passed triaxial record obtained in real-time by geophone sensors installed in two boreholes. The proposed methodology has been tested using available data of blast records and strong-motion free- field records. Proposals for future system expansion include magnitude estimation and hypocentral location using parameters τc and Pd.

The records of strong earthquakes were kindly provided by Prof. Carlos Ventura, Department of Civil Engineering, University of British Columbia, and by Dr.Andreas Rosenberger, Geological Survey of Canada, Pacific Geoscience Centre. Support provided by BC Ministry of Transportation is highly appreciated.

Alessandrini, B Cattaneo, M Demartin, M Gasperini, M & Lanza, V 1994, ‘A Simple P-wave polarization Analysis: its Application to Earthquake Location’, Annali di Geofisica, Vol. XXXVII, N.5, Sept.

Allen, RM 1978, ‘Automatic Earthquake Recognition and
Timing from Single Traces’. Bull. Seism. Soc. Am. 68:5, pp.

Allen, RM & Kanamori, H 2003, ‘The potential for earthquake early warning in Southern California’, Science,

300, pp. 786–789

Cassidy, JF & Rogers, GC 1999, ‘Seismic Site Response in the Greater Vancouver, British Columbia, area: Spectral Ratios from Moderate Earthquakes’, Can. Geotech. J., 36, pp.195-209

Christoffersson, A Husebye, ES & Ingate, SF 1988, ‘Wavefield Decomposition Using ML-probabilities in Modelling Single-site Three-component Records’, Geophys. J. Int., 93, pp.197-213

Cramer, H & Leadbetter, MR 1967, ‘Stationary and related stochastic processes’, Wiley, New York

Esmersoy, C 1984, ‘Polarization Analysis, Orientation and Velocity Estimation in Three Component VSP’, Vertical Seismic Profiling - Part B: Advanced Concepts, eds.: Toksoz, MN & Stewart, RR (Geophysical Press)

Flinn, EA 1965, ‘Signal Analysis Using Rectilinearity and Direction of Particle Motion’, Proc. IEEE, 53, pp. 1874-1876

Fujinawa, Y Rokugo, Y Noda, Y Mizui, Y Kobayashi, M & Mizutani, E 2008, ‘Efforts of Earthquake Disaster Mitigation Using Earthquake Early Warning in Japan’, 14th WCEE, Beijing, China, Oct 12-17

Jackson, GM Mason, IM & Greechalgh, SA 1991, ‘Principal Component Transforms of Triaxial Recordings by Singular Value Decomposition’, Geophysics, 56(4), pp.528-533

Kanda, K Nasu, T Miyamura, M & Koide, E 2008,
‘Development of site-specific earthquake early warning system for hazard mitigation’, 14th WCEE, Beijing, China, Oct 12-17

Kuyuk, HS & Motosaka, M 2008, ‘Spectral Forecasting of Earthquake Ground Motion Using Regional and National Earthquake Early Warning Systems for Advanced Engineering Application Against Approaching Miyagi-Ken Oki Earthquakes’, 14th WCEE, Beijing, China, Oct 12-17

Masaki, K Kurahashi, S & Irikura, K 2008, ‘Development of Alarm Network Using Earthquake Early Warning’, 14th WCEE, Beijing, China, Oct 12-17

Meguro, K 2008, ‘Strategy for Taking Full Advantage of Earthquake Early Warning System for Earthquake Disaster Reduction’, 14th WCEE, Beijing, China, Oct 12-17

Montalbetti, JR & Kanasewich, ER 1970, ‘Enhancement of Teleseismic Body Phases with a Polarization Filter’, Geophys. J. R. Astron. Soc., 21, pp. 119-129

Motosaka, M 2008, ‘Application of Earthquake Early Warning Systems for Disaster Prevention in Schools’, 14th WCEE, Beijing, China, Oct 12-17

Nakamura, Y 2008, ‘First Actual P-wave Alarm Systems and Examples of Disaster Prevention by Them’, 14th WCEE, Beijing, China, Oct 12-17

Onur, T Ventura, CE & Tao, S 2001, Strong Motion Records from the February 28, 2001 Nisqually Earthquake Recorded by Three Structural Arrays in BC, Report EQ 01-02, UBC Earthquake Engineering Research Lab

Ruud, BO Husebye, ES Ingate, SF & Christoffersson, A 1988, ‘Event Location at Any Distance Using Seismic Data from a Single, Three-component Station’, Bull. Seismol. Soc. Am., 78, pp. 308-325

Saita, J Sato, T Nakamura, Y 2008, ‘What is the useful application of the earthquake early warning system?’ 14th WCEE, Beijing, China, Oct 12-17

Shanyou, L & Jindong, S 2008, ‘A New Magnitude Estimation Method Based on Predominant Period and Peak Amplitude’, 14th WCEE, Beijing, China, Oct 12-17

Wu, YM & Kanamori, H 2005, ‘Experiment on an onsite early warning method for the Taiwan early warning system’, Bull. Seism. Soc. Am., 95, pp. 347–353.

Wu, YM & Zhao, L 2006, ‘Magnitude estimation using the first three seconds P-wave amplitude in earthquake early warning’, Geophys. Res. Lett., 33, L16312, doi:10.1029/2006GL026871.

Zaicenco, A & Alkaz, V 2008, ‘Numerical Solution Of An Elastic Wave Equation Using The Spectral Method’, In: Harmonization of Seismic Hazard in Vrancea Zone, Editors: Zaicenco, A Craifaleanu, I Paskaleva, I, Springer, NATO Science for Peace and Security Series C: Env. Security, pp.319-327