* Digital Object Identifier 10.1109/MSP.2012.2182949*
* Date of publication: 9 April 2012*

**T**

he retrieval of
information
con-v e y e d i n d a t a
recorded by
seis-mic arrays plays a
key role in seismology and
geophysical exploration.
Accurate localization and
reli-able detection of seismic
events are major tasks in
seis-mic monitoring systems. The
nonstationarity, low
signal-to-noise ratio (SNR), and weak
signal coherence of seismic
data remain challenging
issues for signal processing
algorithms. The maximum
likelihood (ML) approach that
performs well in such critical
conditions is one of the best solutions for simultaneous detec-tion and localizadetec-tion of seismic events. This article will discuss the methodology of ML for estimation and detection of seismic data and its extension to geoacoustic model selection.

**INTRODUCTION**

Since 1996, the overwhelming majority of countries have signed the Comprehensive Nuclear Test Ban Treaty, which bans

nuclear testing underground, in oceans, and in the atmo-sphere. One important aspect of this treaty is the seismic monitoring regime. The abili-ty to accurately detect and locate worldwide all seismic events above a local magni-tude of four within a circle of about 18 km is a formidable scientific and technical chal-lenge [1]. Recently, the 2011 magnitude nine earthquake in Tohoku (also known as the Great East Japan Earthquake) triggered extremely devastat-ing tsunami waves that struck Japan’s east coast. In addition to tremendous loss of life and destruction, nuclear accidents caused by the tsunami at the Fukushima power plant was a serious health and safety concern to countries around the globe. Hence, the task of reliable, effi-cient, and accurate detection of seismic events is of great impor-tance to early warning systems and disaster mitigation.

To extract information contained in propagating waves, one has to automatically process signals measured by a network of seismic arrays and single stations with robust algorithms [2]– [4]. In this article, we address the aspect of signal detection and parameter extraction from a single regional seismic array [5]. In particular, we will investigate a statistically motivated ML

**[**

### Pei-Jung Chung and Johann F. Böhme

**]**

**[**

### Estimation, detection, and exploration of seismic events

**]**

IMAGE COURTESY OF U.S. DEPARTMENT OF COMMERCE/NOAA/NESDIS/NATIONAL GEOPHYSICAL DATA CENTER

approach for seismic data analysis. Among existing array pro-cessing methods, the ML approach is characterized by optimal statistical properties and robustness against critical scenarios involving low SNRs, closely located sources, and signals coher-ence. In addition, the ML approach is applicable to both narrow band and broadband data. Unlike subspace methods, which require an additional focusing step in the broadband case, the ML approach combines broadband data in a statistically justified manner by exploiting asymptotic normality and independence of Fourier transformed data [5]. Based on the likelihood princi-ple, a multiple test procedure can be formulated to detect seis-mic events with high accuracy and good time resolution [6]. The detection capability is of particular importance when multi-ple events happen within short time intervals. In [7], the experi-mental comparison with subspace methods [8] and the beamforming-based f-k analysis [9] shows that the ML approach provides the best estimation accuracy and resolution in space and time.

In the following, we first describe the array data model, their statistical properties in frequency domain, and the underlying parametric model. Next, the ML estimates for wave parameters are derived and the multiple hypothesis testing for signal detec-tion is addressed. The applicadetec-tion of the ML approach to geo-acoustic model selection is briefly discussed. Finally we show experimental results obtained by processing data recorded by the German Experimental Seismic System (GERESS) array and make a comparison between the proposed ML method and f–k analysis in terms of their performance.

**DATA MODEL**

Seismic events generate waves propagating through the interior
of the earth and along its surface layer. For propagation distance
much larger than the array aperture, the wave fronts lie
approx-imately on a flat plane perpendicular to the propagation
*direc-tion. The wave field is observed by sensors located at N distinct *
*positions pi, i5 1, c, N. Suppose M wave types are present. *

The array outputs are sampled temporally by a properly chosen sampling frequency and short-time Fourier-transformed

*Xl*1v2 5 1

*"T* a

*T*21

*t*50

*wl1t2x1t2e2jvt, l5 1, c, L, (1)*

*where wl1t2 s are orthonormal window functions [10], [11]. *

Assuming the plane wave model, the steering vector associated
*with the ith wave is given by *

*di*1v2 5 3e*2jvji*

*T _{p}*

1_{,}c,e*2jvjiTpN*4*T*_{, (2)}
*where the superscript T defines matrix transpose. The *
expres-sion (2) considers only sensitivity in the vertical direction.
More details on polarization sensitive arrays can be found in
[12]. The slowness vector j*i* is related to the propagation

*velocity Vi*, azimuth a*i*, and elevation f*i* through the following

equation:
j*i*5

1

*Vi* 3cos f

*i *sin a*i*, cos f*i *cos a*i*, sin f*i*4*T*. (3)

In the frequency domain, array outputs can be approximately described through a nonlinear regression model

*Xl1v2 5 H1v,q2 Sl1v2 1 Ul*1v2, _{(4)}

*where the transfer function H1v, q2 5 3d*1*1v2 cdM*1v24

*con-sists of M steering vectors. The nonlinear wave parameters are *
summarized in q 5 3j1,c, j*M4. Sl1v2 and Ul*1v2 denote the

Fourier-transformed signal and noise vectors, respectively. The advantages of using orthonormal windows, for example, prolate spheroidal sequences suggested in [10] and [11] in (1) include the following:

1) The energy in the observation interval can be concentrated in a prespecified frequency band.

2) The orthogonality of window functions decouples the
*dependence of Xl1v2, 1l 5 1, c, L2 and leads to a reduced *

variance in estimation of the array spectral matrix.
The ML approach relies on proper probabilistic modeling of
the data. According to asymptotic theory of Fourier transform
*[13] [14], the array output Xl*1v2 is characterized by the following

sta-tistical features:

■ *For large T, Xl*1v2 is complex normally distributed.
■ *When the signal Sl*1v2 is considered as deterministic,

*the distribution of Xl*1v2 is completely specified by the

*mean H1v,q2Sl1v2 and covariance matrix C*

*Ul*1v2, denoted
*by Xl*1v2 |N*c1H1v, q2Sl1v2, C*

*Ul*1v22.

■ For different frequency bins, v*i*2 v*j, Xl*1v*i2 and Xl*1v*j*2

are mutually independent.

In contrary to analysis based on time domain observation

*x1t2, these properties do not assume normal distribution in *

time domain. In fact, they only require regularity conditions on
*the moments of x1t2. The noise covariance matrix is defined as *
* CUl1v2 5 E3Ul1v2Ul*1v2

*H4, where the superscript H denotes*Hermitian transpose. For uncorrelated and identical sensor

*noise, CUl1v2 5 n1v2I where the noise level n1v 2 is constant for*

*l5 1, c, L. *

The problem of central interest is to detect seismic events
and estimate the wave parameters using the array
observa-tions. The frequencies v1,c, v*J* chosen for processing

should cover most seismic energy. The unknown signal
*parameters Sl*1v

*j*2 and noise spectral parameters n1v*j*2 are

summarized in the vectors **S and ****n, respectively. **

**ML APPROACH**

The ML approach for parameter estimation for sensor array pro-cessing is well known to have excellent statistical performance and robustness. More importantly, the statistical properties of Fourier transformed data allow a natural combination of differ-ent frequency bins. These features make the ML method an attractive solution in processing broadband data. For geophysical applications, determination of arrival times of seis-mic waves plays a crucial role. In the ML framework, localiza-tion and deteclocaliza-tion of seismic activities can be carried out simultaneously. A generalized likelihood ratio test can be easily constructed if an estimate for wave parameters is available.

**PARAMETER ESTIMATION**

Following the properties of normal distribution and indepen-dency discussed in the previous section, the broadband log-like-lihood function is given by

*L1q, S, n2 5 2L a*
*L*
*l*51a
*J*
*j*51*cN*
logn1v*j*2 1
1
n1v*j*2
*1Xl*1v
*j*2
*2 H1vj*, *q2Sl*1v*j*22*H1Xl*1v*j2 2H1vj*, *q2Sl*1v*j*2d . (5)

Direct maximization of (5) is computationally prohibitive due to
a high dimension of parameter space. Fortunately, the signal and
*noise parameters are separable and the dependence on Sl*1v*j*2

and n1v*j*2 can be removed through replacing the unknown

sig-nal and noise parameters by their ML estimates at fixed and unknown wave parameters q [13]. Omitting the constant term, the broadband concentrated likelihood function is given by

*L*1q2 5 2 a
*J*

*j*51

log trS*1I 2 P1vj*, *q22 C^x*1v*j*2T, (6)

*where P*1v*j*, q2 is the projection matrix onto the column space

*of the transfer matrix H*1v*j*, q2 and the sample spectral matrix

*is given by C^x*1v*j*2 5

1

*L*g*Ll*51* Xl*1v*j2 Xl*1v*j*2*H*. The ML estimate q^

is obtained from maximizing (6)

q^ 5 arg max_{q} * L*1q2. (7)
It is worth noting that the noise level n1v*i*2 is estimated by

n^1v*j2 5 1/N tr31I 2 P1vj*, *q2C *^*x*1v*j*24. The likelihood function

(6) suggests that the optimizing parameter minimizes the geo-metric mean of estimated noise power over frequencies.

**SIGNAL DETECTION **

The detection of seismic events can be formulated as a multiple
*hypothesis test. Suppose the maximal number of signals to be M*max.

*The detection procedure discovers the mth signal by testing the *
*null hypothesis Hm against the alternative Am*.

*For m*_{5 1, }

*H*1 : Data contains only noise.

*X*1v*j2 5 U1vj*2

*A*1 : Data contains at least one signal.

*X*1v*j2 5 H*11v*j*; q1*2S*11v*j2 1 U1vj*2 (8)

*For m5 2, c, M*max

*Hm *: Data contains at most *1m 2 12 signals.*
*X*1v*j2 5 Hm*211v*j*; q*m*21*2Sm*211v*j2 1 U1vj*2
*Am : Data contains at least m signals.*

*X*1v*j2 5 Hm*1v*j*; q*m2Sm*1v*j2 1 U1vj*2 (9)

*Starting from the noise only hypothesis H*1, the test decides if a

signal is present. If no signal is present, the procedure stops. If a signal is detected, the procedure goes to the next step and decides if a second signal is present in the observation. This

*pro-cess continues until the maximum number of signals M*max is

reached. The subscript 1#2*m* in (8) and (9) emphasizes the matrix

dimension and the number of parameters associated with the assumed model.

*The test statistic for testing Hm against Am* is constructed

by the likelihood ratio principle [5]

*Tm*5 sup

q*m*

*L*1q*m*2 2 sup

q*m*21

*L*1q*m*212. (10)

Inserting the ML estimate q^*m* and q^*m*21 into the concentrated

likelihood function, the resulting test statistic is given by

*Tm*5 a
*J*
*j*51
*Tm*1v*j*2, (11)
*Tm*1v*j*2 5 loga1 1
*n*1
*n*2
* Fm*1v*j*; q^*m*2b, (12)
*Fm*1v*j*; q^*m*2 5
*n*2
*n*1
tr31P*m*1v*j*; q^*m2 2 Pm*211v*j*; q^*m*21*22C^x*1v*j*24
tr*31I 2 Pm*1v*j*; q^*m22C^x*1v*j*24
.
(13)
*The statistic Fm*1v*j*; q^*m2 is asymptotically Fn*1*, n*2-distributed
*under null hypothesis with degrees of freedom n*1*, n*2 [15].

Taking the estimated nonlinear parameters into account, the degrees of freedom are given by [16]

*n*1*5 L12 1 rm2, n*2*5 L12n 2 2m 2 rm*2 (14)

*with rm*5 dim1q*m*2 denoting the dimension of the nonlinear

*parameter vector associated with the mth signal. *

*At the mth stage, the test statistic is compared with a *
*thresh-old tm to decide whether to reject the hypothesis Hm* or retain it. The

*rejection of Hm* implies that a (further) signal is discovered;

other-wise, no signal is detected. More specifically,

*Tm$ tm1 reject Hm*, (15)
*Tm, tm1 retain Hm*. (16)
**CALCULATION OF TEST THRESHOLD**

*The threshold tm* is chosen to keep a prespecified significance level

(or false alarm rate) at a*m*. The probability that a signal is wrongly

detected should be kept less than a*m*. Typical values for a*m* are

1%, 5%, and 10%. Let F*m*1#*2 denote the null distribution of Tm*.

The threshold is determined by its inverse F*m*211#2 as follows:
*tm*5 F*m*211a*m*2. (17)

*In the narrow band case with J*5 1, the test (9) is equivalent to
*the test [17], the threshold can be easily derived from the *
F-distribution. However, in the broadband case, the null
distribution does not have a closed-form expression. To determine
*the threshold t*a*m*, methods such as normal approximation was
*proposed by [18]. Note that the statistics Tm*1v*j2, 1 j 5 1, c, J2 *

are independent, identically distributed with mean
m*m5 ETm*1v*j*2 and variance s*m*2 *5 Var1Tm*1v*j*22. For large

*sam-ples J, the test statistic Tm* is approximately normally distributed.

*tm*< m*m*1

s*m*

*"J*F

21_{1a}

*m*2, (18)

where F211a*m*2 is the inverse of the standard normal

distribu-tion funcdistribu-tion evaluated at a*m*. The formulas for computing m*m*

and s*m*2 can be found in [19]. A more accurate approximation

can be obtained from the Cornish-Fisher expansion that takes
*higher-order moments of Tm* into account [20] or

simulation-based bootstrap techniques [16], [21].

For simplicity, the test level a*m* is kept constant in each

detec-tion stage. The problem of controlling of the global test level can be achieved by controlling the familywise error-rate (FWE) through the sequentially rejective Bonferroni-Holm procedure [22] or the false discovery rate (FDR) through the Benjamini-Hochberg procedure [23]. The former is discussed in [18] and the latter is addressed in [24].

Given the Fourier transformed data, the ML approach for esti-mation and detection of seismic waves is summarized in Table 1.

**IMPLEMENTATION OF ML ESTIMATION**

The maximization of the likelihood function (6) involves
multi-dimensional search of nonlinear functions. The global maxima
can be computed by stochastic optimization methods such as
the genetic algorithm [25] or simulated annealing in a
straight-forward manner. A computationally more attractive
implemen-tation is to exploit the nature of sequential detection. Recall
*that the likelihood ratio (10) for testing Hm against Am* involves

the estimates q^*m*21 and q^*m*. Assuming that the components in

q^*m*21 for the first *1m 2 12 signals do not deviate from the *

cor-responding components in the estimate q^*m*, we may employ the

knowledge about the discovered signals and fix their values at
q^*m*21. Then the likelihood function is optimized globally over

the slowness parameter j*m associated with the mth signal *

j^*m*5 arg max_{j}

*m*

* L*1q^*m*21, j*m*2. (19)

The computational complexity for global optimization is signifi-cantly reduced. The so-obtained estimate can be refined by carry-ing out local optimization via local optimization methods such as Newton-type algorithms using the initial estimate

q^*m*3045 3q^*m*21, j^*m*4 (20)

with a smaller search interval. This implementation shows excel-lent performance in both simulation and experimental results in seismic and sonar applications [18], [16], [19].

**EXPLORATION OF GEOPHYSICAL STRUCTURES**

The geophysical structure of Earth’s surface has been of great interest to many areas including geology, environmental and civil engineering, and industrial applications. The geophysical inver-sion aims to construct a geophysical model that best fits experi-mental measurements. The majority of current inversion techniques assumes a constant model structure and computes the best model within this class of models [26], [27]. Contrary to

these methods, the ML approach considers models of various complexity and finds the best model in a hierarchy of models. The best model within each class maximizes the likelihood func-tion. The overall best model is selected by hypothesis testing.

Consider the following classes of models of increasing orders
M1( c( M*m*( cM*M*max. (21)
The major difference between the previously discussed seismic
location and detection is that the propagation waves are caused
by man-made signals in a controlled environment. The array
outputs collected during an experiment can be described by

*Xl*1v

*j2 5 d1vj*; q*m2Sl*1v*j2 1 Ul*1v*j*2, (22)

where q*m* contains the geophysical parameters associated with

the model class M*m*. Note that only one signal source is

*pres-ent; hence, Sl*1v

*j2 is a scalar and d1vj*; q*m*2 is a column vector.

*The transfer function d*1v*j*; q*m*2 is derived from underlying

geophysical laws. In a multilayered half-space model discussed
in [28], q*m* may contain thickness parameters of layers,

P-wave and S-wave propagation velocities. The model com-plexity increases with an increasing numbers of layers.

To select the best model, a sequential generalized
likeli-hood ratio test is suggested in [19] and [28]. Unlike the
detec-tion problem, the dimension of signal subspace remains the
same for all models. To obtain a test statistic with known null
distribution, the two models are tested against a joint signal
subspace spanned by both models. Starting from the smallest
model in (21), M1, we compare two adjacent models M*m* and

M*m*11 by carrying out a three-step test.

Step 1

*H1, m X5 U * No signal in the data.

*A1, m X5 3dm dm*114 c

*S*1

*S*2*d 1 U * M

*m* or M*m*11

gener-ates the data.

**[TABLE 1] ML APPROACH FOR ESTIMATION **
**AND DETECTION OF SEISMIC EVENTS.**

INPUT: FOURIER TRANSFORMED DATA 5X*l*1v

*j2, l 5 1, c, L, j 5 0, c, J 2 16 *

*MAXIMAL NUMBER OF SIGNALS M*max

*INITIAL VALUE FOR ESTIMATED NUMBER OF SIGNALS M*^ 5 0.
*FOR m5 1, c, M*max

1. WAVE PARAMETER ESTIMATION

FIND THE ML ESTIMATE: q^*m*5 arg maxq*mL*1q*m*2.
*2. COMPUTE THE TEST STATISTIC Tm* (10).

3. SIGNAL DETECTION

IF *Tm$ tm, REJECT Hm*, A SIGNAL IS DETECTED

* M^5 M^ 1 1, UPDATE THE NUMBER OF SIGNALS *
* m5 m 1 1, THE PROCEDURE CONTINUES *
ELSE BREAK THE LOOP

END END

OUTPUT: ESTIMATED WAVE PARAMETERS ^_{q ^}*M* AND DETECTED NUMBER OF

Step 2

*H2, m X5 dm1 U * M*m* generates the data.
*A2, m X5 3dm dm*114 c

*S*1

*S*2*d 1 U Some components cannot *

be modeled by M*m*.

Step 3

*H3, m X5 dm*11*1 U * M*m*11 generates the data.

*A3, m X5 3dm dm*114 c

*S*1

*S*2*d 1 U Some components cannot *

be modeled by M*m*.

For simplicity, the frequency and parameter dependence has
*been omitted in the formulation. The vector di5 d1vj*; q^*i*2,
*i5 m, m 1 1, is the transfer vector computed at the ML *

esti-mate obtained under M*i. In Step 1, if H1, m* is rejected, we

con-clude that a signal is present in the data and proceed to Step 2.
*If the hypothesis H2, m* is accepted, then model M*m* is the true

model and the test stops here. Otherwise, we proceed to Step 3.
Step 3 is a cross-check to test if the larger model M*m*11 is a

bet-ter one. The test statistics based on likelihood ratio have the same form as (10)–(13). The null distribution and test threshold can be derived in a similar manner as those for signal detection. The three-step procedure for model selection has been success-fully applied to near surface seismic model reconstruction. More details can be found in [28].

**EXPERIMENTAL RESULTS**

In this section, we demonstrate feasibility and performance of the
ML approach by processing seismic data recorded by the 25
verti-cal seismometers of the GERESS array located in the Bavarian
Forest, Germany. The output of each sensor is sampled with
40 Hz. Details about the GERESS arrays can be found in [29].
Because of the relatively small vertical aperture of GERESS
array, the elevation f*i* is taken to be zero.

A teleseismic event, which occurred on 1 March 1994 in the Eastern Mediterranean, was selected for analysis. We also applied the conventional sliding f–k analysis [30] to the same batch of data for comparison. Theoretical slow-ness values for each event are derived from the AK135 Earth model [31] as ref-erence. Relevant information about the selected events is collected in Table 2.

For the ML-based algorithm, we
employed sliding windows of length 3.2 s
with a shift of 0.5 s. In each window, data is
processed by the method described before
and the estimated parameters are plotted at
the center of the data window (e.g.,
Figure 1). The spectral density matrix is
*estimated with L*5 3 Thomson’s
*orthonor-mal windows. J*5 7 frequency bins
between 0.9 and 3.1 Hz are selected for
estimation and detection. The maximum
number of signals is assumed to be three.
The log–likelihood function is optimized
*over search space Vi* (apparent velocity)

[ 30 c504 km/s, a*i* (back–azimuth)

[ 30 c3604 degree. The corresponding
slowness vector is calculated according to
(3). These constraints avoid nonreasonable
*estimates for vi* and offer better accuracy

for a*i*. A genetic algorithm is used in global

optimization similar to those described in
[25]. The sequential testing is carried out
with test level a*m*5 0.033.

The conventional sliding f–k analysis also uses a window length of 3.2 s and a

**[TABLE 2] EVENT LIST FROM THE NATIONAL EARTHQUAKE **
**INFORMATION CENTER (NEIC).**

**ORIGIN TIME **
**HHMMSS.SS **
**LAT. **
**DEG. N **
**LONG. **
**DEG. E **
**DIST. **
**DEG. **
**THEOR. BAZ **
**DEG. **
**MAG. **
**MB LOCATION **
17:42:50 34.88 32.81 19.86 127.6 3.1 CYPRUS
0 30 60 90 120 150 180 210 240 270 300 330 360
−5,000
0
5,000

GERESS Data: 01.03.1994 17:43 34.9N 32.8E mb = 3.1 Cyprus

C2 (Counts) 0 30 60 90 120 150 180 210 240 270 300 330 360 0 2 4 Detected Sources 0 30 60 90 120 150 180 210 240 270 300 330 360 90 135 180 225 Baz (°) 0 30 60 90 120 150 180 210 240 270 300 330 360 0 5 10 15 20 Velocity (km/s) Time: 17:46:30 – 17:52:30 (s) (d) (c) (b) (a)

**[FIG1] Application of wideband ML method to a weak teleseismic event that occurred on **
**1 March 1994: (a) seismic data recorded at station C2, GERESS array in Germany, (b) **
**number of detected signals, (c) estimated values for back-azimuth (?) and theoretical **
**values for back-azimuth (22), and (d) estimated values for velocity (?) and theoretical **

shift of 0.5 s. The data in each window were first filtered with a Butterworth bandpass filter of 0.722.0 Hz, order three. Next, a wideband frequency-wavenumber spectrum analysis following [9] was applied and the corresponding results (back-azimuth and apparent velocity) were plotted at the center of the data window (e.g., Figure 2). Furthermore, the quality factor of the f–k analysis (ratio of the incoherent noise power to the total power of the sensor output) is displayed in addition to the beam steered in the estimated direction and centered at the main sta-tion C2 of the GERESS array.

It is shown in Figure 1 that earthquakes can be detected with good time resolution by its P phases. No signals are detected for the S phases. Sometimes there are false alarms. The estimates for back azimuth are quite accurate while the estimation of apparent velocity is slightly higher than theoretical values 1 3 , *2. More accurate estimates can be obtained by including frequency bins between 0.6 and 0.9 Hz. It was not done in our analysis because the presence of coherent noise structure in this frequency inter-val severely affects the detection performance. The signals detect-ed between 170 and 210 s after begin of analysis is another local seismic event that can be recognized by slightly different azimuths and different velocities at 180 s

and 190 s. The application of f–k analysis to this event shows that under this critical condition it is very difficult to claim that a seismic signal is present with help of the beamformer output and quality of the esti-mates (see Figure 2). In the absence of a seismic signal quality lies in the same range as in the presence of a seismic sig-nal. The estimates for back azimuth are similar to those given by ML approach and apparent velocities seem slightly better. There is no indication for the regional event detected by ML algorithm in f-k analysis.

In another analysis, an earthquake originated from Gulf of Aqaba in the Middle East is contaminated by a one magnitude-unit smaller preshock, locat-ed about 37 km from the main event. For the first time, both events can be detect-ed and localizdetect-ed accurately by applying the ML approach. The f-k analysis was only able to identify the strong event and could not discover the weak one. More details can be found in [6].

Experimental results show that the ML method provides not only reasonable esti-mates for wave parameters but also a reli-able indication about the presence of weak signals or multiple sources for low SNRs and long propagation paths. The detection ability is a significant advantage over routinely used f–k method. Thus the

ML method is a promising alternative to the conventional method in seismic application.

**CONCLUSIONS AND FUTURE WORK**

We discussed the ML approach for seismic parameter estima-tion and signal detecestima-tion, and its extension to geophysical model reconstruction. The ML approach has the advantage of excellent performance, high-resolution capability, and robust-ness against low SNRs and small samples. Exploiting the asymptotic normality of Fourier transformed data, the ML method allows an optimal combination of various frequency components. More importantly, localization and detection of seismic events can be carried out simultaneously. The multi-ple hypothesis test provides a statistically justified framework for signal detection and model fitting.

Experimental results show that the ML approach provides both accurate estimates for velocity and location and reliable indication about the presence of seismic events in critical sce-narios and long propagation paths. The detection ability will significantly enhance the power of modern seismic monitoring systems in minimizing the impact of natural disasters.

−150 0 150

GERESS Data: 01.03.1994 17:43 34.9N 32.8E mb = 3.1 Cyprus

BEAM (Counts) 0.2 0.4 0.6 0.8 1 Quality 0 30 60 90 120 150 180 210 240 270 300 330 360 45 90 135 180 225 Baz (°) 0 30 60 90 120 150 180 210 240 270 300 330 360 0 30 60 90 120 150 180 210 240 270 300 330 360 0 30 60 90 120 150 180 210 240 270 300 330 360 0 5 10 15 20 Velocity (km/s) Time: 17:46:30 − 17:52:30 (s) (d) (c) (b) (a)

**[FIG2] Application of sliding f-k analysis to a weak teleseismic event that occurred on **
**1 March 1994: (a) beamformer output, (b) quality of estimates, (c) estimated values for **
**back-azimuth (?) and theoretical values for back-azimuth (22), and (d) estimated **
**values for velocity (?) and theoretical values for velocity (3). **

In this article, we have considered ML methods under the plane wave model. Since the ML method is a model-based approach, it can be expected that further improvement in estimation accuracy can be achieved if a more realistic seismic model is incorporated into the algorithm. Another important issue is computational complexity. Although the processing speed of modern computers has grown rapidly, computational efficiency still plays an important role in practical systems. The recursive algorithms suggested in [32] and [33] facilitate computationally efficient implementation of ML methods. Their online processing capability can be very benefi-cial to seismic processing systems. As mentioned in the beginning of the article, the ML approach discussed here addresses seismic data processing from the view of a regional seismic array. How to integrate ML techniques into global seismic monitoring networks would be a challenging and very important topic in the future.

**AUTHORS**

* Pei-Jung Chung (peijung.chung@gmail.com) received the *
Dr.-Ing. with distinction in 2002 from Ruhr-Universität
Bochum, Germany. From 2002 to 2004, she held a postdoctoral
position at Carnegie Mellon University and the University of
Michigan, Ann Arbor, respectively. From 2004 to 2006, she was
an assistant professor with National Chiao Tung University, Hsin
Chu, Taiwan. In 2006, she joined the Institute for Digital
Communications, School of Engineering, the University of
Edinburgh, United Kingdom as a lecturer. Currently, she is an
associate member of the IEEE Signal Processing Society Sensor
Array Multichannel Technical Committee and serves on the
IEEE Communications Society, Multimedia Communications
Technical Committee as vice chair of the Interest Group on
Acoustic and Speech Processing for Communications. Her
research interests include array processing, statistical signal
pro-cessing, wireless multiple-input and multiple-output
communi-cations, and distributed processing in wireless sensor networks.

* Johann F. Böhme (Johann.Boehme@rub.de) received the *
Diplom in mathematics in 1966, the Dr.-Ing. in 1970, and the
Habilitation in 1977, both in computer science, from the Technical
University of Hannover, Germany, the University of Erlangen,
Germany, and the University of Bonn, Germany, respectively. From
1967 to 1974, he was with the Sonar-Research Laboratory of
Krupp Atlas Elektronik in Bremen, Germany. He was with the
University of Bonn until 1978 and the FGAN in Wachtberg,
Werthhoven. He has been a professor of signal theory in the
Department of Electrical Engineering and Information Sciences,
Ruhr-Universität Bochum, Germany, since 1980. His research
interests are in the areas of statistical signal processing and its
applications. He is an IEEE Life Fellow and recipient of the 2003
IEEE Signal Processing Society Technical Achievement Award.

**REFERENCES**

*[1] E. S. Husebye and A. Dainty, Eds., Monitoring a Comprehensive Test Ban *

*Treaty. Dordrecht, The Netherlands: Kluwer, 1996.*

[2] T. C. Bache, S. R. Bratt, H. J. Swanger, G. W. Beall, and F. K. Dashiell,
*“Knowledge-based interpretation of seismic data in the intelligent monitoring system,” Bull. Seism. *

*Soc. Amer., vol. 83, no. 5, pp. 1507–1526, 1993.*

[3] L. M. Haikin, A. F. Kushnir, and A. M. Dainty, “Combined automated and off-line computer processing system for seismic monitoring with small aperture arrays,”

*Seism. Res. Lett., vol. 69, no. 3, pp. 235–247, 1998.*

[4] G. Ekstrom, “Global detection and location of seismic sources by using surface
*waves,” Bull. Seism. Soc. Am. A, vol. 96, no. 4, pp. 1201–1212, 2006.*

*[5] J. F. Böhme, “Retrieving signals from array data,” in Monitoring a *

*Comprehensive Test Ban Treaty, E. S. Husebye and A. M. Dainty, Eds. Alvor, *

Portugal: NATO ASI, Jan. 23–Feb. 2, 1995, pp. 587–610.

[6] P.-J. Chung, M. L. Jost, and J. F. Böhme, “Seismic wave parameter estimation and
*signal detection using broadband maximum likelihood methods,” Comput. Geosci., *
vol. 27, no. 2, pp. 147–156, Mar. 2001.

[7] P.-J. Chung, A. B. Gershman, and J. F. Böhme, “Comparative study of
two-dimen-sional maximum likelihood and interpolated root-MUSIC with application to
*teleseis-mic source localization,” in Proc. IEEE Signal Processing Workshop Statistical *

*Sig-nal and Array Processing, Pocono Manor, PA, Aug. 14–16, 2000, pp. 68–72.*

[8] D. V. Sidorovich and A. B. Gershman, “Two-dimensional wideband interpolated
*root-MUSIC applied to measured seismic data,” IEEE Trans. Signal Processing, vol. *
46, no. 8, pp. 2263–2267, Aug. 1998.

*[9] J. Capon, “High-resolution frequency–wavenumber spectrum analysis,” Proc. *

*IEEE, vol. 57, no. 8, pp. 1408–1418, 1969.*

*[10] D. J. Thomson, “Spectrum estimation and harmonic analysis,” Proc. IEEE, vol. *
70, no. 9, pp. 1055–1096, Sept. 1982.

*[11] D. B. Percival and A. T. Walden, Spectral Analysis for Physical Applications, *

*Multitaper and Conventional Univariate Techniques. Cambridge, U.K.: Cambridge *

Univ. Press, 1993.

[12] D. Maiwald, D. V. Sidorovitch, and J. F. Bohme, “Broadband maximum likelihood
*wave parameter estimation using polarization sensitive arrays,” in Proc. IEEE Int. *

*Conf. Acoustics, Speech, and Signal Processing, vol. 4, 1993, pp. 356–359.*

*[13] J. F. Böhme, “Array processing,” in Advances in Spectrum Analysis and Array *

*Processing, S. Haykin, Ed. Englewood Cliffs, NJ: Prentice-Hall, 1991, pp. 1–63.*

*[14] D. R. Brillinger, Time Series: Data Analysis and Theory, expanded ed. San *
Francisco, CA: Holden-Day, 1981.

*[15] N. L. Johnson and S. Kotz, Continuous Univariate Distributions-2. New York: *
Wiley, 1970.

[16] D. Maiwald, “Breitbandverfahren zur Signalentdeckung und -ortung mit Sen-sorgruppen in Seismik- und Sonaranwendungen,” Dissertation, Faculty Elect. Eng., Ruhr Univ. Bochum, Germany, May 1995.

[17] R. H. Shumway, “Replicated time-series regression: An approach to signal
*esti-mation and detection,” in Handbook of Statistics, vol. 3, D. R. Brillinger and P. R. *
Krishnaiah, Eds. New York: Elsevier Science Publishers B.V., 1983, pp. 383–408.
[18] D. Maiwald and J. F. Böhme, “Multiple testing for seismic data using bootstrap,”
*in Proc. IEEE Int. Conf. Acoustics, Speech, and Signal Processing, Adelaide, South *
Australia, 1994, vol. 6, pp. 89–92.

[19] C. F. Mecklenbräuker, P. Gerstoft, J. F. Böhme, and P.-J. Chung, “Hypothesis
*testing for geoacoustic environmental models using likelihood ratio,” J. Acoust. Soc. *

*Amer., vol. 105, no. 3, pp. 1738–1748, Mar. 1999.*

[20] P.-J. Chung, M. Viberg, and C. F. Mecklenbräuker, “Broadband ML estimation
*un-der model orun-der uncertainty,” Signal Process., vol. 90, no. 5, pp. 1350–1356, May 2010.*
*[21] A. M. Zoubir and D. R. Iskander, Bootstrap Techniques for Signal Processing. *
Cambridge, U.K.: Cambridge Univ. Press, 2004.

*[22] S. Holm, “A simple sequentially rejective multiple test procedure,” Scand. J. *

*Stat-ist., vol. 6, pp. 65–70, 1979. *

[23] Y. Benjamini and Y. Hochberg, “Controlling the false discovery rate: A practical
*and powerful approach to multiple testing,” J. Roy. Statist. Soc. B, vol. 57, no. 1, pp. *
289–300, 1995.

[24] P.-J. Chung, J. F. Böhme, C. F. Mecklenbräuker, and A. O. Hero, “Detection of the
*number of signals using the Benjamini–Hochberg procedure,” IEEE Trans. Signal *

*Processing, vol. 55, no. 6, pp. 2497–2508, June 2007.*

*[25] D. E. Goldberg, Genetic Algorithms in Search, Optimization and Machine *

*Learning. Reading, MA: Addison-Wesley, 1989.*

[26] M. Sambridge, “Geophysical inversion with a neighbourhood algorithm. I.
*Searching a parameter space,” Geophys. J. Int., vol. 138, pp. 479–494, 1999. *
[27] N. Qiu, Q.-S. Liu, Q.-Y. Gao, and Q.-L. Zeng, “Combining genetic algorithm and
generalized least squares for geophysical potential field data optimized inversion,”

*IEEE Geosci. Remote Sensing Lett., vol. 7, no. 4, pp. 660–664, Oct. 2010.*

[28] M. Westebbe, J. F. Böhme, and H. Krummel, “Model fitting and testing in near
*surface seismics using maximum likelihood in frequency domain,” in Proc. 32nd *

*Asilomar Conf. Signals, Systems and Computers, Pacific Grove, CA, Nov. 1998, pp. *

1311–1315.

*[29] H.-P. Harjes, “Design and sitting of a new regional array in Central Europe,” Bull. *

*Seism. Soc. Amer. B, vol. 80, no. 6, pp. 1801–1817, 1990.*

[30] S. Mykkeltveit and H. Bungum, “Processing of regional seismic events using data
*from small-aperture arrays,” Bull. Seism. Soc. Am., vol. 74, no. 6, pp. 2313–2333, *
1984.

[31] B. L. N. Kennett, E. R. Engdahl, and R. Buland, “Constraints on seismic velocities
*in the earth from traveltimes,” Geophys. J. Int., vol. 122, pp. 108–124, 1995. *
[32] P.-J. Chung and J. F. Böhme, “EM and SAGE algorithms for towed array data,” in

*The Applications of Space–Time Adaptive Processing, R. Klemm, Ed. London: IEE *

Publishers, 2004, pp. 733–753.

[33] P.-J. Chung and J. F. Böhme, “Recursive EM and SAGE-inspired algorithms with
*application to DOA estimation,” IEEE Trans. Signal Processing, vol. 53, no. 8, pp. *