- Email: [email protected]
Automatic Travel Time Determination from a Frequency-domain Transfer Function: The Sompi Event Analysis
C H A P T E R 20
A UTOMATIC T RAVEL T IME D ETERMINATION FROM A F REQUENCY- DOMAIN T RANSFER F UNCTION : T HE S OMPI E VENT A NALYSIS Yoko Hasada a,b , Mineo Kumazawa b,c , Kayoko Tsuruga b,d , Takahiro Kunitomo b,c , and Junzo Kasahara b,c
Contents 1. Introduction 2. Theory 2.1. Concept of the model 2.2. Model of pulse sequence 2.3. Modeling dispersive waves 2.4. Modeling observed data 2.5. Determination of the complex travel times 2.6. Determination of the complex amplitudes 2.7. Maximum likelihood approach using linearized model 2.8. Determination of the model order 2.9. Error estimation by bootstrap method 3. Numerical Experiment 4. Discussion and Conclusion References
381 382 382 382 385 386 387 388 388 389 391 392 393 395
1. I NTRODUCTION ACROSS (Accurately-Controlled Routinely-Operated Signal System) is a subsurface exploring/monitoring system whose concept of data acquisition is distinctive: data obtained by the system constitute a complex a b c d
Graduate School of Environmental Studies, Nagoya University, Japan Japan Atomic Energy Agency, Japan Faculty of Science, Shizuoka University, Japan Geophysical Services Department, JGI Inc., Japan
E-mail address: [email protected] (Y. Hasada). Handbook of Geophysical Exploration: Seismic Exploration, Volume 40 ISSN 0950-1401, DOI: 10.1016/S0950-1401(10)04025-5
c 2010 Elsevier Ltd.
All rights reserved.
381
382
Y. Hasada et al.
discrete frequency series within a limited bandwidth (Kumazawa et al., 2000). To utilize the ACROSS system, we must establish a series of analysis methods, by which the subsurface structure is estimated from the observed frequency-domain data. As the first step, we developed a method to extract “events” localized in a time domain from frequency-domain data; this is called Sompi event analysis (Hasada et al., 2001). The data observed by ACROSS is a transfer function between the source input and the receiver output sampled at discrete frequencies in a limited bandwidth. It carries information along various paths between the transmitter and the receiver. In order to distinguish the information about each path from other paths, we attempt to translate a transfer function into a pulse sequence in the time domain, which is equivalent to an impulse response function if the frequency range of the transfer function is infinitely wide. The Sompi event analysis method has been developed to extract “events” localized in a time domain from a transfer function in the frequency domain, which is the basic analysis method of ACROSS (Hasada et al., 2001). Although the validity of Sompi event analysis has been confirmed through various numerical tests, especially for analysis of dispersive waves, there are some problems in its practical application. Here, we introduce a revised procedure, by means of weighted least-squares fitting, based on maximum likelihood estimation.
2. T HEORY 2.1. Concept of the model In Sompi event analysis, we assume that the complex frequency sequence to be analyzed is a transfer function or Green’s function between the source input and the receiver output sampled at discrete sequential frequencies. Through a kind of autoregressive modeling in the frequency domain, we obtain a set of “events” characterized by complex travel time and complex amplitude, in which the former expresses travel time and attenuation, whilst the latter expresses amplitude and phase angle.
2.2. Model of pulse sequence Using ACROSS, we measure the transfer function at discrete sequential frequencies in a limited band. A transfer function, which indicates how much the sinusoidal wave of each frequency is amplified and delayed in phase during transmission through a linear system, is obtained by dividing the complex frequency spectrum of the output signal by that of the input signal. The inverse Fourier transform of the transfer function yields an impulse response function.
383
The Sompi Event Analysis
Assuming a certain medium with neither dispersion nor attenuation, we consider that the waveform f (t) observed at a certain distance from the source is represented as a superposition of delayed pulses similar to the source time function f 0 (t), M X
f (t) =
Am f 0 (t − τm ),
(1)
m=1
where Am and τm are the amplification factor and the delay time of the m-th wave arrival, respectively. In the frequency domain, the Fourier transform of f (t) is described using the source spectrum F0 (ω) as F(ω) =
M X
Am e−iωτm F0 (ω).
(2)
m=1
The transfer function between the source spectrum and the observation spectrum is obtained as H (ω) =
M X F(ω) = Am e−iωτm , F0 (ω) m=1
(3)
which corresponds to the impulse response function h(t) =
M X
Am δ(t − τm ).
(4)
m=1
The transfer function in Eq. (3) contains the exponentials of imaginary numbers, whose real and imaginary parts are cosine and sine functions, respectively. Figure 1(a) and 1(b) show impulse response functions composed of one impulse and the corresponding transfer functions. The oscillation “frequency” in the complex transfer function is equal to the delay time in the impulse response function. To consider attenuation and phase shift, it is useful to introduce complex values for both amplification factor and delay time. We define the complex amplitude as αm = Am eiφm
(5)
and the complex travel time as qm = τm + iυm .
(6)
384
Y. Hasada et al.
(a) 10 8 6 4 2 0 0
1
2
3
4
5 6 time
7
8
9
10
1
2
3
4
5 6 time
7
8
9
10
8 6 4 2 0 1
2
3
4
5 6 time
7
8
9
10
8 6 4 2 0 1
2
3
4
5 6 time
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
1 0.5 0 –0.5 –1
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 frequency
(d) 10
0
1 0.5 0 –0.5 –1
frequency
(c) 10
0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 frequency
(b) 10 8 6 4 2 0 0
1 0.5 0 –0.5 –1
real part imaginary part
7
8
9
10
1 0.5 0 –0.5 –1
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 frequency
Figure 1 The impulse response functions (left) and the corresponding transfer functions in the frequency domain (right) for a nondispersive medium without attenuation (a and b) and with attenuation (c and d).
Then we have a model for transfer functions as H (ω) =
M X
αm e−iωqm ,
(7)
m=1
with the corresponding impulse response function being M 1 X iαm h(t) = R , π m=1 t − qm
(8)
where R indicates the real part. Figure 1(c) and (d) show examples for Eq. (7) and (8). The transfer function defined by Eq. (7) is a superposition of decaying or growing sinusoids. Figure 2 shows the schematic image of our model containing plural decaying sinusoids. We refer to each decaying or growing sinusoid as a “wave element” representing a pulse in the time
385
The Sompi Event Analysis Transfer function
Impulse response function
(a)
(b)
FT IFT
τ1 τ2
frequency
τ3 time
(c)
τ1 , υ1, A1 , φ1
(d)
τ2 , υ2, A2 , φ2 FT τ3 , υ3, A3 , φ3
IFT
frequency
time
Figure 2 A schematic image of a model in Sompi Event Analysis: (a) an example of a transfer function composed of the three decaying sinusoids shown in (c); (b) the corresponding impulse response function. The three waves in (c) correspond to the three pulses in (d). The impulse response functions in (b) and (d) are calculated by a complex Lorentzian model and are not Fourier transforms of (a) and (c).
domain. From Eq. (7), the growth rate is given by the imaginary part of the complex travel time, υm , which is required to be negative for the energy not to diverge to infinity. The constraint for real h(t) is H (ω) = H ∗ (−ω) for negative ω, where ∗ denotes the complex conjugate. Each pulse in Eq. (8) has the form of a linear combination of the Lorentz function and its Hilbert transform (π/2 radian phase shift); we refer to the content of the square bracket as a “complex Lorentzian” (Hasada et al., 2001). When αm is a real number, or φm in Eq. (5) is zero, the pulse becomes a Lorentz function whose half-value width is −2υm . This function has long tails on both sides of the peak, which seems to violate causality, since an attenuating medium necessarily involves dispersion, or frequency dependence. To satisfy causality and how to deal with dispersive transfer functions will be shown in the following section.
2.3. Modeling dispersive waves Here, for convenience, we consider a homogeneous dispersive media. The transfer function is written using a frequency dependent wave number k(ω)
386
Y. Hasada et al.
and propagating distance x, H (ω) = αe−ik(ω)x .
(9)
Expanding k(ω) into a Taylor series around ω0 up to the first order as ∂k ω0 ω − ω0 k(ω) = k(ω0 ) + (ω − ω0 ) = + , ∂ω ω=ω0 c(ω0 ) u(ω0 )
(10)
we obtain the transfer function −iω0 c(ωx
H (ω) = αe
0)
e
−i(ω−ω0 ) u(ωx
0)
= αe ˜ −i(ω−ω0 )q˜ ,
(11)
where c(ω0 ) and u(ω0 ) denote the phase velocity and group velocity at the frequency ω0 , respectively. Equation (11) is similar to Eq. (7), with the complex amplitude and the complex travel time defined as α˜ = αe
−iω0 c(ωx
0)
,
q˜ =
x . u(ω0 )
(12)
This indicates that the complex travel time obtained by our method is not related to the phase velocity, but to the group velocity. Thus, our model, as defined by Eq. (7), is approximately adaptable for a dispersive wave if the frequency range of the data to be analyzed is narrow enough. It is effective in dividing the data into segments, with bandwidth suitable for the frequency dependency inherent in the data when we turn to analyzing highly dispersive waves.
2.4. Modeling observed data We assume that the observed transfer function X j consists of the signal sequence H j and the noise sequence E j . X j = Hj + E j
( j = 1, . . . , N ),
(13)
where the real and the imaginary parts of E j have an independent and identical normal distribution N (0, σ j2 ), whose variance depends on the frequency. In the case of ACROSS observation, the errors of the measured transfer function can be evaluated by taking root mean squares of the data at the vicinal frequencies of the transmitted signals. The estimated error of the transfer function ε j is assumed to be proportional to the noise standard deviation σ j .
387
The Sompi Event Analysis
The signal sequence H j is the superposition of decaying sinusoids, as stated above Hj =
M X
αm e−2πi( j−1)1 f qm .
(14)
m=1
The parameters to be determined are the complex travel times qm and the complex amplitudes αm . Because the order of model M is also unknown, we change M in an arbitrary range and attempt to find the optimum model.
2.5. Determination of the complex travel times Hasada et al. (2001) proposed a procedure to determine complex travel times and complex amplitudes from observed transfer functions by means of autoregressive modeling. The signal sequence composed of M decaying sinusoids should satisfy the following homogeneous autoregressive (AR) equation of order M. M X
al H j−l = 0
( j = M + 1, . . . , N ),
(15)
l=0
where a = {a0 , . . . , a M } is a complex AR coefficient vector. We determine a by minimizing the fitting error 1 N−M
2 N M X X al X j−l → min, j=M+1 l=0
(16)
under a constraint to obtain a non-zero a, a · a∗ = 1.
(17)
Minimization of Eq. (16) with the condition Eq. (17), using a Lagrangian undetermined multiplier λ, leads to an eigenvalue problem, Pa = λa.
(18)
Here P is a non-Toeplitz Hermitian autocovariance matrix Pml =
1 N−M
N X j=M+1
∗ X j−l X ∗j−m = Plm
(l, m = 0, . . . , M). (19)
388
Y. Hasada et al.
Solving Eq. (18) yields M eigenvectors. We adopt the eigenvector a corresponding to the minimum eigenvalue, since λ reflects the noise power. In order to convert the AR coefficient vector a to a set of the complex travel times, we should solve the characteristic equation M X
al z −l = 0,
(20)
l=0
where z is the unit-time-shift operator (zH j = H j+1 ). The M characteristic solutions z 1 , . . . , z M of Eq. (20) are related to the complex travel times by z m = e−2πi1 f qm ,
(21)
which results in M complex travel times, qm =
i ln z m 2π1 f
(m = 1, . . . , M).
(22)
2.6. Determination of the complex amplitudes The complex amplitudes αm are determined by least-squares fitting. If we have an estimated error of X j , weighted least-squares fitting provides a more accurate estimation 2 M X w j X j − αm e−2πi( j−1)1 f qm → min . j=1 m=1
N X
(23)
The weight w j is defined in terms of a maximum likelihood procedure, as follows, 1 wj = 2 εj
X N k=1
1 , εk2
(24)
where ε j is the estimated error of X j .
2.7. Maximum likelihood approach using linearized model Hasada et al. (2004a,b) has proposed a procedure to improve estimation by using a maximum likelihood approach through linearization of the model.
389
The Sompi Event Analysis
Assuming that the estimated complex travel times and complex amplitudes are the initial values of the model parameters, H j = H j0 +
M X
0 0 δαm − 2πiαm δqm e−2πi( j−1)1 f qm ,
(25)
m=1 0 where H j0 is the transfer function calculated using the initial values, αm and qm0 are the initial values, and δαm and δqm are the small variations in parameters, respectively. These small variations are determined by linearized weighted least-squares fitting, 2 N M X X 0 0 0 −2πi( j−1)1 f qm w j X j − H j − δαm − 2πiαm δqm e → min . j=1 m=1
(26) Therefore, we have the revised model parameters qm = qm0 − δq m and 0 + δα . αm = αm m However, the revised parameters do not always yield a larger likelihood with the original or unlinearized model, particularly when the residual of the original model is large. Therefore we suggest that the revised parameters be adopted only when the AIC (Akaike’s Information Criterion) of the new model is less than the original one. Adding this procedure to every stage of analysis where complex travel times and complex amplitudes have been determined enables a more accurate estimation.
2.8. Determination of the model order In the original Sompi method for spectral analysis, the model order, or the number of wave elements, is determined by means of a two-parameter AIC proposed by Matsu’ura et al. (1990). One of the two parameters to specify the model is the AR order, the other is the number of wave elements. We describe each model as MODEL(M, M0 ), where M is the AR order and M0 is the number of wave elements, in which the largest M0 wave elements are adopted from the M candidates in terms of the mean power densities (MPD), defined as N 2 1 X −2πi( j−1)1 f qm MPDm = (27) αm e . N j=1 The complex amplitudes of MODEL (M, M0 ) are recalculated by leastsquares fitting for the M0 wave elements.
390
Y. Hasada et al.
The definition of the AIC is ˆ + 2k, AI C = −2l(θ)
(28)
ˆ is the maximum logarithmic likelihood and k is the number of where l(θ) free parameters. The two-parameter AIC of MODEL (M, M0 ) is given by AI C(M, M0 ) = 2N ln 2π σˆ 02 −
N X
ln w j + 2N + 8M0 ,
(29)
j=1
where σˆ 02 is the maximum likelihood estimate of the noise variance per unit weight calculated by 2 M0 N X X 1 w j X j − αm e−2πi( j−1)1 f qm . σˆ 02 = 2N j=1 m=1
(30)
We can determine the optimum MODEL (M, M0 ) by searching the minimum AIC for a certain range of both M and M0 . Using the AIC sometimes leads to a model with an improperly high order, in terms of practical application to analysis of observed data, because the parametric model is not fully suitable to a natural system. To avoid false solutions, we suggest screening of the solutions or wave elements by the growth/decay rates and the peak heights. The growth/decay rate of a wave element represents the pulse width in the time domain, it is expressed by the imaginary part of the complex travel time. In the analysis presented here, we adopt the wave elements whose imaginary parts are between −1/(N 1 f ) and 1/(N 1 f ). The peak height can be calculated from the growth/decay rates and the complex amplitude as |αm | 1 − e2π N 1 f υm , N 1 − e2π 1 f υm
(31)
where υm denotes the imaginary part of the complex travel time. We suggest that the threshold of peak height be set at the noise level in the time domain, which can be calculated from errors in the frequency domain or estimated by the time-domain waveform. For both screening procedures reasonable thresholds have to be adjusted according to the purpose. Figure 3 shows the flow of analysis, including the process by which to search for the optimum model order.
391
The Sompi Event Analysis
Transfer function with error
Determination of complex travel times Screening by decay rate
Determination of complex amplitudes
Screening by peak height Calculation of MPD and sorting by MPD
Recalculation of complex amplitudes
Calculation of AIC All candidate models
Finding minimum AIC model Bootstrap error estimation Any inappropriate estimates?
No
Yes Finding next minimum AIC model Optimum model
Figure 3
Schematic flow of analysis.
2.9. Error estimation by bootstrap method A bootstrap method is generally used to evaluate estimation errors. In our analysis, the residuals have to be weighted before resampling if we deal with data whose uncertainties depend on frequency. To identify wave elements obtained from bootstrap replications, we adopt the closeness of the complex travel times. We calculate the distance of all combinations of wave elements from the original dataset and bootstrap replications in the complex plane. Then, in descending order of MPD for the original wave elements, we choose the closest wave element from the bootstrap replication. As with the ordinal bootstrap method, the estimated standard error is given by the standard deviation of all replicated estimates.
392
Y. Hasada et al.
8 6 4 2 0 –2 –4 –6 15
0.2 0.15 0.1 0.05 0 –0.05 –0.1 –0.15
0
real imaginary
15.5
16
16.5
2
4
6
17 17.5 18 frequency (Hz)
8
10 12 time (s)
18.5
19
19.5
20
14
16
18
20
standard deviation of noise
Figure 4 The synthetic transfer function (top) and time-domain waveform (bottom) obtained by Fourier transformation. 1.5
1
0.5
0 15
16
17
18
19
20
frequency (Hz)
Figure 5
The frequency dependency of noise amplitude.
3. N UMERICAL E XPERIMENT To demonstrate automatic Sompi event analysis, we analyze a synthetic transfer function consisting of three decaying sinusoids. We synthesize the transfer function by three complex decaying sinusoids and a frequencydependent random noise sequence. The sampled frequency range is 1519.95 Hz with an interval of 0.05 Hz. Figure 4 shows the synthetic data, composed of three wave elements at 0.1 s, 1.0 s, and 4.0 s. The frequency dependency of the noise sequence is shown in Figure 5. We scanned the AR order from 1 to 20, and found the optimum model, MODEL(20, 5). The AIC of all the models are shown in Figure 6. In
393
The Sompi Event Analysis
AIC 2
1000
4 800
AR order (M)
6 8
600
10 12
400 14 16
200
18 20
0 5 10 number of wave elements (M0)
Figure 6
15
Two-parameter AIC of all models. The minimum AIC model is MODEL(20, 5).
Figure 7, we plot all the resultant wave elements in (a) the time domain and (b) the frequency domain, with Table 1 showing the estimated values of the complex travel times and the complex amplitudes in comparison with the input values. For the three largest wave elements, the estimated values recover the input values within the range of error. The two small redundant wave elements are considered to be created by noise. This result indicates that automatic processing using the Sompi event analysis method provides good estimates of the input complex decaying sinusoids, although it outputs some excessive wave elements, caused by noise.
4. D ISCUSSION AND CONCLUSION We have developed an automatic travel time determination by Sompi event analysis of transfer functions in the frequency domain. In Sompi event analysis, the observed transfer function is modeled as an autoregressive process and decomposed into wave elements. Given that an autoregressive model is based on the solution of a wave equation, the physical significance of the model is clear. However, extended modeling is necessary in cases of multiple reflections or resonance. Model order determination using a two-parameter AIC has previously been difficult, because this method tends to estimate improperly high order models. To overcome this difficulty and realize automatic processing, we proposed two screening criteria; growth/decay rates and peak heights, which
394
Y. Hasada et al.
(a)
(b)
original
original
#1
#1
#2
#2
#3
#3
#4
#4
#5
#5
residual
residual 0
5
10 time (s)
15
15
16
17
18
19
frequency (Hz)
Figure 7 The decomposed wave elements and residual in the (a) time domain and (b) frequency domain, calculated with the parameters by Sompi event analysis. The time-domain waveforms are calculated by Fourier transformation. The largest three wave elements recover the input values well. The smallest two are not input wave elements, but are rather created by noise. The residual has a similar frequency dependency to input noise amplitude. Table 1 Comparison of the estimated complex travel times and complex amplitudes with the input values. The errors are estimated by the bootstrap method. For the largest three wave elements, the estimated values recover the input values within the range of errors. The two smallest solutions are attributed to noises
Estimate
Error
Input
q1
0.11
−0.02i
q2
1.00
−0.09i
0.03
1.00
−0.10i
q3
4.01
+0.00i
0.03
4.00
+0.004i
q4
14.2
+0.1i
0.3
q5
10.8
+0.1i
0.1
α1
2.3
−1.1i
0.5
2.0
−1.0i
α2
3.9
+0.7i
0.8
4.0
+1.0i
α3
1.3
+0.1i
0.6
1.3
+0.0i
α4
−0.01
−0.00i
0.07
α5
−0.03
+0.10i
0.1
0.02
0.10
−0.01i
The Sompi Event Analysis
395
are confirmed, through numerical testing of synthetic data, to be effective in reducing redundant solutions. Repeated least-squares fitting by means of the linearized model increases the estimation accuracy. For more direct model optimization, some nonlinear optimization techniques should be considered.
REFERENCES Hasada, Y., Kumagai, H., Kumazawa, M., 2001. Autoregressive modeling of transfer functions in frequency domain to determine complex travel times. Earth Planets Space 53, 3–11. Hasada, Y., Kumazawa, M., Tsuruga, K., Kunitomo, T., 2004a. Traveltime estimation from the data acquired in frequency domain by ACROSS: Sompi Event Analysis. In: Special issue “Active Monitoring”. Chikyu Monthly 61–68 (in Japanese). Hasada, Y., Kumazawa, M., Tsuruga, K., Kunitomo, T., 2004b. Travel time estimation from a transfer function in frequency domain: The revised Sompi event analysis. In: Fujii, N., Kasahara, J., Higashihara, H., and Ogawa, K. (Eds.), The Proceedings of 1st International Workshop on Active Monitoring in the Solid Earth Geophysics (IWAM04), Task Group for Active Monitoring, Mizunami, Japan. pp. 241–244. Kumazawa, M., Kunitomo, T., Yokoyama, Y., Nakajima, T., Tsuruga, K., 2000. ACROSS: Theoretical and technical developments and prospect to future applications. Technical report, Japan Nuclear Cycle Development Institute, 9. pp. 115–129 (in Japanese). Matsu’ura, T., Imanishi, Y., Imanari, M., Kumazawa, M., 1990. Application of a new method of high-resolution spectral analysis, Sompi, for free induction decay of nuclear magnetic resonance. Applied Spectroscopy 44, 618–626.
A UTOMATIC T RAVEL T IME D ETERMINATION FROM A F REQUENCY- DOMAIN T RANSFER F UNCTION : T HE S OMPI E VENT A NALYSIS Yoko Hasada a,b , Mineo Kumazawa b,c , Kayoko Tsuruga b,d , Takahiro Kunitomo b,c , and Junzo Kasahara b,c
Contents 1. Introduction 2. Theory 2.1. Concept of the model 2.2. Model of pulse sequence 2.3. Modeling dispersive waves 2.4. Modeling observed data 2.5. Determination of the complex travel times 2.6. Determination of the complex amplitudes 2.7. Maximum likelihood approach using linearized model 2.8. Determination of the model order 2.9. Error estimation by bootstrap method 3. Numerical Experiment 4. Discussion and Conclusion References
381 382 382 382 385 386 387 388 388 389 391 392 393 395
1. I NTRODUCTION ACROSS (Accurately-Controlled Routinely-Operated Signal System) is a subsurface exploring/monitoring system whose concept of data acquisition is distinctive: data obtained by the system constitute a complex a b c d
Graduate School of Environmental Studies, Nagoya University, Japan Japan Atomic Energy Agency, Japan Faculty of Science, Shizuoka University, Japan Geophysical Services Department, JGI Inc., Japan
E-mail address: [email protected] (Y. Hasada). Handbook of Geophysical Exploration: Seismic Exploration, Volume 40 ISSN 0950-1401, DOI: 10.1016/S0950-1401(10)04025-5
c 2010 Elsevier Ltd.
All rights reserved.
381
382
Y. Hasada et al.
discrete frequency series within a limited bandwidth (Kumazawa et al., 2000). To utilize the ACROSS system, we must establish a series of analysis methods, by which the subsurface structure is estimated from the observed frequency-domain data. As the first step, we developed a method to extract “events” localized in a time domain from frequency-domain data; this is called Sompi event analysis (Hasada et al., 2001). The data observed by ACROSS is a transfer function between the source input and the receiver output sampled at discrete frequencies in a limited bandwidth. It carries information along various paths between the transmitter and the receiver. In order to distinguish the information about each path from other paths, we attempt to translate a transfer function into a pulse sequence in the time domain, which is equivalent to an impulse response function if the frequency range of the transfer function is infinitely wide. The Sompi event analysis method has been developed to extract “events” localized in a time domain from a transfer function in the frequency domain, which is the basic analysis method of ACROSS (Hasada et al., 2001). Although the validity of Sompi event analysis has been confirmed through various numerical tests, especially for analysis of dispersive waves, there are some problems in its practical application. Here, we introduce a revised procedure, by means of weighted least-squares fitting, based on maximum likelihood estimation.
2. T HEORY 2.1. Concept of the model In Sompi event analysis, we assume that the complex frequency sequence to be analyzed is a transfer function or Green’s function between the source input and the receiver output sampled at discrete sequential frequencies. Through a kind of autoregressive modeling in the frequency domain, we obtain a set of “events” characterized by complex travel time and complex amplitude, in which the former expresses travel time and attenuation, whilst the latter expresses amplitude and phase angle.
2.2. Model of pulse sequence Using ACROSS, we measure the transfer function at discrete sequential frequencies in a limited band. A transfer function, which indicates how much the sinusoidal wave of each frequency is amplified and delayed in phase during transmission through a linear system, is obtained by dividing the complex frequency spectrum of the output signal by that of the input signal. The inverse Fourier transform of the transfer function yields an impulse response function.
383
The Sompi Event Analysis
Assuming a certain medium with neither dispersion nor attenuation, we consider that the waveform f (t) observed at a certain distance from the source is represented as a superposition of delayed pulses similar to the source time function f 0 (t), M X
f (t) =
Am f 0 (t − τm ),
(1)
m=1
where Am and τm are the amplification factor and the delay time of the m-th wave arrival, respectively. In the frequency domain, the Fourier transform of f (t) is described using the source spectrum F0 (ω) as F(ω) =
M X
Am e−iωτm F0 (ω).
(2)
m=1
The transfer function between the source spectrum and the observation spectrum is obtained as H (ω) =
M X F(ω) = Am e−iωτm , F0 (ω) m=1
(3)
which corresponds to the impulse response function h(t) =
M X
Am δ(t − τm ).
(4)
m=1
The transfer function in Eq. (3) contains the exponentials of imaginary numbers, whose real and imaginary parts are cosine and sine functions, respectively. Figure 1(a) and 1(b) show impulse response functions composed of one impulse and the corresponding transfer functions. The oscillation “frequency” in the complex transfer function is equal to the delay time in the impulse response function. To consider attenuation and phase shift, it is useful to introduce complex values for both amplification factor and delay time. We define the complex amplitude as αm = Am eiφm
(5)
and the complex travel time as qm = τm + iυm .
(6)
384
Y. Hasada et al.
(a) 10 8 6 4 2 0 0
1
2
3
4
5 6 time
7
8
9
10
1
2
3
4
5 6 time
7
8
9
10
8 6 4 2 0 1
2
3
4
5 6 time
7
8
9
10
8 6 4 2 0 1
2
3
4
5 6 time
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
1 0.5 0 –0.5 –1
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 frequency
(d) 10
0
1 0.5 0 –0.5 –1
frequency
(c) 10
0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 frequency
(b) 10 8 6 4 2 0 0
1 0.5 0 –0.5 –1
real part imaginary part
7
8
9
10
1 0.5 0 –0.5 –1
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 frequency
Figure 1 The impulse response functions (left) and the corresponding transfer functions in the frequency domain (right) for a nondispersive medium without attenuation (a and b) and with attenuation (c and d).
Then we have a model for transfer functions as H (ω) =
M X
αm e−iωqm ,
(7)
m=1
with the corresponding impulse response function being M 1 X iαm h(t) = R , π m=1 t − qm
(8)
where R indicates the real part. Figure 1(c) and (d) show examples for Eq. (7) and (8). The transfer function defined by Eq. (7) is a superposition of decaying or growing sinusoids. Figure 2 shows the schematic image of our model containing plural decaying sinusoids. We refer to each decaying or growing sinusoid as a “wave element” representing a pulse in the time
385
The Sompi Event Analysis Transfer function
Impulse response function
(a)
(b)
FT IFT
τ1 τ2
frequency
τ3 time
(c)
τ1 , υ1, A1 , φ1
(d)
τ2 , υ2, A2 , φ2 FT τ3 , υ3, A3 , φ3
IFT
frequency
time
Figure 2 A schematic image of a model in Sompi Event Analysis: (a) an example of a transfer function composed of the three decaying sinusoids shown in (c); (b) the corresponding impulse response function. The three waves in (c) correspond to the three pulses in (d). The impulse response functions in (b) and (d) are calculated by a complex Lorentzian model and are not Fourier transforms of (a) and (c).
domain. From Eq. (7), the growth rate is given by the imaginary part of the complex travel time, υm , which is required to be negative for the energy not to diverge to infinity. The constraint for real h(t) is H (ω) = H ∗ (−ω) for negative ω, where ∗ denotes the complex conjugate. Each pulse in Eq. (8) has the form of a linear combination of the Lorentz function and its Hilbert transform (π/2 radian phase shift); we refer to the content of the square bracket as a “complex Lorentzian” (Hasada et al., 2001). When αm is a real number, or φm in Eq. (5) is zero, the pulse becomes a Lorentz function whose half-value width is −2υm . This function has long tails on both sides of the peak, which seems to violate causality, since an attenuating medium necessarily involves dispersion, or frequency dependence. To satisfy causality and how to deal with dispersive transfer functions will be shown in the following section.
2.3. Modeling dispersive waves Here, for convenience, we consider a homogeneous dispersive media. The transfer function is written using a frequency dependent wave number k(ω)
386
Y. Hasada et al.
and propagating distance x, H (ω) = αe−ik(ω)x .
(9)
Expanding k(ω) into a Taylor series around ω0 up to the first order as ∂k ω0 ω − ω0 k(ω) = k(ω0 ) + (ω − ω0 ) = + , ∂ω ω=ω0 c(ω0 ) u(ω0 )
(10)
we obtain the transfer function −iω0 c(ωx
H (ω) = αe
0)
e
−i(ω−ω0 ) u(ωx
0)
= αe ˜ −i(ω−ω0 )q˜ ,
(11)
where c(ω0 ) and u(ω0 ) denote the phase velocity and group velocity at the frequency ω0 , respectively. Equation (11) is similar to Eq. (7), with the complex amplitude and the complex travel time defined as α˜ = αe
−iω0 c(ωx
0)
,
q˜ =
x . u(ω0 )
(12)
This indicates that the complex travel time obtained by our method is not related to the phase velocity, but to the group velocity. Thus, our model, as defined by Eq. (7), is approximately adaptable for a dispersive wave if the frequency range of the data to be analyzed is narrow enough. It is effective in dividing the data into segments, with bandwidth suitable for the frequency dependency inherent in the data when we turn to analyzing highly dispersive waves.
2.4. Modeling observed data We assume that the observed transfer function X j consists of the signal sequence H j and the noise sequence E j . X j = Hj + E j
( j = 1, . . . , N ),
(13)
where the real and the imaginary parts of E j have an independent and identical normal distribution N (0, σ j2 ), whose variance depends on the frequency. In the case of ACROSS observation, the errors of the measured transfer function can be evaluated by taking root mean squares of the data at the vicinal frequencies of the transmitted signals. The estimated error of the transfer function ε j is assumed to be proportional to the noise standard deviation σ j .
387
The Sompi Event Analysis
The signal sequence H j is the superposition of decaying sinusoids, as stated above Hj =
M X
αm e−2πi( j−1)1 f qm .
(14)
m=1
The parameters to be determined are the complex travel times qm and the complex amplitudes αm . Because the order of model M is also unknown, we change M in an arbitrary range and attempt to find the optimum model.
2.5. Determination of the complex travel times Hasada et al. (2001) proposed a procedure to determine complex travel times and complex amplitudes from observed transfer functions by means of autoregressive modeling. The signal sequence composed of M decaying sinusoids should satisfy the following homogeneous autoregressive (AR) equation of order M. M X
al H j−l = 0
( j = M + 1, . . . , N ),
(15)
l=0
where a = {a0 , . . . , a M } is a complex AR coefficient vector. We determine a by minimizing the fitting error 1 N−M
2 N M X X al X j−l → min, j=M+1 l=0
(16)
under a constraint to obtain a non-zero a, a · a∗ = 1.
(17)
Minimization of Eq. (16) with the condition Eq. (17), using a Lagrangian undetermined multiplier λ, leads to an eigenvalue problem, Pa = λa.
(18)
Here P is a non-Toeplitz Hermitian autocovariance matrix Pml =
1 N−M
N X j=M+1
∗ X j−l X ∗j−m = Plm
(l, m = 0, . . . , M). (19)
388
Y. Hasada et al.
Solving Eq. (18) yields M eigenvectors. We adopt the eigenvector a corresponding to the minimum eigenvalue, since λ reflects the noise power. In order to convert the AR coefficient vector a to a set of the complex travel times, we should solve the characteristic equation M X
al z −l = 0,
(20)
l=0
where z is the unit-time-shift operator (zH j = H j+1 ). The M characteristic solutions z 1 , . . . , z M of Eq. (20) are related to the complex travel times by z m = e−2πi1 f qm ,
(21)
which results in M complex travel times, qm =
i ln z m 2π1 f
(m = 1, . . . , M).
(22)
2.6. Determination of the complex amplitudes The complex amplitudes αm are determined by least-squares fitting. If we have an estimated error of X j , weighted least-squares fitting provides a more accurate estimation 2 M X w j X j − αm e−2πi( j−1)1 f qm → min . j=1 m=1
N X
(23)
The weight w j is defined in terms of a maximum likelihood procedure, as follows, 1 wj = 2 εj
X N k=1
1 , εk2
(24)
where ε j is the estimated error of X j .
2.7. Maximum likelihood approach using linearized model Hasada et al. (2004a,b) has proposed a procedure to improve estimation by using a maximum likelihood approach through linearization of the model.
389
The Sompi Event Analysis
Assuming that the estimated complex travel times and complex amplitudes are the initial values of the model parameters, H j = H j0 +
M X
0 0 δαm − 2πiαm δqm e−2πi( j−1)1 f qm ,
(25)
m=1 0 where H j0 is the transfer function calculated using the initial values, αm and qm0 are the initial values, and δαm and δqm are the small variations in parameters, respectively. These small variations are determined by linearized weighted least-squares fitting, 2 N M X X 0 0 0 −2πi( j−1)1 f qm w j X j − H j − δαm − 2πiαm δqm e → min . j=1 m=1
(26) Therefore, we have the revised model parameters qm = qm0 − δq m and 0 + δα . αm = αm m However, the revised parameters do not always yield a larger likelihood with the original or unlinearized model, particularly when the residual of the original model is large. Therefore we suggest that the revised parameters be adopted only when the AIC (Akaike’s Information Criterion) of the new model is less than the original one. Adding this procedure to every stage of analysis where complex travel times and complex amplitudes have been determined enables a more accurate estimation.
2.8. Determination of the model order In the original Sompi method for spectral analysis, the model order, or the number of wave elements, is determined by means of a two-parameter AIC proposed by Matsu’ura et al. (1990). One of the two parameters to specify the model is the AR order, the other is the number of wave elements. We describe each model as MODEL(M, M0 ), where M is the AR order and M0 is the number of wave elements, in which the largest M0 wave elements are adopted from the M candidates in terms of the mean power densities (MPD), defined as N 2 1 X −2πi( j−1)1 f qm MPDm = (27) αm e . N j=1 The complex amplitudes of MODEL (M, M0 ) are recalculated by leastsquares fitting for the M0 wave elements.
390
Y. Hasada et al.
The definition of the AIC is ˆ + 2k, AI C = −2l(θ)
(28)
ˆ is the maximum logarithmic likelihood and k is the number of where l(θ) free parameters. The two-parameter AIC of MODEL (M, M0 ) is given by AI C(M, M0 ) = 2N ln 2π σˆ 02 −
N X
ln w j + 2N + 8M0 ,
(29)
j=1
where σˆ 02 is the maximum likelihood estimate of the noise variance per unit weight calculated by 2 M0 N X X 1 w j X j − αm e−2πi( j−1)1 f qm . σˆ 02 = 2N j=1 m=1
(30)
We can determine the optimum MODEL (M, M0 ) by searching the minimum AIC for a certain range of both M and M0 . Using the AIC sometimes leads to a model with an improperly high order, in terms of practical application to analysis of observed data, because the parametric model is not fully suitable to a natural system. To avoid false solutions, we suggest screening of the solutions or wave elements by the growth/decay rates and the peak heights. The growth/decay rate of a wave element represents the pulse width in the time domain, it is expressed by the imaginary part of the complex travel time. In the analysis presented here, we adopt the wave elements whose imaginary parts are between −1/(N 1 f ) and 1/(N 1 f ). The peak height can be calculated from the growth/decay rates and the complex amplitude as |αm | 1 − e2π N 1 f υm , N 1 − e2π 1 f υm
(31)
where υm denotes the imaginary part of the complex travel time. We suggest that the threshold of peak height be set at the noise level in the time domain, which can be calculated from errors in the frequency domain or estimated by the time-domain waveform. For both screening procedures reasonable thresholds have to be adjusted according to the purpose. Figure 3 shows the flow of analysis, including the process by which to search for the optimum model order.
391
The Sompi Event Analysis
Transfer function with error
Determination of complex travel times Screening by decay rate
Determination of complex amplitudes
Screening by peak height Calculation of MPD and sorting by MPD
Recalculation of complex amplitudes
Calculation of AIC All candidate models
Finding minimum AIC model Bootstrap error estimation Any inappropriate estimates?
No
Yes Finding next minimum AIC model Optimum model
Figure 3
Schematic flow of analysis.
2.9. Error estimation by bootstrap method A bootstrap method is generally used to evaluate estimation errors. In our analysis, the residuals have to be weighted before resampling if we deal with data whose uncertainties depend on frequency. To identify wave elements obtained from bootstrap replications, we adopt the closeness of the complex travel times. We calculate the distance of all combinations of wave elements from the original dataset and bootstrap replications in the complex plane. Then, in descending order of MPD for the original wave elements, we choose the closest wave element from the bootstrap replication. As with the ordinal bootstrap method, the estimated standard error is given by the standard deviation of all replicated estimates.
392
Y. Hasada et al.
8 6 4 2 0 –2 –4 –6 15
0.2 0.15 0.1 0.05 0 –0.05 –0.1 –0.15
0
real imaginary
15.5
16
16.5
2
4
6
17 17.5 18 frequency (Hz)
8
10 12 time (s)
18.5
19
19.5
20
14
16
18
20
standard deviation of noise
Figure 4 The synthetic transfer function (top) and time-domain waveform (bottom) obtained by Fourier transformation. 1.5
1
0.5
0 15
16
17
18
19
20
frequency (Hz)
Figure 5
The frequency dependency of noise amplitude.
3. N UMERICAL E XPERIMENT To demonstrate automatic Sompi event analysis, we analyze a synthetic transfer function consisting of three decaying sinusoids. We synthesize the transfer function by three complex decaying sinusoids and a frequencydependent random noise sequence. The sampled frequency range is 1519.95 Hz with an interval of 0.05 Hz. Figure 4 shows the synthetic data, composed of three wave elements at 0.1 s, 1.0 s, and 4.0 s. The frequency dependency of the noise sequence is shown in Figure 5. We scanned the AR order from 1 to 20, and found the optimum model, MODEL(20, 5). The AIC of all the models are shown in Figure 6. In
393
The Sompi Event Analysis
AIC 2
1000
4 800
AR order (M)
6 8
600
10 12
400 14 16
200
18 20
0 5 10 number of wave elements (M0)
Figure 6
15
Two-parameter AIC of all models. The minimum AIC model is MODEL(20, 5).
Figure 7, we plot all the resultant wave elements in (a) the time domain and (b) the frequency domain, with Table 1 showing the estimated values of the complex travel times and the complex amplitudes in comparison with the input values. For the three largest wave elements, the estimated values recover the input values within the range of error. The two small redundant wave elements are considered to be created by noise. This result indicates that automatic processing using the Sompi event analysis method provides good estimates of the input complex decaying sinusoids, although it outputs some excessive wave elements, caused by noise.
4. D ISCUSSION AND CONCLUSION We have developed an automatic travel time determination by Sompi event analysis of transfer functions in the frequency domain. In Sompi event analysis, the observed transfer function is modeled as an autoregressive process and decomposed into wave elements. Given that an autoregressive model is based on the solution of a wave equation, the physical significance of the model is clear. However, extended modeling is necessary in cases of multiple reflections or resonance. Model order determination using a two-parameter AIC has previously been difficult, because this method tends to estimate improperly high order models. To overcome this difficulty and realize automatic processing, we proposed two screening criteria; growth/decay rates and peak heights, which
394
Y. Hasada et al.
(a)
(b)
original
original
#1
#1
#2
#2
#3
#3
#4
#4
#5
#5
residual
residual 0
5
10 time (s)
15
15
16
17
18
19
frequency (Hz)
Figure 7 The decomposed wave elements and residual in the (a) time domain and (b) frequency domain, calculated with the parameters by Sompi event analysis. The time-domain waveforms are calculated by Fourier transformation. The largest three wave elements recover the input values well. The smallest two are not input wave elements, but are rather created by noise. The residual has a similar frequency dependency to input noise amplitude. Table 1 Comparison of the estimated complex travel times and complex amplitudes with the input values. The errors are estimated by the bootstrap method. For the largest three wave elements, the estimated values recover the input values within the range of errors. The two smallest solutions are attributed to noises
Estimate
Error
Input
q1
0.11
−0.02i
q2
1.00
−0.09i
0.03
1.00
−0.10i
q3
4.01
+0.00i
0.03
4.00
+0.004i
q4
14.2
+0.1i
0.3
q5
10.8
+0.1i
0.1
α1
2.3
−1.1i
0.5
2.0
−1.0i
α2
3.9
+0.7i
0.8
4.0
+1.0i
α3
1.3
+0.1i
0.6
1.3
+0.0i
α4
−0.01
−0.00i
0.07
α5
−0.03
+0.10i
0.1
0.02
0.10
−0.01i
The Sompi Event Analysis
395
are confirmed, through numerical testing of synthetic data, to be effective in reducing redundant solutions. Repeated least-squares fitting by means of the linearized model increases the estimation accuracy. For more direct model optimization, some nonlinear optimization techniques should be considered.
REFERENCES Hasada, Y., Kumagai, H., Kumazawa, M., 2001. Autoregressive modeling of transfer functions in frequency domain to determine complex travel times. Earth Planets Space 53, 3–11. Hasada, Y., Kumazawa, M., Tsuruga, K., Kunitomo, T., 2004a. Traveltime estimation from the data acquired in frequency domain by ACROSS: Sompi Event Analysis. In: Special issue “Active Monitoring”. Chikyu Monthly 61–68 (in Japanese). Hasada, Y., Kumazawa, M., Tsuruga, K., Kunitomo, T., 2004b. Travel time estimation from a transfer function in frequency domain: The revised Sompi event analysis. In: Fujii, N., Kasahara, J., Higashihara, H., and Ogawa, K. (Eds.), The Proceedings of 1st International Workshop on Active Monitoring in the Solid Earth Geophysics (IWAM04), Task Group for Active Monitoring, Mizunami, Japan. pp. 241–244. Kumazawa, M., Kunitomo, T., Yokoyama, Y., Nakajima, T., Tsuruga, K., 2000. ACROSS: Theoretical and technical developments and prospect to future applications. Technical report, Japan Nuclear Cycle Development Institute, 9. pp. 115–129 (in Japanese). Matsu’ura, T., Imanishi, Y., Imanari, M., Kumazawa, M., 1990. Application of a new method of high-resolution spectral analysis, Sompi, for free induction decay of nuclear magnetic resonance. Applied Spectroscopy 44, 618–626.