A Method to Estimate Spatial Resolution in 2-D Seismic Surface Wave Tomographic Problems ()
1. Introduction
In the recent years, the general aspects of the lithospheric shear wave velocity structure of the South America continent has been studied via surface waves [1] -[4] . In some cases, the investigation was concentrated in specific areas of the continent, such as central and western South America [1] . In this study, they used waveform of Rayleigh waves to image the 3-D S wave velocity structure of the region. The general features of four slices (at 100, 150, 200 and 300 km depth) of the proposed model for the region were presented and discussed. The spatial resolution of this study was obtained through resolution tests (based essentially on synthetic data with random noise) and it was qualitatively classified into very good (to depths above 300 km) and none to good (to depths below to 300 km). In other study in the South America continent [2] , they used both Rayleigh (from 10 to 150 s) and Love waves (from 20 to 70 s) group-velocities to analyze the 3-D S-velocity structure of the region. In this case, the region was discretized into cells of 2˚ × 2˚ and the problem was solved via conjugate-gradient method. A regionalized dispersion curve at each cell of the grid was then used to estimate the 1-D S-velocity in depth. The spatial resolution was obtained through both checkerboard and realistic structure tests. According to [2] , the spatial resolution at 30 km is around 4˚, while at 100 km it is 6˚ and at 150 km it is 8˚. The general features of the S wave velocities anomalies at those depths were presented and discussed. More recently, [3] used joint inversion of waveforms and fundamental mode Rayleigh wave group-velocities to investigate the South America continent. The methodologies used by both [1] [2] were used in this study to get group-velocities and for waveforms inversions, respectively. The spatial resolution was obtained with checkerboard tests and according to the authors it was 6˚ at both 100 and 300 km depths. The S wave velocities anomalies to these depths were presented and discussed. The most recent study about the S wave velocity structure beneath South America continent was presented by [4] . He used a group of 128 earthquakes occurred at western South America continent and recorded by 71 seismic stations distributed throughout continental part of the region. Rayleigh wave group-velocities varying from 5 to 250 s were used to obtain the source-station dispersion curves. The region was divided into cells of 4˚ × 4˚ and a regionalized Rayleigh wave dispersion curve was obtained for each cell of the grid. Each regionalized dispersion curve was inverted to get a 1-D S wave velocity in each cell. None quantitative evaluation about spatial resolution of the estimated shear wave velocity in the region was presented. Only a qualitative evaluation of the resolution (i.e., the terms bad resolution and poor resolution were applied to classify the resolution in depth) was presented and discussed. In the present study, a mathematical approach to estimate quantitatively the spatial resolution in 2-D seismic surface wave problems is proposed and it is applied in the northeastern Brazilian region.
2. The Two-Dimensional Surface Wave Tomographic Problem
Let’s consider an area, defined by a pair of geographycal coordinates, which is discretized into square cells. Suppose that there are 
source-station surface wave paths crossing this area and also assume that for a particular path, the total source-station surface wave travel time is the sum of the travel time in each cell. Thus, for an 
path and a given period
, one can write [5]
(1)
where 
is an attributed phase and/or group velocity for the 
cell, 
is the source-station phase and/or group velocity for the 
path calculated from attributed
, 
is the observed source-station phase and/or group velocity for the 
path, 
is the phase and/or group velocity for the 
cell, 
is the length of the 
path in the 
cell, and 
is the source-station distance of the 
path.
The result of the procedure aforementioned is a dispersion curve (phase and/or group velocity as a function of period) for each cell of the target area. It is also important to remember that only the cells visited by at least one surface wave path are considered in such computation. Equation (1) can be rewritten in a matrix form as
(2)
where 
is a residual vector (m × 1), 
is a vector with the unknown parameters (n × 1) and 
is a mxn matrix relating model parameters to observations. In terms of damped least-squares [6] [7] , Equation (2) is given by
(3)
where 
is the transpose of matrix
, 
is the Levenberg-Marquardt parameter, and 
is the identity matrix. In this context, the expressions for resolution 
and covariance 
matrices are the following
(4)
(5)
where 
is the variance in the observations. The lines of the resolution matrix 
are commonly called of resolving kernels. The standard deviations of the estimated group-velocities at each cell are obtained in the diagonal elements of the covariance matrix. Furthermore, several factors can contribute to uncertainties in group-velocity estimation [8] . These include origin time error, epicenter mislocation error, source finiteness error, rise time error, etc. These errors were estimated and incorporated to the 
term of the covariance matrix
.
3. Estimation of the Spatial Resolution
The best and more accurate way of getting the resolution of the estimated phase and/or group velocities in any place of the earth’s planet, i.e., the spatial extent of the smallest phase and/or group velocity anomaly that can be imaged, is through the computation of resolution matrix. However, the computation of this matrix is not a simply task, because it involves operations with very large matrices. For this reason, the major part of these kinds of studies does not compute the resolution matrix. In general, those studies use a different approach, called checkerboard tests, which are based on the inversion of synthetic data constructed as from subjective assumptions about input and output models [9] According to [9] , in such analysis the quality criterion is the similarity between both final and initial models used to calculate the synthetics.
In this study, we calculated the exact resolution matrix (Equation (4)), by using some efficient computational techniques into a Personal Computer [10] . It is fundamental to observe that the computation of the resolution matrix through Equation (4) is based only on the surface wave paths distributions of the target area.
According to [11] , if the resolution has a single sharp maximum centered about the main diagonal, thus the data are well resolved. However, if the resolution is very broad, thus the data are poorly resolved (Section 4.2 of [11] ). In spite of these important definitions, none way of quantifying the resolution is presented in the seismological literature.
In this study, a novel methodology to quantify the spatial resolution in 2-D seismic surface wave tomographic problems is proposed. It is based on both the resolving kernels obtained via computation of the full resolution matrix (Equation (4)) and a basic property of a Gaussian function, because a resolving kernel shape can be approximated to a Gaussian function. In fact, as 
in Equation (4) is always a positive number, then the lines of 
assume a maximum value in the position of the unknown parameter under evaluation (i.e. phase and/or group velocity at a particular cell) and the neighborhoods points tend to zero. A graphical representation of such resolving kernel shows that it follows a Gaussian function. Thus, by using this assumption and the concept of Full Width at Half Maximum (FWHM) of a Gaussian function, which represents the difference between the two values of the independent variable at which the dependent variable is equal to half of its maximum value [12] . In spectroscopy, the concept of FWHM of a Gaussian function is applied in the definition of spectrometers’ resolution [13] -[14] .
According to [15] , the FWHM of a Gaussian function is given by
(6)
where 
is the standard deviation. Thus, let’s consider a Gaussian function described by the following expression
(7)
where 
is the independent variable and 
was defined previously. Then, rewriting Equation (7) in a condensed way
(8)
where 
is equal to 
and
(9)
and
(10)
Applying the natural logarithm in both sides of Equation (8)
(11)
The Equation (11) is a second order polynomial function or a quadratic function, which is represented graphically by a parabola. In this case, a quadratic least-squares fit can be used to get the three coefficients of the quadratic function and, consequently, the parameter
.
Extracting s from Equation (10), we get
(12)
Substituting Equation (12) in Equation (6), one gets an expression for the width of the Gaussian at half maximum 
as a function only of the parameter c, which is given by
(13)
Equation (13) provides a quantitative spatial resolution measurement of the unknown parameter at any cell of the grid under investigation. Here, the unknown parameter is the estimated phase and/or group velocity either at a specific cell or at a particular period (Equation (1)). As the periods of the surface waves dispersion curves are related to depths (i.e. short periods are associated with shallow depths while long periods are related to deep parts of the earth’s interior), then the Equation (13) also allows the evaluation of the spatial resolution in depth.
The Equation (4) shows that the resolution matrix is a square matrix, i.e., if the number of unknown parameters is
, then 
has 
elements. As each element of a resolving kernel is related to the phase and/or group velocity at a given cell, then each element of 
matrix has a spatial location, because it is represented by the geographical coordinates. Furthermore, if a resolving kernel has a shape similar to a Gaussian function [11] , thus the maximum value of it should be located at the cell (i.e. the center of it) under investigation and it should decay to zero around a central value. I am interested in analyzing the spatial shape of the resolving kernel around a specific cell, i.e., the spread of the resolution or spatial decay around a central point, which is located in the diagonal of the resolution matrix [11] .
4. Application to the Northeastern Brazilian Region
The region limited by the geographical coordinates (−125, 45) and (25, −75) is, initially, divided into square cells of 2˚ × 2˚ (Figure 1). To have a good path coverage of the target area (i.e., northeastern Brazil), a large quantity of earthquakes recorded at twenty three seismic stations (from 1988 to 1998 [16] ) of the IRIS Consortium [17] , located either in or around South America continent, were selected in this study.
The digital seismograms, related to seismic events with focal depth shallower than 120 km and magnitude higher than 5.0 mb, were requested and received from IRIS. The Seismic Analysis Code (SAC) package [18] was used to the identification of the surface waves in the seismic signals. This analysis produced 3134 digital seismograms, corresponding to 1082 earthquakes (Figure 1). The vertical components of the motion were used here and they represent 3134 source-station Rayleigh wave paths in the area displayed in Figure 1. To remove the effect of the seismograph on the seismic signal, an instrumental correction was applied in the seismograms (via SAC). Correction for source group delay was not applied, because it is considered a small source of error [19] . The Multiple-Filter Technique [20] was applied in the digital seismograms selected for this study. Twenty two periods, ranging from 10 to 102 s (i.e., 10.04, 12.05, 14.03, 16.00, 18.29, 20.08, 24.38, 28.44, 32.00, 36.57, 42.67, 46.55, 51.20, 56.89, 60.24, 64.00, 68.27, 73.14 78.77, 85.33, 93.09, 102.40 s), were selected to compose the source-station dispersion curves. This computation produced 3134 source-station Rayleigh wave group velocity dispersion curves and all details about the surface wave data processing can be found in [16] .
The surface wave tomography methodology, described in Section 2, was applied in data set aforementioned to get the Rayleigh wave dispersion curves of all cells (visited by the surface waves) in Figure 1. As my goal is to
Figure 1. Seismic stations (red triangles) and epicenters of the earthquakes (blue circles) used in the present study. The area is limited by the geographical coordinates (−125, 45) and (25, −75). The heavy square box represents northeastern Brazilian region.
investigate northeastern Brazilian region, thus we selected only the Rayleigh wave dispersion curves associated with each cell of such region (Figure 2). For each period, its source-station Rayleigh wave path distributions were plotted into maps to give an idea about their spatial distribution (Figures 3-7).
5. Discussion of the Results
The width of the Gaussian at half maximum (W—Equation (13)) was used to get the spatial resolution of the estimated Rayleigh wave group-velocities at each cell of the northeastern Brazil. As the region under investigation is formed by 49 cells (Figure 2) and each cell has a dispersion curve (which is composed of 22 periods), I got 1078 resolution estimates (49 × 22) for the northeastern Brazilian region. All those results are condensed in Table1 For a given cell of the northeastern Brazilian region, I got five elements of the resolution matrix (Equation (4)). The central point (third element) is related to the center of the cell under evaluation and the others four points are located in the centers of the neighborhood cells (two on left and two on right). Here, the separation between each central point is equal to 2˚, because the studied area was divided into cells of 2˚ × 2˚, thus the separation between the five points is equal to 2˚.
A few cases displayed in Table 1 are used in Figures 8-10 to show the procedure to computation of the spatial resolution at each cell in the northeastern Brazil. The five elements of the resolving kernel (centered at a particular cell and for a specific period) are submitted to a Gaussian fit to calculate the corresponding coefficients. By using the parameter c (from Equation (11)), the spatial resolution at a given cell is calculated through Equation (13). The black circles are the five elements of the resolution matrix (i.e., the resolving kernel computed from Equation (4)) for a specific cell and period. They are located at 2, 4, 6 (the third element or central point), 8 and 10 positions. As mentioned previously, they are separated by 2˚, because of the cell size in the 2-D inversion
Figure 2. The grid of the Figure 1 was divided into cells of 2˚ × 2˚ and each cell on it was numbered (from top to bottom and from left to right) so that any region inside northeastern Brazil must be associated with the numbers shown on the map.
Figure 3. Source-station Rayleigh wave path distributions (for periods of 10.04 and 12.05 s) in the area of the northeastern Brazil. The black lines are the great-circle paths connecting epicenters and stations. The target area is relatively well covered by the Rayleigh wave path distribution. A clear decrease in the quantity of the source-station paths, when the period increases, is observed in these maps.
process. The red square is the theoretical Gaussian function calculated with the Gaussian fit coefficients. In the upper right side is shown the width of the Gaussian at half maximum (W—Equation (13)) and the correlation coefficient (CC) between the resolving values (computed via Equation (4)) and the theoretical Gaussian function calculated with the Gaussian fit coefficients. An analysis of the correlation coefficients values displayed in Table 1
Figure 4. Source-station Rayleigh wave path distributions (for periods of 16.00 and 20.08 s) in the area of the northeastern Brazil. The black lines are the great-circle paths connecting epicenters and stations. The target area is relatively well covered by the Rayleigh wave path distribution. A clear decrease in the quantity of the source-station paths, when the period increases, is observed in these maps.
Figure 5. Source-station Rayleigh wave path distributions (for periods of 24.38 and 32.00 s) in the area of the northeastern Brazil. The black lines are the great-circle paths connecting epicenters and stations. The target area is relatively well covered by the Rayleigh wave path distribution. A clear decrease in the quantity of the source-station paths, when the period increases, is observed in these maps.
shows that 70% of them are equal and higher than 0.7. The large quantity of high CC values indicates that Gaussian function seems to be a good representation of the resolving kernels’ shape.
The results displayed in Table 1 also show some interesting features. First, the spatial resolution varies from 3.2˚ to 12.7˚, but the major part of them is concentrated between 4.0˚ and 8.0˚, i.e., from two to four cells. This spatial resolution range is in agreement with all previous tomographic and 3-D shear wave velocities studies developed in the South American continent ([21] —8˚; [22] —from 6˚ to 9˚; [23] —5˚; [2] —from 4˚ to 8˚; [3] —6˚).
Figure 6. Source-station Rayleigh wave path distributions (for periods of 46.55 and 68.27 s) in the area of the northeastern Brazil. The black lines are the great-circle paths connecting epicenters and stations. The target area is relatively well covered by the Rayleigh wave path distribution. A clear decrease in the quantity of the source-station paths, when the period increases, is observed in these maps.
Figure 7. Source-station Rayleigh wave path distributions (for periods of 85.33 and 102.40 s) in the area of the northeastern Brazil. The black lines are the great-circle paths connecting epicenters and stations. The target area is relatively well covered by the Rayleigh wave path distribution. A clear decrease in the quantity of the source-station paths, when the period increases, is observed in these maps.
Table 1 also shows that the group-velocity anomalies smaller than 3.2˚ cannot be resolved in northeastern Brazil and, therefore, those anomalies should not be considered in the spatial resolution analysis. It is important to emphasize that the procedure developed here allows analyzing any group velocity anomaly in any particular point of the area studied.
Despite the fast reduction in the Rayleigh wave paths quantity, when the period increases (Figure 3), the spatial resolution value, for a given cell, varies smoothly and it is clearly observed in Table1 In general, the spatial resolution values (for a particular cell) show either a slightly increase in their magnitude (which suggest a
Figure 8. Examples of spatial estimates (at cells 1764 and 1842; and periods 10.04 and 18.29 s, respectively) displayed in Table1 The black circles are the five elements extracted from the resolution matrix (Equation (4)) at a particular cell (central element— position 6) and period. The red squares are the theoretical values computed with the coefficients of the Gaussian fit. In the upper left side, the legends of the symbols are: res-resolving values from resolution matrix; and gau-theoretical Gaussian values. In the upper right side (and for a specific cell), one can find both the estimated spatial resolution value (W) and the Correlation Coefficient (CC) between the estimated resolving kernels via resolution matrix (Equation (4)) and the theoretical values calculated with the coefficients of the Gaussian fit. In all cases, the theoretical Gaussian curve is very close to the estimated resolution value obtained with the resolution matrix (Equation (4)). This is also confirmed by the large quantity of high CC values in Table1 It shows that a Gaussian function represents very well the resolving kernel’s shape.
Figure 9. Examples of spatial estimates (at cells 1916 and 1993; and periods of 32.00 and 51.20 s, respectively) displayed in Table1 The black circles are the five elements extracted from the resolution matrix (Equation (4)) at a particular cell (central element— position 6) and period. The red squares are the theoretical values computed with the coefficients of the Gaussian fit. In the upper left side, the legends of the symbols are: res-resolving values from resolution matrix; and gau-theoretical Gaussian values. In the upper right side (and for a specific cell), one can find both the estimated spatial resolution value (W) and the Correlation Coefficient (CC) between the estimated resolving kernels via resolution matrix (Equation (4)) and the theoretical values calculated with the coefficients of the Gaussian fit. In all cases, the theoretical Gaussian curve is very close to the estimated resolution value obtained with the resolution matrix (Equation (4)). This is also confirmed by the large quantity of high CC values in Table1 It shows that a Gaussian function represents very well the resolving kernel’s shape.
Figure 10. Examples of spatial estimates (at cells 2070 and 2142; and periods of 73.14 and 102.40 s, respectively) displayed in Table1 The black circles are the five elements extracted from the resolution matrix (Equation (4)) at a particular cell (central element— position 6) and period. The red squares are the theoretical values computed with the coefficients of the Gaussian fit. In the upper left side, the legends of the symbols are: res-resolving values from resolution matrix; and gau-theoretical Gaussian values. In the upper right side (and for a specific cell), one can find both the estimated spatial resolution value (W) and the Correlation Coefficient (CC) between the estimated resolving kernels via resolution matrix (Equation (4)) and the theoretical values calculated with the coefficients of the Gaussian fit. In all cases, the theoretical Gaussian curve is very close to the estimated resolution value obtained with the resolution matrix (Equation (4)). This is also confirmed by the large quantity of high CC values in Table1 It shows that a Gaussian function represents very well the resolving kernel’s shape.
decrease in a spatial resolution) or an almost constant throughout the entire period range (or depth range) displayed in Table1 In some cases (cells 1840, 2141 and 2142, Table 1), the spatial resolution values increase when the period increases. In these cases, we can observe two main regions on Table 1 (i.e., periods lower than 46.55 s and periods higher than 46.55 s). In the first region (periods £ 46.55 s), the spatial resolution values are higher than those for the second region (periods > 46.55 s). As shown in Table 1 and expected, it seems that azimuthal distribution is a parameter more important than the quantity of the surface wave paths. These results indicate that it is possible to get a better spatial resolution by using a reduced quantity of seismic data.
Another important point observed in the proposed methodology is related to the improvement of the spatial resolution in surface wave tomography studies. In general, in the major part of the surface wave studies in the seismological literature, the improvement of spatial resolution is associated with the increase of data quantity. In fact, the improvement of the spatial resolution should be related to the cell size. Let’s consider the case of a photo. In order to increase the resolution of a photo we must decrease the size of the discretization element, i.e., in this case, we increase the number of dots per inch (or increase the density of dots and, consequently, decrease the dot size). Thus, to get a better spatial resolution we must either decrease the cell size or increase the number of source-station paths with a good azimuthal coverage. The previous discussions show that the proposed methodology describes the main features of the seismological problem.
6. Conclusion
A novel methodology to quantify the spatial resolution in 2-D seismic surface wave tomographic problems is proposed. It is based on both the resolving kernels computed via full resolution matrix and the concept of Full Width at Half Maximum (FWHM) of a Gaussian function. The method allows estimating quantitatively the spatial resolution at any cell of a gridded area. The spatial resolution range estimates with the proposed method in northeastern Brazil is in agreement with those obtained with all previous surface wave investigations in the South America continent.
Acknowledgements
This research was developed by using data from GSN and GTSN seismic networks available at IRIS. The author would like to thank these institutions for providing all seismic data. He also thanks Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ) for financial support during the development of this study (Process #: E-26/170.407/2000) and MCTI/Observatório Nacional for all support. Several figures were prepared with the software Generic Mapping Tools [24] .