- Email: [email protected]
Plane wave fitting method for a plane, small aperture, short period seismic array: a MATHCAD program
Computers & Geosciences 28 (2002) 59–64
Short Note
Plane wave fitting method for a plane, small aperture, short period seismic array: a MATHCAD program$ Edoardo Del Pezzo*, Flora Giudicepietro INGV-Osservatorio Vesuviano, Via Diocleziano 328, 80124 Napoli, Italy Received 8 January 2001; received in revised form 26 June 2001; accepted 30 June 2001
1. Introduction A seismometer array is a dense network of seismic sensors deployed in the field with a geometrically regular configuration. It conceptually acts like a directional antenna for radio waves, in the sense that a seismometer array measures the directional properties of the wavefield (slowness or wavenumber vector) radiated by a seismic source. Arrays are used in seismology in different aspects of seismic monitoring, spanning from the detection and identification of underground nuclear explosions to the wavefield study at short distance of weak seismic events like, for example, volcanic tremor or volcanic earthquakes. The recent improvement of the seismic instrumentation and of the digital data loggers allowed a fast development of portable arrays, which actually can be easily deployed for experiments of short time duration like the study of seismic noise for seismic engineering microzonation purposes, or for experiments aimed at the detection and tracking of the tremor and long period (LP) events on active volcanoes. In particular, tremor and LP’s are considered as possible precursors of impending eruptions, and for this reason the array techniques are becoming important tools in volcano monitoring (Chouet, 1996). The array techniques can be easily divided into two main categories: those working in time domain of the seismic signals and those in frequency (or spectral) domain. Among the spectral methods, MUSIC approach (Goldstein and Chouet, 1994) is able to resolve multiple closely spaced sources, being selectively sensi$
Code on server at http://www.iamg.org/CGEditor/ index.htm *Corresponding author. Fax: +390816108351. E-mail addresses: [email protected] (E. Del Pezzo), [email protected] (F. Giudicepietro).
tive to the strongest source. On the other hand MUSIC, as all the other approaches based on Fourier transforming the time signal recorded at the array sensors, fails when it is necessary to study short duration phases composed of a few periods. In this situation, the methods based on the time-domain estimate of the cross-correlation function between the signal pairs are the most appropriate: they can be, in principle, applied to signals of every duration, having the only defect of being more time-consuming than those in frequency domain. Among the methods in time domain the so-called ‘‘zero lag cross-correlation‘‘ (ZLCC) method (Frankel et al., 1991) seems appropriate to the use for small time window signals, as recently shown by Saccorotti and Del Pezzo (2000). As ZLCC is based on a grid search of the optimal solution, it is quite time-consuming. For a quicker use in the field during short duration field experiments, a method based on the search of the minimum of a misfit function between the observed time delays at the array station pairs and the theoretical delays expected for a plane wavefront of given slowness vector, can be used as a simplification of the ZLCC technique. This technique is generally called ‘‘plane wave fitting technique’’ (PWF). The observed time delays can be obtained by a visual picking of the phase of interest, or, more precisely, by an estimate of the time value for which the cross-correlation function of signal pairs takes its maximum. In the present paper, we outline the PWF method; then we will describe the MATHCADTM professional program we have developed (PWFP) and finally show an example of the application to real data taken during an experiment carried out at Deception volcano, Antarctica. Due to the ease of understanding MATHCADTM programs, even by a non-specialized public, the example could also be used for didactic or illustrative purposes of this kind of technique.
0098-3004/02/$ - see front matter r 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 3 0 0 4 ( 0 1 ) 0 0 0 7 6 - 0
60
E. del Pezzo, F. Giudicepietro / Computers & Geosciences 28 (2002) 59–64
2. The method
3. The Mathcad program PWFP
Given an array of seismic stations in a flat area located at points with Cartesian coordinates xi ; yi having an almost regular geometrical configuration. In the hypothesis that a wave packet of a given slowness vector p ðpx ; py Þ crosses the array, with the inter-distances between the station pairs sufficiently small to permit a good degree of correlation among the waveforms, the travel time differences between station i and j are given by the dot product
The program developed to perform the calculations described previously is reported in Appendix 1, available on the IAMG server. It contains several plots, internally generated by the MATHCAD graphical system, used to illustrate the calculations, and is widely commented. PWFP reads and plots the true array coordinates, import the signals recorded at the array stations as ASCII files, each containing the digitized data samples for a single array station. Graphically the user can select the portion of the signal which needs to be treated. The evaluation of the numerical cross-correlation between the station pairs is carried out using P the internal MATHCAD routine correlðx; yÞi ¼ k xk yiþk : The
tij ¼ p Dxij ;
ð1Þ
where Dxij is the two-element vector of the differences between the coordinates of station i and j: Dxij ðxi xj ; yi yj Þ as also explained in the comment beside the definition of DXX array in Appendix 1 (available on the IAMG server). tij can be estimated through the cross-correlation function between the signal pairs ai ðtÞ and aj ðtÞ Z þN ai ðtÞaj ðt þ tÞ dt ð2Þ Cij ðtÞ ¼ N
whose first maximum occurs at a time, tmax ; which corresponds to the phase delay between signals at station i and j: If we denote as t the one-column matrix of the NðN 1Þ=2 travel time differences between the possible station pairs (tij ¼ tji ) we can write t ¼ p . Dx;
ð3Þ
where Dx is the two-columns matrix containing the NðN 1Þ=2 row elements Dxij : Eq. (3) can be viewed as an overdetermined system of NðN 1Þ=2 linear equations in the two unknowns px and py and solved using the least-squares approach: p ¼ ðDxT DxÞ1 DxT t
ð4Þ
where superscript T indicates the transpose of a matrix. Azimuth f and apparent velocity v ¼ 1=jpj can be obtained by Eq. (5) v ¼ ðp2x þ p2y Þ1=2 ð5Þ f ¼ tan1 ðpx =py Þ: Errors can be obtained by the covariance matrix of the data, covðtÞ; expressed by Menke (1984) T covðpÞ ¼ ðDxT DxÞ1 Dx covðtÞ ðDxT DxÞ1 Dx : ð6Þ It is interesting to note that, assuming an a priori reasonable covariance matrix of data, Eq. (6) can easily be used for experimental design of the array geometry prior to the deployment of the stations in the field.
Fig. 1. (a) Plot of synthetic signal (in arbitrary units) used for testing MATHCAD procedure. Horizontal axis represents time (s); (b) plot of slowness components (x and y axes represent respectively, x and y components of slowness vector in s/m) calculated for test with synthetic data. Cross represents the solution of PWF technique with error bars, and box is the true slowness value.
E. del Pezzo, F. Giudicepietro / Computers & Geosciences 28 (2002) 59–64
61
Fig. 2. Map of Deception volcano with spatial configuration of two arrays deployed during experiment. Geometrical configuration of ‘‘Fumarolas’’ array, represented by A in upper panel, is represented enlarged in right panel.
output is ordered using the internal function recenter which puts the lag-0 element in the middle of the sequence. Even though in principle a precision greater than the uncertainty introduced by the time sampling of the signal can be reached using a robust interpolator of the numerical array generated by the correl algorithm, we assumed an uncertainty of 1 sampling time interval for the evaluation of the time corresponding to the maximum of the cross-correlation function. Matrix operations are internally optimized in MATHCAD. In particular, the inversion of a matrix is carried out using the singular value decomposition algorithm (Menke, 1984).
4. A test with synthetic data Synthetic signals SðtÞ having the characteristics of the true signals which can be recorded on active volcanoes have been generated by Almendros et al. (1999) using the following formula: B t SðtÞ ¼ A et=t0 sinð2pf0 tÞ; t0
ð7Þ
where A is the maximum amplitude at time t0 ; f0 is the predominant frequency and B takes into account the sharpness of the signal envelope from the start time to the time corresponding to the maximum amplitude. In Fig. 1(a), a plot of SðtÞ is reported, calculated using A ¼ 10; B ¼ 1:6; t0 ¼ 0:3 s; f0 ¼ 3 Hz. The test signal is contaminated with a Gaussian noise with zero mean and standard deviation of 0.2. Then the data set is built making the single trace be shifted according to a test slowness vector psynth ; given as input parameter. The synthetic traces are finally used as input for the program PWFP . The solution pest is plotted together with psynth in Fig. 1(b). Appendix 2 (on the server) contains the run of the PWFP program applied to the synthetic signal.
5. Application to real data Two seismic arrays operated at Deception Island, an active volcano located in Antarctica, in 1998–1999 austral summer during a 3-month experiment (Fig. 2) aimed at the study of wavefield and source properties of the short period seismic radiation generated by volcanic and hydrothermal activity. During the experiment an intense swarm (up to 3000 earthquakes of Magnitude ranging from -1 to 3.0) of volcano-tectonic events was recorded at the two seismic arrays, besides the background hydrothermal activity which characterizes this active volcano. A detailed description of this experiment can be found in Saccorotti et al. (2001). The volcanotectonic earthquakes of the swarm have been studied using the zero-lag cross-correlation (ZLCC) technique to calculate apparent slowness and back-azimuth of the seismic phases related to the first P-wave pulse. A description of this method is in Frankel et al. (1991) and Del Pezzo et al. (1997). Almendros et al. (1999) modified the technique allowing circular wave fronts, and Saccorotti and Del Pezzo (2000) studied the error function associated to the grid search estimate. The ZLCC technique provides a complete 3-D image of the error function, but is slow, being the search of the best solution based on a grid search. Del Pezzo et al. (1997) and La Rocca (1999) compare ZLCC with the techniques working in the frequency domain as MUSIC (Goldstein and Chouet, 1994) or high resolution and beam-forming (Capon, 1979), showing that ZLCC produce the same results if applied to high frequency (> 1 Hz) seismic signals of short duration. In the present paper we, apply PWF technique to five earthquakes belonging to the seismic swarm described previously in order to test PWFP with real data. We check the results obtained with the PWF technique comparing them with those obtained with ZLCC. In Fig. 2 a sketch map of Deception Island with the array positions is reported together with the configuration of the array we used for the present example. An
62
E. del Pezzo, F. Giudicepietro / Computers & Geosciences 28 (2002) 59–64
Fig. 3. Plot of four examples: two vertical lines confine time window selected for analysis.
example of the recorded data is reported in Fig. 3, where we plot the seismograms of 4 of the 5 selected earthquakes. In Appendix 1 (on server) the MATHCAD sheet shows the application of the PWF method to the last example. In Fig. 4 the PWF solution is plotted together with the ZLCC solution. There is agreement between the two solutions which have superimposed error bounds. The length of the error interval (not reported in the plot) for ZLCC technique is of the same order of that obtained using PWF. It must be also noted
that due to the high apparent velocity of the wavefront incoming to the array, the examples reported in the present paper represent the minimum resolution for both techniques.
6. Discussion and conclusions The PWF method is the fastest technique working in time domain. It has the advantage of being based on the
E. del Pezzo, F. Giudicepietro / Computers & Geosciences 28 (2002) 59–64
63
Fig. 4. Plot of slowness components (x and y axes represent, respectively, the x and y components of slowness vector in s/m) calculated for four earthquakes recorded at array A whose seismograms are plotted in Fig. 3. Crosses are PWF solutions with error bars. Squares are ZLCC solutions.
solution of a linear overdetermined system of equations which can be performed using a least-squares approach. In this approach, the unit-covariance of the solution permits the a priori experimental design of the array geometry, being independent of the data covariance. Conversely, the advantage of using ZLCC is that the average cross-correlation function is obtained and plotted for the whole slowness space, giving a complete picture of the error field. MATHCADTM allows the use of symbolic routines (based on MAPLE V) together with numerical procedures. It is based on Windows 98 operating system and takes the advantage of Microsoft’s OLE 2 object linking to work with other Windows 98 applications. Its graphic system is easy to use and at the same time sufficiently powerful to visualize all the program steps allowing a satisfactory control of the program functioning. Conversely, MATHCADTM
sometimes suffers from crashing problems related to the Windows 98 operating system and for many executions may be slower than other programming languages. The advantage of MATHCADTM ; a standard for technical professionals, in seismological use lies in its user-friendliness, which favors the diffusion and the understanding of the details of many algorithms in a wide non-specialist community, including students. In the seismological community the analysts of the seismic laboratories may take advantage of using MATHCADTM for solving routine problems as for example magnitude evaluation, spectral analysis and other applications. In the field, the easy use of PWFP in MATHCADTM environment by the field personnel (including non-expert students) may be of great help in experiments carried out using seismic arrays. To conclude, the easy understanding through the
64
E. del Pezzo, F. Giudicepietro / Computers & Geosciences 28 (2002) 59–64
MATHCADTM interface of the structure of the PWFP program may find favor even among non-expert users in making improvements or changes. For example, even though it has been written in the approximation of a plane wavefront, the program structure can be easily modified to allow for a circular wavefront (Almendros et al., 1999). Note. Appendix 1 and 2 are available from IAMG server at http://www.iamg.org/CGEditor/index.htm.
Acknowledgements The authors thank G. Saccorotti, J.M Ibanez, J. Almendros, M. La Rocca and M. Martini for continuing support and discussion. J. Ibanez kindly provided the data used as an example of application together with the ZLCC solutions. The paper has been partly supported by the Antartic ANT98-1111 project of the Spanish Government and partly by Civil Protection of Italy, Progetti Speciali legge 74.
References Almendros, J., Ibanez, J., Alguacil, G., Del Pezzo, E., 1999. Array analysis using circular-wave-front geometry: an application to locate the nearby seismo-volcanic source. Geophysical Journal International 136, 159–170.
Capon, J., 1979. Maximum likelihood spectral estimation. In: Kaykin, S. (Ed.). Nonlinear Methods of spectral analysis. Springer Verlag, Berlin. Chouet, B., 1996. New methods and future trends in seismological volcano monitoring. In: Scarpa, R., Tilling, R.I. (Eds.), Monitoring and Mitigation of Volcano Hazards. Springer, Berlin, pp. 23–97. Del Pezzo, E., La Rocca, M., Ibanez, J., 1997. Observation of high-frequency scattered waves using dense arrays at Teide volcano. Bulletin Seismological Society of America 87 (6), 1637–1647. Frankel, A., Hough, S., Friberg, P., Busby, R., 1991. Observations of Loma Prieta aftershocks from a dense array in Sunnyvale, California. Bulletin of the Seismological Society of America 80, 1900–1922. Goldstein, P., Chouet, B., 1994. Array measurements and modeling of sources of shallow volcanic tremor at Kilauea volcano, Hawaii. Journal of geophysical Research 99, 2637–2652. La Rocca, M., 1999. Multichannel array techniques applied to volcanic hazard. Ph.D. Dissertation, Faculty of Sciences, University of Salerno, Salerno, Italy. Menke, W., 1984. Geophysical Data Analysis: Discrete Inverse Theory. Academic Press, New York. Saccorotti, G., Almendros, J., Carmona, E., Ibanez, J., Del Pezzo, E., 2001. Slowness anomalies from two dense arrays at Deception Island, Antarctica. Bulletin of the Seismological Society of America, 91, 557–571. Saccorotti, G., Del Pezzo, E., 2000. A probabilistic approach to the inversion of data from a seismic array and its application to volcanic signals. Geophysical Journal International 143, 1–18.
Short Note
Plane wave fitting method for a plane, small aperture, short period seismic array: a MATHCAD program$ Edoardo Del Pezzo*, Flora Giudicepietro INGV-Osservatorio Vesuviano, Via Diocleziano 328, 80124 Napoli, Italy Received 8 January 2001; received in revised form 26 June 2001; accepted 30 June 2001
1. Introduction A seismometer array is a dense network of seismic sensors deployed in the field with a geometrically regular configuration. It conceptually acts like a directional antenna for radio waves, in the sense that a seismometer array measures the directional properties of the wavefield (slowness or wavenumber vector) radiated by a seismic source. Arrays are used in seismology in different aspects of seismic monitoring, spanning from the detection and identification of underground nuclear explosions to the wavefield study at short distance of weak seismic events like, for example, volcanic tremor or volcanic earthquakes. The recent improvement of the seismic instrumentation and of the digital data loggers allowed a fast development of portable arrays, which actually can be easily deployed for experiments of short time duration like the study of seismic noise for seismic engineering microzonation purposes, or for experiments aimed at the detection and tracking of the tremor and long period (LP) events on active volcanoes. In particular, tremor and LP’s are considered as possible precursors of impending eruptions, and for this reason the array techniques are becoming important tools in volcano monitoring (Chouet, 1996). The array techniques can be easily divided into two main categories: those working in time domain of the seismic signals and those in frequency (or spectral) domain. Among the spectral methods, MUSIC approach (Goldstein and Chouet, 1994) is able to resolve multiple closely spaced sources, being selectively sensi$
Code on server at http://www.iamg.org/CGEditor/ index.htm *Corresponding author. Fax: +390816108351. E-mail addresses: [email protected] (E. Del Pezzo), [email protected] (F. Giudicepietro).
tive to the strongest source. On the other hand MUSIC, as all the other approaches based on Fourier transforming the time signal recorded at the array sensors, fails when it is necessary to study short duration phases composed of a few periods. In this situation, the methods based on the time-domain estimate of the cross-correlation function between the signal pairs are the most appropriate: they can be, in principle, applied to signals of every duration, having the only defect of being more time-consuming than those in frequency domain. Among the methods in time domain the so-called ‘‘zero lag cross-correlation‘‘ (ZLCC) method (Frankel et al., 1991) seems appropriate to the use for small time window signals, as recently shown by Saccorotti and Del Pezzo (2000). As ZLCC is based on a grid search of the optimal solution, it is quite time-consuming. For a quicker use in the field during short duration field experiments, a method based on the search of the minimum of a misfit function between the observed time delays at the array station pairs and the theoretical delays expected for a plane wavefront of given slowness vector, can be used as a simplification of the ZLCC technique. This technique is generally called ‘‘plane wave fitting technique’’ (PWF). The observed time delays can be obtained by a visual picking of the phase of interest, or, more precisely, by an estimate of the time value for which the cross-correlation function of signal pairs takes its maximum. In the present paper, we outline the PWF method; then we will describe the MATHCADTM professional program we have developed (PWFP) and finally show an example of the application to real data taken during an experiment carried out at Deception volcano, Antarctica. Due to the ease of understanding MATHCADTM programs, even by a non-specialized public, the example could also be used for didactic or illustrative purposes of this kind of technique.
0098-3004/02/$ - see front matter r 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 3 0 0 4 ( 0 1 ) 0 0 0 7 6 - 0
60
E. del Pezzo, F. Giudicepietro / Computers & Geosciences 28 (2002) 59–64
2. The method
3. The Mathcad program PWFP
Given an array of seismic stations in a flat area located at points with Cartesian coordinates xi ; yi having an almost regular geometrical configuration. In the hypothesis that a wave packet of a given slowness vector p ðpx ; py Þ crosses the array, with the inter-distances between the station pairs sufficiently small to permit a good degree of correlation among the waveforms, the travel time differences between station i and j are given by the dot product
The program developed to perform the calculations described previously is reported in Appendix 1, available on the IAMG server. It contains several plots, internally generated by the MATHCAD graphical system, used to illustrate the calculations, and is widely commented. PWFP reads and plots the true array coordinates, import the signals recorded at the array stations as ASCII files, each containing the digitized data samples for a single array station. Graphically the user can select the portion of the signal which needs to be treated. The evaluation of the numerical cross-correlation between the station pairs is carried out using P the internal MATHCAD routine correlðx; yÞi ¼ k xk yiþk : The
tij ¼ p Dxij ;
ð1Þ
where Dxij is the two-element vector of the differences between the coordinates of station i and j: Dxij ðxi xj ; yi yj Þ as also explained in the comment beside the definition of DXX array in Appendix 1 (available on the IAMG server). tij can be estimated through the cross-correlation function between the signal pairs ai ðtÞ and aj ðtÞ Z þN ai ðtÞaj ðt þ tÞ dt ð2Þ Cij ðtÞ ¼ N
whose first maximum occurs at a time, tmax ; which corresponds to the phase delay between signals at station i and j: If we denote as t the one-column matrix of the NðN 1Þ=2 travel time differences between the possible station pairs (tij ¼ tji ) we can write t ¼ p . Dx;
ð3Þ
where Dx is the two-columns matrix containing the NðN 1Þ=2 row elements Dxij : Eq. (3) can be viewed as an overdetermined system of NðN 1Þ=2 linear equations in the two unknowns px and py and solved using the least-squares approach: p ¼ ðDxT DxÞ1 DxT t
ð4Þ
where superscript T indicates the transpose of a matrix. Azimuth f and apparent velocity v ¼ 1=jpj can be obtained by Eq. (5) v ¼ ðp2x þ p2y Þ1=2 ð5Þ f ¼ tan1 ðpx =py Þ: Errors can be obtained by the covariance matrix of the data, covðtÞ; expressed by Menke (1984) T covðpÞ ¼ ðDxT DxÞ1 Dx covðtÞ ðDxT DxÞ1 Dx : ð6Þ It is interesting to note that, assuming an a priori reasonable covariance matrix of data, Eq. (6) can easily be used for experimental design of the array geometry prior to the deployment of the stations in the field.
Fig. 1. (a) Plot of synthetic signal (in arbitrary units) used for testing MATHCAD procedure. Horizontal axis represents time (s); (b) plot of slowness components (x and y axes represent respectively, x and y components of slowness vector in s/m) calculated for test with synthetic data. Cross represents the solution of PWF technique with error bars, and box is the true slowness value.
E. del Pezzo, F. Giudicepietro / Computers & Geosciences 28 (2002) 59–64
61
Fig. 2. Map of Deception volcano with spatial configuration of two arrays deployed during experiment. Geometrical configuration of ‘‘Fumarolas’’ array, represented by A in upper panel, is represented enlarged in right panel.
output is ordered using the internal function recenter which puts the lag-0 element in the middle of the sequence. Even though in principle a precision greater than the uncertainty introduced by the time sampling of the signal can be reached using a robust interpolator of the numerical array generated by the correl algorithm, we assumed an uncertainty of 1 sampling time interval for the evaluation of the time corresponding to the maximum of the cross-correlation function. Matrix operations are internally optimized in MATHCAD. In particular, the inversion of a matrix is carried out using the singular value decomposition algorithm (Menke, 1984).
4. A test with synthetic data Synthetic signals SðtÞ having the characteristics of the true signals which can be recorded on active volcanoes have been generated by Almendros et al. (1999) using the following formula: B t SðtÞ ¼ A et=t0 sinð2pf0 tÞ; t0
ð7Þ
where A is the maximum amplitude at time t0 ; f0 is the predominant frequency and B takes into account the sharpness of the signal envelope from the start time to the time corresponding to the maximum amplitude. In Fig. 1(a), a plot of SðtÞ is reported, calculated using A ¼ 10; B ¼ 1:6; t0 ¼ 0:3 s; f0 ¼ 3 Hz. The test signal is contaminated with a Gaussian noise with zero mean and standard deviation of 0.2. Then the data set is built making the single trace be shifted according to a test slowness vector psynth ; given as input parameter. The synthetic traces are finally used as input for the program PWFP . The solution pest is plotted together with psynth in Fig. 1(b). Appendix 2 (on the server) contains the run of the PWFP program applied to the synthetic signal.
5. Application to real data Two seismic arrays operated at Deception Island, an active volcano located in Antarctica, in 1998–1999 austral summer during a 3-month experiment (Fig. 2) aimed at the study of wavefield and source properties of the short period seismic radiation generated by volcanic and hydrothermal activity. During the experiment an intense swarm (up to 3000 earthquakes of Magnitude ranging from -1 to 3.0) of volcano-tectonic events was recorded at the two seismic arrays, besides the background hydrothermal activity which characterizes this active volcano. A detailed description of this experiment can be found in Saccorotti et al. (2001). The volcanotectonic earthquakes of the swarm have been studied using the zero-lag cross-correlation (ZLCC) technique to calculate apparent slowness and back-azimuth of the seismic phases related to the first P-wave pulse. A description of this method is in Frankel et al. (1991) and Del Pezzo et al. (1997). Almendros et al. (1999) modified the technique allowing circular wave fronts, and Saccorotti and Del Pezzo (2000) studied the error function associated to the grid search estimate. The ZLCC technique provides a complete 3-D image of the error function, but is slow, being the search of the best solution based on a grid search. Del Pezzo et al. (1997) and La Rocca (1999) compare ZLCC with the techniques working in the frequency domain as MUSIC (Goldstein and Chouet, 1994) or high resolution and beam-forming (Capon, 1979), showing that ZLCC produce the same results if applied to high frequency (> 1 Hz) seismic signals of short duration. In the present paper we, apply PWF technique to five earthquakes belonging to the seismic swarm described previously in order to test PWFP with real data. We check the results obtained with the PWF technique comparing them with those obtained with ZLCC. In Fig. 2 a sketch map of Deception Island with the array positions is reported together with the configuration of the array we used for the present example. An
62
E. del Pezzo, F. Giudicepietro / Computers & Geosciences 28 (2002) 59–64
Fig. 3. Plot of four examples: two vertical lines confine time window selected for analysis.
example of the recorded data is reported in Fig. 3, where we plot the seismograms of 4 of the 5 selected earthquakes. In Appendix 1 (on server) the MATHCAD sheet shows the application of the PWF method to the last example. In Fig. 4 the PWF solution is plotted together with the ZLCC solution. There is agreement between the two solutions which have superimposed error bounds. The length of the error interval (not reported in the plot) for ZLCC technique is of the same order of that obtained using PWF. It must be also noted
that due to the high apparent velocity of the wavefront incoming to the array, the examples reported in the present paper represent the minimum resolution for both techniques.
6. Discussion and conclusions The PWF method is the fastest technique working in time domain. It has the advantage of being based on the
E. del Pezzo, F. Giudicepietro / Computers & Geosciences 28 (2002) 59–64
63
Fig. 4. Plot of slowness components (x and y axes represent, respectively, the x and y components of slowness vector in s/m) calculated for four earthquakes recorded at array A whose seismograms are plotted in Fig. 3. Crosses are PWF solutions with error bars. Squares are ZLCC solutions.
solution of a linear overdetermined system of equations which can be performed using a least-squares approach. In this approach, the unit-covariance of the solution permits the a priori experimental design of the array geometry, being independent of the data covariance. Conversely, the advantage of using ZLCC is that the average cross-correlation function is obtained and plotted for the whole slowness space, giving a complete picture of the error field. MATHCADTM allows the use of symbolic routines (based on MAPLE V) together with numerical procedures. It is based on Windows 98 operating system and takes the advantage of Microsoft’s OLE 2 object linking to work with other Windows 98 applications. Its graphic system is easy to use and at the same time sufficiently powerful to visualize all the program steps allowing a satisfactory control of the program functioning. Conversely, MATHCADTM
sometimes suffers from crashing problems related to the Windows 98 operating system and for many executions may be slower than other programming languages. The advantage of MATHCADTM ; a standard for technical professionals, in seismological use lies in its user-friendliness, which favors the diffusion and the understanding of the details of many algorithms in a wide non-specialist community, including students. In the seismological community the analysts of the seismic laboratories may take advantage of using MATHCADTM for solving routine problems as for example magnitude evaluation, spectral analysis and other applications. In the field, the easy use of PWFP in MATHCADTM environment by the field personnel (including non-expert students) may be of great help in experiments carried out using seismic arrays. To conclude, the easy understanding through the
64
E. del Pezzo, F. Giudicepietro / Computers & Geosciences 28 (2002) 59–64
MATHCADTM interface of the structure of the PWFP program may find favor even among non-expert users in making improvements or changes. For example, even though it has been written in the approximation of a plane wavefront, the program structure can be easily modified to allow for a circular wavefront (Almendros et al., 1999). Note. Appendix 1 and 2 are available from IAMG server at http://www.iamg.org/CGEditor/index.htm.
Acknowledgements The authors thank G. Saccorotti, J.M Ibanez, J. Almendros, M. La Rocca and M. Martini for continuing support and discussion. J. Ibanez kindly provided the data used as an example of application together with the ZLCC solutions. The paper has been partly supported by the Antartic ANT98-1111 project of the Spanish Government and partly by Civil Protection of Italy, Progetti Speciali legge 74.
References Almendros, J., Ibanez, J., Alguacil, G., Del Pezzo, E., 1999. Array analysis using circular-wave-front geometry: an application to locate the nearby seismo-volcanic source. Geophysical Journal International 136, 159–170.
Capon, J., 1979. Maximum likelihood spectral estimation. In: Kaykin, S. (Ed.). Nonlinear Methods of spectral analysis. Springer Verlag, Berlin. Chouet, B., 1996. New methods and future trends in seismological volcano monitoring. In: Scarpa, R., Tilling, R.I. (Eds.), Monitoring and Mitigation of Volcano Hazards. Springer, Berlin, pp. 23–97. Del Pezzo, E., La Rocca, M., Ibanez, J., 1997. Observation of high-frequency scattered waves using dense arrays at Teide volcano. Bulletin Seismological Society of America 87 (6), 1637–1647. Frankel, A., Hough, S., Friberg, P., Busby, R., 1991. Observations of Loma Prieta aftershocks from a dense array in Sunnyvale, California. Bulletin of the Seismological Society of America 80, 1900–1922. Goldstein, P., Chouet, B., 1994. Array measurements and modeling of sources of shallow volcanic tremor at Kilauea volcano, Hawaii. Journal of geophysical Research 99, 2637–2652. La Rocca, M., 1999. Multichannel array techniques applied to volcanic hazard. Ph.D. Dissertation, Faculty of Sciences, University of Salerno, Salerno, Italy. Menke, W., 1984. Geophysical Data Analysis: Discrete Inverse Theory. Academic Press, New York. Saccorotti, G., Almendros, J., Carmona, E., Ibanez, J., Del Pezzo, E., 2001. Slowness anomalies from two dense arrays at Deception Island, Antarctica. Bulletin of the Seismological Society of America, 91, 557–571. Saccorotti, G., Del Pezzo, E., 2000. A probabilistic approach to the inversion of data from a seismic array and its application to volcanic signals. Geophysical Journal International 143, 1–18.