Thi Phuong

Khanh T P Nguyen

Mitra Fouladirad
email: mitra.fouladirad@utt.fr

Antoine Grall
email: antoine.grall@utt.fr

New methodology for improving the inspection policies for degradation model selection according to prognostic measures

Keywords: Prognostic prediction, RUL prediction, Non-periodic inspection, Model selection, Degradation process

New methodology for improving the inspection policies for degradation model selection according to prognostic measures

I. INTRODUCTION

Nowadays, the implementation of Prognostic and Health Management (PHM) solutions leads to improve the competitiveness of companies by increasing production availability, avoiding contract penalties and reducing maintenance costs. Indeed, PHM is concerned with estimating component health and predicting Remaining Useful Life (RUL) of a given system. Prognostic methods can be classified into three main approaches: model-based prognostics, data-driven prognostics and hybrid. Data-driven approaches uses very often machine learning algorithms to transform monitoring data before failure into compatible models [START_REF] Liu | A novel dynamic-weighted probabilistic support vector regression-based ensemble for prognostics of time series data[END_REF], [START_REF] Torabi | A Survey on Artificial Intelligence-Based Modeling Techniques for High Speed Milling Processes[END_REF]. However, it is difficult to acquire such large amount of training data to built suitable machine learning algorithms. Model-based approaches allow us to overcome this drawback by taking into account the physical degradation phenomenon of the system and deriving appropriate mathematical models [START_REF] Daigle | Model-Based Prognostics With Concurrent Damage Progression Processes[END_REF], [START_REF] Wang | An additive wiener process-based prognostic model for hybrid deteriorating systems[END_REF], [START_REF] Zheng | A nonlinear prognostic model for degrading systems with three-source variability[END_REF]. Furthermore, prior knowledge on the degradation phenomenon helps us to limit the list of candidate models. It can improve PHM procedure. A hybrid approach allows benefiting from the strengths of both categories [5], [START_REF] Cheng | A fusion prognostics method for remaining useful life prediction of electronic products[END_REF]. It can be based on a study of physical degradation phenomena to propose a list of models and an estimation of their parameters based on the acquired data. In this framework, one important issues is to select in the competing candidate list the model which is the best to describe the underlying degradation phenomenon in the presence of deteriorating data [START_REF] Ye | Stochastic modelling and analysis of degradation for highly reliable products[END_REF], [START_REF] Zhou | Degradation modeling applied to residual lifetime prediction using functional data analysis[END_REF].

For model selection, numerous papers in literature propose decision criteria that depend mainly on the specific purpose for which the model can be used, see for instance [START_REF] Bousquet | Bayesian gamma processes for optimizing conditionbased maintenance under uncertainty[END_REF], [START_REF] Freitas | Reliability assessment using degradation models: bayesian and classical approaches[END_REF], [START_REF] Zhou | Degradation modeling applied to residual lifetime prediction using functional data analysis[END_REF]. Authors in [START_REF] Saxena | Metrics for Offline Evaluation of Prognostic Performance[END_REF] introduce a set of metrics to assess the prognostic ability of a given model when degradation data and true failure times are known. [START_REF] Baussaron | Reliability assessment based on degradation measurements: How to compare some models?[END_REF] proposes a method to compare degradation models based on the quality of reliability assessment method. [START_REF] Hu | Online performance assessment method for a model-based prognostic approach[END_REF] develops an approach for online assessing the performance of a prognostic method in situations of very poor knowledge on the degradation process. To our knowledge, no existing research address the issue of the relevant number of monitoring data for model selection. Continuous monitoring data provide a large amount of information and therefore can help improving the model selection accuracy. However, continuous monitoring can be very costly and is sometimes impossible to implement in practice. Limits due to the monitoring technology or to human resources can lead to some constraints on the time between two successive inspections. Among different inspection policies, the policy with more frequent inspections can capture more information on the degradation evolution. Nevertheless, it is also very costly. This paper aims to address the gap in the literature by investigating how to find a balance between the monitoring cost and the necessary information level for degradation modeling. The competing models are assumed to follow the Wiener and Gamma processes that have been largely addressed in the framework of degradation modeling [START_REF] Gorjian | A review on degradation models in reliability analysis[END_REF], [START_REF] Tamura | Reliability analysis based on a jump diffusion model with two wiener processes for cloud computing with big data[END_REF], [START_REF] Van Noortwijk | A survey of the application of gamma processes in maintenance[END_REF]. In fact, the Wiener process is very popular in engineering and deterioration modeling when the health indicator evolves non-monotonically. The statistical properties of the failure time of a Wiener process are studied in [START_REF] Chhikara | The inverse Gaussian distribution as a lifetime model[END_REF], [START_REF] Cox | Analysis of survival data[END_REF]. Since the 1970s, this process has been widely studied in reliability and lifetime analysis. Authors in [START_REF] Doksum | Models for variable-stress accelerated life testing experiments based on wiener processes and the inverse gaussian distribution[END_REF] used the Wiener process with drift to model accelerated life testing data. In [START_REF] Whitmore | Estimation degradation by a wiener diffusion process subject to measurement error[END_REF], [START_REF] Whitmore | Modelling accelerated degradation data using wiener diffusion with a time scale transormation[END_REF] the impact of measurement errors on the degradation of self-regulating heating cables is analyzed. Authors in [START_REF] Padgett | Inference from accelerated degradation and failure data based on gaussian process models[END_REF], [START_REF] Park | Stochastic degradation models with several accelerating variables[END_REF], [START_REF] Whitmore | Recurrent first hitting times in wiener diffusion under several observation schemes[END_REF] also focused on the stopping time (failure time) of Wiener process and expanded the existing theoretical results in this domain. Recently, the Brownian motion with non-linear drift has attracted more attention in engineering problems and residual lifetime estimation, see for example [START_REF] Si | Remaining useful life estimationa review on the statistical data driven approaches[END_REF], [START_REF] Wang | A model for residual life prediction based on brownian motion with an adaptive drift[END_REF], [START_REF] Wang | Wiener processes with random effects for degradation data[END_REF]. When the degradation path is monotonically increasing, the Gamma process which a is a positive stochastic process with independent increments is a suitable candidate to model the degradation. Hence very often the Gamma process is used to describe the model a deterioration caused by the accumulation of wear [START_REF] Avramidis | Efficient simulation of gamma and variance-gamma processes[END_REF]. This process is a jump process which can be roughly considered as a succession of the frequent arrival of tiny increments. This rough description makes it relevant to model gradual deterioration such as corrosion, erosion, wear of structural components, concrete creep, crack growth [START_REF] Bažant | Stochastic process for extrapolating concrete creep[END_REF]. Moreover, the Gamma process has a probability distribution function which permits feasible mathematical developments. This process has been widely used in deterioration modeling for condition-based maintenance (see [START_REF] Van Noortwijk | A survey of the application of gamma processes in maintenance[END_REF]).

Considering a list of model candidates including the Wiener and Gamma processes, the paper proposes a new methodology to investigate optimal inspection policies that are less expensive but ensure information quality for model selection. Therefore, it could help companies to improve their actual monitoring policy for reducing operation and maintenance costs.

The remainder of the paper is as follows. Section II describes the problem statement. Section III presents the proposed methodology to optimize inspection policies for model selection. In Section IV, a set of models under consideration is exposed. Section V is dedicated to illustration of the performance of our methodology through numerical examples. Finally, Section VI presents the conclusion and further research work.

II. PROBLEM STATEMENT

A regular inspection policy provides a large amount of information for degradation modeling and ensures the prognostic precision of system failure. However, it can be very costly and is sometimes impossible to implement in practice. Therefore, the question is how to find a balance between the inspection policy cost and the required information level for degradation modeling. This paper aims to address this issue. A methodology is proposed to improve the inspection policy, i.e to reduce the inspection cost while ensuring the information quality.

The system state is assumed to be observed through periodic or non-periodic inspections. Each inspection instantaneously reveals the system state. For periodic inspections, the length between two successive inspections ∆ Ti is a constant. For non-periodic inspections, ∆ Ti is a function of the degradation state observed at the previous inspection.

In practice, due to technological and logistic constraints, the inspection interval is always greater than ∆ min , which is the minimum achievable time between two successive inspections. Therefore the periodic inspection policy with inspection interval equal to ∆ min is the policy which can collect the maximum information about the system degradation evolution. This latter may also be the most costly. Hereafter, this inspection policy is called the Benchmark policy and is considered as a reference to compare the efficiency of the other policies.

It is necessary to optimize inspection policies to find a balance between the monitoring cost and the minimum information required for relevant degradation modeling. A "relevant" degradation model is obtained when the model selection results are as good as for the Benchmark policy. In detail, data is firstly acquired following different inspection policies, including the periodic or non-periodic inspections. Corresponding to every set of data, a degradation model selection procedure is then performed. The results of the model selection procedure are compared to the results obtained by the Benchmark policy. If the number of inspections of a given policy is small while the selection result based on this policy is consistent with the one based on the Benchmark policy, this policy is considered as better than the Benchmark policy. The optimal inspection policy is the one having a minimal number of observations and ensuring the best consistency with the Benchmark policy.

III. METHODOLOGY TO IMPROVE THE INSPECTION POLICY FOR DEGRADATION MODEL SELECTION

We firstly define an investigation set of different types of inspection policies, see subsection III-A. Afterward, corresponding to each inspection policy, the model selection procedure is performed, see subsection III-B. Finally, in subsection III-C we propose a criteria set to optimize the inspection policies such that the selection model results based on these policies are acceptable compared to the Benchmark policy.

A. Inspection policies for acquiring degradation information

An inspection policy is considered as relevant when the number of inspections is small, but the collected data is informative enough. For instance, for a new system, it is not necessary to often inspect the degradation state. On the opposite it is preferable to observe frequently the degradation state for an old system. Therefore, the function characterizing the length between two successive inspections (∆ Ti = T i -T i-1 ) is generally a decreasing function of the observed degradation state at time T i-1 denoted by X i-1 . Consider a fixed failure threshold L. The duration ∆ Ti before the next inspection time T i is assumed to be one of the functions described hereafter. 1) Case 1: a linear function,

∆ Ti = max b -b • X i-1 a • L , ∆ min ; if a > 0 = b; if a = 0 (1)
2) Case 2: a concave function, 3) Case 3 : a convex function,

∆ Ti = max b ln(aL) ln(aL -X i-1 ), ∆ min (2)
∆ Ti = max b -b X i-1 L a , ∆ min ; (3)
Considering ∆ min is given, an inspection policy is characterized by 3 parameters (a, b, c) where:

• The parameter c characterizes the case number i.e. the shape of the decreasing function which can be linear (c = 1), concave (c = 2) or convex (c = 3). • The parameter b represents the first inspection time,

T 1 = ∆ T1 = b, ∆ min ≤ b ≤ L.
• The parameter a adjusts the decreasing rate of ∆ Ti , a ∈ [0, 1]. When the degradation state exceeds the threshold a • L, X i-1 ≥ a • L, the periodic inspection with ∆ Ti = ∆ min will be applied. Figure 1 illustrates different inspection functions. When c = 1, see Figure 1(a), ∆ Ti constantly decreases with system degradation state. When a = 0, a periodic inspection policy is considered and the period length is defined by b i.e. ∆ Ti = b. When c = 2, see Figure 1(b), ∆ Ti is slowly decreasing in early stage (when the degradation state is still insignificant) while it rapidly declines later. Contrarily, when c = 3, see Figure 1(c), ∆ Ti rapidly decreases in early stage while it slowly declines when system is more degraded.

B. Procedure of model selection

The model selection procedure has the following steps: 1) Data acquisition. The degradation data are recorded at the inspection times

T i , i = 1, 2, • • • n until it exceeds a given threshold τ • L, where τ ∈ [0, 1].
For the Benchmark inspection policy, the length between two inspections are constant an equal to ∆ Ti = ∆ min .

For other periodic or non-periodic policies, ∆ Ti is a function of the degradation state observed at the previous inspection time. The considered function ∆ Ti have been introduced in subsection III-A. 2) Parameter estimation. In reliability engineering, parameter estimation refers to the process of using sample data (degradation records at inspection times) to find the model parameters allowing the best data fitting. Among numerous parameter estimation methods, the Maximum Likelihood Estimation (MLE) is one of the most widely used in statistics [START_REF] Myung | Tutorial on maximum likelihood estimation[END_REF]. It proposes a general approach to determine the values of the model parameters through the maximization of the likelihood function. Suppose that n is the number of collected degradation records, x j is the degradation records at the jth observation time (1 ≤ j ≤ n). Let p(x j , t j |x (j-1) , t (j-1) ; θ) denote the transition probability density function of the degradation process from the state x (j-1) at time t (j-1) to the state x j at time t j knowing that θ is the vector of model parameters. The log-likelihood function for an uncensored sample dataset is given by:

L(θ) = n j=1
log(p(x j , t j |x (j-1) , t (j-1) ; θ))

The optimal parameter vector is given by θ * = arg max θ L(θ). In this paper, we propose to use Nelder-Mead (NM) method as a numerical method to maximize the likelihood or the log-likelihood function. The NM algorithm was firstly proposed by Nelder-Mead in 1965 to solve an unconstrained optimization problem. As it is only based on the evaluation of the objective function values, but no derivatives, this algorithm could be used for objective functions whose closed form derivatives are difficult to obtain. Although it generally lacks rigorous convergence properties it can reach satisfactory results after the first few iterations in practice [START_REF] Singer | Nelder-mead algorithm[END_REF]. The choice of this generic algorithm is to emphasize that the optimization problem can be solved numerically easily. 3) Evaluation of selection criteria. The candidate models are compared according to different criteria focusing on the prognostic measures. They are described hereafter:

• Prognostic accuracy criterion (P AC) is a measure that evaluates the precision of Residual Useful Lifetime (RUL) estimation corresponding to a given amount of observation data within a fixed observation period. Let L be the failure threshold of the system and t F the failure time i.e. the first time the degradation level crosses L. Let t τ be the first time the degradation level exceeds the threshold

L τ = τ • L with τ ∈ [0, 1].
The prognostic accuracy criterion (P AC) is defined in this paper as a prognostic measure for small values of . It evaluates the probability that the RUL of component calculated at time t τ and denoted by RU L tτ remains between the -bounds:

+ = (1 + )t F -t τ and -= (1 -)t F -t τ .
The P AC is defined as follows:

P AC = P[ -≤ RU L tτ ≤ + ] (5)
The choice of , ζ, and τ depend on the operation requirements. The approach to evaluate the predicted RUL distribution is presented later in section IV. Among candidate models, the one having the maximal PAC will be selected. • Hybrid Criterion (HyC) is a measure that allows to take into account the goodness-of-fit information of observed data in the evaluation of the RUL estimation accuracy. The HyC is calculated as the sum between the logarithm of PAC indicator based on RU L tτ and the mean log-likelihood of the observed degradation increments without exceeding the threshold L τ :

HyC = -log(P[ -≤ RU L tτ ≤ + ) - 1 n log p(x 1 , t 1 |X 0 , 0; θ) + n|(xn≤τ •L) i=2 log p(x i , t i |x i-1 , t i-1 ; θ) ( 6
)
where n is the number of observation points. For the HyC criterion, the model with the minimal value of HyC will be chosen 4) Normalization of selection criteria. According to each selection criterion, the best model will be chosen. For model selection, we are interested in the ranking of candidate models according to the identified selection criterion more than its measured values. Moreover, in order to be able to compare the consistency of ranking results of different selection criteria, normalized values need to be considered. The latter allows us to adjust values measured on different scales to a common scale, for example between 0 and 1. Let CR i be the measured values corresponding to the Hybrid Criterion of the ith model in a set of M considered models, CR i ∈ R + . The normalized values is given by:

CR i = min i∈M (CR i ) CR i (7)
For the Prognostic Accuracy criteria, let CR i be the value of the P AC indicator for the ith model in a set of M candidates. The normalized value is defined as:

CR i = CR i max i∈M (CR i ) (8)
After normalization, according to all criteria, the value of

CR i falls into [0, 1]. The model i with CR i = 1 is the best one.

C. Optimization criteria of inspection policies for model selection

For model selection, the optimal inspection policy is the one having a small number of observations but still offers sufficient information for the model selection procedure. In detail, considering Benchmark monitoring as a reference policy, an optimal inspection policy satisfies the following conditions:

1) It has the minimal ratio between its inspection numbers and the inspection numbers of the Benchmark policy. Let n j and n s be respectively the observation numbers of the jth inspection policy and the Benchmark monitoring policy, the first optimality criterion is defined by:

OC1 j = n j n s (9)
For an identified degradation process, the number of observations for each inspection policy is always smaller that the one of the Benchmark policy. Therefore, the ratio between them falls into [0, 1]. Its value is the proportion of the sample size obtained from the jth inspection policy by comparison with the one of the Benchmark policy.

2) It has the minimal inconsistency score obtained from its resulting model selection when compared with the standard result based on information acquired from the Benchmark policy. The model selection result are characterized by the vector CR in space R M where M is the number of considered models. The M components of CR are evaluated by Eq.( 7) or ( 8). Let CR s and CR j be respectively the standard model-selection-result obtained from Benchmark policy and the model-selection-result based on jth inspection policy. The inconsistency score between them is given by:

OC2 j = 1 M CR s -CR j p (10)
where . p is a given norm. Note that we have to divide the "distance" between CR s and CR j by M in order to avoid the impact of the size of model candidate set, M .

IV. PRESENTATION OF CANDIDATE MODELS

In the framework of this paper, the set of candidate models is based on Brownian and Gamma family, which has widely been considered for degradation modeling, [START_REF] Tamura | Reliability analysis based on a jump diffusion model with two wiener processes for cloud computing with big data[END_REF], [START_REF] Van Noortwijk | A survey of the application of gamma processes in maintenance[END_REF].

A. Wiener process and its extensions

To consider a general degradation modeling framework and take into account the possible existing physical models, it is possible to introduce stochastic differential equations (SDE) based on a standard Brownian motions B t :

dX t = m(X t , t)dt + σ(X t , t)dB t ,
(m, σ) : R × R + → R are respectively the drift and the diffusion coefficient. These equations appear at the beginning of 20th century in statistical mechanics and have been thoroughly formulated by Itô [START_REF] Itô | Stochastic processes[END_REF], [START_REF] Tseng | Stochastic diffusion modeling of degradation data[END_REF]. Such equations can be considered as derived from existing physical models by adding Gaussian "white noises" on measurements. They permit a wide range of degradation modeling due to the flexibility of the structure and functional parameters.

In this section, some specific Generalized Wiener processes and extensions derived from SDE are presented. For each case, the differential equation, the related distribution functions, and some stochastic properties are exposed.

1) Diffusion processes:

A Diffusion process has the following properties. X t is solution of the SDE dX t = µ(t)dt + σdB t , where µ(t) is a function of t and B t is a standard Brownian motion. The transition probability to X t = x knowing that X s = y is given by:

p(x, t|y, s) = 1 4πσ(t -s) exp - (x + M (t, s) -y) 2 4σ 2 (t -s) , (11)
where M (t, s), γ(t, s) are given by:

M (t, s) = - t s µ(u)du, γ(t, s) = t s σ 2 (u) 2 du (12)
The mean and variance values of X t are given by:

E[X t ] = -M (t, 0), Var[X t ] = σ 2 t (13)
a) M 1 : Wiener process with linear drift: This process is the special case of a Diffusion process when the drift and the variance are not time dependent (µ, σ are constant). This Wiener process (or drifted Brownian motion) which is also a Lévy process is suitable for fluctuating degradation records linearly increasing in time. It will be referred as M 1 in Section V.

The RUL cumulative distribution function (cdf) for a drifted Brownian motion given the observation value X t = x at the observation time t are given as follows [START_REF] Baussaron | Reliability assessment based on degradation measurements: How to compare some models?[END_REF]:

F RU L (x,t) (u) = Φ -L + µu + x σ √ u +e 2 µ σ 2 (L-x) Φ -L -µu + x σ √ u (14)
where Φ(•) is the standard normal cumulative distribution function. b) M 2 : Wiener process with time-dependent drift: This process is the particular case of a diffusion process when the degradation process is exponentially increasing in time (µ(t) = at b ). It will be referred as M 2 in Section V. In this case, the ratio between its drift and diffusion is not a constant and also depends on the time. Therefore, it is difficult to derive the explicit expression of the RUL distribution. Its evaluation requires solving a non-singular Volterra Integral Equation. It can be done numerically, see e.g. [START_REF] Deng | Residual Useful Life Estimation Based on a Time-Dependent Ornstein-Uhlenbeck Process[END_REF].

2) Diffusion process with purely time-dependent drift and diffusion (M 3 ): This process is the particular case of a diffusion process when µ(t) and σ(t) are time dependent functions independent of X t ). This is suitable for a degradation process including random walks with time-dependent drift and diffusion terms. In this paper, we consider a special case of the purely time dependent drift and diffusion Brownian motion:

µ(t) = cat b and σ(t) = √ 2at b .
It will be denoted as M 3 in Section V. As the power-law drift is proportional with drift, according to [START_REF] Molini | First passage time statistics of Brownian motion with purely time dependent drift and diffusion[END_REF], the RUL CDF of the process at time t given a degradation level at time t, X t = x, is derived:

F RU L (x,t) (u) = Φ -L -cγ(u + t, t) + x 2γ(u + t, t) +e c(L-x) Φ -L -cγ(u + t, t) + x 2γ(u + t, t) (15)
where γ(u + t, t) is given by Eq.( 12).

B. Gamma process

The increments of a gamma process (X t ), are independent and such that for s < t X t -X s ∼ Ga(α(t) -α(s), β) with the transition probability density function from X s = y to X t = x as follows:

p(x, t|y, s) = β α(t)-α(s) Γ(α(t) -α(s)) (x -y) α(t)-α(s)-1 e -(x-y)β
(16) where the shape function α(t) is an increasing function defined on R + . Γ is the Euler's Gamma function. and β the scale parameter. The mean and variance values of X t are given by:

E[X t ] = α(t) β , Var[X t ] = α(t) β 2 (17)
The choice of α(.) and β allows to model various deterioration behaviors from almost deterministic to very chaotic. Note that, based on the form of α(t), the Gamma process can be:

• Homogeneous if α(t) is a linear function in t: α(t) = at. This process is denoted M 4 . • Non-homogeneous if α(t) is a non-linear function: for instance α(t) = at b , a > 0, b > 1 for the process referred as M 5 in the following. Given a degradation level X t = y at time t, the cumulative distribution of the remaining useful life of a gamma process is given by [START_REF] Paroissin | Failure Time of Non Homogeneous Gamma Process[END_REF]: These values are considered as the reference values and used to assess the efficiency of other inspection policies characterized by different decision variables (a, b, c). Note that, for P AC criterion, we assess the probability that the estimated failure time belongs to the interval [ -, + ] where = 0.01.

F RU L (y,t) (u) = Γ(α(u + t) -α(t), (L -y t )/β) Γ(α(u + t) -α(t)) ; Γ(•,
To eliminate the biased results induced by the random feature of degradation processes, we repeat the above model A. Discussion on the convergence for optimal results with respect to the number of simulation processes In order to ensure the robustness of our results, we study the convergence for optimal results (objective function) with respect to the number of simulated processes N . However, since our problem is a bi-objective optimization problem, the solution is not given by a single point but rather by a set of efficient solutions, called the Pareto set. Therefore, the issue is to measure the difference between the Pareto fronts or to assess its convergence when the simulation number, N , is increasing. The following approach is proposed to address this problem. In order to measure the distance between two sets having different size, we propose to use the Hausdorff distance d H which is well-known and commonly used in the field of evolutionary multi-objective optimization (EMO). It defines a metric in the mathematical sense on the set of compact subsets of the R n , therefore returns one value characterized information about the relation between two sets having different sizes, refer to [START_REF] Schutze | Using the Averaged Hausdorff Distance as a Performance Measure in Evolutionary Multiobjective Optimization[END_REF].

In this purpose, first let use define GD(A N ) and IGD(A N ) respectively the Generational Distance and the Inverted Generational Distance from the candidate set A N to the reference set P Nmax as follows:

GD(P Nmax ) = 1 n A N n A N i=1 d p i 1/p IGD(A N ) = 1 n Nmax   n Nmax j=1 dp j   1/p (19)
where d p i denotes the minimal p-norm distance from a i to P Nmax , i.e. d p i = min yj ∈P Nmax ||a i -y j || p-norm and dp j denotes the minimal p-norm distance from y j to A N , i.e. dp j = min ai∈A N ||a i -y j || p-norm . The Hausdorff distance is evaluated as follows:

d H (A N , P Nmax ) = max(GD(P Nmax ), IGD(A N )) (20)
After evaluating the Hausdorff distance (d H ) from a Pareto front A N corresponding to N numbers of simulations to the reference Pareto front P Nmax , we find that this Hausdorff distance could be considered as a convex decreasing function in N . In fact, it is slightly constant when N > 80, see Figure 2 for illustration. More over its standard deviation (σ) is also decreasing when simulation number N increases. For degradation data generated by different models M 1 , M 2 , ..., M 5 with coefficient of variation vc = 30%, Figure 2 illustrate the variation of d H with respect to N . Based on Figure 2 and other simulations carried out and not exposed in the paper, we consider that d H will converge when N is sufficiently large. Therefore, for numerical experiments, based on obtained results, N is set to 100 in order to ensure the robustness of our results. 10)) between selection result based on the Benchmark inspection policy and the one based on other inspection policies. The vertical axis represents the proportion between inspection numbers (Eq.( 9)) of such inspection policy when compared with the observation points of the Benchmark policy. The coordinate of every black point in subfigures characterizes the bi-dimensional objective function's value for every considered inspection policy. Each point is corresponding to an inspection policy characterized by the set of parameters (a, b, c). Note that all types of linear, concave, convex inspection policies are investigated at the same time. The Pareto front, characterized by red circles, represents the problem trade-offs of an objective function, that can not be improved in one dimension without being worsened in another. We find that the Pareto fronts of the optimal inspection policies have convex form, see Figure 3 for example. When the inconsistency score increases from 0 to 1, the proportion of observation numbers also decreases from 1 to 0. In other words, reducing the inconsistency score will require more inspection numbers.

B. Pareto front of optimization of the condition based inspection policies for model selection

For more details, we consider Figures 4 and5 which illustrate the optimal results for data generated by different models Mi, i E {1, • • • , } and different variance coefficients considering the two decision criterions.

For data generated by diffusion type degradation processes, we find that:

• According to PAC criterion, it is difficult to ensure an acceptable inconsistency score for the model selection results based on RUL estimation accuracy when reducing the number of observations. Considering Figure 4(a) for example, it requires at least 60% of observation points of the Benchmark policy to guarantee that the inconsistency score of the result when comparing with the Benchmark policy one is less than 0.2.

• In the H yC criterion, taking into account the goodnessof-fit information of observation data in the evaluation of the RUL estimation accuracy allow us to guarantees the robustness of mode l selection results when reducing the number of observation, especially for a high coefficient variance. Indeed, considering Figure 4(b) for example when vc = 0.5 it requires less than 10% of observation points of the Benchmark policy to ensure that the inconsistency score of the result when comparing with the Benchmark policy one is less than 0.2. However, for data generated by Gamma degradation processes (see Figure 5), we find that the H yC do not allow us to guarantee the robustness of model selection results when reducing the number of observation, especially for a high coefficient variance. Indeed, it requires more observation numbers than the PAC criterion whe n vc is great. In detail, for data generated by Gamma homogeneous process (M 4) •,. . with high variance coefficient (vc = 0.5), it requires less than 40% of observation points according to the PAC criterion for the threshold of 0.2 for the inconsistency score. However, according to the H yC criterion, more than 60% of inspection numbers is necessary to satisfy this threshold

d •, • •• ~ . ..... ••• ••••••••~••••••••••••••••• • ' 0,
•, o d ' •. •, •••, Q. i : : : ..__• •_ • ._ : i- :: • .c.... • a... • • • a... • • ;c. • • =-c • • =-c • • =-c • • ;-c • • =-c •• ... 0.

C. Discussion on the optimal inspection policies

In the previous section, for each configuration, the optimal inspection policy with regards to number of inspections and closeness to the best choice is proposed. Nevertheless, the information on the type of this policy is not discussed. Indeed, it is not discussed that in each configuration, the proposed optimal policy corresponds to which percentage of optimal policies obtained in the overall N simulations. It is possible that the number of times that this policy is considered to be the best is very close to the number of times a second candidate is optimal. Moreover, it is possible that the second candidate can be implemented much easier than the first candidate without enormous loss of optimality. To explore these matters, in this section, the percentage of inspection policy types are given and analyzed.

In this section, we consider under every selection criteria, which inspection policies will be preferred for a given data set. In this purpose, numerous numerical examples are carried out: model selection is performed for data generated by 5 models M i , i ∈ {1, 2, ..., 5} with different values of vc = [0.1, 03, 0.5]. For this numerical example, we consider 50 couples of the inconsistency score threshold (IST ) and the proportion of observation points threshold (P OT ) where IST ∈ {0.1, 0.2, ..., 0.5} and P OT ∈ {0.1, 0.2, ..., 1}. Thus, we obtain 50 sets of inspection policies. Each policy in these sets ensures that the corresponding inconsistency score and proportion of observation points are inferior to the coupling threshold of IST and P OT .

To illustrate these conclusions, we illustrate in Figure 6 the number of chosen inspection policies and the proportion of their type (Periodic, Linear, Concave, Convex) corresponding to 50 cases of couples (IST, P OT ) when data is generated by M 1 with vc = 0.3. For the first 10 cases (Case 1 to Case 10), P OT is respectively increasing from 0.1 to 1 with IST = 0.1. Next, for Case 11 to Case 20, P OT is respectively increasing from 0.1 to 1 with IST = 0.2 and similar for next cases. We find that according to P AC, among chosen inspection policies, the linear and concave form are the two best types (see sub-figure 6.a) while according to HyC, the periodic inspections and concave form are preferred (see sub-figure 6.b). For other cases, the characteristics of the chosen policies are then summarized in Table I.

The results in Table I are summed up as follows.

• In general, according to each model selection criteria, each data set admits its optimal inspection policies. For example, according to P AC, the optimal non-periodic linear and concave inspection policies are approved for data generated by a Brownian process while the optimal periodic and non-periodic concave polices are preferred for the Gamma non-homogeneous model. On the other hand, according to HyC, the optimal concave and linear policies are favored for the Gamma model family while the periodic inspections are chosen for Wiener model family. Therefore, when we do not have the prior information about data type, periodic inspections are recommended.

• When IST and P OT are small, the criteria lead to totally different results:

according to the P AC criterion, non-periodic inspections are generally preferred for degradation 8.

• Un

T 1 ↓ Periodic ↑ T 1 ↑↓ Periodic ↑ T 1 ↑↓ Periodic ↑ T 1 ↓ Periodic ↓↑ T 1 ↑ Periodic ↓↑
small, periodic inspections are generally preferred for degradation modeling of data generated by a Wiener process (see sub-figure 6.b for an example) while non-periodic inspections are favored for data generated by Gamma processes (see sub-figure 7.(b) for an example).

• When the variance coefficient vc is increasing, the first inspection time T 1 arrives very early for data generated by a homogeneous Gamma process or a Wiener process, while no distinguished features can be derived for T 1 for data generated by other processes. for Wiener or diffusion processes, periodic inspections are more often used according to HyC while non-periodic ones are more used according to P AC. according to both P AC, periodic inspections are more preferred for Gamma homogeneous process while for the non-homogeneous Gamma process the periodic inspections not always favored.

D. Discussion of quality-index-based-filter

In the previous section, we have discussed the features of the chosen optimal inspection policies according to a couple of the inconsistency score and the proportion of observation point threshold (IST, P OT ). However, there could exist cases where the prognostic metrics (before normalization) are not good enough but the inconsistency scores of the selection results are acceptable. To avoid these unusual cases and to ensure their elimination in decision making, we propose in this section a quality-index-based-filter. This quality index is also the ratio between the prognostic accuracy metrics given by an inspection policy when comparing with the Benchmark policy.

Recall that given observation data acquired by the Benchmark policy, the CR s i is the selection criteria value before normalization for the ith candidate model. When considering Hybrid criterion (HyC), the kth model having the minimal value CR s i is the best candidate: k = arg min (i∈M ) (CR s i ). When considering Prognostic Accuracy criterion (P AC), the kth model having the maximal value of CR s i is chosen: k = arg max (i∈M ) (CR s i ). For model selection, the value of CR s k could be considered as a standard value. When we use an inspection policy (for example j-th policy), we hope that its value CR j k is close to the standard value, CR s k . Based on this concept, the quality index corresponding to jth inspection policy, QI j is defined by :

QI j = CR s k CR j k ; if using HyC criterion = CR j k CR s k ; if using P AC criterion ( 21
)
QI j ∈ [0, 1]
, where 1 is the best value.

For a quality-index-based-filter, only inspection policies whose quality index is superior or equal than a threshold QIT are kept. Considering the inconsistency score threshold IST = 0.5, the minimal proportions of observation numbers according to different values of quality index threshold, QIT ∈ [0.05, 1], are represented in Figure 8 and9. We find that when QIT is increasing, the minimal required observation number is also increasing. • For data generated by Wiener or diffusion processes, taking into account the goodness-of-fit information of observation data (HyC) in the evaluation of the RUL estimation accuracy allows us to guarantee a high-quality index of prognostic metric even if the proportion of observation numbers is less than 0.2. Indeed, considering Figure 8(a) for example, when QIT = 0.8, according to P AC the minimal proportion of observation points when comparing with the Benchmark policy is more than 0.6 while it is less than 0.2 according to HyC. for both criterions P AC and HyC, the effect of high variance coefficient is not significant.

• For data generated by Gamma processes, according to both P AC and HyC, taking into account the goodness-of-fit information of observation data in the evaluation of the RUL estimation accuracy do not allow us to guarantee a highquality index of prognostic metric when reducing the inspection numbers. high variance coefficient strongly effects on the requirements of observation points to satisfy a quality index threshold. For example, see Figure 8(b), when QIT = 0.8, the proportion of observation points when comparing with the Benchmark policy is less than 0.4 at vc = 0.1 while it is equal to 1 at vc = 0.5. Next, we discuss the type of the chosen optimal inspection policies after using the quality-index-based filter. For numerical examples, QIT is chosen as : QIT = 0.8. Similar to Section V-C, we consider 50 couples of (IST, P OT ). The numerical results are given in Figures 10 and9. The results can be summed up as follows.

• For Wiener or diffusion processes, according to P AC, in order to satisfy at the same time the requirements of the quality index and the inconsistency score, it requires at least 50% of the observation numbers of the Benchmark policy. For example, consider Figure 10.a, no inspection policy and help the manager to reduce more than 50% of the observation numbers. according to HyC, when the requirement of proportion of observation points is small, P OT ≤ 0.5, the periodic inspection policies is often chosen. When P OT is high, for example in Figure 10.b when P OT ≥ 0.5, the non-periodic inspection with concave form is the best policy.

t .. •.••• J ;; •, '• j ~ \~ ~ ' ' \.\ ••:•• . • • • • • • ~ • •• i •• •. , 0.2 " o.e 0.8 '0 l g~--,-----. •-•-• ~••••••r• ~••~• 0. 2
• For Gamma process, when the variance coefficient is high, there is no inspection policy which satisties in the same time the requirements of the quality index, the inconsistency score and the proportion of observation points. Whereas, when the variance coefficient is small, according to P AC, when the requirement of observation number is low, periodic inspection policies and non-periodic inspections with concave form are favored, see Figure 9.(a) for example. according to HyC, non-periodic inspection policies (linear and concave function) is widely used, see Figure 9.(b). In order to have a better understanding of the chosen optimal policies, we consider the case study QIT = 0.8, IST = 0.2, P OT = 0.5. According to P AC criterion, there exists no inspection policy that satisfies at the same time (QIT = 0.8, IST = 0.2, P OT = 0.5). Table II represents the chosen optimal inspection policies for model selection according hybrid criterion. The values of decision variables a, b and c are respectively characterized Among the chosen optimal inspection policies, there exist:

• one periodic inspection policy (a = 0, c = 1) with the constant period length ∆ Ti = 2, • two non-periodic policies with linear inspection function.

For example, considering the policy characterized by a = 0. VI. CONCLUSION In this paper, we have presented a new methodology that allows decision-makers to improve an inspection policies for making it less expensive while ensuring the quality of information for model selection which can e.g. be used for failure prognostic. The selection procedure is based on the prognostic accuracy criterion (P AC) and the hybrid criterion (HyC) that takes into account the goodness-of-fit information of the observation data when evaluating the prognostic measure. Besides, we have investigated different types of inspection policies including periodic or non-periodic ones by considering multiples functions that are linear, concave or convex in the degradation state of the system. Among different competing policies, the optimal ones are chosen based on the multi-criteria optimization process that takes into account minimal requirements for the observation point ratio, and the inconsistency score with the quality index filter. The advantages and disadvantages of this methodology are discussed through numerous numerical examples for different types of degradation process, particularly diffusion and Gamma processes.

The paper also provides instructions to improve the selection of appropriate inspection policies according to the degradation process family as well as the evaluation criteria. In general, in case of total absence of prior information about the degradation process family (Wiener or Gamma for example), the periodic inspections could be encouraged to acquire data. We discussed, in addition, the impact of the data variance coefficient on the inspection policy optimization results. In details, high variance ratio significantly effects on the minimum requirements for observation number when considering Gamma family while its impacts are not important when examining Wiener models.

In the scope of this paper, we utilized the grid search algorithm to solve the multi-objective optimization problem. It can lead to long implementation times when investigating large search spaces for real work problems. In further works, an efficient method will be developed to reduce the optimization time. Moreover, the set of models can also be extended, and some other models such as Inverse Gaussian, Variance-Gamma or jump diffusion models can be included to capture the various characteristic of the system degradation process in practice.

Fig. 1 .

 1 Fig. 1. Illustration of different forms of inspection policy based on the degradation state when b = 5

  •) is the upper incomplete Gamma function.[START_REF] Myung | Tutorial on maximum likelihood estimation[END_REF] V. NUMERICAL EXPERIMENTS For simulation of a model selection procedure presented in subsection III-B, we generate a stochastic degradation path by one of the five models presented in Section IV (M 1: drifted Brownian motion, M 2: time-dependent-drifted Brownian motion, M 3: time-and-diffusion-dependent drifted Brownian motion, M 4: homogeneous gamma process, M 5: non-homogeneous gamma process). Its parameters are chosen in order to have the mean degradation level 100 at time t = 100 : x = 100 with different values of coefficient of variation vc ∈ [10%, 30%, 50%] at time t = 100. The evolution of a degradation process from 0 to T m = 250 is generated with the time step ∆t = 0.1. Considering the failure threshold L = 80, the degradation data is observed until threshold L τ = τ • L where τ = 0.7. Based on data inspected by the Benchmark policy, we estimate parameters for every candidate models and then evaluate their corresponding measured values of the selection criterion (CR s ) and its normalized values, CR s .

Fig. 2 .

 2 Fig. 2. Illustration of convergence of the optimal results with respect to number of simulations N

•

  We firstly consider a large number N max as the maximal simulation numbers due to the limits of simulation time. Then, we use its corresponding Pareto set as a reference set, noted P Nmax = {y 1 , ..., y n Nmax } where n Nmax is the number of points on reference Pareto front P Nmax .• According to every different value of N < N max , we obtain a candidate Pareto set, noted A N = {a 1 , ..., a n A N } where n A N is the number of points on the Pareto front A N . • The distance between reference Pareto set P Nmax and candidate Pareto set A N is evaluated and used to measure the convergency of the optimization results.

Figure 3 presents

 3 Figure3presents Pareto fronts of inspection policy optimization. In every sub-figure, the horizontal axis represents the inconsistency score (Eq.[START_REF] Deng | Residual Useful Life Estimation Based on a Time-Dependent Ornstein-Uhlenbeck Process[END_REF]) between selection result based on the Benchmark inspection policy and the one based on other inspection policies. The vertical axis represents the proportion between inspection numbers (Eq.(9)) of such inspection policy when compared with the observation points of the Benchmark policy. The coordinate of every black point in subfigures characterizes the bi-dimensional objective function's value for every considered inspection policy. Each point is corresponding to an inspection policy characterized by the set of parameters (a, b, c). Note that all types of linear, concave, convex inspection policies are investigated at the same time.

  ~--• • _ :•: • ~----...

2 OA

 2

( c )Fig. 4 .

 c4 Fig. 4. Requirement for proportion of observation points according to maximal inconsistency threshold of model selection result for diffusion type degradation data i § \ •• I : -: t i O • , • • • • § :;J ...__ • ~•• _•:_: , _, • ~•• -•• cc •• :..; • •,..: • •c.:.: • •c.:.: • •:.:.; • •c., ••=., • • f_

Fig. 6 .

 6 Fig. 6. Number of chosen inspection policies and their type proportions when data are generated by M 1-Brownian process with vc = 0.3: Periodic (Per), Linear (Lin), Concave, Convex.

Fig. 7 .

 7 Fig. 7. Number of chosen inspection policies and their type proportions when data are generated by M 4-Gamma process with vc = 0.3: Periodic (Per), Linear (Lin), Concave, Convex.

Fig. 8 .

 8 Fig. 8. Requirement for proportion of observation points according to quality index of model selection result for Wiener degradation data

Fig. 9 .

 9 Fig. 9. Requirement for proportion of observation points according to quality index of model selection result for Gamma degradation data

Fig. 10 .

 10 Fig. 10. According to QIT = 0.8, number of chosen inspection policies and their type proportions when data are generated by M 1-Brownian process with vc = 0.3: Periodic (Per), Linear (Lin), Concave, Convex.

Fig. 11 . 8 "•

 118 Fig. 11. According to QIT = 0.8, number of chosen inspection policies and their type proportions when data are generated by M 4-Gamma process with vc = 0.3: Periodic (Per), Linear (Lin), Concave, Convex.

'

  nsiste ncysc,:;f'il lnc,;-m.stEn,c;.ys,;;om (b) Data generated by J\,f(i: Gamma non-

  2, b = 41, c = 1, the first inspection time is T 1 = 41. At the inspection time, if the observed state is higher than 0.2 • 41 = 8.2, the duration until the next inspection is ∆ min . If not, it is evaluated by a linear function defined by Eq. (1).• five non-periodic policies with concave inspection function. For example, considering the policy characterized by a = 0.5, b = 25, c = 2, the first inspection time is T 1 = 25. At the inspection time, if the observed state is higher than 0.5 • 25 = 12.5, the duration until the next inspection is ∆ min . If not, it is evaluated a concave function defined by Eq. (2).

  ......•...........•...

		PAC criterion		0 i :; l . § :; l		HyC criterion
	2	0.4	0.6	0.8	i ::: ..	0 .2	0 .4	0 .6	0 .8
	(b) Data generated by M 2: Time-dependent-drifted-Brownian process
					i I : ;		HyC criterion
							0 .4	0 .6	0.8
		1 neonslstet1cy score				1noon• tency score

\

:.::::::,:,,

o d .\ .• . ' • . . i ~ . . . _ •_ • . _

. -: : -' -: : . ; ; . .; • . ; ; . .: : . . . . : :-=-• ' .:..; • •.:..; • •.;.. • •.:..; • . .,_ • .:..; • •.;.. • •.L:.LJ o.z

z Case


Clement Papadacci

T Mirault

B Dizier

M Tanter

E Messas

M Pernot

Non-invasive evaluation of aortic stiffness dependence with aortic blood pressure and internal radius by Shear Wave Elastography and Ultrafast Imaging

  L'archive ouverte pluridisciplinaire

Introduction

Arterial Stiffness is of major interest in the understanding of the cardiovascular function [START_REF] Lee | Vascular mechanics for the cardiologist[END_REF]- [START_REF] Cavalcante | Aortic Stiffness: Current Understanding and Future Directions[END_REF] and is widely recognized as a powerful prognostic factor of cardiovascular events. Many techniques have been developed to measure arterial stiffness non-invasively. However, arterial stiffness indexes are still not used routinely in clinical practice, in part because available techniques are not yet suitable for widespread clinical practice. Indeed, all these techniques rely on indirect estimates of the arterial stiffness through measurements of arterial pressure, distensibility or a combination, which at the end makes the arterial stiffness less accurate and more complex to estimate. Therefore, there is still a need for an easy to use, reliable and noninvasive method to quantify arterial stiffness.

The pulse wave velocity (PWV) is the most widely available technique [START_REF] Asmar | Assessment of arterial distensibility by automatic pulse wave velocity measurement. Validation and clinical application studies[END_REF]. PWV is related to the elasticity of the arterial wall because the propagation of the arterial pulse wave is driven by the ability of the arterial wall to deform when the pressure increases. Applanation tonometry [START_REF] Nelson | Noninvasive Measurement of Central Vascular Pressures With Arterial Tonometry: Clinical Revival of the Pulse Pressure Waveform?[END_REF] and oscillometric method have been used to detect the time delay of pressure waveforms at two arterial sites (generally at the carotid and femoral arteries). The PWV is then retrieved by dividing the time delay to the distance between the two sites [START_REF] Asmar | Assessment of arterial distensibility by automatic pulse wave velocity measurement. Validation and clinical application studies[END_REF]. This method is non-invasive, provides an average velocity over a long arterial distance, but suffers from inaccuracies due to the challenging measurement of the vascular length especially in the potential arterial tortuosity in older patients. Recently, real-time ultrasound imaging of the pulse wave propagation with a reduced number of transmitted beams has been proposed by Kanai et al [START_REF] Kanai | Propagation of spontaneously actuated pulsive vibration in human heart wall and in vivo viscoelasticity estimation[END_REF] and later, extensively studied by other groups to estimate the local PWV [START_REF] Eriksson | Arterial pulse wave velocity with tissue Doppler imaging[END_REF]- [START_REF] Luo | Pulse wave imaging of the human carotid artery: an in vivo feasibility study[END_REF]. However, although the PWV has been used successfully for decades as an arterial stiffness index, it remains an indirect measurement of the arterial stiffness which depends strongly on the size of arteries (diameter and wall thickness). For example, the PWV of arteries of different size should not be used to compare their arterial stiffness.

Local arterial stiffness can also be assessed by measuring the distensibility of the arterial wall using real-time ultrasound imaging. The technique known as echotracking has been developed to measure finely the change of diameter of the artery during the cardiac cycle [START_REF] Stadler | The application of echo-tracking methods to endothelium-dependent vasoreactivity and arterial compliance measurements[END_REF]. The arterial wall distensibility can then be computed by using the diameter and the blood pressure. As a matter of fact, this approach requires the knowledge of the intraluminal arterial pressure which cannot be measured directly. Pressure of superficial arteries is usually estimated using aplanation tonometry, but for deep arteries such as aorta, a transfer function must be applied to estimate the blood pressure. These transfer functions are established on average populations, however, there is a large scatter in measured transfer functions so that central pressure may differ considerably from the direct measurement [START_REF] Segers | Individualizing the aorto-radial pressure transfer function: feasibility of a model-based approach[END_REF], [START_REF] Hope | Use of arterial transfer functions for the derivation of aortic waveform characteristics[END_REF].

In another hand, Shear wave elastography (SWE) has been proposed to measure directly, locally, noninvasively and in real time the elasticity of soft tissues [START_REF] Bercoff | Supersonic shear imaging: a new technique for soft tissue elasticity mapping[END_REF]. Shear waves are generated by the acoustic radiation force of a focused ultrasonic beam emitted by an ultrasonic probe. By imaging the propagation of the shear wave at very high frame rate using ultrafast ultrasound imaging, the shear wave speed can be quantified locally, and the local shear modulus can be derived. This technique has been successfully implemented to many organs such as the breast [START_REF] Bercoff | Supersonic shear imaging: a new technique for soft tissue elasticity mapping[END_REF], the liver [START_REF] Muller | Quantitative viscoelasticity mapping of human liver using supersonic shear imaging: preliminary in vivo feasibility study[END_REF], the cornea [START_REF] Nguyen | Monitoring of cornea elastic properties changes during UV-A/riboflavin-induced corneal collagen cross-linking using supersonic shear wave imaging: a pilot study[END_REF], the tendon [START_REF] Brum | In vivo evaluation of the elastic anisotropy of the human Achilles tendon using shear wave dispersion analysis[END_REF] and the heart [START_REF] Couade | In vivo quantitative mapping of myocardial stiffening and transmural anisotropy during the cardiac cycle[END_REF]- [START_REF] Papadacci | Shear Wave Imaging of the heart using a cardiac phased array with coherent spatial compound[END_REF]. We previously proposed to apply this technique for the estimation of the arterial stiffness in Couade et al. [START_REF] Couade | Quantitative assessment of arterial wall biomechanical properties using shear wave imaging[END_REF]. We showed that SWE has the potential to measure quantitatively, locally and non-invasively the stiffness of the arterial wall and its variation during the cardiac cycles. However, the technique was not yet compared against other gold standard measurements and no index was proposed to be used in the clinic.

In this study, we propose to measure arterial stiffness using SWE over a large range of arterial pressure (40 mmHg-190 mmHg) and compare it with gold standard measurements. The aortic stiffness of rats (N=5) was measured non-invasively by SWE as well as the artery diameter by ultrafast imaging. Simultaneously, the blood pressure was assessed using invasive pressure catheter. Arterial blood pressure was varied using phenylephrine and the incremental Young's modulus was obtained using both invasive and non-invasive methods. Finally, to validate the qualitative mechanical behavior of the artery, shear wave speed versus diameter curves obtained non-invasively were compared to the incremental elastic modulus and the pressure elastic modulus obtained with invasive measurements.

Material and methods

1)

Stiffness evaluation in artery with SWE SWE has been extensively studied in the past years. It relies on the creation of an acoustic ultrasonic burst emitted from a focused beam of an ultrasonic probe which produces pressure radiation acoustic force. The pressure radiation force locally pushes the tissue and its relaxation produces a shear wave propagating at a speed 𝑐 in the medium. In soft tissues, it has been demonstrated that the elastic modulus (Young's modulus) 𝐸 which represents stiffness was only a function of shear modulus µ in soft isotropic tissues as: 𝐸 ≈ 3µ

(1) and shear modulus µ related to shear wave speed c and density 𝜌 of the soft tissue as: µ ≈ 𝜌 𝑐 2 (2) The shear wave speed can then be assessed by tracking the shear wave propagation with ultrafast plane wave imaging [START_REF] Bercoff | Supersonic shear imaging: a new technique for soft tissue elasticity mapping[END_REF]. While arterial wall is a soft tissue like the liver or the breast, it is rather more complicated to assess its stiffness. Firstly, because the elastic modulus is a function of stress and strain. In static soft tissues such as the breast or the liver, the stress and strain remains the same which implies a unique measurement under one strainstress condition. However, in artery the stress applied on the arterial wall varies naturally and dynamically in a cardiac cycle due to the internal blood pressure variation. The strain follows the stress variation by modifying the arterial wall thickness. It is known that the relationship between stress and strain in biological tissue is a non-linear relationship [START_REF] Fung | Biomechanics: Mechanical Properties of Living Tissues, Second Edition[END_REF]. Therefore, an assessment of elasticity at different blood pressure values (or different stress values) is necessary to assess the mechanical behavior of the artery and to characterize it under different stress-strain conditions. Secondly, because the arterial wall is an anisotropic soft tissue, the elastic modulus is a tensor with different components. The configuration shown in figure 1 (a),(b),(c), was the configuration of the shear wave elastography acquisition in this study. The stiffness along the longitudinal direction û 𝑧 (figure 1 (a)) could be different from the stiffness assessed in the circumferential direction û 𝜃 . Finally, because the arterial wall is small compared to the shear wavelength, the shear waves are guided in the arterial wall. The guided wave propagation model to retrieve the quantitative elastic modulus is complex and depends on the shear wave central frequency, the cylindrical geometry of the artery and the ratio between the wall thickness and the radius of the artery. In this study, shear wave speeds were used as a stiffness index with no attempts to derive the shear modulus.

2) Experimental Setup

A high frequency ultrasound linear probe (15 MHz, 128 elements, 0.125 mm pitch, Vermon, France) connected to an ultrasonic scanner (Aixplorer, SuperSonic Imagine, Fance) was positioned non-invasively, on the abdominal aorta of anesthetized rats (N=5). 2 pressure catheters were set invasively in the left common carotid artery and in the left femoral artery. Blood pressure at the two sites were recorded in real time (figure 2, green and black line). ECG was recorded via 4 needle electrodes implanted in the legs (figure 2, blue line). The echograph was trigged on the R-wave from ECG. Trig out of the echograph was also recorded to verify that SWE acquisitions and ultrafast acquisitions occurred at the desired time (figure 2, red line). Pressure from the two sites were delayed and averaged to estimate the blood pressure at the location of the SWE acquisition. The blood pressure was varied from about 40 mmHg to 190 mmHg using a vasoconstrictor (phenylephrine) infused via a venous catheter.

3) Data acquisition a) SWE sequence

Longitudinal shear wave speed was assessed during SWE sequence. One pattern of the SWE acquisition consisted in the emission of one pushing beam focused at the superior arterial wall (between 4 mm and 6 mm depth depending on the animal) during 300 µs, at the probe center with an F over D (focal distance over aperture) of 1.5. 80 tilted plane waves (with angles of [-2° 2°]) were transmitted at the frame rate of 31,000 images/s in order to coherently sum 40 compounded plane waves with a frame rate of 15,500 images/s according to [START_REF] Montaldo | Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography[END_REF]. This pattern was repeated 10 times every 15 ms to cover an entire cardiac cycle. This SWE sequence was trigged on the R-wave of the ECG (i.e. figure 3 a).

b) Ultrafast plane wave acquisition

Non-invasive evaluation of diameter variation was performed using ultrafast plane wave acquisition data.

After SWE acquisitions, independent ultrafast plane wave imaging acquisitions were performed for identical value of blood pressure (i.e. figure 3 b). The acquisitions were also triggered on the R-wave of the ECG. The ultrafast imaging acquisition consisted in the emission of 4500 tilted plane waves (with angles of [-2° 0 2°]) at a frame rate of 27,000 images/s in order to achieve 1500 compounded plane waves at a frame rate of 9,000 images/s. The ultrafast acquisition lasted 166 ms which corresponded to one rat cardiac cycle.

Under incompressibility assumption, wall thickness ℎ as a function of time was derived by resolving the equation:

( 𝑑 2 + ℎ) 2 - 𝑑 4 2 = 𝐶 0 (5)
With the constant 𝐶 0 :

𝐶 0 = ( 𝑑 0 2 + ℎ 0 ) 2 - 𝑑 0 4 2 (6)
Finally, the external radius 𝑏 was found as: 𝑏 = 𝑎 + ℎ (7)

5) Comparison with validated techniques

We compared the longitudinal elastic modulus values estimated by SWE to the circumferential elastic modulus assessed by two validated methods: the pressure elastic modulus [START_REF] Kamenskiy | Nonlinear Mechanical Behavior of The Human Common, External, and Internal Carotid Arteries In Vivo[END_REF] and the incremental elastic modulus [START_REF] Cox | Arterial wall mechanics and composition and the effects of smooth muscle activation[END_REF]- [START_REF] Bussy | Intrinsic Stiffness of the Carotid Arterial Wall Material in Essential Hypertensives[END_REF] derived from pressure and distensibility (relative cross-sectional diameter change for a given pressure) measurements. These two methods are based on stress and strain measurements.

a)

The pressure elastic modulus Pressure Young's modulus 𝐸 𝑃𝜃 was defined by [START_REF] Kamenskiy | Nonlinear Mechanical Behavior of The Human Common, External, and Internal Carotid Arteries In Vivo[END_REF]:

𝐸 𝑃𝜃 = 𝑃 𝑎 2 ℎ 𝛿𝑎 (8)
It was estimated by measuring invasively the blood pressure 𝑃 with the catheters, the internal radius 𝑎, the arterial radius variation 𝛿𝑎 and the wall thickness ℎ from ultrafast acquisitions data.

b) The incremental elastic modulus

This method also assumes homogeneity and isotropy of the material. The incremental elastic modulus was calculated from the pressure-radius data as:

𝐸 𝑖𝑛𝑐𝜃 = 3 2 𝑎 2 𝑏 (𝑏 2 -𝑎 2 ) ∆𝑃 ∆𝑏 (9)
Where ∆𝑃 ∆𝑏 represents the slope of the pressure-external radius curve at a specific point (i.e., pressure or strain). Incremental elastic modulus at a point was determined by calculating the derivative of the pressure-external radius curve considering the data points below and above.

Results

1)

Shear wave speed over one cardiac cycle At each SWE acquisition, shear wave propagation could be tracked in the superior arterial wall of the aorta with high temporal resolution. Figure 6 displays snapshots of the shear wave propagation superimposed on the B-mode image at different time. Shear waves appeared to propagate faster in the arterial wall than in the surrounding tissues.

SWE acquisition had the sensitivity to measure the variation of shear wave speed associated to the natural variation of arterial blood pressure during one cardiac cycle. Figure 7 a b, displays an example of shear wave speed estimations at ten different moment of one cardiac cycle of a rat as a function of time. In this example, the blood pressure varied between 92 mmHg in diastole to 130 mmHg in systole (i.e. figure 7 d). The shear wave speeds were found to follow the blood pressure variations from 5.1 m/s in diastole to 10.2 m/s in systole.

The internal radius variation was also found to follow the blood pressure variation from 811 µm to 863 µm (figure 7 c).

2) Stiffness as a function of pressure

Arterial blood pressure was gradually varied from 40 mmHg (minimum in diastole) to 190 mmHg (maximum in systole) for the N=5 rats. For each set of pressure, shear wave speeds were assessed simultaneously. Figure 8 displays the results of shear wave speed as a function of invasive measure of blood pressure for the 5 rats.

The shear wave speeds were averaged over the 5 rats over 10 mmHg step from 40 mmHg to 190 mmHg, error bars were computed as the standard deviation. Arterial stiffness was found to increase with pressure and presents two different behavior depending on the pressure range. First, a small increase of shear wave speed (slope 1.9 10 -2 𝑚/𝑠. 𝑚𝑚𝐻𝑔 -1 ) from (5.06±0.82) m/s to (5.97±0.90) m/s was found for pressure variation from (40±5) mmHg to (95±5) mmHg respectively. For higher pressure range ((105±5) mmHg-(190±5) mmHg), a large increase (slope 8.8 10 -2 𝑚/𝑠. 𝑚𝑚𝐻𝑔 -1 ) of shear speeds from (6.90±1.45) m/s to (13.85±5.51) m/s was shown. These two parts evidenced the non-linear behavior of the artery's elasticity under a radial stress. Two linear fits were applied to the two curve parts showing a large slope difference.

3) Comparison of shear wave speed against internal radius with validated techniques.

Shear wave speeds measured non-invasively with the SWE method were compared against internal radius. The results are shown in figure 9 for N=3 rats. Once again, the curve is divided in two parts. First, a small increase (slope of 18.2 𝑚/𝑠. 𝑚𝑚 -1 ) of shear wave speed (5 m/s to 5.9 m/s) was found for internal radius variation from 798 µm to 835 µm. For larger radius range (838 µm-862 µm), a large increase (slope of 174 𝑚/𝑠. 𝑚𝑚 -1 ) of shear wave speed was found (5.8 m/s to 9.9 m/s).

To compare these results, incremental elastic modulus and pressure elastic modulus estimations were performed. Measure of internal radius, external radius and invasive measurement of blood pressure were combined in order to estimate the incremental elastic moduli with respect to equation 9 and to estimate pressure elastic moduli (equation 8).

Results are shown in figure 9 for N=3 rats. The incremental and pressure elastic moduli are shown as a function of internal radius (curve blue and black respectively). Similarly, to the SWE results, two parts can be separated. The first parts of the two curves corresponded to a small increase from 166 kPa to 309 kPa for the incremental elastic moduli and from 154 kPa to 246 kPa for the pressure elastic moduli, for radii values between 798 µm and 835 µm respectively. A larger increase is then observed for internal radius variation from 798 µm to 835 µm. Interestingly, the three methods gave very close qualitative results. Although, the values of the three methods could not be compared, the mechanical behavior of the artery remained the same at small and large radius.

Discussion

In this study, SWE was performed in the aortas of N=5 rats. Shear wave propagation was observed on the anterior wall of the rats' aorta (figure 6 a b c) at ten different moments during one cardiac cycle. Shear wave speed variation due to the physiological blood pressure change was measured (figure 7 a b). Two intravenous catheters were inserted on the two edge of the abdominal aorta to monitor and record blood pressure variations in real-time (figure 7 d). UF imaging was also performed to evaluate the internal radius variation (figure 7 c) which was found to follow shear wave speed (figure 7 b) and blood pressure (figure 7 d).

Shear wave speed in the abdominal aortic wall of N=5 rats was then measured over a large range of arterial pressure by increasing the mean blood pressure using phenylephrine injection. With the invasive pressure measurements, the relationship between shear wave speed and blood pressure was assessed. An increase of shear wave speed with the pressure was observed (figure 8). This increase is due to the fact that the blood pressure level is proportional to the stress induced on the arterial wall which is also proportional to the elastic modulus. Therefore, when the blood pressure increases, the shear wave speed increases as well.

However, we observed a non-linear relationship between shear wave speed and blood pressure level indicating of a non-linear mechanical behavior. The curve can be separated in two parts. A first part can be observed at low stress between 40 mmHg and 100 mmHg as shear wave speed increases linearly and slowly with the blood pressure level with a slope of 1.9 10 -2 𝑚/𝑠. 𝑚𝑚𝐻𝑔 -1 . A second part at higher stress between 105 mmHg and 195 mmHg, still displays a linear variation but with a much larger increase of shear wave speed with blood pressure level (slope of 8.8 10 -2 𝑚/𝑠. 𝑚𝑚𝐻𝑔 -1 ). These results indicate two different mechanical behaviors of the artery which depends on the level of stress. This particular behavior could be associated with the elastin-collagen model. At low stress (BP<100 mmHg), elastin fibers are recruited to stiffen the arteries while stress is increasing [START_REF] Callaghan | Relationship between pulse-wave velocity and arterial elasticity[END_REF]. At higher stress (BP>100 mmHg), when elastin fibers are completely stretched, collagen fibers are recruited for the stiffening of the arteries [START_REF] Ong | Effect of celiprolol on prevention of cardiovascular events in vascular Ehlers-Danlos syndrome: a prospective randomised, open, blinded-endpoints trial[END_REF]. These two fibers have different mechanical properties explaining that their elasticity does not respond identically to an increase of stress. This model could explain the difference in slopes observed on the plot of the figure 8. In addition, the slope corresponding to the collagen fibers was found larger than the one corresponding to the elastin fibers which was in a good agreement with literature.

In the third part of the study, shear wave speed against arterial internal radius was assessed and compared to the incremental elastic modulus and the pressure elastic modulus (figure 9). As expected, shear wave speed which is related to stiffness, was found to increase with internal radius. Similarly to the previous part of the study, a non-linear behavior was found and two parts of the plot could be separated. For a radius of 798 µm to 835 µm, a small linear increase was found (slope of 18.2 𝑚/𝑠. 𝑚𝑚 -1 ). For a larger radius from 838 µm to 862 µm, a larger linear increase was noticed (slope of 174 𝑚/𝑠. 𝑚𝑚 -1 ). This non-linear mechanical behavior was also associated to the elastin-collagen model. Small radius corresponding to the recruitment of elastin fibers while at larger radius elastin fibers are completely stretched and collagen fibers are used. Interestingly, we found a good agreement with the two other techniques. Although, the values cannot be directly compared as the two invasive methods give an estimation of the elastic modulus, the same behavior can nevertheless be identified. A small increase was measured for the three methods for radius from 798 µm to 835 µm indicating of an elastin fibers stiffening. For the same value of radius (838 µm), the mechanical behavior of the artery changed with the recruitment of collagen fibers.

In this study, the elastic modulus was not derived with shear wave speed. However, although the proposed method does not evaluate quantitatively the elasticity via the elastic modulus, the shear velocity is an interesting stiffness index to investigate the mechanical behavior of the artery. It is also completely noninvasive compared to the incremental and pressure elastic modulii that need an invasive measurement of the blood pressure to assess the elasticity accurately. With SWE and UF sequences the proposed method can obtain the two slopes of the curve which could be used as an index to characterize the state of health of an artery. In Ehlers-Danlos genetic disease for instance, the artery has a collagen failure. As a consequence, the artery does not stiffen efficiently at high pressure resulting in the creation of aneurysms and the death of the patients. There is currently no imaging methods that can detect and diagnose these patients. The proposed method could be of a great interest to detect and characterize the Ehlers-Danlos patients by evaluating the slope of the second part of the shear wave speed/radius curve corresponding to the mechanical properties of the collagen.

In this study, the two measurements of shear wave speed and radius evaluation were performed separately over two successive cardiac cycles, but it may be possible to combine these two sequences into one and provide the stiffness-diameter relationship in real-time.

However, the proposed method still needs the measurement of the internal and external radius at a reference time to assess radius variations. This measurement performed on the focused B-mode image can be inaccurate for several reasons. First of all, the probe can be not perfectly set across the artery, inducing an underestimation of the diameter. An ultrasonic 3D ultrafast image of the artery could be used to overcome this issue [START_REF] Provost | 3D ultrafast ultrasound imaging in vivo[END_REF]- [START_REF] Apostolakis | Feasibility and Validation of 4-D Pulse Wave Imaging in Phantoms and In Vivo[END_REF]. Finally, the method could provide two index to characterize the mechanical properties in arteries of patients in the clinic.

Conclusion

In this study, we measured shear wave speed in the arterial walls of rats' abdominal aortas with shear wave elastography within entire heart cycles. Shear wave elastography was able to follow the shear wave speed variation with the natural arterial blood pressure and diameter variations. After increasing blood pressure with phenylephrine injections, the elasticity of the artery increased as a response to higher stress. Shear wave speed was measured over a large range of pressure (from 40 mmHg to 190 mmHg), and a non-linear relationship between pressure and shear wave speed was found. We also proposed a method totally noninvasive, which relies on the estimation of shear wave speed and the measurement of arterial diameter with ultrafast plane wave imaging technique. This method enabled the estimation of the arterial diameter variation with high temporal resolution through an entire cardiac cycle. The shear wave speeds and the radius were coupled to highlight the non-linear properties of the abdominal aorta and calculate the slopes associated with elastin and collagen stiffening. The proposed technique was successfully compared to two other validated techniques which relies on invasive blood pressure measurements. The slopes derived from the proposed method could be a useful index to characterize arteries completely non-invasively in clinics without the need to use blood pressure measurements.

  signals from the 80 tilted plane waves were stored in memory. Coherent compounding was performed to get 40 compounded images. B-mode images were retrieved from absolute value of compounded IQ signals (figure 4 a). Tissue velocities (figure 4 b) were obtained using a per pixel frame to frame 1D cross-correlation on demodulated compounded IQ images with an axial kernel size of 3 pixels (150 µm) to obtain images of tissue frame-to-frame axial displacements. The superior arterial wall was selected manually (figure 4 a, blue line). Tissue velocity was averaged from -1 pixel to 1 pixel in depth around the selected line and displayed as a function of time (figure 4 b). Shear wave speed was assessed by fitting the maxima of tissue velocity of the shear wave front with a linear fit (figure 4 b, white line). The slope of the fit gave access to the longitudinal speed of the shear wave. b) UF data IQ signals from the 4500 tilted plane waves were stored in memory. Coherent compounding was performed on the three angles to get 1500 compounded images. Similarly to SWE post processing, tissue velocities (figure 4 b) were obtained using a per pixel frame to frame 1D cross-correlation on demodulated compounded IQ images with an axial kernel size of 3 pixels (150 µm) to obtain images of tissue frame-to-frame axial displacements. Internal diameter 𝑑 0 and wall thickness ℎ 0 of the artery at a reference time 𝑡 0 were measured on the focused B-mode images of the SSI echograph (i.e. figure 5 a). Arterial wall radial displacement of the anterior wall 𝛿𝑑 1 and posterior wall 𝛿𝑑 2 were obtained with ultrafast sequence by integrating tissue velocities over time of the anterior wall and posterior wall respectively (i.e. figure 5 b). Internal diameter variation was therefore computed as: 𝛿𝑑 = 𝛿𝑑 1 -𝛿𝑑 2 (3) Internal diameter 𝑑 and internal radius 𝑎 as a function of time was calculated as: 𝑑 = (𝑑 0 + 𝛿𝑑) ; 𝑎 = 𝑑 2


Clement Papadacci
email: clement.papadacci@espci.fr

Camille Zordan

Mickaël Tanter

Mathieu Pernot

Mickael Tanter

Quantitative cardiac output assessment using 4D ultrafast Doppler imaging: an in vitro study

I. INTRODUCTION

Echocardiography has the advantages of being fast, low-cost, and to enable real-time measurements on patients' beating heart. It is routinely used in the clinic to evaluate the cardiac function by assessing various cardiac indexes. Anatomical indexes such as ventricular volume measurements or functional indexes such as Cardiac Output are performed using standard two dimensional (2D) echocardiography. However, some limitations are inherent to 2D echocardiography measurements. Probe positioning between the ribs of the patient can be challenging and finding the standard cross-sectional view of the heart requires an important training of the cardiologists or radiologists. The positioning induces inter-operator variability which is often responsible for low result reproducibility [START_REF] Sahn | Recommendations regarding quantitation in M-mode echocardiography: results of a survey of echocardiographic measurements[END_REF]. In addition, some approximations are necessary to quantify three dimensional (3D) volumes from 2D cross-sectional images. For instance, Cardiac Output is calculated from Left ventricular outflow tract velocity time integral (LVOT VTI) which requires the measurement of the aorta cross section and the assumption of a cylindrical geometry. This approximation is responsible for an error of more than 10% on the value of the Cardiac Output [START_REF] Huntsman | Noninvasive Doppler determination of cardiac output in man. Clinical validation[END_REF]. 3D echocardiography could be the answer to increase repeatability and reduce inter-operator dependency [START_REF] Dorosz | Performance of 3-Dimensional Echocardiography in Measuring Left Ventricular Volumes and Ejection Fraction: A Systematic Review and Meta-Analysis[END_REF]. It is actually used in the clinic to perform anatomical measurements such as the ejection fraction [START_REF] Dorosz | Performance of 3-Dimensional Echocardiography in Measuring Left Ventricular Volumes and Ejection Fraction: A Systematic Review and Meta-Analysis[END_REF]. However, due to its low volume rate, standard 3D echocardiography is not able to map flow velocity over 3D volumes due to low frame-rate of the line-byline ultrasonic wave transmit scheme used in current clinical 3D ultrasound systems.

In research, ultrafast imaging has been studied for more than two decades [START_REF] Tanter | Ultrafast imaging in biomedical ultrasound[END_REF]. It enables to increase the frame rate in 2D to reach few kilo images per second. The method relies on the emission of unfocused wave to insonifiate the entire medium in few transmits. Recently, ultrafast imaging was extended to 4D ultrasound imaging [START_REF] Provost | 3D ultrafast ultrasound imaging in vivo[END_REF] for different applications: to image blood flow in small vessels [START_REF] Provost | 3-D ultrafast doppler imaging applied to the noninvasive mapping of blood vessels in Vivo[END_REF], shear waves [START_REF] Gennisson | 4-D ultrafast shear-wave imaging[END_REF], pulse wave [START_REF] Apostolakis | Feasibility and Validation of 4-D Pulse Wave Imaging in Phantoms and In Vivo[END_REF], quasi-static [START_REF] Papadacci | 3D Quasi-Static Ultrasound Elastography With Plane Wave In Vivo[END_REF] and dynamic strain [START_REF] Papadacci | 3D Myocardial Elastography In Vivo[END_REF]. In the heart, 4D ultrasound ultrafast imaging was performed to image blood flow [START_REF] Provost | 3D ultrafast ultrasound imaging in vivo[END_REF] and fibers [START_REF] Papadacci | High-contrast ultrafast imaging of the heart[END_REF] in the left ventricle of a human volunteer and tissue displacement in canine's hearts [START_REF] Papadacci | 3D Myocardial Elastography In Vivo[END_REF]. One very interesting property of ultrafast imaging is that it enables quantification. In the case of Doppler blood flow, one spectrogram can be retrieved from each voxel which allows quantification as if pulsed Doppler was performed at each location of the volume. Another advantage of 4D ultrafast Doppler is that, flow rate evaluation can become independent from the hypothesis on the geometry and the angle of insonification. Indeed, the flow speed can be evaluated over an entire cross-section and the normal to the cross-section can be chosen to be collinear to the ultrasonic beam axis. In this study, we developed this flow rate evaluation method based on 4D ultrafast Doppler independent of assumptions and beam angle. The method was evaluated in vitro on a pipe with a controlled flow rate (N=7).

II. MATERIAL AND METHODS

A. Experimentation

A soft plastic pipe (internal/external diameter 1.8/2.7 cm) was connected to a peristatic pump. The pipe was maintained vertically in a water tank with small nylon thread attached to the side of the tank. The ultrasonic probe was placed on top of the pipe's end to mimic the in vivo aortic flow measurement in apical view condition. A solution of water with corn starch (6%) was infused through the pipe. Ultrasonic measurements were performed for N=7 different flow rates. The flow was monitored with a flowmeter. A scheme of the experiment is displayed in figure 1.

Figure 1 : Scheme of the experiment. A pipe filled with water and corn starch, connected to a peristatic pump, is immerged in a water tank. A matrix array probe connected to a 3D ultrafast scanner is set vertically on top of the pipe to image the outgoing flow.

B. Ultrasonic acquisition

An ultrasonic matrix array probe (central frequency 2.5MHz, 1024 elements, pitch 0.3 mm, bandwith 60%, Vermon, France) connected to a 1024 channels ultrasound scanner prototype was used to image the pipe in three dimensions. 2000 diverging waves (angular aperture 80° [START_REF] Papadacci | High-contrast ultrafast imaging of the heart[END_REF]) were emitted at a volume rate of 2000 volumes/s during 1s. Radio Frequency data were received from each transmitted wave at a sampling frequency of 10MHz and stored in memory.

C. Post processing a. Beamforming

The 500 volumes (dimensions 80°x80°x50mm, sampling 2°x2°x0.3mm) were reconstructed. A conventional delay and sum algorithm was implemented on GPU (Titan Xp, nvidia) using cuda language. High pass filter (cut-off frequence 20Hz) was applied to remove low frequency motion and keep only the signal from the flow. Power Doppler image was computed by summing the 500 filtered images over time. In each voxel, Doppler spectrums were computed by performing the Fourier transform over time. Color Doppler volumes (quantitative flow speed volumes) were computed by calculating the first moment of the Doppler spectrums in each voxel. The Color Doppler volumes were averaged over time to decrease noise. Power velocity integral (PVI) was calculated by multiplying the power Doppler and the Color Doppler volumes. The flow volume was obtained by thresholding the PVI while keeping the same surface S for each flow rate acquisition.

b. Flow rate

The pipe flow rate (Q) at a depth z is given by the following expression:

= . → ( )
is the flow speed (m/s). Here we chose an arbitrary surface S perpendicular to the ultrasonic beam axis, we obtain the integration of the axial velocity over S. In other word, the flow rate is the integral of the Doppler velocity estimates over the cross section of the flow. and are collinear as the normal to the surface is always along the beam axis before scan conversion process.

c. Data representations

3D volumes were then scan-converted and represented using Amira Software. Calculation was performed using Matlab. Beamforming was performed on a TITAN XP from Nvidia Corporation.

III. RESULTS

A. Power Doppler

The power Doppler volume was successfully observed allowing 3D representations in top view (figure 2 A) and side view (figure 2 B). In each voxel of the volume, a Doppler spectrum was assessed (i.e. figure 2 C). The first momentum was associated to the flow speed. In other word, in each voxel of the volume the flow speed was assessed.

B. Color Doppler

Color Doppler volume was successfully observed in side view (figure 3 A) and top view (figure 3 B). The 3D visualization enabled to visualize flow going in two different directions: blue color going towards the probe, red color going away from the probe (see white arrows). The flow going towards the probe was associated to the outgoing flow while the flow going away from the probe was interpreted as the fallout of the flow along the pipe. Power velocity integral method allowed us to obtain a composite image between the power Doppler and the Color Doppler (figure 3

C. Evaluation of flow rates in the pipe

The pipe output was calculated after application of the mask and flow speed integration at the output of the pipe. The comparison between output assessed with the flowmeter and the output calculated from the 4D ultrafast method is presented in figure 5. A linear regression was performed and gave an excellent correlation between the measured ouputs ( = 0.96).

IV. CONCLUSION

In study, we developed a new method to evaluate the flow rates in a pipe based on 4D ultrafast Doppler imaging. The measured flow rates were found to be in a good agreement with the flowmeter values used as a gold standard ( = 0.96). The 4D approach enabled to evaluate a spectrogram at each voxel of the volume from which a blood flow speed was assessed. Flow rates were assessed on a cross section of the tube using the estimated blood flow speed at each voxel of the cross section. This method presents the advantage of being independent of the direction of the ultrasonic beam. No assumptions were made on the pipe geometry. Limitation includes the incapability for the pump to deliver higher flow rate (1L/min max). However, it would still be possible to perform the method at higher flow rate as long as no aliasing is present in the Spectral Doppler signal. In case of aliasing the method could be coupled to dealiasing scheme [START_REF] Posada | Staggered multiple-PRF Ultrafast Color Doppler[END_REF]. The method would be extremely valuable to estimate Cardiac Output in acute heart failure for therapeutic optimization for instance. The four dimensional nature of the acquisition has the potential to enable the evaluation of the Cardiac Output in vivo in patients without the need of making any assumptions on the valve geometry or the direction of the ultrasonic beam usually responsible for errors.

Figure 2 :

 2 Figure 2 : 3D representation of Power Doppler in the pipe with a Doppler spectrum associated with one voxel of the volume. Cross section of the Power Doppler in top view (A) and side view (B) are presented. A Doppler spectrum can be derived for each voxel of the volume (C).

  C D). A threshold was applied on the PVI to obtain the volume of outgoing flow. (figure 4 A B).

Figure 3 :

 3 Figure 3 : 3D representation of Color Doppler flow and Power velocity integral in the pipe. Cross section of the Power Doppler in top view (A) and side view (B) are presented and displays flows in two directions. Cross section of the Power velocity integral in top view (A) and side view (B) are displayed and lobes from Color Doppler flow are reduced.

Figure 4 :

 4 Figure 4 : 3D mask associated with the flow in the pipe in side view (A) and top view (B). It was computed by thresholding the Power velocity integral.

Figure 5 :

 5 Figure 5 : Flow rate measured with 4D ultrafast Doppler against flow rate measured with the flowmeter is presented. A linear regression was applied to the data. A high correlation was found ( = 0.96)

V. ACKNOWLEDGEMENTS

This study was supported by the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n° 311025 and the ANR-10-IDEX-0001-02 PSL* Research University. We acknowledge the ART (Technological Research Accelerator) biomedical ultrasound program of INSERM. The Titan X Pascal used for this research was donated by the NVIDIA Corporation.

VI. BIBLIOGRAPHY