Your Paper's Title Starts Here: Please Center Estimation Algorithm for Low Frequency Seismic Motions Based on Extended Semblance Function Lianming Sun Faculty of Environmental Engineering, The University of Kitakyushu 1-1 Hibikino, Wakamatsu, Kitakyushu 808-0135, Japan, sun@kitakyu-u.ac.jp Abstract. Monitoring low frequency and long period ground motions is very important for earthquake detection and alarm systems, and signal processing techniques help to investigate the mechanism of deep ground motions and earthquake occurrence by detecting effective seismic information from the seismograms. In this paper a frequency domain algorithm is presented to estimate the travel time difference of seismic waves with multi-path interferences caused by the refraction effect. The spectra of the seismograms are calculated through the fast algorithm, and the multi- path parameters as well as the difference of travel time between a reference position and the seismographic stations are given by an optimization algorithm in the frequency domain. The new approach can work under conditions of refraction interference, and improve the estimation performance by using an extended semblance function. Its effectiveness is demonstrated through some numerical examples, and it is shown that the proposed algorithm is applicable to the analysis of low frequency seismograms. Keywords: Spectral analysis, low frequency seismograms, semblance function, frequency domain. 1. Introduction It has been revealed that the fault activities often accompany by low frequency underground motions, and long period of their accumulation may induce earthquakes. Since many geological phenomena are so complicated that the mechanisms of deep ground tremors and fault slips still remain mysteries, it is difficult to construct an effective physical model for the low frequency motions. Consequently, it is an essential task to monitor the low frequency seismic motions and then to give the alarm of earthquake occurrence. With the significant progress of instrumentation technology, several seismographic networks have been applied to directly monitor the seismic activities in the urgent earthquake detection and alarm systems. Incorporated Research Institutions for Seismology (IRIS) in America, International Data Center managed by Comprehensive Test Ban Treaty (CTBT), Hi-net, Kyoshin-net (K-net) in Japan, are the typical monitoring networks [1]. Using the monitoring networks, the low slip or low-frequency earthquakes, which are considered as the characteristic mode of moment released at a deep structure of a subduction plate interface, are detected from the seismograms [2-5]. Moreover, revealing their mechanisms may help to explore the intrinsic characteristics of seismic phenomena. Nevertheless, the detected signals of low frequency seismograms do not have obvious P-wave and S-wave arrivals; hence it is not easy to locate the epicenters through the ordinary hypocenter determination method [2]. Some array signal processing methods were proposed to maximize a semblance function using a grid search algorithm, then to maximize the cylindrical wave index using steepest descent algorithm [5, 6], or to synchronize the peak of correlation functions of the detected seismic waves in the array records [7]. These methods were performed in time domain, and had been applied in analysis of seismic tremor in Bungo channel region [3], Tokachi channel region [5], Kyoto basin [7], and the long period ground motion of deep ground structure in Yokohama city [8]. Furthermore, the large array scale is often desired in order to acquire much seismic information, whereas a fine grid is required to avoid local minimum in the existing time domain methods. As a result, they lead to considerably heavy computational load. On the other hand, the representation of semblance coefficient has been investigated in the frequency domain, and a frequency domain algorithm has been developed to reduce the computational complexity since some fast frequency algorithm can be utilized [10, 11]. It could achieve the similar accuracy as that of the time domain; moreover, the performance of the presented frequency approach can be improved by compensating the discretization error at low sampling rate. However, when the seismic waves have refraction effects, the multi-path interferences occur and they will shift the maximum of the conventional semblance function due to the side-lobe effects of multi-path. As a result, the estimation accuracy of travel time as well as epicenter location degrades significantly; therefore the effective techniques to reduce the multi- path effects are necessary in the processing algorithms. An extended semblance function is investigated in this paper. By performing subspace decomposition of the spectral data space, the multi-path parameters can be estimated, and the optimization for the extended semblance function yields the estimation of travel time difference between a reference position and tiltmeter stations. It is seen that the multi-path effects can be handled easily in the frequency domain using the estimated parameters, and the estimation of fine travel time difference can improve the location accuracy of epicenter from the seismograms spectra in the proposed algorithm. Some numerical simulations demonstrate the effectiveness of the proposed approach. 2. Problem Statement Consider the seismograms observed by the tiltmeter arrays, which are composed of several stations in a tiltmeter network. If the separation distance between these stations is smaller than the wavelength of the interested low frequency seismograms, the sampling of the seismic waves holds most of the original information with little aliasing. In the seismograms of low frequency earthquake, though the records have not explicit P-wave and S-wave arrivals, they might be analyzed by array processing techniques to detect and locate the low frequency earthquakes, and to utilize the information to reveal the occurrence mechanism and activities in subduction zones. Fig.1 Tiltmeter array in monitoring network Fig.1 illustrates a tiltmeter array in the monitoring network. It has L stations, and each station is equipped with a sensor unit at the bottom of a borehole deeper than 100m [3]. Assume that the separation between stations is smaller than the wavelength of seismic waves to reduce the affection of heterogeneity in crust and mantle. The conventional semblance function )( ref tC [5, 6] is defined by                  K k L l ll K k L l ll tkttaL tktta tC 1 1 2 trv,smpref 1 2 1 trv,smpref ref )( )( )( , (1) where k is the sample number within the time window Kk 1 , tref indicates a reference time, tsmp and tl, trv represent the sampling interval, the travel time difference of the lth station from a reference position, respectively, while the later is a function of the propagation path, wave velocity or the horizontal apparent slowness vector, which is a reciprocal vector of velocity. In [5], the semblance function is maximized by finding an appropriate apparent slowness vector using a grid search algorithm, where the grid is set as small spacing in latitude and longitude, and the estimate of epicenter is determined from the estimated vector of apparent slowness. When the grid search algorithm in the time domain is used simultaneously for multiple arrays, the estimation efficiency will be very poor especially for smaller spacing in both latitude and longitude. An efficient frequency domain approach has been considered for the estimation of travel time difference by using the seismic spectra with low computational load [9]. Define the relation function of the observed signal at stations l1 and l2 as follows: .)()(),( ,    K k llllllll tnkttatnktta K nnR 1 smpsmprefsmpsmpref 22112121 1 (2) Then the numerator and denominator of the semblance function can be rewritten as , , , )( ,, , ,, ,                                                 L l ll ll L l L ll ll ll t t t t RL t t t t R L tC 1 smp trv smp trv 1 1 1 smp trv smp trv ref 1 12 21 21 2 1 (3) where [x] is the nearest integer to the real number x. It is easily seen that maximizing the semblance function (1) is equivalent to maximizing the second term of (3). On the other hand, let the Fourier transform of the observed record be defined by ,)()(     K k ik l i l ff ekttaeA 1 smpref  (4) where Ff f /  is the f th frequency grid, FfF  1 . Then the inverse Fourier transform of )()( * ff i l i l eAeA   21 can be given by ).,()()()( , * , kkkReeAeA KF kX ll F Ff iki l i lll fff     11 1 212121 1  (5) Therefore, the semblance function can be approximated by            L l ll L l L ll llll XL kkX L tC 1 1 1 1 ref 0 2 1 1 12 2121 , , , )( , (6) then )( ref tC in (3) reaches a maximum if the values kl1, l2= kl1- kl2 maximize the following function     ,,maxarg,, 1 1 1 , ,, ,12,1 1 12 2121 ,12,1         L l L ll llll kk LL kkXkk LL  (7) and smptrvtrv 2121 tktt llll ,,,  . If the integer F is the power of 2, both )( fi l eA  and   2121 llll kkX , , can be computed through the algorithm of fast Fourier transform, therefore the algorithm can be performed efficiently [10]. Furthermore, a genetic algorithm based approach can be used to give the estimation of epicenter location [11]. 3. Extended Semblance Function Assume that the observed signal at tiltmeter station l can be approximated as follows:  )()()( ,,,, 1smpref10smpref0smpref lllll kttagkttagktta  , (8) where gl,0, gl,1 and 10 ,, , ll  are the gain and delay time of the direct wave, delayed wave, respectively. )( smpref ktta  is the source wave comes from the epicenter. Then the spectral representation of )( smpref ktta l  can be given by ).()( )()()( , , ff fmlfff ii l i M m i ml K k ik l i l eAeG eAegekttaeA            01 smpref (9) It is seen that 0,trv, ll t  if no delay waves in the observed signal )( smpref ktta l  . Let the extended semblance function be given by                             L l F Ff i l i l i l i l L l L ll F Ff i l i l i l i l f f f f f f f f eG eA eG eA L eG eA eG eA L tC 1 1 * * 1 1 1 1 * * ref )( )( )( )( )( )( )( )( 2 1 )( 1 12 2 2 1 1         , (10) then the parameters of gl,m and ml , might be estimated by optimizing the semblance function in (10). Nevertheless, the direct optimization of (10) is difficult since it is a severe nonlinear problem, and the optimization may converge to local solution easily. In order to improve the processing efficiency, a simplified semblance function will be used in the estimation algorithm. 4. Estimation Algorithm Without loss of the generality, let the reference position be the tiltmeter station 1, and fix g1,0 as g1,0=1, then for l = 2, . . ., L, the travel time difference tl,trv of direct waves from the reference position is 010trvl, ,,   lt . Let the sub- semblance function )( , ref1 tC l be given by     . )()()()()()()()( )()()()()()()()( )( **** **** ,           F Ff ii l i l ii l iii l F Ff i l ii l iii l ii l l ffffffff ffffffff eGeAeAeGeGeAeAeG eGeAeAeGeGeAeAeG tC 1 1111 1 1111 ref1 2 2 1   (11) Commonly only a few frequency components are contained in the seismograms, hence the optimization of (11) will be an ill-conditioned problem, and is influenced by noise easily. Instead the direct optimization, the optimization will be performed in the following procedures. A.Initial Values.Set the initial values gl,0=1, gl,1= gl,2=…=0, then the initial estimates of 010trvl, ,,   lt can be obtained by (7). B.Estimation of Multi-path Parameters.These parameters will be estimated by a subspace decomposition technique in a grid research way. Construct a spectral data matrix ),,,,,( ,,,,  212111 lll kkkk containing both the main lobe and the side-lobe of seismic waves as follows: , )()( )()( )()( )()( ),,,,,( 2, 1, 1 )( 1 )( 1 )( 1 )( )()( 2,1,2,11,1 trv,1,1trv,1,1 trv,0,1trv,0,1 trv,1,11trv,1,11 1 T l l T itkiitki itkiitki i l tkii l tki i l i l lll FllFFllF FllFFllF FlFFlF FF eAeeAe eAeeAe eAeeAe eAeA kkkk                                     Ω Ω Ω            (12) where tl,trv is the estimate given in Step A,  is a grid interval. In order to reduce the influence of noise, the subspace decomposition of (12) is performed as follows: ),,,,,(),,,,,( ,,,,,,,,  212111212111 lllll H l H lll kkkkkkkk ΩΩVSU  , (13) then the noise subspace vectors  V , which are the last column vectors in Vl, can be obtained. Denote 1,V and 2,V as the parts corresponding to 1,l Ω , 2,l Ω , respectively, where 1, V is associated with the parameters g1,m, whereas 2,V is associated with the parameters gl,m. Therefore the optimal values of gl,m , kl,m will be given by   .)(maxarg ,,,,, , ,,,,,,,, ,,,,,,,, ,, ,,,, , ,,        VΩΩVVΩΩVVΩΩV VΩΩVVΩΩV l H l H l H l H l H l H l H l H l H l H kk MlM Ml kkgg g 1 1 1 22221111 11222211 11111 01 11    (14) The first term in (14) indicates an approximation of the extended semblance function )( , ref1 tC l , while the second one is the criterion for multi-path parameter estimation. On the other hand, instead of the extended semblance function, an appropriate pre-filter can be utilized to compensate the multi-path effect by using the estimates of gl,m and kl,m [12]. C.Optimization of Extended Semblance Function.Substitute the estimates of gl,m and kl,m into the extended semblance function )( , ref1 tC l , and leave trvl, t as a variable, then the optimization of )( , ref1 tC l with respect to trvl, t yields fine estimate of the difference of travel time between the tiltmeter stations and reference position. D.Estimation of Epicenter Location.Using the estimates of difference travel time trvl, t , the epicenter location can be estimated by the methods given in [5, 11]. Moreover, in order to decrease the influence of noise, the value of semblance function )( , ref1 tC l can be used as weight function for trvl, t in the epicenter estimation. 5. Numerical Examples Some simulation examples are considered to investigate the effectiveness of the proposed algorithm. The main frequency components are within 0.03~1.5Hz, and the components beyond the frequency band, the interferences caused by refraction or reflection are considered as noise. The seismic wave travels at velocity of 3km/s. 10 stations whose separation are within 18.1km~20.5km are utilized in the tiltmeter array, and 10 simulation runs be performed in the example. The true locations of epicenter in 10 runs are randomly distributed as normal distribution as N(39.0622', (3.2909') 2 ) west in longitude, N(18.3011', (2.6903') 2 south in latitude, as uniform distribution as 10~15km in the depth from the reference station 1. There is a refraction wave in the layer depth of 21km. The seismograms are sampled at the sampling rate 20Hz, and the records in a time window of 105s are used for the data analysis. The simulations are performed under the signal to noise ratio of 12.5, 15, 20dB, respectively. The estimation errors of difference travel time and epicenter location are summarized in Tab.1 and 2. For comparison, the results of conventional methods without compensation of multi-path effects are also given in the tables. It is seen that the proposed algorithm can reduce the estimation errors under the situations with multi-path effects, and it works well especially if a great deal of data can be used in estimation even if the signal to noise ration is low. Tab.1 Standard deviation of difference travel time estimation in 10 simulation runs Method Difference travel time(  s) 12.5dB 15dB 20dB Conventional algorithm 0.03716 0.02809 0.02998 Proposed algorithm 0.02518 0.02325 0.01158 Tab.2 Standard deviation of epicenter location in 10 simulation runs Method Epicenter location (  %) 12.5dB 15dB 20dB Conventional algorithm 3.443 3.119 3.934 Proposed algorithm 2.975 1.751 1.754 6. Conclusions The extended semblance function based algorithm has been proposed to estimate the difference of travel time and the epicenter location for low frequency seismograms. It has been shown that the proposed algorithm can perform optimization in the frequency domain, where the side-lobe can be treated by the shift property of signal spectra. Compared with the conventional methods, the proposed algorithm improves the estimation accuracy of the difference travel time when the seismograms have multi-path effects; therefore it can help to locate the epicenter more precisely. The approach to the case with several epicenters, and the effectiveness validation of the real seismograms will be investigated in the future work. References [1] M. Kikuchi. Realtime Seismology, University of Tokyo Press, 2003 [2] K. Obara. “Nonvolcanic deep tremor associated with subduction in southwest Japan,” Science, vol.269, pp.1679- 1681, 2002. [3] H. Hirose and K. Obara.“Repeating short- and long-term slow slip events with deep tremor activity around the Bungo channel regions, southwest Japan,”Earth Planets Space, vol.57, pp.961-972, 2005. [4] Committee for Earthquake Prediction Studies, Seismological Society of Japan.The Science of Earthquake Prediction, University of Tokyo Process, 2007. [5] Y. Asano, K. Obara and Y. Ito. “Spatiotemporal distribution of very-low frequency earthquakes in Tokachi-oki near the junction of the Kuril and Japan trenches revealed by using array signal processing, ” Earth Planets Space, vol.60, pp.871-875, 2008. [6] N.S. Neidell and M.T.Taner.“Semblance and other coherency measures for multichannel,”Geophysics,vol.36, pp.482-497, 1971. [7] C. Shirakawa and T. Iwata.“Ground motion characteristics of South-East Kyoto Basin site using three dimensional small aperture seismic array at Uji Campus, Kyoto University,”Annuals of Disaster Prevention Research Institute, Kyoto University, no.50(B), pp.251-258, 2007. [8] H. Miura and S. Midorikawa.“Effects of 3-D deep underground structure on characteristics of rather long-period ground motion,” Earthquake, no.54, pp.381-395, 2001. [9] W. Zhu, L. Sun and X. Zhu.“Analysis of low frequency seismograms in frequency domain,”ICIC Express Letters, vol.5, no.9(A), pp.3249-3254, 2011. [10] L. Sun and A.Sano.“System identification and error analysis in frequency domain,”International Journal of Innovative Computing, Information and Control, vol.3, no.5, pp.1201-1218, 2007. [11] W. Zhu, L. Sun and X. Zhu.“Epicenter location via genetic algorithm for low frequency seismograms, ICIC Express Letters, vol.6 (to appear) [12] G. Jacovitti and G. Scarano. “Discrete time techniques for time delay estimation,” IEEE Trans. Signal Processing, vol.41, no.2, pp.525-533, 1993.