Ultrasound computed tomography on standing trees: accounting for wood anisotropy permits a more accurate detection of defects

Abstract

Key message

Considering anisotropy in image reconstruction algorithm for ultrasound computed tomography of trees resulted in a more accurate detection of defects compared to common approaches used.

Context

Ultrasound computed tomography is a suitable tool for nondestructive evaluation of standing trees. Until now, to simplify the image reconstruction process, the transverse cross-section of trees has been considered as quasi-isotropic and therefore limiting the defect identification capability.

Aims

An approach to solve the inverse problem for tree imaging is presented, using an ultrasound-based method (travel-time computed tomography) suited to the anisotropy of wood material and validated experimentally.

Methods

The proposed iterative method focused on finding a polynomial approximation of the slowness in each pixel of the image depending on the angle of propagation, modifying the curved trajectories by means of a raytracing method. This method allowed a mapping of specific elastic constants using nonlinear regression. Experimental validation was performed using sections of green wood from a pine tree (Pinus pinea L.), with configurations that include a healthy case, a centered, and an off-centered defect.

Results

Images obtained using the proposed method led to a more accurate location of the defects compared to the filtered backprojection algorithm (isotropic hypothesis), considered as reference.

Conclusion

The performed experiments demonstrated that considering the wood anisotropy in the imaging process led to a better defect detection compared to the use of a common imaging technique.

1 Introduction

Due to its biological origin, the properties of wood have a high degree of variability. Trees are indeed complex lignified plants with a morphology and internal structure marked by a high heterogeneity that is structured by the genetic origin of the trees but also by their conditions of management and environmental factors, including the climate. Characterizing wood properties and inner heterogeneities of trees brings information useful for clonal selection, grading of logs, but also for tree risk management in urban areas. Standing trees in urban areas present an important ecological and health role (Pokorny, 2003). The tree population requires control and management processes to assess potential risks. For example, falling trees pose a danger to pedestrians, houses, and cars. Aiming to minimize the risk associated with tree failure, several significant advances in decay detection equipment, formulas, and guidelines for assessing hazardous trees have been done (Johnstone et al. 2010; Leong et al. 2012).

The most common methods to evaluate the inner state of standing trees consist of drilling measures or in acoustical or electrical spot measurements (Pellerin and Ross 2002). However, these techniques remain limited, providing only one-dimensional information. To overcome this limitation, nondestructive acoustic and ultrasonic imaging methods have been presented as an alternative to carrying out the assessment of the inner structures of trees, without altering their condition (Bucur 2003; Arciniegas et al. 2014a).

The most known commercial devices use acoustic waves, where a hammer is used as emitter. Different studies have shown that these acoustic devices present different drawbacks, such as a low spatial resolution, images difficult to interpret, detection precision not always optimal, and a relatively long time of execution (Rabe et al. 2004; Gilbert and Smiley 2004; Wang et al. 2007; Deflorio et al. 2008).

The quality of the obtained image (time-of-flight tomography) is influenced by different aspects, resulting from the interaction between transmitted waves and heterogeneous and anisotropic media. This includes the frequency (sound or ultrasound), the signal to noise ratio, the phenomena of attenuation and diffraction, the total number of sensors, the modeling and the assumption of the propagation, and the image reconstruction algorithm (Arciniegas et al. 2014b). The ultrasonic approach aims to increase the frequency of excitation to obtain images with higher resolution.

The basic principle for ultrasound computed tomography (USCT) is that decay inside wood affects the propagation of elastic waves, leading to a decreasing velocity and an increasing attenuation. The Fermat’s principle states that the first-arriving wave that of the compressional wave will tend to travel along the fastest path; therefore, the presence of pronounced compressional wave velocity contrasts will tend to curve the rays (Maurer et al. 2005). Decay regions will be avoided, as they slow down waves, resulting in larger time-of-flight (TOF) measurements compared to healthy wood.

USCT of wood was first presented by Tomikawa et al. (1986), aiming to perform non-destructive testing of wooden poles. TOF measurements were obtained considering a fanbeam geometry and the filtered backprojection algorithm for image reconstruction. It was found that the proposed USCT algorithm and device were adequate to detect rotten areas, but some drawbacks were highlighted, such as a weak image resolution and long computation time. Afterward, different approaches have been presented (Martinis et al. 2004; Lin et al. 2008; Brancheriau et al. 2008; Brancheriau et al. 2011). It was pointed out that ultrasonic techniques are suitable for standing trees quality evaluation, using frequencies ranging from 22 to 1 MHz, enabling to detect decay, knots, fungal attack, and other defects. Up to now, these approaches have not considered the anisotropy property of wood in a mechanical modeling of wave propagation, affecting the image reconstruction. To build the tomographic image, it is necessary to know the trajectories (rays) followed by the waves from the transmitter to each receiver. As shown in different studies, wood anisotropy leads to deformed wavefronts, compared to the spherical wavefronts obtained for isotropic media (Payton 2003; Schubert et al. 2008; Gao et al. 2014; Espinosa et al. 2019a). This deformation leads to curved trajectories (Maurer et al. 2006), compared to the straight-line trajectories under isotropic condition. These curved trajectories should be considered for a suitable reconstruction.

This study aimed to consider the anisotropy property of wood in the tomography image reconstruction process, by developing an iterative method that takes into account the curved ray paths. This inversion procedure allowed determining locally the specific elastic constants of the wood material. We proposed an experimental configuration to evaluate the precision of defect detection using the proposed image reconstruction technique, compared to a common method that considers a straight-line hypothesis (filtered backprojection algorithm). The proposed inversion procedure was tested in a pine sample, considering three cases: a healthy case, a simulated centric defect case, and a simulated eccentric defect case. Results will be discussed according to several criteria with the aim of applying the method with a view to in situ operations: the computing time with a common computer, the number of probes and sensors, and the possibility to implement the proposed method in the existing commercial devices.

2 Materials and methods

2.1 Wood samples

From a healthy trunk of a pine tree (Pinus pinea L.), 3 wood disks were obtained. The disk diameter was 24 cm, and the thickness of each disk was 3 cm. The tree age was 55 years. To reduce water loss during experiments, logs were sealed and stored in a temperature-controlled room (+ 4 °C). The average moisture content was 92% (moisture reduction during the tests was less than 3%).

To simulate the presence of a defect, a circular hole was drilled with a diameter of 7.6 cm. Two defect positions were tested: centered and off-centered (halfway between the center and the bark). Even when the tree section used for testing did not present any decay, there was a presence of reaction wood and juvenile wood that affects the physical and mechanical properties (Timell 1986; Ross 2010; Gardiner et al. 2014), and consequently the propagation of waves (Bucur 2006; Brancheriau et al. 2012a).

2.2 Ultrasonic testing

Ultrasonic measurements were obtained at 16 different positions around the wood disks, distributed uniformly. To validate experimentally the image reconstruction procedure, it was decided to debark the disks in order to improve the coupling of the sensors (to increase the energy transfer and the quality of the signals); this avoided bias not due to the inversion procedure itself. The ultrasonic chain of measurement is shown in Fig. 1. The ultrasonic (US) sensors had a main resonant frequency at 60 kHz (Physical Acoustics Corporation R6α). To improve energy transfer, a fluid couplant was used. The electro-acoustical set-up consisted of a signal generator and an oscilloscope (Picoscope 2000), connected to a PC for data acquisition. Additionally, two electronic amplifiers were used, one located after the signal generator (Single Channel High Voltage Linear Amplifier A800, FLC Electronics) to increase the energy sent to the US transmitter, and other at the US receiver output (Physical Acoustics Corporation AE2A/AE5A), both with an amplification of 40 dB. The excitation signal was a chirp signal (central frequency of 60 kHz, bandwidth of 48 kHz), leading to a concentrated power spectrum around the sensor central frequency that allows a time-of-flight (TOF) measurement using a cross-correlation (Espinosa et al. 2018b). The full set of TOF measurements (i.e., projections) was obtained by placing the transmitter at the 16 different points along the wood disk circumference. For every transmitter position, the receiver position changed in the remaining 15 points.

Fig. 1

Ultrasonic chain of measurement, including the ultrasonic pair of sensors, electrical signal generator and oscilloscope, with input and output amplifiers

2.3 Inversion procedure

It is necessary to establish the mathematical model allowing to express the relation between the ultrasonic measurements and the inner mechanical wood parameters. The inverse problem consists in finding the model parameters from a set of measurements.

Considering the transversal section of a trunk (2D case), wood mechanical properties are determined by two perpendicular axes: the radial axis (R) aligned from the bark to the pith and the tangential axis (T), tangent to the growth rings. The stress (σ) and strain (ε) relations are given by Hook’s law, with the rigidity matrix C for an orthotropic material. For the RT plane, the elastic parameters are Young’s moduli E_R (radial direction) and E_T (tangential direction), the shear modulus G_RT, and Poisson’s coefficient ν_RT.

As time-of-flight (first arrival) tomography was considered, therefore only pure compression waves were considered (hereinafter refers to as ultrasonic waves), omitting all other second-order phenomena (refraction on the defect, shear waves regardless of the incidence conditions, mode conversion, amplitude modification, and dispersion). The ultrasonic wave velocity V in wood depends on the direction of propagation θ, the elastic parameters, and the wood density as expressed by the Christoffel equation (Espinosa et al. 2018):

$$ {\displaystyle \begin{array}{c}{\Gamma}_{RR}={C}_{11}{\cos}^2\theta +{C}_{66}{\sin}^2\theta \\ {}\begin{array}{c}{\Gamma}_{TT}={C}_{66}{\cos}^2\theta +{C}_{22}{\sin}^2\theta \\ {}\begin{array}{c}{\Gamma}_{RT}=\left({C}_{12}+{C}_{66}\right)\cos \theta \sin \theta \\ {}V=\sqrt{\frac{\Gamma_{RR}+{\Gamma}_{TT}+\sqrt{{\left({\Gamma}_{TT}-{\Gamma}_{RR}\right)}^2+4\cdotp {\Gamma}_{RT}^2}}{2\rho }}\end{array}\end{array}\end{array}} $$

(1)

The θ angle is formed between the direction vector n (normal to the wavefront) and the radial (R) direction (Fig. 2). Then, higher velocities are obtained for waves propagating in the radial direction (θ = 0°) compared to the tangential direction (θ = 90°). Therefore, wavefronts are not strictly circular as for the case of isotropic materials, resulting in curved trajectories between the transmitter and each receiver (Espinosa et al. 2019a).

Fig. 2

Geometrical definition of θ angle and of vector n normal to the wavefront

For image reconstruction, it is necessary to know the real trajectories (or ultrasonic ray paths). This process considers a space divided into cells (or pixels), that are traversed by the rays. The objective is to estimate a local slowness α for every cell. The addition of individual slowness values along the ray leads to the TOF measurement at the receiver (Eq. 2). For a given pair transmitter-receiver, which trajectory m crosses k cells, the time-of-flight t_m can be written as:

$$ {t}_m={\sum}_{k\ \mathrm{along}\ m}{l}_{mk}{\alpha}_k $$

(2)

with l_mk as the length of the ray segment. Solving all the equations for α results in the reconstruction of the inner parameter (inverse problem). Two considerations arise at this point: first, the matrix formulation should be adapted to consider the slowness dependency on the angle of propagation, and second, the estimation of the curved trajectories, a priori unknown, affected by the wood anisotropy and the presence of defects (slow propagation regions).

In many algebraic reconstruction techniques, the l_mk terms are replaced by 1’s and 0’s, depending on whether or not the image cell is traversed by the m ray (more details in Kak and Slaney 2001). The total length of a ray L_m is then considered as the sum of individual ray segments l_mk, assuming these segments to be uniform. We have L_m = D_ml_m, with D_m as the total number of ray segments in the path m. Dividing both sides of Eq. 2 by L_m we have:

$$ {A}_m=\frac{t_m}{L_m}=\frac{1}{{\mathrm{D}}_m}{\sum}_{k\ \mathrm{along}\ m}{\alpha}_k $$

(3)

with A_m as the total slowness of the ray. So, for every pair of transmitter-receiver, an equation can be formulated, resulting in a system of linear equations. To adapt the matrix formulation to consider the slowness dependency on the angle of propagation, the Christoffel equation was approximated by a polynomial function of the angle of propagation (5th degree). This polynomial was selected from an R² evaluation using different degrees, allowing the determination of the specific elastic parameters of the Christoffel equation (Espinosa et al. 2020). Then, the slowness for a pixel k crossed by a ray m was written as:

$$ {\alpha}_k={\beta}_{5,k}{\theta}_{k,m}^5+{\beta}_{4,k}{\theta}_{k,m}^4+{\beta}_{3,k}{\theta}_{k,m}^3+{\beta}_{2,k}{\theta}_{k,m}^2+{\beta}_{1,k}{\theta}_{k,m}^1+{\beta}_{0,k}{\theta}_{k,m}^0, $$

(4)

The inverse problem solution corresponded to find the polynomial coefficients β for all the pixels. Combining Eqs. 3 and 4 for every pixel, we obtain a set of linear equations, as presented in Eq. 5. The SIRT (Simultaneous Iterative Reconstruction Technique) method was used to solve this matrix problem (Kak and Slaney, 2001). To increase the number of equations, an interpolation of the TOF measurements was done. Virtual sensors were located between the original ones, passing from 16 sensors to 32 sensors, by linear interpolation of the TOF measurements and the positions of the sensors. The choice of an interpolation to 32 sensors was linked to a previous research showing that reconstruction algorithms such as SIRT (used in the proposed method) converged from 30 transducer positions (Arciniegas et al. 2014b). The sinogram interpolation has shown to be an alternative to increase the defect identification capacity when using a limited number of sensors (Espinosa et al. 2019c).

$$ \left[\begin{array}{cccccccccc}{\theta}_{1,1}^5& \cdots & {\theta}_{1,1}^0& 0& \cdots & 0& \cdots & {\theta}_{N,1}^5& \cdots & {\theta}_{N,1}^0\\ {}0& \cdots & 0& {\theta}_{2,2}^5& \cdots & {\theta}_{2,2}^0& \cdots & 0& \cdots & 0\\ {}\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ {}{\theta}_{1,M}^5& \cdots & {\theta}_{1,M}^0& 0& \cdots & 0& \cdots & {\theta}_{N,M}^5& \cdots & {\theta}_{N,M}^0\end{array}\right]\left[\begin{array}{c}{\beta}_{5,1}\\ {}\vdots \\ {}{\beta}_{0,1}\\ {}{\beta}_{5,2}\\ {}\vdots \\ {}{\beta}_{0,2}\\ {}\vdots \\ {}{\beta}_{5,N}\\ {}\vdots \\ {}{\beta}_{0,N}\end{array}\right]=\left[\begin{array}{c}{D}_1{A}_1\\ {}{D}_2{A}_2\\ {}\vdots \\ {}{D}_m{A}_m\end{array}\right], $$

(5)

For the estimation of the curved trajectories, an initial guess corresponded to straight-line paths. With this initial guess, it was possible to perform a first inversion process, leading to the polynomial coefficients β of the slowness for all the pixels. Using these coefficients, we can perform a numerical simulation to obtain an estimation of the TOF and the trajectories. The model used for the simulation of the forward problem corresponded to a raytracing approximation (Espinosa et al. 2019a, Espinosa et al. 2019b). This process was repeated until the TOF and trajectory difference with respect to the model was minimized.

From the polynomial approximation obtained for every pixel, the values of the specific elastic constants could be estimated using a nonlinear least square regression (Seber and Wild 1989). The result is an image of an elastic parameter of wood. In this case, it was decided to determine the Young’s modulus in the radial direction (E_R) for a given density value (E_T and G_RT can also be computed with the proposed procedure). The proposed inversion schema is shown in Fig. 3.

Fig. 3

Flowchart for the proposed reconstruction method presenting an iterative schema

2.4 Assessment of decay detection

To evaluate the improvement in the image reconstruction by considering the wood anisotropy, a comparison was made with images obtained using an algorithm that considers an isotropic and straight-line rays’ assumptions. The analytical method was based on the sum of the inverse Fourier transforms of the filtered backprojections (so-called FBP method) (Kak and Slaney, 2001). The sinograms used for the FBP method were also interpolated. From each image, a profile is presented, which crosses the center of the disk.

The images for the defective cases were segmented using a threshold. It has been previously shown that using a threshold value to separate healthy and defective regions is not always the most suitable way of interpreting tomograms (Espinosa et al. 2017). In this specific case, the thresholding segmentation was chosen because it allowed comparing quantitatively the two image reconstruction methods. The threshold values corresponded to 30%, 50%, 70%, and 90% of the mean value from the reconstructed image. To assess the defect discrimination capability, two statistical measurements were obtained: the sensitivity or true positive rate (TP) and the fall-out or false-positive FP rate. The first is computed as the ratio of correctly identified pixels inside the defective area; for correct identification, expected values should be as near as possible to 100%. The second corresponds to the ratio of incorrectly identified pixels, i.e., the pixels classified as a defect, that are outside the defective area; for correct identification, expected values should be as near as possible to 0%.

3 Results

3.1 Healthy case

First, a healthy disk was tested. The cross-section is shown in Fig. 4a, drawing attention to the presence of compression wood (right area in the image). After performing the tomographic reconstruction using the proposed method and the reference FBP method, the corresponding E_R parametric image, and the reference velocity image are presented in Fig. 4b and c respectively. For the E_R image, the mean value was 1258 MPa, with a standard deviation of σ = 126 MPa. For the FBP image, the mean velocity was 1551 m/s with σ = 315 m/s. Even when both images presented relatively flat surfaces, the FBP image presented more artifacts near to the borders, where the precision of the TOF measurements is lower (Espinosa et al. 2018b). Two horizontal profiles traversing the center of the image are shown in Fig. 4d and e, for the E_R parametric image and the reference velocity image respectively. The profile obtained in Fig. 4d presents three different parts: in the left, E_R values corresponding to normal wood, then around the pith, lower E_R values that could be associated with the presence of juvenile wood, and in the right, higher E_R values for the part with compression wood.

Fig. 4

a Oak healthy trunk, b E_R image obtained with the proposed method, c FBP image for comparison, corresponding horizontal profiles from d the E_R image and e the FBP image. The dashed line in b and c images indicates the profile position

Juvenile wood density is lower than the mature one (Ross 2010). This can be associated with the presence of thinner cell walls and shorter tracheids in juvenile wood. Another characteristic of juvenile wood is a large microfibril angle that might be associated with low stiffness and a low Young’s modulus (Barnett and Bonham 2004). Juvenile wood has also been associated with a decreased velocity of ultrasonic waves (Brancheriau et al. 2012a; Palma et al. 2018). This wood is often considered as low-quality for many industrial uses as it typically presents low strength properties (Senft et al. 1985).

In the case of compression wood, density is higher compared to normal wood, given that the cell wall is much thicker (Timell 1986; Kollmann and Côté 2012). With the proposed method, regions affected by these density changes could lead to a biased value, considering that a single density value was fixed for the whole disk. Compression wood also has a greater microfibril angle compared to normal wood, leading to a higher stiffness in the RT plane and lower stiffness in the longitudinal direction (Brancheriau et al. 2012b; Gardiner et al. 2014). Then, the velocity of propagation of ultrasonic waves in compression wood is higher compared to normal wood (Saadat-Nia et al. 2011; Brancheriau et al. 2012b). This effect has been used for the detection of compression wood using ultrasound in previous approaches (Bucur 2006).

3.2 Centered and off-centered defects

For the centered defect position, the tested disk is shown in Fig. 5a. The hole drilled to simulate the defect was centered with respect to the pith that was not located in the geometric center of the disk, due to the presence of compression wood. In Fig. 5b, the trajectories computed iteratively in the inversion procedure are shown, and it is possible to observe how after convergence of the method these trajectories avoid the defective area. The reconstructed images obtained with the proposed method and the FBP method are shown in Fig. 5 c and d respectively. The proposed method allowed a more detailed defect identification, as the void region is distinct compared to the smooth variation obtained with the FBP image. This enhanced contrast between the defect and the healthy wood using the proposed method is linked to the use of curved rays, considering this kind of trajectories avoid the low-speed propagation areas. In contrast, the FBP method tends to smooth the output image (Kak and Slaney, 2001). The corresponding horizontal profiles are shown in Fig. 5 e and f. These profiles corroborate that the defect shape is better approximated with the proposed method, resulting in sharper defect borders.

Fig. 5

a Pine with a centric defect, b rays obtained after the reconstruction procedure, c E_R image obtained with the proposed method, d FBP image for comparison, corresponding horizontal profiles from e the E_R image and f the FBP image. The dashed line in c and d images indicates the profile position. The dashed line in e and f images indicates the defect position

To compare quantitatively the ability of defect detection, the thresholding process was applied. Figure 6 presents the thresholding results, based on the mean values in the reconstructed images. The mean velocity in the FBP image was 1527 m/s (σ=553 m/s). For the E_R image, the mean value was 1361 MPa (σ=289 MPa). In the case of the FBP image, using the first two threshold values, no defect was found. For the threshold values of 70% and 90%, the true positive rate was 33% and 96%, respectively, with false-positive rates of 7% and 29%. Improvements were obtained using the proposed method (E_R image), with true-positive rates of 61%, 66%, 75%, and 95% as the threshold value increased, and false-positive rates of 5%, 8%, 11%, and 25% respectively. The E_R image with a threshold of 70% presented the best balance.

Fig. 6

Thresholding of the FBP images (top) and the ER images (bottom) for the centered defect case. Dashed circle represents the true geometry of the defect. White regions correspond to the pixels classified as defective

For the off-centered defect position, the tested disk is shown in Fig. 7a. Iteratively estimated trajectories are shown in Fig. 7b, and there is a convergence towards the result. The images obtained are presented in Fig. 7c for the reconstruction using the proposed method and in Fig. 7d using the FBP method. As for the centered defect case, the E_R image presented a more contrasted defect region compared to the smooth variation obtained with the FBP image. Figure 7 e and f present a vertical profile for each reconstruction, with a better-defined defect boundary for the E_R image.

Fig. 7

a Pine with an off-centered defect, b rays obtained after the reconstruction procedure, c E_R image obtained with the proposed method, d FBP image for comparison, corresponding horizontal profiles from e the E_R image and f the FBP image. The dashed line in c and d images indicates the profile position. The dashed line in e and f images indicates the defect position

Figure 8 presents the segmented images for the off-centered case. For the E_R image, the mean value was 1250 MPa (σ = 812 MPa), and for the FBP image, the mean velocity was 1330 m/s (σ = 251 m/s). Using the threshold value of 30% for the FBP image, defective pixels were not detected. For the other threshold values, the true-positive rates were 7%, 62%, and 91%, and the respective false-positive rates were 5%, 9%, and 21%. For the E_R images, true-positive rates were improved compared to the FBP results, with values of 73%, 84%, 94%, and 95% as the threshold increased, and false-positive rates of 3%, 7%, 26%, and 50% respectively. Segmentation using the E_R image and a 50%threshold led to the best compromise. Compared to the centered defect case, here we obtained higher false-positive rate values. It has been shown that off-centered defects resulted in lower TOF variations compared to the centric defect cases, leading to less contrasted regions in the tomographic images compared to centric cases, which increases the difficulty in identifying defect location (Espinosa et al. 2019a).

Fig. 8

Thresholding of the FBP images (top) and the E_R images (bottom) for the off-centered defect case. Dashed circle represents the true geometry of the defect. White regions correspond to the pixels classified as defective

4 Discussion

The evaluation process of standing trees, especially in an urban context, involves subjective assessments of various parameters to be measured in decay detection and loss of resistance from external signs. However, the use of tools such as tomography in this process provides greater support and the ability to have detailed and accurate information about the internal state of standing trees. This allows making better decisions based on a deeper understanding of risk factors for failure in trees. The proposed method may facilitate the interpretation of the results and increase the accuracy of the extension and location of decomposition, avoid invasive tests, and enable timely intervention of health affectations of trees in the technical evaluation in accident prevention and silvicultural management of trees. In this study, we tried to design the reconstruction methodology taking into account the following criteria: a fast computing time with a common computer, a low number of sensors used in situ, and the possibility to implement the proposed method in the existing commercial devices.

The common reconstruction techniques for acoustic or ultrasonic imaging of standing trees are based on the hypothesis of quasi-isotropic behavior in the cross-section. The propagation velocity of waves is supposed to be independent of the angle of propagation between the direction vector and the radial axis; it is thus possible to construct an image representing a velocity value for each pixel of the image. However, the hypothesis of quasi-isotropic behavior is not true, and it was proposed in this study a theoretical method allowing computing a specific elastic constant for each pixel. The specific elastic constant can be Young’s modulus in the radial direction, in the tangential direction, or the shear modulus, divided by the density. Poisson’s coefficient was set to a constant (ν_RT = 0.3) to improve the convergence of the algorithm. Furthermore, it was chosen to present the results with the radial modulus of elasticity and a constant value of the density equal 661 kg/m³ (Ross 2010). To compute more precise values of elastic properties using the proposed method, it would be necessary to determine the variation of the density within the cross-section of the tree using another non-destructive technique than elastic waves propagation analysis.

The sinogram was interpolated to virtually add sensors and increase image resolution while limiting the number of sensors really used (which allows reducing the in situ time inspection). A sinogram grouped the set of experimental times-of-flight obtained for all the emitters (projections) in a matrix form (the rows corresponded to the emitters and the columns to the receivers). Assuming that the number of emitters used was sufficient to detect the presence of a defect and that the physical quantity measured (time-of-flight) was continuous between two sensors, it was thus possible to determine this quantity for a virtual sensor placed between two real adjacent sensors by interpolation. In this case, it was chosen to use a linear interpolation procedure, and the number of sensors was duplicated, passing from 16 to 32 sensors. The final number of sensors (interpolated sinogram) was only limited by the computing time. Sinogram interpolation allowed to obtain more accurate results than interpolating the reconstructed image because the final image includes additional errors from the reconstruction process.

It was previously shown that the presence of bark significantly contributed to signal attenuation (Brancheriau et al. 2012a). Furthermore, the irregularity of the bark induces a problem of coupling between the ultrasonic sensor and the wood. In this study, it was decided to remove the bark because the aim was to validate experimentally the reconstruction method of the inner part of a tree. For in situ ultrasonic testing, it would be necessary to adapt or to design specific sensors that optimize the coupling with the wood. A possibility would be to screw the sensor into the wood through the bark.

The experimental tests were performed on green wood with an average moisture content of 92%. The protocol was notably designed to avoid wood drying, and the time between the tree harvesting and the end of the ultrasonic tests was approximatively 1 month. Maintaining wood above the fiber saturation point allowed mimicking the state of the matter in standing tree, while having the possibility to create controlled defects in the wood discs. It was reported that variation of moisture content above the fiber saturation point had an impact on the propagation of ultrasonic waves in wood (Sakai et al. 1990; Unterwieser and Schickhofer 2011; Yamasaki et al. 2017). The harvesting and the period of storage would have softened the gradient of moisture content within the wood disks and produced a quasi-homogeneous hygroscopic state. Taking into account this last hypothesis, the influence of the moisture content variation was not considered in the proposed model. In addition, the internal moisture content gradient is difficult to assess experimentally during in situ field-testing on standing trees. An extension of this study would focus on the effect of the internal variation of moisture content on tomographic image reconstruction and on the correction of the computed characteristic parameters.

The proposed reconstruction method was based on the bi-dimensional hypothesis (transverse cross-section of a tree). This hypothesis is encountered in all commercial tomographic devices. Therefore, three-dimensional imaging is a result of the interpolation of bi-dimensional images (Martinis et al. 2004). Due to the difference of wave velocities between the longitudinal axis of a tree and the transverse axes, the shape of the wavefront is more altered in the longitudinal direction. The consequence is that a time-of-flight measured between two sensors placed at a given height of a tree can be the result of a wave traveling with a 3D path. This is why wood disks were used in this study and not tree trunks. In order to continue improving the tomographic reconstruction procedure of standing trees, it would be necessary to develop a full 3D inversion technique (Goncharsky et al. 2014). A full 3D setup would require a distribution of sensors at different heights of a tree (not necessarily around cross-sections). The reconstruction method would consider the cylindrical orthotropic condition of wood in the three directions to model the wave propagation. The solution of this inverse problem does not exist up to now but would need high computing capacities, as it would use a large amount of data by considering all the possible transmitter-receiver combinations.

5 Conclusions

With the purpose of improving the detection of the presence of pests and diseases in standing trees, an alternative procedure for the tomographic imaging of wood by elastic waves was proposed. This procedure was developed for ultrasonic waves but is also valid for acoustic waves. Common tomographic devices mainly use acoustic waves, and the informatics algorithms developed for this work can be integrated into these commercial devices. The proposed procedure gave a solution to the inverse problem of tomographic imaging by considering the orthotropic behavior of the wood matter. However, the reasoning shown here can also be applied for buildings structural health monitoring or medical imaging, in the case of anisotropic behaviors, as examples. The orthotropic behavior of wood had the consequence that the local wave velocity was not an intrinsic parameter of the matter (despite this, images of local velocities are however used by common tomographic devices) because the wave velocity was a function of the propagation angle between the acoustic ray and the local orthotropic axes. The proposed procedure allowed determining the local values of the specific elastic parameters (Young’s moduli in radial or tangential direction, or shear modulus, divided by density) that are intrinsic parameters of the matter, and to determine the propagation trajectories in the same computing process. It mainly consisted of an iterative process focused on finding a polynomial approximation of the slowness at each pixel in the image as a function of the propagation angle, modifying the curved trajectories by a raytracing approach. The proposed inversion executed fast considering in situ testing. A specific interpolation procedure was implemented allowing limiting the number of sensors used, which can be less time consuming for field expertise. Based on the hypothesis of the continuity of the measured physical property (time-of-flight) between sensors, the sinogram interpolation increased the spatial resolution of the image, which can improve the quantification of the geometry and size of defects. It was demonstrated by experimental tests that the imaging process led to a better defect detection compared to the use of a common imaging technique. The value of the intrinsic parameter (specific modulus in radial direction) in the experimental tomographic image was modified by the presence of juvenile wood and by the presence of compression wood. It was thus possible to detect the compression wood by this imaging procedure. Further work should consider two main research axes to obtain a better representation of the tree inner state using elastic wave tomography. On one hand, the effect of the moisture content above the fiber saturation point, and its variability within the cross-section of the tree, on the physical parameters measured and reconstructed should be taken into account in the imaging process. On the other hand, the wave propagation in the three directions of the space should be considered by developing a specific experimental methodology and by developing a 3D inversion procedure including orthotropy.