April 22, 2004

10:41

WSPC/130-JCA

00223

Journal of Computational Acoustics, Vol. 12, No. 2 (2004) 1–22 c IMACS

ON THE USE OF THE REASSIGNED WAVELET TRANSFORM FOR MODE IDENTIFICATION

MICHAEL I. TAROUDAKIS∗ and GEORGE TZAGKARAKIS Department of Mathematics, University of Crete, 714 09, Heraklion, Crete, Greece ∗ [email protected]

Received 00 Month 2002 Revised 00 Month 2003 This paper is concerned with the use of the reassigned wavelet transform for mode identification in shallow water acoustic propagation. Mode identification is important for inverse procedures in underwater acoustics. An efficient way to recognize the modal structure of the acoustic field when a single hydrophone is available is to refer to the time frequency analysis of the recorded signal using wavelet transform. However, the standard wavelet transform in some cases may result in an obscure representation of the dispersion curves. Thus, a reassigned process is proposed which brings important improvements in the time frequency representation of the signal. This is achieved by moving the calculation point of the scalogram in the center of gravity of the energy concentration, associated with each one of the propagating modes. This argument is supported by two illustrative examples corresponding to propagation of low frequency tomographic signals, in shallow water. Keywords:

1. Introduction Among the most interesting inverse problems in underwater acoustics is that of determining the sea-bed or water column parameters from acoustic measurements obtained in the water column. Ocean observatory systems monitoring the change of oceanographic processes by acoustic means are based on inversion procedures exploiting the information included in the acoustic signals, which are received at a single hydrophone or an array of hydrophones. A vast literature is devoted to mathematical methods and the associated technology for treating the relative problems, yet the field is open for novel techniques of signal processing and inversions, due to the complexity of the inverse problems, being by nature ill-posed. The exploitation of the signal in the time domain was the basis of the traditional Ocean Acoustic Tomography, which is concerned with the problem of estimating the sea water parameters and more specifically the sound speed structure using measurements of the acoustic field at long ranges from a known source. 1 The first inversion procedures in the ∗

Also at the Institute of Applied and Computational Mathematics, FORTH, P.O.Box 1527 711 10 Heraklion, Crete, Greece. 1

April 22, 2004

2

10:41

WSPC/130-JCA

00223

M. I. Taroudakis & G. Tzagkarakis

framework of Ocean Acoustic Tomography used the concept of ray arrivals associating the peaks of the received signal with ray paths before inverting the arrival time information for obtaining the sound speed profile. These applications were limited to deep water situations in order to avoid introducing complications into the inversion procedure due to unknown interactions with unknown sea-bed properties. The input data for the inversion procedure in these cases are arrival times of the signal peaks, which are recognized as “ray arrivals”. Thus, information in the time domain only is exploited and the amplitudes of the peaks are not significant for the inversion procedure. Moreover, a single hydrophone is sufficient as a receiving device and the association of the signal peaks with ray arrivals is done by means of appropriate techniques. As an alternative to ray tomography, modal tomography exploits the modal structure of the acoustic field either in the frequency or in the time domain. The measurable quantity in the frequency domain is the modal phase, provided that an array of hydrophones is available. Examples of modal phase tomography applications can be found in Refs. 2–4. In the time domain, it is frequently possible to identify modal arrivals and base the inversion procedure on the measurements of the modal arrivals in either a linear 5,6 or a nonlinear way.7 Finally methods based on the identification of the local maxima of the signal magnitude (peaks) have also appeared in the literature and have been used with real data. 8,9 Focusing our attention to modal travel time tomography, identification of the modal arrivals can be made using travel time information only, by comparing location of local maxima of the arrival pattern with predicted arrival times of modal packets through an inversion procedure10 based on the central frequency of the broad-band signal. This procedure may result in some cases to false identifications, which can affect seriously the inversion results. The reason is that sometimes there exist local maxima of the arrival pattern not corresponding to a modal packet, which appear in locations very close to a predicted modal packet and they trap towards them the identification process. The identification may be improved by exploiting the time-frequency analysis of the signal, which has the advantage that consistency of the arrival structure over the frequency bandwidth is retained. Time frequency analysis for geoacoustic inversions has been used by several scientists including Dosso and Brooke11 who studied dispersion characteristics of the shear waves and most recently by by Miller and Potty 12,13 who used a hybrid inversion approach, in connection with an optimization process based on a cost function defined over the dispersion properties of the waveguide. Time frequency analysis can be applied by means of Fourier transform or wavelet transform. The latter has several advantages over the Fourier transform when a very wide-band signal is used due to the flexibility in the determination of the time-frequency window. However, even with the wavelet analysis, there exist cases where the energy dispersion for higher order modes is not well reproduced. When these modes are of high importance (case of geoacoustic inversions) their identification is critical and therefore an improvement in the time-frequency analysis is requested. The paper presents an improved wavelet analysis, based on the so called “reassignment” method which has several advantages over the traditional analysis and may lead to identi-

April 22, 2004

10:41

WSPC/130-JCA

00223

Reassigned Wavelet Transform for Mode Identification

3

fication of the higher order modes by refining the representation of the energy dispersion curves. This analysis may be used in connection with an inversion procedure based on the central frequency only10 or on the full bandwidth.12,13 The presentation is supported by two illustrative examples showing a good performance of the method in two simulated shallow-water environments. 2. Modelling the Source-Channel-Receiver System using Normal-Mode Approach The acoustic signals to be considered in the present study are simulated broad-band tomographic signals modelled according to normal-mode theory applied in shallow water environments. It is well known that acoustic energy carried on by acoustic pulses propagates in wave packets each one of which is characterized by individual dispersion properties. The dispersion character of the sound channel can be studied by applying normal-mode theory. A usual way to treat the problem for simulation purposes is to consider the source-channel-receiver system as a linear filter fully described by its impulse response function h SR (t), in the time domain, or, equivalently, by its transfer function H SR (ω), which is the Fourier transform of hSR (t), in the frequency domain: HSR (ω) = F[hSR (t); t → ω].

(2.1)

Actually, the transfer function is referred to a specific source/receiver location and thus it is a function of their spatial coordinates x and x 0 respectively: HSR (x; x0 ; ω). Denoting the source excitation function by S(ω), the pressure field at the receiver has the form: p(x; x0 ; ω) = HSR (x; x0 ; ω)S(ω) .

(2.2)

The transfer function is the solution of the boundary value problem defined by the Helmholtz equation: 1 δ(r)δ(z − z0 ) , (2.3) 2πr supplemented by the appropriate boundary and radiation conditions. Here k(z) is the wavenumber. Using normal-mode theory, the solution at long ranges from the source is given by the formula: ∇2 HSR (r, z; z0 ; ω) + k 2 (z)HSR (r, z; z0 ; ω) = −

π N X ei 4 un (z0 )un (z) ikn r √ HSR (r, z; z0 ; ω) = √ e , kn r 8πρ(z0 ) n=1

(2.4)

where un (z) are the eigenfunctions of the eigenvalue problem defined by the ordinary differential equation d2 un (z) + (k 2 (z) − kn2 )un (z) = 0 , dz 2

(2.5)

April 22, 2004

4

10:41

WSPC/130-JCA

00223

M. I. Taroudakis & G. Tzagkarakis

supplemented by the appropriate boundary conditions, which normally incorporate a pressure release surface, continuity of the function u n (z) and the product ρ−1 (z)dun (z)/dz across the interfaces of the horizontally stratified medium and vanishing of u n (z) as z goes to infinity. ρ(z) is the density which is considered constant in each one of the layers of the medium and kn is the eigenvalue of order n. When a semi-infinite substrate is considered, expression (2.4) implies that only the discrete part of the eigenvalue spectrum is taken into account, the continuous spectrum being assumed negligible at long ranges. Note also that the normalization condition for the eigenfunctions u n (z) involves the weight ρ−1 (z). Thus, the pressure field is: π N X ei 4 un (z0 )un (z) ikn r √ p(r, z; z0 ; ω) = S(ω) √ e kn r 8πρ(z0 ) n=1

=

N X

pn (r, z; z0 ; ω) .

(2.6)

n=1

The pressure field in the time domain is given by the inverse Fourier transform: p(r, z, z0 ; t) = F −1 [p(r, z; z0 ; ω); ω → t] Z ∞ 1 = HSR (r, z; z0 ; ω)S(ω)eiωt dω 2π −∞ Z ∞X N 1 = pn (r, z; z0 ; ω)eiωt dω . 2π −∞

(2.7)

n=1

Interchanging integration with summation, we come up to the result: N Z ∞ X 1 p(r, z; z0 ; t) = pn (r, z; z0 ; ω)eiωt dω 2π n=1 −∞ =

N X

p˜n (r, z; z0 ; t) .

(2.8)

n=1

This expression shows that the acoustic energy is propagating in energy packets each one corresponding to a specific mode number. Taking the absolute value of the pressure field one identifies the location of the maximum energy in each packet among the peaks of this function. This approach has been shown to be appropriate for the simulation and analysis of the propagation of broad-band signals. The dispersion character of the channel is studied by introducing the notion of “group velocity” which is the speed of propagation of the energy packets: ∂ω vg = (2.9) , ∂α ω0 where α is the admissible horizontal wavenumber. 14

April 22, 2004

10:41

WSPC/130-JCA

00223

Reassigned Wavelet Transform for Mode Identification

5

For the eigenvalue problem defined above in Eq. (2.5), the horizontal wavenumbers correspond to the eigenvalues kn which means that a certain mode propagates with a group velocity as defined in Eq. (2.9), where the horizontal wavenumber is the corresponding eigenvalue. ∂ω vg,n = (2.10) . ∂kn ω0

The sound channel is dispersive. This means that the signal as it propagates in the water column changes its shape and characteristics. As regards the presence and behavior of the distinct peaks corresponding to the modal packets, we easily conclude that for a specific source-receiver configuration their amplitude depends on the degree of interference between the individual modal packets. Thus, narrow shapes of modal packets are expected in cases where the group velocity of a specific mode does not change a lot with frequency within the frequency bandwidth of the source. This is easily explained by the fact that if the group velocity does not change a lot, the various frequency components travel with the same speed and therefore there is little pulse spreading. 15 In a general broad-band signal mode dispersion is significant and can be represented by calculating the group velocity as a function of frequency across the frequency bandwidth of the signal. On the other hand, mode dispersion is characteristic for a specific acoustic waveguide and can be used for its recognition provided that acoustic measurements are appropriately processed to obtain the group velocity vs frequency relations for each one of the propagating modes (dispersion curves). Time-frequency analysis is the appropriate means for obtaining the energy distribution over the propagating modes in time and frequency. The arrival times of the modal packets over the frequency bandwidth are directly associated with the corresponding group velocities and thus the dispersion curves are reproduced. Information obtained at the central frequency only, 6 may be not adequate for the characterization of the waveguide especially in cases where the waveguide is characterized by strong dispersion for two reasons: • Association of the peaks of the signal magnitude with modal peaks is not always reliable due to the strong interference characteristics of the signal. Thus we are exposed to the risk of assigning wrong peaks to the modal packets. • As an inverse problem is ill posed, one should exploit the maximum of the available information to obtain reliable estimates of the recoverable parameters. By focusing the inversion at a single frequency, the broad-band character of the signal is not exploited. It remains to decide the best method to apply time-frequency analysis in the measured signal. 3. Wavelets and Computation of the Dispersion Curves A natural way to obtain the information carried by a signal in the time-domain is to map it into the time-frequency plane. The Short-Time Fourier Transform (STFT) is, probably,

April 22, 2004

6

10:41

WSPC/130-JCA

00223

M. I. Taroudakis & G. Tzagkarakis

the most common tool for this purpose. The main drawback of the STFT is the fixed timefrequency resolution over the entire time-frequency plane, since it uses the same window at all frequencies. On the other hand, frequency is directly proportional to the number of cycles per unit time, which means that it takes a narrow window in the time domain to locate high-frequency components and a wide window in the time domain to capture the low-frequency behavior. The Continuous Wavelet Transform (CWT) 16 overcomes the resolution limitation of the STFT, providing a flexible time-frequency window which automatically narrows when observing high-frequency components and widens when observing low-frequency components. The CWT of the signal pR (t) = Real{p(r, z; z0 , t)} relative to a basic wavelet ψ(t) is defined by Z ∞ t − b 1 (Wψ pR )(b, α) = p pR (t)ψ dt , (3.11) α |α| −∞ where pR , ψ ∈ L2 (R), α ∈ R, α 6= 0 is the scale parameter and b ∈ R is the time-location parameter. For the definition of the CWT to be valid, the basic wavelet ψ must satisfy the “admissibility condition” Z ∞ |Ψ(ω)|2 Cψ = dω < ∞ , (3.12) |ω| −∞ where Ψ(ω) is the Fourier transform of ψ(t), and ω denotes circular frequency. By setting t − b 1 , ψb, α (t) = p ψ α |α|

(3.13)

the CWT coefficients becomes

(Wψ pR )(b, α) = hpR , ψb, α i ,

(3.14)

which yields that the CWT gives a measure of similarity between the wavelets ψ b, α and the signal itself. The similarity is in the sense of similar frequency content. The CWT coefficients measure how close to the wavelet, at the current scale, the signal is. In fact the CWT is a time-scale rather than a time-frequency representation. A general relation between scale α and frequency ω is given by α=

ω0 , ω

(3.15)

where ω0 is the central frequency of the basic wavelet’s spectrum, which corresponds to a bandpass filter. The important property of the CWT is that at high frequencies it provides a high time resolution and a low frequency resolution, while at low frequencies the time resolution is low and the frequency resolution is high. This paper demonstrates the effectiveness of this property for the identification of the propagation modes of an underwater acoustic signal. For this purpose it seems to be more

April 22, 2004

10:41

WSPC/130-JCA

00223

Reassigned Wavelet Transform for Mode Identification

7

convenient to use the energy distribution associated with the CWT. This energy distribution is defined by17 (SCψ pR )(b, α) = |(Wψ pR )(b, α)|2

(3.16)

and is called scalogram. Compared to the spectrogram, the scalogram retains the timefrequency resolution property of the CWT. In both, the spectrogram and the scalogram, time and frequency resolutions are bounded by the Heisenberg’s Uncertainty Principle 17 1 , 4 where σt and σω are the standard deviations for time and frequency, respectively. σt2 σω2 ≥

(3.17)

4. The Reassigned Scalogram The scalogram, as every other bilinear distribution, faces the problem of the presence of interference terms (or cross-terms) to those regions of the time-frequency plane where the signal components’ scalograms overlap. These cross-terms may lead to a wrong visual interpretation of the signal’s time-frequency representation as they may overlap with the searched time-frequency pattern (in this study the searched pattern are the modes). Kodera, Gendrin and de Villedary proposed a method to reduce the cross-terms improving the readability, which corresponds to an increase of the signal components concentration.18,19 Their method remained unused because of implementation difficulties. A new formulation of this method, which is called the reassignment method, was presented by Auger and Flandrin,20 leading to its practical use because of its computational efficiency. The new formulation is generalized to any bilinear time-frequency or time-scale distribution. The critical point of the reassignment method is that the modified version of a representation is created by moving its values away from where they are computed, so as to improve the sharpness of the localization of the signal components by reallocating its energy distribution in the time-frequency plane. The scalogram belongs to the class of time-scale distributions. Such representations result performing an affine smoothing of the Wigner–Ville distribution (WVD), as follows 20 : Z ∞Z ∞  u dω , ω0 − αω W V (pR ; t − u, ω) du , (4.18) φT F T SR(pR ; t, α) = α 2π −∞ −∞ where ω0 is the central frequency of the bandpass filter. The last expression shows that the value of a time-scale representation at a point (t, α = ω 0 /ω) is the average of the signal energy located in a domain centered on (t, ω) and bounded by the essential support of φT F . Hence, the time-scale representation can be nonzero on a point (t, ω) where the WVD indicates no energy, if there are some nonzero WVD values at the neighborhood of this point. One way to avoid this, is the assignment of the point on which the average is computed, to the center of gravity ( tˆ, ω ˆ ) of these energy contributions. Auger and Flandrin 20 give analytical expressions for the coordinates of the center of gravity and for the reassigned time-scale representation on a point (t 0 , α0 ).

April 22, 2004

8

10:41

WSPC/130-JCA

00223

M. I. Taroudakis & G. Tzagkarakis

Taking as smoothing kernel φT F the WVD of a window function h(u) yields the scalogram. The coordinates of the center of gravity in this case are given by " # α (W p )(t, α) (W p )(t, α) R R T h h tˆ(pR ; t, ω) = t − Re , (4.19) |(Wh pR )(t, α)|2 " # ω0 (WDh pR )(t, α) (Wh pR )(t, α) ω ˆ (pR ; t, ω) = + Im , (4.20) α α |(Wh pR )(t, α)|2 where T h = t · h(t) and Dh = dh/dt. These coordinates can be computed easily using only two additional scalograms. The reassigned scalogram is defined by Z ∞Z ∞ dt dα (RSCh pR )(t0 , α0 ) = (α0 )2 (SCh pR )(t, α) δ(t0 − tˆ(pR ; t, α)) δ(α0 − α(p ˆ R ; t, α)) α2 −∞ −∞ (4.21) and retains all the properties of the scalogram except the bilinearity, causing the attenuation of the interference terms. The following example demonstrates the difference between the scalogram and its reassigned version. Consider the signal  cos(2π20t), 0 ≤ t ≤ 0.8 y(t) = cos(2π50t), 0.8 ≤ t ≤ 1 with duration equal to 1 sec, which consists of two frequency components, the first is a 20 Hz component in the time interval [0, 0.8] sec and the second is a 50 Hz component in the time interval [0.8, 1] sec. For the calculation of the scalogram and its reassigned version the Morlet’s wavelet will be used. The Morlet’s wavelet is given by the formula ψ(t) = e−t

2 /2

cos 5t

(4.22)

and it has an essential support of [−4, 4]. Figure 1(a) represents the scalogram of the signal y(t). We observe that the time resolution is quite good, since the two components occur at the right time intervals. However the frequency resolution is bad, since the scalogram indicates energy in two frequency bands and not at the two certain frequencies (the first frequency band is around the 20 Hz component and the second is around the 50 Hz component). An interference term appears at the time instant 0.8 sec, because of the sudden change of the frequency at that point. Figure 1(b) shows the reassigned scalogram of the signal y(t) using the Morlet’s wavelet. The time and frequency resolution is almost perfect. The energy is concentrated in the frequencies of 20 and 50 Hz and the interference term has been attenuated. 5. Applications The following test cases demonstrate the effectiveness of using the reassigned scalogram for mode identification in shallow water acoustic transmissions. It is noted that the analysis

April 22, 2004

10:41

WSPC/130-JCA

00223

Reassigned Wavelet Transform for Mode Identification

(a) Scalogram (a) Scalogram.

(a) Scalogram

(b) Reassigned Scalogram.

(b) and Reassigned Scalogram Fig. 1. Scalogram reassigned scalogram for the function y(t). Figure 1. Scalogram and reassigned scalogram for the function y(t). (b) Reassigned Scalogram Figure 1. Scalogram and reassigned scalogram for the function y(t).

9

April 22, 2004

10

10:41

WSPC/130-JCA

00223

M. I. Taroudakis & G. Tzagkarakis

(0,0)

Sea surface

Water

r

!w

,

cw (z)

source

receiver

!

" (r,z) sound speed profile

D Sediment

!b , cb (z)

Substrate

!sb , csb

D+ h

z Fig. 2. The test environment.

Figure 2. The test environment

presented here, corresponds to the forward problem only, as the signals considered are based on simulations and the dispersion characteristics are known. Figure 2 presents the sea environment of the test cases, which is assumed range independent and axially symmetric. It consists of a shallow water layer, a single sediment layer and a semi-infinite bottom (the substrate). All layers are considered fluid and the sound speed profile may vary with depth in the water layer and the sediment, while it is constant in the substrate. For simplicity, the density of each layer is assumed to be constant. In both cases a low frequency tomographic signal will be considered. 5.1. Test case I In the first test case the tomographic signal will be modelled in the frequency domain using a Gaussian pulse of central frequency f 0 = 50 Hz and a bandwidth ∆f = 20 Hz placed at the depth of 30 m. The environmental parameters appear in Table 1. The simulated data correspond to measurements at the depth of 30 m and range of 100 km from the source. The data are calculated using the Normal-Mode program MODE1 developed at F.O.R.T.H. These data are provided as input to the discrete inverse Fourier transform to yield the signal in the time domain which is displayed in Fig. 3. The magnitude of the received signal is presented in Fig. 4 and we observe several peaks. The identification of the peaks that correspond to modal arrivals is in principle possible

April 22, 2004

10:41

WSPC/130-JCA

00223

Reassigned Wavelet Transform for Mode Identification

11

Fig. 3. The simulated signal of Test Case I.

Figure 3. The simulated signal of Test Case I.

Fig. 4. The arrival pattern of the simulated signal of Test Case I.

Figure 4. The arrival pattern of the simulated signal of Test Case I.

using the identification procedure proposed by Taroudakis. 10 This identification algorithm based on the determination of the group velocity at the central frequency only is easily applied and it is very fast as there is no need to calculate the field in the time domain.

April 22, 2004

12

10:41

WSPC/130-JCA

00223

M. I. Taroudakis & G. Tzagkarakis

Table 1. The shallow water environment (Test Case I). Water Depth Sound speed profile in the water: cw (0) cw (30) cw (200) Compressional velocity in the sea-bed Sea-bed density

200 m 1500 1498 1500 1800 1500

m/sec m/sec m/sec m/sec kg/m3

Fig. 5. The dispersion curves of the signal Test Case I.

Figure 5. The dispersion curves of the signal of Test Case I

However, by closer inspection of the simulated arrival pattern it becomes clear that the presence of many peaks renders the mode identification based at the central frequency only, very tedious and may not be accurate. For this reason the identification approach based on the scalogram or the reassigned scalogram seems to be promising, since the received signal is analyzed at a whole frequency band. The feature to be exploited is the dispersion characteristics of the signal and not the arrival pattern. Figure 5 presents the dispersion curves of the simulated signal as recorded at the receiver, in the time-frequency plane. The dispersion curves are usually presented in the group velocity-frequency plane but it is very easy to display them in the time-frequency plane using

April 22, 2004

10:41

WSPC/130-JCA

00223

Reassigned Wavelet Transform for Mode Identification

13

MISC f0T = 2.3, Nf = 400, Threshold = 0.04%

Fig. 6. The scalogram and corresponding dispersion curves for the signal of Test Case I.

the expression

Figure 6. The scalogram and corresponding dispersion curves for the signal of Test r tn = Case , I (5.23)

vg, n

where tn is the arrival time, r is the range and v g, n is the group velocity of a modal packet. The numbers in the dispersion curves correspond to mode order. Our purpose is to demonstrate the effectiveness of using the reassigned scalogram for the identification of the propagating modes. We have used the Morlet’s wavelet since it has an explicit formula and it is symmetric, which means that beginning the computations for the analysis of the signal from the left-hand side, yields the same results with those obtained when starting the computations from the right-hand side. Figure 6 presents the scalogram of the analyzed signal, relative to the Morlet’s wavelet, together with the analytical dispersion curves. The computation has been performed using 800 frequency bins and a time-bandwidth product of the Morlet’s wavelet equal to 2.3. Table 2 presents the arrival time of each mode at the central frequency. The time differences between the modes are large enough, so we expect that the modes will be easily identified. Indeed, the first six modes are clearly visible. Around the 7th mode the scalogram indicates no energy, which means that this mode is not excited. This is easily explained by looking at the shape of the mode 7 at the central frequency, showing a null at the depth of the receiver (see Fig. 7). An interesting observation is that the scalogram indicates energy for frequencies above 50 Hz and times greater than 74 sec, which corresponds to traces of the last three modes (8th, 9th, 10th). These modes are not excited at the central frequency and thus they cannot

April 22, 2004

14

10:41

WSPC/130-JCA

00223

M. I. Taroudakis & G. Tzagkarakis

Table 2. Arrival time and arrival time differences between two adjacent modes, for the identified modes of Test Case I. Central frequency f0 = 50 Hz. Mode no.

Arrival time (sec)

Arrival time difference (sec)

1 2 3 4 5 6 7

66.8664 67.3104 68.0889 69.2459 70.8430 72.9440 75.5011

0.4439 0.7784 1.1569 1.5971 2.1009 2.5571

Fig. 7. Amplitude of mode 7 of the signal corresponding to Test Case I.

be exploited for inversion purposes the7 inversion is based onCase information on the Figure 7. Amplitude of if mode of the signalscheme corresponding to test I central frequency only. Although there is a very good agreement between the experimental (wavelet analysis) and the analytical (dispersion curves calculation) results, the time-frequency resolution is not so good as someone would expect. Thus, there is still some uncertainty in any inversion process based on the comparison of the predicted dispersion curves with the energy distribution provided by the scalogram. It is interesting to note that the standard Fourier Transform gives a representation which is very similar to the scalogram provided by the original wavelet transform (see Fig. 8), which is expected in our case, as the signal used is narrow-band and for these types of signals the benefits of the wavelet transform are not fully exploited.

April 22, 2004

10:41

WSPC/130-JCA

00223

Reassigned Wavelet Transform for Mode Identification

15

Fig. 8. The spectrogram (STFT) and corresponding dispersion curves for the signal of Test Case I.

Figure 8. The spectrogram and= corresponding dispersion RMSC f0T(STFT) = 2.3, Nf 400, Threshold = 0.004% curves for the signal of Test Case I

Fig. 9. The reassigned scalogram and corresponding dispersion curves for the signal of Test Case I.

Figure 9. The reassigned scalogram and corresponding dispersion curves for the signal of Test case I

One way to improve the time-frequency resolution providing a better visual interpretation, is to use the reassigned scalogram. As it is presented in Fig. 9 the reassigned scalogram gives much better results than the scalogram. Now, the first six reassigned experimental

April 22, 2004

16

10:41

WSPC/130-JCA

00223

M. I. Taroudakis & G. Tzagkarakis

Fig. 10. The reassigned scalogram for the signal of Test Case I when attenuation of 0.015 nep/m KHz is introduced in the bottom.

Figure 10. The reassigned scalogram for the signal of Test Case I when attenuation of 0.015 nep/m KHz is introduced in the bottom

modes are perfectly localized to the corresponding analytical dispersion curves. The 7th mode remains invisible as before, something natural as the reassignment procedure does not reassign a point where the representation value is zero. The identification of the last three modes although still poor is significantly better after their reassignment. By introducing some attenuation in the bottom, higher order modes are in general damped. The reassigned scalogram is capable of identifying higher order modes even in this case, provided that the absorption is of course reasonable. Figure 10 presents the reassigned scalogram of Test Case 1 when an attenuation coefficient of 0.015 nep/m Khz is introduced in the bottom properties. Modes 8 and 9 are still visible although with low energy. By comparison with Fig. 11 representing the original scalogram, we observe the improvement in the representation of all the modes including the higher order ones. 5.2. Test case II The second test case also corresponds to a shallow water environment. The geoacoustic parameters of the environment appear in Table 3. The signal emitted by a source, placed at the depth of 75 m, is simulated using a Gaussian pulse of central frequency f 0 = 150 Hz and a bandwidth ∆f = 40 Hz. Figure 12 presents the arrival pattern at a receiver placed at the depth of 75 m and range of 20 km from the source. The magnitude of the signal is shown in Fig. 13, while Fig. 14 presents the dispersion curves for this environment (numbers in the curves correspond to mode order). The travel time of each mode at the central frequency appear in Table 4. The travel time differences for the identified modes indicate that it will

April 22, 2004

10:41

WSPC/130-JCA

00223

Reassigned Wavelet Transform for Mode Identification

17

Fig. 11. The scalogram for the signal of Test Case I when attenuation of 0.015 nep/m KHz is introduced in the bottom.

Figure 11. The scalogram for the signal of Test Case I when attenuation of 0.015 nep/m KHz is introduced in the bottom

Fig. 12. The simulated signal of Test Case II.

Figure 12. The simulated signal of Test Case II

April 22, 2004

18

10:41

WSPC/130-JCA

00223

M. I. Taroudakis & G. Tzagkarakis

Table 3. The shallow water environment (Test Case II). Water Depth Sound speed profile in the water: cw (0) cw (100) Compressional velocity in the sea-bed Sea-bed density

100 m 1500 1520 1590 1200

m/sec m/sec m/sec kg/m3

Table 4. Arrival time and arrival time differences between two adjacent modes, for the identified modes of Test Case II. Central frequency f0 = 150 Hz. Mode no.

Arrival time (sec)

Arrival time difference (sec)

1 2 3 4 5 6

13.3042 13.2986 13.3528 13.4488 13.5690 13.6919

0.0056 0.0542 0.0959 0.1201 0.1229

Fig. 13. The arrival pattern of the simulated signal of Test Case II.

Figure 13. The arrival pattern of the simulated signal of Test Case II be difficult to distinguish the first three modes, while the fact that the dispersion curves are relatively vertical indicates sharp modal arrivals. The time-frequency content of the analyzed signal is extracted using the scalogram and

April 22, 2004

10:41

WSPC/130-JCA

00223

Reassigned Wavelet Transform for Mode Identification

19

Fig. 14. The dispersion curves for the signal of Test Case II.

Figure 14. The dispersion curves for the signal of Test Case II MISC f0T = 2.7, Nf = 256, Threshold = 0.04%

Fig. 15. The scalogram and corresponding dispersion curves for the signal of Test Case II.

Figure 15. The scalogram and corresponding dispersion curves for the signal of test Case II its reassigned version. Figure 15 shows the scalogram, which has been calculated using 512 frequency bins and the Morlet’s wavelet with time-bandwidth product equal to 2.7. The first two modes is impossible to be distinguished by the scalogram, since they are very close to each other, while the third mode is slightly distinguished from the first two modes. The scalogram indicates no energy around the 4th and the 7th analytical dispersion curves,

April 22, 2004

20

10:41

WSPC/130-JCA

00223

M. I. Taroudakis & G. Tzagkarakis

Fig. 16. Amplitude of mode 4 of the signal corresponding to Test Case II.

Figure 16. Amplitude of mode 4 of the signal corresponding to test Case II RMSC f0T = 2.3, Nf = 256, Threshold = 0.04%

Fig. 17. The reassigned scalogram and corresponding dispersion curves for the signal of Test Case II.

Figure 17. The reassigned scalogram and corresponding dispersion curves for the signal of test Case II which means that the corresponding modes are not excited. Note that as it is indicated by Fig. 16 mode 4 is not excited due to the position of the source at a null of the corresponding function, whereas mode 7 is not excited at all at the central frequency. Finally, it is possible to identify the 5th and the 6th modes. Once again the time-frequency resolution provided

April 22, 2004

10:41

WSPC/130-JCA

00223

Reassigned Wavelet Transform for Mode Identification

21

by the scalogram is not good, although there is a relatively good agreement between the experimental and the analytical results for the modes 3, 5, 6. The reassigned scalogram is used to improve the time-frequency resolution, as in the previous test case. The reassigned modes are shown in Fig. 17. We observe that the first two modes are still impossible to be distinguished, since their travel time difference is almost zero. On the other hand, the third mode is clearly identified after its reassignment and it can be distinguished from the first two modes. The reassigned scalogram also indicates that the modes 4 and 7 are not excited. Although there are some interference terms between the 5th and the 6th mode, their time-frequency resolution has been improved. In general, the reassigned scalogram provides a perfect agreement between the experimental modes and the analytical curves for this environment. This test case also supports our statement that the reassigned scalogram is capable of distinguishing closely spaced modes of a multimode acoustic signal and may be proven a useful tool for inversion schemes based on the identification of the signal modal structure. 6. Conclusions The use of the reassigned scalogram is shown to improve the time-frequency representation of an acoustic signal recorded in shallow water. Improvement is evident for the identification of all the propagating modes. However for geoacoustic inversions which in principle are based on higher order modes, the use of the reassigned scalogram is more beneficial, as these modes are normally the more difficult to be identified using standard time-frequency analysis or identification procedures based on the central frequency only. As it was illustrated by means of two characteristic examples based on simulated signal recordings, the reassigned scalogram resulted in a more clear representation of the dispersion curves for these modes. Future work includes the development of an automatic inverse procedure based on the use of the reassigned scalogram and study its efficiency for a wide range of different signals and shallow water environments. References 1. W. Munk and C. Wunsch, Ocean acoustic tomography: a scheme for large scale monitoring, Deep Sea Res. 26a (1979) 123–161. 2. S. Rajan, J. F. Lynch, and G. V. Frisk, Perturbative inversion methods for obtaining bottom geoacoustic parameters in shallow water, J. Acoust. Soc. Am. 82 (1987) 998–1017. 3. E. C. Shang, “Ocean acoustic tomography based on adiabatic mode theory, J. Acoust. Soc. Am. 85 (1989) 1531–1537. 4. M. I. Taroudakis and J. S. Papadakis, A modal inversion scheme for ocean acoustic tomography, J. Comp. Acoustics 1 (1993) 395–421. 5. E. C. Shang and Y. Y. Wang, On the possibility of monitoring the El Ni˜ no by using modal ocean acoustic tomography, J. Acoust. Soc. Am. 91 (1992) 136–140. 6. M. I. Taroudakis, A comparison of the modal-phase and modal- travel time approaches for ocean acoustic tomography, in Proc. 2nd European Conf. Underwater Acoustics, ed. L. Bjørnø, 1994, pp. 1057–1062.

April 22, 2004

22

10:41

WSPC/130-JCA

00223

M. I. Taroudakis & G. Tzagkarakis

7. M. I. Taroudakis and M. Markaki, Tomographic inversions in shallow water, using modal travel time measurements, Acta Acustica United with ACUSTICA 87 (2001) 647–658. 8. E. K. Skarsoulis, G. A. Athanassoulis and U. Send, Ocean acoustic tomography based on peak arrivals, J. Acoust. Soc. Am. 100 (1996) 797–813. 9. G. A. Athanassoulis, J. S. Papadakis, E. K. Skarsoulis and M. I. Taroudakis, A comparative study of modal and correlation arrival inversion in ocean acoustic tomography, in Full Field Inversion Methods in Ocean and Seismic Acoustics, eds. O. Diachok, A. Caiti, P. Gerstoft and H. Schmidt (Kluwer Academic Publishers, 1995), pp. 127–132. 10. M. I. Taroudakis, Identifying modal arrivals in shallow water for bottom geoacoustic inversions, J. Comp. Acoustics 8 (2000) 307–324. 11. S. E. Dosso and G. H. Brooke, Measurement of seismo-acoustic ocean-bottom properties in the high Arctic, J. Acoust. Soc. Am. 98 (1995) 1657–1666. 12. G. R. Potty and J. H. Miller, Geoacoustic tomography: range dependent inversions on a single slice, J. Comp. Acoustics 8 (2000) 325–346. 13. G. R. Potty, J. H. Miller, J. F. Lynch and K. B. Smith, Tomographic inversion for sediment parameters in shallow water, J. Acoust. Soc. Am. 108 (2000) 973–986. 14. I. Tolstoy and C. S. Clay Ocean Acoustics. Theory and Experiment in Underwater Sound (American Institute of Physics, New York, 1987). 15. C. T. Tindle, K. M. Guthrie, G. E. J. Bold, M. D. Johns, D. Jones, K. O. Dixon and T. G. Birdsall, Measurements of the frequency dependence of normal modes, J. Acoust. Soc. Am. 64 (1978) 1178–1185. 16. C. K. Chui, An Introduction to Wavelets (Academic Press, 1992). 17. S. Mallat, A Wavelet Tour of Signal Processing (Academic Press, 1998). 18. K. Kodera, C. de Villedary and R. Gendarin, A new method for the numerical analysis of time varying signals with small BT values, Phys. Earth Planet Interiors 12 (1976) 142–150. 19. K. Kodera, R. Gendarin and C. de Villedary, Analysis of time-varying signals with small BT values, IEEE Trans. Acoust. Speech Signal Processing ASSP-34 (1986) 64–76. 20. F. Auger and P. Flandrin, Improving the readability of time-frequency and time-scale representations by the reassignment method, IEEE Trans. Signal Processing 43 (1995) 1068–1089.