GEOPHYSICAL RESEARCH LETTERS, VOL. 27, NO. 20, PAGES 3433-3436, OCTOBER 15, 2000 Capabilities of 3-D Wavelet Transforms to Detect Plume-like Structures from Seismic Tomography. StephenY. Bergeronand David A. Yuen, Dept. of Geology and Geophysics,University of Minnesota, Minneapolis, MN $$415-1227, USA and Minnesota Supercomputer Institute, University of Minnesota, Minneapolis, MN $$415-1227, USA. Alain P. Vincent D•partement de Physique, Universit• de Montreal, C.P. 6128, Succ. Centre-Ville, Montreal, QC H3C 3J7, Canada and CERCA, 5160 Boul. D•carie, Bureau 400, Montreal, QC H3X 2H9, Canada. Abstract. The wavelet transform methods have been applied to viewing 3-D seismic tomography by casting the transformed quantities into two proxy distributions, E-max, the maximum of the magnitude of the local spectra about a local point and the associatedlocal wavenumber, k-max. Using a stochasticbackground noise, we test the capability of this procedure in picking up the coherent structures of upper-mantle plumes. Plumes with a Gaussian shape and a characteristic width up to 2250 km have been tested for var- iousamountsof the signal-to-noise ratios (SNR). We have found that plumes can be picked out for SNR as low as 0.08 db and that the optimal plume width for detection is around 1500 km. For plume width ranging between 700 km and 2000 kin, the SNR can be lower than 1 db. This lengthscale falls within the range for plume-detection based on the signal-to-noiselevels associatedwith the current global tomographical models. intrinsically noisy and have errors due mainly to the resolution, the choice of inversion algorithm and the parametriza- tion used[Boschiand Dziewonski,1999].We havenot been investigated the effect of noise in the dataset and how this may influence the ability of wavelets for detecting coherent structures in the mantle. The coherent structures we have in mind are mantle plumes, which supposedlyare portrayed by the sharp lateral temperature gradients and the associated elastic and anelastic properties. Recently, there have been some successesin detecting plumes in the mantle with to- mography[Bijwaardand $pakman,1999;$hen et at., 1998; Ni et at., 1999; Atten et at., 1999]. Thus, it is important to separate out the lower- and upper-mantle components of these plume-like objects. For this reason, regional tomography should be used in conjunction with global tomographic models. Wavelets can be used to cross-correlate these differ- ent scalesin the different tomographic models. The purpose of this paper then is to illustrate the capabilities of wavelets Introduction in discriminating plumes from the presenceof various levels of noise in the background mantle. This assessmentis conIn the last severalyears, waveletshave begun to take hold ducted by means of numerical experiments on the synthetic in the scientificand engineering community[Burke-Hubbard, data. One of the critical points we wish to communicate is 1998; Wickerhauser, 1994]as a truly viabletool for investi- that there is an optimal way for extracting the most valuable gating a wide variety of phenomenawith multiple-scalesand information to find plumes. to extract features from complex data-sets. Wavelets have enabled fluid dynamicists to delineate coherent structures Description of Wavelet Technique and in turbulent flows [Fargeet at., 1996], to analyzespatialthe Extraction Scheme temporalpattern in seismicmeasurements[Bethouxet at., 1998]and to determinedifferentregimesin atmospheric and The wavelet technique is a localized transform in both climate dynamics [Kumar and Foufouta-Georgiou, 1993]. physical and wavenumber space that has been formalized The idea of local 1-D spectral wavelet analysis has also been in the 1980'sby [Grossmannand Morte•, 1984;Daubechies, used to trace the sourcesof postglacial rebound in a spher- 1988; Matta•, 1989] among others. The dual approachin ical geometry[Simonsand Hager,1997]. physical and spectral space makes it possible to extract sigRecently[Bergeronet at., 1999,2000]haveemployed3-D nificant but weak, higher-order features that would be dif- wavelet transforms to study seismicanomalies derived from ficult to discern with a conventional Fourier or windowed tomographicinversion[Zhou,1996]. Thesestudiescalled Fourier approach[Daubechies, 1992]. With the classical for the usage of two proxy quantities derived from applying the wavelet transform to the seismic data set. Some inter- esting results, such as the detection of the Icelandic plume [Bijwaardand Spakman,1999], came out from this initial exploration with wavelets. However, the seismicmodels are Copyright2000 by the AmericanGeophysicalUnion. Papernumber 1999GL011243. 0094-8276/00/1999GL011243505.00 Fourier transform the scale is global and imposed a priori, whereas with the wavelets one can chooselocally the scale of interest. This formalism leads naturally to the concept of local spectra[Perr/er e• at., 1995]. Here we review briefly the wavelet local spectra extraction technique in 3-D. We also succinctly describethe extraction procedure of the two proxy quantities, E-max and k-max, in order to synthesize this large amount of information. In contrast to our pre- viouspresentation[Bergetone• at., 2000],all mathematical expressionsare given in the physical space. We employ the 3433 3434 BERGERON ET AL.- CAPABILITIES OF 3-D WAVELETS TO DETECT PLUMES whether the local anomaly is slow (negative)or fast (positive) relative to a 1-D backgroundEarth model suchas the PREM [Dziewonskiand Anderson,1981]. Obviously,this information is lost by using the L2-norm. Consequently, we chooseto incorporate the sign of the signal in the representation of E-max. [Bergeronet al., 2000] have employedthis waveletapproach to detect a simplified subducting slab and a simplified mantle plume in the presenceof a noisy background. In the case of the slab, the boundaries are sharp, in the sense D that the scalar field associated with the structure varies abruptly from 0 to 1. In the caseof the plume-like structure, < ........................................ D ............................................................ the boundaries are smoother since they are varied accordingly to a Gaussian function. They have demonstrated that the method can be used with a high noise level of ldb and Figure 1. Dimensionusedfor the synthetic mantle plume. The that sharp boundaries are easier to detect by a simple visual box width D and height H are fixed to 10,000 km and 700 km respectively. The plume is a Gaussian function centered in the inspection. Here we are concernedwith the sensitivity of our horizontal plane with a diameter of 2rrn varying from 250 km to method for detecting plume-like structures. The main pur2250 km. pose of this paper is to determine at which noise level and the width of the structure this proxy approach can be used efficiently. secondderivative of the Gaussian, also known as the MexiThe plume-like structure is independent of the vertical can hat as describedby [Murenziand Antoine,1996]for the direction and decays as a Gaussian function along the median vertical axis with an horizontal standard deviation an: wavelet •b' •p•'-g= 3-•'-•• exp -• a a a , Splume( • = exp_x 2rr• _!/2 )' (1) (3) whereF- (x, y, z) is the vectorpositionand b is the location The width of the structure is given by two times the stanwhere the center of the wavelet is positioned. The scalinga is isotropic and the wavelet is dilated equally in each direction. We also imposed that the L2-norm, a mathematical measure of the magnitude, of the wavelet is scale independent. The wavelettransformof the signalf (here the seismicvelocity anomalies)is given by f(r-') •, a (2) where Lx, Ly and Lz are the lengths of the periodic box in Cartesianspace. The local wavenumberk• is related to the inverseof the scale a by k' - 1/a. Waveletscan thus be viewed as an operation involving a 3-D spatial convolution integral, which is set to zoom in both physical space at a given location and in spectral space at a given scale. The best way for handling this transformation is to go to the Fourier space, since the convolution involved in the previous equation becomes then a multiplicative operation. The dard deviation. The typical size of interest here is a plume with a height H of 700 km and a width D varying from 225 km to 2250 km. Fig. I depicts sucha typical plume. We note thatother forms forSplum e(r-)arealso likelyandshould be examined in the future. The noise is a Gaussian stochastic background generated by usingthe techniquedescribedin [ Vincentand Meneguzzi, 1991]. It followsa power-lawdecayin the Fourier space. Here we have fixed the value of the power-law exponent to -2, since recent estimates of 3-D tomographic data show that the value shouldlie between-2 and-2.5 [Cadek et al., 1998;Passlet and Snieder,1995]. This stochasticbackground,noted Bnoise,can be viewedhas a noisyfield in which the plume-like signal is drown out. The purpose of this construction is to pick out the plume-like structures Splum e fromthenoisy dataSmeasured , where $measured (0 = $plume (0 q-C Bnois e(0- (4) detailednumericalprocedurecan be found in [Bergeronet al., 2000]. The local power spectrum of the signal at a given and fixed location b is given by the L2-norm of the wavelet trans- format thislocation: II](a,g)11 Thetypeof information produced by this analysis is difficult to visualize, since at each location we need to render a whole spectrum depend- ing on the waveletwavenumberk• = 1/a. In order to synthesize and digest this large amount of information more efficiently, we extract two proxy quantities: the maximum of the local energy, E-max, and the related local wavenum- ber, k-max [Bergetonet al., 1999]. We note that E-max is sensitive to the variations in the background, while k-max is very usefulin pickingup boundariesof structures[Bergeton et al., 2000]. Seismicvelocityanomaliesare characterizedby Figure 2. Horizontal cut in the center of the box of the analytic (plate a) and noisy (plate b) signal at a SNR of 0.02 db. The width of the cylinder is fixed to 1250 km. The elongated structures visible in the noisy data are artifacts of the method used to generate the noise. The k-max proxy, shown in plate c, does not detect the plume-like structure. The image dimension of each panel is 256x256 pixels coded on a grey scale. BERGERON ET AL- CAPABILITIES OF 3-D WAVELETS The parameter C determines the actual signal-to-noise ra- TO DETECT PLUMES 3435 6 tio (SNR). The measureof the SNR in decibelis given by [Starcketal., 1998]' 5 Zhou 1996 etimated 4 2) SNR (5) SNR(db) = 10log m C; • ilB•--iois•il2 . A Cartesian box Lx x Ly x Lz covers an horizontal area of 10,000km x 10,000km and a height of 700 km. In order to compare this synthetic signal with the resolution of current and near future tomographical models, we used a smaller gridconsisting of 128x 128evenlyspaced gridpointshorizontally and 32 points vertically. -1 • 0 • 500 • 1000 Results of synthetic experiments 1500 • 2000 2500 width (km) We illustrate the method described in section 2 by fixing the width of the cylinder imposing a relatively low SNR ratio of 0.02 db. Fig 2 showsthe three main steps involved in the Figure 4. Sensitivity curve for plume detection. Minimal SNR in decibels versus the plume width. The optimal detection occurs at a width of 1500 km and 0.08 db. procedure: First,thepuresignal (Splume) witherh=625km is defined and displayedon Fig. 2a. Second,the noisy version $measured is shownon Fig. 2b. It is computedby adjust- also find that the optimal detection occurs for a width of ing the C parameter accordingly. The elongated structures visible in the noisy data are artifacts of the method use to generate the noise. Finally, the k-max proxy value is ex- tracted (Fig. 2c). Obviously,the plume-likestructure can not be detected for this SNR value of 0.02 db. We then slightly increase the SNR until k-max detects the large structure. This is accomplished by decreasing the value of the parameter C in equation 4 and by determining whether k-max can detect the large structure in the center of the box. For a synthetic plume with a width of 1250 km, detection occurs for a signal-to-noise ratio of 0.1 db and the whole procedure is summarized on Fig 3. It appears clearly that the large structure(Fig. 3a) can be detectedin k-max space(Fig. 3c) eventhoughit cannotbe detectedby a direct visual inspectionof the noisydata (Fig. 3b). Moreover,it is very difficult to differentiate between the two SNR level as shown on Fig 2b and Fig 3b, since at such low SNR, the noise dominates over the measured signal. Those two examples illustrate our notion of optimal detection. This whole procedure can be repeated for different values of erh. The results represent the sensitivity curves, 1500 km at 0.08 db. This result can be compared with the estimatedSNR levelof 5 db of [Zhou,1996]tomographical model, also given in Fig. 4 as a comparison. It follows that isolated plume-like structure with a typical diameter lying between 250 km and 2250 km, if present in this database, should be detectable. For an implausibly large plume diameter greater than 2250 km, the k-max proxy cannot detect the plume-like structure, even if there is no noise in the signal. This is due to the fact that the horizontal computational domain of 10,000km x 10,000km is too small to encompass fully the convolution of the Gaussian shape with the analyzing wavelet. Thus we give caution to the use of this technique for detecting very large scale structures, which may be present in tomographical data. Larger structures would require a bigger computational domain, which would mean going to spherical geometry. Concluding Remarks Our study with plume-like structures in the upper mantle shows that the sensitivity level obtained with the Mexican displayingthe minimal SNR (maximum noise)and achievhat lies well below the typical SNR ratio estimated for curing detection in k-max space as a function of the structure width 2an. rent tomographicalmodel [Boschiand Dziewonski,1999]. Such a sensitivity curve is plotted in Fig. 4 for plume widths varying between 250 km and 2250 km. We find that km to 2250 km (Fig. 4). Our study alsoshowsthat there is detection is achieved for SNR level of 2.5 db or lower. We This observation holds for plume widths varying from 250 currently no need to employ higher-order Gaussian derivatives for the wavelet function. We have explicitly verified this by using higher-order Gaussian derivatives up to degree 8, without producing any fundamental changes in the results. The same procedure can be applied to other type of 3-D structures, e.g. a slab-like structure, or to different resolu- tion [Bergeronet al., 2000]. From theseresults,we believe Figure 3. Horizontal cut in the center of the box of the analytic (plate a) and noisy (plate b) signal at a SNR of 0.02 db. The width of the cylinder is fixed to 1250 km. The elongated structures visible in the noisy data are artifacts of the method used to generate the noise. The k-max proxy, shown in plate c, detects the plume-like structure. The image dimension of each panel is 256x256 pixels coded on a grey scale. that this method will prove to be even more valuable with larger datasets coming from future high-resolution seismic network(USA net) with a resolution10 timesgreaterthan today's resolution. Acknowledgments. This researchhas been supported by geosciencesprogram of D.O.E. and geophysics program of N.S.F. 3436 BERGERON ET AL.: CAPABILITIES OF 3-D WAVELETS TO DETECT PLUMES References Allen, R. M., G. Nolet, W. J. Morgan, K. Vogfjord, B.H. Bergsson, P. Erlendsson, G.R. Foulger, S. Jakobsdottir, B.R. Julian, M. Pritchard, S. Ragnarsson, and R. Stefansson,The thin hot plume beneath Iceland, Geophys.J. Int., 137, 51-63, 1999. Bergeron, S.Y., Vincent, A.P., Yuen, D.A., Tranchant, J.S. and C. Tchong, Viewing seismicvelocity anomalies with 3-D continuous Gaussian wavelets, Geophys. Res. Lett., 26 No. 15, 2311-231J, 1999. Bergeron, S.Y., Yuen, D.A. and A.P. Vincent, Looking at the Inside of the Earth with 3-D Wavelets: A Pair of New Glasses For Geoscientists,Electronic Geosciences, (œ000)5:3, edited by J.D. Clemens, Springer-Verlag, Heidelberg, Germany, 2000. Bethoux, N. Ouillon, G. and M. Nicolas, The instrumental seismicity of the western Alps: spatial-temporal patterns analyzed with the wavelet transform, Geophils. J. Int., 138, 177-194, Mallat, S., A theory for multiresolution signal decomposition: the wavelet representation, IEEE Tran. on Pattern Anal. and Mach. Intel., 11 (7), 674-693, 1989. Murenzi, R., Antoine, J.P., The continuous wavelet transform, from I to 3 dimensions. A contributed paper in the bookSubband and Wavelet Transforms: Design and Application, Mark J.T. Smith and All Akansu, eds. Kluwer, Dordrecht, 1996. Ni, S., Ding, X., Helmberger, D.V. and M. Gurnis, Lowvelocity structure beneath Africa from forward modelling, Earth Planet. Sci. Left., Vol. 170, 497-507, 1999. Passier, M.L. and R. K. Snieder, On the presenceof intermediatescale heterogeneity in the upper-mantle, Geophys.J. Int., 1œ3, 817 - 837, 1995. Perrier, V., Philipovitch, T. and C. Basdevant, Wavelet spectra compared to Fourier spectra, Journal of Mathematical Physics, 36, No. 3, 1506-1519, 1995. Shen, Y., Solomon, S.C., Bjarnason, I.T. and C.J. Wolfe, Seis1998. mic evidence for a lower-mantle origin of the Iceland plume, Bijwaard, H. and W. Spakman, Tomographic evidence for a narNature, 395, 6œ-65, 1998. row whole mantle plume below Iceland, Earth Planet. Sci. Simons, M. and B.H. Hager, Localization of the gravity field and Lett., 166, 1œ1-1œ6, 1999. the signature of glacial rebound, Nature, 390, 500-504, 1997. Boschi, L. and A.M. Dziewonski, High- and low-resolution images Starck, J.-L., Murtagh, F. and A. Bijaoui, Image Processingand of the Earth's mantle: Implications of different approachesto Data Analysis: The Multiscale Approach, Cambridge Univertomographicmodeling,J. of Geophils.Res., 104(Bll), œ5567sity Press, Cambridge, 1998. œ5594, 1999. Vincent, A.P. and M. Meneguzzi, The spatial structures and staBurke-Hubbard, B., The World according to Wavelets, second tistical properties of homogeneousturbulence, Journal of Fluid edition, A.K. Peters, Wellesley, Ma, 1998. Mechanics, œœ5,1-œ0, 1991. Cadek, O.P., Yuen, D.A., Cizkova, H., Kido, M., Zhou, H.W., Wickerhauser, M.V, Adapted Wavelet Analysis from Theory to Brunet, D. and Ph. Machete], New perspectiveson mantle dySoftware, IEEE Press, Piscatawail, N.J., 1994. namics from high-resolution seismictomographic mode] P1200, Zhou, H, A high resolution P-wave mode] for the top 1200 km of Pure and Applied Geophilsics, 151, 503-5œ5, 1998. the mantle, J. Geophys. Res., 101, œ7791-œ7810,1996. Daubechies, I., Orthonormal bases of compactly supported wavelets, Comm. on Pure and Appl. Math, XLI, 906-996, S. Y. Bergeron and D. A. Yuen, Minnesota Super1988. Daubechies, I., Ten Lectures on Wavelets, SIAM, Philadelphia, computer Institute, University of Minnesota, Minneapo1992. Dziewonski, A.M. and D.L. Anderson, Preliminary reference earth model (PREM), Phys. Earth Planet. Int., œ5,œ97.356, 1981. lis, MN 55415-1227, USA. (e-mail: davey@ krissy.msi.umn.edu) bergeron@msi.umn.edu; A. P. Vincent, D•partement de Physique, Universit• de Montreal, C.P. 6128, Succ. Centre-Ville, Montreal, Quebec, CA) Farge,M., Kevlahan,N., Perrier,V. and]•. Goirand,Wavelets CanadaH3C 3J7. (e-mail: Vincent@ASTRO.UMontreal. and turbulence, Proceedingsof the IEEE, 84, No. 4, 639-669, 1996. Grossmann, A. and J. Morlet, Decomposition of Hardy functions into square integrable wavelets of constant shape, SIAM J. Math. Anal., 15 (4), 7œ3-736,1984. Kumar, P. and E. Foufoula-Georgiou, A multicomponent decomposition of spatial rainfall fields: 2. Self-similarity in fluctua- tions, Water ResourcesResearch,œ9(8), œ533-œ544,1993. (ReceivedNovember18, 1999; revisedMay 22, 2000; acceptedJune 26, 2000.)