
Citation: | Chunyan Zhang, Sidao Ni (2011). Wave separation for the great Sumatra-Andaman earthquake with regional seismic array. Earthq Sci 24(1): 127-132. DOI: 10.1007/s11589-011-0776-4 |
We analyzed the seismic waveforms from the December 26, 2004 Sumatra-Andaman earthquake recorded at broadband seismic stations in western Europe. Previous studies involving of the beam-forming technique and high frequency analysis suggest that the earthquake ruptured with a duration of around 500 s. This very long duration makes P wave overlap with later arrivals such as PP wave, which follows P in about 200 s. Since P waves are crucial for modeling earthquake processes, we propose an iterative method to separate P and PP waveforms. The separated P waveform confirms a second large energy release around 300 s after the initial rupture. The iterative signal separation technique is particularly useful for mixed signals that are not independent and the number of recording stations far exceeds number of mixed signal sources.
The devastating great Sumatra-Andaman earthquake (December 26, 2004) is the strongest earthquake in the past four decades, and is well observed with various instrument arrays such as seismometers, hydroacoustic stations, infra-sonic stations, tidal gauges, and satellites (GPS for ground displacement, altimetry for sea surface variation) (Ammon et al., 2005; Banerjee et al., 2007; Catherine et al., 2005; Guilbert et al., 2005). Particularly, hundreds of modern broadband seismometers all over the world provided unprecedented records for the study of how great earthquakes occur and its effects on tsunami excitation. Various methods have been applied to study the source processes of this great earthquake, and some of the fundamental parameters have been obtained by array techniques. Ni et al. (2005) used global seismic network (sparse array) to study the high frequency behavior of P waves from the earthquake, and they found that the quake started from northwestern tip of Sumatra (solid circle in Figure 1) and ruptured towards Andaman Islands (around latitude of 13 degrees north) with a length of 1 300 km and lasted about 500 s, which is one of the longest ones in recorded history. Their rupture model can explain the extended distribution of aftershocks (open circles in Figure 1). Kruger and Ohrnberger (2005) obtained rupture length and duration of the earthquake by analyzing broadband seismic data recorded by German regional seismic array with beam-forming method. For the case of curved wave front, similar technique is applied to the short period data from the dense Japanese seismic array (Hi-net), again leading to compatible estimates of rupture length and duration (Ishii et al., 2005).
All these important results argue for importance of array technology in observing and modeling earthquake, thus substantially contributing to hazard mitigation. However, details of this huge earthquakes are far from being understood and one major obstacle is due to its extra long duration (500 s). P wave, which is always the first arrival, is routinely analyzed to study earthquake processes because P wave is very little affected by Earth structure, thus providing the most direct and reliable information of earthquakes (Kikuchi and Kanamori, 1982). All the major earthquakes (except the great Sumatra earthquake) in the past four decades have been successfully modeled with P waves. But for this great earthquake, P waveform is severely contaminated by PP wave which is about 200 s after P wave for epicentral distances around 80 degrees. For example, for the seismic data recorded by seismic stations in Germany, Switzerland and Austria (referred as European array later on)(Figure 2), it is very hard to tell how long the earthquake lasted, because seismic signal does not appear to lose strength even at 650 s after the earthquake occurred. This is because P waves after 200 s are overlapped by PP waves (indicated by the thick dashed line in Figure 2, according to theoretical arrival time by IASP91 Earth model). Therefore remaining 300 s of P wave is not directly available for modeling the details of the earthquake processes, and signal separation has to be performed to extract P wave from the contaminated seismic records.
Generally speaking, there are three goals for array signal processing, i.e., to determine (1) number of sources or model order, (2) direction of arrivals (DOA) and (3) extract of waveforms (Vorobyov et al., 2005). And various algorithms have been proposed to achieve one or two of the goals. For example, Ng et al. (2005) proposed a method to simultaneously determine number of sources and DOA, based on Markov chain Monte Carlo methods, and their methods work well for both broadband and narrow band signals. As for DOA, they are numerous methods based on MUSIC, CBF, constant modulus, maximum likelihood, information theory (Weinstein et al., 1993; Amindavar and Reza, 2005; Ozcetin 2003; Vorobyov et al., 2005). Impressively, an efficient method that estimates DOA and channel parameters at the same time is described in the work by Amindavar and Reza (2005). As for signal separation, the decorrelation algorithm is developed to separate signals mixed in two channels (Weinstein et al., 1993).
However, most of these methods can not be applied to studies on great earthquakes because of the following reasons: (1) seismic signals are wideband, and large earthquakes typically have more low frequency engery; (2) large earthquake sources are spatially and temporally distributed. For example, the great Sumatra earthquake ruptured almost a 1 300 km length; (3) mixed signals are usually not independent of each other. (4) multipathing effect is more or less well known. For example the travel time of PP can be calculated fairly accurate according to some reference Earth models, so we do not have to estimate the channel properties if we view the Earth as the transmission channel for seismic waves. It is the first three reasons (very wideband and spatially distributed, and correlation between mixed signals) that make many of the above mentioned methods not immediately applicable for seismic modeling.
For example, the method of multi-channel signal separation requires the two mixed sources are uncorrelated (Weinstein et al., 1993). And from this assumption, two mixed sources can be recovered from two outputs. Signal separation based independent component analysis (ICA) also requires that the mixing sources are independent. For the problem of separating P waves from PP waves for large earthquakes, we can not assume P and PP are uncorrelated, actually P and PP correlates fairly well since both are excited by the same earthquake source though they are radiated with different strength because of the radiation pattern. However, seismologists have many more observing stations than the number of mixing signals, and we can take advantage of this and develop an efficient, physics-based algorithm for separating mixing seismic signals, instead of just assuming independence and performing purely blind signal separation.
Low frequency P wave (< 0.5 Hz) propagates in the Earth with little modification and the transfer function is usually modeled as delta function (Kikuchi and Kanamori, 1982). PP wave is a little bit more variable, but still can be modeled as traveling wave whose arrival time can be affected by heterogeneity in the Earth with little medication on waveforms, as far as recording seismic array is not very large (Kruger and Ohrnberger, 2005). So we can assume, the seismic signals recorded at the k-th station can be modeled as, for k=1 to N:
|
(1) |
In equation (1), N is the number of seismic stations, and N=27 for this study (Figure 1). Qk(t) is noise which could be the real noise in the Earth or scattered waves. SP(t) and SPP(t) are the P and PP signals from the earthquake respectively, Tk is the differential time at the k-th station between PP wave and P wave. Because P and PP waves travel at different speeds (or slowness, P wave needs about 5 s to travel 1.0 degree distance, while PP wave needs about 8 s, at epicentral distance of 84 degrees). For example, 1D ray tracing for PREM predicts that PP and P are separate by 201.34 s at station SENIN, and by 188.81 s at station RUE. It is the difference in traveling slowness that makes separation of P and PP possible.
In frequency domain, equation (1) can be written as
|
(2) |
It is apparent that equation (2) is very similar to principle component analysis (PCA). The idea is that X1(ω), X2(ω), ···, XN(ω), these N≫1 functions of ω can be expressed as linear combination of two functions SP(ω) and SPP(ω). This also leads to the possibility of identifying how many signals are mixed in the N seismic signals, based on PCA.
However, the separation is very difficult to do directly in frequency domain, because the assumption in equation (3) in the work by Weinstein et al. (1993) usually does not meet because realistic Xk(t) consists of other signals such as PPP, and some scattered waves.
To derive a time domain algorithm, we can start from equation (2). If we sum Xk(ω) for k=1 to N, we have
|
(3) |
For large enough N, the standard deviation σ[∑Qk(ω)]~Sqrt(N) caused by the noise spectrum at the k-th station Qk(ω) are typically assumed to be independent at different stations. And it is expected that σ[∑exp(iωTk)] is about Sqrt(N) if Tk is randomly distributed. Now we can get a decent estimate of SP(ω) = ∑Xk(ω)/N with an accuracy of Sqrt(N).
Then we can estimate SPP(ω) from equation (2) again with an accuracy of Sqrt(N)
|
(4) |
From equation (3) we can estimate SP(ω) from the updated SPP(ω) with a better approximation, and then with SP(ω) we can estimate SPP(ω) with equation (4) better. Thus we have an iterative algorithm.
In time domain, the algorithms look like this:
Initially set SPP(t) = 0;
step 1:
for each k
shift SPP(t) back by Tk
Xk(t) = Xk(t) − SPP(t)
end
SP(t)=average of Xk(t)
step 2:
For each k
Xk(t) = Xk(t) − SP(t)
Shift Xk(t) back by −Tk
End
SPP(t)=average of Xk(t)
Go to step 1 until minimum mismatch is obtained.
We implemented the algorithm with Matlab, and applied the algorithm to 27 traces of seismograms, the differential arrival time between P and PP for each station is calculated with IASP91 model. The result is shown in Figure 3. SPP(t) (the second panel) is almost zero before 200 s, which should be the case because PP does not arrive until 200 s after P according to theoretical travel time calculations. Then SPP(t) between 200 and 300 s shows similar waveform as SP(t) between 0 and 100 s, confirming that SP(t) and SPP(t) correlate, which should be the case because both waves are radiated from the same source. SP(t) also shows stronger signal between 300 and 400 s, and loses strength after 450 s, consistent with results in Ni et al. (2005) and Ishii et al. (2005).
To verify that the separation of P and PP is successful, we compare original seismic record (dashed line in Figure 3c) with reconstructed signal SP(t)+SPP(t − Tk) for station UBBA. The match is much better than just a simple average of all the Xk(t) (Figure 3d).
The algorithm converges (Figure 4) fairly fast because for each time the estimate values approach the true values with an error of the order of 1/Sqrt(N). So this algorithm only works for the case where number of stations far exceeds the number of mixed signals.
Though the great Sumatra-Andaman earthquake is recorded with unprecedented number of modern broadband seismic stations, the earthquake source processes are still not well understood because P wave is not available for earthquake source modeling before it is separated from PP waves. Here we proposed an iterative method to separate P from PP wave with data from a regional seismic network in western Europe. This separated P wave confirms that there is a second large energy radiation around 300–400 s, and the earthquake lasted at least 450 s (Ni et al., 2005; Ishii et al., 2005; Kruger and Ohrnberger, 2005). We will apply the algorithm to other regional seismic work to get separated P wave in order to invert much better earthquake processes.
We also note that our algorithm does not consider the distributed nature of great earthquake sources. However, the differential time between PP and P does not change much (3 s) for the regional array we used here, as compared with the dominant period of the seismic data (20–30 s) (Kruger and Ohrnberger, 2005). Therefore, our algorithm should provide a fairly good estimate of P and PP waves. In future studies, we will take into account that SP(t) and SPP(t) are radiated from different parts of the earthquake rupture zone. This correction can be straightforwardly performed since our algorithm is in time domain.
In general, we proposed an iterative method to separate mix signals where the source durations overlap each other substantially and the sources correlate with each other. Our method works well if number of observation far exceeds number of mix signals, which are usually true for geophysical problems, for example, modern global seismic networks feature hundreds of broadband stations.
We are grateful that IRIS and ORFEUS provide broadband seismic data for the great Sumatra earthquake. The study is supported by CAS fund (KZCX2-YW-116-1), National Natural Science Foundation of China (40821160549 and 41074032) and China Earthquake Administration fund (200808078).
Amindavar H and Reza A M (2005). A new simultaneous estimation of directions of arrival and channel parameters in a multipath environment. IEEE Trans Signal Processing 53(2): 471-483. doi: 10.1109/TSP.2004.840673
|
Ammon C J, Ji C, Thio H K, Robinson D, Ni S D, Hjorleifsdottir V, Kanamori H, Lay T, Das S, Helmberger D, Ichinose G, Polet J and Wald D (2005). Rupture process of the 2004 Sumatra-Andaman earthquake. Science 308(5725): 1133-1139. doi: 10.1126/science.1112260
|
Banerjee P, Pollitz F, Nagarajan B and Bürgmann R (2007). Coseismic slip distributions of the 26 December 2004 Sumatra-Andaman and 28 March 2005 Nias earthquakes from GPS static offsets. Bull Seismol Soc Am 97(1A): S86-S102. doi: 10.1785/0120050609
|
Catherine J K, Gahalaut V K and Sahu V K (2005). Constraints on rupture of the December 26, 2004, Sumatra earthquake from far-field GPS observations. Earth Planet Sci Lett 237: 673-679. doi: 10.1016/j.epsl.2005.07.012
|
Guilbert J, Vergoz J, Schissele E, Roueff A and Cansi Y (2005). Use of hydroacoustic and seismic arrays to observe rupture propagation and source extent of the MW9.0 Sumatra earthquake. Geophys Res Lett 32: L15310, doi: 10.1029/2005GL022966.
|
Ishii M, Shearer M, Houston H and Vidale J (2005). Extent, duration and speed of the 2004-Sumatra-Andaman earthquake imaged by the Hi-net array. Nature 435: 933-936. doi: 10.1038/nature03675
|
Kikuchi M and Kanamori H (1982). Inversion of complex body waves. Bull Seismol Soc Am 72: 491-506.
|
Kruger F and Ohrnberger M (2005). Tracking the rupture of the MW=9.3 Sumatra earthquake over 1150 km at teleseismic distance. Nature 435: 937-939. doi: 10.1038/nature03696
|
Ng W, Reilly J P, Kirubarajan T and Larocque J (2005). Wideband array signal processing using MCMC methods. IEEE Trans Signal Processing 53(2): 471-483. doi: 10.1109/TSP.2004.840673
|
Ni S D, Helmberger and Kanamori H (2005). Energy radiation from the Sumatra earthquake. Nature 434: 582. doi: 10.1038/434582a
|
Ozcetin A (2003). Music, CBF and differential algebraic constant modulus algorithms for direction of arrival estimation in passive coherent locators. Turk J Elec Engin 11(2): 95-115.
|
Vorobyov S A, Gershman A B and Wong K M (2005). Maximum likelihood direction-of-arrival estimation in unknown noise fields using sparse sensor arrays. IEEE Trans Signal Processing 53(1): 471-483.
|
Weinstein E, Feder M and Oppenheim A V (1993). Multichannel signal separation by decorrrelation. IEEE Trans Signal Processing 1(4): 405-413.
|