Fractal analysis with applications to seismological pattern recognition of underground nuclear explosions

Signal Processing 80 (2000) 1849}1861

Fractal analysis with applications to seismological pattern recognition of underground nuclear explosions夽 Daizhi Liu*, Ke Zhao, Hongxing Zou, Juan Su Section 103, Xi'an Research Institute of High Technology, Hongqing Town, Baqiao, Xi'an, Shanxi Province 710025, People's Republic of China Received 18 January 1996; received in revised form 15 March 2000

Abstract The problem of discriminating between underground nuclear explosions and natural earthquakes is considered in this paper. The feature of statistical self-a$ne fractal yielded from the logarithmic power spectrum is presented, and "ve novel features obtained from the evolution of the signal energy with scale are also included, together with the relationships between the features advocated. Seismological pattern recognition results are shown to demonstrate the potential of the proposed methodologies of feature extraction.  2000 Elsevier Science B.V. All rights reserved. Zusammenfassung In dieser Arbeit betrachten wir die Problematik bei der Unterscheidung zwischen nuklearen Untergrundexplosionen und natuK rlichen Erdbeben. Wir stellen hierzu die Eigenschaften statistischer selbst-a$ner Fraktale vor, die durch das logarithmische Leistungsspektrum erzeugt wurden und betrachten fuK nf neue Merkmale, die mit Hilfe der Evolution der skalierten Signalenergie gewonnen wurden, sowie die Beziehungen der Merkmale untereinander. Wir zeigen Resultate der seismologischen Strukturerkennung, um die FaK higkeiten der vorgestellten Methoden zur Feature Extraktion zu demonstrieren.  2000 Elsevier Science B.V. All rights reserved. Re2 sume2 Le proble`me de la discrimination entre les explosions nucleH aires souterraines et les tremblements de terre naturels est abordeH dans cet article. L'attribut de fractaliteH auto-a$ne statistique tireH du spectre de puissance logarithmique est preH senteH , et cinq attributs nouveaux obtenus de l'eH volution de l'eH nergie du signal selon l'eH chelle sont eH galement inclus, et les relations entre ces attributs sont deH battus. Des reH sultats de reconnaissance de formes seH ismologiques sont deH crits a"n de mettre en eH vidence le potentiel des meH thdologies proposeH es pour l'extraction d'attributs.  2000 Elsevier Science B.V. All rights reserved. Keywords: Seismic signal; Underground nuclear explosion; Natural earthquake; Pattern recognition; Fractal analysis; Wavelet analysis; Signal processing

夽 This work was supported in part by the Scienti"c Bureau of National Defense under Grant No. EP94037 and by the National Natural Science Foundation under Grant No. 49974029. * Corresponding author. Tel.: #86-29334-4090; fax: #86-29-334-4290. E-mail address: liudiaizhi}[email protected] (D. Liu).

0165-1684/00/$ - see front matter  2000 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 0 0 ) 0 0 0 9 3 - 1

1850

D. Liu et al. / Signal Processing 80 (2000) 1849}1861

1. Introduction The wave "eld f (x, y, z, t) created by natural earthquake or underground nuclear explosion is mathematically a function of time and spatial variables, where f can be the amplitude, power, or other quantities, x, y, and z denote the spatial coordinates, and t is the time. If x, y and z are set, then the wave "eld can be expressed as f (t), which is the seismographic record of a seismic station, and a kind of representation of seismic signal in a time domain. There are basically two approaches to studying the seismic signals (seismograms): one is phase analysis, such as measurement of arrival time and travel time, and inversion, etc.; the other is waveform and spectral analysis. The arrival time and travel time which can be utilized to explore the interior structure of the earth relate only to the propagation path and the characteristics of propagation media, whilst the waveform and spectral analysis involve the whole wave group, and the seismic waveform relates very much to the source characteristics, propagation media (path characteristics), the characteristics of seismograph, etc. [2]. Virtually, the phase determination is for analyzing the kinetic characteristics of seismic wave, whereas the waveform and spectral analysis are for investigating the dynamic characteristics. In seismic detection of underground nuclear explosion, the former is better suited for locating the epicenter; the latter, in many respects, is suited to identify the explosions from earthquakes. Fractal theory is without doubt an active and important branch in the area of nonlinearity, and constitutes as natural a mathematical tool for the analysis and processing of geometrical occurrences of many nonlinear problems as the Fourier transform for stationary stochastic processes. Fractal, the term coined by Mandelbrot [21] some two decades ago, has come into frequent use in a rather short time [13,22], and been brought into use in a multitude of diverse "elds of science and technology in the early 1990s [12]. Wavelet analysis, which is commonly known as `mathematical microscopea, is the latest development of harmonic analysis, and "nds its applications in many scenarios [6,7]. Wavelet-based fractal analysis has been

presented in the literature [14,30,31]. There is a substantial body of techniques for discriminating underground nuclear explosions and natural earthquakes (see, e.g., [1}6,8,10,16}19,23,24,26,27] and references therein), and consequently, di!erent techniques result in di!erent extracted features. In this paper, based on the wavelet transform and the analysis of seismic signals from the fractal perspective, three kinds of features are extracted accordingly. This paper is a sequel to [9] and presents novel results on the wavelet-based fractal features. It is illustrated that the extracted/selected features are feasible for seismological pattern recognition of underground nuclear explosions.

2. Underlying ideas and principles It is well known that a fractal is by de"nition a set with fractal dimension, and the fractal dimension D is the main feature for characterizing the fractal set, no matter whether the set is discrete, or continuous, whether in the sense of pure mathematics, or of statistics. The essence of fractal analysis is to calculate the corresponding fractal dimension D of di!erent fractal sets using di!erent dimension estimation methods [15,20,25,32]. The seismic signals in time domain can, in general, be regarded as a set and the self-a$ne signals, just like most geophysical signals [29]. The basic concept behind the fractal approach to seismological pattern recognition of underground nuclear explosions is to extract some features via the analysis and processing of seismic waveform, then construct some classi"cation discriminants, and ultimately distinguish the explosion events from earthquake events automatically. Two major problems in the study of fractality are: (1) whether the object of interest possesses the fractal properties; (2) under a certain resolution, how much the fractal dimension D is. In fractal analysis, the key step is to "nd the scale-invariant interval so as to evaluate the corresponding fractal dimension. In order to determine the scale-invariant interval, i.e., to evaluate the fractal dimension D, the scale decomposition should be carried out. The ideas underlying the fractal-based approach to analyzing the seismic signal are as follows:

D. Liu et al. / Signal Processing 80 (2000) 1849}1861

(1) Fourier analyze the seismic signal within a "nite time duration, and determine whether the seismic signal in time domain has the properties of selfa$ne fractal in accordance with the relationship between power spectrum and frequency (whether the relationship "ts a power law or not), and then, calculate the fractal dimension D from the slope of a straight line which is "tted from a log}log plot of power spectrum and frequency; (2) perform the wavelet decomposition for seismic signal in time domain, "nd the signal `energya (in a broad sense) at di!erent scales, investigate the relationship between the signal energy at each scale, then, evaluate the relevant features; (3) use these features as the classi"cation discriminants, perform the recognition tests, and "nally, select the adequate features for seismological pattern recognition applications. 2.1. Self-azne fractal dimension D To "x our notation, let x(t)3¸(R) denote the seismic signal in time domain (vertical component as involved in our work). For 0)t)¹, the Fourier transform of x(t) is de"ned by



X( f, ¹)"

2

x(t)e\p DR dt

(1)



and the power spectral density function is given by 1 S( f )" "X(f, ¹)". ¹

(2)

If x(t) has fractal properties, the relationship between the power spectral density and the frequency will obey the power law in a certain frequency band such that S( f )"Kf\@.

(3)

Indeed, via Mandelbrot's work we can express the fractal dimension D by [29,31] D"(5!b)/2.

(4)

Take the logarithm of both sides of Eq. (3), "t the log}log plot with a straight line (the slope is just !b), and the fractal dimension D will be easily computed via Eq. (4).

1851

2.2. Multiresolution energy fractal features For the signal whose power spectrum possesses the characteristics of self-a$ne fractal, our research indicated that its log}log plot of the energy of continuous wavelet transform with respect to scales has a scale-invariant interval, from which a new feature as what we call the energy fractal dimension can be derived. Since in such a continuous case the scales can take continuous values, the scale-invariant interval can therefore be determined precisely, but at the expense of greatly increased computations. Unlike its continuous counterpart, the discrete wavelet transform is much more computationally e$cient, and is therefore well suited for extracting features from the large dataset. In regard to this consideration, the present paper is just based on the discrete wavelet transform. As we know, the multiresolution wavelet transform can be considered as decomposition of a signal of interest into its coarse and detail components, respectively. Approximation of the signal at resolutions 2\H, j"1,2,2, can be seen as viewing the signal at decreasing levels of detail, and this is analogous to the upward continuation of geophysical potential "eld in which the continuation height corresponds to the wavelet decomposition scale. Both cases have the e!ect of low-pass "ltering. Take dyadic wavelet decomposition as an example. The scale factor is a"2H, where j represents the scaling index. When j"1, the scale is small, the time window is narrow, but the time resolution is high, this case corresponds to the lower continuation, and consequently, the "ner details of the original signal can be re#ected. With the increase of j, the time window becomes wider, the frequency resolution becomes higher, this case is similar to the higher continuation, and as a result, the lower frequency component of the original signal can be revealed. When a signal is scaled from a resolution of 2\H\}2\H (corresponding to upward continuous one), some signal details are discarded, and this discarded information is retained as the detail component, namely, the detail signal, at resolution 2\H. Needless to say, the coarse approximation signals and detail signals at di!erent scales are closely related to the scale. In view of this consideration,

1852

D. Liu et al. / Signal Processing 80 (2000) 1849}1861

we regard the energy of the coarse approximation signals and detail signals at di!erent scales as functions of decomposition scale, therefore, relevant features can be obtained from this kind of function. The concrete computation procedure is as follows. Again, let x(t)3¸(R) denote the original seismic signal in time domain. Certainly, the real observations of x(t) are discretized. Let us denote by c "+c ,, n"1,2,2, N, the samples of x(t). Ac L cording to multiresolution analysis, the coarse approximation signal and detail signal at di!erent scale ( j"1,2,2, J) can be obtained recursively from +c ,. At each scale 2H, let us denote by L c "+c , Z , the coarse approximation signal, by H HI IZ d "+d , Z the detail signal, such that H HI IZ , c " h(n!2k)c , (5) HI H\L L , d " g(n!2k)c . (6) HI H\L L Two pairs of mirror "lters used in this paper are plotted in Fig. 1, where Fig. 1(a) corresponds to biorthogonal wavelets, Fig. 1(b) to the orthogonal wavelets. Let E ( j) be the energy function of coarse  approximation signal, and E ( j) of detail signal,  respectively, viz., E ( j)" c ,  HI IZ8

(7)

E ( j)" d , (8)  HI IZ8 then we have the energy E and E at a certain scale   j. Within a certain scale interval (scale-invariant interval), if E ( j) and E ( j) satisfy the power law   with respect to j, i.e., E ( j)&a" , E ( j)&a" , (9)  H  H then D and D can be calculated from the slopes of   straight lines "tted from log}log plots of E ( j) and  E ( j). By analogy with the fractal dimension D, we  could refer to D and D as the energy fractal   dimensions. Since the feature D characterizes ap proximately the scale-invariant interval in the low resolutions, for di!erent samples of the same nature, the scale-invariant intervals of D may be  di!erent, and occasionally, this kind of di!erence

Fig. 1. Waveforms of "lter pair h and g.

may be signi"cant. In [9], we have investigated the discriminative power of D and D , but we are now   concerned with the following problem: in which common scale-invariant interval, the coarse approximation signals and detail signals of all training samples carry the most discriminative information. This deserves attention as it concerns new methods of feature extraction. In this paper, since the number of sample points N for each truncated seismic signature is 1024, the maximum decomposition scale J of discrete dyadic wavelet transform is therefore set to 10 (taking into account the fact that log N"10). It is manifest  that D and D are closely related to the scale   intervals within which the straight lines are "tted. If we choose arbitrarily two scales to form a scale interval, we will have C groups of scale intervals.  Within each scale interval, calculating D and  D for each signal yields two sets of features, then  the two scale intervals corresponding to the two best recognition results achieved by a dichotomy decision rule [28] will be used as the `besta scale intervals for D and D , respectively. Although the   `besta scale interval is essentially di!erent from the scale-invariant interval, we still treat D and D as   the energy fractal feature. We attribute this to the following two reasons: (1) feature D re#ects the  variation trend of the energy of coarse approximation signals versus the scales; therefore, D  can characterize the trend of log}log plot in the

D. Liu et al. / Signal Processing 80 (2000) 1849}1861

scale-invariant interval (i.e., the longest linear part in log}log plot); (2) feature D re#ects the unsym metry with respect to the apex in the log}log plot of detail signals; hence, it can also characterize the scale-invariant interval to some extent. In the above procedure, we observed a some kind of apex appearing on the curve of E ( j), both sides  of the apex can be "tted with a straight line; hence another two slopes (denoted by D  and D  ) are   obtained. In fact, the features D  and D  are in   essence of generalized fractal dimensions. Meanwhile, we noticed that, the scale position aH corresponding to the apex can also be used as a feature for pattern recognition applications. (Note that aH is found from the quadratic curve "tted from the same plot of E ( j), hence it may not be an integer.)  In summary, after processing and analysis of the seismic signal x(t), we have extracted the following three kinds of features: (1) self-a$ne fractal dimension D; (2) energy fractal dimensions D , D , D     and D  ; (3) scale position aH of the apex of detail  signal energy. 2.3. Adopted pattern classixers Two approaches to pattern classi"cation, namely, the dichotomy decision rule for single features and the nearest-neighbor rule [28] for multifeatures, are used in the present work. Central to the dichotomy decision rule is to "nd a threshold that can almost correctly delineate the point of demarcation between `largea and `smalla features, while the idea behind the nearest-neighbor rule is to "nd two `besta prototypes which assigns a pattern of unknown classi"cation to the class of its nearest neighbor. In this paper, we use the criterion whether the average correct recognition rate is the highest to determine the `besta threshold and the `besta prototypes. The rationale of this criterion is to treat the false alarm and false dismissal as equally important factors, or speci"cally, to assign them

 The word `besta appears in quotes when it is used in this less restrictive sense. Since, strictly speaking, both the dichotomy decision rule and the nearest-neighbor rule are suboptimal procedures; their use will usually lead to error rates greater than the minimum possible, the Bayes rates.

1853

with the same risk probabilities; otherwise, if we emphasize the high recognition rate for the nuclear explosion alone, the probability of false alarm will increase accordingly. Since, as is often the case, large dataset of underground nuclear explosions is not yet available for evaluating and benchmarking the existing algorithms and pattern recognition systems we use all the samples on hand as training samples, so the correct recognition rate we obtained is actually the one for the training samples, which, as we may anticipate, is higher than that obtained from the Bayesian classi"er.

3. Computation results and discussions 3.1. Feature extraction We extracted the self-a$ne fractal dimensions and the multiresolution energy fractal dimensions from the samples of short-period dataset (consisting of 80 vertical components of explosion records and 80 vertical components of earthquake records, with the sampling frequency f  "40 Hz), and the sam ples of intermediate-period dataset (consisting of 50 vertical components of explosion records and 50 vertical components of earthquake records, with the sampling frequency f  "20 Hz), respectively.  One thousand and twenty four sample points starting with the onset were truncated from each event. Fig. 2(a) depicts an earthquake waveform (intermediate period), Fig. 2(b) depicts the logarithmic power spectrum corresponding to Fig. 1(a) with b+!0.3558 and D+2.6779. Fig. 2(c) illustrates an explosion waveform (short-period), and its corresponding logarithmic power spectrum is shown in Fig. 2(d) with the features b+!4.8643 and D+4.9321, respectively. Throughout the computation we found that the scale-invariant intervals of all the samples are located in the low-frequency parts of the logarithmic power spectra. To be more speci"c, the scale-invariant intervals for intermediate-period samples are located in the beginning of the spectra, akin to that of Fig. 2(b). For the shortperiod samples, barring a few exceptions, most of the scale-invariant intervals are located in the slopes ahead of the apices, which are similar to that of Fig. 2(d). In general, the scale-invariant intervals

1854

D. Liu et al. / Signal Processing 80 (2000) 1849}1861

Fig. 2. (a) Waveform of a natural earthquake (intermediate-period) and (b) its corresponding logarithmic power spectrum; (c) waveform of an underground nuclear explosion (short-period) and (d) its corresponding logarithmic power spectrum.

Fig. 3. Histograms of feature D for (a) intermediate-period earthquake samples, (b) intermediate-period explosion samples, (c) short-period earthquake samples, and (d) short-period explosion samples. (Note: The number of bins is 15.)

for short-period samples di!er signi"cantly as compared with those for intermediate-period samples, therefore, they are di$cult to be determined automatically. Figs. 3(a)}(d) show the histograms of feature D for intermediate period earthquake and

Fig. 4. Energy fractal features for a certain short-period explosion sample.

explosion samples, and for short-period earthquake and explosion samples, respectively. Fig. 4 presents typical energy fractal features of a short-period explosion sample. Fig. 4(a) indicates the relationship between the logarithmic energy log E ( j) of coarse approximation signal and the  logarithmic scale log 2H. Fitting the curve within the interval j"1}7 yields D +!1.9200. Fig. 4(b) 

D. Liu et al. / Signal Processing 80 (2000) 1849}1861

1855

illustrates the relationship between the logarithmic energy log E ( j) of detail signal and the logarithmic  scale log 2H. Fitting this curve within the interval j"4}8 yields D +!2.6772. It is easily found  from the curve in Fig. 4(b) that an apex corresponds to the scale s"j"5, therefore, curve "tting within the intervals j"1}s and s}8, respectively, will yield two features D  +2.4133 and D  +!3.6481.   Moreover, a quadratic curve "tting within j"3}10 yields aH+2.6896. It is worthwhile noting that all the scale-invariant intervals involved in the aforementioned computations are best intervals which correspond to the highest correct recognition rates. Note also that for di!erent scale intervals, the correct recognition rates may be the same, largely due to the small size of training sets. Di!erent samples have di!erent s values as evidenced by Fig. 5. Conceptually, the curve, as illustrated in Fig. 4, can be interpreted as the `energy spectruma. It is clear that the energy of coarse approximation signal, E ( j), decreases with the increase of scale 2H. As  we know, when the scale increases, the components decomposed (forming the detail signal) increase, and consequently, the energy of the remainder of the signal (forming the coarse approximation signal) decreases. One of the main reasons for an apex

to occur on the `energy spectruma of the detail signal is that, the energy of seismic signal at di!erent frequencies are di!erent. Since the amplitude spectra of seismic signals created by natural earthquakes and underground nuclear explosions have the so-called `peak frequenciesa (refer to Fig. 2), two features (i.e., the peak frequency and steepness) extracted from the amplitude spectrum can be used to identify the source nature [2]. Recall that the dyadic wavelet function is by nature a constant Q "lter, its center frequency and bandwidth will decrease by a factor of 2 with the increase of j, in this sense, we say that the peak frequency F of power spectrum must appear in the passband of the wavelet function at a certain scale. Therefore, it is straightforward to show that the peak frequency F corresponds to a scale aH at which the highest energy of detail signal appears. Fig. 6 shows the histogram of feature F. Since the underlying physical signi"cance of F and aH are the same, and evaluating the steepness of power spectrum in general is not easy, we shall ignore the notational distinction, treating D  and D  as the steepness of   power spectrum on both sides of peak frequency F. Pattern classi"cation experiments have shown the feasibility of the features aH, D  and D  . Figs. 7}9   show respectively the histograms of D , D and  

Fig. 5. Histograms of s values for (a) intermediate-period earthquake samples, (b) intermediate-period explosion samples, (c) short-period earthquake samples, and (d) short-period explosion samples. (Note: The number of bins is 7.)

Fig. 6. Histograms of peak frequency F for (a) intermediateperiod earthquake samples, (b) intermediate-period explosion samples, (c) short-period earthquake samples, and (d) shortperiod explosion samples. (Note: The number of bins is 20.)

1856

D. Liu et al. / Signal Processing 80 (2000) 1849}1861

Fig. 7. Histograms of feature D for (a) intermediate-period  earthquake samples, (b) intermediate-period explosion samples, (c) short-period earthquake samples, and (d) short-period explosion samples. (Note: The number of bins is 15.)

Fig. 9. Histograms of feature aH for (a) intermediate-period earthquake samples, (b) intermediate-period explosion samples, (c) short-period earthquake samples, and (d) short-period explosion samples. (Note: The number of bins is 20.)

aH for the training samples. In fact, the energy fractal features D , D , D  and D  characterize     approximately the signal in di!erent frequency bands, while the self-a$ne fractal feature D characterizes comprehensively the signal in the scale-invariant interval. From the perspective of signal energy, however, the energy fractal features are analogous to the self-a$ne fractal feature to some extent, and the wavelet decomposition scale is inherently related to frequency. Since they have substantial similarities, we shall elaborate this subject after introducing the recognition results. 3.2. Classixcation using single features

Fig. 8. Histograms of feature D for (a) intermediate-period  earthquake samples, (b) intermediate-period explosion samples, (c) short-period earthquake samples, and (d) short-period explosion samples. (Note: The number of bins is 20.)

How discriminative are the single features? A pragmatic answer is to assign the discrimination power based on the recognition performance using each single feature. Table 1 presents the correct recognition rates using the dichotomy decision rule from the previously discussed features D, D , D , D  , D  , aH and F. In addition, the best     scale intervals for linear and quadratic curve "tting

D. Liu et al. / Signal Processing 80 (2000) 1849}1861 Table 1 Recognition results using single features Feature

Short period

Table 2 Recognition results using multifeatures Intermediate period Features

D D  D  F aH D  D 

68.75% j"1}7, 83.75% j"4}8, 84.375% 70% j"3}10, 79.375% j"1}s, 56.25% j"s}8, 66.875%

1857

78% j"1}7, 82% j"3}5, 82% 65% j"2}8, 83% j"1}s, 63% j"s}8, 69%

are also given in Table 1. By reference to Table 1 we see that the features D, D , D and aH can provide   more informative discriminations, and thus are preferred to the features D  and D  . The discrimina  tion power of feature D for intermediate-period samples is higher than that for short-period samples, while the discrimination power of feature F is inferior to that of feature aH for both cases. The recognition results may ultimately prove features D and D to be equally valuable, and to be su  perior to feature D.

D, D , D , aH   D, D , D   D, D , aH  D, D , aH  D , D , aH   D, D  D, D  D, aH D ,D   D , aH  D , aH 

Short period (%)

Intermediate period (%)

85.625 85 83.75 84.375 86.25 83.75 85 83.125 84.375 84.375 86.25

86 86 87 87 84 88 84 87 82 83 83

3.3. Classixcation using multifeatures In this part, we use the nearest-neighbor rule to investigate the discrimination power of di!erent combinations of D, D , D and aH. The highest   correct recognition rate for short-period samples can be achieved by the combination of D and aH,  and that for intermediate-period samples can be achieved by the combination of D and D , as evid enced by Table 2. Since most of the power spectra of short-period samples di!er signi"cantly, and are more complicated than that shown in Fig. 2(d), the feature D may adversely a!ect the performance of the combined features. In order to visualize the discrimination power of the combined features, we devised a novel surrogate feature}feature distribution, as illustrated in Figs. 10}15 (where, `
a stands for explosion samples, `#a for earthquake samples). In the upperright corners of Figs. 10}12, an enlarged portion of the distribution is also illustrated. These distributions are obtained from the intermediate-period samples, the ones obtained from the short-period

Fig. 10. D versus aH. 

samples are not included in this paper bec-ause of space limitation. 3.4. Further discussions of relationships between features The aforementioned features D, D , D , D  , D      and aH reveal physically how the signal energy evolves with frequency or with scale, thus they warrant special mention because subtle ties may exist among them. Fortunately, this kind of tie can be exploited by means of "gures rather than by rigorous mathematical derivations.

1858

D. Liu et al. / Signal Processing 80 (2000) 1849}1861

Fig. 11. D versus aH. 

Fig. 13. D versus D .  

Fig. 12. D versus aH.

Fig. 14. D versus D . 

As was mentioned earlier, the peak frequency of power spectrum corresponds to the scale at which the highest energy of detail signal appears. The dotted line in Fig. 2(d) is "tted within the predetermined scale-invariant interval, its slope b is related to the fractal dimension D via Eq. (4). It is convenient to regard approximately b as the steepness of the spectrum on the low-frequency side of peak frequency. Since the coarse approximation signals at lower scales and detail signals at higher scales can better characterize a signal, the feature D  should be evaluated from the energy at lower

scales, and correspondingly, the feature D should  be evaluated from the energy at higher scales. As expected, the searched best scale intervals are consistent with this qualitative analysis. From Fig. 4(b) we can see that the feature D  corresponds to the  steepness of the power spectrum on the high-frequency side of the peak frequency, evidently it is not a good choice, as against the feature D  .  The above explanations in general hold true for the short-period samples, but not valid for the intermediate-period samples. It was not clear until recently what caused this di!erence, but it could be

D. Liu et al. / Signal Processing 80 (2000) 1849}1861

1859

the unique value. Because aH is evaluated by quadratic curve "tting on the energy of multiresolution signals, aH reveals the information of scales, on one of which the highest energy of detail signal appears. Furthermore, these scales involve the information of steepness of power spectrum in the neighborhood of peak frequency. In one word, aH consists of two kinds of information, one is the peak frequency, the other is the steepness, hence aH is more discriminative than the peak frequency feature F. 3.5. Ewect of diwerent wavelets

Fig. 15. D versus D . 

that for the latter ones, important information is contained in the lower parts of the power spectra, then the corresponding scale-invariant intervals are located automatically at the beginning of the frequency band, and as a result, the yielded feature D is more discriminative than D  and D  .   The features D and D are yielded within the   searched best scale intervals that guarantee the best recognition results for both explosion and earthquake samples. However, looking more closely at the underlying physics of D and D , which was our   main motivation that led to our discovery of the other fractal features, we observed that "nding the best scale intervals is equivalent to "nding the best frequency bands of the power spectra in the vicinity of peak frequencies; therefore, their discrimination power are almost the same, as evidenced by Table 1. Surprisingly, the discrimination power of scale feature aH is relatively high. It is worth noting, however, that the correspondence between peak frequency of power spectrum and scale to the highest energy of detail signal depends on the choice of wavelets. If a certain wavelet other than the dyadic one is used, the resulting aH will be certainly di!erent. Since the resolution of dyadic wavelet per se is "nite, the scale feature aH is an approximate value. But if the continuous wavelet is adopted, then we can obtain a unique value for the scale feature aH. Therefore, aH can be viewed as an approximation to

Generally speaking, the choice of wavelet is application dependent, where, typically, but not necessarily, we use the biorthogonal wavelet (recall Fig. 1(a)). In this paper, we shall not delve into details of the choice of wavelet, but we have to bear in mind that if we choose the wavelet judiciously, fascinating new insights into the signal will be highlighted. Table 3 presents the recognition results using Daubechies wavelet (refer to Fig. 1(b)). From Tables 2 and 3 we can see that there are no dramatic di!erences between the two results. The only signi"cant di!erence may be due to the feature D, where, notwithstanding the discrimination, power of D (refer to Table 1) is inferior to that of D ,  D and aH for the intermediate period samples, the  feature D does provide informative discrimination when combined with other features. Since the Table 3 Recognition results via Daubechies wavelets

Features D  D  aH D, D , D , aH   D, D , D   D, D , aH  D, D , aH  D , D , aH   D, D  D, D  D, aH D ,D   D , aH  D , aH 

Short period (%)

Intermediate period (%)

83.125 85.625 82.5 86.88 86.25 85 86.88 86.88 84.375 85.625 83.125 85.625 85 86.88

84 83 82 87 87 86 86 84 85 86 86 84 83 84

1860

D. Liu et al. / Signal Processing 80 (2000) 1849}1861

scale-invariant interval for evaluating D can be automatically searched for the intermediate-period samples, in this case, D is recognized as a valuable feature.

4. Concluding remarks The increasing demand for improving the accuracy of discrimination between underground nuclear explosions and natural earthquakes has introduced a problem of limited reliability and e$ciency in the Comprehensive (Nuclear) Test Ban Treaty (CTBT) verixcation research. Although there are many discrimination techniques currently available, there still exists a need to develop more e$cient and convincing methods to further alleviate this problem. In this paper, several methods of feature extraction from a fractal perspective for seismological pattern recognition have been proposed. The recognition results provided suggest that the features D, D , D and aH are good reme  dies for seismological pattern recognition pathology. It should be remarked that the goal here is not to demonstrate the power of the features advocated or compare their performances with conventional techniques for discrimination, but rather to use these features as a vehicle for characterizing the nature of source when applied to the same task. One of the major issues in the "eld of pattern recognition has always been which feature, if any, is absolutely `besta. Given these developments, together with many other successful features documented in the literature, it appears that at this stage of our knowledge, trying to prove which feature is `besta is premature, to say the least. Because the wavelet-based fractal analysis approach to the topic is relatively new, there remains, quite naturally, many unanswered questions. Nevertheless, how to determine the scale-invariant interval in which the fractal phenomenon holds is a question of taste.

5. Further reading The following reference is also of interest to the reader: [11].

Acknowledgements The authors would like to express their gratitude to Prof. Yinkang Wei and Dr. Guofu Sun for their valuable input to this work. Thanks also are due to the anonymous associate editor and the reviewers for their constructive remarks and helpful suggestions in the preparation of this paper.

References [1] V.V. Adushkin, I.O. Kitov, Distinguishing between underground and surface (contact) explosions, Trans. (Dokl.) USSR Acad. Sci. Earth Sci. Sect. 325 (6) (1994) 80}83. [2] M. Bath, Spectral Analysis in Geophysics, Elsevier, Amsterdam, 1977. [3] B.A. Bolt, Nuclear Explosions and Earthquakes, The Parted Veil, Freeman, San Francisco, 1976. [4] C.H. Chen, Feature extraction and computational complexity in seismological pattern recognition, Proceedings of the Second International Joint Conference on Pattern Recognition, IEEE, Piscataway, 1974. [5] C.H. Chen, On statistical and structural feature extraction, Pattern Recognition and Arti"cial Intelligence, Academic Press, New York, 1976. [6] C.H. Chen, I.C. Lin, Pattern analysis and classi"cation with the new ACDA seismic signature data base, Report AD. A015, 925, Arms Control and Disarmament Association, Washington, DC, 1975. [7] C.K. Chui, Wavelets: A Tutorial in Theory and Applications, Academic Press, New York, 1992. [8] O. Dahlman, H. Israelson, Monitoring Underground Nuclear Explosions, Elsevier, Amsterdam, 1977. [9] Daizhi Liu, Hongxing Zou, Yinkang Wei, Ke Zhao, Juan Su, Weimin Jia, Fractal analysis with applications to seismic pattern recognition of nuclear explosion, in: Baozong Yuan, Xiaofang Tang (Eds.), 1996 Third International Conference on Signal Processing Proceedings of the ICSP'96, Publishing House of Electronics Industry, Beijing, China, 1996. [10] G.R. Dargahi-Noubary, Identi"cation of seismic events based on stochastic properties of the short-period records, Soil Dynamics Earthquake Eng. 17 (2) (1998) 101}115. [11] I. Daubechies, Ten Lectures on Wavelets, SIAM, Philadelphia, PA, 1992. [12] K.J. Falconer, Fractal Geometry: Mathematical Foundations and Applications, Wiley, New York, 1990. [13] J. Feder, Fractals, Plenum Press, New York, 1988. [14] P. Flandrin, Wavelet analysis and synthesis of fractional Brownian motion, IEEE Trans. Inform. Theory 38 (2) (1992) 910}917. [15] R. Fox, M.S. Taqqu, Large-sample properties of parameter estimates for strongly dependent stationary Gaussian time series, Ann. Stist 14 (2) (1986) 517}532.

D. Liu et al. / Signal Processing 80 (2000) 1849}1861 [16] R.B. Ives, Dynamic spectral ratios as features in seismological pattern recognition, Conference on Computer Graphics, Pattern Recognition and Data Structure, IEEE, Piscataway, 1975. [17] G.S. Jang, F. Dowla, V. Vemuri, Application of neural networks for seismic phase identi"cation, IEEE International Joint Conference on Neural Networks, IJCNN'91, 1991, pp. 899}904. [18] G.S. Jang, F. Dowla, V. Vemuri, Performance comparison of some neural network paradigms for solving the seismic phase identi"cation problem, Proc. SPIE 1721 (1992) 265}276. [19] H.H. Liu, A rule-based system for automatic seismic discrimination, Pattern Recognition 18 (6) (1984) 459}463. [20] T. Lundahl, W.J. Ohley, S.M. Kay, R. Si!ert, Fractional Brownian motion: a maximum likelihood estimator and its application to image texure, IEEE Trans. Med. Imaging MI-5 (9) (1986) 152}161. [21] B.B. Mandelbrot, Fractals: Form, Chance and Dimension, Freeman, San Francisco, CA, 1977. [22] B.B. Mandelbrot, The Fractal Geometry of Nature, Freeman, San Francisco, CA, 1982. [23] C.L. Mason, Cooperative seismic data interpretation for nuclear test ban treaty veri"cation, Appl. Artif. Intell. 9 (4) (1995) 371}400. [24] M. Myers, Computer-aided seismic analysis and discrimination, Computer 10 (10) (1977) 64}71.

1861

[25] H.E. Schepers, J.H.G.M. van Beek, J.B. Bassingthwaighte, Four methods to estimate the fractal dimension from self-a$ne signals, IEEE Eng. Med. Biol. Mag. 11 (2) (1992) 57}64. [26] S.R. Taylor, H.E. Hartse, Evaluation of generalized likelihood ratio outlier detection to identi"cation of seismic events in Western China, Bull. Seismolog. Soc. Am. 87 (4) (1997) 824}831. [27] T. Tiira, Discrimination of nuclear explosions and earthquakes from teleseismic distances with a local network of short period seismic stations using arti"cial neural networks, Phys. Earth Planet. Inter. 97 (1}4) (1996) 247. [28] J.T. Tou, R.C. Gonzalez, Pattern Recognition Principles, Addison-Wesley, Amsterdam, 1974. [29] D.L. Turcotte (Chinese translators: Chen Yong et al.), Fractal and Chaos with Applications to Geology and Geophysics, Beijing, Seismic Press, 1993, pp. 79}96 (Chapter 7). [30] G.W. Wornell, A Karhunen}Loe`ve-like expansion for 1/f processes via wavelets, IEEE Trans. Inform. Theory 36 (7) (1990) 859}861. [31] G.W. Wornell, A.V. Oppenheim, Estimation of fractal signals from noisy measurements using wavelets, IEEE Trans. Signal Process. 40 (3) (1992) 611}623. [32] Y. Yajima, On estimation of long-memory time series models, Austral. J. Stist. 27 (3) (1985) 321}325.