Emmanuel Prados

Patrick Criqui

Constantin Ilasca

Mathematical formulation of REDEM algorithm for National soft landing CO2 trajectories under global carbon budgets

Keywords: global warming, climate change, international climate negotiations, benchmark, national GHG emission curves

Global warming may be one of the greatest threats facing the human civilization. It is now widely shared that it is necessary to reduce quickly and significantly the greenhouse gas emissions to avoid uncontrolled and irreversible evolutions of climate. It has now become urgent to develop a legal instrument addressing the post-2020 period and to achieve a successful outcome in the international climate negotiations. In this paper we propose a new computational tool which provides elements of benchmarking for the climate negotiations. The model and algorithm we propose is designed on rationale elaborated by energy and climate policy experts. We detail how to estimate the parameters of the model and how this benchmarking tool could be used.

Formulation mathématique de l'algorithme REDEM, pour des trajectoires de référence des émissions nationales de CO2 sous contrainte de budget carbone global

Résumé :

Le changement climatique est probablement l'une des plus grandes menaces à laquelle la civilisation humaine doit faire face. Il est désormais largement partagé qu'il est nécessaire de réduire rapidement et de manière très importante les émissions de gaz à effet de serre afin d'éviter un emballement irréversible du système climatique. Il est aujourd'hui urgent de développer les instruments politiques et juridiques pour la période post-2020 et de faire aboutir les négociations internationales. Dans ce rapport, nous proposons un nouvel outil permettant de calculer des courbes de références des émissions nationales de gaz à effet de serre. L'algorithme et le modèle proposés se basent sur une logique développée par des économistes spécialistes de l'énergie et du changement climatique. Nous montrons comment les paramètres du modèles peuvent être estimés et comment ces outils peuvent être utilisés.

Mots-clés : réchauffement climatique, négociations internationales, émission de gaz à effet de serre, courbes de référence

Introduction

Climate change and global warming may be one of the greatest threats facing the human wellbeing. The world is not on track to stay within the 2°C limit of temperature rise. The consequences would be catastrophic [START_REF] Stocker | Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment[END_REF]. Howsoever bad the situation may seem, the latest scientific research indicates that keeping below the 2°C limit is challenging but feasible. Nevertheless only an internationally coordinated, goal oriented approach to operationalizing the 2°C limit will allow humanity to avoid dangerous climate change [START_REF] Guerin | Pathways to deep decarbonization: Interim 2014 report[END_REF].

Compared to expectations, the 16th Conference of the Parties to the UNFCCC at Copenhagen was a failure. The main outcome of COP-17 is the decision to launch a new negotiation process in order to develop a "legal instrument, or agreed outcome" addressing the post-2020 period and "applicable to all Parties". To break the deadlock in international climate negotiations is more than ever crucial and urgent.

The installation of an international climate regime is a complex process. One of the main difficulties of the next climate accords is that they must be able to combine a bottom-up process and a top-down rationale. In effect, since the beginning, the distinction made in the international negotiations -between the Annex I (developed countries) and non-Annex I (developing countries)-drives to a greater emphasis on differentiated responsibilities than on common responsibilities [START_REF] Aldy | Climate negotiations open a window: Key implications of the durban platform for enhanced action[END_REF]. Then this drift has been significantly increasing when developing countries stepped forward to join the main emitters in the limelight. In particular, the conclusions of the COP-16 at Copenhagen clearly show that the will of most Parties is to propose their own pledges to emission reductions. Continuing to follow the bottom-up approach is feasible but not desirable : when added up, the actual national pledges lead to a 3-4°C warming of the temperature [START_REF] Vieweg | Climate shuffle climate action tracker update[END_REF]. On the other side, the outcomes of COP-17 at Durban recognize this emission gap and mention the 2°C target which clearly illustrate the importance of the top-down rationale. Whatever the protocol will come from the negotiations by 2015, this new accord has to include a differentiation criterion that doesn't violate basic equity principles in emission reduction and a regulatory approach which may provide a reasonable chance to reach the 2°C target. This top-down rational is desirable but for the moment hardly feasible.

In this context, we propose a tool which attempt to reconcile bottom-up aspirations and top-down scheme. More precisely, this paper proposes a harmonized approach that aims at developing a common benchmarking tool for national emission reduction policies. Indeed, the decision to keep the 2°C temperature increase as a reference and the mention to the legal form need to create tools which allow to compare GHG emission reductions of the various countries and put into coherence.

This technical report describes the mathematical model and the associated algorithms of the computational tool we propose. For a more complete introduction, a comparison with the state of the art, and a deep analysis of the evolution and actual situation of international climate negotiations, we refer the reader to our previous report [START_REF] Criqui | National Soft Landing CO2 trajectories under global carbon budgets[END_REF].

A benchmarking tool which reconciles bottom-up aspirations and top-down scheme

In that perspective, our goal is to identify practical solutions for reaching the Convention's main objectives. In order to reconcile bottom-up aspirations and top-down scheme, the solutions have to verify the features expressed in the Parties' positions as well as in the latest developments of the discussions within the ADP 1 process. The key features, that the benchmarking system for national pledges has to verify, can be synthesized as follows:

1. The objectives and the benchmarking profiles must be differentiated, based on the countries' responsibility and capability (and not on the affiliation to different "clubs" or stages);

2. The reasoning should be based on the simplest, yet the most uncontested variables (population, GDP, level of emissions) in order to develop differentiation indicators which could be least subject to controversy;

3. The national pledges, formulated in a bottom-up manner, should not be leading to an imposed emission level. Transparency, visibility and common accounting rules could be provided through an international benchmark system. 4. The scheme should be submitted to a measurement and verification process in order to insure the reference objectives and the convergence of the emission reduction efforts.

5. The benchmark process should be able to ensure the connection between the two approaches (top-down and bottom-up). More precisely, it should show either the coherence or the discrepancy between the global and the national emission profiles.

6. The differentiating system should be based on an economical perspective, ensuring that no shock is imposed to the national economies through emission reduction targets.

7. The scheme should be submitted to a transparent measurement and verification process in order to ensure the convergence of the emission reduction efforts towards the 2°C target.

In the sequel of this report, we propose a tool that is designed for the benchmarking of national emission reduction trajectories and which fulfills these underlying principles. The idea is to develop a tool which generates curves which allow the countries to evaluate and compare themselves to others. The comparisons must be done in a framework ensuring the consistency with the global objective of GHG emissions. This framework must verify some common equity rules which are the simplest and the most objective possible and which allow to differentiate the countries according to their capacity and responsibility. The obtained curves are benchmark curves (of future GHG emissions until 2100). They are not forecasting curves based on simulations of complex economic phenomena.

To differentiate the countries, we propose here to use a differentiation indicator which is a simple combination of each country's per capita emission and per capita GDP. The per capita emission can be considered as a reasonable proxy of a country's responsibility, and the per capita GDP a proxy of its capacity for action. In the following, we call this indicator, the CRI (Capacity-Responsibility Indicator).

One of the characteristics of our method is that it principally focuses on the rate of variation of GHG emissions and on the total emission budget for the period considered. In addition, it is based on the fundamental observation that when emission curves follow a pattern of type "Peak, plateau and decline" [START_REF] Yawitch | South africa's long term mitigation scenarios: Process and outcomes[END_REF] with a final stabilization at a low level, the rates of emission growth 1) decrease, then 2) reach a minimum and then 3) slowly increase and converge to zero.

More precisely, in our model, we focus here on the peak in the rate of reduction of emissions, or decarbonization rate. This peak -indeed rather a trough -takes place later than the emission peak, and can be interpreted as a peak of effort. We propose to parameterize the maximum effort for each country according to the above-mentioned capacity-responsibility indicator in a "proportional" way. By doing so, we then obtain a common benchmarking method which proposes differentiated objectives based on a simple indicator. In addition, the algorithm is designed so as to compute emission curves for all the countries in such a way that the global aggregated Inria emission budget is consistent with the chosen temperature target (for example, the 2°C climate target).

In [START_REF] Criqui | National Soft Landing CO2 trajectories under global carbon budgets[END_REF], the reader can find much more details and discussions about the underlying logic to our proposition. We also describe all interface parameters in a simple and intuitive manner. In the following, we detail the mathematical model we propose and associated algorithms which allow to compute the parameters.

Mathematical formulation

Capacity-Responsibility Index

To be consistent with the first two rules mentioned above, we then propose to enforce the effort of reduction of GHG of states to directly depend on a Capacity-Responsibility Index (CRI) which combines: 1) the per capita emissions of GHG and 2) the per capita PIB for each state. These indicators are fixed once and for all (for the whole considered period). They are computed at a reference year that has to be negotiated. Here we fix this reference year to 2010.

Let i ∈ {1 ... n} be the indices corresponding to the states or groups of states. Let p i ≥ 0 be the per capita PIB for the state i (or group of states i), expressed in k$(in thousands of US Dollars). Let e i ≥ 0 be the per capita emissions of GHG for the state i (or group of states i), expressed in tons (CO 2 and other GHG).

We denote N i the CRI of state i (or region i). We propose to define the CRI indexes by using the classical 1-norms used in R 2 :

N i = rp i + e i δ ,
where δ = max i (rp i + e i ). For each state i, we have then 0 ≤ N i ≤ 1. Let us also note that, by construction, there exists a specific i such that N i = 1.

In some sense, the CRI allows to project the 2-dimensional distribution of country indicators (p i , e i ) to a 1-dimensional distribution included in the range [0, 1]. The parameter r is an important parameter and must be carefully considered. It directly and significantly impacts the repartition of the indicator.

When the isolines are horizontal lines, then they are orthogonal to the trajectories (in the space (p, e)) that would follow a state which would only reduce his GHG emission (without changing his GDP per cap.). Those trajectories would then be collinear to the gradient of the CRI. For such states, a reduction of the GHG emissions would be widely rewarded, when a variation of their GDP per capita would not be penalized. Inversely, when the isolines are vertical then variations in the GHG emissions do not change CRI whereas GDP variations modify in a maximal way this index.

Model formulation and constraint specification

Let E i be the total emissions of GHG for the state or region i (in million of tons). For all states and regions, the emissions E i of GHG evolve over time; in other words, E i is a function of time E i : R → R : t → E i (t). Their derivatives E i (t) measure the variations of the emissions E i (t); more exactly they measure their growing. Now let us consider the decarbonization rate R i (t) = -100

E i (t)
Ei(t) . This ratio is a percentage that can be interpreted as an indicator of efforts to make. The effort is maximal, when R i (t) is maximal.

We suggest to use the CRI indexs to introduce equity in the reduction efforts of GHG emissions. We propose to perform this by introducing the CRI indexs in the peak of effort. To simplify, we propose that the date of peak of effort is the same for all the states and regions (if it is relevant and useful, an eventual relaxation of this constraint will be the subject of a future paper). We call t max this date that has to be negotiated. A reasonable choice for t max could be 2050. At t max , we propose to force where µ = 1 n i N i is the mean of the CRIs (of all the states and regions) and σ their standard deviation. σ is a parameter which has to be negotiated. In a sense, here we have rescaled and positioned the CRIs in an "effort" range via a simple affine transformation of the CRIs; by changing their mean and their standard deviation: see Figure 1. Parameter σ gives more or less weight to this equity desire. It allows to weaken or to magnify the differences of treatment imposed to the states and regions. It gives more or less importance and effects to the CRIs. When σ = 0, there is no differentiation ; this would correspond to a one-size-fits-all approach. In the following, we are going to denote

R i (t max ) = μ + σ σ (N i -µ), Inria
Ñi = μ + σ σ (N i -µ).
In the sequel we will explain how to get a relevant value for μ.

The reduction efforts being assumed to be maximal at t max , we have then

R i (t max ) = 0.
We can also have at our disposal some approximations of the actual efforts D i made by all the states. So we have to enforce R i to verify

R i (t 0 ) = D i ;
i.e. R i (t 0 ) is fixed and known for a reference year t 0 . In our experiments, we have fixed t 0 = 2010. In other respects, to ensure the 2°C global objective, we have to force the global GHG emission budget to do not exceed a certain bound. In other words, we have to enforce the following constraint:

i t budget t0 E i (t)dt ≤ B budget , (1)
where B budget is an adequate bound (in million of tons) and [t 0 , t budget ] is the budget period. In our experiments, we have chosen t budget = 2100.

Mathematical formulation in brief

The reference emission curves E i of the benchmarking tools have then to verify the following equation:

100 E i (t) + R i (t)E i (t) = 0, (2)
with E i (t 0 ) = E 0i which are the (known) emissions of GHG at the reference year t 0 . Equation ( 2) is a classical homogeneous linear differential equation of order 1. In other respects, the R i functions have to verify

R i (t 0 ) = D i , (3)
R i (t max ) = 0.
It is also natural to impose that the R i are regularly increasing between t 0 and t max , and that they are regular and monotonic after t max .

In addition, the global emission curves E i have to verify

i t budget t0 E i (t)dt ≤ B budget . ( 4
)
In the following we show how to compute emission curves E i that verify these equations and constraints.

Computation of emission curves

Resolution of the differential equation (2)

While the whole problem formulated by the set of the equations ( 2), ( 3) and ( 4) is not well-posed (in the sense that it has in general an infinite number of solutions as it will be illustrated in the following), the ordinary differential equation ( 2)

100 E i (t) + R i (t)E i (t) = 0, E i (t 0 ) = E 0i , (5)
is completely well-posed. Also it is quite well-known that this differential equation has a unique solution:

E i (t) = E 0i e -1 100
t t 0
Ri(s)ds , (6) see for example [START_REF] Robinson | An Introduction to Ordinary Differential Equations[END_REF]. The issue here is then to find some adequate decarbonization ratio curves R i .

Decarbonization ratio curves R i

The goal here is then to propose some adequate decarbonization ratio curves R i . Clearly some ambiguities appear because the three conditions (3) are not sufficient to completely determine these curves. We then propose to select some simple and reasonable solutions. On the interval [t 0 , t max ], we then pose

R i (t) = D i -Ñi t max -t t max -t 0 γ + Ñi ,
where c can be interpreted as an acceleration coefficient. This last parameter also indicates the degree of convexity of the curves of the rate of variation. We fix this parameter to be identical for all countries. The value of γ is exogenous and have to be negotiated between all the Parties. By default, we fix γ = 2 in our experiments. To simplify notations, we pose

a i := D i -Ñi (t max -t 0 ) γ Inria and b i := Ñi . So R i (t) = a i (t max -t) γ + b i .
In the long term, emissions have to be stabilized. It is then necessary that R i (t) regularly converges towards zero when t ≥ 0. We then propose the following function for all t ≥ t max :

R i (t) = Ñi e -(t-tmax ) 2 2θ 2
.

θ is another acceleration term, which plays a similar role as γ .

Emission curves E i before the peak of effort

The reader can easily verify that for t in the interval [t 0 , t max ], we have

t t0 R i (s)ds = -a i (t max -t) γ+1 γ + 1 + b i t + c i ,
where

c i = a i (t max -t 0 ) γ+1 γ + 1 -b i t 0 = D i -Ñi t max -t 0 γ + 1 -Ñi t 0 . (7)
Then, for t ∈ [t 0 , t max ], we have

E i (t) = E 0i e 1 100 [ai (tmax-t) γ+1 γ+1 -bit-ci] . (8)

Emission curves E i after the peak of effort

We have:

t tmax R i (s)ds = t tmax Ñi e -(s-tmax ) 2 2θ 2 ds = θ π 2 Ñi erf t -t max √ 2θ , (9)
where the Gauss error function is the (non-elementary) sigmoid function

erf(t) = 2 √ π t 0 e -s 2 ds.
So, for t ≥ t max , we have

E i (t) = E 0i e -1 100 bitmax+ci+θ √ π 2 Ñi erf t-tmax √ 2θ . ( 10
)

Enforcing the 2°C global objective

The issue here is to find a value for the parameter μ such that constraint (4) reminded below is verified

n i=1 t budget t0 E i (t)dt ≤ B budget .

RR n°8601

In practice the states will always try to minimize their efforts, so the constraint is de facto

n i=1 t budget t0 E i (t)dt = B budget (11)
instead of equation ( 4). Since Ñi are increasing with respect to μ and R i are increasing with respects to Ñi , we obtain R i increasing with respects to μ. With similar monotonic arguments we can conclude that i t budget t0 E i (t)dt is decreasing with respect to μ. So there exists at most one μ which verifies equation (11). Let us remark that t budget t0 E i (t)dt is an exponential integral that cannot be written as a combination of elementary functions. So we propose to solve equation (11) by using the Newton method. Let us note

f (μ) = n i=1 t budget t0 E i (t, μ)dt -B budget .
Starting from a first approximation μ0 of the solution, Newton method consists in computing a sequence of approximations of the solution by performing the following iterations until the convergence:

μk+1 = μk - f (μ k ) f (μ k ) .
To apply this method, we need to calculate the derivative of f with respect to μ. Here we have:

f (μ) = n i=1 tmax t0 ∂ μE i (t, μ)dt + t budget tmax ∂ μE i (t, μ)dt
where, for t ≤ t max , we have

∂ μE i (t, μ) = E i (t, μ) 100 a i (μ) (t max -t) γ+1 γ + 1 -b i (μ)t -c i (μ) .
Since

a i (μ) = - 1 (t max -t 0 ) γ , b i (μ) = 1 , c i (μ) = - t max + γt 0 γ + 1 ,
then, for t ≤ t max , we have

∂ μE i (t, μ) = 1 100 E i (t, μ)[- (t max -t) γ+1 (γ + 1)(t max -t 0 ) γ -t + t max + γt 0 γ + 1 ].
For t ≥ t max , we have

∂ μE i (t, μ) = - 1 100 E i (t, μ) t max - t max + γt 0 γ + 1 + θ π 2 erf t -t max √ 2θ .

Inria 6 Examples of use of the algorithm

Our tool aims at illustrating the conditions of convergence between the national pledges and the 2°C target.

In practice, the parameters r, t 0 , t max , σ, θ and γ are fixed by users. They result of a negotiation by all the Parties. Then the algorithm computes the emission curves E i for all countries (whose curves of rate of variation verify the parameterization defined above) such that the objective for the expected increase in temperature is verified (the temperature target being also indicated as parameter). The actual version of the program produces a set of curves: first it displays the curves of the rates of variation of emissions (in % per year), the emission trajectories, the emissions per capita, and a chart illustrating the distribution of the budget of emissions (budget between 2000-2100) for all countries (see Figure 2). Finally, it is also possible to get the rate of variation curves and emission curves for one given country as obtained for several budgets on the same graph, as illustrated in Figure 3 (all the other parameters being fixed). This tool would thus allow to easily visualize "national decarbonisation corridors" which are globally consistent for all the countries and verifying the rules mentioned in the previous sections.

This program can be used to examine different "effort profiles". For example, various scenarios can allow to compare an "early" and a "delayed" action, the impacts of low and high levels of probability of attaining the 2°C target and finally low and high levels of differentiation across countries. Such scenarios are characterized by changes in the values of the parameter t peak , of the emission budget and of the differentiation parameter σ . The results obtained with several examples of scenario are detailed in our technical report [START_REF] Criqui | National Soft Landing CO2 trajectories under global carbon budgets[END_REF]. We refer interested reader to this report because of space constraint. The current process of climate negotiation is focused on the framing of an outcome supposed to be applicable to all the countries. Meanwhile, the sum of the mitigation contributions -which are to be based on national circumstances and policies, should be compatible with the 2°C target. In this context our approach provides a rationale to produce mitigation trajectories that are consistent in a cross-country analysis with the control of the increase in world temperature.

Our tool shows a practical way to guide the national trajectories or pledges through a convergence mechanism into a comprehensible framework. The algorithm starts from the consideration of a generic emission profile (peak, plateau and decline) and a consistent way to relate the carbon budget in 2100 to the temperature increase. At the same time, it considers the equity matter, which is provided through the Capacity Responsibility Indicator and afterwards by the use of a standard deviation for the differentiation of the national decarbonisation trajectories. The main feature of the algorithm consists in the simulation of the decarbonisation rate for each individual country while it also fully accounts for total emissions and emission budgets over the period. We detail in this report the mathematical model we propose and how we compute the benchmarking emission curves and the associated parameters.

We finally rapidly illustrate how the tools could be used. The simulations describe different settings for national and global emission profiles. They illustrate different options and tradeoffs for the international community so as their potential consequences in terms of national trajectories. The logic of providing national decarbonisation corridors, representing different stringency efforts, may provide a help for implementing a benchmarking system for national decarbonisation policies. Thus the use of our tool or of its results may be a help in designing a better alignment of national contributions towards the 2°C target.

Figure 1 :

 1 Figure 1: Parameterization the maximum effort for each of the countries with its CRI.

Figure 2 :

 2 Figure 2: Examples of output curves for all countries: rate of variation of emissions, emission trajectories, emissions per capita, distribution of the budget

Figure 3 :

 3 Figure 3: Example of national decarbonisation corridors

3 . 2 3 . 3 7 45

 32337 -Responsibility Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Model formulation and constraint specification . . . . . . . . . . . . . . . . . . . 6 Mathematical formulation in brief . . . . . . . . . . . . . . . . . . . . . . . . . . Computation of emission curves 8 4.1 Resolution of the differential equation (2) . . . . . . . . . . . . . . . . . . . . . . 8 4.2 Decarbonization ratio curves R i . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 4.3 Emission curves E i before the peak of effort . . . . . . . . . . . . . . . . . . . . . 9 4.4 Emission curves E i after the peak of effort . . . . . . . . . . . . . . . . . . . . . . 9 Inovallée 655 avenue de l'Europe Montbonnot 38334 Saint Ismier Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

RR n°8601

Prados & Criqui & Ilasca


Amandine Dubois

Francois Charpillet

A Gait Analysis Method Based on a Depth Camera for Fall Prevention

This paper proposes a markerless system whose purpose is to help preventing falls of elderly people at home. To track human movements, the Microsoft Kinect camera is used which allows to acquire at the same time a RGB image and a depth image. Several articles show that the analysis of some gait parameters could allow fall risk assessment. We developed a system which extracts three gait parameters (the length and the duration of steps and the speed of the gait) by tracking the center of mass of the person. To check the validity of our system, the accuracy of the gait parameters obtained with the camera is evaluated. In an experiment, eleven subjects walked on an actimetric carpet, perpendicularly to the camera which filmed the scene. The three gait parameters obtained by the carpet are compared with those of the camera. In this study, four situations were tested to evaluate the robustness of our model. The subjects walked normally, making small steps, wearing a skirt and in front of the camera. The results showed that the system is accurate when there is one camera fixed perpendicularly. Thus we believe that the presented method is accurate enough to be used in real fall prevention applications.

I. INTRODUCTION

The ageing of the population is a major problem for the next years. One of the problems generating a loss of autonomy is fall. Each year a third of the persons aged 65 years and older fall [START_REF] Bourdessol | Prévention des chutes chez les personnes âgées à domicile[END_REF]. Falls have multiple risk factors: age, muscular weakness, gait modification, obstacles in the environment (carpet)... Fall risk can be detected by measuring gait modifications.

Several systems exist for computing gait parameters. There are systems like accelerometers, gyroscopes which measure the angular velocity of the segment on which they are put. There exist optical systems (as the VICON system) for movement tracking based on markers fixed on certain key points of the body. Other systems without wearable sensors exist as for example force platforms, actimetric carpets, on which the person walks and which detect the exact moment of the foot-ground contact and calculate the velocity and the length of steps. There are also telemeters computing velocity. We would like to develop a system based on a low cost and not imposing ambiant sensor (without wearable markers). In the litterature, some authors begin to be interessed in developing ambiant systems allowing to analyze the gait of a person. In 2009 Hagler et al [START_REF] Hagler | Unobtrusive and Ubiquitous In-Home Monitoring: A Methodology for Continuous Assessment of Gait Velocity in Elders[END_REF] installed four PIR motion sensors on the ceiling to obtain the speed of the gait. Auvinet et al [START_REF] Auvinet | Gait analysis with multiple depth cameras[END_REF] in 2011 are interessed in detecting symmetrical gait troubles with three RGB-D cameras. In 2011 Stone et al [START_REF] Stone | Evaluation of an Inexpensive Depth Camera for Passive In-Home Fall Risk Assessment[END_REF] developped a system to evaluate a fall risk with an inexpensive depth camera. The idea is to measure some gait parameters of the person with a RGB-D camera.

The aim of our system is to record and to analyze the gait parameters of a person regularly enough to note a possible evolution over time. On some period, the aim is to know if the gait parameters tend to remain stable or if they are evolving which would suggest an increased risk of fall. In the last case, the results could be transmited to the physician.

The proposed system uses a RGB-D camera chosen for the advantage of the depth reconstruction, its low cost and its invariance to visible lighting. We imagine to fix the system high in a corner at home. Several stages are necessary to achieve our goal. The system must be able to detect the person in the field of view, to detect when she walks and finally to analyze her gait. The algorithm that we use to extract the person and to detect her activity (as walking) is presented in [START_REF] Dubois | Human Activities Recognition with RGB-Depth Camera using HMM[END_REF]. In this article, we present an algorithm for analysing the gait. Our method consists in extracting gait parameters, of interest for fall prevention, from the trajectory of the center of mass.

This paper is organized as follows. Section II presents background elements on activity recognition that we use for detecting walking situations. Section III presents the algorithm allowing to analyze the gait of a person. Section IV presents the evaluation of our algorithm.

II. BACKGROUND

This work is the continuation of previous work on fall detection with RGB-D camera [START_REF] Dubois | Human Activities Recognition with RGB-Depth Camera using HMM[END_REF]. The detection of the falls was included in the more general problem of activity recognition. The method consists in tracking a person in the field of the view of the camera. To extract the person from the background, a standard method in computer vision is used. A Hidden Markov Model (HMM) is used to analyze the activity of the person. Each state represents an activity that the person can do in her daily life. Eight activities were defined: walking, lying (on a bed, on a couch for example), sitting, falling, lying down, squatting, going up on an obstacle (a chair, a footboard for example) and bending. These activities are distinguished using 3 parameters: the position of the center of mass of a person on vertical axis, the vertical speed and standard deviation of the person's body.

To know in what state the person is, the Forward-Backward algorithm is used. 26 healthy subjects, of 20 at 53 years old, performed the 8 activities, which represented a training database of 208 sequences. The observation probabilities of the three parameters for each of the 8 situations were learnt from of 128 sequences. The validity of the algorithm was been tested with 80 sequences. The results showed that the algorithm could distinguished the different activities. The only error was for the state "bending" which was replaced by state "sitting" or "squatting".

The HMM algorithm provides the state of the person and the general idea is to analyze the gait of the person when she is in the state "walking". The following stage consists in analyzing the gait.

III. METHOD

The background extraction is realized using the "running average" method as presented in [START_REF] Dubois | Human Activities Recognition with RGB-Depth Camera using HMM[END_REF]. In this part, the method for tracking the person and analyzing the gait is presented.

A. Tracking of the person

To obtain the displacement of the walking body, our approach consists in analyzing the displacement of her center of mass. The horizontal and vertical displacements of the center of mass when a person walks is described in Figure 1(a The background extraction allows us to have the mobile pixels belonging to the person. To calculate the center of mass of the person, the method consists in averaging the location of all the mobile pixels of the person. The displacement of the center of mass of a person walking perpendicular to the camera, obtained with our method, is drawn in Figure 2. The horizontal trajectory of the Figure 2(b) is smoother than the trajectory described in litterature. The vertical trajectory of the Figure 2(a) is similar to the theoretical sinusoidal curve described by Saunders in Figure 1(b). Only the displacement on the vertical axis is discussed in the following section. The center of mass trajectory is smoothed with a Kalman filter [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF].

B. Analysis of the gait 1) Selection of gait parameters:

To analyze the gait of a person, some parameters can be cited according the purpose of the analysis. A first category corresponds to the spatiotemporal parameters, as the length of steps, the duration of steps, the speed of the gait, the time propulsion, which are widely studied because they illustrate the characteristics of the gait. Then, the kinematic parameters, as knee angle, are the indicators reflecting the movement of the lower limbs during walking. And we can also cite the dynamic parameters, as the reaction forces from the ground and the center of pressure, allowing an evaluation of the foot. Our aim is to develop a system for fall prevention through a global gait analysis, so the first category (the spatiotemporal parameters) proves to be more suitable. Several articles show that some gait parameters are indicators of an increased fall risk. Hausdorff et al [START_REF] Hausdorff | Gait Variability and Fall Risk in Community-Living Older Adults: A 1-Year Prospective Study[END_REF], in 2001, on a one year experiment showed that the stride time variability, the swing time variability and the inconsistency of the variance is associated with an increased risk of future fall. Logistic regression also showed that stride time variability predicts falls. Auvinet et al [START_REF] Auvinet | Gait Abnormalities in Elderly Fallers[END_REF] showed in 2003 that the speed, the frequency and the length of the steps had reduced values for fallers. The irregularity of the length of steps is considered by them as a reliable variable to predict falls. We are interested in the gait parameters as indicators for the prevention of falls. Therefore we focus, according to the litterature, on measuring the length and duration of the steps and the speed of the gait 2) Extraction of the selected gait parameters: The gait parameters (length and duration of steps and the speed of the gait) are extracted from the vertical displacement of the center of mass. In medicine the definition of the length of steps is the distance separating the two heels at the time of double support, as illustrated in Figure 4. Our method consists in searching the position of the supporting foot. The center of mass is located at the vertical of the supporting foot when it is highest as shown in Figure 3. The position of the supporting foot is estimated by the projection on the ground of the center of mass (at the time of the maxima). The length of step is calculated as the distance between the projection of two local vertical maxima of the trajectory of the center of mass. In the same way the duration of the steps, representing the time separating two consecutive double support instants, is calculated as the time between the projection of two consecutive local vertical maxima of the trajectory of the center of mass. The speed of the gait is calculated, on a straight line, between the moment of appearance of the person in the field of view of the camera and its last moment of appearance. The extraction of the three gait parameters is shown in Figure 3.

IV. EVALUATION OF THE ALGORITHM

A. Description of the experiment

Eleven subjects took part in the experiment among which four women and seven men aged from 22 to 53. The experiment consisted in using an actimetric carpet (GAITRite) of 5 meters length, as a reference for evaluating our system. The error of measurement is 1.27 cm which is the distance between the sensors in the carpet. The experiment was held in a smart room where the sujects were filmed by a camera RGB-D. The field of view of the camera is 3m wide and 6m deep. The persons walked perpendiculary to the camera (except for a situation that we describe in the following paragraph), at a distance of 3m from the camera.

The instruction given to the subjects was to make 4 situations. Each subject had to walk three times in every situation to have a significant amount of data (165 sequences were analysed). The situations were to:

• walk in a straight line "normally";

• walk straight making small steps;

• walk with a long skirt to know if the clothes disturb the gait (for all the subjects);

• walk in front of the camera (the other situations are performed with the camera perpendicular at the subject). The gait parameters obtained with the carpet were compared to the values provided by the camera algorithm. In this experiment we compared three gait parameters: length and duration of the steps and the speed of the gait.

B. Results

1)

Results based on the error on each sequence: The mean of the length, the duration of the steps and the speed of gait are calculated for the carpet and our algorithm for each sequence. A sequence is a walk composed of 2 to 5 steps, depending on the subject and the situation. The error between the mean obtained with the carpet and the camera is calculated. Then the average error is obtained for each situation. The results are shown in Table I.

For the length of steps the average error for situations "Normal", "Small" and "Skirt" (all exept "Front") is less than 4.7%.

For the duration of steps, the average error for all situations (exept "Front") is less than 6.6%.

For the speed of the gait, the average error for situation "Normal" and "Small" is less than 3.69%. For the "Skirt" situation, the average error is more important, 4.37%. The situation with the camera in front of the subject is the worst for the accuracy of the speed.

Only the error on the situation with the camera in front is important. When the camera is in front of the subject the results are degraded. The system gives accurate enough

  ) and 1(b).

Fig. 1 .

 1 Fig.1. Displacement of the center of mass as described in the literature[START_REF] Saunders | The major determinants in normal and pathological gait[END_REF].

Fig. 2 .

 2 Fig. 2. Displacement of the center of mass extracted by our algorithm (before filtering) for a person walking in straight, perpendicularly of the camera.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Extraction of the gait parameters from the center of mass.

Fig. 5 .

 5 Fig. 5. Gait acquisition screen capture.

Inria, Villers-lès-Nancy, 54600, France.

Université de Lorraine, LORIA, UMR 7503, Vandoeuvre-lès-Nancy, 54506, France.

CNRS, LORIA, UMR 7503, Vandoeuvre-lès-Nancy, 54506, France.

ACTIMETRIC CARPET (1: NORMAL STEPS, 2: SMALL STEPS, 3:SITUATION WITH SKIRT, 4: FRONT CAMERA).

values on all the other situations. The speed of gait is the most accurate of the gait parameters. 2) Results based on the mean of each gait parameter: Another way to analyze the results is to calculate the mean of each parameter of the gait. The results are shown in Table II. In this table the mean for the three gait parameters is calculated for the actimetric carpet and for our camera system.

First for the steps length, we can see that the mean obtained with the carpet and with the camera are similar for all the situations. The maximum error between carpet and camera is 2 cm. The situation with the small steps is as accurate as the normal steps.

For the duration of the steps, the mean maximum error is 0.04 s for all the situations.

For the speed of the gait, the error is more important in the "Front" situation than in the other situations. In "Front" situation, the mean error is 6.3 cm/s. In the other situations, the mean error is less than 3.7 cm/s.

In this part, we can see that our system provided values similar to those provided by the carpet in all the situations when we are interested in the mean of the length and duration of the steps and in the speed of the gait.


Florent Margnat

Xavier Gloerfelt

On compressibility assumptions in aeroacoustic integrals: a numerical study within subsonic mixing layers

Keywords: numbers: 43, 28, Ra, 43, 20, Wd

Two assumptions commonly made in predictions based on Lighthill's formalism are investigated: a constant density in the quadrupole expression, and the evaluation of the source quantity from incompressible simulations. Numerical predictions of the acoustic field are conducted in the case of a subsonic spatially evolving two-dimensional mixing layer at Re = 400. Published results of the direct noise computation (DNC) of the flow are use as reference and input for hybrid approaches before the assumptions on density are progressively introduced. Divergence free velocity fields are obtained from an incompressible simulation of the same flow case, exhibiting the same hydrodynamic field as the DNC. Fair comparisons of the hybrid predictions with the reference acoustic field valid both assumptions in the source region for the tested values of the Mach number. However, in the observer region, the inclusion of flow effects in the Lighthill source term is not preserved, which is illustrated through a comparison with the Kirchhoff wave-extrapolation formalism, and with the use of a convected Green function in the integration process.

I. INTRODUCTION

The overall problem area addressed here is that of estimating the noise radiated by unsteady flows. We aim at improving the class of methods in two steps, usually described as hybrid, which use experimentally or numerically generated flow data as an input for solutions of wave propagation equations, in particular those solutions written with an integral formalism such as Lighthill's analogy and the Kirchhoff wave extrapolation method. Because the two-step approach allows cost reduction in the estimation of the radiated field, hybrid methods are frequently chosen for the prediction of flow generated noise in many applications such as jet, high-lift devices and landing-gear on planes, car mirror and window vortex, fans, propellers, etc. That prediction consists in extracting among the unsteady motions in the flow those able to excite outward propagating sound waves. Consequently, once it is obtained, the analysis of the acoustic field and its correlation to the flow events can be conducted in order to identify the noisy vortical motion. In that sense, hybrid methods may also help the physical understanding of aeroacoustic issues.

Any analogy approach relies on the explicit knowledge of the source quantity that excites the wave operator. Since Lighthill's initial derivation 1 of an inhomogeneous wave equation from the equations of motion, with a source term defined by the flow quantities, a constant a) Electronic address: florent.margnat@univ-poitiers.fr research effort has been devoted to improve its potential, from theoretical works on the source term and surrounding flow effect modelling, to numerical works increasing quickness and accuracy of the solution. In particular, Crow 2 investigated analytically the distinct role of the compressibility in the source term and in the wave operator. By the method of matched asymptotic expansion, he concluded that, for low Mach number and compact source region, Lighthill's solution for the density is adequate if the latter is assumed constant in the quadrupole. Because Lighthill's equation is an exact rearrangement of the flow equations, it contains flow effects on acoustics such as wave convection and refraction. The left-hand side being the wave operator without flow, those effects are included in the source side. Consequently, the source term thus obtained is not the true 3 acoustic source, because it accounts for other phenomena in addition to acoustic energy production. Questing for such a reduced source expression, Lilley [START_REF] Lilley | On the noise from jets[END_REF] and Ribner [START_REF] Ribner | Effects of jet flow on jet noise via an extension to the Lighthill model[END_REF] derived an analogy equation for jet noise modelling, through the specification of unidirectional, transversely sheared, mean flow in the Lighthill source term, so that its solution displays sound refraction by velocity gradients. This was numerically applied to a mixing layer by Colonius et al [START_REF] Colonius | Sound generation in a mixing-layer[END_REF] . Recently, Suzuki & Lele 7 derived approximate Green functions for a source in a mixing layer. It should be the appropriate propagator to associate, in the convolution integral, with a source expression that would be purified from convection and refraction effects.

The full knowledge of the source term with fine space and time resolutions brought by the numerical simulations allows to check the validity of assumptions Compressibility assumptions in aeroacoustic integrals 1 made in theoretical or experimental studies. Moreover, building reduced-order models of the flow may reduce the computational cost of closed-loop control systems. In that context, Samanta et al 8 studied the sensitivity of several analogy formulations to errors intentionnally introduced into the source quantity, namely removing the highest modes from the decomposition or perturbating the most energetic ones.

The present contribution builds on these efforts to characterise the content of the Lighthill source term, in order to provide background for current source term evaluation in analogies. The assumptions regarding density are investigated numerically, through the association of a constant density with compressible velocity fields, as well as the evaluation of the quadrupoles from fully incompressible data. We also analyse how those assumptions are sensitive to Mach number variations and how they affect the inclusion of flow effects on acoustics in Lighthill's formalism. For the latter, the convected wave equation and the Kirchhoff formalism are used as reference points.

Numerical evaluations of aeroacoustic integrals are thus conducted in the frequency domain in the case of a subsonic spatially evolving two-dimensional mixing layer at Re = 400. Such flow is free from surface effect such as reflexion or diffraction, thus focusing our analysis on the original source term. Moreover, a two-dimensional configuration allows the integration of many source quantities over many volume extents with an affordable computational effort. The integrands are read from databases generated by the numerical solution of the Navier-Stokes equations. A direct noise computation (DNC) provides the compressible velocity, pressure and density fields, and the reference result for the acoustic field at the same time. In addition, an incompressible simulation provides the source quantity using divergence free velocity field. To the best of the authors' knowledge, the resolution of Lighthill's equation with source data from a DNC and an incompressible simulation of the same flow was not conducted before.

The paper is organised as follows. In section II, Lighthill's and Kirchhoff's formalisms are recalled in the Fourier space, for either a convected wave equation or a propagation in a medium at rest. In section III, the mixing-layer parameters are presented along with the numerical methods applied for the flow simulations and integral evaluations. The coherence between the several prediction strategies is checked in section IV, and the inclusion of the convection effect in the Lighthill source term is illustrated. The influence of compressibility in the evaluation of the quadrupoles is then adressed in section V. Eventually, the main results of the paper are summarised and further discussed in section VI.

II. MATHEMATICAL FORMULATION

A. Kirchhoff's formalism

For a moving medium, the acoustic pressure p ′ at any time t and location x in a volume V is related to the distribution Q of sources within V and the distribution of the pressure and its derivative on the boundary of V , noted Σ, by the generalised Green formula [START_REF] Goldstein | Compressibility assumptions in aeroacoustic integrals 11[END_REF] . For a 2-D configuration, with U ∞ = (U ∞ , 0) in the observer domain, it can be written as

p ′ (x, t) = ∫ ∞ -∞ ∫ ∫ V Q(y, τ ) G(x, t|y, τ )dydτ + ∫ ∞ -∞ ∫ Σ { G ∂p ′ ∂y i -p ′ ∂ G ∂y i } n i dΣ(y)dτ + U ∞ c 2 0 ∫ ∞ -∞ ∫ Σ { p ′ D ∞ G Dτ - G D ∞ p ′ Dτ } n 1 dΣ(y)dτ (1)
where n = (n 1 , n 2 ) is the unit normal vector pointing inward the observer domain V , c 0 is the sound speed of the ambient fluid, D ∞ /Dt = ∂/∂t + U ∞ i ∂/∂x i , and G is the time-domain Green function solution to the uniformly moving medium wave equation. By taking the following Fourier transform:

F [f (x, t)] = f (x, ω) = ∫ ∞ -∞ f (x, t)e -iωt dt ( 2
)
where ω is the angular frequency and i 2 = -1, formula (1) reduces to the form

p ′ (x, ω) = ∫ ∫ V Q(y, ω)G(x|y, ω)dy + ∫ Σ { G(x|y, ω) ∂p ′ (y, ω) ∂y i -p ′ (y, ω) ∂G(x|y, ω) ∂y i } (n i -M 2 ∞ n 1 )dΣ(y) + 2iωM ∞ c 0 ∫ Σ p ′ (y, ω)G(x|y, ω)dΣ(y) (3)
where the time factor exp(-iωt) has been omitted. The Green function for the propagation in a uniform flow is given in the 2-D frequency domain by [START_REF] Lockard | An efficient, two dimensional implementation of the Ffowcs Williams and Hawkings equation[END_REF][START_REF] Gloerfelt | Direct computation of the noise radiated by a subsonic cavity flow and application of integral methods[END_REF][START_REF] Morfey | Green's functions in computational aeroacoustics[END_REF] :

G c (x|y, ω) = i 4β exp ( iM ∞ kr 1 β 2
) H

(2) 0

( kr β β 2 ) (4)
where y = (y 1 , y 2 ) and x = (x 1 , x 2 ) are the source and observer locations respectively, and

r i = x i -y i . H (m) ν
is the Hankel function of order ν and kind m, k = ω/c 0 ,

β 2 = 1 -M 2
∞ is the Prandtl-Glauert factor, and

r β = √ (x 1 -y 1 ) 2 + β 2 (x 2 -y 2 ) 2 .
The Green function for the uniform medium at rest, and the related pressure formula, are recovered when M ∞ is set to 0.

In our configuration sketched in figure 1, the extrapolation of the acoustic waves from the control surface Σ = (y 2 = y s 2 ) ∪ (y 2 = -y s 2 ) is obtained by the following formula:

p ′ (x, ω) = n 2 i 4β ∫ Σ [ ∂p ′ ∂y 2 H (2) 0 ( kr β β 2 ) -p ′ kr β β 2 H (2) 1 ( kr β β 2 )] exp ( iM ∞ kr 1 β 2 ) dy 1 (5)
Kirchhoff integral (5) theoretically predicts the acoustic pressure field provided that no source phenomenon is left Compressibility assumptions in aeroacoustic integrals 2 inside V , that the pressure distribution on Σ contains the acoustic information generated by the sources outside V , and convection is the only flow effect acting inside V .

B. Lighthill's formalism

Lighthill's equation combines the mass and momentum conservation equation to form an inhomogeneous wave equation for the density. Consequently, for an observer located in the acoustic (isentropic) region, the solution for the 2-D acoustic pressure field is given, in the frequency domain, by:

p ′ (x, ω) = i 4 ∫ ∫ D

S(y, ω)H

(2) 0 (kr) dy (6) where D is the source region and r = ||x -y||, for a propagation in an unbounded medium at rest. In (6), S is the Lighthill source term, defined by:

S = -∂ 2 T ij ∂y i ∂y j = -∂ 2 ∂y i ∂y j [ ρu i u j + (p ′ -c 2 0 ρ ′ )δ ij + τ ij ] (7)
where u i are the components of the velocity field u = (u, v), ρ ′ is the density fluctuation from its value in the medium at rest and τ ij stands for the components of the viscous stress tensor τ .

Expression (6) can be transformed taking the space derivatives over the Green function instead of the source term, leading to:

p ′ (x, ω) = -i 4 ∫ ∫ D [ k 2 r i r j r 2 H (2) 2 (kr) -k δij r H (2) 1 (kr) ] T ij (y, ω) dy (8)
In the case of a propagation in a uniform flow at

U ∞ = (U ∞ 1 , U ∞ 2 )
, the derivation of a convected wave equation 11 yields:

p ′ (x, ω) = ∫ ∫ D S c (y, ω)G c (x|y, ω) dy ( 9
)
where

S c = ∂ 2 T c ij ∂y i ∂y j (10) = ∂ 2 [ ρ(u i -U ∞ i )(u j -U ∞ j ) + (p ′ -c 2 0 ρ ′ )δ ij + τ ij ] ∂y i ∂y j
Note that it is not harmless for the physical interpretation that the source expression depends on the propagation medium governing equation. It is indeed an intrinsic consequence of the analogy formalism, as explained by several authors [START_REF] Goldstein | On identifying the true sources of aerodynamic sound[END_REF][START_REF] Samanta | Robustness of acoustic analogies for predicting mixing-layer noise[END_REF][START_REF] Cabana | Identifying the radiating core of Lighthill's source term[END_REF] .

Similarly to the non-convected case, the space derivatives can be analytically transferred to the Green function instead of numerically estimated on the source term. The second-order derivatives of the convected Green function ( 4) are given in [START_REF] Gloerfelt | Direct computation of the noise radiated by a subsonic cavity flow and application of integral methods[END_REF] .

III. FLOW CONFIGURATION AND NUMERICAL TOOLS A. Reference flow and source data

As a reference, the acoustic field radiated by a spatially evolving subsonic mixing layer is obtained by directly solving the compressible Navier-Stokes equations on a spatial domain containing both the mixing region and the acoustic far field. The flow configuration and solver are the same as in ref. [START_REF] Moser | Numerical study of Mach number and ther-mal effects on sound radiation by a mixing layer[END_REF] .

The frame origin is at the inflow on the centerline, and the mean inflow velocity profile, illustrated in figure 1, has a hyperbolic tangent shape

ū(0, y 2 ) = U c + ∆U 2 tanh ( 2y 2 δ ω ) (11)
where

U c = U1+U2 2
and ∆U = U 1 -U 2 , noting U 1 and U 2 the high-and low-speed flow respectively. The vorticity thickness, which varies in the streamwise direction, is defined as

δ ω (y 1 ) = ∆U max y2 ( ∂ ū ∂y 2 ) (12)
Its value at the inflow is noted δ ω and is taken as the reference length throughout the paper. Two subsonic isothermal cases with U 1 = 2U 2 are considered here, with different value of M = ∆U/c 0 : 0.25 and 0.40. The Reynolds number is defined by Re = ρ 0 δ ω ∆U/µ, where ρ 0 is the density of the surrounding fluid (set to unity in those adimensionalized equations) and µ is the dynamic viscosity, and its value is fixed at Re = 400. In order to obtain a periodic and spatially fixed vortex pairing phenomenon, the flow is pertubed just downstream of the inlet, at the frequency f 0 = 0.132U c /(δ ω,f ) (where δ ω,f = 2δ ω is the local vorticity thickness at x ≈ 80), and its first subharmonic f 0 /2, using a divergence-free forcing function with low amplitude, as presented in ref. [START_REF] Moser | Numerical study of Mach number and ther-mal effects on sound radiation by a mixing layer[END_REF] .

The space derivatives are evaluated using sixth-order compact finite difference schemes, while a third-order Compressibility assumptions in aeroacoustic integrals 3 Runge-Kutta scheme provides the time marching. Finally, the boundary conditions are written using a characteristic formulation, and terms are set inflow and outflow. The computational domain lays over 800δ ω in both directions. The grid size is (2071 × 785). It is stretched from the centerline in y 2 direction. In y 1 direction, it is almost constant below 400δ ω , and streched after, combined with a dissipation function to create the sponge zone. The flow solver and numerical parameters are described in details in ref. [START_REF] Moser | Numerical study of Mach number and ther-mal effects on sound radiation by a mixing layer[END_REF] .

These simulations provide both the reference acoustic field and the input data for the Kirchhoff solution (5) and the computation of the Lighthill source quantities ( 7), ( 10), ( 15) and ( 16).

B. Hybrid method implementations

Thanks to the periodic nature of the flow, the computation of the aeroacoustic integrals can be conducted in the frequency domain thus avoiding the integration over the whole history of the time series that is implied by 2D Green function in the time domain . Standard Fast Fourrier Transform routines are applied to the source time series. The latter are recorded about 500 times a period from the direct computation. The Hankel functions are computed using Amos library [START_REF] Amos | A portable package for bessel functions of a complex argument and nonnegative order[END_REF] . A mid-panel quadrature is implemented for the integration of the Lighthill source term, while the Kirchhoff surface integral are evaluated with a trapezoidal rule. Streamwise, the integration domain extends over the whole DNC domain. The observer grid covers the DNC domain and is uniform with 150 points in each direction, leading to ∆xf 0 /c 0 ≈ 8. For the formulations based on a convected wave equation, one sets U ∞ = (U 1 , 0) and M ∞ = U 1 /c 0 for the high-speed flow, and U ∞ = (U 2 , 0) and M ∞ = U 2 /c 0 for the low-speed flow.

The truncation of the source domain, typically at the outflow, can be treated either by a dissipation region, a spatial weighting, or the addition of a residual term based on the Reynolds transport theorem [START_REF] Martínez | Efficient implementation of equivalent sources from CFD data in curle's analogy[END_REF] accounting for the missed region. Here, the two first techniques are used. Indeed, the dissipation function used in the direct computation prevents the paired vortices to generate spurious acoustic waves as they pass outflow, but their damping is not sufficient to avoid a strong contribution to the integral. Thus a weighting function based on a Tukey windowing is applied to the source term (7). With respect to the spatial top-hat window which is de facto applied when integrating (6) over the computational domain, such weighting function reduces the spectral leakage [START_REF] Obrist | The influence of spatial domain truncation on the prediction of acoustic far-fields[END_REF] , that is spatial spurious noise directly translated into temporal noise through the dispersion relation. The source term expression in the convected form (10) appeared more sensitive to truncation effects, that is why a weighting of the following Gaussian form is applied:

W (y 1 ) = exp ( - [ y 1 -y 0 1 σ(y L 1 -y 0 1 ) ] 4 ) (13)
where y 0 1 is the position before which W is set to 1., y L 1 and σ being adjusted by the user for a given computational configuration in order that no radiation is emitted from the outflow. Those weighting functions have a streamwise extent from 500δ ω to the outflow boundary. The numerical singularity in the Hankel function, appearing when the observer is located inside the integration domain, can be treated by removing a small disk around it when computing the integral [START_REF] Schram | A boundary element extension of curle's analogy for non-compact geometries at low-mach numbers[END_REF] , or by setting the observer grid in a staggered way with respect to the source grid. That latter solution suits better when the spatial derivatives are evaluated on the Green function itself.

IV. COHERENCE OF THE PREDICTION STRATEGIES

A. Validation of the hybrid methods

The present bi-harmonically excited mixing layer experiences a periodic vortex-pairing phenomenon, which occurs around y 1 = 200δ ω , and which generates acoustic waves at that pairing frequency. The resulting instantaneous acoustic pressure field, obtained by the DNC, is plotted in figure 2a). Both Kirchhoff's (5) and Lighthill's (9) formalisms, based on a convected Green's function and fed by source data from the DNC, accurately predict that acoustic field, as visible in figure 2b) and 2c). The comparison is fair with the field directly computed solving the compressible Navier-Stokes equations plotted in figure 2a). In particular, the wavefront pattern, with a maximum radiation around θ = 50 o on both layer sides, is well recovered (here and henceforth, θ is the angle counted counterclockwise from the downstream direction and centered on the pairing location).

In figure 2d), the acoustic field is plotted resulting from the numerical evaluation of Lighthill's solution (6) with the Green function for an observer medium at rest integrated over the full DNC domain. The fair comparaison to the field directly computed illustrates well that the Lighthill analogy accounts for the convection effects, yet this is the case here because the source data come from the DNC itself and the integration domain extends over the whole observer region. The integration of a source term estimated from a compressible solver, over a domain including the observer, is then found equivalent to the application of a convected Green function to a localised source region, as firstly described by Goldstein [START_REF] Goldstein | Compressibility assumptions in aeroacoustic integrals 11[END_REF] and numerically illustrated by Bogey et al [START_REF] Bogey | An illustration of the inclusion of sound-flow interactions in Lighthill's equation[END_REF] . As shown in section V, that does not hold any more when a constant density is assumed in the quadrupole.

The mean square value of the fluctuating pressure is plotted in figure 3a) on a circle at r = 150δ ω away from the apparent source location (x 1 = 200δ ω , x 2 = 0). This quantity approaches the acoustic intensity in the far field, and provides here a more quantitative comparison between the four methods. They agree well for all radiation angles. In particular, the level fall for increasing |θ| is perfectly reproduced by the hybrid methods. For |θ| higher than 90 o , some directivity lobes are visible, at very low levels. The method using the integration of the Compressibility assumptions in aeroacoustic integrals 4 source term over the whole domain is the most affected by this artefact, especially on the low-speed side.

In figures 2 and 3a), the control surface used in the Kirchhoff method is defined by y s 2 = 40δ ω , whereas the convected Lighthill source term is integrated over the volume inside this surface. How the predicted acoustic field depends on the control surface position and on the source extent is analysed through figure 3b), where the error with respect to the DNC is plotted, being defined as:

E 1 = ||I Hybrid (x) -I DN C (x)|| ||I DN C (x)|| ( 14
)
where I is the acoustic intensity and ||f (x)|| is the mean value of |f | over a set of observer points. Here, that set is defined as the line x 2 = 300δ ω for 100 ≤ x 1 /δ ω ≤ 500.

For both hybrid methods, the error converges towards a constant value of about 5% for y s 2 > 40δ ω , but the Compressibility assumptions in aeroacoustic integrals 5 behavior is different for lower values. The convected Lighthill computation reaches the steady error level from y s 2 ≈ 10δ ω after a regular decrease from unity. For the Kirchhoff computation part, the error is much higher when the control surface is located very close to the layer, and it catches the Lighthill curve only for y s 2 30δ ω . From this result, the conclusion is that the acoustic source mechanism is confined in |y 2 | ≤ 10δ ω , and that below y 2 ≈ 30δ ω , the fluctuations associated with the vortex dynamics dominate in the pressure field and prevent the Kirchhoff integral from a correct extrapolation of the acoustic waves.

B. Uniform flow effect

A decomposition of the integration domain D is now introduced, in order to analyse convection effect inclusion in the Lighthill source term from a spatial point of view. Using the notations presented in figure 1, the region defined by |y 2 | < y S 2 is noted D SOU RCE and hereafter called 'the source region', while its complement with respect to D is noted D OBS and called 'the observation region'. Hybrid computations are then performed using the Green function for a propagation medium at rest, that is using (5) with M ∞ = 0 for the Kirchhoff method, and (6) with an integration over D SOU RCE for the Lighthill analogy.

The resulting instantaneous acoustic pressure fields for y s 2 = 40δ ω , plotted in figure 4a) and 4b), are very similar between the two methods, but the prediction does not agree with the DNC reference shown in figure 2a). However, if the contribution of the, still not convected, Lighthill source term with an integration over D OBS , shown in figure 4c), is added to both, the correct pattern is recovered in figure 4d) for the Kirchhoff method, and is the same as in figure 2d) for the Lighthill analogy. The lobes for |θ| ≥ 90 o are again visible in figure 4d), suggesting they are an effect from the acoustic/propagation domain (ie: from the boundary conditions in the DNC), and not from the source mechanism, otherwise they should have been present in the convected Kirchhoff prediction in figure 2c).

The convection effect can thus be provided to the wave extrapolation from a Kirchhoff surface by the integration of the Lighthill source term over the observer domain, which plays here exactly the same role as the convected Green function, that is, schematically:

∫ |y2|=y s 2 { G c ∂p ′ ∂n -p ′ ∂G c ∂n } dy 1 ≈ ∫ |y2|=y s 2 { G ∂p ′ ∂n -p ′ ∂G ∂n } dy 1 + ∫ ∫ |y2|>y s 2 S • G dy
Note that the Kirchhoff formalism leads to the same flow quantity in the integrand, namely the pressure and its normal derivative to the control surface, for both the convected and static medium wave equations, unlike the Lighthill formalism, for which the source term expression is dependent on the specific choice of the propagation operator [START_REF] Goldstein | On identifying the true sources of aerodynamic sound[END_REF] . The fully compressible solution was needed for that illustration, which may not improve the efficiency of hybrid methods, though it confirms the robustness of the interpretation of Lighthill's equation as an implicit equation for (aero)acoustic propagation within a uniformly moving flow. Moreover it is of theoretical in-Compressibility assumptions in aeroacoustic integrals 6 terest that, within the integration process, the convected action of the Ligthill source term is spatially limited to the observer region (provided that the Kirchhoff surface accounts for the whole source process). Now that the hybrid approaches have been validated when the source quantity is evaluated from the DNC, whether the density variations should be included in the Lighthill source term to ensure a correct prediction is investigated in the following.

V. DENSITY VARIATIONS IN THE QUADRUPOLES

Attempting to separate a purely radiating part from flow induced effects, Cabana et al. [START_REF] Cabana | Identifying the radiating core of Lighthill's source term[END_REF] suggested a decomposition of the Lighthill source term, featuring density gradient, dilatation and vorticity from the double divergence. Subterms including the two formers were identified as a compressible reaction to subterms identified as driving terms based on vorticity and kinetic energy.

Margnat & Fortuné 20 applied that decomposition to the spatially evolving mixing-layer and noticed the importance of the dilatation transport term to account for the mean flow effect. Moreover, applying the Reynolds decomposition to density and velocity components of the Lighthill tensor, Moser et al. [START_REF] Moser | Numerical study of Mach number and ther-mal effects on sound radiation by a mixing layer[END_REF] reported a significant effect of the term made up with density fluctuations and mean streamwise velocity when the Mach number is increased. Those studies pose the question of the preservation of the source and propagation mechanisms included in the Lighthill source term when the density field is submitted to assumptions, as it may occur in both theoretical works and hybrid predictions.

A. Two assumptions about compressibility

The analogy approach assumes that the flow drives the acoustics with no feedback from the latter, so that Lighthill's equation can be solved explicitely once the source distribution is known. This is theoretically valid for weakly compressible flows, say, low Mach numbers. In the absence of heating, once a low Mach number is assumed, it is tempting to assume a constant density in the evaluation of the source term, thus replacing ρ by ρ 0 in the source expression, that is:

S ≈ ρ 0 ∂ 2 u i u j ∂y i ∂y j = S 0 (15)
Similarly, we define

S c,0 = ρ 0 ∂ 2 ∂y i ∂y j [ (u i -U ∞ i )(u j -U ∞ j ) ] ( 16
)
This assumption is usually made when the source quantity is evaluated from experimental data in flows without significant density variations (subsonic and isothermal flows): the compressible content of the velocity field is preserved, while density gradients are neglected. It is also often assumed in theoretical works, e. g. [START_REF] Huerre | Sound generation by instability waves in a low Mach number jet[END_REF] . For the present mixing layer, figure 5a) shows that neglecting density variations does not alter the acoustic prediction using Lighthill's analogy in the convected form (10). However, for the form (7) that includes the convection effect in the source term, it does, as visible in figure 5b), yielding an incorrect wave pattern.

Note that in acoustic motions, the fluctuating activity involves the pressure and the velocity together with the density, meaning that a constant density assumption will not remove the acoustic part of the Lighthill source term contained in the velocity field. An incompressible simulation, however, will lead to a source quantity evaluated with a divergence-free velocity field, that is:

S c ≈ ρ 0 ∂ 2 ∂y i ∂y j [ (û i -U ∞ i )(û j -U ∞ j ) ] = Ŝc,0 (17)
where ∇ • û = 0. Alternately, Ŝ0 = ρ 0 ∂ 2 ûi ûj ∂y i ∂y j . Such assumption on the source term, which is common in the use of hybrid methods, is investigated in details in the following subsections.

B. Compressible and incompressible near-fields

For the purpose of evaluating (9) under assumption (17), an incompressible simulation of the mixing layer at the same Reynolds number is conducted. The incompressible flow solver is that presented in [START_REF] Laizet | High-order compact schemes for incompressible flows: A simple and efficient method with quasi-spectral accuracy[END_REF] . The perturbation of the inflow profile (11) is only applied to the Compressibility assumptions in aeroacoustic integrals 7 transverse velocity component. The perturbation amplitude at f 0 is 0.0027U 2 , twice that of f 0 /2, while the phase between the two frequencies is 0.4π. These settings are chosen in order to obtain the same longitudinal evolution as obtained in the compressible simulations. The computational domain is limited to 100δ ω on each side of the mixing layer since the acoustic field can not be simulated. The dynamic viscosity is linearly increased on the second half of the computational domain, providing a dissipation of paired vortices before the outflow, in order to prepare the treatment of the truncation in the hybrid acoustic prediction.

The resulting vorticity field is compared in figure 6 to that from the DNC at different Mach numbers. As shown by Moser et al [START_REF] Moser | Numerical study of Mach number and ther-mal effects on sound radiation by a mixing layer[END_REF] there are almost no observable differences between the vorticity field of the simulation at M = 0.25 and that of the simulation at M = 0.40. No more differences appear when the comparison includes the incompressible simulation. In the following, the velocity and pressure statistics in the mixing region are quantitatively compared between the incompressible simulation and the two compressible cases. The Reynolds decomposition is used noting f = f + f ′ where f is the mean value of f over a pairing period.

The spreading of the mixing layer is characterised with the help of the streamwise evolution of a thickness. In addition to the vorticity thickness (12), two other definitions are used. How far the constant velocity profiles must be moved to deliver the same mass flow leads to the displacement thickness:

δ * (y 1 ) = ∫ 0 -∞ ( ū U 2 -1 ) dy 2 + ∫ +∞ 0 ( 1 - ū U 1
) dy 2 (18) Finally, the summation of the half thickness of each stream yields a third definition:

δ 1/2 (y 1 ) = δ 1 (y 1 ) + δ 2 (y 1 ) (19)
where δ 1 and δ 2 are defined by ū (y 1 , δ 1 ) = U 1 -∆U/4 and ū (y 1 , δ 2 ) = U 2 + ∆U/4, respectively.

The evolutions of the three thicknesses, of the streamwise and transverse velocity RM S fluctuations and of the instantaneous pressure fluctuation are plotted along the mixing-layer axis in figure 7. All the quantities show excellent agreement between the three simulations. In particular, the saturation of the first unstable mode occurs at the same position. The vortex pairing causes a global maximum of the thicknesses and the velocity fluctuations, and a second saturation for the pressure fluctuation. The fluctuations from the incompressible simulation have a slower amplitude decrease at the end of the plotted domain, likely because the dissipation region is not designed in the same way. The vorticity thickness exhibits spurious local variations dowstream of the pairing region. Indeed, because it is defined using the maximum slope of the mean velocity profile, its computation may be affected by the presence of more than one slope maximum, which happens after the pairing.

Concerning the transverse decays of fluctuations, it is shown through figure 8 that for |y 2 | < 25δ ω , the fluctuations of both velocity components and of the pressure, once normalized by the low-speed flow velocity, have exactly the same level for the three simulations. That region may then be referred to as hydrodynamic. Above |y 2 | = 25δ ω , the fluctuations from the incompressible simulation go on decaying exponentially, while the acoustic fluctuations start to be dominant in the compressible simulations. The latter depend on the Mach number, since they scale with M 7/2 in 2D as observed in ref. [START_REF] Moser | Numerical study of Mach number and ther-mal effects on sound radiation by a mixing layer[END_REF] .

At this point, it is clear from figures 6, 7, and 8, that the hydrodynamic part of the flow is the same between both compressible simulations and the incompressible one. On that basis, the assumption of an incompressible velocity field can be tested in the Lighthill prediction.

Compressibility assumptions in aeroacoustic integrals 8

C. Results for the convected Lighthill predictions

In figure 9, the directivity of the acoustic intensity at r = 150δ ω obtained using either S c,0 or Ŝc,0 is compared to the convected Lighthill computation using the fully compressible source term S c . The DNC result is shown too, for reference. Three extents of the source region are considered, in order to visualise any phenomenon that could be missed where the hydrodynamic and acoustic parts coexist. A first worthnoting result appears through the intensity fall for |θ| > 50 o , which is well captured by all the Lighthill predictions. Below these angles, the two sides of the mixing layer do not exhibit the same trends. On the high-speed side, the modelling of the source term does not significantly affect the prediction, even for the smallest source extent. The incompressible source term tends to yield a slight overprediction at small |θ|, however. On the low-speed side, a strong overprediction at aft angles is observed for the three hybrid computations when the source region is limited to |y 2 | ≤ 20δ ω . For larger source extents, the constant density modelling yields nearly the same result as the fully compressible source term, while the incompressible assumption still leads to an overprediction at small |θ|, though reduced.

That influence of the source extent suggests a weakly compressible phenomenon in the near-field, which affects the acoustic radiation while being accounted for by the Lighthill source term if evaluated with a compressible velocity field. Regarding this, the pressure fluctuation near-field is plotted in figure 10 for both flow simulations. The DNC makes visible a pattern of two lines of spots where the fluctuation level is strongly lower than around, like sinks, downstream after the pairing and around y 2 = ±30δ ω . The incompressible flow simulation returns only one such sink, in the vicinity of the pairing. Moreover, an asymmetry between the two streams is noticed, the high-speed flow exhibiting more similarities between the two simulations. Such near-fields may influence strongly the propagation at small angles (for which acoustic rays stay longer in the perturbed region), what could explain why any of the convected Lighthill computations miss the acoustic level there when the source domain do not include it, and why the incompressible source modelling misses this in the low-speed region even for a larger source domain.

Compressibility assumptions in aeroacoustic integrals 9

D. Increasing the Mach number

As investigated by Moser et al [START_REF] Moser | Numerical study of Mach number and ther-mal effects on sound radiation by a mixing layer[END_REF] , increasing the Mach number of mixing layers with U 1 = 2U 2 leads to a reinforcement of the directivity peak around θ = 50 o especially for the radiation towards the high-speed flow. This is illustrated in figure 11a) and figure 11c) with the instantaneous acoustic pressure field as given by the DNC for M = 0.25 and M = 0.4 respectively. That reinforcement is perfectly captured by the incompress-ible source term modelling, as qualitatively visible in figure 11b) and figure 11d), and quantitavely confirmed by the intensity plots at r = 150δ ω in figure 12 showing well the intensity fall after θ ≈ 50 o . In spite of the aforementioned flaw at aft angles, such a successful association of the incompressible source modelling with the convected Green function, even for the highly subsonic case at (M 1 = 0.8, M 2 = 0.4), is worth to notice.

Furthermore, it can be deduced that the effect of the Mach number on the directivity is brought by the acous-Compressibility assumptions in aeroacoustic integrals 10 tic wave convection only, because the three vortical fields do not show any significant difference while the convolution to the convected Green function yields that directivity evolution. This is consistent with Ffowcs Williams' theoretical work [START_REF] Williams | Sound production at the edge of a steady flow[END_REF] , according to which there may be no refraction effect of shear at the interface between two uniform flows of different speed and density, in the limit of a small wavelength with respect to the transverse extent of the flow and large with respect to the interface thickness. The latter conditions are satisfied in the present mixing layers, where the wavelength is about 50 -80δ ω and the transverse extent is indeed infinite due to the boundary conditions (in any case, the computational domain is 400δ ω on each side of the mixing-layer).

VI. CONCLUDING REMARKS

Based on several expressions of the source term, comprehensive hybrid predictions were carried out in order to bring facts about which content of the Lighthill source quantity should be taken into account, regarding the inclusion of density fluctuations and the flow effects on acoustics. The mixing-layer flow case was selected in order to focus on free shear flows while avoiding wall effects. This apparently simple 2D case still requires careful numerical implemention. Instantaneous pressure fields were preferred to point spectra in order to visualise wavefront patterns.

Strong background is provided about the modelling of the quadrupole term, the principal conclusions being:

• The constant density assumption is numerically validated.

• The use of incompressible velocity fields in Lighthill source term is numerically validated, except for the radiation at low angles.

• The validity of both assumptions is limited, however, to the source region, because they cancel out the inclusion of convection effects in the Lighthill source term. Those effects are correctly predicted when the source quantity is fully compressible (e. g. from DNC) and when a convected Green function is associated to either constant density or incompressible quadrupoles.

The two following peripheral points are also emphasized:

• While the hydrodynamic pressure field from the incompressible simulation matches perfectly with those from the DNC, it yields an incorrect acoustic prediction when put into the Kirchhoff formalism. This is an expected result when the control surface is far from the vortical region, but it shows that the pressure field generated by the hydrodynamic waves is not appropriate as such to model the aeroacoustic excitation by the (2D, low-Re) mixing-layer.

• Finally, because the reinforcement of the directivity around |θ| ≈ 50 o due to the Mach number increase is well captured with incompressible sources and a convected Green function, it is concluded that such effect of the Mach number is not refraction in the sheared region but convection in the observer region.

FIG. 1 .

 1 FIG. 1. Configuration and notations for the application of aeroacoustic integrals to the mixing layer flow.

FIG. 2 .FIG. 3

 23 FIG. 2. Instantaneous acoustic pressure field of the mixing layer, predicted by a) DNC b) Lighthill's analogy with convected Green's function (9), c) the convected Kirchhoff method, and d) Lighthill's analogy (6) with the full source domain. Levels are the same for each field, linearly increasing from -4.0 × 10 -5 c 2 0 (black) to 4.0 × 10 -5 c 2 0 (white).

FIG. 4 .

 4 FIG. 4. Hybrid computations with the non-convected Green function. a) Lighthill's analogy (6) integrated over DSOURCE: |y2/δω| < 40; b) Kirchhoff's method (5) with M∞ = 0 for y s 2 = 40δω; c) Lighthill's analogy (6) integrated over DOBS: 40 < |y2/δω| < 400; and d) sum of b)and c). The sum of a) and c) is identical as figure2d). Levels are the same for each field, linearly increasing from -4.0 × 10 -5 c 2 0 (black) to 4.0 × 10 -5 c 2 0 (white).

FIG. 5 .

 5 FIG. 5. Hybrid predictions based on Lighthill's analogy with a constant density in source subterms. The acoustic pressure field radiated by a) Sc,0 with the convected Green function; b) S0 integrated over the full DNC domain. Levels are the same for each field, linearly increasing from -4.0 × 10 -5 c 2 0 (black) to 4.0 × 10 -5 c 2 0 (white). DSOURCE and DOBS are defined as inside and outside |y s 2 | = 40δω respectively.

FIG. 6 .

 6 FIG. 6. Vorticity snapshots for DNC's at M = 0.40 (top) and M = 0.25 (middle) and for the incompressible simulation (bottom). Levels are from -0.8U2/δω (black) to 0.0 (white).

FIG. 7 .

 7 FIG. 7. Evolution of statistic quantities along the mixing-layer axis. Top: mixing-layer thicknesses; middle: fluctuating longitudinal and transverse velocities; bottom: instantaneous pressure fluctuation and its extreme values in time. M = 0.00 stands for the incompressible simulation.

2 FIG. 8 .FIG. 9 .

 289 FIG.8. Transverse evolution of the velocity and pressure fluctuations at x/δω = 300. M = 0.00 stands for the incompressible simulation.

FIG. 10 .

 10 FIG. 10. Isocontours of the RM S pressure fluctuation in the near-field of the compressible (top) and incompressible (bottom) mixing layer, for M = 0.25. The level of a contour is twice the previous one, starting from 0.0032U 2 2 .

FIG. 11 .FIG. 12 .

 1112 FIG. 11. Instantaneous acoustic pressure field of the mixing layer predicted by DNC (a) and c)) and convected Lighthill's analogy with incompressible source data and y s 2 = 40δω (b) and d)), for M = 0.25 (a) and b), levels from -4.0 × 10 -5 c 2 0 (black) to 4.0 × 10 -5 c 2 0 (white) ) and M = 0.4 (c) and d), levels are five times larger).

Compressibility assumptions in aeroacoustic integrals 12

Acknowledgments

The authors gratefully acknowledge Eric Lamballais and Veronique Fortuné for courtesy share of the flow solvers and compressible data.