I would like to express my deep gratitude to my research supervisors Flavio Prieto and Loïc Brancheriau, for their patient guidance, enthusiastic encouragement and useful critiques of this research work. I would also like to thank Phillipe Lasaygues, whose expertise was invaluable for the success of this thesis. I would like to offer my special thanks to the members of the thesis committee in France throughout these years: Yann Monerie, Alfredo Napoli, Joseph Gril, Claire Atger, Philippe Gallet, Yves Caraglio and Emmanuel Le Clézio. In Colombia, I would also like to thank Jan Bacca and Claudia Patricia Pérez for participating in the evaluation of my thesis proposal and qualifying exam.

I wish to acknowledge Sylvain Lotte and Daniel Guibard, for his capable and willing technical assistance in the processing of wood samples. I would also like to extend my thanks to Pierre Aversenq, Sylvain Dujardin and Yolima Cortes, arboricultural experts, for letting us accompany them in the technical assessment of standing trees.

I would like to thank my fellow doctoral students in Colombia and France for their feedback, cooperation and of course friendship. In addition, I would like to express my gratitude to the staff of CIRAD; assistance provided by Isabelle Chalon was greatly appreciated.

Finally, I wish to thank my family: my parents and my brothers for supporting me throughout writing this thesis. Some special words of gratitude go to my friends who have always been a major source of support. i

SUMMARY

The tree plays a major ecological role in modern cities. The management of the plants is the subject of requests from urban operators: the diagnosis is essentially visual, even when the extent of internal damage and the associated hazard cannot be precisely evaluated by simple observation. Ultrasonic imaging methods allow answering biological questions related to the adaptation of the tree to exogenous constraints, such as pathogenic attacks, presence, and type of internal damages, the extent of degraded or traumatized areas.

The major scientific issues are linked to the image production (reconstruction of the intrinsic parameter from a set of measurements) and to the image interpretation (discrimination for detection of alterations and its positioning). The overall aim of this thesis was to develop an ultrasonic imaging method for the diagnosis of the internal condition of urban trees. The scientific objectives were to develop numerical models to study the factors of influence on the propagation of ultrasonic waves in the cross-section of a tree and to propose an image reconstruction solution, suited for orthotropic materials, allowing the discrimination and positioning of decay. The development of a protocol for the acquisition, processing, analysis and interpretation of ultrasound tomography signals and images is of great importance for wood science. Obtaining reliable and interpretable images is a recurring demand from urban operators.

Initially, to set-up the ultrasonic chain of measurement, a comparative experimental study was done to choose the excitation signal parameters, such as shape, temporal duration, and frequency response, and then the choice of a suitable time-of-flight determination technique. Then, we were concerned on evaluating the influence of the orthotropic condition of wood on the propagation of ultrasonic waves, by performing a time-of-flight (TOF) estimation using a raytracing approach, a method used in the field of exploration seismography to simulate wavefronts in elastic media. The anisotropy of wood in the radial-tangential plane influenced the wave velocity depending on the direction of propagation, that led to deformed wavefronts compared to the perfectly circular wavefronts for an isotropic case. The paths from each receiver to the transmitter in the wood presented a curvature, therefore the trajectories differed from the straight-line distance obtained for an isotropic case. A numerical comparison was made using the Finite Elements Method (FEM); the TOF estimates and wavefronts agreed with those of the raytracing approach. A similar experimental validation was performed. Wood sections from two species were SUMMARY vi tested. Defects in the wood were simulated by drilling holes. The shape of TOF curves computed using the raytracing algorithm and those obtained from the experiments were in good agreement. Defects located in the center of the trunk presented larger TOF variations compare to defects located in off-centered positions. Thus, off-centered defects would be more difficult to determine and characterize by tomographic inversion.

Then, we were interested in the influence of the wood orthotropic condition on the tomography image reconstruction process (inverse problem) and how it should be adapted to the standing tree constraints. For wood, the ray paths between the ultrasonic transmitter and the receivers are not straight as for isotropic media; therefore, the image reconstruction method should be adapted to deal with curved rays. The proposed method considers the orthotropy property of wood material, performing an iterative process that approximated the curved rays. A slowness function was defined for every pixel and a nonlinear regression allowed the mapping of the inner elastic constants. Initially, four numerical configurations were tested representing real cases usually found in standing tree monitoring. The reconstructed images using the proposed method were compared with a straight-ray reconstruction method (filtered back projection algorithm), highlighting a more detailed identification and quantification of the inner state of the anisotropic structure of the trunk.

Then, the inversion procedure was tested using wood samples from two species for three different configurations: a healthy case, a centered defect case, and an off-centered defect case. As for the numerical study, the proposed method resulted in a more accurate defect representation when compared to a straight-ray reconstruction, especially for the case of

RESUMÉ viii

Les méthodes classiques de diagnostic interne des arbres sont basées soit sur des mesures de perçage, soit sur des mesures acoustiques ou électriques ponctuelles. Ces méthodes restent limitées, les informations fournies sont principalement locales. Les méthodes d'imagerie permettent d'analyser la structure interne des arbres (reconstruction bi ou tridimensionnelle) sans altérer leur état. L'imagerie du bois peut être réalisée en utilisant des ondes électromagnétiques (rayons X) mais la mise en oeuvre sur le terrain est délicate.

Une alternative aux ondes électromagnétiques consiste à utiliser des ondes élastiques (tomographie acoustique ou ultrasonore) présentant comme avantages un coût réduit et une mise en oeuvre plus facile. Les approches d'imagerie utilisant des ondes élastiques consistent en une tomographie quantitative (bidimensionnelle dans la très grande majorité des cas). L'image est appelée quantitative car il s'agit d'une cartographie d'un paramètre mécanique intrinsèque du matériau (déduit de la mesure de la vitesse de propagation des ondes dans la section transverse de l'arbre au niveau du sondage). Les dispositifs commerciaux les plus connus utilisent des ondes acoustiques, tels que Arbotom, Picus et Fakkop. Des sondes acoustiques sont positionnées en contact tout autour de l'arbre et un marteau est utilisé pour impacter le tronc, produisant des ondes à des fréquences audibles (inférieures à 20 kHz). Ces ondes se propagent dans le tronc et sont captées par les sondes acoustiques. De nombreuses études ont évalué les performances de ces dispositifs acoustiques, soulignant différents problèmes : une faible résolution spatiale et des images difficiles à interpréter. L'approche par ultrasons vise à augmenter la fréquence des ondes analysées pour obtenir des images avec une résolution spatiale plus élevée.

La faisabilité de la tomographie par ultrasons pour la détection de la décomposition du bois dans les arbres vivants a été évaluée dans plusieurs études ; les auteurs ont indiqué que les techniques ultrasonores sont appropriées pour l'évaluation de la qualité des arbres sur pied. Cependant, cette technique telle qu'utilisée actuellement présente quelques inconvénients : l'effet de l'anisotropie du bois dans la reconstruction de l'image n'est pas pris en compte (l'image obtenue est biaisée), la mesure de la vitesse de propagation est imprécise (nécessité de répéter les essais), les signaux ultrasonores sont fortement atténués (le dispositif expérimental est plus complexe à mettre en oeuvre que pour la tomographie acoustique).

Afin d'améliorer la tomographie par ultrasons, il est nécessaire de prendre en compte la complexité du matériau bois et de développer des techniques de traitement du signal et de reconstruction d'image adaptées à cette complexité. La première question ix abordée est liée au problème direct. Il s'agit de modéliser la propagation des ondes élastiques dans un milieu orthotrope à géométrie cylindrique (la section transverse du tronc d'un arbre). La deuxième question scientifique concerne l'algorithme de construction de l'image tomographique, appelé problème inverse. Cet algorithme doit prendre en compte des rayons courbes dans un milieu anisotrope pour déterminer les propriétés intrinsèques du matériau en tout point de la section transverse du tronc. Les images ainsi obtenues doivent être comparées à celles obtenues par la technique classique utilisant une hypothèse d'isotropie.

L'équation de Christoffel est généralement utilisée pour décrire la relation entre la vitesse de l'onde ultrasonore et les paramètres mécaniques intrinsèques du matériau bois.

La propagation des ondes dans un solide élastique anisotrope homogène peut être décrit par un ensemble d'équations différentielles. L'équation de Christoffel conduit à une solution pour cet ensemble d'équations sous la forme d'ondes planes, reliant la vitesse de propagation aux constantes élastiques du matériau et à la direction de propagation des ondes. L'anisotropie du bois dans le plan radial-tangentiel a une influence directe sur la vitesse de l'onde en fonction de la direction de propagation. L'évolution de la vitesse en fonction de la direction de propagation dépend de l'essence de bois considérée, et une différence a été constatée entre les résineux et les feuillus. La direction radiale correspond à la vitesse d'onde la plus rapide. L'influence des paramètres mécaniques du plan radialtangentiel du bois sur le calcul de la vitesse de propagation utilisant l'équation de Christoffel a été évaluée par une étude de sensibilité. Les paramètres mécaniques de six essences de bois ont été choisis. L'analyse de sensibilité a montré une influence plus grande des modules de Young, suivis du coefficient de Poisson et enfin du module de cisaillement.

Cependant, ces deux derniers paramètres ne doivent pas être négligés lors de l'utilisation de l'équation de Christoffel pour résoudre le problème inverse en tomographie (négliger ces paramètres revient à introduire un biais non négligeable).

Une étude a été réalisée pour déterminer les paramètres du signal ultrasonore d'excitation, tels que la forme, la durée et la réponse en fréquence ; puis pour sélectionner une technique de détermination du TOF. La mesure du temps de propagation est en effet une étape critique dans l'exécution d'essais non destructifs par ultrasons des arbres sur pied, avec une influence directe sur la précision de la détection des défauts. Un arbre vivant a été testé en plaçant des capteurs (avec des fréquences de résonance de 36 kHz et 60 kHz) dans quatre positions différentes, avec cinq signaux d'excitation différents et trois méthodes de RESUMÉ x détection du TOF. Parmi toutes les configurations, celle qui présentait le moins de variations sur les mesures de TOF était la combinaison d'un signal « chirp » (signal modulé en fréquence autour d'une fréquence porteuse) avec la méthode de corrélation croisée. Le signal « chirp » a été ajusté à la réponse du capteur (gamme utile autour de la fréquence de résonance). La position des capteurs a eu un effet sur la mesure : lorsque l'angle de position du capteur approchait de la direction radiale, les valeurs de TOF présentaient moins de variabilité (la mesure était plus stable).

Un modèle numérique a ensuite été développé, avec l'équation de Christoffel, pour simuler la propagation des ondes dans le bois et déterminer le temps de propagation (ou temps de vol, TOF) de l'onde (paramètre physique permettant la construction de l'image tomographique). La méthode de « raytracing » a été utilisée pour ce modèle. Cette méthode est notamment employée dans la géophysique pour simuler les fronts d'onde dans les milieux élastiques. Le « raytracing » construit les fronts d'onde successivement à partir des précédents avec un pas de temps constant. L'anisotropie dans le plan radial-tangentiel du bois modifie la forme des fronts d'onde par rapport au cas d'un matériau isotrope. Les rayons entre émetteur et récepteur sont courbes. Par conséquent, les distances parcourues sont différentes des distances en ligne droite obtenues pour le cas isotrope. Afin de comparer et de valider les résultats obtenus avec l'approche « raytracing », la méthode des éléments finis (FEM) a été utilisée pour modéliser la propagation des ondes élastiques dans le bois. Le modèle FEM a abouti à des estimations des TOF très proches de celles obtenues avec l'approche « raytracing ». Les fronts d'onde dans les deux cas concordent et les valeurs de TOF estimées avaient une différence relative inférieure à 2% entre les deux modèles. La méthode « raytracing » est moins complexe que le modèle FEM, ce qui entraine des temps de calcul plus courts.

Une validation expérimentale du modèle « raytracing » a été effectuée sur des disques de deux essences (chêne et frêne). Des défauts dans le tronc ont été créés en perçant des trous. Ces défauts ont été testés dans deux positions (centrée et excentrée) et avec trois diamètres différents pour chaque position. Les expériences réalisées ont permis d'obtenir des profils de temps de propagation similaires à ceux obtenus par le modèle numérique. La présence des défauts centriques a eu pour effet des variations de TOF plus importantes que celles des défauts excentrés. Un défaut centré et de grande taille est plus facilement détectable.

xi Les approches classiques de tomographie des arbres sur pied utilisent des techniques de reconstruction prenant en compte des trajectoires rectilignes (rayons), comme la technique de rétroprojection filtrée (FBP) et les méthodes algébriques, comme par exemple la technique de reconstruction itérative simultanée (SIRT). L'hypothèse des trajectoires rectilignes produit une image biaisée (le matériau est supposé isotrope). De plus, ces approches fournissent une cartographie de vitesse, mais étant donné que ce paramètre varie en fonction de l'angle de propagation par rapport au repère local d'orthotropie, le sens physique de l'image résultante n'est pas clair.

Une méthode de reconstruction d'image tomographique bidimensionnelle adaptée au bois a été développée dans cette thèse. La méthode proposée prend en compte l'orthotropie locale du matériau avec une géométrie cylindrique ; c'est un processus itératif qui reconstruit à la fois les rayons de propagation et les propriétés intrinsèques locales du matériau. Dans ce processus, une fonction de lenteur (inverse de la vitesse de propagation fonction de l'angle de propagation) est définie pour chaque pixel de l'image et une régression non linéaire permet de déterminer les constantes élastiques spécifiques (constantes divisées par la densité) à partir des paramètres de cette fonction de lenteur.

RESUMEN

Título en español: Tomografía ultrasónica para la gestión sostenible y racional de árboles en pie en áreas urbanas Los árboles urbanos desempeñan un importante papel ecológico, sanitario y estético en las ciudades modernas. Estos beneficios incluyen reducir la contaminación del aire, brindar refugio y alimento a muchas especies de aves y animales pequeños, proteger contra el calor y mejorar el paisaje. Sin embargo, algunos riesgos están asociados con los árboles urbanos: la caída de un árbol (o rama) puede causar daños a los transeúntes, vehículos y casas. La evaluación del riesgo de los árboles en las ciudades a menudo está a cargo de operadores urbanos. El diagnóstico es esencialmente visual, mientras que la extensión del daño interno y el peligro asociado no pueden evaluarse con precisión mediante una simple observación. Con frecuencia, los árboles con descomposición interna parecen sanos. La respuesta a la gestión de las plantaciones urbanas requiere, por lo tanto, el uso de métodos de diagnóstico interno in situ para los árboles. Se pueden utilizar técnicas modernas para minimizar los riesgos asociados con la caída de árboles; ha habido un progreso significativo en los equipos de control y diagnóstico, así como en las fórmulas y pautas para evaluar árboles peligrosos.

Tuvimos la oportunidad de participar en un proceso de inspección de árboles in situ, para evaluar su nivel de peligro, con un experto forestal. Los árboles en cuestión pertenecían al género Platanus, ubicado a lo largo de una carretera en la localidad de Lodève, diagnóstico solicitado por el consejo del departamento de Hérault (Francia). El diagnóstico se divide en dos fases: evaluación visual y evaluación especializada con instrumentos. Dependiendo del entorno local y su importancia patrimonial, el diagnóstico de un árbol puede implicar ambas fases o solo la evaluación visual. Por ejemplo, los árboles que inspeccionamos, que podríamos definir como árboles ornamentales, normalmente se evalúan visualmente, a menos que muestren signos significativos de degradación. Si ese es el caso, se utilizan técnicas usando instrumentos, como por ejemplo medidas de perforación. El diagnóstico debe ser rápido en su ejecución. Algunos árboles urbanos con gran importancia patrimonial pueden requerir un análisis más profundo y demorado.

Técnicas como la tomografía acústica y la ultrasónica se convierten en herramientas indispensables.

RESUMEN xvi

Los métodos estándar de diagnóstico interno de árboles se basan en mediciones de perforación o en mediciones acústicas o eléctricas. Estos métodos son limitados, la información proporcionada es principalmente local. Los métodos de obtención de imágenes tomográficas pueden analizar la estructura interna de los árboles (reconstrucción en dos o tres dimensiones) sin alterar su estado. En el caso de la madera, estas imágenes se pueden obtener utilizando ondas electromagnéticas (rayos X), pero la implementación en campo es difícil. Una alternativa a las ondas electromagnéticas es utilizar ondas elásticas (tomografía acústica o ultrasónica) con las ventajas de un costo reducido y una implementación más sencilla. Los métodos que utilizan ondas elásticas consisten en una tomografía cuantitativa (bidimensional en la gran mayoría de los casos). La imagen se llama cuantitativa porque es un mapeo de un parámetro mecánico intrínseco del material (deducido de la medición de la velocidad de propagación de las ondas en la sección transversal del árbol en el nivel de sondeo). Los dispositivos comerciales más conocidos utilizan ondas acústicas, como por ejemplo Arbotom, Picus y Fakkop. Los sensores acústicos se colocan en contacto alrededor del árbol y se usa un martillo para impactar el tronco, produciendo ondas en frecuencias audibles (por debajo de 20 kHz). Estas ondas se propagan en el tronco y son captadas por los sensores acústicos. Numerosos estudios han evaluado el desempeño de estos dispositivos acústicos, destacando varios problemas: baja resolución espacial e imágenes difíciles de interpretar. El enfoque de ultrasonido tiene como objetivo aumentar la frecuencia de las ondas analizadas para obtener imágenes con una resolución espacial más alta.

La viabilidad de la tomografía ultrasónica para la detección de la descomposición de la madera en árboles vivos se ha evaluado en varios estudios; los autores han indicado que las técnicas ultrasónicas son apropiadas para evaluar la calidad de los árboles en pie.

Sin embargo, esta técnica en su estado actual tiene algunas desventajas: el efecto de la anisotropía de la madera en la reconstrucción de la imagen no se tiene en cuenta (la imagen obtenida está sesgada), la medición de la velocidad de propagación es imprecisa (es necesario repetir las pruebas), las señales ultrasónicas están muy atenuadas (el dispositivo experimental es más complejo de implementar que en el caso de la tomografía acústica).

Para superar las desventajas de la tomografía ultrasónica, es necesario tener en cuenta la complejidad del material de madera y desarrollar técnicas de procesamiento de señales y reconstrucción de imágenes adaptadas a esta complejidad. El primer problema está relacionado con el problema directo. Se trata de modelar la propagación de las ondas xvii elásticas en un medio ortotrópico con geometría cilíndrica (la sección transversal del tronco de un árbol). La segunda pregunta científica se refiere al algoritmo de construcción de la imagen tomográfica, llamado problema inverso. Este algoritmo debe considerar los rayos curvos en un medio anisotrópico para determinar las propiedades intrínsecas del material en cualquier punto de la sección transversal del tronco. Las imágenes obtenidas de este modo deben compararse con las obtenidas mediante la técnica convencional utilizando una hipótesis de isotropía.

La ecuación de Christoffel se usa generalmente para describir la relación entre la velocidad de la onda ultrasónica y los parámetros mecánicos intrínsecos del material de madera. La propagación de la onda en un sólido elástico anisotrópico homogéneo se puede describir mediante un conjunto de ecuaciones diferenciales. La ecuación de Christoffel conduce a una solución para este conjunto de ecuaciones en forma de ondas planas, relacionando la velocidad de propagación a las constantes elásticas del material y a la dirección de propagación de las ondas. La anisotropía de la madera en el plano radialtangencial tiene una influencia directa sobre la velocidad de la onda en función de la dirección de propagación. La evolución de la velocidad en función de la dirección de propagación depende de la especie de madera considerada, y se observó una diferencia entre las maderas blandas y las maderas duras. La dirección radial corresponde a la velocidad de onda más rápida. La influencia de los parámetros mecánicos del plano radialtangencial de la madera en el cálculo de la velocidad de propagación utilizando la ecuación de Christoffel se evaluó mediante un estudio de sensibilidad. Se eligieron los parámetros mecánicos de seis especies de madera. El análisis de sensibilidad mostró una mayor influencia de los módulos de Young, seguido por el coeficiente de Poisson y finalmente el módulo de cizalla. Sin embargo, estos dos últimos parámetros no deben descuidarse al usar la ecuación de Christoffel para resolver el problema inverso en tomografía (descuidar estos parámetros introduce un sesgo considerable).

Se realizó un estudio para determinar las características de la señal de excitación ultrasónica a utilizar, como lo son la forma, la duración y la respuesta de frecuencia; a continuación, se debía seleccionar una técnica de determinación de TOF. La medición del tiempo de propagación es un paso crítico en la realización de pruebas ultrasónicas no destructivas de árboles en pie, con una influencia directa en la precisión de la detección de defectos. Se realizaron pruebas en un árbol en pie colocando sensores (con frecuencias de resonancia de 36 kHz y 60 kHz) en cuatro posiciones diferentes, con cinco señales de RESUMEN xviii excitación diferentes y tres métodos de detección de TOF. De todas las configuraciones, la que tuvo la menor variación en las mediciones de TOF fue la combinación de una señal "chirp" (una señal modulada en frecuencia alrededor de una frecuencia portadora) con el método de detección usando correlación cruzada. La señal "chirp" se ajustó a la respuesta en frecuencia del sensor (bandas alrededor de la frecuencia de resonancia). La posición de los sensores tuvo un efecto en la medición: cuando el ángulo de posición del sensor se acercaba a la dirección radial, los valores de TOF mostraban una menor variabilidad (la medición era más estable).

Luego se desarrolló un modelo numérico utilizando la ecuación de Christoffel, para simular la propagación de las ondas en la madera y determinar el tiempo de propagación (o tiempo de vuelo, TOF) de la onda (parámetro físico que permite la construcción de la imagen tomográfica). Se utilizó el método conocido como "raytracing" para este modelo.

Este método se usa particularmente en geofísica para simular frentes de onda en medios elásticos. El método raytracing construye frentes de onda sucesivamente con un paso de tiempo constante. La anisotropía en el plano radial-tangencial de la madera modifica la forma de los frentes de onda en comparación con el caso de un material isotrópico. Las trayectorias entre el transmisor y los receptores son curvas. Por tanto, las distancias recorridas son diferentes a las distancias en línea recta obtenidas para el caso isotrópico.

Para comparar y validar los resultados obtenidos con el método raytracing, se utilizó el método de elementos finitos (FEM) para modelar la propagación de ondas elásticas en la madera. El modelo FEM dio como resultado estimaciones de TOF muy cercanas a las obtenidas con el enfoque de trazado de rayos. Los frentes de onda en ambos casos concuerdan y los valores de TOF estimados tuvieron una diferencia relativa de menos del 2% entre los dos modelos. El método raytracing es menos complejo que el modelo FEM, lo que resulta en tiempos de cálculo más cortos.

La validación experimental del modelo raytracing se llevó a cabo en discos de madera de dos especies (roble y fresno). Los defectos en el tronco se crearon perforando agujeros. Estos defectos se probaron en dos posiciones (centrada y excéntrica) y con tres diámetros diferentes para cada posición. Los experimentos realizados permitieron obtener perfiles de tiempo de propagación similares a los obtenidos con el modelo numérico. Los defectos ubicados en el centro causaron mayores variaciones en los valores de TOF que los xix defectos excéntricos en comparación con el caso sano (antes de la perforación). Un defecto más grande y céntrico es detectado con más facilidad.

Los enfoques tradicionales de tomografía de árboles en pie utilizan técnicas de reconstrucción que toman en cuenta trayectorias rectilíneas (rayos), como la técnica de retroproyección filtrada (FBP) y métodos algebraicos, como la técnica de reconstrucción iterativa simultánea (SIRT). La hipótesis de trayectorias rectilíneas produce una imagen sesgada (se supone que el material es isotrópico). Además, estos enfoques proporcionan mapas de velocidad, pero como este parámetro varía con el ángulo de propagación en relación con el punto de referencia ortotrópico local, el significado físico de la imagen resultante no está claro.

En esta tesis se ha desarrollado un método de reconstrucción de imágenes tomográficas bidimensionales adaptado a la madera. El método propuesto tiene en cuenta la ortotropía local del material con una geometría cilíndrica; es un proceso iterativo que reconstruye las trayectorias de propagación (rayos) y las propiedades locales intrínsecas del material. En este proceso, se define una función de lentitud (inversa de la velocidad de propagación en función del ángulo de propagación), para cada píxel de la imagen y una regresión no lineal permite determinar las constantes elásticas específicas (constantes divididas por la densidad), de los parámetros de esta función de lentitud.

Primero debemos considerar una imagen dividida en N píxeles atravesados por los rayos ultrasónicos. El objetivo es estimar un parámetro local para cada píxel (la lentitud α), que sumado a lo largo del radio permite determinar el TOF entre el transmisor y el receptor.

Para un par de sensores cuya trayectoria m pasa a través de los píxeles k, el tiempo de vuelo tm se puede escribir de la siguiente manera:

𝑡 𝑚 = ∑ 𝑙 𝑚𝑘 𝛼 𝑘 𝑘 𝑎𝑙𝑜𝑛𝑔 𝑚
Con lmk la longitud de un segmento del rayo de propagación. La resolución de todas las ecuaciones de los pares de sensores conduce a la reconstrucción del parámetro interno αk (problema inverso). Sin embargo, surgen dos problemas: en primer lugar, la formulación matricial debe adaptarse a la dependencia de la lentitud en el ángulo de propagación (linealización de la ecuación de Christoffel), y, en segundo lugar, la estimación de las trayectorias verdaderas (desconocidas en principio), se ve afectada por la anisotropía de la madera y la presencia de defectos.

RESUMEN xx

La ecuación de Christoffel se linealizó mediante una aproximación polinómica (grado 5). La lentitud del píxel atravesado por un rayo m se expresa de la siguiente manera: El proceso iterativo comienza con trayectorias rectas. Con estas trayectorias, es posible realizar una primera inversión que conduzca a los coeficientes polinómicos de lentitud β para todos los píxeles. A partir de estos coeficientes, se resuelve el problema directo para obtener una estimación de los valores de TOF y las nuevas trayectorias. Este proceso se repite hasta que la diferencia entre los valores de TOF en relación con los valores reales se minimiza y las formas de las trayectorias calculadas ya no varían. Finalmente, la aproximación polinomial obtenida para cada píxel se usa para estimar los valores de las constantes elásticas mediante una regresión no lineal por el método de mínimos cuadrados.

𝛼 𝑘 = 𝛽
El resultado es una imagen paramétrica de los módulos específicos ER, ET y GRT considerando constante el coeficiente de Poisson νRT (esto permite obtener parámetros estables mediante la regresión no lineal).

Se realizó una validación numérica de este método de inversión. Se probaron cuatro configuraciones numéricas que representan casos reales que generalmente se encuentran en el campo: (1) un tronco en estado saludable, (2) con un orificio céntrico, (3) con un orificio excéntrico y (4) con un defecto excéntrico caracterizado por un gradiente de propiedades mecánicas. Este último caso podría estar asociado con una etapa temprana de degradación de la madera por el ataque de hongos. Las imágenes reconstruidas utilizando el método propuesto se compararon con el método de reconstrucción clásico con la hipótesis de isotropía (FBP, rayos rectos).

La comparación de las imágenes obtenidas reveló una identificación y cuantificación más detallada del estado interno del tronco con el método propuesto. El xxi método propuesto reconstruyó correctamente la curvatura de los rayos; las áreas defectuosas están mejor definidas y son más visibles. La identificación de defectos céntricos fue más precisa que los casos con un defecto excéntrico y un gradiente de propiedad, independientemente del método utilizado. El método propuesto requiere solo unos minutos de cálculo con una computadora portátil comercial; por lo tanto, se adapta al contexto de las pruebas in situ.

El método de inversión propuesto se probó experimentalmente en muestras de madera de dos especies (roble y pino), para tres configuraciones diferentes: un caso saludable, un defecto céntrico y un defecto excéntrico. Al igual que con la validación numérica, el método desarrollado permitió obtener una representación más precisa de los defectos en comparación con una reconstrucción clásica de rayos rectos, en particular en el caso de defectos céntricos.

Otros factores asociados con la variabilidad de la madera pueden incluirse en el modelado para mejorar la precisión del método. Por ejemplo, las variaciones en el contenido de humedad por encima del punto de saturación de las fibras o los cambios en la densidad afectan la velocidad de propagación de las ondas. En esta tesis, hemos considerado solamente las ondas de compresión y no se tuvieron en cuenta todos los demás fenómenos de segundo orden (refracción causada por los defectos, conversión de modo, atenuación y dispersión). Para la continuación de este trabajo, todos estos fenómenos físicos deben tenerse en cuenta para caracterizar la madera con mayor precisión. Se podría considerar el uso del método conocido como "Full Wave Inversion" (FWI). El término "completo" se refiere al uso de toda la señal acústica, es decir, sin rechazar información potencialmente útil. Hoy en día, estos métodos todavía no están adaptados a la madera.

Considerando una implementación 3D del método de reconstrucción, es necesaria una instalación de transductores a diferentes alturas y la adaptación del problema directo a los frentes de onda que se propagan en la dirección longitudinal, radial y tangencial.

xxiii Trees play a major ecological and sanitary role in modern cities. Some benefits are a reduction of air pollution, energy savings from solar shading, wildlife shelter, and food, among multiple others [START_REF] Pokorny | Urban tree risk management: A community guide to program design and implementation[END_REF]. However, some risks are associated with urban trees. For instance, trees (or tree parts) fall can cause damages to passers-by, vehicles, and houses (Figure 1). Trees planted in urban areas are exposed to conditions that can weaken its structure, such as soils that are not appropriate for plants, drought, or traffic shocks that wounds the trunk. Trees can be weakened because of these conditions and deteriorated by the attack of insects or fungi. The risk evaluation associated with trees in cities is often in charge of urban operators. Sustainable and rational management of the tree patrimony is based on the census and expertise of "all" the individuals of a composition. As an example, in Bogotá, the Secretary of Environment and the Botanical Garden oversee the tree patrimony (more than 1.2 million trees in public space), which includes silvicultural practices such as planning and monitoring of the phytosanitary state and physiological condition of trees. As another example, the city of Paris, one of the wooded capitals in Europe, has the Tree and Wood Department of the town hall in charge of the monitoring and maintenance of 200.000 trees. Each tree is monitored by his own "digital ID card" gathering all the information regarding his date of planting, successive watering, pruning, state of health (physiological state, wounded areas, fungi), to help on the diagnosis of dangerous trees. Rather than considering security fell, the tree maintenance is preferred until they are integrated into the renovation plan. This is made possible through monitoring, and constant maintenance.

LIST OF SYMBOLS

Economical aspects related to the management of the plants are the subject of requests from urban operators. Today, many communities in France or Colombia, cannot afford an extensive diagnosis with an expert by visiting the whole tree population.

Moreover, the diagnosis is essentially visual, even when the extent of internal damage and the associated danger cannot be precisely evaluated by simple observation. Frequently, trees presenting inner decomposition seem healthy in appearance; failure symptoms are not always seen from the outside. Offering answers to the management of urban plantations, therefore, require the use of methods for internal diagnosis for trees in situ.

Modern techniques can be used to minimize the risk associated with tree failure.

There have been significant advances in decay-detection equipment, and formulas and guidelines for assessing hazardous trees (Johnstone et al. 2010a;[START_REF] Leong | A purposeful classification of tree decay detection tools[END_REF], as detailed in Appendix A.

Non-Destructive evaluation of wood

Classic methods for trees inner diagnostic are based either in drilling measures or in punctual acoustical or electrical measures [START_REF] Pellerin | Nondestructive evaluation of wood[END_REF]. Intrusive methods, such as drilling (Figure 2a), remain limited due to the incapability of obtaining information other than one dimensional (profile of drilling resistance versus drilling depth) [START_REF] Nicolotti | Application and comparison of three tomographic techniques for detection of decay in trees[END_REF]. Also, there exists the risk of propagating diseases from one tree to another by using the same tools. For punctual methods, as for intrusive, the information provided is mainly local, corresponding to the measurement location [START_REF] Drénou | Vitalité et solidité de l'arbre : choisir les méthodes de diagnostic[END_REF][START_REF] Larsson | Nondestructive detection of decay in living trees[END_REF]).

Nondestructive imaging methods allow analyzing the inner structures of trees (bi or three-dimensional reconstruction), without altering their condition (Bucur 2003a). Wood imaging can be performed using electromagnetic waves (Figure 2b): gamma and X rays (Zhu et al. 1991;Habermehl and Ridder 1992;[START_REF] Giudiceandrea | A high speed CT scanner for the sawmill[END_REF], infrared thermography [START_REF] Catena | Thermography reveals hidden tree decay[END_REF][START_REF] Catena | Overview of thermal imaging for tree assessment[END_REF], radar [START_REF] Nicolotti | Application and comparison of three tomographic techniques for detection of decay in trees[END_REF][START_REF] Hagrey | Geophysical imaging of root-zone, trunk, and moisture heterogeneity[END_REF]) and microwaves [START_REF] Kaestner | Microwave polarimetry tomography of wood[END_REF][START_REF] Hislop | Microwave radar for detection of resin defects in Pinus elliottii Engelm var elliottii[END_REF]. One alternative to the electromagnetic waves is the use of elastic waves (acoustic or ultrasonic tomography) (Bucur 2003b;Johnstone et al. 2010b), presenting as advantages a reduced cost and simplified setting up.

Imaging approaches using elastic waves consist of quantitative tomography; it can be used in a qualitative way by analysis and later segmentation of the produced image. The image is called quantitative because it is a mapping of a parameter: the acoustic velocity (related to the specific elasticity of the material). This parameter changes during the development of the tree and it is strongly influenced by the action of pathogens.

Tomography using elastic waves can be extended to a three-dimensional case by interpolation, making cuts at different heights [START_REF] Socco | Feasibility of ultrasonic tomography for nondestructive testing of decay on living trees[END_REF][START_REF] Martinis | Tomographie ultrasonore pour les arbres sur pied[END_REF][START_REF] Bucur | Ultrasonic techniques for nondestructive testing of standing trees[END_REF]).

The most known commercial devices use acoustic waves, such as Arbotom [START_REF] Rinn | Influence of wood anatomy on stress wave tomography[END_REF], Picus [START_REF] Rust | A new tomographic device for the non-destructive testing of trees[END_REF] and Fakkop [START_REF] Divos | Tree evaluation by acoustic tomography[END_REF][START_REF] Divos | Resolution of stress wave based acoustic tomography[END_REF], where a hammer is used to excite the trunk, producing waves in audible frequencies (below 20 kHz). Multiple studies have evaluated the performance of these acoustic devices, highlighting different problems: a low spatial resolution, images difficult to interpret (Figure 3), detection precision not always optimal and a relative long time of execution [START_REF] Gilbert | Picus Sonic tomography for the quantification of decay in white oak (Quercus alba) and hickory (Carya spp.). Journal of arboriculture Gillis PP (1972) Orthotropic elastic constants of wood[END_REF][START_REF] Rabe | Detection of decay in trees with stress waves and interpretation of acoustic tomograms[END_REF][START_REF] Wang | Acoustic tomography for decay detection in red oak trees[END_REF][START_REF] Wang | Effect of sensor quantity on measurement accuracy of log inner defects by using stress wave[END_REF][START_REF] Deflorio | Detection of incipient decay in tree stems with sonic tomography after wounding and fungal inoculation[END_REF]Johnstone et al. 2010a). The quality of the image obtained is a function of several factors: the frequency of transmitted wave (sound or ultrasound), the signal to noise ratio (attenuation phenomenon), the number of probes used, the physical model of propagation (isotropic, orthotropic, moisture content and density gradient) and the image reconstruction algorithm (Arciniegas et al. 2014b).

Meeting with experts

Inspecting hazardous trees with classic methods

We had the opportunity to accompany an on-site tree inspection process, to evaluate its hazardousness level. An arboricultural expert, was the operator in charge of the diagnostic requested by a departmental organization in the Herault department, France. The concerned trees were from the Platanus genus, located all along a highway in the Lodève locality.

Before starting the testing procedures, we talked with the operator to learn about the diagnostic process and his experience working on the field. First, we were stated that the diagnostic is divided into two phases: visual evaluation and specialized instrumental evaluation. Depending on the tree's local environment, pathogenic state, and patrimonial importance, the procedure could consist of both phases or only one of them. For example, those trees along the highway, defined as ornamental trees, are normally evaluated by visual means, unless they present significant signs of disease (additional tools as drilling would be used in this case). Then, there are some urban trees that present high patrimonial importance or heritage value where deeper analysis can be required; this case is where methods as acoustic and ultrasonic tomography become more valuable tools. For the visual case, methods are expected to be fast in execution; for the case using instruments, a long-term study can be allowed. There are also economic factors: specialized tools present high costs that are difficult to make a profit considering that one ornamental tree diagnostic is estimated to cost 10 euros.

From the operator experience, tomography has not a very good precision on estimating defects (as confirmed by several authors), so usually, a drilling test is also required to validate. However, drawbacks for resistograph (the most common technique)

testing include the fact of creating a hole into the tree, even if it is only 2.5 mm, that will allow fungi propagation (frequently this hole is not closed, as it is expected the tree to heal).

Tree diameters tested for this operator range between 30 and 50 cm, limited by the length of the drill.

A visual inspection was expected for our study case. To start talking about how the visual diagnostic is approached, the first question is related to how the tree response to fungi attack. Tree response to wounding and invasion by decay organisms is described by the CODIT system (Compartmentalization Of Decay In Trees) [START_REF] Shigo | Ultrasonic goniometry immersion techniques for the measurement of elastic moduli[END_REF], detailed in Appendix B. Face to wounding, tree establishes a defense system known as compartmentalization. This system is constituted by three reactions zones (walls 1, 2 and 3 for longitudinal, tangential and radial directions) that limits pathogens propagation in wood at the wounding moment, and by a barrier zone (wall 4) formed by the cambium, that protects new wood from infested wood. Gradually, depending on the severity of the injury and the reaction potential of the tree, the wound will close, thanks to the new wood created annually, covering the necrotic zone that may continue hollowing.

As commented by the operator, two main types of decomposition are searched:

those located at the base of the tree and those aerial (when the branches were cut). In the case of decomposition on the base of the tree, fungi will tend to go up and reach the center of the trunk. The type of fungi will depend on species and environmental factors. To identify the type of fungi attack, a fructification process must take place on the tree, making evident the lignivore mushroom on the trunk surface.

Visual assessment of standing trees goes through two main stages. First, dendrometry data is collected. Most frequent measures are the diameter, the height of the tree and crown volume. Data are mainly used to compute failure risk associated with wind loads. A security factor can be obtained, for example, by using a tool for structural tree analysis like TreeCalc, (arbosafe GmbH, www.treecalc.com). Then, an evaluation of the tree structural state is performed using the standard method known as VTA (Visual Tree Assessment), detailed in Appendix A, introduced by professor Mattheck [START_REF] Mattheck | Field Guide for Visual Tree Assessment (vta)[END_REF]. The objective is to determine hazardous trees using exterior symptoms.

Vitality criteria are used, as for example root problems, growth defects, crown and leaves problems (low foliage density, loss of branches), among others. Other discussed symptoms were cavities presence (for example, when branches were cut), and the presence of fungal fruiting bodies (as observed in one of the trees in front of us).

We made some testing on one of the standing Platanus tree, visually designated as defective. First, a hammer blow was used to evaluate the "sonority" of the tree. This means we looked for a different sound when comparing a defective tree to a healthy one. This procedure implies some experience from the operator: as you have compared a lot of trees, you can increase your accuracy detecting defects. Then, a resistograph from iML (PD500)

was used to evaluate tree status (Figure 4). Drilling measurements were obtained at different points, showing an important inner hollow area. The device worked fast (less than one minute for every direction), and the results could be seen directly in the digital screen of the tool. Defects position was briefly discussed. The most common pattern consists of a centered defect. At late decay stages, heart-rot corresponds to cavities within the center of the tree. Following the CODIT system (Appendix B), the strongest barrier is the wall 4, the cambium's wall, and frequently this is the only barrier that can completely stop the spread of infection by closing the wound with new wood. Finally, one important remark was that tree evaluation using instruments will not deliver a fully automatic diagnostic, this meaning that the expertise of operators will be necessary even in the case of specialized approaches as acoustic tomography.

Inspecting hazardous trees using acoustic methods

Then, we had the opportunity to discuss with an arboricultural expert, who works at the ONF (Office National des Forêts), related to his expertise in the diagnostic process for standing trees by means of acoustic methods.

When we talk about acoustic testing, the acoustic hammer appears as the first approach to be considered (Figure 5). This device is designed to measure the velocity of propagation of mechanical waves in standing trees, as this measurement is directly related to the wood dynamic modulus of elasticity and density, allowing the user to evaluate the stiffness and related properties in standing trees. Advantages of using this technology include a fast test execution and reduced affectation to tree state. However, as commented by the expert, the acoustic hammer is not reliable to evaluate some common fungi attacks.

Additionally, both defect positioning and size estimation are difficult tasks to be addressed with this device as only one information is obtained. In France, its use has been reduced in the last 10 years, being replaced by alternatives techniques. One of the recent techniques used for standing tree defects detection and localization is tomography (cross-sectional image from an inner parameter of the tree). At ONF, they are using both acoustic (Figure 6) and electrical tomography as non-destructive evaluation techniques, both devices from the same manufacturer (Picus, Argus Electronic Gmbh). For the tomography image acquisition, the process includes hitting with a hammer several times the measuring point indicated by a nail (manufacturer recommend 3 times, but expert assures by hitting 6 times); time-of-flight is estimated as the mean value. For the tree contour definition, a caliper is used (Argus Electronic provides its own solution: Picus Calliper); the geometry is defined by a triangulation method. The maximum number of sensors used by this operator is 24. Testing time could take 45 minutes to one hour. When looking for fissures, the trunk is divided into two sections, and after the two tomographies are combined.

The number of cases where tomography is required has been increasing these years.

The most common case for acoustic testing is trees with high patrimonial value. Last year, that operator performed 15 inspections. The ONF counts with 4 tomographic devices and 8 trained operators. In terms of costs, testing 4 trees in a day could reach 700 euros. Buying the device is about 25.000 euros. The arboricultural expert insisted on the fact that this device is intended to be a support in the diagnostic process, not a definitive answer. A 3D option is available, by performing tomography at different heights.

After acoustic tomography testing, the security factor is obtained using the dendrometry data and the tomographic image via the TreeCalc software. Also, the combination of Electric (Electrical Resistance Tomography) and Acoustic tomography leads to a better interpretation of the tree state (Figure 7). The electric resistance of the wood is influenced by the water content and chemical elements which change according to the status of wood. Therefore, this combination could address the questions of the type of defect (hollow, decay or crack), stage of decay (incipient or advanced) and a more detailed decay size estimation. In Colombia, there was also possible to participate to an on-site tree inspection process. In this case, the tree was in Bogotá, and the Secretary of Environment oversaw the testing. The operators used the acoustic tomography unit from Rinntech, Arbotom (Figure 8). Testing procedure was like the previous case with the ONF operator, except for the tree contour definition: no caliper was used, the sensors position was approximated using small variations from a perfect circle. They included drilling profiles using a Resistograph from Rinntech to complete the diagnosis procedure. Testing time was about one hour. The maximum number of sensors used was 24. The basic consideration is decay inside wood influences the propagation of elastic waves: the velocity decreases, and the attenuation increases. The Fermat's principle states that the first arriving wave will tend to travel along the fastest path, therefore the presence of pronounced velocity contrasts will tend to curve the rays [START_REF] Maurer | Application of nonlinear acoustic tomography for non-destructive testing of trees[END_REF].

The first approach to ultrasonic tomography of wood was presented by [START_REF] Tomikawa | Nondestructive Inspection of a Wooden Pole Using Ultrasonic Computed Tomography[END_REF]. He used propagation time of ultrasound to perform non-destructive testing of wooden poles, in tomographic configuration (fan-beam geometry, backprojection). The proposed USCT technique was found to be acceptable to identify rotten areas, however, a weak image resolution and long computation time were presented as the main drawbacks [START_REF] Tomikawa | Nondestructive Inspection of a Wooden Pole Using Ultrasonic Computed Tomography[END_REF]). Then, several scientific approaches have been presented with frequencies ranging from 22 kHz to 1 MHz [START_REF] Martinis | Tomographie ultrasonore pour les arbres sur pied[END_REF][START_REF] Brancheriau | Ultrasonic tomography of green wood using a non-parametric imaging algorithm with reflected waves[END_REF][START_REF] Brancheriau | Ultrasonic imaging of defects in standing trees: development of an automatic device for plantations[END_REF][START_REF] Lin | Application of an ultrasonic tomographic technique for detecting defects in standing trees[END_REF].

The feasibility of ultrasonic tomography for detecting wood decay in living trees has been evaluated in multiple studies, comparing the basic approach from Tomikawa with the results from other nondestructive techniques: ionizing radiation; thermal techniques; microwaves (Bucur 2003a); electric; georadar [START_REF] Nicolotti | Application and comparison of three tomographic techniques for detection of decay in trees[END_REF]; radar [START_REF] Sambuelli | Ultrasonic, electric and radar measurements for living trees assessment[END_REF]. Also, evaluations have been presented using cross-sections cuts from decayed trees [START_REF] Socco | Feasibility of ultrasonic tomography for nondestructive testing of decay on living trees[END_REF]) and drilling profiles [START_REF] Martinis | Tomographie ultrasonore pour les arbres sur pied[END_REF], where the authors stated that ultrasonic techniques are suitable for standing trees quality evaluation, allowing to find knots, decay and fungal attack. However, some common drawbacks using this technique are presented: the anisotropy effect, linked to the apparition of a bias in the image; the relation between frequency and smallest detectable defect; a low signal-to-noise ratio (without debarking the tree, it can be difficult to determine the first arrival of the signal); long data acquisition.

Variability in living trees presents multiple challenges to ultrasonic imaging.

Factors such as the grain orientation, the tree species and the moisture content could affect the ultrasonic wave propagation [START_REF] Schafer | Ultrasound for defect detection and grading in wood and lumber[END_REF]. Also, a study was carried out to evaluate the variation of ultrasonic properties in wood at meso-structural level (1 mm or less), having into account the fact that wave propagation is influenced by the anatomical structure at microscopic scale [START_REF] Feeney | Meso-Structural Considerations of Ultrasonic Propagation in Wood[END_REF]. They indicated that some structures like annual rings, early-wood, and late-wood present a significant variation in the velocity of ultrasonic waves, and how the variation in density values from the bark to the pith will affect the wave speed measurements. As well, it is expected a variation in velocity values passing from sapwood to heartwood, due to the change in density properties. Beall [START_REF] Beall | Overview of the use of ultrasonic technologies in research on wood properties[END_REF] summarized some of these effects as follows: the velocity of acoustic waves decreases with moisture content up to fiber saturation point; an increasing grain angle related to an increment of the wave attenuation and a reduction of the wave velocity; in softwoods, the growth ring angle affects the attenuation, with maximum values in angles from 45 to 60 degrees; the effect of density on wave velocity depends on the species and the propagation angle; in the presence of decay, the wave speed decreases, and the attenuation increases.

To overcome the general drawbacks in ultrasonic tomography, it is necessary to well understand and model the complexity of wood material, to develop techniques for the signal processing and the image reconstruction stages adapted to these conditions.

A crucial step to perform ultrasonic non-destructive testing of standing trees is the time-of-flight measurement, influencing the precision of defect detection. Aiming to increase the accuracy on the estimation, the characteristics of the ultrasonic chain of measurement should be adapted to the constraints of wood testing in living condition.

Considering conventional ultrasonic testing, the object is excited with a pulse, and TOF measurement relies on the estimation of the signal instantaneous power by determining the first arrival above a noise threshold, defined by the user [START_REF] Bucur | Acoustics of Wood[END_REF][START_REF] Loosvelt | A Wavelet-Based Processing method for simultaneously determining ultrasonic velocity and material thickness[END_REF]. Also, a pulse train can be used to boost the transmitted energy for a specific frequency (Brancheriau et al. 2012a). Automatic methods for detecting first arrivals have been proposed, including pickers based on the Akaike information criteria (AIC) [START_REF] Sleeman | Robust automatic P-phase picking: an on-line implementation in the analysis of broadband seismogram recordings[END_REF][START_REF] Zhang | Automatic P-Wave Arrival Detection and Picking with Multiscale Wavelet Analysis for Single-Component Recordings[END_REF]) and the Hinkley criteria [START_REF] Kurz | Strategies for reliable automatic onset time picking of acoustic emissions and of ultrasound signals in concrete[END_REF].

Alternatives include the transmission of encoded waveforms, such as the chirp-coded excitation method, where a recognizable signature is sent through the media and the TOF is estimated using a cross-correlation function [START_REF] Pedersen | Clinical evaluation of chirp-coded excitation in medical ultrasound[END_REF][START_REF] Lasaygues | Use of a Chirp-coded Excitation Method in Order to Improve Geometrical and Acoustical Measurements in Wood Specimen[END_REF][START_REF] Rouyer | The benefits of compression methods in acoustic coherence tomography[END_REF]. Bearing in mind the wide range of signals and TOF detection techniques, the choice of parameters for standing tree ultrasonic testing demands an evaluation of the methods accuracy.

The wood anisotropy has been studied considering ultrasonic testing. The anisotropic behavior leads to wavefronts propagating in the media with deformed shapes, compared to the spherical wavefronts in isotropic media, as shown in different studies, either by analytical [START_REF] Payton | Wave Fronts in Wood[END_REF], by simulation using finite differences [START_REF] Schubert | Acousto-ultrasound assessment of inner wood-decay in standing trees, possibilities and limitations[END_REF][START_REF] Schubert | Acoustic wood tomography on trees and the challenge of wood heterogeneity[END_REF], finite elements [START_REF] Sebera | Numerical simulation of elastic wave propagation in wood with defined tree rings[END_REF]) and by experimental setting [START_REF] Zhang | Experimental investigation of stress wave propagation in standing trees[END_REF][START_REF] Gao | Application of an ultrasonic wave propagation field in the quantitative identification of cavity defect of log disc[END_REF]. This affects directly the image reconstruction process, due to the curvature in the wave ray paths from the emitter to the receiver, being longer than the classic consideration of straight ray's paths. In terms of the ray theory, the trajectories of a wave traveling through an anisotropic media do not follow a straight path for all directions of propagation. This wave path curvature directly affects the image reconstruction process: classic inversion methods consider straight rays; when they are used for anisotropic materials the resulting image is blurry [START_REF] Maurer | Application of nonlinear acoustic tomography for non-destructive testing of trees[END_REF]. Therefore, modeling the ray paths for a wood cross-section is an essential step to address the inverse problem for standing tree tomography.

The ray-tracing algorithm has been widely used for numerical simulation of wave propagation for anisotropic media, particularly in the field of seismography exploration [START_REF] Cerveny | Seismic Ray Theory[END_REF][START_REF] Bóna | Fermat's principle for seismic rays in elastic media[END_REF]. The approach known as wavefront construction has been presented previously in applications for the geophysics field [START_REF] Vinje | Traveltime and amplitude estimation using wavefront construction[END_REF][START_REF] Vinje | Estimation of multivalued arrivals in 3D models using wavefront construction-Part I1[END_REF][START_REF] Vinje | 3-D ray modeling by wavefront construction in open models[END_REF][START_REF] Lambaré | Two-dimensional multivalued traveltime and amplitude maps by uniform sampling of a ray field[END_REF][START_REF] Ettrich | Wave front construction in smooth media for prestack depth migration[END_REF][START_REF] Lucio | 3D multivalued travel time and amplitude maps[END_REF][START_REF] Coman | Traveltime computation by wavefront-orientated ray tracing[END_REF][START_REF] Chambers | A practical implementation of wave front construction for 3-D isotropic media[END_REF]. This method propagates a ray-field using a constant time step, with new wavefronts constructed from the previous ones. Wave velocity is computed by the solution of the Christoffel equation using the material elastic constants.

Travel times are estimated by following individual paths from the transmitter to each receiver.

Ideally, the effect of anisotropy should be included in the image reconstruction process [START_REF] Pratt | Traveltime tomography in anisotropic media-II. Application[END_REF]. Few approaches have dealt with curved rays, and there is no known exact solution. Inversion techniques using curved rays for USCT have proved to be an efficient way to handle with strong scattering in the case of isotropic media [START_REF] Schechter | Ultrasonic tomography using curved ray paths obtained by wave propagation simulations on a massively parallel computer[END_REF]. Jackson [START_REF] Jackson | MIGRATOM: Geophysical tomography using wavefront migration and fuzzy constraints[END_REF] Then, residuals are computed by estimating the difference between the measured and the calculated travel times. These residuals are used to calculate an incremental correction factor (slowness correction) for every pixel crossed by the corresponding rays. After each iteration, velocity values for each pixel are obtained using the SIRT algorithm. With respect to the anisotropy, it is limited to an elliptical dependence of velocity (elliptic anisotropy), defined as a ratio of the velocities in the horizontal and vertical direction. The vertical coordinates are stretched by a factor equal to the anisotropy ratio, then the calculations are performed without any additional modification. Experiments performed with synthetic data showed that for the case of an unconstrained reconstruction, some spurious anomalies appeared. Then, to improve the reliability of the reconstructed images, it was necessary to include constraints in the reconstruction process, like maximum and minimum velocity values or known boundary conditions.

Maurer et al. [START_REF] Maurer | A simple anisotropy correction procedure for acoustic wood tomography[END_REF]) proposed a correction procedure to remove anisotropy effects in the image reconstruction for non-destructive testing of trees. They assumed a "weak anisotropy" simplification, by considering that this effect is governed primarily by a constant tangential-radial velocity ratio. This ratio was used to obtain an anisotropy correction factor to apply to the travel time data before using an isotropic reconstruction method. As this method does not consider the curvature of the rays, its applications were limited to the reconstruction of low contrasted defects. Liu and Li (Liu and Li 2018) presented a method for tree acoustic tomography that performs a 2D velocity map iterative reconstruction following a similar strategy to the one used by the MIGRATOM software, using a propagation model based on the Huygens' principle to find the rays' curvature due to the presence of defects (strong velocity contrast regions). Images were obtained from tree sections and standing trees using this method and comparing it to a straight-line inversion, showing an improvement in the estimation of the size of the defect.

However, the ray curvature due to the wood anisotropy has not been yet considered for the image reconstruction.

Research questions

How to increase the ultrasound image quality for nondestructive testing of living trees, by considering the wood complexity as material in the tomography image construction? This question is closely related to the hypothesis presented down below.

The first question is associated with the direct problem, this is, how elastic waves propagate in an orthotropic and heterogeneous material such as wood. The propagation of stress waves in anisotropic media can be approximated by analytical and numerical models.

When multiples factors of the material variability are considered on the model, its complexity will increase but the time-of-flight measure estimation would be closer to the real value. By developing a model that relates the inner mechanical properties of wood to the wave propagation velocity, the inverse problem can be approached.

The second scientific question is related to the algorithm for imaging reconstruction, known as the inverse problem, having into account that this algorithm must consider curved rays in an anisotropic media. Therefore, it is necessary to evaluate if there is a solution for the tomographic image reconstruction considering orthotropic behavior and therefore curved rays. Furthermore, how different is the image using this approach if we compare with respect to the classic techniques using an isotropic assumption? The image interpretation as a tool to help in the discrimination process of decay and other defects inside the trees should be adapted to the limitations of the image reconstructed and must consider the variability of trees. How precise can be the decay detection using the ultrasound tomography images obtained with the proposed method?

To answer these questions, this thesis project presents a scientific approach to perform ultrasonic imaging of standing trees in urban areas, based on the considerations for wood variability involved in each stage of the imaging process.

Methodology and organization of the thesis

The thesis was divided into four main parts described below. It was developed in the modality of joint supervision between the Universidad Nacional de Colombia (Research Group GAUNAL, Engineering Faculty, Bogotá, Colombia) and the University of Montpellier (Research Unit BioWooEB, CIRAD, Montpellier, France); the activities were developed alternately in the two countries.

Wood and ultrasonic waves: theoretical aspects

The basic considerations of the wood as a material and the governing mechanical equations to describe the propagation of ultrasonic waves are presented in Chapter 2. The Christoffel equation is presented, relating the wave propagation velocity with the material elastic constants and the wave direction of propagation. An evaluation of the sensitivity of the mechanical parameters in the computation of the compression wave velocity using the Christoffel equation was performed, using data from several wood species.

Partial results of these topics were published at:

Configuration of the ultrasonic measurement system

Time-of-flight (TOF) measurement is a critical step to perform ultrasonic nondestructive testing of standing trees. For the experimental set-up, an electro-acoustic measurement system was implemented. Aiming to increase the accuracy on the TOF estimation, the characteristics of the ultrasonic measurement chain should be adapted to the constraints of wood testing in living condition. A characterization of the ultrasonic sensors has been done. Then, the excitation signal parameters were defined through an experimental study, presented in Chapter 3.

Several configurations were tested in a standing tree, including different signal shapes with distinct frequency responses. To perform the TOF measurement automatically, three detection methods were tested. Data analysis included signal-to-noise ratio (SNR)

and root mean square voltage (RMS) measurements and time-frequency analysis. To reduce uncertainty on the TOF measurement, tests were repeated several times, and the best combination for the signal parameters and the TOF determination method was chosen.

Partial results of these topics were published at:

• Espinosa L, Bacca J, Prieto F,

Wave propagation model by a raytracing approach

In Chapter 4, a numerical study was performed to model the effects of wood anisotropy and the presence of defective areas in the propagation of elastic waves. This study allowed us to observe these effects on the time-of-flight measurements and to develop an approximated model, with low computational cost, that could be used in the inversion process. A "raytracing" model was implemented, a method based on the ray's theory, used in geophysics, that traces the ray's propagation inside the studied object. The raytracing model is based on a formulation with the Christoffel equation. Simplifying hypotheses were considered: for example, a 2-dimensional approximation of the trunk geometry. The model using ray-tracing was tested for isotropic and orthotropic media, with and without defects.

A comparison of this model to the solution obtained with a finite element method was done. Tested configurations included a healthy and a defective case. The comparison evaluated both approaches by studying the wavefronts and time-of-flight measurements.

An experimental study for laboratory validation was performed on this stage. The time-of-flight measurements were compared to the model approximations. Configurations included a healthy case, a centered defect case, and an off-centered defect case. Several defect diameters were tested.

Partial results of these topics were published or presented at:

Algebraic solution of the inverse problem for an orthotropic material

In the context of the acoustical tomography imaging of standing trees, the wave velocity values are determined for each local area (pixel of the resulting map) of the scanned cross section by solving an ill-conditioned inverse problem with a low number of acoustic measurements. Up to now, this problem has been solved by assuming that the transverse cross-section of trees is quasi-isotropic. The hypothesis of isotropy blurs the image and makes it difficult to characterize the mechanical state of wood and the presence of a defect. A way to overcome this problem is to consider the cross-section of a standing tree as being cylindrically orthotropic in the process of inversion, such that the elastic constants of wood for each pixel in the radial-tangential plane could be determined.

In Chapter 5, the aim was to develop a reconstruction method adapted to the wood anisotropy. An inversion method based on the algebraic solution was proposed and implemented. The algebraic resolution allowed the use of curved ray models as ray-tracing and is well adapted to heterogeneous and orthotropic materials. This inversion method was tested on numerical simulations and then on laboratory specimens, comparing to the images obtained with an algorithm considering straight-line rays.

Partial results of these topics were published or presented at:

• Espinosa L,

WOOD AND ULTRASONIC WAVES:

THEORETICAL ASPECTS

Introduction

The wood material is presented in this chapter from its anatomical and mechanical characteristics, drawing attention to its orthotropic condition. Then, the equations describing the motion of elastic waves for an orthotropic medium are presented, which solution is known as the Christoffel equation, linking the wave velocity to the wood mechanical parameters and the direction of propagation. For this equation, a sensitivity analysis is performed to evaluate the influence of every mechanical parameter in the wave velocity variation.

Trees and wood

For forest inventory, a tree can be considered as a woody perennial of a species presenting an elongated stem that supports branches and leaves [START_REF] Gschwantner | Common tree definitions for national forest inventories in Europe[END_REF]).

Trees cover a large portion of land ecosystems: The Global Forest Resources Assessment 2010 (FRA 2010) indicated that the world's total forest area is over 4 billion hectares, this is 31% of total land area (FAO 2010).

Anatomy

At a macroscopic level, a growing, living tree has two mains parts: the shoot and the roots. The shoot includes the trunk, the branches, and the leaves. The trunk is composed of multiple layers, from outside to inside: out bark, inner bark, vascular cambium, sapwood, heartwood, and pith (Figure 9a). The main component in the trunk is therefore wood.

Sapwood is the active wood that conducts water from the roots to the leaves while the heartwood is a zone of inactive cellular tissue. From a wood segment, three planes can be distinguished (Figure 9b): the cross or transverse section (face that is exposed when a tree is cut down), the radial section (pith-to-bark direction) and the tangential section (perpendicular to the radial plane). Associated to these planes, three mutually perpendicular axes are defined: the longitudinal axis (L) parallel to the fibers, the radial axis (R) orientated from the bark to the pith and the tangential axis (T), tangent to the growth rings and perpendicular to the grain.

Trees are mainly classified in softwoods and hardwoods. Softwoods come from gymnosperms (essentially conifers) and hardwoods come from angiosperms (flowering plants). Softwoods are generally needle-leaved evergreen trees (spruces and pines, for example), while hardwoods are commonly broad-leaf trees (oaks and maples, for example).

The fundamental differences between woods are related to the characteristics at the microscopic level. Hardwoods present a more complicated microscopic structure than softwoods (Figure 9c): softwoods are mainly composed of axial tracheids, aligned parallel with the trunk, and rays, in the radial direction; for hardwoods, the axial system is composed of fibers, vessels, and axial parenchyma, the radial direction is composed of rays. The cell wall possesses a structure that will establish many wood mechanical and physical properties. It is composed of three chemical components: cellulose, lignin, and (a) Ek et al. 2016). The cell wall is divided into three main regions: the middle lamella, the primary wall and the secondary wall (Figure 10). The main function of the middle lamella is to provide adhesion between the wood cells. The primary wall is characterized by randomly orientated cellulose microfibrils; this layer is thin and often is difficult to differentiate from middle lamella. The secondary cell wall is composed of three layers, S1, S2, and S3, being S1 and S3 thin compared to the central S2 layer. The first layer S1 is adjacent to the primary wall and is characterized by a large microfibril angle (50° to 70°), with respect to the axial direction. The second layer possesses more importance in determining the wood characteristics, due to their thickness and low microfibril angle (5° to 30°). S3 layer is thin and presents a high microfibril angle (> 70°) as the S1 layer, but the lignin concentration is lower. The orientation of the cellulose microfibrils of the S2 layer is highly relevant for determining the physical properties of the wood cell and therefore the wood structure.

(b) (c) R T L hemicellulose (

Mechanical behavior

The cellular structure of wood and the disposition of cellulose in the cell wall explains the anisotropic property of wood [START_REF] Gibson | Cellular Solids: Structure and Properties[END_REF]. At a macroscopic level, wood can be considered as an orthotropic and heterogeneous material. In orthotropic materials, the mechanical properties will change depending on three mutually perpendicular axes. For wood, these are the longitudinal axis (L), the radial axis (R) and the tangential axis (T).

Considering wood as an elastic medium, the constitutive law that relates the applied forces with the associated deformations is the Hook law, defined in form of stress (σ) and strain (ϵ) as:

𝜎 = 𝐶𝜖, ( 1
)
where C is known as the rigidity matrix. The stress vector describes the stress conditions at any point x in the space. For small deformations, the strain vector can be expressed in terms of the displacement vector (u) as:

𝜖 𝑖𝑗 = 1 2 ( ∂𝑢 𝑖 ∂𝑥 𝑗 + ∂𝑢 𝑗 ∂𝑥 𝑖
) , with 𝑖, 𝑗 = 1,2,3.

(2) Nine elastic constants define the linear relation between the stress vector σ and the strain vector s for orthotropic materials [START_REF] Berthelot | Matériaux composites[END_REF]:

( 𝜎 𝑅𝑅 𝜎 𝑇𝑇 𝜎 𝐿𝐿 𝜎 𝑇𝐿 𝜎 𝑅𝐿 𝜎 𝑅𝑇 ) = ( 𝐶 11 𝐶 12 𝐶 13 0 0 0 𝐶 21 𝐶 22 𝐶 23 0 0 0 𝐶 31 𝐶 32 𝐶 33 0 0 0 0 0 0 𝐶 44 0 0 0 0 0 0 𝐶 55 0 0 0 0 0 0 𝐶 66 ) ( 𝜖 𝑅𝑅 𝜖 𝑇𝑇 𝜖 𝐿𝐿 2𝜖 𝑇𝐿 2𝜖 𝑅𝐿 2𝜖 𝑅𝑇 ) . ( 3
)
With the symmetry C12 = C21, C13 = C31 and C23 = C32.

The elastic relations are expressed as a function of the elasticity moduli, that can be determined by mechanical tests. Using the engineering constants notation, the inverse relation using the compliance matrix S (s = Sσ = C -1 σ) is presented as:

( 𝜖 𝑅𝑅 𝜖 𝑇𝑇 𝜖 𝐿𝐿 2𝜖 𝑇𝐿 2𝜖 𝑅𝐿 2𝜖 𝑅𝑇 ) = ( 1 𝐸 𝑅 - 𝜈 𝑇𝑅 𝐸 𝑇 - 𝜈 𝐿𝑅 𝐸 𝐿 0 0 0 - 𝜈 𝑅𝑇 𝐸 𝑅 1 𝐸 𝑇 - 𝜈 𝐿𝑇 𝐸 𝐿 0 0 0 - 𝜈 𝑅𝐿 𝐸 𝑅 - 𝜈 𝑇𝐿 𝐸 𝑇 1 𝐸 𝐿 0 0 0 0 0 0 1 𝐺 𝑇𝐿 0 0 0 0 0 0 1 𝐺 𝐿𝑅 0 0 0 0 0 0 1 𝐺 𝑅𝑇 ) ( 𝜎 𝑅𝑅 𝜎 𝑇𝑇 𝜎 𝐿𝐿 𝜎 𝑇𝐿 𝜎 𝑅𝐿 𝜎 𝑅𝑇 ) . ( 4
)
Where the parameters ER, ET and EL correspond to the Young modulus in each direction (from a tensile or compressive test), the parameters νRT, νTR, νRL, νLR, νTL, and νLT correspond to the Poison's ratios (ratio of the transverse to axial strain) and finally GTL, GLR, and GRT correspond to the shear modulus (from a shear test). To illustrate the effect of anisotropy in wood, Table 1 presents reference values for the Young modulus ratio in the three directions [START_REF] Guitard | Mécanique du matériau bois et composites. Cepadues Editions Habermehl A, Ridder HW[END_REF]. For all known species, the order of magnitude for these moduli are EL >> ER > ET. The mechanical behavior of wood is linked to its physical and chemical state. For example, in the longitudinal direction, the fibers constitute a strengthening element, resulting in the highest Young modulus value. In the radial direction, the appearance of rays constitutes the strengthening factor with respect to the tangential. In the case of density, the value will depend on multiple factors such as the cell wall density, the porosity, and moisture content. The moisture content in the wood also affects the Young modulus values [START_REF] Guitard | Mécanique du matériau bois et composites. Cepadues Editions Habermehl A, Ridder HW[END_REF].

Wave propagation in anisotropic media

The fundamental relation between spatial variations of the stress vector (σ) and the time variations of the displacement vector (u) is known as the elastodynamic equation. It is defined as follows:

∂𝜎 𝑖𝑗 ∂𝑥 𝑗 + 𝑓 𝑖 = 𝜌 ∂ 2 𝑢 𝑖 ∂𝑡 2 , 𝑖, 𝑗 = 1,2,3.
(5)

Where fi represents the external forces acting on the body in Cartesian coordinates.

This equation relates stress acting on a point in a solid to the motion of the particles in the solid. Combining the elastodynamic equation (Equation 5) and the Hook's law (Equation 1) for a continuous medium, we obtain:

∂ ∂𝑥 𝑗 [𝐶 𝑖𝑗𝑘𝑙 𝜖 𝑘𝑙 ] + 𝑓 𝑖 = 𝜌 ∂ 2 𝑢 𝑖 ∂𝑡 2 , 𝑖, 𝑗, 𝑘, 𝑙 = 1,2,3. (6)
Combining the last two equations, and considering that exterior forces to be zero, the resulting expression is:

𝐶 𝑖𝑗𝑘𝑙 ∂ 2 𝑢 𝑘 ∂𝑥 𝑗 ∂𝑥 𝑙 = 𝜌 ∂ 2 𝑢 𝑖 ∂𝑡 2 , 𝑖, 𝑗, 𝑘, 𝑙 = 1,2,3. (7)
As we are interested in the displacements in the cross-section of a trunk, the expression is simplified to the 2-dimensional case in the plane 1,2 (or RT), written as:

A solution to this equation is presented in the form of a plane wave:

𝑢(𝑥 → , 𝑡) = 𝑢 0 𝑝 → 𝑒 𝑖(𝜔𝑡-𝑘 → 𝑥 → ) , (9)
where 𝑢 0 is the wave amplitude, 𝑝 ⃗ is the polarization vector, 𝜔 is the angular frequency and 𝑘 ⃗⃗ is the wave vector. The wave vector is related to the wave number as 𝑘 ⃗⃗ = 𝑘𝑛 ⃗⃗, where 𝑛 ⃗⃗ is the wave front propagation direction, and the wave number 𝑘 is equal to 𝑘 = 𝜔/𝑉, with 𝑉 the wave phase velocity. If the polarization vector is parallel to the wave front normal direction, the wave is known as compression wave or longitudinal wave (not to be confused with the longitudinal axis). In the case of perpendicularity between the two vectors, the wave is called shear (or transverse). Waves representation for both cases is presented in Figure 12. If the angle formed is different to the two previous cases, the wave is said to propagate in a quasi-mode (quasi-longitudinal, quasi-transverse). Using this wave expression in the Equation 8it leads to the equation system:

known as the Christoffel equation [START_REF] Royer | Elastic waves in solids[END_REF]. To obtain the phase velocity value, the equation can be solved as an eigenvalue problem:

det ([ Γ 11 Γ 12 Γ 12 Γ 22 ] -𝜌𝑉 2 [ 1 0 0 1 ]) = 0, ( 11
)
where Γ is known as the Christoffel tensor. Therefore, with the values of the rigidity matrix for a given orthotropic material the equation can be solved to compute the phase velocities with a wave normal 𝑛 ⃗⃗. The eigenvector associated to the problem corresponds to the polarization vector 𝑝 ⃗. Then, phase velocity in a cross-section (plane RT) is a function of the mechanical properties of the material (ER, ET, GRT, νRT, ρ). We are only interested in the longitudinal and quasi-longitudinal waves, since these waves travel faster than shear waves, and therefore are the first to arrive to the ultrasonic receivers.

The velocity is a function of the wavefront normal direction, that can be simply described by an angle 𝜃 between the wave direction vector (𝑛 ⃗⃗) and the direction of the radial axis (Figure 13). This is:

𝒏 𝒙 = 𝒄𝒐𝒔(𝜽), 𝒏 𝒚 = 𝒔𝒊𝒏(𝜽). ( 12
)
Figure 13: RT axis and wave normal direction to define the θ angle.

In that case, solving Equation 11 results in an expression for the phase velocity (compression wave), as follows:

𝑉 = √ Γ 11 +Γ 22 +√(Γ 22 -Γ 11 ) 2 +4⋅Γ 12 2 2𝜌 . (13)
The Christoffel coefficients Γ11, Γ22, and Γ12 are a function of the elements of the rigidity matrix Cij for an orthotropic material, computed using the elastic constants of the RT plane ER, ET, GRT, and νRT. The direction of propagation corresponds to the angle θ between the vector normal to the wavefront and the radial direction, as:

𝛤 11 = 𝐶 11 𝑐𝑜𝑠 2 𝜃 + 𝐶 66 𝑠𝑖𝑛 2 𝜃.

(

) 14
𝛤 22 = 𝐶 66 𝑐𝑜𝑠 2 𝜃 + 𝐶 22 𝑠𝑖𝑛 2 𝜃.

(

) 15
𝛤 12 = (𝐶 12 + 𝐶 66 ) 𝑐𝑜𝑠 𝜃 𝑠𝑖𝑛 𝜃.

(16)

Sensitivity analysis of the Christoffel equation

The Christoffel equation leads to a solution in the form of plane waves, relating the propagation velocity with the material elastic constants and the wave direction of propagation. Therefore, this relationship can be used to pass from a set of measured velocities to an estimation of the mechanical parameters, a procedure known as inverse problem [START_REF] Bucur | Elastic constants for wood by an ultrasonic method[END_REF][START_REF] Castagnede | Optimized Determination of Elastic Constants of Anisotropic Solids from Wavespeed Measurements[END_REF][START_REF] Bucur | Acoustics of Wood[END_REF][START_REF] Dahmen | Elastic constants measurement of anisotropic Olivier wood plates using air-coupled transducers generated Lamb wave and ultrasonic bulk wave[END_REF][START_REF] Longo | Wood elastic characterization from a single sample by resonant ultrasound spectroscopy[END_REF][START_REF] Gonçalves | Elastic constants of wood determined by ultrasound using three geometries of specimens[END_REF][START_REF] Alves | Application of Acoustic Tomography and Ultrasonic Waves to Estimate Stiffness Constants of Muiracatiara Brazilian Wood[END_REF]. For instance, ultrasonic goniometry relies on the principle of inverse problem to determine the elastic constants of the stiffness matrix for the characterization of different composite materials (Siva Shashidhara Reddy et al. 2005;Zhao et al. 2016), including wood [START_REF] Preziosa | Méthode de détermination des constantes élastiques du matériau bois par utilisation des ultrasons Qian S, Chen D[END_REF].

The mentioned procedure works by using small samples with cubic, prismatic, polyhedral geometry, or multifaceted discs, to determine the 9 elastic constants of wood, by assuming the homogeneity of the specimens, and thus that the mechanical properties are constants within the specimens.

Here we examined the sensitivity of the mechanical parameters in the computation of the compression wave velocity using the Christoffel equation. Published data were used to study the effect of anisotropy, according to the orientation of the wavefront relative to the radial-tangential plane. Additionally, fluctuations in each value of mechanical parameters were introduced in the velocity computation according to the orientation of the wavefront. The results made it possible to examine the consequences of proposing various simplified hypotheses based on an inversion process for standing tree tomography.

Mechanical parameters and sensitivity equations

Table 2 presents the tree species selected from previously published data [START_REF] Ross | Wood handbook: wood as an engineering material[END_REF]). The species were chosen to cover a wide range of transverse anisotropy ratio, mechanical parameters, and density. ER/ET is the anisotropy ratio between the stiffnesses in the radial-tangential directions. The variation range of the anisotropy ratio is from 1.36

for Douglas fir to 2.30 for Sweetgum. The species in A set of angles θ ranging from 0° (radial direction) to 90° (tangential direction)

were used in Equation 13to evaluate the influence of the selected species (Table 2) on the wave velocity. The velocity values for each species were then computed by introducing a variation of ±10% for the mechanical parameters (ER, ET, GRT, and νRT) to evaluate the sensitivity of these parameters in Equation 13. For example, Young's modulus in the radial direction was increased as ERsup = 1.1 x ER and decreased as ERinf = 0.9 x ER. Variations such as 4% and 6% on the velocity measurement for wood testing, within the same species, has been previously reported [START_REF] Bucur | Acoustics of Wood[END_REF][START_REF] Chauhan | Assessment of variability in morphological and wood quality traits in Melia dubia Cav. for selection of superior trees[END_REF]. Thus, the influence on the velocity for a variation in the mechanical parameters of 10% is noteworthy.

The corresponding modified velocities values were named Vsup (velocity when the parameters were increased by 10%) and Vinf (velocity when the parameters were decreased by 10%). The variation of velocity (in percentage) for each parameter was obtained as follows:

%𝑉 𝑠𝑢𝑝 = 𝑉 𝑠𝑢𝑝 -𝑉 𝑉 * 100. ( 17
)
%𝑉 𝑖𝑛𝑓 = 𝑉-𝑉 𝑖𝑛𝑓 𝑉 * 100.

(

) %𝑉 = 𝑉 𝑠𝑢𝑝 +𝑉 𝑖𝑛𝑓 2 . ( 18
) 19

Results and discussion

Figure 14 presents the wave velocity values depending on the angle θ for the selected species. The relationship between the wave velocity and the direction of propagation (angle θ), is a direct consequence of the wood anisotropy in the RT plane. For all species, higher velocities were obtained in the radial direction (θ = 0°) since this direction is stiffer than the tangential direction. The anisotropy between ER and ET can be related to the cellular microstructure of wood, which consists mainly of hollow tubular cells leading to an approximated honeycomb structure (Gillis 1972; [START_REF] Kahle | The influence of cell geometry on the elasticity of softwood[END_REF][START_REF] Gibson | Cellular Solids: Structure and Properties[END_REF]. From this approach, several aspects have been linked to the anisotropic behavior of wood. First, the effect of the cell geometry: the cell walls are highly aligned in the radial direction, while the tangential direction follows an irregular pattern [START_REF] Kahle | The influence of cell geometry on the elasticity of softwood[END_REF], resulting in a higher Young's modulus in the radial direction. Second, the mechanical properties change within the annual growth rings. Earlywood exhibits a marked anisotropic behavior (with large thin-walled cells aligned in the radial direction), while latewood exhibits a roughly isotropic behavior (with smaller and thicker-walled cells) [START_REF] Boutelje | The Relationship of Structure to Transverse Anisotropy in Wood with Reference to Shrinkage and Elasticity[END_REF]. Therefore, the proportion between earlywood and latewood affects the relationship between the radial-tangential moduli. Third, the presence of the ray cells reinforces the radial direction, depending on the width of the rays, height, and area fraction [START_REF] Burgert | The Influence of Rays on the Transverse Elastic Anisotropy in Green Wood of Deciduous Trees[END_REF]. For hardwoods, an additional factor to be considered is the vessel distribution, with higher ER/ET ratios for diffuse porous species, such as Yellow poplar, than ring-porous species, such as Oak [START_REF] Beery | Quantitative Wood Anatomy-Relating Anatomy to Transverse Tensile Strength[END_REF].

Differences between maximum and minimum velocities were higher for hardwoods, as they presented a higher anisotropy ratio (difference of 596 m/s for Sweetgum). In contrast, lower velocity differences were found in the trees with a lower anisotropy ratio (270 m/s for Douglas fir). For hardwoods, the minimum values of velocity were found in the pure tangential direction (θ = 90°). For softwoods, the minimum values were not found directly in the tangential direction, but in an angle ranging from 55° to 66°

(Figure 14 and Table 3). This angle depends mainly on GRT and νRT parameters, as they affect the off-diagonal parameter Γ12 in the Christoffel's equation. The Γ12 coefficient has a higher influence on the velocity computation when the terms Γ11 and Γ22 are equal (the term inside the inner root square will only depend on Γ12), which occurred for angles ranging from 50° to 60°.

Each mechanical parameter was increased and decreased by 10%, and the velocity values were computed using Equation 13. Table 4 displays the maximums of the velocity variations (Equation 19) after changing each parameter. For example, the variations of velocity values (Equations 17 and 18) for Sitka spruce are shown in Figure 15. The velocity variation increased as the angle approached 0° when ER was altered (maximum variation of circa 4%, Table 4). On the contrary, the variation was at its maximum when the angle reached 90° for ET (maximum variation of 6%, Table 4). This was explained by the fact that the C11 element is predominant in Eq. 1 (axis 1 is associated with the radial direction), when the angle θ approaches to 0°. The same reasoning can be applied to the element C22

when the angle θ tends to 90° (tangential direction). GRT presented the lowest variation with a maximum of 0.34%, at an angle of 53° for Spruce in Figure 15 (overall variation of 0.8%, Table 4). The maximum variation was found in angles ranging from 49° to 58° when GRT and νRT were changed to be almost equal for both parameters in each species (variation of 2.7% for νRT, Table 4). These angles increased as the ratio of anisotropy (ER/ET) increased (49° for Douglas fir and 58° for Sweetgum). This phenomenon was already explained by the effect of the off-diagonal parameter Γ12 in the Christoffel's equation. The velocity variations for changes in ER, ET, and νRT did not reach zero as it did for GRT. This was explained by the fact that the velocity is computed using the stiffness constants Cij, which was modified by ER, ET, and νRT. Only the stiffness constant C66 was affected solely by GRT. Velocity values were more affected by the ET and ER parameters than by νRT and GRT. The order of influence, from biggest to smallest was ET, ER, νRT, and GRT with a maximum variation of 5.9%, 4.1%, 2.7%, and 0.8%, respectively. Figure 16 shows the variation of velocity %V divided by a variation λ in the mechanical parameters, with λ ranging between 10% to 50%. Even when these variations were not linear for ET, ER, and νRT, the conclusions from the sensitivity study using λ=10% remain valid as the variations were monotonous.

Figure 16: Velocity variation divided by the variation λ in the mechanical parameters, with λ ranging between 10% to 50%.

Considering the inverse problem in the case of standing tree acoustical tomography imaging, when passing from the velocity value to the mechanical parameters, the initial problem counts for 5 parameters (4 elastic parameters and the density) associated with each pixel of the tomogram. Bearing in mind the low sensitivity of νRT and GRT in Equation 13, these two parameters would be determined with low precision. A first approximation would be to set these two parameters to zero to find an initial solution, only for the two Young's moduli, and then to use this solution to attempt again the inversion, but this time with all variables. To establish the velocity variation, in this case, the νRT and GRT parameters were set to zero, and the corresponding velocity was compared to the velocity obtained using all the parameters. Figure 17 presents the variations obtained for Spruce and Oak species. As expected, a higher variation of velocity was obtained when νRT was nil compared to GRT.

When both parameters were nil, the variation of velocity was maximized for angles rating between 50° to 60°. However, even when the sensitivity of νRT and GRT on the velocity computation was low compared to Young's moduli, the maximum velocity variation was circa -20% for Spruce and -30% for Northern red oak. Thus, it was concluded that νRT and GRT cannot be neglected, even in first approximation, for Equation 13. When considering the anisotropy of wood in the transverse cross section of a tree, the propagation paths of acoustic waves are curved, and not straight rays as they are for an isotropic material (Espinosa et al. 2017b). As a result, the notion of wave velocity (considered as an intrinsic parameter of the material) associated with one pixel of the tomogram has no physical sense for anisotropic material, because the velocity depends on the direction of propagation. The only intrinsic parameters to be considered should be the elastic parameters and the density (5 intrinsic parameters). Under consideration that the knowledge on the specific stiffness (stiffness matrix divided by the density) is enough to allow for a tree health assessment, the number of unknowns for each pixel is reduced to four. In this case, the inversion process will lead to 4 tomograms associated with the 4 parameters.

Synthesis

The anisotropy of wood in the radial-tangential plane directly influences wave velocity depending on the direction of propagation. The evolution of wave velocity according to the direction of propagation depends on the considered species, with a difference between softwoods and hardwoods. The radial direction, θ = 0°, corresponded to the fastest wave velocity. The shear modulus and Poisson's ratio determined the angle for the minimum velocity of softwoods, ranging from 55° to 66°. For hardwoods, the minimum velocity was in the tangential direction (θ = 90°). From the sensitivity analysis of the Christoffel equation, it was found that the order of influence of the mechanical parameters on the velocity variation, from largest to smallest was: ET, ER, νRT, GRT.

Considering an initial variation of 10% for each parameter, the maximum of the resulting velocity variations was 7 times higher for ET than for GRT. Young's moduli influence was maximized when the direction of propagation was close to the tangential or radial axis.

Poisson's ratio and shear modulus influences were maximized in directions ranging from 50° to 60°. Even if the influence of the Poisson's ratio and shear modulus was low, the νRT and GRT parameters cannot be neglected in the Christoffel equation to solve the inverse problem of standing tree tomography.

CONFIGURATION OF THE ULTRASONIC

MEASUREMENT SYSTEM

Introduction

The aim of this chapter was to compare several signal shapes and TOF detection methods, for setting up an ultrasonic chain of measurement to perform a nondestructive evaluation of standing trees. Impulsive and encoded signals were tested, combined with three different methods for TOF estimation: Threshold, AIC method, and cross-correlation.

First, the experimental setting is presented, including an electrical specification for the ultrasonic chain, the excitation signal parameters, and a description of the TOF detection methods. Then, energy and signal-to-noise ratios are computed for all configurations. A time-frequency analysis using the Gabor transform is performed, aiming to inspect energy distribution. Lastly, wave transit times are reported, computing dispersion among experiments repetition, to establish which setting leads to the highest accuracy.

Methodology

A standing plane tree (Platanus × acerifolia (Aiton) Willd) was tested (Figure 18).

Probes distance above the ground was 120 cm. The trunk diameter was 23 cm, with a regular cross-section. Tests were conducted in the dormancy period (winter). Two ultrasonic pairs of sensors were used: Physical Acoustics Corporation R3α and R6α. Sensor R3α has a main resonant frequency at 36 kHz and two secondary resonant frequencies at 22 kHz and 95 kHz; operating frequency range indicated by the manufacturer is from 25 to 70 kHz. Sensor R6α has a main resonant frequency at 60 kHz and two secondary resonant frequencies at 37 kHz and 97 kHz; operating frequency range indicated by the manufacturer is from 35 to 100 kHz. These sensors are intended for general purpose ultrasonic testing, presenting a solid stainless-steel body with a flat ceramic face. A fluid couplant was used.

The position of the sensor acting as the transmitter was fixed; the receiver position changed in 4 equidistant points at angles found along half the trunk circumference: 45°, 90°, 135°, and 180°.

Ultrasonic measurements

Ultrasonic chain of measurement is presented in Figure 19. Electrical signal This chain of measurement acts as a continuous linear stationary causal filter, then the input signal s(t) and the output signal y(t) are related by a convolution function:

𝑦(𝑡) = ((ℎ 𝑡 * * 𝑠) * ℎ 𝑚 )(𝑡), ( 20
)
where hm is the response of the tree, s(t) is the electrically generated signal, and ht*(t) is the equivalent electro-acoustic pulse response. The electro-acoustic pulse response ht*(t) is the auto-convolution of the transducers impulse response ht(t), including the response of the amplifier, and considering the transmitter and receiver transducers responses with coupling to be identical.

The five signal shapes tested were an impulse (short duration rectangular pulse), pulse train, Gaussian pulse, half-Gaussian pulse and chirp (Figure 20). For every sensor position and signal shape, ultrasonic measurement was repeated 10 times, removing and replacing the transducers. For the signal amplitude measurements, the root mean square voltage (RMS) and the signal-to-noise ratio (SNR) were computed.

RMS voltage was obtained as:

𝑉 𝑅𝑀𝑆 (𝑦) = √ 1 𝑁 ∑ |𝑦 𝑛 | 2 𝑁 𝑛=1 , (21)
with N as the signal length. SNR was computed as:

𝑆𝑁𝑅(𝑦) = 20log 𝑉 𝑅𝑀𝑆 (𝑦) 𝑉 𝑅𝑀𝑆 (𝜂) (22)
where η is the noise, estimated by selecting the first signal portion before the arrival time.

Time-of-flight detection methods

Threshold

The threshold level for the received signal had to be defined above the noise level (Arciniegas et al. 2014a). The threshold level is defined to be m times the standard deviation of the noise, with m as a user-defined parameter. For the experiments, this value was constant and fixed by trial and error to 8. TOF is then selected to be the first time point where the signal is above the threshold level.

AIC method

This method assumes that the signal can be divided into two local stationary segments, before and after the onset time, each one modeled as an autoregressive process.

The time instant where the Akaike information criterion (AIC) is minimized, corresponds to the optimal separation between noise and signal, this is, the onset time (Brancheriau et al. 2012a). For a signal divided at point k into two segments y1 (before k) and y2 (after k), the AIC criterion is computed as:

𝐴𝐼𝐶[𝑘] = 𝑘log(𝜎 2 (𝑦 1 )) + (𝑁 -𝑘)log(𝜎 2 (𝑦 2 )), (23)
where σ 2 corresponds to the variance. TOF value is obtained by founding the time point where the AIC criterion reaches the global minimum.

Cross-correlation

When a recognizable signature is sent through the media, such as chirp signal, input and output signals delay time can be obtained using cross-correlation [START_REF] Pedersen | Clinical evaluation of chirp-coded excitation in medical ultrasound[END_REF][START_REF] Lasaygues | Use of a Chirp-coded Excitation Method in Order to Improve Geometrical and Acoustical Measurements in Wood Specimen[END_REF][START_REF] Rouyer | The benefits of compression methods in acoustic coherence tomography[END_REF]. The maximum value for the cross-correlation function between two signals indicates their delay time. Normalized cross-correlation function is:

𝑟 𝑠𝑦 [𝑙] = 1 √ 𝐸 𝑠 𝐸 𝑦 ∑ 𝑠[𝑘]𝑦[𝑘 -𝑙] 𝑁 𝑘=0 (24)
where Es and Ey correspond to the signal's energy and N is the signal length. consider that the transducer impulse response will change the signal applied to the tree.

Results

Signal amplitude measurement

Using the chirp signal resulted in a lower RMS ratio for both sensors, and signals such as the half Gaussian pulse and the impulse resulted in the larger ratios.

Time-frequency analysis

As the frequency contents of the received signals varied over time, we used a timefrequency analysis to obtain a representation of the input and output signals behavior for the ultrasonic chain of measurement. From several alternatives to perform the timefrequency analysis, the Gabor transform was used [START_REF] Carmona | Practical Time-Frequency Analysis[END_REF]Qian and Chen 1999). For this study, resolution in time was set to 0.1 ms and resolution in frequency was set to 5 kHz. The receiver angle selected for the analysis was 135°, considering it presents the most energetic signals, with higher SNR ratios. Chirp is the only signal able to concentrate the energy around the central frequency for both sensors on the output signal. Gaussian pulse presented power concentration at frequencies near to the excitation central frequencies only for sensor R3α; mean power frequencies did not correspond for sensor R6α where energy dissipated at different frequencies from 60 kHz (mainly 37 kHz and 97 kHz). The other signals presented energy concentration mainly on the other sensor resonant peaks: for R3α at the third resonant peak (95 kHz), and for R6α in first and third resonant peaks (37 kHz and 97 kHz).

TOF determination

Time-of-flight was obtained for all the experiment configurations, using the Threshold and AIC method. Cross-correlation was used exclusively for the chirp signal,

given that is the only excitation signal with a similar shape on the output for both sensors, and therefore, chirp signal results are studied separately. For the sensor R3α, Figure 25 shows the mean and standard deviation values for the TOF estimated using the threshold Error bars present ±σ [START_REF] Espinosa | Accuracy on the Time-of-Flight Estimation for Ultrasonic Waves Applied to Non-Destructive Evaluation of Standing Trees: A Comparative Experimental Study[END_REF].

To obtain a clearer view of variability on the mean TOF estimation, Figure 26 presents the relative standard deviation (coefficient of variation), computed as the standard deviation divided by the corresponding mean value, and presented as a percentage. TOF values for chirp signal were obtained using the three detection methods, including cross-correlation. Figure 29 presents the mean and standard deviation values for both sensors. Mean TOF values for R3α ranged between 85 μs to 152 μs using crosscorrelation and 120 μs to 160 μs for the other two methods; for R6α ranged between 94 μs to 150 μs with cross-correlation and 90 μs to 150 μs with the other two methods. The standard deviation for R3α ranged between 0.48 μs to 0.79 μs using cross-correlation and 5.7 μs to 33 μs for AIC and threshold methods; for R6α ranged between 0.31 μs to 3.69 μs using cross-correlation and 3.34 μs to 19 μs for AIC and threshold methods. Chirp signal presents small amplitude variations at the beginning, an ill-favored condition when using AIC and Threshold methods, where a first energetic arrival is expected; therefore, the method presenting less variation is the cross-correlation method. Figure 30 presents the relative standard deviation values, where the large difference for cross-correlation compared to the other two methods is clearly observed: for R3α sensor the coefficient of variation using cross-correlation was smaller than 1% while for the other two methods ranged between 3.8% to 27%; similarly for R6α, using cross-correlation resulted in a coefficient of variation ranging between 0.2% to 3.9% compared to a range going from 3% to 12.7% for AIC and threshold methods.

Discussion

Signal energy received in angle 45° was significantly lower than those obtained for the other angles, even if this position implies the shorter distance between transmitter and receiver tested. The transmitter placed at 135° resulted generally in the larger signal energy received. Ultrasonic beams for these sensors are affected by the transducer directivity pattern, resulting in higher radiation intensity in the frontal direction of the sensor, that is orientated in the radial direction in the experiments. Another effect is related to the propagation of waves in wood: wood anisotropy affects wave propagation, including curvature of ray paths from the transmitter to receivers, with respect to straight line paths for an isotropic case [START_REF] Schubert | Acoustic wood tomography on trees and the challenge of wood heterogeneity[END_REF][START_REF] Gao | Application of an ultrasonic wave propagation field in the quantitative identification of cavity defect of log disc[END_REF].

Signals with an initial impulsive response (impulse, pulse train, and half-Gaussian pulse), resulted in larger energy received, but this energy was spread over several frequency bands, as seen on the time-frequency analysis, where the only signal able to concentrate the energy around the sensor central frequency was the chirp, the same one that presented a lower received energy. So, the compromise implies higher received energy but widely spread frequency spectrum or lower received energy but well-concentrated frequency spectrum.

Threshold and Akaike methods for TOF detection presented highly similar results, as observed in a previous study [START_REF] Arciniegas | Tomography in standing trees: revisiting the determination of acoustic wave velocity[END_REF], where it was shown that those two methods performed in agreement when the received signals presented SNR ratios above 20 dB. However, Akaike method presents as advantage that it does not need user-defined parameters, like the m value in threshold case (paragraph 3.2.2), which variation will result in a different TOF estimation. Inaccuracy increases using the AIC method when the SNR is very low, i.e. below 10 dB.

For the chirp signal, the method that presented the lower variations was the crosscorrelation. Among the other signals, the combination AIC-Impulse presented the best results. Figure 31 presents the comparison between the relative standard deviation values, for the Impulse-AIC setting and the chirp-cross-correlation. For most cases, the chirp-crosscorrelation setting resulted in lower variation for TOF estimation. The only case where chirp-cross-correlation combination was inferior to Impulse-AIC corresponded to the sensor R6α located at 45°. In that case, the signal-to-noise ratio was the lower for all configurations, near to 10 dB, while impulse presented an SNR with a mean value of 25 dB. When comparing the difference between the TOF mean values obtained with the R3α and R6α sensors, the AIC-Impulse combination resulted in a lower difference, as presented in Table 9. A dispersion effect became noticeable when using the chirp signal, that could affect the TOF measurements. When the medium is dispersive, wave propagation velocity depends on the frequency, resulting in an output signal that spreads out in time. To visualize this effect, the peaks of the Gabor transform were obtained for both input and output chirp signals, giving an idea of instantaneous frequency for different time instants, as shown in Figure 32 for the case of the sensor R6α located at 135°. Input frequencies present a linear distribution on time, however, the instantaneous output frequencies delayed more for higher frequencies.

Synthesis

Several factors influence the accuracy on the time-of-flight determination: the excitation signal characteristics in energy and frequency, the transducer frequency response, the wood inner variability, the coupling between the sensor and the tree including the bark influence, the effect of the SNR on the TOF estimation, among others. In situ testing was performed comparing five different excitation signals, two different transducers with resonant frequencies at 36 kHz and 60 kHz, 4 different receiver positions around the

WAVE PROPAGATION MODEL BY A

RAYTRACING APPROACH

Introduction

This chapter describes the influence of wood anisotropy condition on the ultrasonic wave propagation, by time-of-flight estimation using the raytracing approach. The wavefront construction method was adapted to the cross-section tree characteristics.

Healthy and defective cases were simulated. Circular defects were tested in centered and off-centered positions. The raytracing approach was compared with a FEM model, contrasting wavefronts and TOF estimations. Experimental validation was performed, using oak and ash trunk sections, simulating defects by drilling circular holes. A disk from isotropic material was used for comparison. Time-of-flight measurements were obtained for all the configurations and for simulated and experimental data. The effect on the tomographic image reconstruction is discussed.

Raytracing modeling

Method description

A circular array of transducers surrounding the tree was considered. The imaging plane (the plane of the array) is orthogonal to the longitudinal axis of the trees considering, therefore, a 2D-problem. This is a conventional hypothesis for USCT (Bucur 2003b), which in practice relies on the use of cylindrically-focused transducers. In the same way, pure compression waves were considered, neglecting all other second-order phenomena (refraction on the defect, mode conversion, amplitude modification, and dispersion). No shear waves were studied in the sample regardless of the incidence conditions. Therefore, only the time-of-flight (TOF) of the acoustic wave was taken into consideration and processed.

In the ray-tracing approximation, using the wavefront construction method, the current wavefront for a travel time T is calculated from the previous wavefront with a displacement magnitude equal to V*Δt as a continuous envelope (Figure 33). The velocity value V was obtained using the Christoffel equation (Equation 13), as presented in Chapter 2. Time step Δt was a fixed value. The algorithm creates new points when the distance between two subsequent points is larger than a predefined maximum distance, set to be twice the distance between two points in the initial wavefront. When the wavefront presents a concavity, a crossing of the normal vectors can occur in the next iteration and a cusp appears, so the algorithm must delete crossing points. To identify crossing points, the algorithm goes through every point in the curve looking if the angle formed with the emitter is a monotonic function. If the monotonicity is no more verified at the current point, this point is deleted until the function becomes monotonous. Defective areas slow down the wave velocity. For this method, in the presence of a defect, the velocity value was set to be a constant value (isotropic behavior), computed as a fraction of the velocity in the healthy area. Depending on the decay degree, different percentages of reduction can be used. For example, in the case of a hole, velocity was set to be the sound speed in air, 343 m/s.

The algorithm for the wavefront construction method, developed in MATLAB (v9.5 2018, The MathWorks, Inc., Natick, Massachusetts, United States), follows the next steps:

I. Load elastic constants for the material: ER, ET, GRT, νRT, ρ.

II.

Compute rigidity matrix C.

III.

Define simulation parameters as the number of points for wavefront discretization, sensor positions, the sensor to act as an emitter, defect radius, defect position.

IV.

Compute the time step Δt.

V. While (not the final wavefront: iterate over the wavefronts) a. Iterate over the points in the wavefront:

i. Find the direction of the normal vector.

ii. Find the direction of the R axis.

iii. Compute the θ angle.

iv. Compute the phase velocity using the Christoffel equation (Equation 13).

v

Find the rays and Time of flight associated

To determine the time step Δt for every simulation, the smallest distance between sensors and therefore the smallest TOF to be computed were considered. For instance, the distance d1-2 between sensors 1 and 2 corresponded to the smallest TOF. The angle formed between these two sensors is near to the tangential direction (for 16 sensors, this angle is near to 80°), then the smallest TOF can be approximated by using the velocity for the tangential direction VT. For an error below 10%, the time step was defined as:

Δ𝑡 = 0.1 𝑑 1-2 𝑉 𝑇 (25)

Methodology for the numerical testing

To evaluate how anisotropy affects the wavefront shape compared to the isotropic case, the same wood elastic constants used for the sensitivity analysis of the Christoffel equation (Section 2.4, Table 2) were selected. For comparison, elastic constants for an isotropic material (PVC) were also defined [START_REF] Cardarelli | Materials handbook: a concise desktop reference[END_REF]. Table 10 sums up the wood species selected for this study, ordered by density, showing their elastic constants and their density, including the isotropic material. For softwood and hardwood, one species with a low, medium, and high anisotropy ratio were chosen, for a total of 6 species. Trunk geometry consisted on a cylinder, with a diameter of 30 cm. The number of sensors was fixed to 16.

Time step Δt for the selected species ranged from 4 µs to 5 µs. To assess the effect in the TOF estimation when a defect is present in the trunk, and the influence of the position and size of the defect, defects were simulated by defining a circular region with three diameters (2 cm, 4 cm, and 6 cm) and in different positions in the trunk, in horizontal, vertical and diagonal offsets from the center. To try different decay stages, the velocity was reduced by a percentage going from 30 to 70% of radial velocity.

Results

Wavefronts shape: orthotropic vs. isotropic

Figure 34 shows the wavefronts obtained for the isotropic material compared to the orthotropic (spruce in this case). As expected, the wavefronts in the isotropic case were perfectly circular, corresponding to concentric circles around the emitter. For the orthotropic case, a deformation was clearly seen, presenting a triangle shape for the first wavefronts, showing a faster velocity value in the radial direction; after passing the pith, waves presented a spherical shape. The propagation angles with slower velocities corresponded to values closer to the tangential direction. When the wavefront passed the pith, the direction of propagation coincided with the radial direction, which means higher velocity values in this area.

Rays and TOF estimation for a healthy case

For wood, the ray paths between the emitter and the receivers were not straight.

Figure 35 presents the corresponding ray paths for the simulation of isotropic and orthotropic materials using 16 sensors. First consideration with curved rays is that the distances traveled differed from the straight-line distance between emitter and receivers. A faster path may result in a longer way. To compare, the distance from each receiver to the emitter is presented for isotropic and orthotropic cases in Figure 36. The distances for the orthotropic case were longer than for the isotropic case, particularly for sensors 6 and 7 and their symmetric pair 11 and 12. Figure 37 shows the ratio between the curved distances traveled by the ray and the straight-line distances from the emitter to the receivers 2 to 9 (considering symmetry, sensors 10 to 16 are equivalents). For the six species, the sensor 6 and 7 (thus 11 and 12)

presented in average variations of 5% and 6%. The maximum variation was obtained for sensor 6 in the Spruce species, equal to 8%. The TOF values estimated for every receiver are presented in Figure 38, for all species, and the corresponding mean values for hardwoods and softwoods. For the isotropic case, the TOF value was constantly increasing until arriving at their maximal in the radial direction (sensor 9). For the orthotropic cases, TOF increased up to the fifth sensor but time values for sensors 6 to 9 were almost the same, creating a plane region in the TOF curve.

Average TOF values for softwood were larger than values for hardwood, considering that maximum velocity values for hardwoods were higher than the values for softwoods. For the isotropic case, the velocity values remained constant for all the sensors, since waves propagate with the same velocity in all directions. For the orthotropic case, the velocity values were higher in the radial direction, with the maximum value found at the sensor located on the opposite side of the trunk (sensor 9). Maximum values were higher for hardwoods. The minimum value changed its position for hardwoods and softwoods: for hardwoods, this value was located at sensors orientated closer to the tangential direction (sensors 2 and 16); for softwoods, this value was in sensors 3, 4 or 5 (depending on the species) and their symmetric pairs 15, 14 or 13, which presented an angle between 45° and 68° with respect to the emitter. Considering the velocity values obtained with the curved rays as the correct value, Figure 40 First, the defect size was changed. Circular defects were chosen, with diameters 2, 4 and 6 cm. The defect was in the center of the trunk, and the velocity reduction in the decay region was fixed to 50% of the velocity value in the radial direction. Figure 41 presents the wavefronts and ray paths for the different defect diameters in the anisotropic and isotropic case. As the defect size increased, the ray paths needed more wavefronts to arrive at the receivers, meaning higher TOF values. In the orthotropic case, the sensors affected were those labeled 7 to 11; for the isotropic case, the changes were observed mainly in the sensor 9. Using these TOF measurements, Figure 43 presents the velocity estimation for both materials, considering the distance from curved rays (as for the following results of this paragraph). The change in velocity values was larger for the orthotropic case. The largest variation was obtained for sensor 9 in the radial direction, with a difference of velocity of 227 m/s in the orthotropic case and 35 m/s for the isotropic case.

To compare the reduction in velocity as the defect size increases, Figure 43 also presents the ratio of velocity reduction for sensors 6 to 9, computed as the velocity value V of these sensors, for a given defect size, divided by the velocity value in the healthy case (Vref). For orthotropic, the velocity decreased almost constantly for those sensors, with a more pronounced slope in the case of sensor 9, that decreased a 16% with respect to the reference velocity value. For the isotropic case, smaller variations were obtained, with a maximum of 3%, mostly affecting the sensor 9.

The second factor studied was the defect position. In this case, the defect location was changed from the center of the trunk in three ways: first in the horizontal axis, second in the vertical axis and finally in diagonal positions, relative to the position of the emitter. For the horizontal shift, the offset from the center was defined as -8 cm and 8 cm, halfway between the trunk center and the bark. Figure 44 presents the wavefronts and ray paths for the different defect positions in the anisotropic and isotropic case. Figure 45 presents the TOF variation for each defect position, for both cases. The defects located in the center of the trunk for the anisotropic material presented a higher variation of TOF measurements compared to the isotropic one. For anisotropic, the TOF variation in the centered case was close to 18%, a larger value with respect to the less than 1% variation in the offset cases. In Figure 46 is presented the velocity variation for both materials, with a corresponding maximum value (in anisotropic) for the centered case of 15%, and for the offset cases a value lower than 1%. Therefore, when the defect position shifted towards the bark, the velocity variations were weaker with respect to the centered case in wood. The defect position was shifted in the vertical direction. The offset from the center was defined again as -8 cm and 8 cm. Figure 47 presents the wavefronts and ray paths for the different defect positions in the anisotropic and isotropic case. were affected, as for example in the case of +8 cm where the only sensor affected was the sensor 9. This effect was also presented in the isotropic case but weaker. In Figure 49 is presented the velocity variation for both materials, with a maximum for the anisotropic case of 15% compared to 5% in the isotropic case. A comparison with offsets in diagonal directions was performed, for the orthotropic case, with respect to the centered case. The defect center was shifted 8 cm in angles from the center of 45°, 135°, 225°, and 315°. Figure 50 presents the wavefronts and ray paths for the different defect positions. The defects with a negative offset in the vertical direction (subfigures c and d) presented less affectation to the ray paths curvature compared to the two cases with a positive offset. To compare the effect, Figure 51 presents the TOF and corresponding velocity variations with respect to the healthy case. The centered defect presented the larger variation, with a TOF change of 18%, followed by the cases with a positive vertical offset, with 12% variation, and finally small variations in the negative vertical offset with less than 1% variation. To evaluate the variations in velocity value depending on the decay affectation on wave velocity, a reduction set to 30, 40, 50, 60 and 70% were tested for wood and isotropic material. Over 70% of reduction, the velocity for Spruce decreased above the sound speed in the air; under 30%, the velocity in the decay region was over the velocity in the tangential direction (1117 m/s). Figure 52 presents the wavefronts and ray paths for the two extreme cases, this is 30% and 70%, in the anisotropic and isotropic case. In the case of wood, the rays passed through the defect even when the velocity was reduced to 50%. In the isotropic case, as the defect appears the rays avoided the region. Figure 53 presents the TOF variation for each decay percentage in both cases. The variations on TOF values for anisotropic were larger, from 7% to 22% for sensor 9 compared to the variation from 2.5% to 3% in isotropic case. Also, in the anisotropic case, multiple sensors were affected, compared to the isotropic cases where the only sensor highly affected was the sensor 9 in the radial direction. Figure 54 presents the velocity variation for all decay cases and presents the relative velocity values with respect to the velocity in the healthy case for sensors 6 to 9. As the decay affectation increases, the velocity values decreased, mostly in sensor 9 in the radial direction.

Finite elements method model comparison

The Finite Elements Method (FEM) has been widely used to study the elastic wave propagation for isotropic and anisotropic media [START_REF] Smith | The Application of Finite Element Analysis to Body Wave Propagation Problems[END_REF][START_REF] Marfurt | Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations[END_REF][START_REF] Serón | Finite-element method for elastic wave propagation[END_REF][START_REF] Lord | Developments in ultrasonic modeling with finite element analysis[END_REF]). Furthermore, specific applications in wood materials have also been presented [START_REF] Mackerle | Finite element analyses in wood research: a bibliography[END_REF], for example for the mechanical characterization of a wooden bar [START_REF] Veres | Wave propagation in a wooden bar[END_REF] or more recently for the ultrasonic evaluation of wooden poles [START_REF] Tallavo | Experimental verification of an orthotropic finite element model for numerical simulations of ultrasonic testing of wood poles[END_REF]). Nowadays, it is possible to find FEM modeling and simulation tools that allow a fast implementation and a variety of solvers adapted to different requirements. Here, we were interested in using a FEM model to study the propagation of elastic waves in a wood 2D cylindrical section to compare and validate the results obtained with the proposed raytracing approach.

Model definition

The FEM model was created using the Abaqus FEA software (v6.14, Dassault Systemes). The geometry was defined in a 2D planar environment as a circle with a diameter of 30 cm. For the defective case, a circular region of diameter 6 cm found at the center of the disk was defined. The material was set to present local orthotropy, with the rotation center defined to be the center of the trunk and the radial direction always pointing from the bark to the trunk center. Unlike the raytracing approach, the FEM simulation considered the presence of shear waves and effects of diffraction, diffusion and mode conversion. Attenuation was not considered (purely elastic).

The elastic constants were taken from Table 10 for the six-wood species. For the defective area, it was considered an isotropic behavior. To obtain the mechanical parameters for the defective region Eiso, νiso and ρiso, two of them were fixed: νiso=0.4 and ρiso=ρ; Eiso was obtained by defining a velocity in the defective region equivalent to 50% percent of the maximum velocity in the healthy region.

Time step was set to 0.1 µs (as for the raytracing case to allow a comparison) and the total simulation time was 300 µs (to observe the first arrival time for all cases). For the excitation, a dynamic load was defined. It corresponded to a square impulse with a short duration of 0.4 µs (for convergence of the simulation). An explicit solver was used, with free boundary conditions. As an impulsive signal was used for excitation, TOF values for every contour node were obtained by using the AIC method.

To define the elements size for the mesh, several experiments were performed, changing the element size, and the corresponding total number of elements. To evaluate the precision of the results obtained for a given mesh, an error measure was obtained, based on the TOF measurements as:

𝑒̅ (𝑛) = ∑ | 𝑇𝑂𝐹 𝑖,𝑛 -𝑇𝑂𝐹 𝑖,𝑛𝑡𝑜𝑡𝑎𝑙 𝑇𝑂𝐹 𝑖,𝑛𝑡𝑜𝑡𝑎𝑙 | 𝑘 𝑖=1 𝑛 𝑡𝑜𝑡𝑎𝑙 , ( 26
)
Where n is the test number and ntotal is the total number of experiments, with experiment number 1 having the least quantity of elements and experiment number 20

having the maximum number of elements; k is the total number of nodes at the trunk contour for the corresponding element size.

Element size was chosen to range from 1.2 mm to 15 mm, with a corresponding total number of elements ranging from 747 to 100.695 in the healthy case and 769 to 101.991 in the defective case. Figure 55 shows the mean error for both cases, healthy and defective, for the corresponding total number of elements. The selected element size was 1.5 mm (around 65.000 elements) that reduces the mean error to 0.25%.

Results

Wavefronts comparison

Wavefronts obtained with both methods were compared. Figure 56 presents the superposition of the simulations performed for the Spruces Sitka species for healthy and defect cases. The border of the wave propagating on the FEM simulation agrees with the corresponding wavefront obtained at the same time step with the raytracing approach.

TOF Comparison

For the specified element size, the number of nodes in the contour was 628 (one sensor every 0.57°). To compare, raytracing simulations were performed using the same number of sensors, and TOF estimations were obtained for every receiver. For healthy case, the comparisons are presented in Figure 57 (a) and (b); for defect case in (c) and (d). The horizontal axis shows the angle in degrees between the receiver and the emitter. The TOF curves obtained with both methods were similar. For the healthy case experiments, a flat region was obtained for receivers whose angle was lower than 40°. For the defective case experiments, an increasing TOF value was obtained starting around 40° and reaching the maximum at 0°. To quantify the difference between the two models, the mean absolute difference was obtained for healthy and defective experiments, as presented in Figure 58. For both cases, the mean absolute difference never exceeded 3 µs, that is a lower value than the time step used for the raytracing simulations. The mean relative difference was lower than 2% when the angle was lower than 85%. Both models resulted in similar TOF estimation under the specified mechanical parameters; however, raytracing may be considered as a less complex technique, resulting in lower processing times. The highest differences were obtained at the nodes found at angles over 85°. In this case, the nodes were really close to the emitter position and the TOF values were comparable to the time step for the simulations. This considering that the number of sensors used for the simulations was high and therefore distances from the emitter were small (in real applications, 32 sensors corresponds to a high value; for simulations, 628 sensors were used).

Experimental validation

Experimental setting

Healthy trunks from Oak (Quercus rubra) and Ash (Fagus sylvatica) were used to obtain 8 wood disks (4 for each species). Oak trunk diameter was 20 cm and Ash trunk diameter was 30 cm. Tree ages were 27 years and 42 years for oak and ash respectively.

Disks thickness was 3 cm for all cases. Disks were debarked. Through the experiments, the trunk was stored with controlled temperature and sealed to reduce water loss. The room where the disks were stored had a temperature regulation fixed to 4°C to reduce water loss and to keep moisture content over the fiber saturation point. Disks were weighted before and after the tests, showing that moisture content reduction was lower than 2%. The moisture content of the disks during the ultrasonic tests ranked between 70% -73% for oak and between 35% -38% for ash. After the ultrasonic tests, the disks were stabilized at 20°C and 65% of relative humidity to determine their density at a moisture content of 12% (577 kg/m3 for the Oak and 690 kg/m3 for the Ash).

Defects were simulated by drilling a circular hole, as shown in Figure 59. Three defect diameters were tested: 2.9 cm, 5.1 cm, and 7.6 cm. The defects were tested in two positions: centered (for oak disks) and off-centered (for ash disks). Ultrasonic tests were performed in 16 positions around the wood disks. The transmitters were placed at an equal distance between each other all around the trunk. To compare with an isotropic material, a PVC disk was used (Figure 60). PVC disk diameter was 30 cm with a thickness of 3 cm. Ash disks tested: (e) Healthy, (f-g) off-centered defects with diameters 2.9, 5.1 and 7.6 cm [START_REF] Espinosa | Effect of wood anisotropy in ultrasonic wave propagation: A ray-tracing approach[END_REF]. The ultrasonic chain of measurement was set up as presented in Chapter 3. TOF measurements were repeated 16 times and the average values were used; relative standard deviation values were around 10%.

For the raytracing simulations, elastic constants and density values for both wood species (at 12% moisture content) were obtained from published data [START_REF] Guitard | Mécanique du matériau bois et composites. Cepadues Editions Habermehl A, Ridder HW[END_REF]) and for PVC from the Materials Handbook [START_REF] Cardarelli | Materials handbook: a concise desktop reference[END_REF], and they are presented in Table 11. For wood, elastic constants ER, ET, and GRT, and density ρ, were corrected to consider the moisture content when green, using the correction equations proposed by Guitard [START_REF] Guitard | Mécanique du matériau bois et composites. Cepadues Editions Habermehl A, Ridder HW[END_REF], the Wood Handbook (Ross 2010), and Sobue [START_REF] Sobue | Simulation study on stress wave velocity in wood above fiber saturation point[END_REF], for moisture fluctuations beyond the fiber saturation point (30%). For the PVC case, ΔT was set to 4.3 µs. For Oak, ΔT was set to 3.3 µs and for Ash was set to 4.9 µs.

Results

Experiments with healthy wood

Figure 61 presents the wavefronts and the corresponding ray paths obtained for the isotropic material compared to the orthotropic. Ash results were used for this comparison, considering that its diameter was like the PVC disk diameter. For the isotropic case, with straight rays, TOF was monotonic increasing as the sensor approached from the tangential (α=90°) to the radial direction (α=0°). For wood, the sensors located after the pith, sensors 6 to 12, resulted in similar TOF estimations. The shape of TOF curves computed using the raytracing algorithm and those obtained from the experiments with the PVC and Ash disks are in good agreement, showing that the proposed model simplification is usable.

For the PVC case, simulated values ranged between 43 µs and 222 µs, while the experimental values ranged between 48 µs and 213 µs. For the Ash case, the simulated values ranged between 48 µs and 190 µs, with a mean TOF value for sensors 6 to 12 of 180 µs; the experimental values ranged between 24 µs and 178 µs, with a mean TOF value for sensors 6 to 12 of 171 µs.

Experiments with centered defects

Figure 62 shows the raytracing simulations for the Oak disks when centered circular defects were created. Also, this figure shows the TOF measurement comparison between the simulated and the experimental data. Ray paths avoided the defective regions, as they slowed down the waves, increasing the time required for the wave to arrive from the transmitter to the receivers located in the direction of the defect. Considering that the defect was in the center, rays traveling to the sensors located after the pith (6 to 12) were affected, presenting a TOF measurement that increased proportionally to the defect size. Figure 63 illustrates the TOF differences when comparing the defective cases to a healthy one. The larger variation was obtained in the radial direction. For this direction, considering simulated data, TOF increased 13 µs for the first defect size (11% increment), 26 µs for the second size (22% increment) and for the largest defect 46 µs (34% increment).

When the receiver angle was lower than 45° (receivers before the pith), the TOF variation was reduced to values under 5 µs. For the radial direction considering the experimental data, TOF increased 18 µs (12% increment) for the smallest defect, 50 µs for the medium size defect (39% increment) and 75 µs for the largest defect (50% increment). Differences for simulated and experimental data followed a similar pattern, with larger variations mostly in the radial direction for the experimental data.

Experiments with off-centered defects

Figure 64 shows the raytracing simulations for the Ash disks when off-centered circular defects were created. The defect position was shifted horizontally from the trunk center. Also, this figure shows the TOF measurement comparison between the simulated and the experimental data.

The TOF curves for the simulated data and the experiments were similar. For the simulated values, minimum TOF was 49 µs and maximum values ranged between 190 µs and 195 µs as the defect size increased; for the experimental values the minimum TOF ranged between 24 µs and 34 µs, and the maximum values ranged between 179 µs and 184 µs as the defect size increased. TOF differences for the off-centered case with respect to the healthy reference were lower than those obtained for the centered case. Figure 65 shows the TOF differences when comparing the defective cases to the healthy one. Considering simulated data, for the receiver number 6 (between angles -40° and -30°), TOF increased 0.62 µs for the first defect size (0.3% increment), 3.4 µs for the second size (1.9% increment) and for the largest defect 7.8 µs (4.2% increment). For the other receivers, the TOF variation was almost zero.

Considering the experimental data, the larger variation was presented in sensor 7 (between angles -30° and -20°), where TOF increased 4.4 µs (2.6% increment) for the smallest defect, 8.9 µs for the medium size defect (5.4% increment) and 17 µs for the largest defect (10% increment). As for the case of centered defects, TOF variations were larger for the experimental data. In the second set of configurations, the defect was shifted from the trunk center in the vertical axis to the bottom (Figure 66a) and to the top (Figure 66d), with respect to the emitter. To obtain these defect positions with the same Ash disk, the sensor acting as the transmitter was changed to positions 5 and 13. The simulation for both defect positions and the TOF estimations for the experimental and simulated data can be seen in Figure 66.

Maximal variations were found in the radial direction (θ=0°). If the defect is located at the bottom of the disk, for the simulated and experimental data, variations were obtained for sensors 7 to 11, with a maximum for the simulated data of 31 µs (15% increment) for the sensor 9, compared to 30 µs (14.9% increment) for the experimental data. When the defect was located at the top of the disk and furthest from the transmitter, the simulation showed only sensor 9 being affected. The largest variation was in the radial direction and corresponded to an increment of 13 µs (6.8% increment). Experimental data variations were present at sensors 8, 9 and 10, with a maximum of 35 µs (18%) in the radial direction. Comparing the TOF curves when the defect was shifted in the two directions, defects shifted vertically were prone to higher variations compared to those shifted horizontally. Moreover, both off-centered cases resulted in lower variations compared to the centered defect cases. It should be noted that off-centered defects will result in less contrasting regions in the tomographic images compared to centered cases, which increases the difficulty in identifying defect location.

Discussion

In the case of healthy wood, for the sensors located on the disk top half (sensors 6 to 12), ray paths tended to pass first for the center of the trunk to follow a faster path, therefore those trajectories were longer with respect to the straight-line distance. In the case of real tomography imaging, there is no possibility to know the real paths, so the common procedure consists of measuring the straight-line distances. Aiming to reduce the affectation in the reconstruction process, a correction for the distance measurement could be performed, using a factor obtained from simulations with a healthy trunk from the species under study, like the procedure presented by Maurer et al. [START_REF] Maurer | A simple anisotropy correction procedure for acoustic wood tomography[END_REF].

TOF estimation for the wood healthy case resulted in similar values for sensors ranging from 6 to 12. For those receivers, trajectories tended to present two ray parts: one part was orientated from the sensor to the pith, and the second part was orientated from the pith to the transmitter. Both parts presented similar velocity values, meaning similar TOF values.

Larger defects could be associated with larger changes in TOF estimation and therefore in velocity variations. This means that larger defects will present more contrast with respect to the healthy areas, allowing an easier identification. This result has been already observed in other studies [START_REF] Schubert | Acousto-ultrasound assessment of inner wood-decay in standing trees, possibilities and limitations[END_REF][START_REF] Zhang | Experimental investigation of stress wave propagation in standing trees[END_REF]. With respect to the defect position, those defects whose center were located on the path between the emitter and one receiver presented a larger TOF variation in that receiver. If the ray path is supposed to pass by the defect center, the ray curvature is going to be maximal. Moreover, defects located in the trunk center tended to cause a larger affectation to the TOF value estimation in orthotropic case. Most of the ray paths above the pith crossed near to the trunk center, so defects in this position are likely to modify the velocity in multiple sensors. Also, the ray paths crossing by the center are privileged directions given that the velocity is maximal near to the radial direction. In the reconstruction process, this will affect the image quality depending on the position of the defect (unknown for real cases).

A combination of centered position and a bigger size will correspond to a higher probability of decay detection using a tomographic image. However, it is necessary to consider the anisotropic behavior of wood when performing the ultrasonic image reconstruction. When classical reconstruction techniques are used, such as the filtered backprojection algorithm, the first assumption is that trajectories are straight lines (isotropic behavior). Then, if the material presents isotropy, and no defect is present, a flat image is obtained corresponding to the wave velocity of propagation. However, if the material presents the wood orthotropy (curved rays), using those methods in a healthy section will result in a gradient of velocities, with the higher velocities in the center and with lower velocity areas near to the border, making difficult the image interpretation (Espinosa et al. 2017b).

Even when the shape of the TOF curves agreed between the simulated and the experimental data, there was an offset between those values. This difference, more appreciable in the Oak case, could be associated to using mechanical parameters for the raytracing simulation from published data, considering that the mechanical parameters for wood species are subject to a large variability and to the effect of factors such as the moisture content.

Synthesis

Ray-tracing approximation using the wavefront construction method allowed to simulate the wave propagation in orthotropic media, and therefore allowed to compute time-of-flight estimations. Anisotropy in the radial-tangential plane in wood resulted in deformed wavefronts with respect to the isotropic case. The paths from each receiver to the transmitter in the wood presented a curvature, therefore the distances traveled differed from the straight-line distance obtained for the isotropic case. Velocity values were higher in the radial direction, with the maximum in the sensor located on the opposite side of the trunk.

Defects located in the center of the trunk presented larger TOF variations compare to defects located in off-centered positions. Off-centered defects located at the top and at the bottom of the disks with respect to the emitter position presented a larger variation than the off-centered defect shifted horizontally from the disk center. Off-centered defects will be more difficult to determine and characterize by tomographic inversion. A comparison with another numerical simulation method, the FEM approach, resulted in similar TOF and wavefronts estimation. An experimental validation of the raytracing method was performed, using wood disks from two species and simulating the presence of defects by drilling holes. The shape of TOF curves obtained with the raytracing algorithm and those obtained from the experimental validation were in good agreement.

ALGEBRAIC SOLUTION OF THE INVERSE

PROBLEM FOR AN ORTHOTROPIC MATERIAL

Introduction

To perform USCT, it is necessary to know the path followed by the wave from the transmitter to the receiver. In Chapters 2 and 4 it has been shown that for anisotropic materials such as wood, the ultrasonic waves velocity depends on the angle of propagation (Espinosa et al. 2017a[START_REF] Espinosa | Effect of wood anisotropy in ultrasonic wave propagation: A ray-tracing approach[END_REF]. Waves propagating in the radial direction (from the bark to the pith) have a higher velocity than those propagating in the tangential direction (perpendicular to the radial direction). This approach to the forward problem showed that the resulting wave paths are curved [START_REF] Espinosa | Effect of wood anisotropy in ultrasonic wave propagation: A ray-tracing approach[END_REF].

Approaches for tree USCT as inverse problem, have used reconstruction techniques that consider straight-line paths (rays), like the filtered backprojection (FBP) method [START_REF] Tomikawa | Nondestructive Inspection of a Wooden Pole Using Ultrasonic Computed Tomography[END_REF]) and algebraic methods such as SIRT (Simultaneous Iterative Reconstruction Technique) [START_REF] Nicolotti | Application and comparison of three tomographic techniques for detection of decay in trees[END_REF][START_REF] Socco | Feasibility of ultrasonic tomography for nondestructive testing of decay on living trees[END_REF][START_REF] Martinis | Tomographie ultrasonore pour les arbres sur pied[END_REF]).

Using a straight-ray approximation results in a biased image. These techniques deliver a velocity map but considering that this parameter is a function of the angle of propagation, the physical sense of the image is unclear.

This chapter is related to the development of a 2D inversion method adapted to wood anisotropy. Aiming to obtain for every pixel a velocity/slowness function depending on the propagation angle, curved rays must be considered. Using this slowness function, it is possible to estimate for every pixel the associated mechanical parameters, such as Young's modulus in the radial and tangential direction, shear moduli and Poisson's ratio.

The curved paths are not known a priori, thus the proposed method iteratively estimates them from an initial straight-ray hypothesis, using a raytracing approach. First, the method was evaluated numerically. Four numerical configurations were tested representing real cases usually found in standing tree monitoring [START_REF] Wang | Acoustic tomography for decay detection in red oak trees[END_REF]): a typical-trunk in a healthy state, with an air-filled centered hole, with an off-centered hole, and with an off-centered defect characterized by a gradient of mechanical properties. The proposed inverse method, tested in these configurations, were compared to a reconstruction method based on straight rays. Then, the inversion method was tested using experimental data from cross-sections of two wood species. Healthy, centered and off-centered cases were tested, simulating the defects by drilling holes and comparing again with straight-line reconstructions.

Proposed inversion method

Algebraic formulation

Let us consider a space divided into N pixels. For every curved path 𝑚 traversing this space, the TOF 𝑡 𝑚 can be obtained by adding individual slowness 𝛼 𝑘 for every pixel k along the path, multiplied by the length of the ray segment 𝑙 𝑚𝑘 . This can be written as:

𝑡 𝑚 = ∑ 𝑙 𝑚𝑘 𝛼 𝑘 𝑘 𝑎𝑙𝑜𝑛𝑔 𝑚 (27)
The total length of a ray 𝐿 𝑚 is equal to the sum of individual ray segments 𝑙 𝑚𝑘 .

Assuming this segment to be uniform (equal pixel length), we have 𝐿 𝑚 = 𝐷 𝑚 𝑙 𝑚 , with 𝐷 𝑚 as the total number of ray segments in the path m. Dividing both sides of Equation 27by 𝐿 𝑚 we have:

𝐴 𝑚 = 𝑡 𝑚 𝐿 𝑚 = 1 𝐷 𝑚 ∑ 𝛼 𝑘 𝑘 𝑎𝑙𝑜𝑛𝑔 𝑚 (28)
With 𝐴 𝑚 as the total slowness considering the ray. Then, an equation can be formulated for every pair of emitter-receiver, creating a system of linear equations. The inverse problem corresponds to the solution of this equation system, leading to the reconstruction of 𝛼 𝑘 for all pixels from its projections. The maximum number of pixels 𝑁, without considering the possibility to perform an interpolation, is linked to the number of sensors M as:

𝑁 = 𝑀 * (𝑀 -1) 2 (29)
A matrix formulation can be proposed considering a positioning of pixels into a grid (Figure 67). In that case, for every curved path m, all the N pixels in the image are arranged in a vector format and a matrix 𝑮 filled with 0's and 1's will define which pixels corresponds to every path and their corresponding total slowness 𝐴 𝑚 , as follows:

[ 1 0 ⋯ 1 0 1 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 1 0 ⋯ 1 ] [ 𝛼 1 𝛼 2 ⋮ 𝛼 𝑁 ] = [ 𝐷 1 𝐴 1 𝐷 2 𝐴 2 ⋮ 𝐷 𝑚 𝐴 𝑚 ], (30)
𝑮𝜶 = 𝑨. ( 31
)
Figure 67: Pixel grid distribution for the algebraic reconstruction [START_REF] Maurer | A simple anisotropy correction procedure for acoustic wood tomography[END_REF]).

However, Equation 30 considers that slowness for every pixel is independent of the propagation angle (algebraic problem [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF]). To consider the effect of this angle, we must modify the 𝑮 and 𝜶 matrices to include the slowness as a function of the propagation angle 𝜃 (using the Christoffel equation). Let us assume that for every pixel we know the propagation angle 𝜃, i.e., we know the ray paths. To linearize the Christoffel equation (Equation 13), a 5 th degree polynomial approximation is proposed, such as:

𝛼 𝑘 = 𝛽
This degree is sufficient to approximate the slowness curve, as observed in Figure 68. This figure shows an estimation of the coefficient of determination R 2 for several polynomial degrees, using as example the slowness curve computed with the elastic parameters of Oak (Table 11). 11).

In that case, the solution of the problem corresponds to finding the polynomial coefficients 𝛽 for every pixel. Replacing Equation 33 in the Equation 30, we obtain:

[ 𝜃 1,1 5 ⋯ 𝜃 1,1 0 0 ⋯ 0 ⋯ 𝜃 𝑁,1 5 ⋯ 𝜃 𝑁,1 0 0 ⋯ 0 𝜃 2,2 5 ⋯ 𝜃 2,2 0 ⋯ 0 ⋯ 0 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ 𝜃 1,𝑀 5 ⋯ 𝜃 1,𝑀 0 0 ⋯ 0 ⋯ 𝜃 𝑁,𝑀 5 ⋯ 𝜃 𝑁,𝑀 0 ] [ 𝛽 5,1 ⋮ 𝛽 0,1 𝛽 5,2 ⋮ 𝛽 0,2 ⋮ 𝛽 5,𝑁 ⋮ 𝛽 0,𝑁 ] = [ 𝐷 1 𝐴 1 𝐷 2 𝐴 2 ⋮ 𝐷 𝑚 𝐴 𝑚 ], (34)
𝚯𝜷 = 𝑨. (35)
So, if we can solve Equation 35to obtain the matrix 𝜷, a reconstruction of the slowness function for every pixel is obtained.

SIRT method [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF] was used to solve Equation 35and to obtain the matrix of Christoffel polynomial coefficients 𝜷. This algebraic iterative algorithm solves the linear equations system from an initial guess. SIRT method has shown better convergence compared to others reconstruction methods, leading to images of higher quality (Arciniegas et al. 2014a).

From the polynomial slowness function for every pixel, a nonlinear regression can be performed to obtain the elastic constants. The Levenberg-Marquardt nonlinear least squares algorithm was used [START_REF] Seber | Nonlinear Regression[END_REF]. This method is used for fitting a nonlinear parameterized function (the Christoffel equation) to a set of measured data points (the estimated polynomial slowness function), through minimization of the sum of the squares of the errors. This minimization process is performed iteratively by updating the parameter values, in this case, the elastic constants. To reduce the error on the estimation of the elastic constants, the Poisson's ratio value was not neglected (Espinosa et al. 2017a) but set to a value of 0.56 [START_REF] Guitard | Mécanique du matériau bois et composites. Cepadues Editions Habermehl A, Ridder HW[END_REF], leaving three variables in the model to be estimated (𝐸 𝑅 , 𝐸 𝑇 and 𝐺 𝑅𝑇 ).

For a given set of TOF values and their corresponding curved trajectories, the elastic constants (from the velocity/slowness function for every point in the space), can be obtained using the proposed inversion method.

Reconstruction schema

Considering a set of TOF measurements, an initial guess for the trajectories can be straight-line paths. This initial guess allows us to perform a first inversion process, obtaining the Christoffel polynomial coefficients (slowness) for every pixel. Using these coefficients, we can perform the direct problem to obtain an estimation of the TOF and the trajectories. Then, we can correct the rays to reduce the TOF difference by applying feedback, as presented in Figure 69. To deal with trajectories, they were approximated by a 3 rd -degree polynomial curve (4 unknown parameters 𝛾 0 , 𝛾 1 , 𝛾 2 and 𝛾 3 ), as a function of the spatial coordinates of the discrete points that constitute the ray. Given that two points were known (emitter and receiver positions), two coefficients can be calculated directly (𝛾 0 and 𝛾 1 ). Thus, using a 3 rd -degree polynomial approximation means that only two variables were used to specify a single trajectory (𝛾 2 and 𝛾 3 ).

Stop criterion is related to the evaluation of the root mean squared error between the real TOF and the estimated one and the differences between the coefficients for ray trajectories (defined by the polynomial coefficients 𝛾 3 and 𝛾 2 ), as:

𝑒 = √ ∑ (𝑇𝑂𝐹-𝑇𝑂𝐹 ̂)2 𝑁 𝑁 * √ ∑ (𝛾 2 -𝛾 2 ̂)2 𝑀 𝑀 * √ ∑ (𝛾 3 -𝛾 3 ̂)2 𝑀 𝑀 . ( 36
)
The algorithm stops iterating when the error stops decreasing or when a maximum number of iterations is reached (set to 10).

The proposed reconstruction method was developed in Matlab (v9.5 2018, The MathWorks, Inc., Natick, Massachusetts, United States). The average processing time for a single iteration was 6 minutes using a desktop PC with Intel Core i5 processor (2x 1.7

GHz) and 8 GB of RAM. Stop criterion was reached with less than 5 iterations in this study.

Numerical validation

Methodology

Elastic constants and density values were obtained from published data [START_REF] Guitard | Mécanique du matériau bois et composites. Cepadues Editions Habermehl A, Ridder HW[END_REF]. Selected wood species was Oak, with: 𝐸 𝑅 =1180 MPa and 𝐸 𝑇 =614 MPa as Young's modulus in the radial and tangential directions respectively, 𝐺 𝑅𝑇 =319 MPa as the shear modulus, 𝜈 𝑅𝑇 =0.56 as the Poisson's ratio and 𝜌=570 kg/m 3 as the density. Trunk geometry consisted on a circular disk, with a diameter of 30 cm. The number of sensors was fixed to 32, considering that this number of sensors corresponds to the maximum for commercial devices, and it was previously shown that several reconstruction methods converge from 30 transducers (Arciniegas et al. 2014a).

Four configurations were tested. First, a healthy case where the elastic parameters were homogeneous in the whole disk. Second, a centered defect was created in the disk, corresponding to a hole with a diameter of 10 cm. Inside the defective region, the velocity value was set to be a constant value (sound speed in air, 343 m/s). Third, an off-centered defect, which diameter was 10 cm shifted horizontally 7.5 cm from the center. Finally, a case corresponding to a gradient of mechanical properties: 𝐸 𝑅 , 𝐸 𝑇 and 𝐺 𝑅𝑇 were reduced linearly from 100% to 50%, in an off-centered circular region (as in the third case), from the border of the region to its center. This last case could be associated to an early stage of wood degradation by fungi attack. For comparison, a reconstruction using the FBP method (straight-ray inversion) was performed.

Results

Tree in a healthy state

Convergence of the method was obtained after 5 iterations. and 1200 MPa. Concerning Young's modulus in the tangential direction, pixels located in the disk center presented higher values than the real value of 614 MPa, with a maximum error in the center value that went over 1000 MPa; however, peripheric pixels were well estimated, ranging from 600 MPa to 700 MPa. For the shear modulus, behavior was similar, with a higher value in the center, and peripheric values ranging between 250 MPa to 300

MPa. Thus, we can observe that the more suitable image for evaluation would be the map of 𝐸 𝑅 values.

Trunk with a centered hole

Convergence of the method was obtained after 3 iterations. Figure 72 presents the reference trajectories and the trajectories estimated with the proposed inversion method.

Again, the estimated trajectories corresponded to a correct approximation of the raytracing reference trajectories, where rays bypassed the defective region. Figure 73 presents the reconstructed images and the corresponding profiles. With the FBP method, the correct size of the defect was difficult to be estimated, due to the velocity gradient highlighted in Figure 71. The images of the three elastic parameters reproduced accurately the hole position and shape with the proposed method.

Trunk with an off-centered hole

Convergence of the method was obtained after 3 iterations. Figure 74 presents the reference trajectories and the trajectories estimated with the proposed method. In this case, trajectories were not all well estimated, as some of them passed through the defect. Figure 75 presents the reconstructed images and the corresponding profiles. As for the centered defect, the hole is difficult to distinguish due to the presence of a gradient of velocities using the FBP algorithm. With the proposed method, the presence of a hole was better defined than for the FBP method, but its actual size was underestimated compared to the central hole case.

Trunk with an off-centered defect

Convergence of the method was obtained after 3 iterations. Figure 76 presents the reference trajectories and the trajectories estimated with the proposed method. The gradient effect on the trajectories was more difficult to capture than the presence of a void region.

Figure 77 presents the reconstructed images and the corresponding profiles. This gradient is almost indistinguishable with the FBP method. Using the proposed method, the linear gradient appeared in the images, but the size was again underestimated.

Discussion

Approximations were made to regularize the problem, using polynomial functions associated to the Christoffel equation and to the trajectories. Regarding the trajectories, the 3 rd -degree polynomial approximation resulted in smooth ray paths that could restrain the possibility of finding complex-shape defective areas. Nevertheless, the proposed method handled with the rays' curvature due to the orthotropy property of wood, that improved the identification of the defects compared to the straight-ray reference method.

The number of sensors used for testing is directly connected to the reconstructed image resolution, as stated in Equation 29. Moreover, the number of sensors and the corresponding total number of pixels influence the processing time of the algorithm, as the size of the matrices to be handled increases fast and algorithms like SIRT are highly timeconsuming. For in situ testing, the duration of an acoustic evaluation of a standing tree is a crucial factor; thus, for a real case, increasing the number of sensors could not be possible.

Another option aiming to increase the image resolution consists of performing an interpolation of the sinogram (interpolation of the projection data before inversion), needing an accurate TOF estimation (Arciniegas et al. 2014b). The pixel size was 1.2 cm, meaning that a defect under that size could not be found. However, the minimum diameter to achieve an identification of a centered defect was 5 cm for a disk with a diameter of 30 cm, corresponding to a ratio defect diameter to disk diameter of 17% (Figure 78). Moreover, off-centered defects proved to be more difficult to identify, with a minimum diameter of 7 cm (ratio of 23%). Central pixels in the reconstruction of Young's modulus in the tangential direction (𝐸 𝑇 ) and the shear modulus 𝐺 𝑅𝑇 resulted in higher values than the real ones. Most of the trajectories passing through these pixels presented a 𝜃 angle near to the radial direction.

Then, there were not enough information of the whole set of angles to obtain the correct slowness function when using the nonlinear regression. To avoid this, constraints could be specified when reconstructing these parameters in the central pixels, such as considering the mean value in the peripheral values. Poisson's coefficient was set as its estimation using the nonlinear regression presented the largest variation when noise (due to measurement uncertainties) was present on the slowness function. To observe this effect, zero-mean normal-distributed noise with a standard deviation equal to k (%) times the function mean value was added to the slowness curve. Relative error for the four elastic parameters was computed as (𝑞 0 -𝑞)/𝑞, with 𝑞 0 as the estimated parameter value when noise was added, and 𝑞 as the real parameter value. When the Poisson's coefficient was set, the relative error was reduced to less than 1%, as presented in Figure 79.

Experimental validation

Materials and methods

Healthy trunks from pine (Pinus pinea) and oak (Quercus rubra) were used to obtain 6 wood disks (3 for each species). Average disk diameters were 20 cm and 24 cm for oak and pine respectively; thickness was 3 cm. Tree ages were 27 years and 55 years for oak and pine respectively. To reduce water loss during the experiments, the trunks were sealed and stored in a room with controlled temperature. The moisture content for oak was in average 72%, and for pine 92% (moisture reduction during the tests was less than 3%).

Defects were simulated by drilling a circular hole of diameter 7.6 cm (above the minimum detectable size obtained in the numerical testing). Centered and off-centered defect positions were tested. The off-centered defect was located halfway between the center and the bark. Ultrasonic chain of measurement was configured as presented in Chapter 3 (Figure 80). 16 sensors were placed around the disks to perform the TOF measurements. For comparison, images considering an isotropic assumption (FBP) were computed. As shown in the numerical validation, the image representing Young's modulus in the radial direction (ER image) was the more suitable for defect detection; then, this image was selected for testing. Profiles from both images were obtained, passing through the disk center. To increase the number of reconstructed pixels, a sinogram interpolation was performed, passing from 16 sensors to 32 sensors, by a linear interpolation of the TOF measurements and the sensors positions. Figure 81 shows an example of sinogram interpolation, for the case of the healthy oak disk TOF measurements. To perform an estimation of the decay detection, the images for the defective cases were segmented using a threshold. The threshold values were fixed using the mean value in the image, as 30%, 50%, 70% and 90% of this value. To evaluate the classification, the true positive rate (sensitivity) and the false positive rate (fall-out) were computed. The first corresponds to 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 is 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%.

Results

Results for Oak

The first case corresponded to a healthy disk of oak, shown in Figure 82a For the centered defect case, shown in Figure 83a., the trajectories were modified by the presence of the hollow region, as observed in Figure 83b. Compared to the FBP approach shown in Figure 83d, the proposed method allowed a more precise defect identification, as observed in the reconstructed image in Figure 83c Figure 84 presents the thresholding of both images to locate the defect. For the FBP images, those using thresholding of 30% and 50% did not identify any pixel as defective, then the sensitivity was 0% and the fall-out was irrelevant. For the thresholding values of 70% and 90%, sensitivity values were 13% and 69% respectively, with fall-out values of 8% and 19%. In the case of the ER images, an improvement on these metrics was obtained, with sensitivity values increasing to 95%, 98%, 98% and 99% for the threshold values of 30%, 50%, 70%, and 90% respectively, and associated fall-out values of 3%, 7%, 11% and 22%. The best balance, in this case, corresponded to the threshold value of 30% for the ER image. In this case, the thresholding of both images to locate the defect is presented in Figure 86. First, for the FBP images, sensitivity was found to be 0%, 13%, 37%, and 60%, as the threshold value increased. The fall-out values were respectively 2%, 7%, 14% and 24%. Improvements were obtained again with the ER image, with sensitivity values 45%, 51%, 60% and 78% and fall-out values of 1%, 4%, 6%, and 19% respectively. Best compromise between sensitivity and fall-out was obtained for the ER images using a threshold of 50% and 70%, but they were lower than those for the off-centered defect case.

Results for Pine

For pine samples, first, a healthy disk was tested, as shown in Figure 87a. The Figure 89 presents the thresholding results. In the case of FBP images, the first two threshold values did not find any defect. For the threshold values of 70% and 90%, sensitivity was 33% and 96% respectively, with fall-out values of 7% and 29%. For the ER image, improvements were obtained, with sensitivity values of 61%, 66%, 75% and 95%

as the threshold value increased, and respectively fall-out values of 5%, 8%, 11%, and 25%.

The best balance was obtained for the ER image with a threshold at 70%. Figure 91 presents the segmented images. For the FBP image, using the threshold value of 30%, no defective pixels were detected. For the other thresholds, the sensitivity values were 7%, 62%, and 91% and the respective fall-out values were 5%, 9% and 21%.

Using the ER images, improvements were obtained mainly for the sensitivity: values were 73%, 84%, 94%, and 95%, with fall-out values of 3%, 7%, 26% and 50%. Best compromise was obtained for the ER image using a threshold of 50%. Outliers for the ER images resulted in higher fall-out values than in the FBP image.

Discussion

For these tests, only one set of TOF measurements were obtained. To reduce uncertainty in the image construction, multiple measurements for every receiver should be necessary to improve the TOF estimation and therefore the image quality. Also, increasing the number of sensors could increase the total number of pixels and the defect detection precision, but this entails a larger implementation time in situ. Here, the sinogram interpolation was used to virtually add sensors and increase image resolution. Considering that the transit time corresponds to physical measurement, we can suppose that the TOF value between two adjacent sensors could be approximated by their mean value. In this case, the number of sensors was duplicated, so passing from 16 to 32 sensors. Similarly, the TOF value for this new sensor was interpolated from the values of the two surrounding sensors.

For the healthy cases, rays were curved with a more complicated pattern compared to the numerical simulation, due to factors that were not considered into the numerical model. In the case of the Pine samples, the profile obtained showed three different zones:

the lowest ER values around the pith, that could be associated to the presence of juvenile wood, the medium ER values in the left part of the profile, corresponding to normal wood, and finally higher ER values for the part with compression wood.

For juvenile wood, density is lower than for mature wood [START_REF] Ross | Wood handbook: wood as an engineering material[END_REF]. The lower wood density is related to the presence of shorter tracheids and thinner cell walls. The large microfibril angle in juvenile wood leads to a low stiffness and a low Young's modulus [START_REF] Barnett | Cellulose microfibril angle in the cell wall of wood fibres[END_REF]. Previous studies have shown a decreasing value of velocity of elastic waves for juvenile wood (Brancheriau et al. 2012a;[START_REF] Palma | Interference from Knots, Wave Propagation Direction, and Effect of Juvenile and Reaction Wood on Velocities in Ultrasound Tomography[END_REF].

Compression wood has a higher density [START_REF] Timell | Compression Wood in Gymnosperms[END_REF][START_REF] Kollmann | Principles of Wood Science and Technology : I Solid wood[END_REF], considering that the cell wall is much thicker in compression wood than in normal wood.

With the proposed method, the ER images were obtained using a fixed density value, for both normal wood and compression wood, thus the elastic modulus was overestimated.

Compression wood also presents a larger microfibril angle compared to than in normal wood, resulting in a lower stiffness in the longitudinal direction, and a higher stiffness in the radial-tangential plane (Brancheriau et al. 2012b;Gardiner et al. 2014). The propagation of ultrasonic waves in wood is therefore affected by the presence of compression wood, resulting in higher velocities in this areas [START_REF] Saadat-Nia | Ultrasonic wave parameter changes during propagation through poplar and spruce reaction wood[END_REF]Brancheriau et al. 2012b).

Other factor influencing the reconstruction is the variation in the wood moisture content along the tree cross-section. For example, in green softwoods, the sapwood presents generally a higher moisture content than the heartwood [START_REF] Ross | Wood handbook: wood as an engineering material[END_REF]. These variations lead to changes in the velocities of the ultrasonic waves, as they influence the mechanical parameters, and that are not considered in the proposed model.

Synthesis

This chapter presented an approach to solve the inverse problem for tree imaging, using a method based on ultrasound travel-time tomography and adapted to the anisotropy of the wood material. The proposed iterative method 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. This method resulted in the mapping of the specific elastic constants using nonlinear regression. The proposed inversion executed fast considering in-situ testing.

Numerical results obtained with four different configurations were presented:

healthy, centered and off-centered defects and a gradient of properties. The method approximated correctly the ray curvature. Defective regions identification and quantification proved to be more precise than the straight-ray inversion (FBP method).

Centered defect identification was more accurate than the cases with an off-centered defect and a properties gradient. Under the proposed geometry, the minimum size of a defect to be detected corresponded to a ratio defect diameter to disk diameter of 17% for the centric case and 23% for the eccentric case. The more suitable image for evaluation would be the map of ER values that presented less variations.

The experimental validation used two wood species, with configurations including a healthy case, a centered and an off-centered defect. To virtually increase the total number of sensors used, a sinogram interpolation was performed. Images obtained with the proposed method resulted in a more precise defect location compared to the FBP approach, as observed using a thresholding segmentation. For the ER images, threshold values of 50% and 70% resulted in the best compromise between sensitivity (true positive rate) and fallout (false positive rate) for segmentation, with a sensitivity from 60% to 98% and a fallout from 4% to 26%. For the healthy case in the pine sample, the presence of compression wood and juvenile wood resulted in variations on the computed image, associated to changes in density and stiffness for these types of wood compared to normal wood.

CONCLUSIONS AND PERSPECTIVES

Considering nondestructive evaluation of standing trees in urban areas, this work focused on the development of a methodology for the reconstruction of 2D ultrasound computed tomography adapted to the wood complexity as material. Aiming to increase the image quality compared to previous tomographic approaches, three main parts were considered: the definition of an electroacoustic measurement system adapted to the in-situ measurement of standing trees, the analysis of the factors influencing the ultrasonic waves propagation in wood by a numerical modeling (direct problem) and the implementation of a reconstruction method adapted to wood anisotropy (inverse problem).

To perform ultrasonic non-destructive testing of standing trees, the time-of-flight The anisotropy of wood in the radial-tangential plane affects the wave velocity depending on the direction of propagation, as observed using the Christoffel equation.

Moreover, from a sensitivity analysis of this equation, it was determined that the order of influence of the mechanical parameters on the velocity variation, from largest to smallest was Young's moduli in tangential direction ET and in radial direction ER, the Poisson's ratio νRT and the shear modulus GRT. Wavefronts propagating in wood present deformed shapes, compared to the spherical wavefronts obtained for isotropic media, leading to trajectories from the emitter to the receivers that do not follow a straight path for all directions of propagation. For the image reconstruction process, this effect was addressed, as the classic inversion methods considering straight rays deliver a biased image.

A ray-tracing approximation using the wavefront construction method was implemented to simulate the wave propagation in orthotropic media, therefore allowing the The raytracing results were compared with a numerical approach, in this case, the FEM method. Unlike raytracing, this method considers second-order phenomena such as refraction on the defect, mode conversion, and dispersion. Both models resulted in similar TOF estimation; however, raytracing may be considered as a less complex technique, resulting in lower processing times. Experimental validation was performed, using two different tree species and simulating the presence of defects by drilling holes of different sizes and positions. The shape of TOF curves computed using the raytracing algorithm and those obtained from the experiments with the wood samples were concordant, showing that the proposed model simplification is admissible.

Using the raytracing model and the information provided about TOF and ray paths, an algebraic inversion method was presented. The method resulted in a 2D mapping of the inner elastic parameters of wood, considering the material anisotropy. For every pixel, a slowness function was obtained depending on the propagation angle. To linearize the Christoffel equation, a polynomial approximation was proposed. Considering the slowness function, it was possible to estimate the associated mechanical parameters via a nonlinear regression. From an initial guess of straight-line trajectories, an iterative process performed the algebraic inversion and a successive raytracing simulation, until a stop criterion related to the TOF and trajectories variation was attained.

To evaluate the performance of the proposed method, a numerical and experimental validation were performed. For the numerical validation, four configurations were tested, including a healthy case, a centered defect, an off-centered defect and finally, a case corresponding to a gradient of mechanical properties. The proposed method approached accurately the ray curvature, resulting in images with a more detailed representation of the inner state compared to the filtered backprojection (FBP) method (isotropic assumption).

The centered defect case resulted in a more accurate identification than the cases with an off-centered defect and a properties gradient. For the experimental validation, two tree species were tested. The tested configurations included a healthy case, a centered defect, and an off-centered defect. Compared to the FBP images, the proposed method resulted in a more accurate defect identification, adapting the curved rays to the defect presence and delivering a parametric image more suitable for the diagnostic process.

Variability factors in the tree mechanical properties may be considered to obtain a better representation of the inner state using ultrasound measurements. For example, changes in the moisture content above the fiber saturation point (FSP) along the wood section have an impact on the propagation of elastic waves [START_REF] Sakai | Effect of moisture content on ultrasonic velocity and attenuation in woods[END_REF][START_REF] Unterwieser | Influence of moisture content of wood on sound velocity and dynamic MOE of natural frequency-and ultrasonic runtime measurement[END_REF][START_REF] Yamasaki | Influence of moisture content on estimating Young's modulus of full-scale timber using stress wave velocity[END_REF]. Above the FSP, the cell wall is saturated with water and microscopic vacant spaces (cell lumina) start keeping the free water. The moisture content variations affect the velocity in two ways: one is the intrinsic effect on the elastic moduli and the other is the effect on the wood density. Similarly, changes in wood density should be addressed. These variations are due to anatomical characteristics such as the ratio of earlywood to latewood and heartwood to sapwood. These interactions may be approached by numerical and experimental settings, and then included in the formulation of the direct problem.

In the longer term, the idea of considering all the physical phenomena involved could make it possible to achieve a much better characterization of the wood material.

Today, the methods developed in geophysics offer great possibilities. Works from this community already inspire the raytracing approach used in this thesis. A potential method to address this strategy is the Full Waveform Imaging (FWI) method (Grechka and Wapenaar 2014; [START_REF] Bernard | Ultrasonic computed tomography based on full-waveform inversion for bone quantitative imaging[END_REF]). The term "full" refers to the use of the full-time series, i.e. without having to discard potentially useful hypothesis in the records. This method relies in an iterative process, reducing the difference between the full recorded data and synthetic waveforms from a numerical model of the wave propagation. In this case, the direct problem of the wave propagation should be efficient, as it is used numerous times during optimization. It was mainly developed in the oil industry and in seismology to obtain maps of the celerity of seismic waves, but recently FWI methods have also started to be applied in medical ultrasound imaging [START_REF] Wang | Waveform inversion with source encoding for breast sound speed reconstruction in ultrasound computed tomography[END_REF][START_REF] Pérez-Liva | Time domain reconstruction of sound speed and attenuation in ultrasound computed tomography using full wave inversiona)[END_REF].

USCT systems aiming for 3D imaging, commonly perform a stacking of 2D images obtained at different heights and apply interpolation [START_REF] Goncharsky | Inverse problems of 3D ultrasonic tomography with complete and incomplete range data[END_REF]. A full 3D setup would require a distribution of transducers around a trunk section covering different heights. In that case, the inverse problem is very computer intensive, as it should be solved using a large amount of data by considering all the possible transmitter-receiver combinations. To model the wave propagation in this case, methods should consider the cylindrical orthotropic condition of wood in the three directions. The resulting wavefronts will be affected, considering that the wave velocities in longitudinal direction are larger than the velocities obtained in the radial or tangential directions.

decay along the wood cross-section. They are cheap and easy to use, but it is an invasive method causing a wound in the tree. Using this device in trees with internal decay, can result in breaking the barrier zone created by the tree to isolate the decay.

3. Penetrometer: These devices record the resistance opposed to a probe driven into the wood. The basic principle is that in the presence of decay, wood density decreases, and wood hardness and drilling resistance decline. One of the most know devices is the Resistograph (Rinntech). As the drill penetrates the wood, resistance to the pressure of the drill is measured and recorded, and the resulting profile is simultaneously printed and immediately available. Drill bit diameter is 1/16 inch (1.59 mm) and the drilling depth is around 50 cm. In the presence of early to early-intermediate stages of decay, the penetrometers may underestimate the extension of decayed areas. Compared to the previous approaches, this device is expensive to purchase and maintain.

4. Sonic and ultrasonic detectors: Stress-wave timers rely on the principle that the speed of sound of wave propagating through the wood is affect by the presence of defective areas. A hammer blow is used to generate the sound waves. Decay regions correspond to slow propagation zones, then the wave transit time is longer than in healthy wood. A limitation to this method is the inability to locate the specific position of the defect or the extent of the internal defect. Similarly, ultrasonic devices work on the same principle, but in this case the waves are produced by an ultrasonic pulse, providing more precise transit time measurements. Limitations are similar to the stress-wave timers, and additionally, they cannot be used in large trunks, due to the attenuation of the ultrasonic waves.

5. Electrical conductivity meters: it consists on electrical resistance (ER) measurements using an electrode inserted in a drilled hole of 3 mm in the trunk to a depth of 30 cm. For the decayed areas, the ER value drops substantially. Some limitations include a difficult interpretation of the results when the moisture content of the wood is below the fiber saturation point or when the drill hole is filled with water, resulting in overestimations of the decay area.

6. Fractometer: This instrument determines the fracture moment and angle of failure of a wood sample (core samples of 5 mm), by placing it in a clamping device and applying an increasing force up to the point of failure. A decrease in fracture moment or an increase in fracture angle could indicate the presence of decay. Limitations to this method includes a lack of reference measurements for diagnostic and the consideration that this is a more invasive method than the previous approaches presented.

  Imagerie ultrasonore pour la gestion durable et raisonnée des arbres sur pied en milieu urbain Les arbres urbains jouent un rôle écologique, sanitaire et esthétique majeur dans les villes modernes. Ces bénéfices incluent entre autres la réduction de la pollution de l'air, l'abri et nourriture à plusieurs espèces d'oiseaux et petits animaux, la protection contre la chaleur, et la mise en valeur du paysage. Cependant, certains risques sont associés aux arbres urbains : la chute d'un arbre (ou d'une branche) peut causer des dommages aux passants, aux véhicules et aux habitations. L'évaluation des risques associés aux arbres dans les villes est souvent à l'initiative des opérateurs urbains. Le diagnostic est essentiellement visuel, alors que l'ampleur des dégâts internes et le danger associé ne peuvent pas être évalués avec précision par la seule observation. Fréquemment, les arbres présentant une décomposition interne semblent sains. Répondre à la gestion des plantations urbaines nécessite donc l'utilisation de méthodes de diagnostic interne des arbres in situ. Des techniques modernes peuvent être utilisées pour minimiser les risques associés à la chute d'arbres ; des progrès significatifs dans le matériel de diagnostic ont été réalisés, ainsi que dans les formules et les directives d'évaluation des arbres dangereux. Nous avons eu l'occasion de participer à un processus d'inspection d'arbres in-situ, afin d'évaluer leur niveau de dangerosité, avec un expert arboricole. Les arbres concernés appartenaient au genre Platanus, situés tout au long d'une autoroute près de la localité de Lodève, diagnostic demandé par le conseil départemental de l'Hérault (France). Le diagnostic est divisé en deux phases : l'évaluation visuelle et l'évaluation spécialisée avec instruments. En fonction de l'environnement local et de son importance patrimoniale, le diagnostic d'un arbre peut impliquer les deux phases ou seulement l'évaluation visuelle. Par exemple, les arbres que nous avons inspectés, et que nous pourrions définir comme arbres d'ornement, sont normalement évalués visuellement, à moins qu'ils ne présentent des signes importants de dégradation. Des techniques instrumentées tels que le perçage sont utilisés dans ce cas. Le diagnostic doit être ici rapide dans son exécution. Certains arbres urbains ayant une importance patrimoniale élevée peuvent nécessiter une analyse plus approfondie et plus longue à mettre en oeuvre. Les techniques telles que l'imagerie acoustique et ultrasonore deviennent alors des outils indispensables.

  Il faut d'abord considérer une image divisée en N pixels traversés par les rayons ultrasonores. L'objectif est d'estimer un paramètre local pour chaque pixel (la lenteur α) qui permet par sommation le long du rayon de déterminer la mesure du TOF entre émetteur et récepteur. Pour un couple émetteur-récepteur, dont la trajectoire m traverse les pixels k, le temps de vol tm peut s'écrire comme suit : 𝑡 𝑚 = ∑ 𝑙 𝑚𝑘 𝛼 𝑘 𝑘 𝑎𝑙𝑜𝑛𝑔 𝑚

Figure 1 :

 1 Figure 1: Fall of a tree affecting a vehicle (source: https://www.paris.fr/arbres).

Figure 2 :

 2 Figure 2: Some methods for nondestructive evaluation of wood: (a) Drilling measure (Source: http://www.resistograph.com/), (b) Thermography imaging (Catena and Catena 2008).

Figure 3 :

 3 Figure 3: Comparison of acoustic tomogram (left) and photograph from cross-section for Red Oak tree (Wang 2007).

Figure 4 :

 4 Figure 4: Standing tree evaluation using the drilling technique.

Figure 5 :

 5 Figure 5: Example of use of acoustic hammer: Hitman ST300 for standing tree evaluation (source: https://www.fibre-gen.com/hitman-st300)

Figure 6 :

 6 Figure 6: Picus Sonic Tomograph from Argus Electronic Gmbh.

Figure 7 :

 7 Figure 7: Example of a combination of electrical and acoustic tomography for different stages of decay evaluation (source: http://www.argus-electronic.de/en/tree-inspection/products/picustreetronic-r-predicting-future-decay).

Figure 8 : 1 .

 81 Figure 8: Acoustic testing of a standing tree in Bogota using an Arbotom device.

  presented MIGRATOM, a software to perform tomographic image reconstruction using measurements of travel time or attenuation, that considers curved rays due to velocity contrasts. Their common application is in geophysics imaging problems, such as fracture detection, fluid monitoring and qualitative stress evaluation. The aim of MIGRATOM is to reconstruct 2D velocity maps, performing an iterative process. An initial velocity model is modified by repeated cycles of three steps: (i) forward computation of travel times, (ii) estimation of the residuals and (iii) application of velocity corrections. For the forward model, this software uses a numerical two-dimensional wavefront propagator based on the Huygens' principle. It considers each point of the wavefronts as an instantaneous point source of wave energy, advancing in a series of time steps, as a circular wavelet with a radius proportional to the local velocity; rays corresponds to the trajectory of an individual point on the wavefront.

Figure 9 :

 9 Figure 9: Wood anatomy. (a) Transverse section of Quercus alba trunk. From outside to inside: out bark (ob), inner bark (ib), vascular cambium (vc), sapwood, heartwood, and pith (p) (Ross 2010). (b) Planes of wood: cross or transverse section, radial section and tangential section (Kollmann and Côté 2012). (c) General structure (A-B) and microscopic structure (C-D) for softwoods (left) and hardwoods (right) (Ross 2010).

Figure 10 :

 10 Figure 10: Wood cell wall (a) simplified structure of a tracheid cell wall showing the middle lamella (ML), primary wall (P), secondary cell wall layers (S1, S2, S3) and warty (W) layer. (b) Micrograph to show the location of the middle lamella (ML) and the middle lamella cell corner (MLcc) (Ek et al. 2016).

  Figure 11 presents the variation of the elastic modulus EL for two species for multiple moisture content values. Before 5%, the elastic constant increases lightly, from 8 to 22% the value decreases almost linearly and after 30% (fiber saturation point) the elastic modulus stabilizes.

Figure 11 :

 11 Figure 11: EL evolution for two wood species (spruce and oak) depending on the moisture content (Guitard 1987).

Figure 12 :

 12 Figure 12: Compression waves (polarization vector parallel to the direction of propagation) and shear waves (polarization vector perpendicular to the direction of propagation).

Figure 14 :

 14 Figure14: Velocity values for all species in function of the angle between the vector normal to the wavefront and the radial direction(Espinosa et al. 2017a).

Figure 15 :

 15 Figure 15: Velocity variations (in percentage) induced by changing each mechanical parameter for Sitka spruce according to the angle between the vector normal to the wavefront and the radial

Figure 17 :

 17 Figure 17: Velocity variation (in percentage) by setting νRT and GRT parameters to zero for (a) Spruce and (b) Oak species (Espinosa et al. 2017a).

Figure 18 :

 18 Figure18: Platanus standing tree tested[START_REF] Espinosa | Accuracy on the Time-of-Flight Estimation for Ultrasonic Waves Applied to Non-Destructive Evaluation of Standing Trees: A Comparative Experimental Study[END_REF]).

  generator and oscilloscope corresponded to a Picoscope 2000 (emission sample rate 1 MHz, reception sample rate 4 MHz), with an interface to a personal computer for data acquisition. Input amplifier reference was FLC Electronics Single Channel High Voltage Linear Amplifier A800 (bandwidth DC to 250 kHz, 40 dB amplification). Output amplifier was Physical Acoustics Corporation AE2A/AE5A wide bandwidth AE amplifier (bandwidth up to 2 MHz, internal 40 dB preamplifier).

Figure 19 :

 19 Figure19: Ultrasonic chain for measurements[START_REF] Espinosa | Accuracy on the Time-of-Flight Estimation for Ultrasonic Waves Applied to Non-Destructive Evaluation of Standing Trees: A Comparative Experimental Study[END_REF]).

  The short duration rectangular pulse, pulse train, and half-Gaussian pulse present a fast-impulsive start, resulting in a large band frequency response, with several resonant lobes in the case of the pulse train and a soft power decay for the half-Gaussian pulse. The Gaussian pulse and chirp signal have a sinusoidal shape, multiplied by a Gaussian window, resulting in a concentrated power spectrum around a central frequency (resonant peaks for the sensors), with a narrower bandwidth for the chirp signal. Parameters fixed for the signals are presented in Table5. Peak voltage for all signals was set to 2V (maximum for signal generator). Signals repetition period was fixed to T=8 ms.

Figure 20 :

 20 Figure 20: Signal shapes tested: (a) impulse, (b) pulse train, (c) Gaussian pulse, (d) half-Gaussian pulse and (e) chirp signal (Espinosa et al. 2018).

Figure 21

 21 Figure 21 presents the root mean square voltage (RMS) mean and standard deviation values, for the received signals, for all the experiment configurations. Correspondingly, Table 6 summarizes the RMS values for the five signals, sorting by the RMS mean value in descending order. Except for the pulse train signal, almost all configurations that used sensor R6α resulted in larger RMS values than the R3α counterpart. Receiver angles with larger RMS values were those located at 90° and 135°. For the R3α sensor, the signals presenting an impulsive behavior (pulse train, half-Gaussian pulse, and impulse) resulted in more energetic received signals. Chirp signal received for both cases ranked in the last positions.

Figure 21 :

 21 Figure 21: Mean values for (a) RMS, and (b) SNR, for all configurations. Error bars present the standard deviation ±σ (Espinosa et al. 2018).

Figure 22 :

 22 Figure 22: Mean values for SNR for all configurations. Error bars present the standard deviation ±σ (Espinosa et al. 2018).

Figure 23 andFigure 23 :Figure 24 :

 232324 Figure 23 and Figure 24 present first the input and output signals on time domain, then their frequency spectrum and finally the input and output spectrograms, for sensors R3α and R6α, respectively.

  technique and AIC methods, for all signals except chirp. Mean TOF values ranged between 65 μs to 143 μs. The difference in mean values estimated with the two methods was always inferior to 1.4 μs. Standard deviation ranged between 0.3 μs to 8.8 μs for the threshold method and 0.2 μs to 6.7 μs for the AIC method.

Figure 25 :

 25 Figure 25: Mean TOF values using the Threshold method (left) and AIC method (right) for R3α.

Figure 27 :

 27 Figure 27: Mean TOF values using the Threshold method (left) and AIC method (right) for R6α.Error bars present ±σ[START_REF] Espinosa | Accuracy on the Time-of-Flight Estimation for Ultrasonic Waves Applied to Non-Destructive Evaluation of Standing Trees: A Comparative Experimental Study[END_REF]).

Figure 28 :

 28 Figure 28: Relative standard deviation for TOF values using the Threshold method (up) and AIC method (down) for R6α (Espinosa et al. 2018).

Figure 29 :

 29 Figure 29: Mean TOF values for the chirp signal using the Threshold, AIC and cross-correlation methods for R3α (left) and R6α (right). Error bars present ±σ (Espinosa et al. 2018).

Figure 30 :

 30 Figure 30: Relative standard deviation for TOF values for the chirp signal using the Threshold, AIC and cross-correlation methods for R3α (up) and R6α (down) (Espinosa et al. 2018).

Figure 31 :

 31 Figure 31: AIC-Impulse and Chirp-Cross-correlation comparison for TOF relative standard deviation values (Espinosa et al. 2018).

Figure 32 :

 32 Figure 32: Chirp dispersion effect for the case of sensor R6α and the receiver located at 135°: instantaneous frequency from Gabor transform for input and output signals.

Figure 33 :

 33 Figure 33: Wavefront construction principle. Arrows correspond to the direction vector normal to the wavefront (Espinosa et al. 2019).

Figure 34 :

 34 Figure 34: Wavefront for (a) isotropic (PVC) and (b) orthotropic (Spruce) materials. For visualization purposes, not all wavefronts are shown.

Figure 35 :

 35 Figure 35: Wave ray paths from the emitter to multiple receivers in (a) isotropic (PVC) and (b) orthotropic materials (Spruce).

Figure 36 :

 36 Figure 36: Distances from receivers to the emitter in isotropic and orthotropic cases.

Figure 37 :

 37 Figure 37: Relation between curved ray path distances with straight line distance for all species.

Figure 38 :

 38 Figure 38: TOF measurements (a) for all species and (b) mean values for softwood, hardwood and isotropic.

Figure 39 :

 39 Figure 39: Estimated velocity value using the real path distance (continuous line) and the straightline distance (line with dots) for (a) softwoods and (b) hardwoods.

Figure 40 :

 40 Figure 40: Relative error for velocity estimation using straight-line distances for all species.

Figure 41 :

 41 Figure 41: Wave fronts using defects with diameters 2 (a-b), 4 (c-d) and 6 cm (e-f), for orthotropic (left) and isotropic (right).

Figure 42 :

 42 Figure42presents the TOF values for each defect size, for both cases. The changes in TOF values were larger for the orthotropic case. For example, the TOF increased 35μs for sensor 9 in orthotropic case comparing the healthy case and the bigger defect case; in the case of isotropic, the variation was near to 6 µs.

  Figure 43: (a-b) Velocity values for orthotropic (left) and isotropic (right), and (c-d) ratio of velocity reduction V/Vref for the sensors 6, 7, 8 and 9.

Figure 44 :

 44 Figure 44: Wave fronts using defects displaced 8 cm to the left (a-b) and displaced 8 cm to the right (c-d) in horizontal axis, for orthotropic (left) and isotropic (right).

Figure 45 :Figure 46 :

 4546 Figure 45: TOF variation in all receivers for the defect changing in horizontal direction for (a) Orthotropic and (b) isotropic.

Figure 47 :

 47 Figure 47: Wavefronts using defects located in the center of the trunk (a-b) and displaced 10 cm in the horizontal axis (c-d), for orthotropic (left) and isotropic (right).

Figure 48

 48 Figure 48 presents the TOF variation for each defect position, for both cases. For the anisotropic material, as the defect position got closer to the receivers, fewer sensors

Figure 48 :Figure 49 :

 4849 Figure 48: TOF variation in all receivers changing the defect position for (a) Orthotropic and (b) isotropic

Figure 50 :

 50 Figure 50: Wave fronts using defects located in the diagonal positions (a) [-8cm, 8cm], (b) [8cm, 8cm], (c) [8cm, -8cm] and (d) [-8cm, -8cm].

Figure 52 :

 52 Figure 52: Wavefronts using defects reducing the velocity by 30% (a-b) and by 70% (c-d), for orthotropic (left) and isotropic (right).

Figure 53 :Figure 54

 5354 Figure 53: TOF variation in all receivers changing the decay percentage for (a) Orthotropic and (b) isotropic

Figure 55 :

 55 Figure 55: Convergence of the mean error depending on the total number of elements

Figure 56 :

 56 Figure 56: Wavefronts comparison for healthy case (a) at 60 µs, (b) at 150 µs; for a defect case (c) at 120 µs, (d) at 180 µs. In black: Raytracing wavefront.

Figure 57 :

 57 Figure 57: TOF comparison. Healthy case: (a) Softwood species and (b) Hardwood species. Defect case: (c) Softwood species and (d) Hardwood species. Dashed lines correspond to raytracing results and continuous lines to FEM results.

Figure 58 :

 58 Figure 58: Mean absolute difference between Raytracing and FEM estimations of TOF. In blue healthy experiments; in red defect experiments.

Figure 59 :

 59 Figure 59: Oak disks tested: (a) Healthy, (b-d) centered defects with diameters 2.9, 5.1 and 7.6 cm.

Figure 60 :

 60 Figure 60: PVC disk for comparison with the isotropic case. Disk diameter of 30 cm and a thickness of 3 cm (Espinosa et al. 2019).

Figure 61 :

 61 Figure 61: Wavefronts (in black) and ray paths (in red) from the transmitter to multiple receivers in (a) isotropic (PVC) and (b) orthotropic materials (Ash). Comparison between TOF for (c) simulated and (d) experimental (Espinosa et al. 2019).

  The shape of TOF curves for the simulated data and the experiments were corresponding. For the simulated values, minimum TOF was 32 µs and maximum values were 120 µs for the healthy reference case, and 133 µs, 146 µs and 161 µs for the three defect sizes respectively; for the experimental values the minimum TOF ranged between 35 µs and 53 µs, and the maximum values were 150 µs for the healthy case, and 168 µs, 209 µs, and 225 µs as the defect size increased.

Figure 62 :

 62 Figure 62: Wavefronts (in black) and ray paths (in red) in (a) healthy case, (b) centered defect of 2.9 cm, (c) centered defect 5.1 cm and (d) centered defect of 7.9 cm. Comparison between TOF for (e) simulated and (f) experimental (Espinosa et al. 2019).

Figure 63 :

 63 Figure 63: TOF difference between centered defects and the healthy case for (a) simulated data and (b) experimental data (Espinosa et al. 2019).

Figure 64 :

 64 Figure 64: Wavefronts (in black) and ray paths (in red) in (a) healthy case, (b) off-centered defect of 2.9 cm, (c) off-centered defect 5.1 cm and (d) off-centered defect of 7.9 cm. Comparison between TOF for (e) simulated and (f) experimental (Espinosa et al. 2019).

Figure 65 :

 65 Figure 65: TOF difference between defective and healthy cases for (a) simulated data and (b) experimental data (Espinosa et al. 2019).

Figure 66 :

 66 Figure 66: Off-centered case located at the bottom: (a) Wavefronts (in black) and ray paths (in red), comparison between TOF for (b) simulated and (c) experimental. Off-centered case located at the top: (d) Wavefronts (in black) and ray paths (in red), comparison between TOF for (e) simulated and (f) experimental (Espinosa et al. 2019).

Figure 68 :

 68 Figure 68: Coefficient of determination R 2 for several polynomial degrees to linearize the Christoffel equation in the case of Oak (Table11).

Figure 69 :

 69 Figure 69: Flow-chart for the proposed reconstruction method: feedback using a raytracing approach

  Figure 70 presents the reference trajectories and the trajectories estimated with the proposed method. Rays obtained with the proposed methodology reproduced correctly the curved behavior expected from anisotropic materials. Figure 71 presents the reconstructed images and the corresponding profiles. FBP method resulted in a gradient of velocity values, with higher values in the disk center, ranging from 1362 m/s to 1912 m/s. With the proposed method, the image representing Young's modulus in the radial direction presented a uniform distribution of reconstructed values around the real value of 1180 MPa, between 1000 MPa

Figure 70 :

 70 Figure 70: Ray paths for a healthy case: (a) reference trajectories obtained with the ray-tracing approach, (b) estimated trajectories obtained with the proposed inversion method.

Figure 72 :

 72 Figure 72: Ray paths for a centered defect case: (a) reference trajectories obtained with the raytracing approach, (b) estimated trajectories obtained with the proposed inversion method.

Figure 74 :

 74 Figure 74: Ray paths for an off-centered defect case: (a) reference trajectories obtained with the ray-tracing approach, (b) estimated trajectories obtained with the proposed inversion method.

Figure 76 :

 76 Figure 76: Ray paths for a case with a gradient of properties: (a) reference trajectories obtained with the ray-tracing approach, (b) estimated trajectories obtained with the proposed inversion method.

  Minimum diameter to identify a (a) centered defect and (b) off-centered defect. Dashed lines indicate the expected profile.

  Figure 79: Relative error for the nonlinear regression of the elastic parameters: (a) with the four parameters, (b) fixing the Poisson's ratio.

Figure 80 :

 80 Figure 80: Experimental setting for the validation of the proposed inversion method for the case of a healthy oak disk.

Figure 81 :

 81 Figure 81: Sinogram interpolation for the case of a healthy oak disk, passing from 16 sensors to 32 sensors.

  . After the inversion procedure, the resulting rays are presented in Figure82b. The two reconstructed images are presented in Figures82c and 82drespectively: the ER parametric image and the reference FBP image. Horizontal profiles (passing from sensor 13 to sensor 5) from the obtained images are presented in Figures82e and 82f. In the case of the ER image, the mean value was 1046 MPa (σ=205 MPa). The mean velocity in the FBP image was 1673 m/s (σ=688 m/s). Both images presented approximated flat surfaces, with a reduced marge of variation for the ER image.

Figure 82 :

 82 Figure 82: (a) Oak healthy trunk, (b) rays obtained after the reconstruction procedure, (c) resulting parametric image, (d) FBP reconstruction for comparison, and horizontal profiles from (e) the parametric image and (f) the FBP image.

Figure 83 :

 83 Figure 83: (a) Oak with a centered defect, (b) rays obtained after the reconstruction procedure, (c) resulting parametric image, (d) FBP reconstruction for comparison, and horizontal profiles from (e) the parametric image and (f) the FBP image. Dashed lines in the profiles (e-f) indicate the defect location.

Figure 84 :Figure 85 :

 8485 Figure 84: Defect segmentation using different threshold values for the FBP images (top) and the ER images (bottom) for the case of Oak with a centered defect. White areas show the pixels classified as defective. Dashed circle represents the real boundary of the defect.

Figure 86 :

 86 Figure 86: Defect segmentation using different threshold values for the FBP images (top) and the ER images (bottom) for the case of Oak with an off-centered defect. White areas show the pixels classified as defective. Dashed circle represents the real boundary of the defect.

  resulting rays after the inversion procedure are presented in Figure87b. The ER parametric image and the reference FBP image are presented in Figures87c and 87d. Two horizontal profiles (sensor 13 to sensor 5) are shown in Figures87e and 87f. In the case of the ER image, the mean value was 1258 MPa (σ=126 MPa). The mean velocity in the FBP image was 1551 m/s (σ=315 m/s).

Figure 87 :

 87 Figure 87: (a) Oak healthy trunk, (b) rays obtained after the reconstruction procedure, (c) resulting parametric image, (d) FBP reconstruction for comparison, and horizontal profiles from (e) the parametric image and (f) the FBP image.

FigureFigure 88 :

 88 Figure 88a presents the pine disk with a centered defect. The trajectories were modified by the presence of the hollow region, as observed in Figure 88b. The proposed method allowed a more precise defect identification, as observed in the reconstructed image in Figure 88c compare to the FBP reconstruction in Figure 88d. Their corresponding horizontal profiles are shown in Figures 88e and 88f respectively. For the ER image, the

Figure 89 :

 89 Figure 89: Defect segmentation using different threshold values for the FBP images (top) and the ER images (bottom) for the case of Pine with a centered defect. White areas show the pixels classified as defective. Dashed circle represents the real boundary of the defect.

Figure 90 :

 90 Figure 90: (a) Oak with an off-centered defect, (b) rays obtained after the reconstruction procedure, (c) resulting parametric image, (d) FBP reconstruction for comparison, and horizontal profiles from (e) the parametric image and (f) the FBP image. Dashed lines in the profiles (e-f) indicate the defect location.

Figure 91 :

 91 Figure 91: Defect segmentation using different threshold values for the FBP images (top) and the ER images (bottom) for the case of Pine with an off-centered defect. White areas show the pixels classified as defective. Dashed circle represents the real boundary of the defect.

(

  TOF) measurement precision is a critical issue. An experimental study comparing several signal shapes and TOF detection methods was carried out, for setting up the ultrasonic chain of measurement having in mind the subsequent testing. Impulsive and encoded signals were tested, for a total of five signals, combined with three different methods for TOF estimation: Threshold, AIC method, and cross-correlation. In situ testing was performed with two different transducers with resonant frequencies at 36 kHz and 60 kHz and 4 different receiver positions around the tree. From all these experiences, it was settled that the one presenting less variation on the TOF measurements corresponded to the combination of an encoded excitation signal, such as chirp signal, with cross-correlation to measure the time delay. Encoded signals such as the chirp signal should be adjusted to the transducer response. Thus, the received signals concentrated energy in frequency bands around the resonant frequency of sensors.

  estimation of TOF values and ray paths. Several considerations were established for the numerical simulation, like defining a 2D geometry and only considering pure compression waves for the estimation of the first arrivals. Simulations included the presence of defects with different diameters and positions. As expected, the estimated trajectories presented a curvature, resulting in larger velocity values in the radial direction. Regarding the defects, those located in the center of the trunk presented larger TOF variations compare to defects located in off-centered positions, thus off-centered defects could be more difficult to determine and characterize by tomographic inversion.

  Avec km l'angle de propagation par rapport à la direction radiale du repère local d'orthotropie du pixel k. En conséquence, la solution du problème inverse correspond à la défaut excentré et un gradient de propriétés quelle que soit la méthode utilisée.Cette méthode ne nécessite que quelques minutes de calcul avec un ordinateur portable du commerce ; elle est donc adaptée au contexte des tests in situ.

	RESUMÉ
	𝛼 𝑘 = 𝛽 5,𝑘 𝜃 𝑘,𝑚 5 + 𝛽 4,𝑘 𝜃 𝑘,𝑚 4 + 𝛽 3,𝑘 𝜃 𝑘,𝑚 3 + 𝛽 2,𝑘 𝜃 𝑘,𝑚 2 + 𝛽 1,𝑘 𝜃 𝑘,𝑚 1 + 𝛽 0,𝑘 𝜃 𝑘,𝑚 0
	recherche des coefficients polynomiaux β pour tous les pixels. Pour résoudre l'ensemble La méthode d'inversion proposée a ensuite été testée expérimentalement sur des
	des équations linéaires, la méthode SIRT (technique de reconstruction itérative simultanée) échantillons de bois de deux essences (chêne, pin) pour trois configurations différentes : un
	est utilisée. Afin d'augmenter le nombre d'équations (résolution spatiale) une interpolation cas sain, avec défaut centré et avec un défaut excentré. Comme pour la validation
	des mesures du TOF a été effectuée en considérant des capteurs virtuels situés entre ceux numérique, la méthode développée a permis d'obtenir une représentation plus précise des
	d'origine (interpolation du sinogramme). défauts par rapport à une reconstruction classique par rayons droits, en particulier dans le
	cas de défauts centrés.
	Le processus itératif démarre avec des trajectoires rectilignes. Avec ces trajectoires,
	formulation matricielle doit être adaptée à la dépendance de la lenteur vis-à-vis de l'angle
	de propagation (linéarisation de l'équation de Christoffel), et deuxièmement, l'estimation
	des vraies trajectoires (inconnues a priori) est affectée par l'anisotropie du bois et la
	présence de défaut.
	xii

Avec lmk la longueur d'un segment du rayon de propagation. La résolution de toutes les équations des couples émetteurs-récepteurs entraîne la reconstruction du paramètre interne αk (problème inverse). Néanmoins, deux problèmes se posent : premièrement, la L'équation de Christoffel est linéarisée par une approximation polynomiale (5 ème degré). La lenteur pour le pixel k traversé par un rayon m est alors exprimée comme suit : il est possible d'effectuer une première inversion conduisant aux coefficients polynomiaux de lenteur β pour tous les pixels. A partir de ces coefficients, le problème direct est résolu pour obtenir une estimation des TOF et de nouvelles trajectoires. Ce processus est répété jusqu'à ce que la différence entre les valeurs des TOF par rapport aux valeurs réellement obtenues soit minimisée et que les formes des trajectoires calculées ne varient plus. Finalement, l'approximation polynomiale obtenue pour chaque pixel est utilisée pour estimer les valeurs des constantes élastiques via une régression non linéaire par la méthode des moindres carrés. Le résultat est une image paramétrique des modules spécifiques ER, ET et GRT en considérant un coefficient de Poisson νRT constant (permettant ainsi d'obtenir des paramètres stables par régression non linéaire). Une validation numérique de cette méthode d'inversion a été réalisée. Quatre configurations numériques ont été testées représentant des cas réels généralement rencontrés sur le terrain : (1) un tronc dans un état sain, (2) avec un trou centré, (3) avec un trou excentré et (4) avec un défaut excentré caractérisé par un gradient de propriétés mécaniques. Ce dernier cas pourrait être associé à un stade précoce de dégradation du bois par une attaque de champignons. Les images reconstruites utilisant la méthode proposée ont été comparées à la méthode de reconstruction classique avec l'hypothèse d'isotropie (FBP, rayons droits). La comparaison des images obtenues a mis en évidence une identification et une quantification plus détaillées de l'état interne du tronc avec la méthode proposée. La méthode proposée a correctement reconstruit la courbure des rayons. Les zones de défaut sont mieux définies et plus visibles avec cette méthode qu'avec la méthode d'inversion classique par rayons droits. L'identification des défauts centrées était plus précise que les xiii cas avec un D'autres facteurs associés à la variabilité du bois peuvent être inclus dans la modélisation afin d'améliorer la fiabilité de la méthode. Par exemple, des variations de la teneur en humidité au-dessus du point de saturation des fibres ou des variations de la densité affectent la vitesse de propagation des ondes. Dans ce travail, nous avons considéré les ondes de compression en négligeant tous les autres phénomènes du second ordre (réfraction sur le défaut, conversion de mode, atténuation et dispersion). Une suite à ce travail serait de prendre en compte l'ensemble de ces phénomènes physiques afin de caractériser le bois plus finement. L'utilisation de la méthode « Full Waveform Imaging » (FWI) pourrait ainsi être envisagée. Le terme « full » fait référence à l'utilisation de la totalité du signal acoustique, c'est-à-dire sans rejeter des informations potentiellement utiles. Aujourd'hui, ces méthodes ne sont pas encore développées pour le cas du bois. Envisageant une implémentation 3D de la méthode de reconstruction, une installation de transducteurs à différentes hauteurs est nécessaire, avec l'adaptation du problème direct aux fronts d'onde se propageant dans la direction longitudinale, radiale et tangentielle. xv

1. INTRODUCTION 1.1. Background and motivation 1.1.1. Urban trees

	LIST OF ABBREVIATIONS
	L AIC	Wood longitudinal axis Akaike Information Criterion
	R ART	Wood radial axis Algebraic Reconstruction Technique
	T CODIT	Wood tangential axis Compartmentalization of Decay in Trees
	σ FBP	Stress Filtered Back Projection
	𝜖 FEM	Strain Finite Element Model
	C ONF	Rigidity matrix Office National des Forêts
	S PVC	Compliance matrix Polyvinyl Chloride
	E RMS	Young modulus Root Mean Square
	G SIRT	Shear modulus Simultaneous Iterative Reconstruction Technique
	ν SNR	Poisson's ratio Signal-to-Noise Ratio
	ρ TOF	Density Time-of-Flight
	p USCT	Wave polarization vector Ultrasound Computed Tomography
	n VTA	Wavefront propagation direction Visual Tree Assessment
	V	Wave phase velocity
	Γ	Christoffel tensor
	𝜃	Angle between the wave direction n and the direction of the radial axis
	Δt	Simulation time step
	𝛼	Wave slowness
	t	Wave transit time
	L	Length of a ray
	𝛽	Coefficients of polynomial regression
		xxv

Table 1 :

 1 Anisotropy in wood: reference values for hardwoods and softwoods[START_REF] Guitard | Mécanique du matériau bois et composites. Cepadues Editions Habermehl A, Ridder HW[END_REF]).

	Type	EL/ER	ER/ET	EL/ET	ρ [g/cm3]
	Hardwood	8	1.7	13.5	0.65
	Softwood	13	1.6	21	0.45

Table 2 :

 2 The shear modulus GRT ranged between 36 MPa and 319 MPa, and the Poisson's ratio νRT ranged between 0.38 and 0.70. Density ρ, ranged from 448 to 706 kg/m 3 . Selected Species and Corresponding Mechanical Parameters. ER and ET: Young's modulus in the radial and tangential directions respectively, GRT: shear modulus, νRT: Poisson's ratio, ρ: density. From published data (Ross 2010).

	Common	Scientific	ER/ET ER	ET	GRT	RT	ρ (kg/m 3 )
	Names	Names		(MPa)	(MPa)	(MPa)
	Douglas fir	Pseudotsuga	1.36	909	668	94	0.39 538
		menziesii

Table 2 are ranked according to the ratio ER/ET. The first half of this table corresponds to softwoods and the second half, to hardwoods. The radial (ER) and tangential (ET) Young's modulus ranged between 909 MPa to 2118 MPa for ER, and between 511 MPa to 1128 MPa for ET.

Table 3 :

 3

Maximum (Vmax)

and Minimum (Vmin) Velocity Values with their Associated Angles (θ).

Table 4 :

 4 Sensitivity values obtained with a variation of 10% for each mechanical parameter.

	Common Names	%V (ER)	%V (ET)	%V (GRT)	%V (νRT)
	Douglas fir	4.4	5.6	0.7	2.2
	Sitka spruce	4.4	5.6	0.3	2.3
	Longleaf pine	4.6	5.4	1.0	1.6
	Northern red oak	4.0	6.0	1.0	2.8
	Yellow poplar	3.5	6.5	0.8	3.9
	Sweetgum	3.7	6.3	1.2	3.1

Table 5 .

 5 Signal parameters. Ts presents the duration of signal portion. Fc indicates the central frequency of every signal. Fco indicates the cut-off frequency range (-3 dB points around central frequency). For the chirp signal, ΔF presents the bandwidth, around the central frequency.

	Signal	Expression		Parameters for R3α	Parameters for R6α
	0 if 𝑡 ≠ 0 Impulse 𝛿(𝑡) = { 1 if 𝑡 = 0	𝑇𝑠: 5 µs	𝑇𝑠: 5 µs
						𝐹𝑐𝑜: [0 90.159] kHz	𝐹𝑐𝑜: [0 90.159] kHz
	Pulse	For one period 𝑇𝑝 = 1/𝑓𝑐:	𝐹𝑐: 36 kHz	Fc: 60 kHz
	train
		𝑝(𝑡) = {	0 if t > 𝑇𝑝/2 1 if t < 𝑇𝑝/2	𝑇𝑠: 83 µs (3 Periods)	𝑇𝑠: 50 µs (3 Periods)
						𝐹𝑐𝑜: [29.39 40.20] kHz	𝐹𝑐𝑜: [49.75 68.01] kHz
	Gaussian	𝑟(𝑡) =				𝐹𝑐: 36 kHz	Fc: 60 kHz
	pulse	𝑠𝑖𝑛(2𝜋𝑓 𝑐 𝑡)𝑒 (𝑡-𝜇) 2 /2𝜎 2	𝑇𝑠: 139 µs (5 periods)	𝑇𝑠: 83 µs (5 periods)
						𝐹𝑐𝑜: [26.46 45.54] kHz	𝐹𝑐𝑜: [44.09 75.90] kHz
	Half-Gaussian	ℎ𝑟(𝑡) = {	0 if 𝑡 < 𝜇 𝑟(𝑡) if 𝑡 > 𝜇	𝐹𝑐: 36 kHz	𝐹𝑐: 60 kHz
	pulse					𝑇𝑠: 69 µs (2.5 periods)	𝑇𝑠: 42 µs (2.5 periods)
						𝐹𝑐𝑜: [24.20 55.28] kHz	𝐹𝑐𝑜: [24.50 81.38]
	Chirp	𝑐(𝑡) = 𝑐𝑜𝑠(2𝜋𝑓(𝑡)𝑡 +	𝐹𝑐: 36 kHz	𝐹𝑐: 60 kHz
		𝜙 0 )𝑒 (𝑡-𝜇) 2 /2𝜎 2		Δ𝐹: 28kHz	Δ𝐹: 48kHz
		𝑓(𝑡) = (𝑓 1 -𝑓 0 )	𝑡 𝑇𝑠 + 𝑓 0	𝑓 0 : 22kHz, 𝑓 1 : 50 kHz	𝑓 0 : 36kHz, 𝑓 1 : 84 kHz
						𝑇𝑠: 45 µs (10 periods)	𝑇𝑠: 27 µs (10 periods)
						𝐹𝑐𝑜: [32.57 40.04] kHz	𝐹𝑐𝑜: [54.55 67.22] kHz

Table 6 :

 6 Mean and standard deviation of RMS values for received signals, sorted from higher to lower.

	Sensor	Signal	𝑹𝑴𝑺 ̅̅̅̅̅̅̅ (mV)	𝝈(𝑹𝑴𝑺) (mV)
	3alpha	Train	90.6	74.4
		Half Gaussian pulse	56.7	37.0
		Impulse	54.8	37.8
		Chirp	30.5	17.1
		Gaussian pulse	29.4	11.1
	6alpha	Half Gaussian pulse	161.3	106.8

Table 7 :

 7

	Sensor	Signal	𝒔(𝒕)	𝒚(𝒕)	Ratio
	3alpha	Impulse	50.0	54.8	-79.2
		Half Gaussian pulse	54.9	56.7	-79.7
		Train	139.6	90.6	-83.7
		Gaussian pulse	78.5	29.4	-88.6
		Chirp	141.9	30.5	-93.6
	6alpha	Half Gaussian pulse	45.0	161.3	-68.9
		Gaussian pulse	60.8	116.6	-74.3
		Impulse	50.0	86.3	-75.2
		Train	109.5	47.9	-87.1
		Chirp	109.4	40.1	-88.6
	Figure 22 presents the signal-to-noise ratio (SNR) mean and standard deviation
	values. Table 8 summarizes the SNR values, sorting by SNR mean value in descending
	order. Average SNR values for all receiver angles ranged between 20 and 40 dB, indicating
	a low presence of noise. The only exception corresponded to chirp signal when using the

Ratio [dB]

between output (y(t)

(VRMS, in [mV]

) before 40 dB amplification) and input (s(t)

(VRMS, in [mV]

) after 40 dB amplification) RMS values for all signals, sorted from higher to lower.

R6α located at 45°, where mean SNR was around 10 dB. As obtained for the RMS measurements, SNR values for the sensor R6α were higher than those obtained for R3α.

Impulsive-like signals, as the pulse train and impulse, presented the highest SNR ratios.

Table 8 :

 8 Mean and standard deviation of SNR values for received signals, sorted from higher to lower.

	Sensor	Signal	𝑺𝑵𝑹 ̅̅̅̅̅̅ (dB)	𝝈(𝑺𝑵𝑹) (dB)
	3alpha	Train	33.11	6.48
		Half Gaussian pulse	29.00	7.09
		Impulse	32.67	5.31
		Chirp	27.58	7.21
		Gaussian pulse	29.77	5.08
	6alpha	Half Gaussian pulse	30.51	6.90
		Gaussian pulse	35.02	11.75
		Impulse	40.52	12.35
		Train	47.9	41.6
		Chirp	40.1	22.4

Table 9 :

 9 Absolute mean differences between TOF obtained with R3α and R6α sensors.

	Angle [°]	Δ Impulse [µs]	Δ Chirp [µs]
	45	1.45	9.08
	90	0.47	25.9
	135	0.54	10.4
	180	0.29	2.50

  . Compute the spatial displacements.

b. Delete crossing points. c. Smooth curve. d. Update wavefront. e. Evaluate if a final front VI.

Table 10 :

 10 Elastic constants for 6 wood species and 1 isotropic material.

	Species/Material Type	Density	𝑬 𝑹	𝑬 𝑻	𝑮 𝑹𝑻	𝝂 𝑹𝑻 𝑬 𝑹 /𝑬 𝑻
			(kg/m 3 )	(MPa)	(MPa)	(MPa)
	Spruce, Sitka	Softwood	448	927	511	36	0.44 1.81
	Douglas-fir	Softwood	538	909	668	94	0.39 1.36
	Pine, Longleaf	Softwood	661	1537	829	181	0.38 1.85
	Yellow-poplar	Hardwood	470	1103	516	132	0.70 2.14
	Sweetgum	Hardwood	582	1429	622	261	0.68 2.30
	Oak, red	Hardwood	706	2118	1128	319	0.56 1.88
	PVC	Isotropic	1550	2100	2100	800	0.31 1

Table 11 :

 11 Elastic constants used for the ray-tracing simulations. Wood parameters at 12% moisture content.

	Type 𝝆 (kg/m3) 𝑬 𝑹 (MPa) 𝑬 𝑻 (MPa) 𝑮 𝑹𝑻 (MPa) 𝝂 𝑹𝑻
	Oak	570	1180	614	319	0.56
	Ash	680	1540	820	280	0.68
	PVC 1550	2100	2100	800	0.31

1.2. Meeting with experts

ACKNOWLEDGMENTS

Concerning the angle, variations were larger when the sensor position angle was 45°, and decreased as this angle approached 180°, this is the sensor located opposed at radial direction. Lower variability was obtained for impulse signal, with coefficients of variation ranging from 0.33% to 1.67%. The Gaussian signal presented the larger variability, considering that the AIC and threshold method works better with an initial impulsive signal.

Figure 26: Relative standard deviation for TOF values using the Threshold method (up) and the AIC method (down) for R3α [START_REF] Espinosa | Accuracy on the Time-of-Flight Estimation for Ultrasonic Waves Applied to Non-Destructive Evaluation of Standing Trees: A Comparative Experimental Study[END_REF]).

In the case of the sensor R6α, Figure 27 presents the mean and standard deviation values for the TOF estimated using the threshold technique and AIC methods, for all signals except chirp. Mean TOF values ranged between 66 μs to 143 μs, equivalent to the values for sensor R3α. The difference in mean values estimated with the two methods was always inferior to 1 μs. Standard deviation ranged between 0.5 μs to 5.2 μs for threshold method and 0.5 μs to 5.7 μs for the AIC method, slightly lower than the difference for sensor R3α.

Figure 28 shows the relative standard deviation presented as a percentage. For the receiver position angle, variations were larger when the sensor position angle was 45° and decreased as this angle approached 180°. Coefficients of variation obtained for half-Gaussian, impulse and pulse train signals were similar, always inferior to 3% for both AIC and threshold approaches. The Gaussian signal presented the larger variability again, reaching 7% when the sensor was located at 180°. tree and three TOF detection methods. Among all configurations, the one presenting less variation on the TOF measurements was the combination of an encoded excitation signal, such as chirp signal, with cross-correlation to measure the time delay. This combination was used for the experimental testing in the following chapters. Chirp signals deserve attention considering that this signal was adjusted to the transducer response and the received signals concentrated energy in frequency bands around the resonant frequency of sensors. The sensor position affected the consistency on time measurements: as the sensor position angle approached to the radial direction, the TOF values presented less variation.

APPENDIX A: TECHNIQUES FOR ASSESSING HAZARDOUS TREES

Visual Tree Assessment

This is a method for the visual inspection of standing trees, composed of three steps [START_REF] Mattheck | Field Guide for Visual Tree Assessment (vta)[END_REF]):

1. Visual inspection aiming to detect symptoms of defects. The symptoms are associated to several factors, such as the tree vitality (dry branches, loose bark, poor growth), the presence of fungal fruiting bodies and open wounds, bulges (symptom for decay), ribs (symptom for radial cracks), 'crackledvarnish' zones in the bark (symptom for brittle fracture). If the tree presents a significant hazard, it follows step 2.

2. Confirmation of defects and determination of residual wall thickness. To perform this confirmation, inspection could include drilling resistance measurement, sound velocity measurement, measurement of the wood strength (Fractometer), annual ring analysis.

3. Assessment of the defects considered to be critical. This include criteria of failure for hollow or decayed trees and criteria of failure for root damages.

For hollow, a criterion using the intact residual wall thickness t and the tree radius R is defined as t/R>0.3.

Decay detection tools

Decay detection assessment can be performed using a variety of devices [START_REF] Pokorny | Urban tree risk management: A community guide to program design and implementation[END_REF]. This include basic devices such as a rubber mallet or an increment borer. More complex devices include penetrometers, sonic and ultrasonic detectors, electrical conductivity meter and fractometer.

1. Rubber mallet: With experience, an operator can hit the bark with this hammer and interpret the resulting sound to determine the presence of sever decay (holes). This method is non-invasive and cheap.

2. Increment borer: this tool allows extracting a core of wood, of approximately 5 mm in diameter, that can be examined for the presence of

APPENDIX B: COMPARTMENTALIZATION OF DECAY IN TREES

When the bark is wounded, the tree will react to stop the spread of pathogens in the healthy wood. Trees respond by compartmentalizing, they attempt to wall off the injured or infected region. To illustrate tree response to wounding and invasion by decay organisms, a model was developed known as CODIT [START_REF] Shigo | Ultrasonic goniometry immersion techniques for the measurement of elastic moduli[END_REF].

Figure A1: CODIT system: 4 different walls to stop the spread of infected wood [START_REF] Shigo | Ultrasonic goniometry immersion techniques for the measurement of elastic moduli[END_REF].

Three walls appear at the time of the injury. Wall 1 resists the vertical spread of infected wood, wall 2 the inward spread and wall 3 the lateral spread. The strongest wall is formed after the wounding and is formed by the cambium, known as wall 4.

Top and bottoms walls 1 are the weakest walls. After wounds, this wall is formed in the vertical vascular system by gums, resins and tyloses. Wall 2 and wall 3 exist before the infection, as they are formed by the annual rings and by the parenchymal rays respectively. After the wounding, these walls are chemically strengthened. Wall 4 is called a barrier zone. It separates wood formed before wounding from wood formed after wounding. Wall 4 is the strongest wall, capable of confining decay to tissues formed prior to wounding. Zhao J, Qiu J, Ji H (2016) Reconstruction of the nine stiffness coefficients of composites using a laser generation based imaging method. Composites Science and Technology 126:27-34. doi: 10.1016/j.compscitech.2016.02.001 Zhu D, Conners R, Araman P (1991) CT image sequence processing for wood defect recognition. IEEE Comput. Soc. Press, pp 75-79