Matthieu Loosvelt

Philippe Lasaygues
email: lasaygues@lma.cnrs-mrs.fr

A wavelet-based processing method for simultaneously determining ultrasonic velocity and material thickness

Keywords: numbers: 43.60.-c, 43.60.+d, 43.35.Ae, 43.35.Bf, 43.35.Cg Ultrasonic wave velocity, thickness, simultaneous measurement, wavelet processing

Methods of measuring ultrasonic wave velocity in an elastic sample require data on the thickness of the sample and/or the distances between the transducers and the sample. The uncertainty of the ultrasonic wave velocity measurements generally depends on that of the data available. Conversely, to determine the thickness of a material, it is necessary to have a priori information about the wave velocity. This problem is particularly hard to solve when measuring the parameters of biological specimens such as bones having a greater acoustical impedance contrast (typically 3 to 5 MRayl) than that of the surrounding soft tissues (typically 1.5 MRayl). Measurements of this kind cannot easily be performed. But obtaining the thickness of a bone structure and/or the ultrasonic wave velocity is a important problem, for example, in biomechanical field for the calculation of elastic modulus, or in acoustical imaging field to parameterize the images, and to reference the grey or color level set to a physical parameter.

The aim of the present study was to develop a method of simultaneously and independently determining the velocity of an ultrasonic wave in an elastic sample and the wave path across the thickness of this sample, using only one acquisition in pure transmission mode. The new method, which we have called the "Wavelet-Based Processing" method, is based on the wavelet decomposition of the signals and on a suitable transmitted incident wave correlated with the experimental device, and the mathematical properties such as orthonormality, of which lend themselves well to the time-scale approach. By following an adapted algorithm, ultrasonic wave velocities in parallelepipedic plates of elastic manufactured material and the apparent thicknesses were both measured using a water tank, a mechanical device and a matched pair of 1 MHz ultrasonic focused transducers having a diameter of 3 mm, a focal length of 150 mm and beam width of 2x2 mm at the focus (mean temperature 22°). The results were compared with those obtained with a conventional Pulse-mode method and with the control values, to check their validity. Measurements performed on bovine and human dry cortical bone samples are also presented to assess the limitations of the method when it is applied to elastic biological samples, including those of an equal-wavelength size (≈ 1.5 mm). The thicknesses and the ultrasonic wave velocities were then measured in this kind of (quasi-)parallelepipedic elastic materials with an mean estimated error ranged from 1% to 3.5% compared to the referenced values.

Introduction

The velocity of ultrasonic waves propagating in a solid elastic material is an important property that can be used, for example, in non-destructive evaluation to calculate the stiffness rigidity matrix, or in the biomedical field to parametrize reconstructed images (in ultrasonic tomography [START_REF] Lasaygues | High Resolution Low Frequency Ultrasonic Tomography[END_REF]).

Conventional techniques used for making velocity measurements in materials require information about the thickness of the material under investigation. These thickness measurements are usually calculated from time-of-flight measurements. Depending on the experimental equipment used, this can involve complicated procedures (several data statements, several acquisition mode, several transducers...), which impose fundamental efficiency limitations (error higher than 5%, reproducibility…) on the velocity measurements. In the case of in vitro ultrasound methods of characterizing bone structures for relevant works in biomechanical field (evaluation of elastic modulus for example) or acoustical field (quantitative ultrasonic imaging for example), the occurrence of physical processes associated with the wave propagation generates complex acoustic signals (weak Signal-to-Noise Ratio (SNR), multiple echoes, dispersion…), which it is often difficult to analyze and to interpret in terms of wave paths or depth-dependencies. The main problem arising here is that it is very difficult to obtain samples with a regular thickness and parallel faces. A large proportion of the incident ultrasonic energy is also scattered (or simply reflected) and the transmitted proportion is refracted (deviated) after undergoing modal conversion at the water/sample interface. To overcome these problems, the use of low ultrasonic frequencies (< 3 MHz) provides an effective alternative approach because some of the zones turn out to be more homogeneous, and the quantity of energy lost during the propagation of the wave therefore decreases. But even with low frequencies, to be able to determine the time-of-flight and thickness, it is necessary to find new methods of recording and processing these signals, and to develop more efficient algorithms.

Several approaches have been used so far for this purpose, such as filtering and spectral analysis methods and a method involving deconvolution based on the use of a transfer function characteristic of the experimental device [START_REF] Lasaygues | High Resolution Low Frequency Ultrasonic Tomography[END_REF]. Deconvolution algorithms improve the quality of the information in the low ( 0) and high ( ∝) frequency ranges, but tend to be over-sophisticated and unstable. Indeed, they are a class of inverse problems associated with regularization process and ill-conditioned systems, and the main drawbacks of the most efficient algorithms are their computational time consuming because they involve often inversion of huge, full or complex matrix, which may not be compatible with an automatic process involving a large number of data.

Although there exist optimized [START_REF] Shin | Sensitivity to point-spread function parameters in medical ultrasound image deconvolution[END_REF], faster [START_REF] Preza | Image reconstruction for 3-D light microscopy with a regularized linear method incorporating smoothness prior[END_REF][START_REF] Haeberlé | Improving the lateral resolution in confocal fluorescence microscopy using laterally interfering excitation beams[END_REF], and adaptive [START_REF] Snoussi | Bayesian blind separation of generalized hyperbolic processes in noisy and under determinate mixtures[END_REF] methods, deconvolution algorithms are sensitive to artifacts, noise and the bias affecting the repeatability of the measurements. Alternative methods based on multi-scale decomposition procedures, such as those based on wavelet transformation of the signals [START_REF] Donoho | Nonlinear solution of linear inverse problems by wavelet-vaguelette decomposition[END_REF], make it possible to process all the local information and much more of the overall information available in terms of the frequency and time parameters. Instead of performing regularization in order to reduce the effects of inverse filter singularities and to restore the original signal, a wavelet transform can be performed in the deconvolution procedure by regularizing the inverse problem [START_REF] Neelamani | Wavelet-based deconvolution for ill-conditioned systems[END_REF][START_REF] Wan | Robust Deconvolution of highfrequency ultrasound images using higher-order spectral analysis and wavelets[END_REF]. Although this method is highly efficient, it involves performing delicate processing steps and having prior knowledge of the singularities of the signal. When used alone, the wavelet transform method lends itself very well to detecting and discriminating between signals during the data pre-processing phase and to extracting information such as the instantaneous frequency and the evanescent properties of the medium [START_REF] Saracco | Propagation of transient waves through a stratified fluid medium: wavelet analysis of a non asymptotic decomposition of the propagator. Part 1. Spherical waves through a two-layered system[END_REF] to remove the tiny echoes embedded in a strong noise background [START_REF] Castillo | Ultrasonic flaw detection in cast iron using a multi-resolution analysis[END_REF] and to abolish speckle noise [START_REF] Abbate | Signal detection and noise suppression using a wavelet transform signal processor: Application to ultrasonic flaw detection[END_REF]. Another efficient wavelet method is also used to assess the spatial changes occurring in a multi-scale structure [START_REF] Grossmann | Decomposition of hardy functions into square integrable wavelets of constant shape[END_REF] and to determine the thickness of irregular plates and granular interfaces [START_REF] Le Gonidec | Multiscale analysis of waves reflected by complex interfaces: Basic principles and experiments[END_REF] from the reflected acoustical waves. The main advantage here is the possibility of performing the optimized shaping of the signal associated with the incident wave propagating through the media [START_REF] Conil | Non-linear synthesis of input signals in ultrasonic experimental setups[END_REF] and applying a matching process, using the wavelet's mathematical properties. At the experimental level, however, the so-called "wavelet response" method is sometimes poorly suited to automatically operating control situations because it requires the use of many probes with different nominal frequencies corresponding to the wavelet scales, and because it focuses entirely on the reflected wave.

Here we present a new method of simultaneously determining the ultrasonic transmitted wave velocity in an elastic material and the associated thickness without any previous knowledge of either parameter. The new algorithm is based on a wavelet decomposition method, which is applied to the acoustical signals, and on a suitable transmitted signal associated with the incident wave, and correlated with the parameters of the experimental device. The "Wavelet-Based Processing" method (WBP method) uses the mathematical properties of these acoustical signals to measure the time-offlight of the transmitted wave through the sample in a one-shot transmission mode, and then to calculate the thickness and the velocity simultaneously. This method, which does not require the use of any specialized equipment, allows discriminate placement of the sample in the ultrasonic path.

The validity of the results obtained was confirmed by measuring the thickness and the wave velocity in rectangular parallelepipedic manufactured plates, the mean thickness of which was closely determined using calipers to within 0.1 mm. The mean control thicknesses and velocities obtained using a more conventional broadband pulse method are also presented, and cases involving thicknesses of the same order as the wavelength (which will be called the "equal-wavelength" case from now on) will be discussed. Lastly, similar experiments and results are presented on bovine and juvenile human bone samples to assess the limitations of the WBP method as a means of biological tissue characterization.

Material and methods

The main working hypothesis adopted in this in vitro experimental study was that the sample was comparable to a parallelepipedic rectangle, and that the ultrasonic incident wave vector was perpendicular to the water/sample interface. Under these assumptions, pure compression waves (ultrasonic waves) were therefore transmitted through water and the sample. Any shear waves, and compression-to-shear wave mode conversion in the sample were neglected. The sample was taken to be homogeneous and isotropic. The wave velocities were taken to be constant and independent of the frequency (non-dispersive). Only the propagation processes were taken into account. The ultrasonic wave attenuation resulting from absorption processes was assumed to be weak and the magnitude of the signals decreased by only a few percent during the propagation. Only the time-offlight (TOF) of the waves was therefore measured.

To have several measurement points on the sample, linear scanning was performed with the transducers over a distance of L mm (16 mm in this study) with a ΔL-mm (1 mm) step on both sides of the sample, using an ultrasonic mechanical scanner. At each measurement, the depth of the wave path versus the mean thickness (e) of the sample in the sounded area (which will be called "the thickness" from now on), and the corresponding ultrasonic wave velocity v b (which will be called "the velocity" from now on) were calculated using the Wavelet-Based Processing (WBP) method and a conventional ultrasonic pulse-mode (Pulse-mode) method as the control method. The thicknesses were also compared with those measured using a caliper. With both methods, a digitized initial signal was obtained without any target in order to measure the TOF t 0 (µs) and the velocity v 0 (m.s -1 ) of the ultrasonic wave, given the distance (d) between the transmitter and the receiver. The data obtained in the first experiment were subsequently used to determine the linear mechanical variations occurring in the ultrasonic mechanical scanner due to the mechanical offset effects.

2.1.Ultrasonic mechanical scanner

Measurements were performed on the ultrasonic mechanical scanner developed at the Laboratory of Mechanics and Acoustics (LMA). The apparatus used consisted of a main symmetrical arm carrying two secondary transversal arms. Each arm supported an aligned transducer, which could be moved linearly. Increments were multiples of 0.75 hundredths of millimeters. The object to be tested was placed in the presumed geometrical center of the bench so that the maximum distance between the transducers and the center was 150 mm. The surrounding fluid medium was water at a temperature of (Θ) (degrees). The matched pair of transducers used for data acquisition purposes were piezocomposite transducers (Imasonic ® ) with a nominal frequency F c of 1 MHz (bandwidth ranged from 0.5 MHz to 1.5 MHz). In this case, the wavelength was around 1.5 mm in water (≈ 1500 m.s -1 ). The diameter of the active interface was 3 mm, and the lateral (angular) and vertical (slice) resolution of the transducer beam was 2 mm in each case with the focus at 150 mm. This distance was corresponded in the center between the both transducers where the interface water/sample was placed (distance l 1 , Figs. 2 and3). The transducers were driven using a waveform generator (TTI ® TGA 1241), and positioned automatically as required. The radio-frequency signals (RF-signals) stored (Agilent ® DSO5014A) were then used to determine the TOF between the source, the object and the receiver. The signal processing algorithms were implemented on a personal computer.

Since the electro-acoustic device (Fig. 1) and the transducers therefore served as a continuous linear stationary causal filter, the input x(t) and output s(t) signals were connected by convolution:

s(t) = (x ⊗ h M )(t) (1) and € x(t) = (h T t ⊗ h T r ⊗ e)(t) = (h T * ⊗ e)(t) (2)
where € ⊗ denotes the operation of convolution. e(t) is the electric signal conveyed to the transmitter via the waveform generator, € h T t (t) and € h T r (t) are the responses of the transmitter and the receiver, which are assumed to be similar, and (t) is the response of the object.

€ h T * (t) = (h T t ⊗ h T r )(t), € h M
In the transmission mode, without the sample, the transmitted response € h M (t) in water depends, when there are no distortions of the propagation, on the time-delay of the wave, which is proportional to the distance (d) between the two transducers and the velocity v 0 . The initial output signal

€ s(t) [ ]
init is therefore equal to the input signal x(t), which is invariant by translation:

€ s(t) [ ] init = x(t) ⊗ δ t - d v 0 % & ' ( ) * = x t - d v 0 % & ' ( ) * = x(t) = h T * ⊗ e ( ) (t) (3)
Note that at this stage, the responses of the transducers do not have to be determined.

2.2.Conventional ultrasonic Pulse-mode method

To obtain control experimental data, the thicknesses and velocities were first calculated using a conventional ultrasonic Pulse-mode method. Three measurements were necessary in the presence of the sample: two in the reflected mode and one in the transmission mode (Fig. 2). The electric signal e(t) was a pulse signal, which was comparable to a Dirac delta function (in terms of the distribution), and

€ s(t) [ ]
init is given by:

€ s(t) [ ] init = h T * (t) (4)
Experiments were carried out in the reflected mode in order to determine TOF t 1 (respectively TOF t 2 ) of the propagating wave, from the left (respectively right) transducer to the left (respectively right) interface of the sample.

The following relation gives the thickness (e):

€ e = d -v 0 t 1 + t 2 2 # $ % & ' ( (5)
The velocity v b was determined in the transmission mode at the normal incidence:

€ v b 2 = v 0 2 1+ v 0 τ e v 0 τ e -2 $ % & ' ( ) (6)
where ( € τ) is a time-delay calculated by cross-correlating the initial signal obtained without the sample

€ s(t) [ ]
init with the output signal obtained with the sample s(t). With this method, the formulas for the thickness and the velocity are directly related because (e) has to be determined before v b .

2.3.Wavelet-Based Processing (WBP) 2.3.1.One-shot transmission mode

With the WBP method, only a one-shot acquisition of the output signal s(t) in the pure transmission mode is necessary, and no reflection mode acquisitions are required (Fig. 3).

A transmission signal is digitized in order to simultaneously measure the first TOF t 1 of the transmitted wave crossing distances l 1 and l 2 in water and the distance (e) in the sample; and the second TOF t 2 of the transmitted wave crossing distances l 1 and l 2 in water and the distance (3e) in the sample (back-wall echo):

€ t 1 = l 1 + l 2 ( ) v 0 + e v b (7)
€

t 2 = l 1 + l 2 ( ) v 0 + 3 e v b (8)
At each measurement, provided (t 0 , t 1 , t 2 ) and (v 0 ) are known, the thickness (e) can be written:

€ e = v 0 t 0 - 3t 1 -t 2 ( ) 2 # $ % % & ' ( ( (9)
and the ultrasonic wave velocity v b can be written:

€ v b = v 0 2t 0 -3t 1 + t 2 ( ) t 2 -t 1 ( ) (10)
In this procedure, the thickness and the ultrasonic wave velocity are calculated independently, and can therefore be obtained simultaneously. The main problem is how to determine TOF t 1 and TOF t 2 , even in the case of an equal-wavelength thickness. An orthogonal wavelet decomposition of the output signal s(t) and an adaptation of the input signal x(t) in term of time and frequency (see paragraph 2.3.5), can be performed to automatically determine these times.

2.3.2.Orthogonal wavelet decomposition (Meyer-Jaffard algorithm)

The wavelet algorithm specifically developed for ultrasonic signal processing purposes [START_REF] Lasaygues | Deconvolution and Wavelet Analysis on Ultrasonic Reflexion Tomography[END_REF] calculates the wavelet coefficients of a signal in the time and frequency bands. Wavelet decomposition of signals is a transformation that depends on discrete parameters [START_REF] Meyer | Orthonormal wavelets[END_REF] noted j and k, (j, k) ∈ Z 2 . Indeed, we associated the series of coefficients C j,k with a signal s(t) as follows (the overline indicates the conjugate complex):

€ C j,k = s(t) -∞ ∞ ∫ Ψ j,k (t) dt (11)
€ Ψ j,k (t) = 2 j Ψ 2 j t - k 2 j Δt % & ' ( ) * + , - . / 0 (12)
where (j, k) belongs to Z 2 , and € Δt is the sampling period. When j < 0, we have dilatation and when j > 0, we have compression of the function

€ Ψ (t).
The function (t) is an orthonormal function meeting the following conditions:

€ Ψ j,k
- € Ψ j,k (t) is a function of the L 2 point space,
-The associated Fourier transform (t) has "physically" compact support, e.g. :

€ ˆ Ψ j,k (υ) has compact support, - € Ψ j,k
€ 20log Ψ j,k (t) max Ψ j,k (t) ( ) ≤ -A dB ( 13
)
where A is a coefficient which can be freely chosen by the user.

The wavelet € Ψ (t) is the so-called "mother" wavelet, which was previously analyzed by Y. Meyer [START_REF] Meyer | Orthonormal wavelets[END_REF] and S. Jaffard [START_REF] Jaffard | Wavelets[END_REF]. Its spectral modulus reaches a maximum at

€ 2 3 F e ( ) ,
where F e (

€ F e = 1/Δt ) is the sampling frequency of the signal, and has a bandwidth equal to 2 octaves 4 3 F e ( ) and a -3 dB bandwidth equal to 1 octave

€ 1 3 F e ,
€ 1 2 F e , F e (
) . In the time domain, it is real and lower than 10 -3 except for an interval 13

€ Δt in length centered on ½ € Δt , which is symmetric with respect to ½ € Δt .
Each wavelet

€ Ψ j,k (t) is focused on € t j 0 = 2 -j t -1 2 Δt
( ) in the time domain and

€ υ j 0 = 2 j+1
3 F e in the frequency domain (with the frequency limits of

€ 2 j 3 F e , 2 j+2 3 F e [ ] )
. The resulting analysis is performed in terms of the time and the frequency. Let the signal s(t) be sampled at frequency F e at N points;

the only non-zero scale lines for s(t) are the lines j ≤ 0, and only dilatations of the "mother" wavelet are performed. In addition, if the sampling frequency is chosen so that the maximum frequency of the signal is

€ F e 3
, the line j = 0 is also equal to zero. The analysis therefore proceeds only up to j < 0. Fig. 4 show an example of a wavelet in the case of j = -5 and k = 0. In this case, the spectral modulus of the wavelet is centered on F c = 0.83 MHz and the bandwidth ranges from 0.41 to 1.66 MHz, which corresponds to the bandwidth of the 1-MHz transducers used.

The wavelet coefficients are defined as projections of s(t) onto the orthogonal basis; they are identical to the cross-correlation between the signal s(t) and the centered wavelet at scale j taken at points € t j,k = k2 -j Δt :

€ C j,k = Ψ j,k (t) s(t) dt = s ⊗ Ψ j,o # [ ] ( k 2 j Δt) -∞ ∞ ∫ = X j t j,k ( ) (14)
where

€ Ψ # (t) = Ψ (-t) (15)
and

€ X j t ( ) = s ⊗ Ψ j,o # ( ) t ( ) (16)
Let s(t) be equal to

€ Ψ J,0 t
( ), which is the wavelet centered on the dilatation scale denoted J (J ∈ Z).

Since orthogonal wavelet decomposition of this particular signal can be performed on the dyadic grid to calculate the cross-correlation X j (t), then:

€ X j (t) = Ψ J,0 ⊗ Ψ j,0 # ( ) (t) (17)
The properties of the orthogonal wavelet decomposition procedure are such that the coefficients, € C j,k are cancelled everywhere except for j = J and k = 0:

€ C j,k = X j t j,k ( ) = 1 if j = J and k = 0 0 elsewhere " # $ (18)
Figs. 5 and6 give an example of the orthogonal wavelet decomposition in the case where J = -5, j ranged from -6 to -4, and the sampling frequency was Fe = 40 MHz. The most interesting term here is

€ X -5 t ( ) = Ψ -5,0 (t) ⊗ Ψ -5,0 # (t)
, which is an autocorrelation term. There is therefore only one valid wavelet coefficient left on the dyadic grid originating from the term C -5,0 (which obviously corresponds to the maximum value of the autocorrelation function).

2.3.3.Algorithm of the TOF estimate

We now take the case of a parallelepipedic rectangle. The response

€ h M (t)
of the sample can be modeled as the sum of two pulse signals, which are comparable to Dirac delta functions (in terms of the distribution):

€ h M t ( ) = A 1 δ t -t 1 ( ) + A 2 δ t -t 2 ( ) (19)
where A 1 (respectively A 2 ) is the amplitude of the first arriving wave located at the TOF t 1

(respectively, of the second arriving wave located at the TOF t 2 ).

Using Eq. ( 1), the output signal s(t) is given by:

€ s t ( ) = x t ( ) ⊗ A 1 δ t -t 1 ( ) + A 2 δ t -t 2 ( ) [ ] (20)
If the input signal x(t) is a wavelet centered on the dilatation scale denoted J (J ∈ Z)

€ x t ( ) = Ψ J,0 (t) (
) , the output signal is given by:

€ s t ( ) = A 1 Ψ J,0 t -t 1 ( ) + A 2 Ψ J,0 t -t 2 ( ) [ ] (21)
The orthogonal wavelet decomposition of this signal on the dyadic grid at scale j is given by:

€ X j (t) [ ] s = A 1 Ψ J,0 t -t 1 ( ) + A 2 Ψ J,0 t -t 2 ( ) [ ] ⊗ Ψ j,0 # (t) = A 1 X j t -t 1 ( ) + A 2 X j t -t 2 ( ) (22)
This equation is the sum of the two cross-correlations products € X j (t) located at TOF t 1 and TOF t 2 , which are attenuated by A 1 and A 2 .

Due to the orthonormal properties of the wavelet decomposition process, TOF t 1 and TOF t 2 are therefore obtained from the coefficients on scale J:

-TOF t 1 is measured by correlating

€ X j (t)
[ ] s with an analyzing pattern P J (t) given by:

€ P J t ( ) = X J t ( ) (23)
-TOF t 2 is measured by correlating the same pattern P J (t) with the set of points belonging to X J (tt 2 ). This set is obtained by using the orthogonal property of the autocorrelation term X J (t), see Eq.( 17).

2.3.4.Numerical example

Here we focus on the input signal x(t) associated with the wavelet (Perspex [START_REF]2.4.1 The speed and attenuation of sound, Kaye & Laby Online[END_REF]). In this case, the wavelength in the water ( λ = 1.74 mm at v 0 = 1485 m.s -1 ) is larger than the thickness of the sample, and it may therefore be difficult to discriminate between the first and second TOF. A 1 is taken to be equal to 1, and A 2 to 0.08, which corresponds to the attenuation of the second arriving wave. The distance (d) between the transducers is taken to be 20 cm, corresponding to a initial TOF t 0 of 134.671 µsec. The attenuation was not based on our experiments, but was chosen to test the validity of the method. . It is impossible to distinguish visually between the two autocorrelation terms X -5 (t-t 1 ) and X -5 (t-t 2 ), and thus between the two TOF. Let us consider only the first term X -5 (t-t 1 ) for the first arriving wave (Fig. 7 middle). Note that since this term is very similar to

€ X -5 (t)
[ ] s , TOF t 1 can be easily obtained by correlating

€ X -5 (t)
[ ] s with an analyzing pattern P -5 (t), as in section 2.3.3, see Eq.( 23). The dyadic grid C -5,k , centered at TOF t 1 is then calculated and the amplitude A 1 is determined at C -5,0 . We now add the second term X -5 (t-t 2 ) (Fig. 7 bottom) centered at TOF t 2 . The overall

€ X -5 (t)
[ ] s is also plotted, for the sake of comparison. It can now be clearly seen that at the former "zeros" on the dyadic grid C -5,k (Fig. 7 middle), we can directly obtain values of X -5 (t-t 2 ) (after eliminating the term C -5,0 , of course): these values constitute the set of points described in section 2.3.3. Correlating € P -5 (t) with this set of points gives an estimate of the position of TOF t 2 . Fig. 8 shows the positions of TOF t 1 and TOF t 2 obtained by performing WBP on the output signal s(t). From these positions, an approximation of the thickness of 1.12 mm (error: 0.63%) and the velocity of 2722 m.s -1 (error: 0.81%) are calculated. This means that in the case of this example, WBP can be used to assess the dimensions of objects, the thickness of which is smaller than the wavelength of the system, with an error lower than 1%.

Now the main problem is how to digitize an initial electric signal e(t), which has to be filtered by the response of the transducers € h T * (t) in order to obtain € x(t) = Ψ J,0 (t) (see Eq. ( 3)).

2.3.5.Generation of a transmitted signal in the form of a wavelet

To digitize the input signal x(t) in the form of wavelets € Ψ J,0 (t), several approaches have been previously tested with various levels of success, such as those based on a simulated annealing method [START_REF] Conil | Non-linear synthesis of input signals in ultrasonic experimental setups[END_REF], which consist in producing the signal iteratively and correcting the quadratic error at each iteration. The simplest method implemented so far is that consisting in directly performing spectral deconvolution of the wavelet by the output signal digitized, using the electric pulse signal e(t) (see Eq. ( 4)) without any immersed sample

€ s(t)
[ ] init (Fig. 9):

€ e(t) = Ψ J,0 (t) ⊗ ⊗ s(t) [ ] init (24)
where

€ • ⊗ ⊗ •
( ) denotes the operation of deconvolution.

In the frequency domain, this equation corresponds to the division of the Fourier transform of each signal:

€ ˆ e (υ) = ˆ Ψ J,0 (υ) ˆ s (υ) [ ] init (25)
This equation is valid only in the frequency bandwidth of the transducers (e.g., with compact support at -3 dB). The inverse Fourier transform gives the electric signal e(t) conveyed to the transmitter via the waveform generator (Fig. 10).

As shown in Fig. 11, the fit between the experimental and theoretical curves is satisfactory, except in the case of small bandwidths around the nominal frequency, probably due to the low stability of the spectral division. But when performing TOF assessments, the fit is perfectly acceptable and the effects of any slight inaccuracy on the results will be negligible. Since the effective transfer function of the apparatus, including the electro-acoustic and mechanical bias, was taken into account, the sensitivity of the solution to this parameter was fairly low. The processing time for wavelet generation was around 4 to 5 seconds on our test configuration (including signal digitalization).

Results and discussion

Rectangular parallelepiped plates

The first experiments were carried out on nine rectangular parallelepipedic plates measuring 10 cm

x 10 cm in size, placed in the vertical position. The constitutive material was Perspex in the case of plate no. 1, Renshape550 in that of plate no. 2, Polystyrene in that of plate no. 3 and PVC in that of plates nos. 4, 5 and 6 obtained from several sources (Fig. 12). Renshape550 is a cream color polyurethane based material, used to simulate the rigidity and durability of the ABS material. Note that the PVC used was the same at a given thickness ("a" and "b"), but differed from one plate to another because the composition of PVC from different manufacturers or batches can vary. The water tank temperature ranged from Θ = 21°5 to 22°7, and the wave velocity in water ranged from v 0 = 1485.3 to 1490.3 m.s -1 . The thicknesses measured using calipers and the mean velocity values published in the literature [START_REF]2.4.1 The speed and attenuation of sound, Kaye & Laby Online[END_REF] are given in Tables 1 and2.

In connection with the experimental protocols detailed in the section 2, paragraphs 2.2 and 2.3, six series of measurements were performed on each sample using the Pulse-mode method (4 measurements) with an electric pulse signal € x(t) = h T * (t) to obtain control values, and WBP (2 measurements) with a wavelet signal € x(t) = Ψ -5,0 (t) as follows:

-Two transmission series without any samples (for Pulse-mode and WBP method);

-Two reflection series on each side of the sample to obtain the thickness (for Pulse-mode method); thickness and the velocity (for WBP method).

Thicknesses and velocities obtained with each of the rectangular parallelepiped plates using the Pulse-mode and WBP methods are given in Tables 1 and2.

The WBP results given in Table 1 show good agreement with the caliper values in terms of the thickness. The error ranged from 0.2 to 3%. (< 3%) The Pulse-mode method gave greater errors ranging from 3.5 to 47 % (> 3%). It is worth noting the case of plates nos. 3, 6-a and 6-b. The thickness (around 1 mm) was below the size of the wavelength (λ ≈ 1.5 mm in the water). The signal obtained in the case of Polystyrene and its WBP are shown in Fig. 13. From these measurements, it may seem difficult at first sight to determine the TOF values: the front and back signals are merged. However, the WBP algorithm makes it possible to evaluate TOF t 1 and TOF t 2 .

Based on these values, the algorithm calculates the thickness, and simultaneously the velocity. The thickness error in comparison with the caliper values is less than 1%. In this case, the Pulse-mode method failed to give thickness values with an error of less than 3%, and in the case of plate no. 6-a, the error amounted to approximately 47%, due to the false value founded for the thickness with this method (0.53 mm vs. 1 mm with the calipers). The Pulse-mode method is based on a detection of the TOF t 1 and the TOF t 2 on signals in reflected mode (see Eq.( 5)). When the echoes of the first and second interface are merged (case of an equal-wavelength thickness), the measurement of this TOF and the calculation of the thickness are false.

As regards the velocities, it can be seen from Table 2 that the results obtained with the WBP method were similar to the previously published mean values. Except for the case of plate no 6-a, the maximum deviation from the published values [START_REF]2.4.1 The speed and attenuation of sound, Kaye & Laby Online[END_REF] was less than 3.5%. In the case of plate no.

6-a, since the thickness (0.53 mm) calculated with the Pulse-mode method was obviously false, the velocity given in Table 2 (2431 m.s -1 ) was that calculated on the basis of the thickness obtained using calipers (1 mm). This value is similar to that calculated using the WBP method. Both velocities differ from the previously published mean values possibly because the velocity in this kind of PVC is nearer to 2400 m.s -1 than 2300 m.s -1 . The second value obtained with the WBP method in the case of plate 6-b was also similar to 2400 m.s -1 .

Velocities calculated with WBP on PVC plates of the same thickness (and presumably having the same composition) were in the same range, and were comparable to the previously published mean values. The mean velocities obtained were 2294 m.s -1 in the case of plates no. 4 (5 mm), 2373 m.s -1 in that of plates no. 5 (1.9 mm), and 2443 m.s -1 in that of plates no. 6 (1 mm). The control values measured using the Pulse-mode method were much less homogeneous ranged from 2400 m.s -1 to 2700 m.s -1 .

Bone samples

Two bovine and one juvenile human femoral long bones were then studied (Fig. 14). The epiphyses of the bone were cut so as to be able to focus only on a cortical area. The age, weight, and sex of each sample were known but were not taken into account in this study. The bones were dry and cleaned for any surrounding soft tissue. They were set in water at room temperature prior to the experiments. The water tank temperature was Θ = 21°7 and the velocities in water v 0 = 1485.7 m.s - 1 . Cortical bones are generally assumed to be elastic material with three anisotropic directions [START_REF] Lasaygues | Ultrasonic characterization of homogeneous, orthotropic elastic bovine bones[END_REF][START_REF] Yoon | Ultrasonic wave propagation in human cortical bone -I. Theoretical considerations for hexagonal symmetry[END_REF][START_REF] Mehta | Bone Elasticity and Ultrasound Velocity are affected by Subtle Changes in the Organic Matrix[END_REF][START_REF] Yoon | Ultrasonic wave propagation in human cortical bone -II. Measurements of elastic properties and hardness[END_REF]. In the axial direction, the axis was assumed to be parallel to the fibers. In the radial direction, the axis was assumed to be perpendicular to the thickness of the area, and in the tangential direction, the axis was assumed to be perpendicular to the fibers. The samples were obtained by cutting the bones in the axial direction and by removing the marrow from each part. They were held on the bench in the horizontal position so that the ultrasonic beam was collinear with the radial axis relative to the thickness. The thicknesses, measured using calipers, in the positions shown Figs. 14 and15, and the previously published mean velocity values are given in Tables 3 and4.

Two heights were tested on the first bovine sample (B1-1 & B1-2), and only one height on the second bovine sample (B2) and the juvenile human sample (H1). The thicknesses and velocities obtained with each bone sample using both the Pulse-mode method and WBP method are given in Tables 3 and4.

The aim of this experiment was to determine the limits of the WBP method when used on biological specimens, and/or when it is not possible to obtain such perfectly parallel interfaces as with manufactured samples. In the case of bone samples, the modeling procedure described in section 2 was therefore not appropriate for dealing with thickness and velocity assessments because of the non-planar shape of the bones. But the focal area of the transducers used was 2x2 mm with the focus set at 150 mm, and at this distance, the wave front was assumed to be plane, without any secondary lobes. The wave reached the interface perpendicularly, which at this scale, was also assumed to be plane. Under these conditions, the above assumptions can be adopted.

Except for the equal-wavelength juvenile human bone sample (H1), as shown in Table 3 as regards the thicknesses (and in Table 4 as regards the velocities), WBP gave values close to the caliper values (as well as close to the previously published values). Indeed, the thickness errors amounted to less than 3% and the velocity errors, to less than 1%. WBP was even found to measure the thickness of the B1-1 sample more accurately (i.e. the calculated error is minimum) than the Pulsemode method (9.18 mm -error = 3% -using the WBP method vs. 8.3 mm -error ≈ 7% -using the Pulse-mode method, while the caliper value was 8.9 ± 0.2 mm). WBP was also found to measure the wave velocities in the bovine bone samples (error < 1%) more accurately (i.e. the calculated error is minimum) than the Pulse-mode method (error > 3%). These results confirm that WBP provides a useful signal-processing tool for obtaining simultaneous longitudinal wave velocity with an error lower than 1%, and thickness values with an error lower than 3%, on a (quasi-) parallelepipedic plate of elastic material.

In the case of a equal-wavelength juvenile human femur (H1, see Fig. 16), the WBP and Pulsemode methods both gave practically the same thicknesses (≈ 1.5 mm), in the rang of values obtained using calipers (≈ 1.6 ± 2 mm). Velocity values differed, however, between the two methods (by 8%), and were lower than the previously published mean velocity (< 3000 m.s -1 ).

Although TOF were detected in the recorded data, some problems arise when measuring the thickness and the velocity, and it is difficult to reach definite conclusions about the accuracy of the WBP method in the particular case of thin bone samples. But it is also worth noting that the surface scan on the samples showed small irregularities of the surface, such as those which can be seen in the zoomed view in Fig. 14, due to the presence of a small remaining porous trabecular part, which may account for the perturbations of the signal digitized and of the values obtained. One of the main hypotheses is that the sample is taken to be homogeneous. This small heterogeneous layer (± 0.1 mm) illustrates the limits of this hypothesis.

The next step consisted in estimating the effects of a lack of flatness of the inner and outer faces, and those of the heterogeneity of bone on the measurements. For this purpose, a linear sweep of ΔL = 1 mm was carried out at various points on a small surface (L = 16 mm) of both B1 (bovine) and H1 (juvenile human) bone samples. This transducer sweep corresponded to an antenna consisting of several components (transducers with an aperture of 2 mm). Position "0 mm", which was the first position of the transducers, corresponded to the previous 1-D experiments (namely, position B1-2 on bovine bone, and H1 on human bone). The corresponding images shown in Figs. 17 and19 give the sounded area (in mm) on the Y-axis vs. time (µs) on the X-axis. The grey level indicates the amplitude of the signals. Linear scanning was first performed without any sample (Fig. 17). The images show the time-shift of the initial signal due to the non-parallel displacements of the two transducers. Since the shape and the amplitude of the signal remained unchanged (Fig. 18), the WBP algorithm is applicable to any displacement of the transducers.

A similar linear scanning procedure was carried out on the bovine and juvenile human bone (Fig. 19) samples.

The localization of TOF t 1 and TOF t 2 can be detected visually (two separate wave fronts can be seen) in the image of the bovine bone (Fig. 19 top). These wavefronts diverged slightly during the scan, which means that the thickness of the sounded area probably increased linearly from top to bottom. However, the image of the human bone (Fig. 19 bottom) gives no information about these TOFs at first sight, due to the small thickness of the bone in comparison with the wavelength of the ultrasonic wave.

WBP was therefore applied and the results obtained are given in Figs. 20 and21.

As shown in Fig. 20 (left), the thickness of the bovine bone was found to increase from top to bottom, which matched the visual interpretation of the bovine bone image (Fig. 19 top). As regards the velocity (Fig. 20 right), the values recorded were fairly consistent over the whole sample (mean value: 3315 m.s -1 , deviation: ± 15 m.s -1 ), which means that this bone segment probably had quite a homogeneous structure, without many irregularities. They also agree with the values published in the literature ( ≈ 3350 m.s -1 , error: 1%). In this experiment, the thickness and velocity showed different patterns of evolution (increasing thickness, stationary velocity), which shows the ability of the WBP method to deal independently with each parameter at each operation.

The case of human bone is slightly more difficult to interpret: both the thickness and the velocity (Fig. 21) values obtained showed some dispersion, even at neighboring positions (± 0.05 mm in the thickness values and ± 100 m.s -1 in the velocity values). This may be attributable to the porous state of the human bone sample surface (Fig. 14). However, the mean thickness (1.55 ± 0.05 mm) and mean velocity (2619 m.s -1 ) were in line with the measurements made with calipers (1.6 ± 0.2 mm) and with the data previously published in the literature (< 3000 m.s -1 ). The porosity of the surface may have affected the ultrasound wave propagation process so drastically (while the thickness remained more or less unchanged) that WBP could not follow these changes exactly, although the frequency was decreased. It is therefore necessary to improve the accuracy of the WBP method in cases of this kind.

Conclusion

The aim of the present study was to present a new ultrasonic method of measuring the thickness of a solid material and the velocity of the wave. This new method involves the use of a Wavelet-Based Processing method to extract the time and frequency information from the signals received, using a suitable transmitted signal. The main advantage of this method is that it requires only a single acquisition in the pure transmission mode to measure the time-of-flight of the wave propagating through the sample under investigation, and to simultaneously calculate the thickness of the sample and the wave velocity, even in the case of an equal-wavelength material thickness. This method was found to be efficient regardless of some imperfections of the devices used, during a linear scanning for example, and of the transducer responses, which can be unknown. The velocity values obtained were quite in line with those previously published in the literature and with those calculated using a more conventional ultrasonic Pulse-mode method. In the case of the manufactured samples with regular flat interfaces, results were obtained on both the thickness and the velocity with an error lower than 3.5%. With the bone samples, although the mean values were satisfactory (the error was less than 1%), the dispersion was greater, especially in a thin (≈ 1.5 mm) human cortical bone sample with a remaining trabecular part. In the case of this natural biological material having a high acoustical impedance (typically 3 to 5 MRayl) contrast with the surrounding soft tissues (or water, typically 1.5 MRayl), the parametric measurements will also depend in the complexity of the wave propagation processes, which generate complex ultrasonic signals forming several packages with different time and frequency signatures. The "Wavelet-Based Processing" method presented above, which makes a compromise between efficiency (error lower than 3.5%) and simplicity (one single measurement), could therefore provide a particularly useful tool for making real-time measurements (≈ 4 to 5 sec. for initial wavelet generation, then 1 to 2 sec. for each data processing time on personal computer, including signal digitalization). However, the theoretical modeling procedure needs to be improved in order to be able to deal with more complex materials (with nonparallelepipedic or circular shapes) and media (inhomogeneous and porous media), and those with irregular, rough or smooth surfaces. It is also proposed to extend this study in order to examine the attenuation of the ultrasonic wave in the material under investigation and to determine whether any dispersion occurs. It is also planned to carry out a more detailed study taking into account the transversal velocity of the wave propagating along the sample. Lastly, we intend to determine whether ultrasonic tomography imaging can be usefully performed on the cross-section of whole circular cylinders and tubes. The main aim of these measurements would be to simultaneously calculate the wall thickness and the velocity of the wave crossing the object, using the "Wavelet-Based Processing" method. ( ) (gray dashed lines) at j = -4 (-top -), j = -5 (-middle -) and j = -6 (-bottom -). Sampling frequency F e = 40 MHz, n= 512 samples. ( ) = Ψ -5,0 (t) ⊗ Ψ j,0 (t) (gray lines) and dyadic coefficients C j,k (black circles) at j = -4 (-top -), j = -5 (-middle -) and j = -6 (-bottom -).

Figure captions

The case j = -5 (middle) corresponds to the autocorrelation function [ ] s (gray dot line), the autocorrelation term X -5 (t-t 1 ) (gray dashed line) and the autocorrelation term X -5 (t-t 2 ) (black line); at the former "zeros" of C -5,k (k ≠ 0), values of X -5 (t-t 2 ) can be measured (black circles).

( ) = A 1 Ψ -5,0 t -t 1 ( ) + A 2 Ψ -5,0 t -t 2 ( ) [ ] ,

5 (

 5 nominal frequency F c = 0.83 MHz, sampling frequency F e = 40 MHz). Let us take the case of a parallelepipedic rectangle with a theoretical thickness of 1.1 mm (1.113 mm in terms of the simulation parameters), and take the known ultrasonic wave velocity to be equal to 2700 m.s-1

Fig. 7 (X - 5

 75 Fig. 7 (top) gives the simulated function

Fig. 1 .

 1 Fig. 1. Overall view of the ultrasonic measurement set-up: mechanical scanner, 8-axis indexer (Rhonax ® ), Waveform Generator (TTI ® TGA 1241), 1-MHz transducers (Imasonic ® ), Sampling oscilloscope (Agilent ® DSO 5014 A), PC computer controlling the scanner and the oscilloscope.

Fig. 2 .

 2 Fig. 2. (-Top drawing -) Conventional ultrasonic Pulse-mode method: signal paths in the reflection (echo) and transmission modes for determining the wave velocities and the wall thickness of a parallelepipedic sample. (-Bottom drawing -) Initial signal path recorded in the transmission mode to determine the initial velocity without any sample.

Fig. 3 .Fig. 4 .

 34 Fig. 3. (-Top drawing -) One-shot transmission mode using the Wavelet-Based Processing method: signal paths recorded in the pure transmission mode to determine the wave velocities and the wall thickness of a parallelepipedic sample. (-Bottom drawing -) Initial signal path recorded in the transmission mode to determine the initial velocity without any sampleFig. 4. (-Top drawing -) Time graph of the wavelet function

Fig. 5 .

 5 Fig. 5. Comparison between the wavelet function

Fig. 6 .

 6 Fig. 6. Comparison between cross-correlation products

Fig. 7 .

 7 Fig. 7. (-Top drawing -) Time graph of the simulated function

  in the case of a parallelepipedic rectangle with a theoretical thickness of 1.1 mm and a wave velocity of 2700 m.s -1 . (-Middle drawing -) Comparison between the simulated time graph of the autocorrelation term X -5 (t-t 1 ) (gray line) and the corresponding dyadic coefficients (black circles) C -5,k . Only C -5,0 is a non-zero coefficient. (-Bottom drawing -) Comparison between the time graphs of the simulated function € X -5 (t)

Fig. 8 .

 8 Fig. 8. Comparison between the simulated output signal s(t) (gray line) in the case of a parallelepipedic rectangle with a theoretical thickness of 1.1 mm, and the wavelet function

Fig. 9 .

 9 Fig. 9. (-Top Drawing -) Time graph of the response

Fig. 10 .

 10 Fig. 10. (-Top drawing -) Time graph of the electric input signal e(t) digitized to obtain a generated wavelet

Fig. 11 .

 11 Fig. 11. (-Top drawing -) Comparison between time graphs of the experimental input signal

Fig. 12 .

 12 Fig. 12. Rectangular parallelepipedic plates (from left to right): 1-Perspex (5.6 mm), 2-Renshape550 (2 mm), 3-Polystyrene (1 mm), 4-PVC (5 mm), 5-PVC (1.9 mm), 6-PVC (1 mm).

Fig. 13 .

 13 Fig. 13. Comparison between the experimental output signal s(t) (black line) obtained in the case of the 1-mm thick rectangular parallelepipedic Polystyrene-plate and the wavelet function

Fig. 14 .

 14 Fig. 14. Bovine bone samples (left: B1, middle: B2) and juvenile human femur (right: H1). Zoomed view of the scanned region of juvenile bone.

Fig. 15 .

 15 Fig. 15. Positions of the calipers used to measure the thickness of the bovine bone samples.

Fig. 16 .

 16 Fig. 16. Comparison between the experimental output signal s(t) (black line) in the case of the 1.6mm juvenile human femur sample: H1, and the wavelet function

Fig. 17 .

 17 Fig. 17. (-Top drawing -) Linear scanning of the transducers without any sample (bovine bone experiment). (-Bottom drawing -) Linear scanning of the transducers without any sample (juvenile human bone experiment).

Fig. 18 .

 18 Fig. 18. Time shift of the wavelet input signal due to the linear scanning at 5 transducer positions.

Fig. 19 .

 19 Fig. 19. (-Top drawing -) Linear scanning of the transducers with bovine bone sample: two separate wave fronts can be visually distinguished (left: position of TOF t 1 , right: position of TOF t 2 ). (-Bottom drawing -) Linear scanning of the transducers with juvenile human bone: no information seems to be available here about TOF t 1 and TOF t 2 . Position 0 mm in this image corresponds to the case shown in Fig.22.

Fig. 20 .

 20 Fig. 20. (-Left drawing -) Thicknesses of the bovine bone; position 0 mm: 8.72 mm (B1-2). (-Right drawing -) Velocities of the bovine bone; position 0 mm: 3317 m.s -1 (B1-2).

Fig. 21 .

 21 Fig. 21. (-Left drawing -) Thicknesses of the juvenile human bone; position 0 mm: 1.5 mm (H1). (-Right drawing -) Velocities of the juvenile human bone; position 0 mm: 2528 m.s -1 (H1).

Acknowledgments

This study was based on research supported by the French National Research Agency (BioGMID -Biological Growth Medium Integrity Diagnosis using bi-modality tomographies -Program) under Grant n°183692 to USAR-CNRS.


Loic Brancheriau
email: loic.brancheriau@cirad.fr

Ghodrati Ashkan

Philippe Gallet

Patrice Thaunay

Philippe Lasaygues

A Ghodrati

Application of ultrasonic tomography to characterize the mechanical state of standing trees (Picea abies)

come

Introduction

Wood is a biological end product generated during cambial growth over successive years [1]. During the growth process, the cambium cells convert most of photosynthesized products and nutriments into various biopolymers for use in the formation of woody tissues (annual ring and cork formation). The wood can thus be considered as a record of responses of trees exposed to different growth conditions but also to changing environments. The properties of an annual ring describe the morphological structure of the species and quantify the trees response to variations of their environment.

In the particular case of plantations, the species are usually selected for their fast growing but the wood quantity produced does not necessarily involve quality. In the objective to improve breeding programs and genetic selections adapted to new forestry practices, non destructive techniques should be directly employed in plantations. The challenge is thus to use accuracy techniques able to quantify a large number of properties in few minutes; the apparatus should also be easy moved through the plantations.

The use of elastic waves (ultrasound or sound waves) to detect decay in trees has been explored and reported by many researchers [2, [START_REF] Bucur | Techniques for high resolution imaging of wood structure: a review[END_REF][START_REF] Divos | Resolution of stress wave based acoustic tomography[END_REF][START_REF] Wang | Acoustic tomography for decay detection in red oak trees[END_REF]. The concept of detecting decay using this method was based on the observation that wave propagation was sensitive to the presence of degradation in wood. In fact, elastic waves appeared to travel more slowly in deteriorated wood than they do in normal wood [START_REF] Socco | Feasibility of ultrasonic tomography for nondestructive testing of decay on living trees[END_REF]. The first generation of equipment used was two-probe systems that measured the wave TOF in a single path. The capability of a single-path approach for tree-decay detection has proven to be limited because stress-wave velocity across tree stems varies substantially even for intact trees, and a standard reference velocity for data interpretation is not readily available [17].

More recently, tomography techniques that were developed for medical and industrial applications or geophysical prospecting have been evaluated for their applicability in standing trees. Different imaging techniques were reported by Bucur in 2003 [START_REF] Brancheriau | Ultrasonic device for the imaging of green wood[END_REF]. In this review article, the author detailed the possibility to use ionizing radiation (X-ray and gamma ray), thermography, microwaves, ultrasonics and nuclear magnetic resonance. A tomograph using ionizing radiation was designed for standing trees [12]. Commercial devices were also developed and were based on the determination of the time of flight (TOF) of elastic waves induced by an instrumented hammer (acousto-ultrasound tools). These devices were presented in details by Schubert in 2007 [START_REF] Schmoldt | CT imaging, data reduction, and visualization of hardwood logs[END_REF]. The three tools were PiCUS®, ARBOTOM® [START_REF] Nicolotti | Application and comparison of three tomographic techniques for detection of decay in trees[END_REF] and FAKOPP3D® [START_REF] Dikrallah | Study of acoustic wave propagation through the cross section of green wood[END_REF]. In contrast to X-and gamma-ray based tomography, tomographic methods based on electrical resistance, microwaves or elastic waves are not dangerous. Nicololli [START_REF] Maurer | A simple anisotropy correction procedure for acoustic wood tomography[END_REF] compared these three approaches by calculating the tomograms of decayed trunks. These tomographic methods consist of inverting electrical resistance or TOF of the electromagnetic and elastic waves by algorithms from the field of geophysical exploitations. Tomography based on elastic waves showed the best results in terms of depicting the geometry and the location of the inner wood-decay. However, the TOF are inverted by filtered back projection which assumes isotropic material behavior and straight rays (PICUS® and FAKOPP3D®). Straight rays are appropriate in isotropic media, otherwise inversion should account for curved rays [START_REF] Deflorio | Detection of incipient decay in tree stems with sonic tomography after wounding and fungal inoculation[END_REF][START_REF] Schubert | Acousto-ultrasound assessment of inner wood decay in standing trees: possibilities and limitations[END_REF]. [START_REF] Socco | Feasibility of ultrasonic tomography for nondestructive testing of decay on living trees[END_REF] used an iterative TOF inversion method (root mean square optimization method) developed for geophysicists. This inversion method can be used with straight or curved rays [START_REF] Schubert | Acousto-ultrasound assessment of inner wood decay in standing trees: possibilities and limitations[END_REF]. Maurer also proposed in 2006 a correction procedure to take into account the anisotropy of wood. This correction assumed a slight anisotropy and took into account a radial to tangential velocity ratio in a standard isotropic reconstruction procedure [START_REF] Lin | Application of an ultrasonic tomographic technique for detecting defects in standing trees[END_REF].

With the challenge to use an accurate and fast nondestructive technique in plantations, the objectives of this study were to design an automatic tomograph and to compare the effect of the presence of bark and drilled holes on the ultrasonic parameters: slowness and attenuation. At the beginning of the development of the tomograph, two ultrasonic probes of 300 kHz were used [START_REF] Beall | Overview of the use of ultrasonic technologies in research on wood properties[END_REF]. The probes consisted of a wheel in which the emitter was placed and the coupling medium was made by an elastomer surrounding the wheel. After several tests on green and dry woods in the transverse and longitudinal directions, these probes were replaced because the transmission depth at 300 kHz was between 3 cm to 10 cm depending on the test configuration (without the presence of bark). The anisotropy of wood was not taken into account and a classic filtered back projection algorithm was used to compute the tomography maps.

Experimental methodology

Materials

Two logs of spruce (Picea abies) were harvested from different trees near the city of Nowshahr (Iran). The sampling height was between 0.80 m and 1.30 m. The associated mean diameter was between 15.0 cm and 16.5 cm. The section shape of the logs was almost circular and the logs had also a straight vertical aspect. The trees were 26 years old. The sampling was made in September 2010 and shipped via France to be tested in October; the samples were wrapped into plastic films to minimize the effect of wood drying.

Experimental protocol

The first log was cut in two disks of 5 cm thick and one disk of 15 cm thick. These disks were sampled in the middle area of the log to avoid cracks due to wood drying and to have high moisture content. The disks of 5 cm thick were used to determine the mean moisture content and also the moisture content profile along the radius (Figure 1). The moisture content was determined before and after the ultrasonic tests. The samples of Figure 1 were weighted then dried in an oven (103°C) during 48 hours then weighted again. The Equation 1 was used to compute the moisture content. The disk of 15 cm thick was tested in tomography with 3 states of the sample: bark, debarked and drilled (2 holes of diameter 1 cm and 2 cm, see Figure 2). The tomography tests were made with two repetitions and an angular step of 30° to estimate the measurement uncertainties and with an angular step of 10° for tomographic imaging. The second log was cut in two disks of thickness 15 cm. These disks were sampled in the middle area of the log. In each disk a hole of diameter 5 cm was drilled, one in the middle and the second offcentered (Figure 3). These two disks were tested in tomography with an angular step of 10°.

Ultrasonic tomography

The ultrasonic tomograph consisted in an aluminium ring supported by a tripod (Figure 4) [START_REF] Beall | Overview of the use of ultrasonic technologies in research on wood properties[END_REF]. To be able to test small disks, a specific support was added and mounted between the tripod (Figure 3). This device was designed to test standing trees. Thus, the height of the ring can be adjusted from 1 m to 1.6 m but the typical height is 1.3 m. If the apparatus is mounted on a trunk (Figure 4) the two parts of the ring are first assembled, and then the tripod is fixed. The apparatus is then placed against the trunk using two narrowing elements. Each probe is mounted on a trolley and these trolleys are slipped within the ring. Finally, the last part of the ring is fixed and the third narrowing element is applied. The probes had a frequency of 55 kHz and were designed for acoustic emission purpose. These probes were mounted on a specific support with a wheel and grease was used as coupling medium. The emission consisted in 5 square periods of 480 V. The output analog signal was amplified of 80 dB. The acquisition was performed by a converter of 16 bits resolution and a sampling frequency of 5

MHz. The duration of acquisition was set to 410 µs. The output signal was previously filtered by a Morlet wavelet at 55 kHz with a bandwidth of 20 kHz for denoising (Figure 5). The ultrasonic parameters used for the image reconstructions were the slowness (s/m) and the attenuation (dB/m at 55 kHz). The slowness Sl was computed using cross correlation between the input and the output signals (Equation 2). The ratio between received energy and emitted energy was used to determine the attenuation At (Equation 3). The tomographic images were computed using an angular step of 10° which led to 900 acquisitions (36 equidistant acquisition points). A complete scan was done in 40 minutes with the apparatus. In comparison, three well trained operators were asked to make their own scan with 18 measurement points (306 acquisitions) and the duration of the scanning process was never less than one hour.

For a projection, the emitter is motionless while the receiver moves around the ring. The next projection is initiated by an angular step motion of the emitter followed by another motion of the receiver all around the ring (an acquisition is done at each angular step). The reconstruction algorithm was based on the inverse Radon transform in fan beam geometry (linear interpolation) with a modified Shepp-Logan filter. Considering the medium as isotropic, the elastic waves were assumed to propagate in a straight line. The velocity (and therefore the attenuation) of waves was then obtained from measurements of travel-time (of energy) differences by inverting the Radon transform.

The theoretical upper limit of the number of computed pixels was 450 (mesh of 39 mm², Equation 4). The physical resolution was however about 26 mm (axial resolution: wavelength) and 31 mm (lateral resolution: Fresnel zone, Equation 5, Figure 6). This resolution led to a mesh of 806 mm²; the images were thus 2D smoothed by a blackman filter of the wavelength size (the errors induced by the use of the modified Shepp-Logan filter was also highly reduced). It was important to highlight that the redundancy of the information (39 mm² compared to 806 mm²) allowed to decrease the effect of erroneous experimental values in the imaging process.

=

With R F : the Fresnel zone radius, S d : the scanned depth (the rotation radius of the emitter or receiver transducer), C: the wave celerity and f: the frequency of the transducer.

Results and discussion

In the next paragraphs, the slowness was presented in 10 3 (s/m) or its inverse: the speed in (m/s). The attenuation in (dB/m) was normalized by the diameter and thus presented in (dB).

Moisture content determination

The mean moisture content was 23% before the tests and 21% at the end of the experiment. The moisture content profile along the radius was presented at the Figure 7. The difference in moisture content from the pith to the bark was 5% before test and 10% after test. Despite the presence of the plastic film on the logs and on the disks, the wood samples tested were not in a green state and the natural drying effect increased while testing (phenomenon important after debarking). Small cracks appeared at the periphery of the disks due to drying (Figure 2, right image).

Effect of bark and holes on the ultrasonic parameters

Mean comparisons were made to study the effect of bark and holes on the ultrasonic parameters (Figure 9). The number of acquisitions was 108 with an angular step of 30°. The statistical tests were firstly a Student's t-test on paired samples (Table 1). 1, were biased by the non normal distribution of the parameters. To highlight this phenomenon, two histograms were plotted (Figure 10) after the tomographic tests with an angular step of 10° (N=900). The bimodal distribution was shown on this figure for speed with two modes between 1200 (m/s). The same type of distribution was also observed concerning the attenuation between -51 (dB). Equivalent observations were made for the debarked and drilled disks. The results of Table 1 were secondly completed by non-parametric Wilcoxon tests (valid when the normality assumption was not satisfied or the sample size was too small). The Wilcoxon tests agreed with those of Student't for attenuation but did not confirm the Student't results for speed (Table 1). The mean attenuation measured on the disk with bark was thus significantly different with the one after debarking. The mean attenuation was higher with bark than without bark. No significant difference was found between the mean attenuation of debarked disk and of drilled disk. Concerning the speed measurements, differences existed between disks with bark, debarked and drilled but the probability values were too high to be absolutely affirmative (Figure 9 and Table 1).

Relationship between slowness and attenuation

The Figure 11 showed a scatter plot between slowness and attenuation for the disk with bark. A significant linear relationship was found with a R² equal to 0.51 (F=933, p < 0.001, N=900). The attenuation increased when the slowness increased or the speed decreased. However, the degree of fit should be read with care because the Figure 11 clearly showed two populations (associated with the distributions of the Figure 10). Outliers were also visible due to the uncertainty on the slowness determination; but they were not associated to particular positions of the probes.

Tomographic images

The computed maps of slowness and attenuation for the first log with and without the bark were presented at Figure 12. The images were displayed in 87x87 pixels and based on a 10° angular step (900 acquisitions). There was no visible effect of the presence of the bark on the outlines of the maps (Figure 12). This phenomenon was due to the low thickness of the bark (3 to 5 mm) compared to the wavelength (26 mm) added with the smoothing effect of the 2D filter. The Figure 12 showed however two different areas from the pith to the periphery. These areas were highlighted at Figure 13 where the paths were drawn if the slowness was lower than 1.2x10 -3 s/m and if the attenuation was upper than -51 dB. The periphery was thus characterized by a higher speed and a lower damping than the ones of the pith.The inverse phenomenon was observed by Socco and interpreted as an effect of the wood anisotropy [START_REF] Schubert | Acousto-ultrasound assessment of inner wood decay in standing trees: possibilities and limitations[END_REF]. The shorter paths (periphery) were more influenced by the tangential behavior of wood, while the radial behavior had a more important effect for longer paths (center). However, the radial velocity of wood is known to be upper than the tangential velocity. In this specific case (Figure 12), this observation was explained by the transition from juvenile to mature wood coupled with the effect of the moisture distribution (Figure 7). The tomographic images obtained for the drilled disks of the two logs were presented at Figure 14 and Figure 15. The effect of the holes was visible on the maps by areas with a high slowness and also a high damping. This observation was valid even for the hole which diameter (10 mm) was the half of the wavelength (26 mm). It seemed that the position of holes in the disk had an effect on the reconstructed pattern; the holes closed to the pith were characterized by two elliptic patterns at the difference of those closed to the periphery. The shape of the patterns associated with the holes was not a circle but an ellipse and several hypotheses were advanced: the anisotropy was not taken into account, some cracks due to natural drying occurred, the reconstruction procedure was not well adapted for missing data (the dimension of the probe supports not allowed to scan the total angular area but only 250° for each projection), the mounting of the probes (Figure 5, right image) induced back waves which altered the received signals. These hypotheses will be investigated in a future study to improve the tomographic process.

Conclusions

An automatic tomograph device specifically designed for standing trees was presented. The scanning process employed constituted an improvement compared to previous works. This improvement was a fast process with only two probes leading to scans with an angular precision of 10 degrees without the intervention of an operator during the process. This tomograph was used on two logs of spruce (Pica abies) to study the effect of the bark and also the effect of defects like holes on the reconstructed images. The tests were not done at a green state with mean moisture content of 22% and a gradient of 7% between the pith and the periphery. The ultrasonic parameters used for imaging of wood were slowness and attenuation. The use of attenuation was not reported in previous works and a better accuracy of determination was found compared to slowness. It appeared that time of flight determination was less accurate even when cross correlation was employed. This phenomenon was due to the dispersion and the dissipation of wood material which caused a high attenuation coupled with modifications of the incident waves. The means of the ultrasonic parameters were not sufficient to discriminate the presence of defects in trees, but the presence of the bark had a significant contribution in the signal attenuations. For the tested samples, it was shown that speed near the bark was higher than in the pith because of the presence of juvenile wood coupled with the moisture gradient (natural drying occurred during the tests). In the same way, the damping near the bark was lower than in the center. A significant relationship between slowness and attenuation was shown but the population was bimodal due to the presence of juvenile wood. It was however possible to observe that a high speed was coherent with a low damping. The presence of the bark was not visible on the tomographic images. This was probably due to its small thickness compared to the wavelength. The presence of the holes was however clearly visible but artifacts did not permit a precise shape reconstruction and localization. Removing these artifacts will be the objective of future developments.

Acknowledgments

The BIOGMID project was supported by a grant from the French National Research Agency (BLAN07-1_183692, "BioGMID").

  With MC: the moisture content, M H : the green weight and M 0 : the oven-dry weight.

Figure 1 .

 1 Figure 1. Sampling pattern for moisture content determination (dimensions of the small samples: 2x3x5 cm 3 ).

Figure 2 .

 2 Figure 2. Tomography tests on the first log (from left to right: bark, debarked and drilled with two holes of diameter 1 cm and 2 cm).

Figure 3 .

 3 Figure 3. Tomography tests on the two disks of the second log (hole diameter of 5 cm).

Figure 4 .

 4 Figure 4. Mounting steps of the ultrasonic tomograph.

Figure 5 .

 5 Figure 5. Ultrasonic signal filtered by a Morlet wavelet at 55 kHz with a bandwidth of 20 kHz (a) and associated Fourier transform (b).

  With Sl: the slowness, At: the attenuation, d: the distance between the emitter and the receiver, e(t): the input signal and s(t) the output signal.

Figure 6 .

 6 Figure 6. Lateral resolution and wavelength as a function of the frequency and according to the sample diameter (the acoustic wave speed was set to 1400 m/s).

  : the number of pixels, N Projections : the number of projections (36), N Sample : the number of acquisitions per projection (25).

Figure 7 .

 7 Figure 7. Moisture content (%) according to the radial position before and after the ultrasonic tests (samples from the first log).

Figure 8 .

 8 Figure 8. Relationship between the two repetitions for (a) speed and (b) attenuation (disk with bark, N=108).

Figure 9 .

 9 Figure 9. Means and standard deviations of speed (a) and attenuation (b) computed on 108 acquisitions for the disk with bark, debarked and drilled.

Figure 10 .

 10 Figure 10. Histograms of speed (a) and attenuation (b) of acoustic waves at 55 kHz computed for the disk with bark (N=900).

Figure 11 .

 11 Figure 11. Relationship between slowness and attenuation for the disk with bark (N=900).

Figure 12 .

 12 Figure 12. Tomographic maps of slowness (10 -3 s/m) (a, b) and attenuation (dB) (c, d) for the disk of the first log, with bark (a, c) and debarked (b, d).

Figure 13 .Figure 14 .

 1314 Figure 13. Paths of emitter to receiver with a slowness lower than 1.2x10 -3 s/m (a) and an attenuation upper than -51 dB (b) for the disk of the first log with bark.

Figure 15 .

 15 Figure 15. Tomographic maps of slowness (10 -3 s/m) (b, e) and attenuation (dB) (c, f) for the two disks of the second log (a, b).

Table 1 .

 1 t-Test and Wilcoxon test on paired samples for speed and attenuation between the disk with bark, debarked and drilled.

				p-value	Wilcoxon
	Data sets	t	ddl	(bilateral)	p (bilateral)
	Bark -Debarked
		Speed -1.934 107	0.056	0.496
	Attenuation -5.888 107	< 0.001	< 0.001
	Drilled -Debarked
		Speed 2.196	107	0.030	0.126
	Attenuation -0.059 107	0.953	0.998

The t-test results, presented in Table