Kailiang Xu

Wang

Pascal Laugier

Dispersion Curves Extraction of Guided Waves 1 Multichannel Processing for Dispersion Curves Extraction of Ultrasonic Axial-Transmission Signals: Comparisons and Case Studies

Keywords: Waveguides, Singular value decomposition, Linear Radon transform, Dispersion curves PACS: 43.60.Fg, 43.60.Gk, 43.80.Vj, 43.80.Jz, 43.80.Ev

Some pioneering studies have shown the clinical feasibility of long bones evaluation using ultrasonic guided waves. Such a strategy is typically designed to determine the dispersion information of the guided modes to infer the elastic and structural characteristics of cortical bone. However, there are still some challenges to extract multimode dispersion curves due to many practical limitations, e.g., high spectral density of modes, limited spectral resolution and poor signal-to-noise ratio (SNR). Recently, two representative signal processing methods have been proposed to improve the dispersion curves extraction.

The first method is based on singular value decomposition (SVD) with advantages of multi-emitter and multi-receiver configuration for enhanced mode extraction; the second one uses linear Radon transform (LRT) with high-resolution imaging of the dispersion curves. To clarify the pros and cons, a face to face comparison was performed between the two methods. The results suggest that the LRT method is suitable to separate the guided modes at low frequency-thickness-product (𝒇 • 𝒉) range; for multimode signals in broadband 𝒇 • 𝒉 range, the SVD-based method shows more robust performances for weak mode enhancement and noise filtering. Different methods are valuable to cover the entire 𝒇 • 𝒉 range for processing ultrasonic axial transmission signals measured in long cortical bones.

I. INTRODUCTION

In the past decade, significant progress has been achieved in quantitative ultrasound (QUS) assessment of cortical bone using axial transmission techniques [START_REF] Talmant | Guided Waves in Cortical Bone[END_REF][START_REF] Moilanen | Ultrasonic guided waves in bone[END_REF] . Such techniques are intended to osteoporotic fracture discrimination [START_REF] Barkmann | A new method for quantitative ultrasound measurements at multiple skeletal sites -First Results of precision and fracture discrimination D-1837-2010 C-9752-2010[END_REF] or to fracture healing monitoring [START_REF] Protopappas | Ultrasonic monitoring of bone fracture healing[END_REF][START_REF] Xu | Axial transmission method for long bone fracture evaluation by ultrasonic guided waves: Simulation, phantom and in vitro experiments[END_REF] . Several axial transmission techniques have been explored using either one transmitter/one receiver [START_REF] Foldes | Quantitative ultrasound of the tibia: a novel approach for assessment of bone status[END_REF] , one transmitter/multiple-receiver [START_REF] Nicholson | Guided ultrasonic waves in long bones: Modelling, experiment and in vivo application[END_REF] , or multiple-transmitter/multiple-receiver configurations [START_REF] Minonzio | Guided wave phase velocity measurement using multi-emitter and multi-receiver arrays in the axial transmission configuration[END_REF] . Signal analysis developed for clinical applications are based on velocity measurement of a single signal, such as the first arriving signal (FAS) [START_REF] Foldes | Quantitative ultrasound of the tibia: a novel approach for assessment of bone status[END_REF] , or a slower second energetic wave which has been interpreted as the fundamental anti-symmetrical A0 Lamb mode based on plate model [START_REF] Nicholson | Guided ultrasonic waves in long bones: Modelling, experiment and in vivo application[END_REF] . If the cortical thickness is much smaller than the longitudinal wavelength, FAS can be seen as the S0 Lamb mode. If the cortical thickness is much larger than the longitudinal wavelength, the FAS corresponds to the non-dispersive lateral compression wave [START_REF] Nicholson | Guided ultrasonic waves in long bones: Modelling, experiment and in vivo application[END_REF] . In between, when the thickness is comparable to the wavelength (e.g., at 1 MHz, the wavelength in cortical bone is 3 to 4 mm), a physical interpretation of the FAS has yet to be provided. Besides FAS, reflection and conversion of ultrasonic body waves have been observed in the relatively thick bovine tibia ex vivo [START_REF] Le | Probing long bones with ultrasonic body waves[END_REF] .

Whereas the signal analysis techniques applied so far in axial transmission meet the need for simplicity and pragmatism, the corresponding biomarkers extracted from a single signal (either FAS or A0 mode) provide an incomplete biomechanical characterization of bone strength. Such a consideration has motivated the research aiming at fully characterizing the response of cortical long bones to an ultrasonic excitation. Particularly, considering that long bones actually support the propagation of multiple guided modes, several studies have shown interest in measuring and interpreting the dispersion curves of guided waves within a wideband frequency-wavenumber range [START_REF] Talmant | Guided Waves in Cortical Bone[END_REF][START_REF] Moilanen | Ultrasonic guided waves in bone[END_REF] . Measurement of dispersion curves, along with suitable waveguide modeling, has been proposed for concurrent estimation of cortical thickness and elastic properties [START_REF] Vallet | In vivo characterization of cortical bone using guided waves measured by axial transmission[END_REF][START_REF] Kilappa | Assessment of the fundamental flexural guided wave in cortical bone by an ultrasonic axial-transmission array transducer[END_REF][START_REF] Tatarinov | Multi-frequency axial transmission bone ultrasonometer[END_REF][START_REF] Song | Analysis of superimposed ultrasonic guided waves in long bones by the joint approximate diagonalization of eigen-matrices algorithm[END_REF][START_REF] Foiret | Combined estimation of thickness and velocities using ultrasound guided waves: a pioneering study on in vitro cortical bone samples[END_REF] . However, the extraction and accurate interpretation of the dispersion characteristics of guided modes propagating along the cortical shell of long bones pose difficulties, such as the high spectral density of modes, limited spectral resolution and poor signal-to-noise ratio (SNR).

Particularly, because of the limited spatial sampling or resolution achievable in clinical measurements with clinical probes (typically a few cm-long array), the traditional two-dimensional (or spatio-temporal) Fourier transform (2D-FT) [START_REF] Chapman | A new method for computing synthetic seismograms[END_REF] cannot achieve a high wavenumber resolution for complete multimode separation [START_REF] Moilanen | Ultrasonic guided waves in bone[END_REF] . From guided signal processing point of view, several methods have been described to distinguish overlapped modes and to extract the corresponding dispersion curves. These include dispersion-based short-time Fourier transform [START_REF] Hong | Dispersion-based short-time Fourier transform applied to dispersive wave analysis[END_REF] , group velocity filtering [START_REF] Moilanen | Assessment of the cortical bone thickness using ultrasonic guided waves: Modelling and in vitro study[END_REF] , time-frequency representation 18 , warped frequency transform [START_REF] Marchi | Ultrasonic guided-waves characterization with warped frequency transforms[END_REF] , time-frequency ridge extraction [START_REF] Xu | Multiridge-based analysis for separating individual modes from multimodal guided wave signals in long bones[END_REF][START_REF] Zhang | Joint spectrogram segmentation and ridge-extraction method for separating multimodal guided waves in long bones[END_REF] , dispersion compensation [START_REF] Xu | Mode separation of Lamb waves based on dispersion compensation method[END_REF] , blind identification [START_REF] Song | Analysis of superimposed ultrasonic guided waves in long bones by the joint approximate diagonalization of eigen-matrices algorithm[END_REF] , generalized warblet transform [START_REF] Yang | Dispersion analysis for broadband guided wave using generalized warblet transform[END_REF] , adaptive Chirplet transform [START_REF] Zeng | Waveform separation and image fusion for Lamb waves inspection resolution improvement[END_REF] , orthogonality relation-based method [START_REF] Ratassepp | Application of orthogonality-relation for the separation of Lamb modes at a plate edge: Numerical and experimental predictions[END_REF] and compressed sensing method [START_REF] Harley | Sparse recovery of the multimodal and dispersive characteristics of Lamb waves[END_REF] . Recently, two methods taking advantage of a multiple-transmitter/multiple-receiver implementation of axial transmission have been proposed, singular vector decomposition (SVD) [START_REF] Minonzio | Guided wave phase velocity measurement using multi-emitter and multi-receiver arrays in the axial transmission configuration[END_REF][START_REF] Minonzio | Impact of attenuation on guided mode wavenumber measurement in axial transmission on bone mimicking plates[END_REF] and linear Radon transform (LRT) [START_REF] Tran | Imaging ultrasonic dispersive guided wave energy in long bones using linear Radon transform[END_REF][START_REF] Tran | Multichannel filtering and reconstruction of ultrasonic guided wave fields using time intercept-slowness transform[END_REF] .

The SVD-based approach is able to significantly enhance the weak-amplitude modes by selecting the singular values and singular vectors corresponding to the signals and by filtering out the small singular values corresponding to the noise [START_REF] Minonzio | Guided wave phase velocity measurement using multi-emitter and multi-receiver arrays in the axial transmission configuration[END_REF] . This method, applied in vitro on radius specimens [START_REF] Foiret | Combined estimation of thickness and velocities using ultrasound guided waves: a pioneering study on in vitro cortical bone samples[END_REF] and in vivo at the forearm [START_REF] Vallet | In vivo characterization of cortical bone using guided waves measured by axial transmission[END_REF] , has shown a good performance to estimate cortical thickness. The LRT method is widely known in seismic data processing [START_REF] Chapman | Generalized Radon transforms and slant stacks[END_REF][START_REF] Sacchi | Statistical and transform methods in geophysical signal processing[END_REF] . Taking advantage of the sparse inversion, the LRT methods are capable of achieving an energy focusing or so-called high-resolution dispersive energy imaging of the wave-packets whose arrival times are linearly dependent on the propagation path [START_REF] Thorson | Velocity-stack and slant-stack stochastic inversion[END_REF][START_REF] Sacchi | High-resolution velocity gathers and offset space reconstruction[END_REF] . The LRT methods have already been applied to study the surface-wave data for Rayleigh wave dispersive energy imaging [START_REF] Luo | Rayleigh-wave dispersive energy imaging using a high-resolution linear Radon transform[END_REF] . Recently, the high-resolution LRT method has been introduced to analyze the ultrasonic guided signals in long bone, and was observed to provide enhanced resolution of the extracted dispersion curves in comparison to the 2D-FT method [START_REF] Tran | Imaging ultrasonic dispersive guided wave energy in long bones using linear Radon transform[END_REF][START_REF] Tran | Multichannel filtering and reconstruction of ultrasonic guided wave fields using time intercept-slowness transform[END_REF] .

While results on bone have been reported with SVD-based and LRT approach by different groups, the pros and cons as well as the applicability conditions of the two different signal processing techniques have not been thoroughly discussed. In this study, a face to face comparison is performed using SVD-based, LRT and the classical 2D-FT method. Synthetic and experimental signals on bone-mimicking phantoms and on ex-vivo human radius are analyzed.

II. BRIEF REVIEW OF THE PROPOSED SIGNAL PROCESSING METHODS

Typically, the frequencies used for the axial transmission measurement of long cortical bone are in the range between 50 kHz and 2 MHz and the typical cortical thickness of human radius and tibia are in the range between 1 to 6 mm approximately. According to the characteristics of the guided signals observed in different ranges of the frequency-thickness product (𝑓 • ℎ), three representatives cases can be distinguished, as shown in

III. THEORY AND METHODS

A. Guided Waves Dispersion

Lamb modes in plates are classified as symmetric (S) and antisymmetric (A) modes, briefly named as 𝑆𝑛 and 𝐴𝑛 modes (𝑛 = 0,1, 2, … ), which are solutions of the Rayleigh-Lamb equations [START_REF] Viktorov | Rayleigh and Lamb waves: physical theory and applications[END_REF][START_REF] Su | Guided Lamb waves for identification of damage in composite structures: a review[END_REF] . The dispersion curves can be expressed as wavenumber 𝑘 versus 𝑓 or 𝑓 • ℎ.

The 2D-FT provides a general relationship between the distance-time space (𝑥, 𝑡) and frequency-wavenumber space (𝑓, 𝑘). 𝑆(𝑘, 𝑓) = ∬ 𝑔(𝑥, 𝑡)𝑒 𝑗(𝑘𝑥-2𝜋𝑓𝑡) 𝑑𝑥𝑑𝑡.

(1),

+∞ -∞

where 𝑔(𝑥, 𝑡 ) is the temporal signals recorded at a series of space positions x. The dispersion curves can thus be obtained by locating the wavenumbers at each frequency where the amplitude of 2D-FT result 𝑆(𝑘, 𝑓) reaches the maxima.

With a given dispersion curve and an excitation signal, spectrum of the dispersive signal at distance 𝑥 0 can be computed by multiplying a phase-spectrum adjustment term 𝑒 -𝑗𝑘(𝑓) 𝑥 0 to the spectrum of excitation.

The temporal waveforms can thus be obtained by using inverse Fourier transform to the phase-adjusted spectrum of excitation. Such a procedure provides us an efficient way to synthesize the temporal signal 𝑔(𝑥, 𝑡) [START_REF] Xu | Mode separation of Lamb waves based on dispersion compensation method[END_REF] .

B. Singular Value Decomposition-based Wavenumber Extraction

The are retained as a basis of the received signal subspace. The wavenumber determination is then achieved by projecting a testing vector onto this basis.

The testing vector is expressed as an attenuated spatial plane wave with a complex wavenumber 𝑘 + 𝑗𝛼

27 𝑣 𝑡𝑒𝑠𝑡 (𝑘, 𝛼) = 1 √ 𝑁 𝑟 2𝛼𝐿 (1 -𝑒 -2𝛼𝐿 ) 𝑒 -𝑗𝑘𝑥-𝛼𝑥 (2a), Denominator term √ 𝑁 𝑟 2𝛼𝐿 (1 -𝑒 -2𝛼𝐿
) is used to normalize 𝑣 𝑡𝑒𝑠𝑡 (𝑘, 𝛼) [START_REF] Minonzio | Impact of attenuation on guided mode wavenumber measurement in axial transmission on bone mimicking plates[END_REF] . 𝛼 is the attenuation coefficient and in the study, a constant was adopted, i.e., 𝛼 = 0.05 𝑁𝑝 • 𝑚𝑚 -1 . Such a value is close to the average attenuation coefficients of the low-order guided modes in the bone mimicking materials [START_REF] Minonzio | Impact of attenuation on guided mode wavenumber measurement in axial transmission on bone mimicking plates[END_REF] . With a constant attenuation coefficient, 𝑣 𝑡𝑒𝑠𝑡 (𝑘, 𝛼) is a function of 𝑘.

At each frequency (𝑘 0 , 𝑓 0 ), the projection of 𝑣 𝑡𝑒𝑠𝑡 onto the 𝑁 𝑟𝑘 first reception singular vector basis 𝑣 𝑖 (𝑖 = 1, … , 𝑁 𝑟𝑘 ) yields the so-called Norm function 𝑆 SVD (𝑘, 𝑓) [START_REF] Minonzio | Guided wave phase velocity measurement using multi-emitter and multi-receiver arrays in the axial transmission configuration[END_REF] 𝑆 𝑆𝑉𝐷 (𝑘 0 , 𝑓

0 ) = ∑ |⟨𝑣 𝑡𝑒𝑠𝑡 (𝑘 0 , 𝛼)| 𝑣 𝑖 (𝑓 0 )⟩| 2 𝑁 𝑟𝑘 (𝑓 0 ) 𝑖=1 (2b),
where ⟨𝑣 𝑡𝑒𝑠𝑡 (𝑘 0 , 𝛼)| 𝑣 𝑖 (𝑓 0 )⟩ is the Hermitian scalar product between the 𝑣 𝑡𝑒𝑠𝑡 (𝑘 0 , 𝛼) and 𝑣 𝑖 (𝑓 0 ). Thus, at each frequency 𝑓 0 , the Norm function is a function of 𝑘 and the maxima correspond to wavenumbers of the guided waves presented in the received signals.

Due to the normalized characteristics of the orthogonal basis, the values of Norm function range from 0 to 1. This value can be interpreted as follows: if a guided mode exists in the signal, the corresponding Norm function value is close to 1; otherwise, the value is close to 0 8 . Details of the SVD-based method and corresponding examples can be learned from [START_REF] Minonzio | Guided wave phase velocity measurement using multi-emitter and multi-receiver arrays in the axial transmission configuration[END_REF][START_REF] Xu | Sparse SVD Method for High Resolution Extraction of the Dispersion Curves of Ultrasonic Guided Waves[END_REF] . Note that the method can be adapted to achieve a more accurate dispersion extraction when the cortical thickness is not uniform but with linear changes in the direction of wave propagation [START_REF] Moreau | Measuring the wavenumber of guided modes in waveguides with linearly varying thickness[END_REF] .

C. Linear Radon Transformation based Wavenumber Extraction

For convenience of the inversion problem, the LRT is usually formulated as the forward relationship between the (𝑥, 𝑡) and the Radon domain. Let 𝑔(𝑥, 𝑡) represent a distance-time matrix measured by single emitter and 𝑁 𝑟 receivers. The LRT can be defined using the following equation [START_REF] Chapman | Generalized Radon transforms and slant stacks[END_REF][START_REF] Sacchi | Statistical and transform methods in geophysical signal processing[END_REF] 𝑔(𝑥, 𝑡) = ∫ 𝑤(𝑝, 𝜏 = 𝑡 -𝑝𝑥)𝑑𝑝 𝑝 (3),

where 𝜏 and 𝑝 denote the intercept time and slope parameter or phase slowness, respectively; 𝑤(𝑝, 𝜏 = 𝑡-𝑝𝑥) designates the signal in the (𝜏, 𝑝) domain. The LRT is commonly named as 𝜏-𝑝 transform or slant stack.

Equation ( 3) can be rewritten in frequency domain as [START_REF] Tran | Imaging ultrasonic dispersive guided wave energy in long bones using linear Radon transform[END_REF][START_REF] Sacchi | Statistical and transform methods in geophysical signal processing[END_REF] 𝐺(𝑥, 𝑓) = ∫ 𝑊(𝑝, 𝑓)𝑒 -𝑗2𝜋𝑓𝑝𝑥 𝑑𝑝 𝑝 (4).

Using matrix notation, we have

𝐺 = 𝐿𝑊 (5),
The 𝑁 𝑟 × 𝑁 𝑝 linear operator 𝐿 is

𝐿 = [ 𝑒 -𝑗2𝜋𝑓𝑝 1 𝑥 1 … 𝑒 -𝑗2𝜋𝑓𝑝 𝑁 𝑝 𝑥 1 ⋮ ⋱ ⋮ 𝑒 -𝑗2𝜋𝑓𝑝 1 𝑥 𝑁 𝑟 … 𝑒 -𝑗2𝜋𝑓𝑝 𝑁 𝑝 𝑥 𝑁 𝑟 ] 𝑝 = 𝑝 1 , 𝑝 2 , … , 𝑝 𝑁 𝑝 (6),
where 𝑁 𝑝 is the number of the slowness sampling points.

A penalized least-squares (LS) solution to Eq. ( 5) has been introduced by minimizing the following cost function [START_REF] Sacchi | Statistical and transform methods in geophysical signal processing[END_REF][START_REF] Thorson | Velocity-stack and slant-stack stochastic inversion[END_REF] 𝐽 = ‖𝐺 -𝐿𝑊‖ 2 2 + 𝜇𝑄(𝑊) (7),

where the Lagrange multiplier 𝜇, also referred as regularization hyperpameter [START_REF] Schultz | Flexible, inversion-based Matlab implementation of the Radon transform[END_REF] , determines the trade-off between the misfit term ‖𝐺 -𝐿𝑊‖ 2 2 describing the data fidelity and the penalty term 𝑄(𝑊) in Radon domain.

A typical implementation of the penalty term is to use the quadratic 𝑙 2 -norm, i.e., 𝑄(𝑊) = ‖𝑊‖ 2 2 . So that, the LS solution can be analytically obtained as

𝑊 ̃= (𝐿 𝐻 𝐿 + 𝜇𝐼) -1 𝐿 𝐻 𝐺 (8).
where 𝐿 𝐻 is the complex-conjugate transpose of 𝐿. It should be noticed that 𝐿 𝐻 𝐿 is not a unitary matrix.

The 𝑊(𝑝, 𝑓) in the Radon field can be readily mapped to the (𝑘, 𝑓) domain via 𝑘 = 2𝜋𝑓𝑝 [START_REF] Tran | Imaging ultrasonic dispersive guided wave energy in long bones using linear Radon transform[END_REF] . The dispersion curves can also be obtained from |𝑆 𝐿𝑅𝑇 (𝑘, 𝑓) | 2 by using LRT method.

If certain non-quadratic terms, which enable to quantify the amount of sparsity of a vector, are adopted as the penalty terms, then we can focus the signal energy on the "best" subspace of the solution spaces. That is actually the sparse solution leading to the high-resolution LRT method. Therefore, the penalty term actually controls the high resolution constraints and also indicates the sparsity of the results. Two typical non-quadratic penalty terms, i.e., 𝑙 1 -norm 𝑄(𝑊) = ‖𝑊‖ 1 and Cauchy norm, are usually adopted for achieving sparsity. The Cauchy norm penalty term is defined as 33

𝑄(𝑊) = ∑ 𝑙𝑛[1 + 𝑊(𝑝 𝑖 , 𝑓)/𝜀 2 ] (9). 𝑁 𝑝 𝑖=1
where the scale factor of the Cauchy distribution 𝜀 2 actually indicates the default power in absence of hyperbolic events [START_REF] Sacchi | High-resolution velocity gathers and offset space reconstruction[END_REF] . According to the discrepancy principle, a proper 𝜀 2 value should be selected to ensure that the misfit matches the power of noise [START_REF] Sacchi | Reweighting strategies in seismic deconvolution[END_REF] . The one-dimension Brent parabolic interpolation method has been used to compute the epsilon [START_REF] Sacchi | Reweighting strategies in seismic deconvolution[END_REF][START_REF] Press | Numerical recipes 2nd edition: The art of scientific computing[END_REF] . It is fixed as a constant 𝜀 = 0.8 in [START_REF] Tran | Multichannel filtering and reconstruction of ultrasonic guided wave fields using time intercept-slowness transform[END_REF] . In the study, 𝜀 = 1 was used for the Cauchy LRT computation. Although there is no analytical solution to high-resolution LRT, it can be solved efficiently by conjugate gradient technique. Details of that can be found in literature [START_REF] Schultz | Flexible, inversion-based Matlab implementation of the Radon transform[END_REF][START_REF] Sacchi | Reweighting strategies in seismic deconvolution[END_REF] .

The trade-off curve, so-called L-curve [START_REF] Engl | Using the L-curve for determining optimal regularization parameters[END_REF] , is usually applied for empirical optimization of the 𝜇 value for the penalized LS-LRT and high-resolution LRT method. A small hyperparameter 𝜇 value leads to a solution with minimized misfit term but with less energy focusing; conversely, a larger 𝜇 value can achieve a high resolution by emphasizing the regularization term. For simplicity and methodology comparison, a fixed value of trade-off 𝜇 = 0.05 were used for the results presented in Section IV. Detailed discussion of trade-off parameter using L-curve can be found in Section V(1).

D. Experiments

Experiments are achieved using an array transducer consisting in 5 emitters and 24 receivers (Vermon, Tours, France) associated with a specific driving electronic device (Althaï s, Tours, France). The pitch of the array transducer is 0.8 mm and the length and width of each rectangular element are 0.8 mm and 8 mm, respectively. The central frequency is 1 MHz and the -6 dB bandwidth goes from 0.5 to 1.6 MHz 10 .

For case II in Table I

IV. RESULTS

A. Synthetic signals, narrow wavenumber-band Lamb modes S0 and A0

This example corresponds to case I in Table I. Two fundamental narrow k-band (0 < 𝑘 ≤ 2 𝑟𝑎𝑑 • 𝑚𝑚 -1 ) Lamb modes A0 and S0 on a 2 mm-thick bone-mimicking plate were synthesized with peak-to-peak amplitudes of 1 and 0.3, respectively (see Fig. 1a). A Gaussian random noise was added with SNR of 30dB. The SNR is defined as the ratio of the power of the signal and that of the noise. Fig. 1b presents the 2D-FT results in the (𝑘, 𝑓) domain. After SVD decomposition, the singular values were normalized in dB scale. Those singular vectors associated with singular values above the threshold 21dB

were remained as the signal subspaces and the rest were filtered out as noise (see Fig. 7a ). The SVD result is depicted in Fig. 1c. Fig. 1d is the (𝜏, 𝑝) result obtained by LRT with 𝑙 2 -norm. Since there is a large velocity difference between the S0 and A0 modes (see Fig. 1a), they are projected as two separate regions in the (𝜏, 𝑝) domain (Fig. 1d). As shown in Figs. 1e-f, the high-resolution LRT results using 𝑙 1 -norm and

Cauchy norm are able to significantly concentrate the 𝑘-𝑓 energy of the narrow wavenumber-band and S0 and A0 modes. The colors of the 𝑘-𝑓 and 𝜏-𝑝 energy distribution present the mode energy with highest values in red and lowest values in blue.

B. Phantom signals (1) Wide wavenumber-band and multiple guided modes

This example corresponds to case II in Table I. Compared with the 2D-FT method, using the SVD-based method, the multimode dispersion curves can be identified with an enhancement of the weak modes, e.g., A0, A1, S0, S4 etc., and also the low-amplitude S0 and A1 mode energy in wavenumber range of 3 < 𝑘 < 4 𝑟𝑎𝑑 • 𝑚𝑚 -1 . In Fig. 2d, the LRT method can obtain a projection in the (𝜏, 𝑝) domain with focused energy points. But as shown in Figs. 2e-f, the concentrated region in the (𝜏, 𝑝) domain only represents the strong modal energy close to the center frequency, which actually cannot lead to an effective reconstruction of the dispersion trajectories of the weak modes. For example, the A0, A1 and S4 modes, whose energy is far from the center frequency of the probe, are not clearly depicted on the LRT results.

(

2) Body waves

This example corresponds to case III in Table I Furthermore, such a slant-stack operator is very suitable to detect the weak components (see Fig. 3d).

However, as shown in Figs. 3e-f, for the signals measured in high 𝑓 • ℎ range, there is still no evidence that the 𝑙 1 -norm and Cauchy norm LRT method can improve the k-f resolution for better dispersive energy imaging. The SVD-based method still provides a result with best dispersion energy extraction in 𝑘-𝑓 domain.

C. Ex vivo guided signals analysis in a human radius

This example corresponds to case II in Table I. The guided signals measured from an ex-vivo 2.5 mm-thick human radius can be seen in Fig. 4a. Figs. 4b-c are the 2D-FT and SVD 𝑘-𝑓 results, respectively. Fig. 4d is the multimode energy distributions in the (𝜏, 𝑝) domain using the 𝑙 2 -norm LRT method. The experimental SNR is around 55dB (see Fig. 7b). The threshold of singular values is 20dB. The resolution improvement of the fundamental A1, S1 and S2 modes can be observed close to center frequency bandwidth in Fig. 4e using the 𝑙 1 -norm LRT method, but for the relatively weak and wideband modes, e.g., A0, A2 and S0, the LRT cannot provide sufficient mode enhancement. As shown in Fig. 4f, due to the improper value of the regularization parameter, it seems that the Cauchy norm LRT method enforces the results with "all-zero" solution. The reason of that will be discussed in Section V.

V. DISCUSSION

In this study, we performed a face to face comparison between two signal processing approaches, namely the SVD-based and LRT method, which have been recently proposed to extract the dispersion curves of guided waves transmitted in long bone. To this goal, the methods were applied to synthetic signals and experimental signals recorded on a bone-mimicking plate and on a human radius ex vivo.

A. Parameter optimization of the LRT and SVD-based method

The hyperparameter 𝜇 of the LRT methods, which controls the trade-off between data fidelity and mode energy concentration (or sparseness), can be heuristically determined by the L-curve. The optimal value of hyperparameter, usually determined on the "elbow" of L-curves, actually corresponds to the maximal curvature point where the misfit and penalty terms are minimized together [START_REF] Tran | Imaging ultrasonic dispersive guided wave energy in long bones using linear Radon transform[END_REF][START_REF] Engl | Using the L-curve for determining optimal regularization parameters[END_REF] . For the SVD-based method, a SNR threshold is used to selectively separate the noise and signal spaces.

(

1) Hyperparameter μ of the LRT method

It has been observed in Section IV that for case I in Table 1, the LRT methods provide similar 𝑘-𝑓 dispersion loci as 2D-FT method; but for cases II and III, the over-sparse solutions of k-f energy distributions are readily to be obtained when using the high-resolution LRT methods. To clarify the reason of that, the L-curves are investigated for case I: synthetic signals of narrow wavenumber-band Lamb modes A0 and S0 on a 2 mm-thick phantom plate (Fig. 5) and for case II: 2.5 mm-thick human bone ex vivo (Fig. 6). Strictly speaking, the hyperparameter μ needs to be optimized as a function of f, slowness, penalty term and misfit term. Therefore, different L-curves at different frequencies are shown at 3D space of μ, penalty term and misfit term. The optimal values of μ parameter are searched from 2 -15 to 2 7 in different bandwidths of the interest.

As shown in Fig. 5, for the narrowband case, the L-curves are computed with the slowness range of 0 < 𝑝 <2.56 𝜇s • 𝑚𝑚 -1 and frequency 0.1 < 𝑓 <0.7 MHz. A group of stable L-curves is obtained with values of the hyperparameter 𝜇 = 0.001, 0.001 and 0.015 for 𝑙 2 -norm, 𝑙 1 -norm and Cauchy norm, respectively.

Thus, for Case I (see Fig. 1), both high resolution and noise filtering can be achieved using LRT methods with a fixed μ value for all frequencies.

The L-curves obtained from the signals of a 2.5 mm-thick human bone ex vivo (case II) are depicted in Fig. 6, with the slowness and frequency ranges of 0 < 𝑝 <2.56 𝜇s • 𝑚𝑚 -1 and 1.1 < 𝑓 <1.6 MHz.

According to the L-curves, small hyperparameter values of 𝜇 < 1 are still preferable for this case.

However, different from Fig. 5, even using such a small 𝜇 value, most of the penalty terms are still obtained with the values below 0.01, which actually indicates the "all-zero" solutions. Such a result explains why in Fig. 4, mode enhancement cannot be achieved by using the LRT methods in bandwidth of 1.1 < 𝑓 < 1.6

MHz.

(

2) SNR threshold of the SVD-based method

The performance of the SVD-based method mainly depends on the singular value selection, which can be optimized by a SNR threshold. Fig. 7 presents the 5 singular values 𝜎 𝑖 as functions of the frequency in dB scale, where the signals are obtained in (a) case I: synthetic signals, narrow wavenumber-band Lamb modes A0 and S0 on a 2 mm-thick Sawbones plate (see Fig. 1), and (b) case II: 2.5 mm-thick human bone ex vivo (see Fig. 4). Our typical experimental signals are recorded with general SNR around 60 dB, but for some low-amplitude modes, the SNR can be less than 10 dB, for example, the 𝜎 4 (𝑓) in dash line between 0.5 MHz and 1 MHz (Fig. 7b). It can be found that in both cases, SNR thresholds can be heuristically determined from the 𝜎 𝑖 functions.

B. Application Condition

For guided signals measured on a relatively large reception length (> 50 mm), the 2D-FT method is able to characterize the dispersion curves of several fundamental Lamb modes, e.g., S0, A0 and S1 in the plate model, or longitudinal guided modes, e.g., L(0,1), L(0,2), and L(0,3) in the cylindrical model [START_REF] Lefebvre | Development of a new ultrasonic technique for bone and biomaterials in vitro characterization[END_REF][START_REF] Ta | Identification and analysis of multimode guided waves in tibia cortical bone[END_REF] . However, the limited spatial sampling achievable with clinical probes (typically a few cm-long transducer arrays) results in a poor wavenumber resolution with the consequence that only high-amplitude and non-overlapped modes can be readily identified using the classical 2D-FT method. Such a limitation of the 2D-FT method can be improved by using the LRT methods and SVD-based method. However, different application conditions of the two methods should be taken into account.

(

1) LRT methods

The merits of the LRT methods originate from the theory of the Radon transform. The linear-path functions, i.e., wave-packets in rays, can be projected to the Radon domain as different energy foci along the linear slant-stack operator. Furthermore, the inversion based sparse technique, so-called high-resolution

LRT method, can be employed to sharpen the 𝜏-𝑝 and 𝑘-𝑓 resolution. As a result, only the less dispersive modes with clear temporal rays can be concentrated as foci by using high-resolution LRT methods with sparsity.

For case I (see Fig. 1), if the wave-packets actually propagate at significantly different velocities and if the temporal overlapping is mainly caused by the short propagation distance instead of dispersion, then the modes can be perfectly separated by the slowness range selection in the Radon domain. Furthermore, considering the reversibility of the LRT method between (𝑘, 𝑓) and (𝑥, 𝑡) fields, the LRT method is capable of providing another good solution to separate some narrowband modes, for example, to extract the slowest fundamental A0 modes in long bone [START_REF] Tran | Multichannel filtering and reconstruction of ultrasonic guided wave fields using time intercept-slowness transform[END_REF] . The extraction of A0 mode in 𝜏-𝑝 field might be more efficient than the temporal wave-packets extraction using the so-called group velocity mask filtering [START_REF] Moilanen | Measuring guided waves in long bones: Modeling and experiments in free and immersed plates[END_REF] .

However, for case II (Figs. 2 and4), i.e., wideband multimodal signals with high attenuation and dispersion, it has been shown that the LRT methods can only enhance the resolution of some high-amplitude modes close to the center frequency of the probe, which actually fails to achieve a wideband dispersion curves extraction.

For case III (Fig. 3), the axial transmission signals, measured from the 25.5 mm-thick PMMA plate, mainly consist of non-dispersive body waves, which are similar to the seismic signals measured from the large-scale media with multiple ray-paths. The LRT methods are able to extract the individual temporal wave-packets propagating along different paths, even for the reflected signals with low amplitudes. The slowness of each wave-packet can be directly read from the τ-p domain. The results suggest that the LRT methods are suitable for extracting FAS [START_REF] Foldes | Quantitative ultrasound of the tibia: a novel approach for assessment of bone status[END_REF] and other multipath body waves [START_REF] Le | Probing long bones with ultrasonic body waves[END_REF] in the long bone.

A beamforming and angle steering strategy at the emission stage can be used to obtain the ultrasonic guided modes with narrowband phase velocity spectrum leading to relatively clear ray contributions in the distance-time diagram [START_REF] Ditri | Excitation of guided waves in generally anisotropic layers using finite sources[END_REF] . For instance, the phase velocities of the multimode signals are approximately in a range of 3 to 5 𝜇s • 𝑚𝑚 -1 28 . In such case, although mode dispersion and temporal overlapping exist, because of the presence of relatively clear rays for the different modes, the high-resolution LRT method can still be used to improve the resolution of 𝑘-𝑓 dispersion curves imaging. However, without enough enhancement of the low-amplitude modes, the LRT methods usually provide identical maxima loci in comparison to 2D-FT. Regarding the identification of those modes with high dispersion and weak amplitude in wide 𝑘-𝑓 ranges, e.g., S1 and S2 modes in Fig. 4, some improvements are still necessary.

We found that (1) the high-resolution LRT methods can enhance the 𝑘-𝑓 resolution of some modes with narrow velocity range, e.g., two fundamental Lamb modes (S0 and A0) in Figs.1 e-f and S0, A0 and A1 and S1 modes in Fig. 4e; (2) in contrast, for wideband highly-dispersive and low-amplitude modes, it is still challenging to concentrate 𝑘-𝑓 trajectories using the LRT methods, e.g., for guided waves signals in

Figs. 2 and4 corresponding to case II. The sparse assumption of the Radon projection of the linear events in the (𝑥, 𝑡) field is well satisfied when different wave-packets propagate at constant velocities (Fig. 3 case III). For the multimode signals with severe dispersion, the assumption is valid, when there are still clear linear events in (𝑥, 𝑡) domain, for instance, signals with only S0 and A0 modes in Fig. 1 (case I) and some seismic data presented in [START_REF] Luo | Rayleigh-wave dispersive energy imaging using a high-resolution linear Radon transform[END_REF][START_REF] Wang | Effect of near-surface topography on high-frequency Rayleigh-wave propagation[END_REF][START_REF] Xia | Generating an image of dispersive energy by frequency decomposition and slant stacking[END_REF] . However, for case II, both the dispersion and short propagation distances of a few centimeters are responsible for modes overlapping for wavenumbers ranging from 0 to 5 𝑟𝑎𝑑 • 𝑚𝑚 -1 10,28,29 . As a consequence, there is no clear linear events observed in the wideband multimode signals (see Figs. 2 and4), so that the low-amplitude signals cannot be effectively enhanced by using the slant-stack operator of the LRT methods. Such a challenge causes the inefficiency of the sparse penalty term, i.e., norm of the 𝑊(𝑝, 𝑓), involved in the LRT methods (Fig. 6). It could explain the difficulties during our application of the LRT methods for extracting the wideband dispersion curves, in particular for the low-amplitude multimode signals under a poor SNR (case II Figs. 2 and4).

(

2) SVD-based method

For both narrowband (case I) and wideband (case II) guided signals, the multi-emitter and multi-receiver configuration combined with the singular vectors selection strategy of the SVD allows achieving a stable performance for noise filtering and extraction of the dispersion curves.

C. Other Potential Approaches and Improvements

Many classical spectra estimation methods [START_REF] Castanié | Spectral Analysis: parametric and non-parametric digital methods[END_REF] , such as the Burg method or the multiple signal classification method (MUSIC), can be used to achieve high-resolution wavenumber estimation. Recently, sparse methods have been introduced for dispersion curves extraction. Harley et al. [START_REF] Harley | Sparse recovery of the multimodal and dispersive characteristics of Lamb waves[END_REF] have proposed a compressed-sensing-based sparse wavenumber method for the recovery of the dispersion curves. These authors showed that the sparse penalty regularization can be directly performed using the wavenumber penalty. Such a sparse strategy in the (𝑘, 𝑓) field may be more efficient and convenient. Currently, the sparse wavenumber extraction method proposed by Harley et al. has been verified on aluminum metal plates [START_REF] Harley | Sparse recovery of the multimodal and dispersive characteristics of Lamb waves[END_REF] . The practical challenge encountered with axial transmission in cortical bone is to efficiently enhance the weak modes under the conditions of severe broadband overlapping (more than 5 modes in some frequency band, see Fig. 4). The sparse singular value decomposition (S-SVD) technique, i.e., an improved SVD-based method by using the sparse strategy in (𝑘, 𝑓) rather than the (𝜏, 𝑝) domain, may significantly overcome the limitation of poor wavenumber resolution [START_REF] Xu | Sparse SVD Method for High Resolution Extraction of the Dispersion Curves of Ultrasonic Guided Waves[END_REF][START_REF] Xu | Sparse inversion SVD method for dispersion extraction of ultrasonic guided waves in cortical bone[END_REF] . Other improvements for 3D multi-emitter and multi-receiver space-time signal processing, e.g., high-dimensional seismic data processing [START_REF] Yu | Interpolation and denoising of high-dimensional seismic data by learning a tight frame[END_REF] , might be also helpful for guided waves dispersion analysis, but the application has not been reported in community of ultrasonic bone evaluation to date. In addition, other guided modes excitation technology, e.g., coded excitation [START_REF] Song | A base-sequence-modulated Golay code improves the excitation and measurement of ultrasonic guided waves in long bones[END_REF][START_REF] Lin | Excitation Waveform Design for Lamb Wave Pulse Compression[END_REF][START_REF] Yucel | Coded waveform excitation for high resolution ultrasonic guided wave response[END_REF] and wideband dispersion reversal method [START_REF] Xu | Wideband dispersion reversal of Lamb waves[END_REF] etc., can also be helpful to enhance the SNR of the ultrasonic axial transmission signals in the long cortical bone.

Generally speaking, in order to interpret the relatively complex guided signals, the signal processing methods should be robust enough to allow the wideband dispersion curve extraction and the low-amplitude mode detection. By retaining the singular values above noise level, the SVD-based method significantly enhances the weak mode extraction when they are poorly detected by the 2D-FT method and LRT methods.

In this sense, the SVD-based method could be more suitable for signal processing of the ultrasonic guided waves in the long bone, especially for highly dispersive wideband signals, in presence of severe attenuation and low SNR.

Figure 2

 2 presents the experimental signals measured in a 4 mm-thick bone-mimicking plate. Fig. 2a is the distance-time diagram of the array-signal. As shown by the 2D-FT and SVD 𝑘-𝑓 results (Figs. 2b-c), the detectable wavenumber dispersion is in the range of 0 < 𝑘 < 4 𝑟𝑎𝑑 • 𝑚𝑚 -1 with more than 5 modes. The experimental SNR is around 60dB. Those singular vectors associated with singular values higher than an heuristic threshold of 20dB were remained as the signal subspaces and the rest were filtered out as noise. Fig. 2d plots the 𝑙 2 -norm-based energy distributions in the (𝜏, 𝑝) domain. Figs 2e-f depict the high-resolution LRT results of the multimode energy distribution in 𝑘-𝑓 field using 𝑙 1 -norm and the Cauchy norm.

  . A 25.5 mm-thick PMMA plate was measured to obtain the signals mainly consisting of ultrasonic body waves. The experimental SNR is as the previous example around 60dB. As ultrasonic body waves have been investigated in a 6.5 mm-thick bovine bone ex vivo[START_REF] Le | Probing long bones with ultrasonic body waves[END_REF] , such a very-thick plate is also prepared to clarify the performance of the two methods in the high 𝑓 • ℎ range. It should be noted that 25.5 mm-thick waveguides are unlikely to be encountered in human cortical bone whose thickness varies from less than a millimeter to a few millimeters at best. Different from the highly dispersive guided modes in the thin plates (TableI, case II), in a thick plate, there mainly exist body waves propagating as temporally separated wave-packets. As shown in Fig.3a, there are five different wave-fronts in the distance-time diagram with different ray paths. Pd and Sd represent the longitudinal and shear waves axially propagating along the plate surface. PrP and SrS correspond to the first reflection of the longitudinal and shear waves on the bottom wall. PrS corresponds to the longitudinal-to-shear wave conversion when the P wave is reflected on the bottom surface of the plate.The markers on Fig.3awere computed according to the distance-velocity relationship. Figs. 3b-c are the 2D-FT and SVD 𝑘-𝑓 results, respectively. After projecting the array-signal from distance-time domain to Radon field, accurate slowness and amplitude of each wave-packet can be obtained from those energy focusing maxima in τ-p domain, which is convenient for detecting the components with different ray paths.

Figs. 4e-

  Figs.4e-fshow the 𝑘-𝑓 energy distribution obtained by 𝑙 1 -norm and Cauchy norm high-resolution LRT methods. Comparing with the LRT and 2D-FT method, SVD-based method is capable of detecting the noise polluted A0 and S0 mode and also part of A1, A2, S1 and S2 modes in wideband frequency-thickness range. Similar to the 4 mm-thick phantom signals (Fig.2), because of high dispersive characteristics, the wideband guided modes in the human radius cannot be readily separated and enhanced in the Radon field.

  methods are necessary to cover the entire cortical bone thickness range of the human long bone. The LRT methods have the advantage of the reversibility between (𝜏, 𝑝) and (𝑥, 𝑡) fields showing a good ability to separate modes with large velocity difference, which is suitable for data analysis of ultrasonic guided waves at low 𝑓 • ℎ range (𝑓 • ℎ < 1 MHz • 𝑚𝑚), e.g., signals consisting of two fundamental modes S0 and A0 with a large velocity difference presented in the study (Fig.1). For the highly dispersive multimode signals in a broadband 𝑓 • ℎ range (0<𝑓 • ℎ < 6 MHz • 𝑚𝑚), which are quite

Fig. 1 .

 1 Fig. 1. (color online) Synthetic signals, narrow wavenumber-band Lamb modes A0 and S0 on a 2 mm-thick Sawbone plate with peak-to-peak amplitudes of 1 and 0.3, and SNR of 30dB, (a) distance-time diagram of the array-signals, (b) 2D-FT 𝑘-𝑓 result, (c) SVD 𝑘-𝑓 result, (d) the 𝜏-𝑝 energy distributions obtained by LRT with 𝑙 2 -norm, and 𝑘-𝑓 results obtained by LRT with high-resolution regularization strategies, i.e., (e) 𝑙 1 -norm, (f) Cauchy norm, respectively.

Fig. 2 .

 2 Fig. 2. (color online) Experimental signals measured in a 4 mm-thick bone-mimicking plate, (a) distance-time diagram of the array-signals, (b) 2D-FT 𝑘-𝑓 result, (c) SVD 𝑘-𝑓 result, (d) the 𝜏-𝑝 energy distributions obtained by LRT with 𝑙 2 -norm, and 𝑘-𝑓 results obtained by LRT with high resolution regularization strategies, i.e., (e) 𝑙 1 -norm, (f) Cauchy norm, respectively.

Fig. 3 .

 3 Fig. 3. (color online) Experimental signals measured in the 25.5 mm-thick PMMA plate, (a) distance-time diagram of the synthetic array-signals, (b) 2D-FT 𝑘-𝑓 mode energy distribution, (c) SVD k-f mode energy distribution, (d) the 𝜏-𝑝 energy distributions obtained by LRT with 𝑙 2 -norm, and 𝑘-𝑓 energy distributions obtained by LRT with high resolution regularization strategies, i.e., (e) 𝑙 1 -norm, (f) Cauchy norm, respectively.

Fig. 4 .

 4 Fig. 4. (color online) Experimental guided signals measured in the 2.5 mm-thick ex-vivo human radius, (a) distance-time diagram of the synthetic array-signals, (b) 2D-FT 𝑘-𝑓 mode energy distribution, (c) SVD 𝑘-𝑓 mode energy distribution, (d) the 𝜏-𝑝 energy distributions obtained by LRT with 𝑙 2 -norm, and 𝑘-𝑓 energy distributions obtained by LRT with high resolution regularization strategies, i.e., (e) 𝑙 1 -norm, (f) Cauchy norm, respectively.

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. (color online) L-curves of the synthetic signals (case I shown in Fig. 1), i.e., narrow wavenumber-band Lamb modes A0 and S0 on a 2-mm-thick Sawbones plate, (a) 𝑙 2 -norm, (b) 𝑙 1 -norm, (c) Cauchy norm. The slowness 𝑝 and frequency f range from 0 to 2.56 𝜇s • 𝑚𝑚 -1 and 0.1 to 0.7 MHz, respectively.

Table

  1: (1) low frequency-thickness product 𝑓 • ℎ < 1 MHz• mm, corresponding to thin cortical bone about 1 mm to 4 mm (mainly for radius), where mainly two fundamental Lamb modes S0 and A0 are observed, (2) intermediate 𝑓 • ℎ range (1 <𝑓 • ℎ<6 MHz• mm) where multiple guided modes overlap without temporally separated wave-packets, and (3) high 𝑓 • ℎ > 6 MHz• mm, corresponding to the thick cortical bone, e.g., tibia, which is mainly in the range of body waves propagation. Various signal processing techniques have been reported to identify and disentangle the overlapping modes for short axial propagation distances. Table II lists the different signal processing approaches proposed for axial transmission, together with a brief description of the methods and conditions of application.

  3D multi-emitter and multi-receiver signals matrix is denoted as 𝑀(𝐸, 𝑥, 𝑡), where 𝑥 is the distance sampled by 𝑁 𝑟 receivers, and 𝐸 denotes different emitters, 𝐸 = 1, 2, … , 𝑁 𝑒 . ℳ(𝐸, 𝑥, 𝑓) is the spectrum of 𝑀(𝐸, 𝑥, 𝑡) in frequency domain. After the SVD decomposition of the 2D response matrix ℳ(𝐸, 𝑥, 𝑓 = 𝑓 0 ) at each frequency 𝑓 0 , we obtain two unitary matrices, i.e., an 𝑁 𝑒 × 𝑁 𝑒 matrix 𝑈 and an 𝑁 𝑟 × 𝑁 𝑟 matrix 𝑉, and an 𝑁 𝑒 × 𝑁 𝑟 rectangular matrix ∑ with the singular values as the diagonal entries 𝜎 𝑖 (𝑖 = 1, … , 𝑁 𝑒 ).The 𝑁 𝑟 columns of 𝑉, i.e., 𝑣 𝑖 (𝑖 = 1, … , 𝑁 𝑟 ), are the reception singular vectors. At a particular frequency, only the first 𝑁 𝑟𝑘 singular vectors 𝑣 𝑖 associated with the most energetic singular values 𝜎 𝑖 (𝑖 = 1, … , 𝑁 𝑟𝑘 )

  , experiments are carried out in two phantoms, i.e., a 4 mm-thick bone-mimicking material (Sawbones, Pacific Research Laboratory Inc, Vashon, WA) and 2.5 mm-thick ex-vivo human radius. A 25.5 mm-thick polymethylacrylate (PMMA) plate is also measured for case III in Table1.Ultrasound gel (Aquasonic, Parker Labs, Inc, Fairfield, NJ) is used to ensure the coupling between the probe and the sample. The characteristics of the specimens are listed in table III. 𝑉 𝑇 is the shear velocity, and 𝑉 𝐿‖ , and 𝑉 𝐿⊥ are the pure compression bulk wave velocities in the direction parallel and normal to the direction of the fibers (bone-mimicking material) or to the long axis of the bone (human radius). The

density and thickness are denoted by ρ and h. Typical values for human radius density, shear and longitudinal velocity, derived from the literature, are used for computation of the theoretical dispersion curves of the 2D transverse isotropic free plate model

[START_REF] Foiret | Combined estimation of thickness and velocities using ultrasound guided waves: a pioneering study on in vitro cortical bone samples[END_REF]

. The average thickness of the human radius specimen was obtained by X-ray computed tomography (XtremCT, Scanco Medical, Bruttisellen, Switzerland).

Table II . Signal processing methods for assessment of the long cortical bone using axial transmission

 II

	frequency-thickness	Case I: Low range	Case II: Middle range	Case III: High range
	product (𝑓 • ℎ) (MHz • 𝑚𝑚)	𝑓 • ℎ < 1			1 < 𝑓 • ℎ < 6	6 < 𝑓 • ℎ
	2D-FT 15	1. Separation of several		N/A	N/A
		low-order modes from
		narrowband	signals	1. Under the wideband	1. Similar to many
		measured in bovine bone		excitation, more than	geophysical
		ex vivo. 43,44			5 guided modes with	applications, mainly
	SVD 8	1. Enhancement of weak modes;	overlapping velocity	body waves with N/A
		ranges; 2. Applicable to extract dispersion curves from	linear ray paths 9 ;
	High-resolution LRT 32	2. For the relatively short propagation distances (a few cm) 8,12,18,44,57 , a complete wideband guided wave signals without a priori knowledge; 3. Extraction of wideband attenuation for each guided mode 27 ; 4. Bone-mimicking phantoms 38,59,60 , ex vivo 14 and in dispersion extraction vivo measurements of human bones 10 . of the overlapping 1. High-resolution solution of the mode trajectories multimode is still basing on a linear slant stack of the distance-time challenging. data 28,29 ;	2. FAS is the lateral wave propagating at the bulk compression velocity 7,58 ; 3. Velocities of S0 and A0 modes converge 1. Thick bones, e.g., to the Rayleigh tibia, have not been velocity. Velocities of tested yet,
		2. Cervine tibiae 28 with narrow k-band angle beam	other guided modes 2. General acceptance in
		excitation; 3. Extraction of the slowest A0 mode from wideband	converge to the shear seismic field for velocity 35 . multichannel
		3. The multimode signals in bovine femur in vitro 29 . dispersion	processing of signals
		information of the		with	ray
		lateral arrival A0 mode		characteristics	in
		can be used to estimate cortical thickness 11,56 .		distance-time domain 39 .
	Group velocity
	mask filtering 45
	Signal processing methods

Table III . Shear and longitudinal velocities, density and thickness of the specimens used in experiments

 III : to the best knowledge of the authors, it has not been reported in the literatures of cortical bone guided waves processing.

	Time-frequency extraction 62 Joint approximate diagonalization of eigen-matrices (JADE) method 63 Bone-mimicking plate 27 PMMA 27 N/ASpecimens Human radius specimen 14	1. Bone-mimicking phantoms 61 , ex vivo human long bones 17,56 ; 2. Extraction of A0 mode component from narrowband guided waves signals. 1. Wideband A0 mode separation without a prior N/A velocity knowledge. 2. Ex-vivo measurement of bovine tibiae 20,44 and of human long bone 11 . N/A 1. Blind identification, has been used to separate the narrowband guided waves for cortical thickness determination of bovine tibiae 13 . ρ (𝑔 • 𝑐𝑚 -3 ) h (mm) 𝑉 𝑇 (𝑚𝑚 • 𝜇𝑠 -1 ) 1.64 4 1.62 1.18 25.5 1.37 1.85 1.58 1.8	N/A N/A N/A 𝑉 𝐿 (𝑚𝑚 • 𝜇𝑠 -1 ) (𝑉 𝐿‖ , 𝑉 𝐿⊥ ) = (3.57, 2.91) 2.7 (𝑉 𝐿‖ , 𝑉 𝐿⊥ ) = (4.0, 3.41)

List of figure captions:

J. Acoust. Soc. Am. 140 (3), September 2016 Xu et al.: Dispersion Curves Extraction of Guided Waves 20

J. Acoust. Soc. Am. 140 (3), September 2016 Xu et al.: Dispersion Curves Extraction of Guided Waves 21

usual in the axial transmission measurement of the human long bone (Figs. 2 and4), the SVD-based method shows more robust performances for weak mode enhancement and noise filtering. Finally, regarding computation time, the SVD-based method can be accomplished efficiently without any iterations, but the 𝑙 1 -norm and Cauchy norm LRT methods are relatively time-consuming due to the reweighting strategy at each frequency. Future work will be to improve the resolution of the SVD method using the sparse strategy, i.e., recently proposed S-SVD method [START_REF] Xu | Sparse SVD Method for High Resolution Extraction of the Dispersion Curves of Ultrasonic Guided Waves[END_REF][START_REF] Xu | Sparse inversion SVD method for dispersion extraction of ultrasonic guided waves in cortical bone[END_REF] , which may achieve a high-resolution extraction of the dispersion curves of ultrasonic guided waves.

ACKNOWLEDGMENT

The authors also acknowledge Dr. Maryline Talmant for her constructive comments and excellent advices to improve the paper. This work was supported by the National Natural Science Foundation of China (11304043, 11327405, 11525416), the NSFC-CNRS-PICS (11511130133) and the CNRS PICS programme n°07032.

List of tables:


Patricia Bouyer

Nicolas Markey

Nicolas Perrin

Philipp Schlehuber-Caissier

Timed-automata abstraction of switched dynamical systems using control funnels

The development of formal methods for control design is an important challenge with potential applications in a wide range of safety-critical cyber-physical systems. Focusing on switched dynamical systems, we propose a new abstraction, based on time-varying regions of invariance (the control funnels), that models behaviors of systems as timed automata. The main advantage of this method is that it allows automated verification of formal specifications and reactive controller synthesis without discretizing the evolution of the state of the system. Efficient constructions are possible in the case of linear dynamics. We demonstrate the potential of our approach with two examples.

Introduction

Verification and synthesis are notoriously difficult for hybrid dynamical systems, i.e. systems that allow abrupt changes in continuous dynamics. For instance, reachability is already undecidable for 2-dimensional piecewise-affine maps [START_REF] Koiran | Computability with low-dimensional dynamical systems[END_REF], or for 3-dimensional dynamical systems with piecewise-constant derivatives [START_REF] Asarin | Reachability analysis of dynamical systems having piecewise-constant derivatives[END_REF].

To enable automated logical reasoning on switched dynamical systems, most methods tend to entirely discretize the dynamics, for example by approximating the behavior of the system with a finite-state machine. Alternatively, early work pointed out links between hybrid and timed systems [START_REF] Maler | From timed to hybrid systems[END_REF], and several methods have been designed to create formal abstractions of dynamical systems that do not rely on a discretization of time. In [START_REF] Frazzoli | Maneuver-based motion planning for nonlinear systems with symmetries[END_REF], a finite maneuver automaton is constructed from a library of motion primitives, and motion plans correspond to timed words. In [START_REF] Liu | Switching control of dynamical systems from metric temporal logic specifications[END_REF][START_REF] Fu | Computational methods for stochastic control with metric interval temporal logic specifications[END_REF], switched controller synthesis and stochastic optimal control are realized via metric temporal logic (MTL) or metric-interval temporal logic (MITL) specifications. In [START_REF] Quottrup | Multi-robot planning : a timed automata approach[END_REF][START_REF] Maler | Approximating continuous systems by timed automata[END_REF], grid-based abstractions and timed automata are used for motion planning or to check timed properties of dynamical systems. In [START_REF] Sloth | Complete abstractions of dynamical systems by timed automata[END_REF], a subdivision of the state space created from sublevel sets of Lyapunov functions leads to an abstraction of dynamical systems by timed automata that enables verification and falsification of safety properties. The same kind of abstraction is used in [START_REF] Sloth | Timed game abstraction of control systems[END_REF] for controller design via timed games, but the update map of the timed games obtained is such that synthesis cannot be realized using existing tools. In [START_REF] David | Application of model-checking technology to controller synthesis[END_REF], the state space of each mode of a piecewise-affine hybrid system is portioned into polytopes, and thanks to control laws that prevent the system from exiting through a given facet, or that force the system to exit through a facet in finite time, reactive control problems can be solved as timed games on timed automata.

Our contribution is a novel timed-automata abstraction of switched dynamical systems based on a particular kind of time-varying regions of invariance: control funnels. Recent results have shown that these invariants are very useful for robust motion planning and control [START_REF] Tedrake | LQR-trees: Feedback motion planning via sums-of-squares verification[END_REF][START_REF] Majumdar | Robust online motion planning with regions of finite time invariance[END_REF][START_REF] Majumdar | Control design along trajectories with sums of squares programming[END_REF], and that funnels or similar concepts related to the notion of Lyapunov stability can be used for formal verification of hybrid systems [START_REF] Julius | Trajectory based verification using local finite-time invariance[END_REF][START_REF] Duggirala | Verification of annotated models from executions[END_REF], and for reactive controller synthesis [START_REF] Decastro | Synthesis of nonlinear continuous controllers for verifiably-correct high-level, reactive behaviors[END_REF].

The paper is organized as follows: Section 2, describes how control funnels, especially for trajectory tracking controllers, can be used to create timed transition systems that abstract the behavior of a given switched dynamical system, and as a result can potentially allow the use of verification tools for motion planning. In Section 3, we show how these timed transition systems can be encoded as timed automata. In Section 4, we consider the case of linear dynamics and introduce the notion of fixed size LQR funnel. In Section 5, we present two examples of application and efficient algorithms that manipulate these LQR funnels. In the first one, a timed game is solved by the tool Uppaal-Tiga [START_REF] Behrmann | UPPAAL-Tiga: Time for playing games! In[END_REF] for the synthesis of a controller that can reactively adjust the phase of a sine wave controlled in acceleration. In the second example, we show that, using our timed-automata abstraction with LQR funnels along constant velocity trajectories, a non-trivial solution to a pick-and-place problem can be computed by the model checker Uppaal [START_REF] Behrmann | Uppaal 4.0[END_REF]. Section 6 concludes and presents avenues for future work.

Graphs of control funnels

Control funnels

Consider a controlled dynamical system governed by the following equation:

ẋ = f (x, t, u(x, t)), (1)
where x ∈ R d is the state of the system (which can contain velocities1 ), t ∈ R + is a real (clock) value corresponding to time (we restrict ourselves to nonnegative time values), u :

R d × R + → R k is the control input function, and f is a continuously differentiable function from R d × R + × R k to R d (
which ensures the uniqueness of the solution for given initial conditions). Assuming that the function u is fixed, we also use the following notation for Equation (1):

ẋ = f u (x, t). (2)
Fig. 1. An example of control funnel for a controller tracking a reference trajectory.

The dashed line is a trajectory of the controlled system in the state space. On the right side, switching transitions between control funnels are depicted.

A control funnel for the above dynamical system is a function F : I → 2 R d such that I ⊆ R + and for any solution x(t) of ( 2), the following property holds:

∀t 1 ∈ I. ∀t 2 ∈ I. (t 2 > t 1 and x(t 1 ) ∈ F(t 1 )) ⇒ x(t 2 ) ∈ F(t 2 ). (3)
It corresponds to time-varying regions of invariance.

Example 1. A typical example of a control funnel based on a trajectory tracking controller (that is, a control funnel asymptotically converging towards a reference trajectory in the state space) is shown in Fig. 1.

Example 2. For a concrete example, consider the simple system whose trajectories are of the form e -t • x 0 . Then any set W ⊆ R d defines a control funnel

F W : t → {e -t • w | w ∈ W }.
The notion of funnel was popularized by Mason [START_REF] Mason | The mechanics of manipulation[END_REF], and it usually specifically refers to operations that eliminate uncertainty (as is the case in the example of Fig. 1) by collapsing a large set of initial conditions into a smaller set of final conditions (see for instance [START_REF] Tedrake | LQR-trees: Feedback motion planning via sums-of-squares verification[END_REF]). In our case, the control funnel may or may not reduce uncertainty, and it is important to note that the set F(t) does not have to decrease in size over time. This more general concept is closer to the definition of viability tubes [START_REF] Aubin | Viability tubes[END_REF], but we nevertheless use the term control funnel as some reduction of uncertainty is often essential to the usefulness of our abstractions. We address the computation of control funnels in Section 4, and leave them as relatively abstract objects for now.

Motion planning

Let us suppose that we have a finite set U of control laws u 1 (x, t), u 2 (x, t), . . . , u n (x, t) that respectively set the dynamics of a given system to ẋ = f u1 (x, t), ẋ = f u2 (x, t), . . . , ẋ = f un (x, t).

We say that the system can switch to the control law u i (x, t) at some state x whenever there is t 0 ∈ R + and a solution x(t) of ẋ = f ui (x, t) with initial condition x = x(t 0 ). Typically, if u i (x, t) corresponds to a trajectory tracking controller, t 0 identifies the point of the trajectory where the tracking is triggered.

Informally, the motion planning problem asks, given a finite set of control laws as above, an initial point x 0 , a target zone T f ⊆ R d , and an obstacle Ω ⊆ R d , whether there exists a sequence of control law switches that generates a trajectory from x 0 to T f while avoiding the obstacle Ω. Several variants of this problem can be considered, that vary on the objective (for instance some tasks can be expressed as ω-regular objectives), but we focus here on a reachability with avoidance objective.

More formally, motion planning asks for a finite sequence of time values t 1 0 < t Intuitively, this means that we can switch conveniently between all the control laws, causing discrete changes in the system dynamics, and ensure the global (reachability with avoidance) objective. The continuous trajectory generated by the solution above is the concatenation of the trajectory portions

{x p (t) | t p 0 ≤ t ≤ t p 1 } for 1 ≤ p ≤ P .

Motion planning with graphs of control funnels

We now explain how the motion planning problem can be abstracted using timed transition systems based on control funnels. For each control law u i (x, t), we assume that we have a finite set of control funnels F 1 i , F 2 i , . . . , F mi i , respectively defined over

I 1 i ⊆ R + , I 2 i ⊆ R + , . . . , I mi i ⊆ R + .
We assume that for every 1 ≤ i ≤ n, for every 1 ≤ j ≤ m i , for every t ∈ I j i , it holds F j i (t) ∩ Ω = ∅, which means that trajectories contained in these funnels always avoid the obstacle Ω.

Consider a control law switch at position x to law u i (x, t) with clock value t 0 . If there exists a control funnel F j i such that t 0 ∈ I j i , and x ∈ F j i (t 0 ), then we know that the state of the system will remain inside F j i (t) for any t > t 0 in I j i (as long as control law u i (x, t) is used). To always keep the system inside one of the control funnels, we can impose sufficient conditions on the switches. For instance, if the state is inside F j i (t 0 ), and if for some future clock value t 1 , there exists a control funnel F l k and t 2 ∈ I l k such that F j i (t 1 ) ⊆ F l k (t 2 ), then when the clock value is t 1 we can safely switch to the control law u k (x, t) while setting the clock to t 2 . Indeed, we know that the state of the system at the switch instant will be inside F l k (t 2 ), and therefore it will remain inside F l k (t) after the switch. Such transitions from a funnel to another are illustrated on the right side of Fig. 1. It is worth noting that similar transitions could be achieved with, instead of control funnels, controller specifications as introduced in [START_REF] Le Ny | Sequential composition of robust controller specifications[END_REF]. For some control funnels F j i and F k i associated to the same control law, it is the case (see Section 4) that when funnel F j i is entered at time t, then at any time

t ≥ t + h j→k i
(where h j→k i is a constant), the state of the system is inside F k i (t ). In that case, we say that the funnel F k i h j→k i -absorbs the funnel F j i . These rules for navigating between control funnels give to the set of control funnels the structure of an infinite graph, or, more precisely, of a timed transition system with real-valued clocks. One of the clocks of this timed transition system is c t , the controller clock. We add two other clocks: a global clock c g , and a local clock c h .

The funnel timed transition system T U,F associated with the family of laws U = (u i (x, t)) 1≤i≤n and the family of funnels F = ((F j i , I j i )) 1≤i≤n,1≤j≤mi is defined as follows. The configurations are pairs (F j i , v) where v assigns a nonnegative real value to each of the clocks c t , c g and c h , with v(c t ) ∈ I j i , and its transition set contains three types of elements:

the time-elapsing transitions:

(F j i , v) → (F j i , v + ∆) whenever [v(c t ), v(c t ) + ∆] ⊆ I j i (where v + ∆ denotes the valuation that maps each clock c to v(c) + ∆); -the switching transitions: (F j i , v) → (F l k , v ) whenever v (c g ) = v(c g ), v (c h ) = 0, v(c t ) ∈ I j i , v (c t ) ∈ I l k , and
F j i (v(c t )) ⊆ F l k (v (c t )); -the absorbing transitions: (F j i , v) → (F k i , v ) whenever F k i h j→k i -absorbs F j i , v(c h ) ≥ h j→k i , v (c h ) = 0, v (c g ) = v(c g ) and v (c t ) = v(c t ).
A run in this timed transition system is a finite sequence of configurations (F j0 i0 , v 0 ), (F j1 i1 , v 1 ), . . . , (F

j P i P , v P ) such that v 0 (c h ) = v 0 (c g ) = 0, v 0 (c t ) ∈ I j0 i0
, and all the transitions (F jp ip , v p ) → (F jp+1 ip+1 , v p+1 ) for 0 ≤ p < P are valid transitions that belong to T U,F . Such a run is of total duration v P (c g ), and it corresponds to a schedule of control law switches that results from the following rules: initially, the control law is set to u i0 (x, t), and the controller clock c t is set to v 0 (c t ). For every timeelapsing transition (F j i , v) → (F j i , v + ∆), the same control law u i (x, t) is kept for a duration of ∆ time units, and for every switching transition (F j i , v) → (F l k , v ), the control law is switched from u i (x, t) to u k (x, t), with an initialization of the controller clock to v (c t ). Absorbing transitions are discarded, as they just correspond to an instantaneous change of funnels for the same control law. Let us denote this sequence of switches by r. Then, it is fundamental to notice that for every x ∈ F j0 i0 (v 0 (c t )), if we follow the schedule of control law switches just described, then the system remains inside control funnels and reaches at the end of the run a unique point of R d , that we denote r(x). The trajectory going from x to r(x) is also uniquely defined.

The funnel timed transition system T U,F satisfies the following property:

Theorem 1. Let r = (F j0 i0 , v 0 ), (F j1 i1 , v 1 ), . . . , (F j P i P , v P ) be a run in T U,F . If x ∈ F j0 i0 (v 0 (c t )), then r(x) ∈ F j P i P (v P (c t )).
In some sense, the funnel timed transition system T U,F is a correct abstraction of trajectories that can be generated by the set of control laws: if

x 0 ∈ F j0 i0 (v 0 (c t ))
Fig. 2. Run of a funnel timed transition system with three control funnels:

r = (F 1 1 , v 1 0 ), (F 1 1 , v 1 1 ), (F 1 2 , v 2 0 ), (F 1 2 , v 2 1 ), (F 1 3 , v 3 0 ), (F 1 3 , v 3 1 ), (F 1 4 , v 4 0 ), (F 1 4 , v 4
1 ) , with:

∀1 ≤ i ≤ 4, v i 0 (ct) = t i 0 , v i 1 (ct) = t i 1 , v i 0 (c h ) = 0, v i 1 (c h ) = t i 1 -t i 0 , v i 1 (cg) = v i 0 (cg)+v i 1 (c h ), and v 1 0 (cg) = 0, and ∀2 ≤ i ≤ 4, v i 0 (cg) = v i-1 1 (cg).
and F j P i P (v P (c t )) ⊆ T f , then such a run witnesses a solution to the motion planning problem. However, this abstraction is obviously not complete. Example 3 (An example with obstacles). The example in Fig. 2 shows a run with three control laws u 1 (x, t), u 2 (x, t) and u 3 (x, t), three control funnels

F 1 1 , F 1 2
and F 1 3 , and an obstacle in the state space. The domains of definition of the control funnels I 1 1 , I 1 2 and I 1 3 are such that for all α ∈ {1, 2, 3} and all t ∈ I 1 α , F 1 α (t) does not intersect the obstacles. With the previous property, any run in the corresponding funnel timed transition system leads to a trajectory that avoids the obstacles. The example of Fig. 2, where reaching F 1 1 (t 4 1 ) from F 1 1 (t 1 0 ) requires a series of switches between the different control funnels, shows the potential interest of automated verification in timed transition systems, as it can result in the generation of obstacle-free dynamic trajectories via appropriate control law switches.

Reduction to timed automata

Timed automata [START_REF] Alur | A theory of timed automata[END_REF] are a timed extension of finite-state automata, with a wellunderstood theory. They provide an expressive formalism for modelling and reasoning about real-time systems, and enjoy decidable reachability properties; much efforts have been invested over the last 20 years for the development of efficient algorithms and tools for their automatic verification (such as the tool Uppaal [START_REF] Behrmann | Uppaal 4.0[END_REF], which we use in this work).

Let C be a finite set of real-valued variables called clocks. A clock valuation over a finite set of clocks C is a mapping v : C → R + . We write R C for the set of clock valuations over C. If ∆ ∈ R + , we write v + ∆ for the clock valuation defined by (v + ∆)(c) = v(c) + ∆ for every c ∈ C. A clock constraint over C is a boolean combination of expressions of the form c ∼ α where α ∈ Q, and ∼ ∈ {≤, <, =, >, ≥}. We denote by C(C) the set of clock constraints over C. We write v |= g if v satisfies g (defined in a natural way). A reset of the clocks is an element res of (Q ∪ {⊥}) C (which we may write R(C)), and if v is a valuation, its image by res, denoted res(v), is the valuation mapping c to v(c) whenever res(c) = ⊥, and to res(c) ∈ Q otherwise.

We define a slight extension of timed automata with rational constants, general boolean combinations of clock constraints and extended clock resets; those timed automata are as expressive as standard timed automata [START_REF] Bouyer | Updatable timed automata[END_REF], but they will be useful for expressing funnel timed transition systems. A timed automaton is a tuple

A = (L, L 0 , L F , C, E, Inv) where L is a finite set of locations, L 0 ⊆ L is a set of initial locations, L F ⊆ L is a set of final locations, C is a finite set of clocks, E ⊆ L × C(C) × R(C) × L is a finite set of edges, and Inv : L → C(C) is an invariant labelling function. A configuration of A is a pair ( , v) ∈ L × R C such that v |= Inv(
), and the timed transition system generated by A is given by the following two rules:

time-elapsing transition: ( , v) → ( , v + ∆) whenever v + δ ∈ Inv( ) for every 0 ≤ δ ≤ ∆; switching or absorbing transition: ( , v) → ( , v ) whenever there exists ( , g, res,

) ∈ E such that v |= g ∧ Inv( ), v = res(v), and v ∈ Inv( ).
A run in A is a sequence of consecutive transitions. The most fundamental result about timed automata is the following:

Theorem 2 ([1]).
Reachability in timed automata is PSPACE-complete.

We consider again the family of control laws U = (u i (x, t)) 1≤i≤n , and the family of funnels F = ((F j i , I j i )) 1≤i≤n,1≤j≤mi , as in the previous section. For every pair 1 ≤ i, k ≤ n, and every 1 ≤ j ≤ m i and 1 ≤ l ≤ m k , we select finitely many tuples (switch, [α, β], (i, j), γ, (k, l)) with α, β, γ ∈ Q such that (i) [α, β] ⊆ I j i , (ii) γ ∈ I l k , and (iii) for every t ∈ [α, β], F j i (t) ⊆ F l k (γ). This allows us to underapproximate the possible switches between funnels. For every 1 ≤ i ≤ n, for every pair 1 ≤ j, k ≤ m i , we select at most one tuple (absorb, ν, (i, j, k)) such that ν ∈ Q and F k i (t) ν-absorbs F j i (t). This allows us to under-approximate the possible absorbing transitions. For every 1 ≤ i ≤ n and every 1 ≤ j ≤ m i , we fix a finite set of tuples (initial, α, (i, j)) such that α ∈ Q and x 0 ∈ F j i (α). This allows us to under-approximate the possible initialization to a control funnel containing the initial point x 0 . For every 1 ≤ i ≤ n and 1 ≤ j ≤ m i , we fix finitely many tuples (invariant, S i,j , (i, j)), where S i,j ⊆ I j i is a finite set of closed intervals with rational bounds. This allows us to under-approximate the definition set of the funnels. Finally, for every 1 ≤ i ≤ and 1 ≤ j ≤ m i , we fix finitely many tuples (target, [α, β], (i, j)), where α, β ∈ Q and [α, β] ⊆ I j i ∩ {t | F j i (t) ⊆ T f }. This allows us to under-approximate the target zone. We denote by K the set of all tuples we just defined.

We can now define a timed automaton that conservatively computes the runs generated by the funnel timed transition system. It is defined by A U,F,K = (L, L 0 , L F , C, E, Inv) with:

-L = {F j i | 1 ≤ i ≤ n, 1 ≤ j ≤ m i } ∪ {init, stop}; L 0 = {init}; L F = {stop}; -C = {c t , c g , c h };
-E is composed of the following edges:

• for every (initial, α, (i, j)) ∈ K, we have an edge (init, true, res, F j i ) in E, with res(c t ) = α and res(c g ) = res(c h ) = 0; • for every (switch, [α, β], (i, j), γ, (k, l)) ∈ K, we have an edge (F j i , α ≤ c t ≤ β, res, F l k ) with res(c t ) = γ, res(c h ) = 0 and res(c g ) = ⊥; • for every (target, [α, β], (i, j)) ∈ K, we have an edge (F j i , α ≤ c t ≤ β, res, stop) in E, with res(c t ) = res(c g ) = res(c h ) = ⊥; • for every (absorb, ν, (i, j, k)) ∈ K, we have an edge (F j i , c h ≥ ν, res, F k i ) with res(c h ) = 0 and res(c t ) = res(c g ) = ⊥; for every (invariant, S i,j , (i, j)) ∈ K, we let Inv(F j i )

[α,β]∈Si,j (α ≤ c t ≤ β). We easily get the following result:

Theorem 3. Let (init, v 0 ) → ( 1 , v 1 ) → • • • → ( P , v P ) → (stop, v P )
be a run in A U,F,K such that v 0 assigns 0 to every clock. Then r = (( 1 , v 1 ), . . . , ( P , v P )) is a run of the funnel timed transition system T U,F that brings x 0 to r(x 0 ) ∈ T f while avoiding the obstacle Ω.

This shows that the reachability of stop in A U,F,K implies that there exists an appropriate schedule of control law switches that safely brings the system to the target zone. Of course, the method is not complete, not all schedules can be obtained using the timed automaton A U,F,K . But if A U,F,K is precise enough, it will be possible to use automatic verification techniques for dynamic trajectory generation.

Remark 1. We could be more precise in the modelling as a timed automaton, if we could use non-deterministic clock resets [START_REF] Bouyer | Updatable timed automata[END_REF]; but we should then be careful with decidability issues. Additionally, non-deterministic resets are not implemented in Uppaal, which is why we have chosen timed automata with deterministic resets only.

Remark 2. As we show with some examples in Section 5, our timed-automata abstraction can be used for other types of objectives than just reachability with avoidance. In particular, the approach can be extended to timed games [START_REF] Asarin | Controller synthesis for timed automata[END_REF], where special uncontrollable transitions model uncertainty in the environment. In that case, the aim is not to synthesize one single run in the system, but rather a strategy that dictates how the system should be controlled, depending on how the environment evolves. It is worth knowing that winning strategies can be computed in exponential time in timed games, and that the tool Uppaal-Tiga [START_REF] Behrmann | UPPAAL-Tiga: Time for playing games! In[END_REF] computes winning strategies. In Section 5.1, we give an example of application where timed games and Uppaal-Tiga are used.

LQR funnels

In this section we consider the particular case of linear time-invariant stabilizable systems whose dynamics are described by the following equation:

ẋ = Ax + Bu, (4)
where A ∈ R d×d and B ∈ R d×k are two constant matrices, and u ∈ R k is the control input. We also consider reference trajectories that can be realized with controlled systems described by Eq. ( 4

To track ẋref , we compute u ∆ as an infinite-time linear quadratic regulator (LQR, see [START_REF] Sontag | Mathematical control theory: deterministic finite dimensional systems[END_REF]), i.e. a minimization of the cost:

J = ∞ 0 x T ∆ Qx ∆ + u T ∆ Ru ∆ dt
, where Q and R respectively are positive-semidefinite and positive-definite matrices. The solution is u ∆ = -Kx ∆ , with K = R -1 B T P and P being the unique positive-definite matrix solution of the continuous time algebraic Riccati equation:

P A + A T P -P BR -1 B T P + Q = 0.
The dynamics can be rewritten ẋ∆ = (A -BK)x ∆ = Āx ∆ , i.e.:

ẋ = ẋref + Ā(x -x ref ), (6)
and the matrix Ā is Hurwitz, i.e. all its eigenvalues have negative real parts. Additionally, V : x ∆ → x T ∆ P x ∆ is a Lyapunov function (V (0) = 0 and for all x ∆ = 0, it holds V (x ∆ ) > 0 and V (x ∆ ) < 0). The solutions of Eq. ( 6) can be written: -t0) x ∆ (t 0 ). Since Ā is Hurwitz, the term e Ā(t-t0) tends to 0 exponentially fast, and the tracking asymptotically converges towards the reference trajectory x ref (t). The Lyapunov function V can be used to define control funnels as follows. For α > 0, we let:

x(t) = x ref (t) + e Ā(t
F α (t) = {x ref (t) + x ∆ | V (x ∆ ) ≤ α} (7)
F is a control funnel defined over R: if x ∆ (t) = x(t) -x ref (t) is a solution of Eq. ( 5) such that x(t 1 ) ∈ F α (t 1 ), then for any t 2 > t 1 , since V (x ∆ ) only decreases,

V (x ∆ (t 2 )) ≤ V (x ∆ (t 1 )) ≤ α, and thus x(t 2 ) = x ref (t 2 ) + x ∆ (t 2 ) ∈ F α (t 2 ). F α (t) is a fixed d-dimensional ellipsoid centered at x ref (t)
. Without going into details, it is possible to get lower bounds on the rate of decay of V (x ∆ ), and effectively compute β > 0 such that, for any solution x ∆ (t) of Eq. ( 5):

∀t ∈ R, ∀δt ∈ R + , V (x ∆ (t + δt)) ≤ e -β.δt V (x ∆ (t)) (8)
This proves that if the system is inside the control funnel F α (t) at a given instant, then after letting time elapse for a duration of δt, the system will be inside the control funnel F αe -β.δt (t). Using the terminology of Section 2.3, this can be equivalently stated as follows: for 0 < α < α, the control funnel F α (t) 1 β log( α α ) -absorbs the control funnel F α (t). Thanks to this property, for a given LQR controller and a reference trajectory x ref (t), we can define a finite set of fixedsize control funnels F α0 (t), F α1 (t), . . . , F αq (t), with α 0 > α 1 > • • • > α q > 0, and absorbing transitions between them in the corresponding timed automaton.

In the remainder of the article, we will only use this kind of fixed size control funnels, which we call "LQR funnels". They are convenient because the larger ones can be used to "catch" other control funnels, and the smaller ones can easily be caught by other control funnels. Figure 3 depicts a typical sequence, where first a large control funnel (in green) catches the system, then after some time, an absorbing transition can be triggered, and finally, a new transition brings the system to a larger control funnel (in blue) on another trajectory. Besides that, testing for inclusion between fixed-size ellipsoids is easy, and therefore LQR funnels allow relatively efficient algorithms for the computation of the tuples needed for the timed-automaton reduction ((switch, [α, β], (i, j), γ, (k, l)), (invariant, S i,j , (i, j)), . . . , see Section 3).

It should be noted that the concepts of fixed size control funnels and absorbing transitions, introduced here for linear systems, are also suitable for general nonlinear systems. Lyapunov functions in general, and quadratic ones in particular, can be computed via optimization, for example with Sum-of-Squares techniques as shown in [START_REF] Majumdar | Control design along trajectories with sums of squares programming[END_REF]. By imposing specific constraints on the optimization, fixed size control funnels with exponential convergence can be obtained inducing the same kind of absorbing transitions as introduced in the last paragraph.

Examples of application

Synchronization of sine waves

In this example, there is a unique reference trajectory: x ref (t) = sin( 2π τ t), for t ∈ [0, τ ] and τ ∈ Q, and a unique LQR controller tracking this trajectory. We define two fixed size LQR funnels F 1 (the large one) and F 2 (the small one) defined over [0, τ ] such that F 2 γ-absorbs F 1 for some γ ∈ Q. The size of F 1 is computed so that an upper bound on the acceleration is always ensured, as long as the state of the system remains inside the control funnel.

The set F 1 (τ /2) contains the smaller control funnel F 2 (t) for a range of time values [α, β] for some α < τ 2 ∈ Q and β > τ 2 ∈ Q. This allows switching transitions from F 2 to F 1 with abrupt modifications of the controller clock c t . Together with the absorbing transition and "cyclic transitions" that come from the equalities F 1 (0) = F 1 (τ ) and F 2 (0) = F 2 (τ ), it results in an abstraction by the timed automaton shown on the left side of Fig. 4. The goal is to synchronize the controlled signal to a fixed signal sin( 2π τ t + ϕ 0 ). The phase ϕ 0 is initially Fig. 4. On the left: the timed automaton for the controlled signal (the system). On the right: the timed automaton used to model the target signal with an initially unknown phase ϕ0. The opponent transition (dashed) is the one used to set ϕ0. unknown, which we model using an adversary: we use a new clock c t , and an opponent transition as in the timed automaton on the right side of Fig. 4.

With these two timed automata, we can use the tool Uppaal-Tiga to synthesize a controller that reacts to the choice of the adversary, and performs adequate switching transitions until c t = c t . It is even possible to generate a strategy that guarantees that the synchronization can always be performed in a bounded amount of time. We show in Fig. 5 a trajectory generated by the synthesized reactive controller. In this example, the phase chosen by the adversary is such that it is best to accelerate the controlled signal. Therefore, the controller uses twice the switching transition from F 2 to F 1 with a reset of the controller clock from α to τ /2 ( 1 and 2 in Fig. 5). Between these switching transitions, an absorbing transition is taken to go back to the control funnel F 2 ( A in Fig. 5). After the first two switching transitions, the remaining gap ε between c t and c t is smaller than τ 2 -α, and therefore the controller waits a bit longer (until τ 2 -ε) to perform the switching transition that exactly synchronizes the two signals ( 3 in Fig. 5).

This example shows that our abstraction can be used for reactive controller synthesis via timed games. The main advantage of our approach over methods based on full discretization is that, since a continuous notion of time is kept in our abstraction, the reactive strategy is theoretically able to exactly synchronize the controlled signal to any real value of ϕ 0 . One of our hopes is that extensions of this result can lead to a general formal approach for signal processing.

A 1D pick-and-place problem

In this second example, we show that timed-automata abstractions based on control funnels can be used to perform non-trivial planning. We propose a one-dimensional pick-and-place scenario. The set-up consists of a linear system controlled in acceleration moving along a straight line. On this line, four positions are defined as lanes (see Fig. 6). On three of these lanes (1, 2 and 3), packages arrive that have to be caught at the right time by the system and later delivered to lane 0. The system has limited acceleration and velocity, and can carry at most two packages at a time.

The LQR funnels in this example are constructed based on 12 reference trajectories. The first four have different constant positive velocities (x i ref with i ∈ {1, . . . , 4}, the fastest one being x 4 ref , and the slowest one x 1 ref ). The next four are the same trajectories with negative velocities. On each of these reference trajectories, five different control funnels of constant size are defined (F j i for j ∈ {0, . . . , 4}, the largest one being F 0 i ). The control funnels with negative constant velocity are the mirror image of those with positive velocity. Additionally, four stationary trajectories x Lk ref (with k ∈ {0, . . . , 3}) at the positions of the lanes are defined. The controllers associated to these trajectories simply stabilize the system at lane positions. For each of these trajectories a small (j = 0) and a large (j = 1) control funnel are constructed. 6). This allows the system to reach a higher or lower velocity without the need of an explicitly defined Fig. 7. In order to define the tuples (switch, [α, β], (i, j), γ, (k, l)) (see Section 3), N regularly spaced points are chosen in x k ref (defining the ellipses F l k (tn) for n ∈ {1, . . . , N }), and for each n, we set γ = tn, and if a point

x i ref (t) such that F j i (t) ⊂ F l k (γ) is found, an in- cremental search is performed to define a range [α, β] such that ∀t ∈ [α, β], F j i (t) ⊂ F l k (γ).
acceleration trajectory. While the abstraction based on these control funnels does not represent all the possible behaviors of the system (it is not complete), switching between different velocity references allows the system to perform a great variety of trajectories with continuous and bounded velocity and bounded acceleration.

To fully specify the timed-automata abstraction, the tuples defining the transition guards must be computed (see Section 3). Here, the regions of invariance defining the funnels are identically-shaped ellipses (only translated along a reference trajectory and scaled), thus the test for inclusion is computationally very cheap. Therefore, many points can be tested for inclusion on each trajectory, as depicted in Fig. 7, which leads to precise ranges for the switching transitions. Since the funnels are fixed sets translated along reference trajectories, knowing velocity or acceleration bounds on these references, and using offsets in the inclusion tests, we can ensure inclusion on the whole range of a switching transition with only a finite number of inclusion tests.

We consider an example where three packages respectively arrive on lanes 3, 2 and 1 at times t 1 arrive = 40, t 2 arrive = 111 and t 3 arrive = 122 (corresponding equality tests on c g can be used to refer to these moments in the timed automaton abstraction). The goal is to find a trajectory that catches all the packages and delivers them to lane 0. At the moment of the catch (c g = t p arrive ), the reference x i ref tracked by the system must be exactly at the correct position (i.e. on the lane of the arriving package). Depending on the reference trajectory, this corresponds to a particular value of c t . We add the following constraints on the catches: an upper bound on velocity such that the system cannot be tracking

x 4 ref , x 3 ref , x -3 ref or x -4
ref when it catches a package, and a bound on uncertainty such that the system must be in a small control funnel to catch a package. Using additional constructions in our timed automaton abstraction (for example a bounded counter that keeps track of the number of packages being carried by the system), it is easy to specify these constraints and the objective as a reachability specification that can be checked by Uppaal. Uppaal outputs a timed word that corresponds to the schedule of control law switches and the trajectory shown on Fig. 8, which successfully catches the packages and delivers them to lane 0. The two upper graphs of Fig. 8 show the evolution of the system in its state space and some of the regions of invariance when taking a switching transition (colored ellipses). The green dots mark positions at which absorbing transitions take place (F j i → F j+1 i ). Purple crosses represent a package. The lower graph compares the evolution of the position of the real system with the reference. One can see that even though the reference velocity can only take seven different values, a relatively smooth trajectory is realized. Before catching the first package, the system switches from F 4 1 to F 0 L3 1 . It then converges to F 1 L3 2 just before the catch. The difference between the real system position and the reference is very small at that point in time. The system then switches to F 0 -1 3 in order to return to lane 0. It is interesting to notice that the system chooses to return to lane 0 after having picked only one package, therefore adopting a non-greedy strategy. This is because it wouldn't have time to perform a delivery to lane 0 between the arrival of the second and third packages.

When the second package arrives on lane 2, the system catches it while being in F 4 -1 4 . This is again a non-trivial behavior: in order to get both the second and the third packages, the system has to first go a little bit further than lane 2 so as to be able to catch the two packages without violating the limit on acceleration. A slight adjustment of the reference position 5 has to be done to catch the third package exactly on time 6 . After that, the system performs a local acceleration 7 to reach lane 0 as soon as possible, and delivers the two packages.

Conclusion and future work

We have presented a timed-automata abstraction of switched dynamical systems based on control funnels, i.e. time-varying regions of invariance. Applying verification tools (such as Uppaal) on this abstraction, one can solve motion planning or more complicated problems with timing requirements. In the example of Section 5.2, we are able to generate a non-trivial solution for a pick-and-place problem. Synthesis of controllers that react to the environment can be done by solving timed games, and in the example of Section 5.1 we use Uppaal-Tiga to generate a controller that can reactively adjust the phase of a signal controlled in acceleration.

To go further and improve our abstraction, as mentioned in Section 3 (Remark 1), we could use non-deterministic clock updates and study the related decidability issues. We could also exploit the specific structure of the timed automata used in our abstraction and design dedicated verification and synthesis algorithms. Indeed, the timed automata of our model have three clocks, and there is non-determinism for only one of them (c t ). This makes us believe that we could potentially outperform the general algorithms of Uppaal and Uppaal-Tiga and solve more complex problems. Finally, in this quest to scale our approach up to larger models and more advanced specifications, we also plan to combine it to numerical and optimization methods.

Fig. 3 .

 3 Fig. 3. An absorbing transition (in green) between two switching transitions.

Fig. 5 .

 5 Fig.5. The reactive controller performs three switching transitions to exactly adjust its phase to that of the target signal.

Fig. 6 .

 6 Fig. 6. The figure on the left shows the set-up. The black dots correspond to the position of the lanes. On the right are shown some LQR funnels along the constant velocity reference trajectories in the state space.

Fig. 8 .

 8 Fig. 8. Execution of a succeeding control strategy given as a timed word.

  ), i.e. trajectories x ref (t) for which there exists u ref (x, t) such that ẋref = Ax ref + Bu ref . We can combine this equation with (4) and get ẋẋref = A(x -x ref ) + B(u -u ref ), which rewrites ẋ∆ = Ax ∆ + Bu ∆ .

This work has been partly supported by ERC Starting grant EQualIS (FP7-308087) and by European FET project Cassting (FP7-601148)

In this paper, we mostly consider state spaces that describe the position and velocity of systems controlled in acceleration. The continuity of trajectories in the state space ensures that the position is always a continuously differentiable function of time.