Natural-neighbor Isosismals

Bulletin of the Seismological Society of America, Vol. 92, No. 5, pp. 1933–1940, June 2002

Natural-Neighbor Isoseismals by Livio Sirovich, Fabio Cavallini, Franco Pettenati, and Muzio Bobbio

Abstract Natural-neighbor (n-n) isoseismals are proposed as a new tool that solves the centennial problem of drawing objective and reproducible isoseismals from earthquake damage sparsely observed in a region (felt reports). The algorithm uses the n-n coordinates for weighting, interpolating, and contouring the felt reports. In our computer implementation, at each step, the surface of irregularly distributed observations is partitioned into a unique set of Voronoi polygons computed on a fine regular grid. The interpolation is local, because the weight of an experimental site brought to a new neighbour point is proportional to the area of the intersection of their Voronoi polygons. In the n-n approach, the interpolant (a) fits the data exactly at the observation sites; (b) is isoparametric and bounded by the data values; and (c) is continuously differentiable at all points, except the data sites. Moreover, the n-n isoseismals do not increase the complexity of the quantitative geophysical interpretation because they do not introduce new (contouring) parameters; and, finally, they may be intersected automatically with geological and topographical information. The new natural-neighbor isoseismals appear as a happy compromise between the crude objectivity of the Voronoi tessellation and the intuitive appeal of the somewhat subjective classical isoseismals. Introduction The problem of drawing objective and reproducible isoseismals is a centennial matter in seismology. It came to light from the earliest trials by Egen (1828) and Mallet (1862; see the chapter “On the forms and areas, of the meizoseismal and isoseismal curves, and the position of the seismic vertical therein,” Vol. 2, 252–258). Even earlier attempts to trace isoseismals date back to an earthquake of the seventeenth century in the Neapolitan kingdom (De Poardi, 1627). Notwithstanding the length of time elapsed since then, no reliable solution has been found; in fact, already in 1921 Davison noticed the contradictory contours obtained by various experts (Davison, 1921; p. 124). Since then, the subjectivity of hand-drawn isoseismals and the absence of an international protocol for drawing these lines automatically were stressed by many authors (e.g., more recent works by Bollinger, 1977; Hanks and Johnston, 1992; Pettenati el al., 1999). The ambiguity of results was also demonstrated by an international experiment (Cecic’ et al., 1996). We also argued that tracing isoseismals from irregularly spaced point observations, using the techniques available at the time, was an ill-posed problem (Pettenati et al., 1999). In particular, we stressed that (a) the sampling obtained from the survey does not often obey the spatial version of the Nyquist principle along the coordinate axes (Press et al., 1992, 494–495); (b) the regional intensity pattern is the result of the superposition of regional continuous effects (radiation from the source, path) and of the local effects of

geological and topographical heterogeneities, which appear to be of discontinuous nature because of the insufficient sampling; and (c) the combined process of sampling plus contouring in automatic procedures constitutes a two-dimensional filter (Wren, 1975), which results in implicit gratuitous physical assumptions about the regional distribution of damage. In spite of all this, in some cases, an easy-to-grasp contouring would help the analyst (for example, in relating semiquantitatively active faults to preinstrumental earthquakes). Thus, the need for a straightforward drawing tool remains. Moreover, if these pictures are based upon a sufficient number of surveyed sites per unit area (frequency completeness), such graphical solutions could help the analyst in grasping some essence of some phenomena intuitively as well. Think, for example, of the meizoseismal areas of strong transcurrent earthquakes, which are elongated in the sense of the fault, or think of the more compact, but often rather irregular, aspect of the meizoseismal areas of earthquakes provoked by dip-slip mechanisms. In our opinion, however, the technique to be used must honor the observed data. To fulfill this condition, consider that contouring is an all-pass filter only if, when resampling the isoseismals at the same locations where intensity was observed in the field, one obtains results that coincide with the observations used. We will show that our new n-n isoseismals meet this condition.

1933

1934

L. Sirovich, F. Cavallini, F. Pettenati, and M. Bobbio

-119 -118

S3 c

b

S2

h

g

X

a

f e

d

S4

S1

34 MM Scale 4

S5

5 6 7

Figure 1.

Local interpolation scheme of the natural-neighbor (n-n) approach for weighting the intensity at the new point X from the values observed at sites Si.

8 9

Figure 2.

The new n-n isoseismals are produced by our program ConVor (Contouring by Voronoi), which is based on the n-n approach (Sibson, 1980; Watson, 1992; Braun and Sambridge, 1995; Sambridge et al., 1995) and on the Voronoi tessellation. Note that we already chose the Voronoi technique when inverting intensity data for the source (Pettenati et al., 1999; Sirovich and Pettenati, 2002). For example, Figure 2 shows the macroseismic data set of the Northridge 1994 earthquake after tessellation by Voronoi polygons. Finally, the possibility of deriving information on the source from the regional intensity patterns has also established itself in seismology (Atkinson and Boore, 1998; for more related references, see Sirovich et al., 2001b). Consider, in this regard, that our n-n isoseismals do not increase the complexity of quantitative source-inversion because they do not increase the number of unknowns by the contouring parameters.

Theory and Algorithms Recall that the Voronoi tessellation solves the problem of proximity in the plane so that 1. Each surveyed site Si, where i ⳱ 1, . . . , N (see Fig. 1), is circumscribed by a convex polygon TSi (thick segments in the figure). 2. Any point p inside each polygon TSi is closer to Si than to any other Sj; that is, TSi ⳱ {p 僆 R : d(p, Si)
d(p, Sj), for all j ⬆ i},

Northridge 1994 earthquake. Observed Modified Mercall; intensities (courtesy of J. Dewey, U.S. Geological Survey, written personal communication, 1994) after tessellation with Voronoi polygons. Dots indicate surveyed sites.

3. The nearest-neighbor search for tracing the segments connecting adjacent sites may be efficiently solved by the Delaunay triangulation (thin lines in Fig. 1). See details and theorems in Preparata and Shamos (1985) and the implementation of the n-n algorithm by Sambridge et al. (1995) at the Website http://rses.anu.edu.au/ jean/PROJECTS/NN_QUICK/index.html. The Voronoi polygons allow one to correlate observed intensities with their geographical information density (Pettenati et al., 1998), but may give a not completely satisfactory sketch of the intensity field in areas with few observations. Note, in particular, the peripheral polygons in Figure 2. For the reader’s convenience, we now briefly recall the basic ideas underlying the n-n method, which is based on a multiple application of the Voronoi tessellation. For a comprehensive treatment, we refer to Sibson (1980) and Sambridge et al. (1995), who gave full details of the theory and analytical solutions of the n-n method to locally interpolate intensity in a new point X (see Fig. 1). In this interpolation scheme, the locally interpolated intensity is a function of observed intensities Ii : I(X) ⳱

2

(1)

where R2 is the two-dimensional Cartesian space and d is the Euclidean distance between points.

兺Ni⳱1wi(X) • Ii ,

(2)

area (TX 艚 TSi) . area (TX)

(3)

where wi(X) ⳱

1935

Natural-Neighbor Isoseismals

Here, referring to Figure 1, we call TX the new Voronoi polygon a-b-c-d-e generated by X. TX intersects some of the polygons TSi. For example, a-e-f-g-h is TX 艚 TS1. The interpolated intensity I(X) is uniquely determined by the data set {(Si, Ii): i ⳱ 1, . . . , N}, and honors the data because wi (Si) ⳱ 1 and wi (Sj) ⳱ 0 if i ⬆ j. This surface is continuous; its derivative is also continuous, except at Si (see the proof in Sambridge et al., 1995). For simplicity, we implemented this interpolation scheme (2) with approximate graphical means, using the screen pixel as the infinitesimal surface element (Okabe et al., 2000, 51). The logical scheme is as follows.

isoseismals of the Northridge 1994 earthquake hand-traced by Dewey et al. (1995) according to USGS practice. Figure 3 is an example of what we are used to seeing in the best international practice. In a previous paper (Pettenati et al., 1999), we noticed that the stochastic contouring technique called inverse weighted-distance average (IWDA) (Eckstein, 1989, equation A.6) is one of the best candidates to substitute for subjective practices in the light of its mean value property (Appendix A in Pettenati et al., 1999); in fact, this property allows IWDA to minimize some drawbacks of other competing interpolants, for example, by separating two neighboring intensity fields leaving as few points as possible on the wrong side. The IWDA technique gave better fits than the cubic spline (Eckstein, 1989) and the kriging interpolants (Philip and Watson, 1986). Indeed, the IWDA was adopted by the USGS for contouring the amplification values within the Los Angeles Basin (Hartzell et al., 1996). Thus, we chose to compare the contour obtained by the n-n scheme with that obtained by IWDA (see Figs. 4a–c and 5). In Figures 4–6 we used all the original USGS data within the figures Ⳳ 0.2 in longitude and Ⳮ 0.3ⳮ0.2 in latitude; the higher northward extension was necessary due to the scarcity of sites in that direction. Thus, a total of 171 data were used. No cosmetic smoothing was used in Figures 4–5. Like most algorithms of this kind, the IWDA technique needs a case-by-case tuning: Figures 4a–c are three examples of the many different isoseismals that could be produced by this technique by tuning the seven contour parameters differently. In particular—according to our subjective ex-

1. Produce the TSi Voronoi polygons and intersections TXi 艚 TSi for all points Xi of a regular grid. 2. Count the pixels of TXi and of TXi 艚 TSi. 3. Compute the interpolated field I(Xi) using equations (2) and (3). 4. Compute and display the 2D density plot of I(X) by cubic piecewise surfaces using the values of I computed at each node. Moreover, one may display a contour plot: isoseismals are back again.

Results Figure 3 shows the macroseismic data points in the Modified Mercalli Intensity scale of 1931, MMI (Wood and Neumann, 1931), obtained from felt reports and from field surveys (courtesy of the U.S. Geological Survey (USGS) at Denver; J. Dewey, written communication, 1994) and the

119.0

118.0

118.5

Palmdale

V

34.5 Fillmore

VII

Sta . C l a r i t a

Sta . Susana Mtns.

San G ab rie

San Fer n a ndo

IX

VIII

l M tns.

V

C h a tswor t h

IX

Oxnard

Gl e n d a l e

Tarazana

VI

Santa Monic a Mt ns.

VIII

34.0

Sta . M o n i c a

VII Los Angeles Whittier

VI Redondo B e a c h

Anaheim

Epicenter

9 Felt at Intensity IX F Felt

Long Beach 20

0 Km

Figure 3. Hand-traced isoseismals of the Northridge, 1994 earthquake according to the USGS practice (Dewey et al., 1995, Fig. 2; courtesy of the USGS).

1936

L. Sirovich, F. Cavallini, F. Pettenati, and M. Bobbio

A.

VII

VII

VI

-118

-118

VI

-119

VIII

B.

-119 VI

VII

VII

VIII

VIII

VII

I VII

VI

IX VI VII

VII

II VI

VII

34

34

5

6

6

VI

5

7

VI

7

VI

MM S c a l e 4

MM S c a l e 4

8

8

9

9

C.

-118 VI

-119VI VII

VIII

VI

VI I

I VII VII

I VII

34 MM S c a l e 4

VII

VI

5 6 7 8 9

Figure 4. Three examples of automatic contours of the data of Fig. 3 by the inverse weighteddistance average (IWDA) stochastic technique. (a) subjectively tuned at best; (b) contour parameters as in Fig. 4a, but with e ⳱ ⳮ1; (c) contour parameters as in Fig. 4a, but with e ⳱ ⳮ2 and r ⳱ 30 km (see text). perience and opinion—Fig. 4A is the best achievable IWDA result when one wants to honor the abrupt local intensity changes and to best reproduce the meizoseismal area of MMI 9 without producing artifacts and spikes. In Figure 4a, the 171 data were automatically interpolated by using the fol-

lowing parameters: maximum search radius r ⳱ 20 km; minimum search radius rmin ⳱ 2 m; distance exponent for interpolation weighting e ⳱ ⳮ4; maximum number of experimental intensities considered in all directions nmax ⳱ 36; minimum number per quadrant nmin ⳱ 2; 200 ⳯ 180 inter-

1937

Natural-Neighbor Isoseismals

-119

VI

-118

VI

VII

VI I

IX

II VI

VIII

VI I

VI

34 VI

MM S c a l e 4 5 7

VI

6 8 9

Figure 5.

Data as in Figs. 3–4, but n-n isoseismals.

polation grid; null smoothing. If nmin
2 in a quadrant, then r overtakes the limit of 20 km to fulfill the nmin ⱖ 2 condition. The dispersion (calculated intensity minus observed) in Figure 4a is low because the data set of this earthquake is relatively simple and interpolation is very local (short r); in general, the IWDA contours of more irregular points give a higher dispersion. Regarding this, note the higher dispersion and some peculiarities of figures 4b–c. Figure 4b was obtained with the same contour parameters as Figure 4a, but with e ⳱ ⳮ1. It is worth commenting, for example, that Figure 4b does not honor many intensity points. Note, in particular, the two anomalous data at Saugus (MMI 4; 34.41138, 118.53917; see the northern black arrow in the figure) and at the Los Angeles–I10/Fairfax site (MMI 9; 34.04500, 118.36000; see the southern black arrow in the figure) which is very close to an MMI 6 site; also consider that the former site is a statistical outlier (Pettenati et al., 1999). Then, in the same figure, the IWDA contour is able to trace only two very small areas of degree IX (the two areas are almost invisible in the figure, being close to the tips of the two white arrows in Fig. 4b). Also, in the southeastern quadrant of the figure, many sites with degree 5 are included in the area of degree VI, some of them being randomly honored by microzones of degree V and some not. Figure 4c was also obtained with the same contour parameters as Figure 4a, but with e ⳱ ⳮ 2 and r ⳱ 30 km; this figure honors the data, but it is unable to circumscribe the six epicentral sites with MMI 9 by one close curve delimiting a continuous meizoseismal area; regarding this, and from the point of view of the spatial density of information, its con-

tinuity is objectively indicated by the fact that in Fig. 2 there are six adjacent Voronoi polygons with MMI 9. Finally, keeping the Northridge data set fixed, some contour parameters give highly unstable IWDA results and more artifacts (data not shown). Figure 5 shows the n-n isoseismals of the same earthquake; it is worth recalling that the isolines shown in the figure are uniquely determined: the algorithm is deterministic and has no adjustable parameters. A 200 ⳯ 180 interpolation grid was used in Figure 5. We stress that in the n-n approach the grid density is the only limit for graphically honoring the most abrupt observations; however, this density is a mere graphical choice external to the contour and, thus, the number of nodes per axis is not a contour parameter, as it is in other algorithms. So, given the all-pass filter nature of the n-n isoseismals, there is no numerical dispersion in Figure 5. In conclusion, this figure was obtained without any explicit or implicit assumption about the observed phenomenon, and no contouring parameters were introduced. We also performed an experiment to analyze the robustness and stability of the n-n isoseismals and of the optimized IWDA intensity isolines in Figure 4a. Consider that (a) IWDA was optimally tuned specifically for the case under study; thus, (b) the test conditions are in its favor; and (c) if one had used Figure 4b instead, the performance of IWDA would have been very poor. For this comparison, we prepared 38 intensity data sets, each containing 164 data; each one of the 38 data sets was obtained by discarding 4% of the original USGS data set of 171 values; the jackknife technique was used for this (Barnett and Lewis, 1978, 47–48). We then performed 38 series of interpolations on regular grids by using the IWDA technique and the optimal contour parameters of Figure 4a. Then, in each node of the grid, we calculated the root-mean-square error (rms) of the intensity residuals between the 38 values obtained in the same node by using the 38 trimmed sets, minus the value predicted in the same node by using the complete USGS set of 171 data (see Fig. 4a). This procedure was also repeated in the case of the n-n algorithm (program ConVor), and the rms values thus obtained are shown in Fig. 6b. We also checked to see whether 38 trimmed sets were sufficient to obtain stable RMS, and whether the RMS were biased by the order of the used trimmed sets. With this in mind, we randomly chose four series of N data sets from the 38 trimmed sets, with N going from 2 to 38, and calculated the chi-square of each series. The results of this analysis are shown in the insets of Figures 6a (IWDA algorithm) and 6b (n-n algorithm); the chi-square values obtained from the four series of trimmed data sets become stable above N ⳱ 10– 20 in ConVor and above N ⳱ 20–30 in IWDA. Thus, we conclude that Figures 6a and 6b are reliable. Next, note that the chi-square values in the inset of Figure 6b are approximately 10 times lower than the corresponding values in the inset of Figure 6a. Note also that in both Figures 6a and 6b the maximum rms values are due to the two aforementioned anomalous data, at Saugus (MMI 4) and at Los Angeles–

1938

L. Sirovich, F. Cavallini, F. Pettenati, and M. Bobbio

B

0.1

0.1

0.2

0.1

0.2

0.3

0.6

0.1

0.4

0.1

0.2

0. 1

A 0.1

0.4

0.1

50

40

6

0.4

5

C hi-square

1 0.

4 3

30

0.1

C hi-square

0.7

0.2

0.1

0.1

0.1

0.1

0.1

2

20 0

10

20

30

40 N

1 0

10

20

30

40 N

Figure 6. (a) geographical distribution of rms intensity given by the IWDA technique on 38 N-7 trimmed data sets; (b) as in (a) but by the n-n approach. (For insets, see text).

T10/Fairfax (MMI 9). The isoseismals of Figure 6b are robust because these two points perturb the n-n isoseismals only locally. On the other hand, the Saugus datum largely perturbs Figure 6a. The rms values in Figure 6b are lower than in 6a and involve smaller areas; thus, the n-n isolines are more stable than IWDA isolines. Notwithstanding our efforts to tune IWDA properly, this algorithm was unable to reproduce the meizoseismal area of degree 9 at best, without losing control of Figure 6a in the northwest corner (where rms is greater than 2 degrees; see also Fig. 4a); this is due to the different densities of observation points in the two zones.

Discussion and Conclusions The Voronoi polygons lead to an objective representation, which depends only on the values of the field data and on their spatial distribution, and can be profitably used for inversion (Sirovich and Pettenati, 1999); however, the graphics so obtained have a somewhat unnatural look, because the Voronoi tessellation obviously does not allow an intensity variation within a given polygon. This feature is more evident where the data are scarce. On the other hand, the n-n representation in Figure 5 is objective, honors the data and—although ultimately based on Voronoi tessellation—has a more natural look (for example, because it is consistent with the fact that, at a regional scale, intensity usually decreases with increasing distance). The IWDA algorithm can give acceptable contours, but it must be tuned for each individual case. Kriging is consid-

ered the best linear unbiased estimator, and is also able to honor data rather well, but the individual fit of the semivariogram is crucial (see the classic text by Ripley, 1981). The triangulation technique honors the data completely, but the contour obtained is unsatisfactory (data not shown). In general, the use of the standard automatic techniques is critical and contradictory in the zone southeast of the meizoseismal area in Figures 3–5, where surveyed sites with MMI 5–6 are distributed in spots (see the details in Figs. 1a and 1b of Pettenati et al., 1999). On the other hand, the n-n isoseismals adapt to the complicated nature of intensity data sets, which are the result of the superposition of regional effects (radiation from the source, and path) and of local effects (geological and topographical heterogeneities). Strictly speaking, both effects result from continuous phenomena, but at a regional scale local heterogeneities appear to be of a discontinuous nature because of the (unavoidably) insufficient sampling provided by macroseismic surveys. The n-n approach is a straightforward technique for getting an objective picture of this kind of data (Sambridge et al., 1995). It makes a totally local interpolation that honors the empirical evidence completely. Strictly speaking, however, it does not overcome all the criticisms that we mentioned in the Introduction about the illposed approach of the current isoseismals. In fact, in case of abrupt intensity changes, the n-n isoseismals give the picture an aspect that is probably smoother than most situations in nature are. This is an unavoidable problem, however. We stress that the n-n isoseismals include or exclude a site on the basis of an exclusively objective criterion of in-

1939

Natural-Neighbor Isoseismals

formation continuity in the field. Look again, for example, at the intensities distributed in spots in the southeastern corners of Figures 2 and 5, and appreciate how inclusions and exclusions of sites from certain isoseismals are simply driven by the continuity or discontinuity between areas determined by the Voronoi polygons in Figure 2. Thus, the n-n isoseismals may also be used to objectively evidence areas with more or less abrupt, positive or negative, increments of MMI; the search for extended site effects, or for errors, can then be addressed there. For this purpose, the n-n isoseismals may also be intersected automatically with geological and topographical information. We stress that the isoseismals in Figure 5 are produced without any subjectivity, any implicit physical assumption, and without the need for any contouring parameter. An alternative would be to draw intensity data sets simply as points, possibly of different colors or diameters, but it is very difficult for the human eye to synthesize a global pattern of intensity from these scattered data alone. In this sense, the Voronoi tessellation is useful, but not entirely satisfactory. The n-n isoseismals allow one to better guess at the shape of the intensity patterns and to direct successive quantitative analyses. In fact, they help guess whether the analyzed data set appears as a more or less continuous field or whether it is made up of a series of data in spots (which, as noted, could be due to local effects, errors, etc.). It is a commonplace that any two-dimensional correlation ought to be preliminarily perceived by the human eye. By using the new n-n isoseismals, the preliminary and qualitative stage of the observation of the intensity patterns is rendered more objective through a reproducible protocol. We argue that the proposed n-n isoseismals have the qualities to be taken into consideration as an international standard for graphically representing regional intensities. Can the n-n isoseismals be used also for other quantitative applications? For example, estimating the magnitude of preinstrumental earthquakes from intensity data is a wellestablished (and unavoidable) practice. Normally, this is done by using the dimensions or mean radii of equal intensity areas. The choice of areas and radii however, is biased by the many kinds of contouring techniques used. Our technique can help define areas and mean radii objectively. We successfully tested ConVor in the case of the 1925 Charlevoix earthquake in Canada (Cavallini et al., 2000) and the 1987 Whittier Narrows earthquake in California (Pettenati and Sirovich, 2000). From a more general point of view, we suggest using n-n isolines of continuous- or discretevalued fields to produce objective and reproducible pictures that are easy to grasp and reliable. In summary, the new n-n isoseismals are a straightforward method obtained from geometry, and no ambiguities or tunable parameters are involved; they appear to be a happy compromise between the crude objectivity of the Voronoi tessellation and the intuitive appeal of the somewhat subjective classical isoseismals. Of course, extensive exper-

imentation with field data is needed to substantiate this view definitely.

Acknowledgments Supported by GNDT of CNR and of the Italian Ministry of Civil Defence, Grants Number 97.00536.PF54 and 98.03227.PF54. We are grateful to James W. Dewey, of the U.S. Geological Survey at Denver, who provided the computer files of intensity data points. L. Sirovich dealt mainly with seismological aspects, F. Pettenati with seismological and contouring matters, F. Cavallini with the n-n approach, and M. Bobbio with the computer programs. We are also grateful to Viviana Castelli for the De Poardi reference.

References Atkinson, G. M., and D. M. Boore (1998). Evaluation of models for earthquake source spectra in Eastern North America, Bull. Seism. Soc. Am. 88, 917–934. Barnett, V., and T. Lewis (1978). Outliers in Statistical Data, Wiley Series in Probability and Mathematical Statistics; Applied Probability and Statistics, Wiley, Chichester, UK, 355 pp. Bollinger, G. A. (1977). Reinterpretation of the intensity data for the 1886 Charleston, South Carolina, earthquake, in Studies Related to the Charleston, South Carolina, Earthquake of 1886: A Preliminary Report, U.S. Geol. Surv. Profess. Pap. 1028, 17–32. Braun, J., and M. Sambridge (1995). A numerical method for solving partial differential equations on highly irregular evolving grids (abstract), Nature 376, 655–660. Cavallini, F., M. Bobbio, F. Pettenati, and L. Sirovich (2000). ConVor, a new-generation methodology for tracing objective and reproducible isoseismals: the Case of the Feb. 28, 1925 Charlevoix Earthquake in Canada. AGU Spring Meet., Washington, 30 May–3 June 2000. EOS, Trans. Am. Geophys. Union 9 May 2000, S311. Cecic⬘, I., R.M.W. Musson, and M. Stucchi (1996). Do seismologists agree upon epicentre determination from macroseismic data? A survey of ESC Working Group “Macroseismology”, Ann. Geofis. (Rome) 39, 1013–1027. Davison, C. (1921). On scales of seismicity intensity and on the construction and use of isoseismal lines, Bull. Seism. Soc. Am. 11, 94–129. De Poardi, G. V. (1627). Nuova relatione del grande e spaventoso terremoto successo nel Regno di Napoli, nella Provincia di Puglia, in Venerd’ li 30 Iuglio 1627, Roma (one map). Dewey, J. W., B. G. Reagor, L. Dengler, and K. Moley (1995). Intensity Distribution and Isoseismal Maps for the Northridge, California, Earthquake of January 17, 1994, U.S. Geol. Surv. Open-File Rept. 95-92, 35 pp. Eckstein, B. A. (1989). Evaluation of spline and weighted average interpolation algorithms, Comput. Geosci. 15, 79–94. ¨ ber das Erdbeben in den Rhein und Niederlanden Egen, P. N. C. (1828). U von 23 Feb. 1828, Ann. Phys. Chem. 13, 153–163. Hanks, T. C., and A. C. Johnston (1992). Common features of the excitation and propagation of strong ground motion for North American earthquakes, Bull. Seism. Soc. Am. 82, 1–23. Hartzell, S., S. Harmsen, A. Frankel, D. Carver, E. Cranswick, M. Meremonte, and J. Michael (1996). Site-Response Maps for the Los Angeles Region Based on Earthquake Ground Motions, U.S. Geol. Surv. Open-File Rept., 96–723 (http://greenwood.cr.usgs.gov/pub/openfile-reports/ofr-96-0723/). Mallet, R. (1862). The Great Neapolitan Earthquake of 1857, 2 vols., Chapman and Hall, London, 431, 399. Okabe, A. B., B. Boots, K. Sugihara, and S. N. Chiu (2000). Spatial Tessellation. Second ed., Wiley, Chichester, 671 pp. Pettenati, F., F. Cavallini, and L. Sirovich (1998). Quantitative treatment of macroseismic intensity and objective evaluation of synthetic results,

1940 in Proc. 11th Europ. Conf. on Earthq. Engng., Paris, Sep. 6–11 1998, Philippe Bisch, Pierre Labbe´ and Alain Pecker (Editors). Balkema, Rotterdam, CD-ROM (10 pp.; file PETQTO.PDF). Pettenati, F., L. Sirovich, and F. Cavallini (1999). Objective treatment, and synthesis of macroseismic intensity data sets using tessellation, Bull. Seism. Soc. Am. 89, 1203–1213. Pettenati, F., and L. Sirovich L. (2000). New-generation objective and reproducible isoseismals, and Tests of Source Inversion of the USGS “Felt Reports” in Proc. 4th Int. Conf. on Recent Advances in Geotechnical and Earthquake Engineering and Soil Dynamics, S. Prakash (Editor), 26–31 March 2001, San Diego California (CD-ROM), Univ. of Missouri-Rolla (ISBN-1-887009-05-1). Philip, G. M., and D. F. Watson (1986). Comments on “Comparing splines and kriging”, Comput. Geosci. 12, 243–245. Preparata, F. P., and M. I. Shamos (1985). Computational Geometry: An Introduction, Springer Verlag, New York, 390 pp. Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery (1992). Numerical Recipes in Fortran: The Art of Scientific Computing, Second Ed., Cambridge University Press, Cambridge, pp. 963. Ripley, B. D. (1981). Spatial Statistics, Wiley, New York, 252 pp. Sambridge, M., J. Braun, and H. McQueen (1995). Geophysical parametrization and interpolation of irregular data using natural neighbours, Geophys. J. Int. 122, 837–857. Sibson, R. (1980). A vector identity for Dirichlet tessellation, Math. Proc. Camb. Philos. Soc. 87, 151–155. Sirovich, L., and F. Pettenati (1999). Seismotectonic outline of SouthEastern Sicily: an evaluation of available options for the scenario earthquake fault rupture, J. Seismol. 3, 213–233.

L. Sirovich, F. Cavallini, F. Pettenati, and M. Bobbio Sirovich, L., and F. Pettenati (2002). Test of source parameters inversion of the intensities of a 54,000-death shock of the XVII Century in SE Sicily, Bull. Seism. Soc. Am. 91, 792–811. Sirovich, L., F. Pettenati, and C. Chiaruttini (2001b). Test of source-parameter inversion of intensity data. Nat. Haz. 24, no. 2, 105–131. Watson, D. F. (1992). Contouring: A Guide to the Analysis and Display of Spatial Data, Pergamon, Oxford, 321 pp. Wood, H. O., and F. Neumann (1931). Modified Mercalli intensity scale of 1931, Bull. Seism. Soc. Am. 21, 277–283. Wren, A. E. (1975). Contouring and the contour map: a new perspective, Geophys. Prospect. 23, 1–17. National Institute for Oceanography and for Experimental Geophysics (OGS) Borgo Grotta Gigante, 42c 34010 Sgonico, Trieste Italy (L.S., F.C., M.B.) National Group for Defence Against Earthquakes (GNDT) National Institute for Geophysics and Vulcanology (INGV), at OGS Borgo Grotta Gigante, 42c 34010 Sgonico, Trieste Italy (F.P.) Manuscript received 9 May 2001.