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 criticalconditions 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
Xl1v2 5 1
"T a
T21
t50
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
di1v2 5 3e2jvji
Tp
1,c,e2jvjiTpN4T, (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 ji is related to the propagation
velocity Vi, azimuth ai, and elevation fi through the following
equation: ji5
1
Vi 3cos f
i sin ai, cos fi cos ai, sin fi4T. (3)
In the frequency domain, array outputs can be approximately described through a nonlinear regression model
Xl1v2 5 H1v,q2 Sl1v2 1 Ul1v2, (4)
where the transfer function H1v, q2 5 3d11v2 cdM1v24
con-sists of M steering vectors. The nonlinear wave parameters are summarized in q 5 3j1,c, jM4. Sl1v2 and Ul1v2 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 Xl1v2 is characterized by the following
sta-tistical features:
■ For large T, Xl1v2 is complex normally distributed. ■ When the signal Sl1v2 is considered as deterministic,
the distribution of Xl1v2 is completely specified by the
mean H1v,q2Sl1v2 and covariance matrix C
Ul1v2, denoted by Xl1v2 |Nc1H1v, q2Sl1v2, C
Ul1v22.
■ For different frequency bins, vi2 vj, Xl1vi2 and Xl1vj2
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 E3Ul1v2Ul1v2H4, 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, vJ chosen for processing
should cover most seismic energy. The unknown signal parameters Sl1v
j2 and noise spectral parameters n1vj2 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 l51a J j51cN logn1vj2 1 1 n1vj2 1Xl1v j2 2 H1vj, q2Sl1vj22H1Xl1vj2 2H1vj, q2Sl1vj2d . (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 Sl1vj2
and n1vj2 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
L1q2 5 2 a J
j51
log trS1I 2 P1vj, q22 C^x1vj2T, (6)
where P1vj, q2 is the projection matrix onto the column space
of the transfer matrix H1vj, q2 and the sample spectral matrix
is given by C^x1vj2 5
1
LgLl51 Xl1vj2 Xl1vj2H. The ML estimate q^
is obtained from maximizing (6)
q^ 5 arg maxq L1q2. (7) It is worth noting that the noise level n1vi2 is estimated by
n^1vj2 5 1/N tr31I 2 P1vj, q2C ^x1vj24. 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 Mmax.
The detection procedure discovers the mth signal by testing the null hypothesis Hm against the alternative Am.
For m5 1,
H1 : Data contains only noise.
X1vj2 5 U1vj2
A1 : Data contains at least one signal.
X1vj2 5 H11vj; q12S11vj2 1 U1vj2 (8)
For m5 2, c, Mmax
Hm : Data contains at most 1m 2 12 signals. X1vj2 5 Hm211vj; qm212Sm211vj2 1 U1vj2 Am : Data contains at least m signals.
X1vj2 5 Hm1vj; qm2Sm1vj2 1 U1vj2 (9)
Starting from the noise only hypothesis H1, 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 Mmax is
reached. The subscript 1#2m 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]
Tm5 sup
qm
L1qm2 2 sup
qm21
L1qm212. (10)
Inserting the ML estimate q^m and q^m21 into the concentrated
likelihood function, the resulting test statistic is given by
Tm5 a J j51 Tm1vj2, (11) Tm1vj2 5 loga1 1 n1 n2 Fm1vj; q^m2b, (12) Fm1vj; q^m2 5 n2 n1 tr31Pm1vj; q^m2 2 Pm211vj; q^m2122C^x1vj24 tr31I 2 Pm1vj; q^m22C^x1vj24 . (13) The statistic Fm1vj; q^m2 is asymptotically Fn1, n2-distributed under null hypothesis with degrees of freedom n1, n2 [15].
Taking the estimated nonlinear parameters into account, the degrees of freedom are given by [16]
n15 L12 1 rm2, n25 L12n 2 2m 2 rm2 (14)
with rm5 dim1qm2 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 am. The probability that a signal is wrongly
detected should be kept less than am. Typical values for am are
1%, 5%, and 10%. Let Fm1#2 denote the null distribution of Tm.
The threshold is determined by its inverse Fm211#2 as follows: tm5 Fm211am2. (17)
In the narrow band case with J5 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 tam, methods such as normal approximation was proposed by [18]. Note that the statistics Tm1vj2, 1 j 5 1, c, J2
are independent, identically distributed with mean mm5 ETm1vj2 and variance sm2 5 Var1Tm1vj22. For large
sam-ples J, the test statistic Tm is approximately normally distributed.
tm< mm1
sm
"JF
211a
m2, (18)
where F211am2 is the inverse of the standard normal
distribu-tion funcdistribu-tion evaluated at am. The formulas for computing mm
and sm2 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 am 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^m21 and q^m. Assuming that the components in
q^m21 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^m21. Then the likelihood function is optimized globally over
the slowness parameter jm associated with the mth signal
j^m5 arg maxj
m
L1q^m21, jm2. (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^m3045 3q^m21, j^m4 (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( Mm( cMMmax. (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
Xl1v
j2 5 d1vj; qm2Sl1vj2 1 Ul1vj2, (22)
where qm contains the geophysical parameters associated with
the model class Mm. Note that only one signal source is
pres-ent; hence, Sl1v
j2 is a scalar and d1vj; qm2 is a column vector.
The transfer function d1vj; qm2 is derived from underlying
geophysical laws. In a multilayered half-space model discussed in [28], qm 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 Mm and
Mm11 by carrying out a three-step test.
Step 1
H1, m X5 U No signal in the data.
A1, m X5 3dm dm114 c
S1
S2d 1 U M
m or Mm11
gener-ates the data.
[TABLE 1] ML APPROACH FOR ESTIMATION AND DETECTION OF SEISMIC EVENTS.
INPUT: FOURIER TRANSFORMED DATA 5Xl1v
j2, l 5 1, c, L, j 5 0, c, J 2 16
MAXIMAL NUMBER OF SIGNALS Mmax
INITIAL VALUE FOR ESTIMATED NUMBER OF SIGNALS M^ 5 0. FOR m5 1, c, Mmax
1. WAVE PARAMETER ESTIMATION
FIND THE ML ESTIMATE: q^m5 arg maxqmL1qm2. 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 Mm generates the data. A2, m X5 3dm dm114 c
S1
S2d 1 U Some components cannot
be modeled by Mm.
Step 3
H3, m X5 dm111 U Mm11 generates the data.
A3, m X5 3dm dm114 c
S1
S2d 1 U Some components cannot
be modeled by Mm.
For simplicity, the frequency and parameter dependence has been omitted in the formulation. The vector di5 d1vj; q^i2, i5 m, m 1 1, is the transfer vector computed at the ML
esti-mate obtained under Mi. 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 Mm 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 Mm11 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 fi 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 L5 3 Thomson’s orthonor-mal windows. J5 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, ai (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 ai. A genetic algorithm is used in global
optimization similar to those described in [25]. The sequential testing is carried out with test level am5 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.
