domain: spi phys

text:
Design methodology of force feedback laws for active side stick interface*

Gemma Prieto Aguilar 1 , Laurent Binet 2 , and Thomas Rakotomamonjy 3

Abstract-SAFRAN Electronics & Defense and the Information and Systems Processing Department (DTIS) of ONERA have begun a cooperation to evaluate the interest and the methods of use of Active Side Stick Units (ASSU) to improve the safety and flight qualities of helicopters. This paper describes the work carried out to model an environment for simulation and evaluation of haptic feedback laws. An experiment, implemented in the simulator PycsHel at ONERA Salon de Provence, has brought some insight about the influence of ASSU's parameters on the detection of specific haptic feedbacks (Softstops). The results obtained will be added to the simulation model in order to allow the specification, optimal if possible, of the haptic cues. In the early days of aviation, aircraft control was based on the use of mechanical linkages between the flight control surfaces and the pilot's commands. The development of civil and military aviation, and the emergence of increasingly larger, faster and more agile aircraft, led to greater efforts on commands and the need of assistance systems. This is when hydromechanical controls appear: the mechanical linkages are now connected to actuators to move the different control surfaces. This system represented a cost of maintenance too important for the civil aviation, what made the mechanical linkages to be replaced by electrical wires, and the actuators by servo-motors. Nowadays this technology, known as "flyby-wire", is used on the most popular commercial transport aircrafts.

Aircraft manufacturers have followed different trends concerning the pilot's commands, offering each one different advantages. SAFRAN Electronics & Defense (E&D) is currently working on the maturation and development of a new thomas.Rakotomamonjy@onera.fr concept of pilot's command: the Active Side Stick Unit (ASSU). This new technology offers:

• A better ergonomics, offering a clear view to flight displays. • An ability to restore static forces, lost with the transition to the "fly-by-wire" commands; and generate dynamically different haptic feedbacks. These haptic cues can, for example, prevent the approach to the pilot of a critical flight parameter like the stall.

The advantages using ASSUs and dedicated haptic cues are numerous:

• A pilot workload reduction and situational awareness improvement.

• An improvement of the flight envelope safety, • Better performances of the aircraft, since the pilot can apply instructions without hesitation • A coupling of the pilot's and copilot's commands. The ASSUs offer the possibility of generating forces in the grip which can be felt by the pilot. These forces can be adjusted to vary with time, angular position of the grip, aircraft state variables, aircraft/helicopter limitations and other parameters that are related to the flight envelope security. This set of forces defines the static characteristics of the ASSU, and can be decomposed by a combination of softstops (SS), detents, gates, friction, vibrations and the gradient of the nominal force-displacement law, or QF (see Fig. 1).

We also speak about dynamic parameters to refer to the damping ratio and the response frequency of the ASSU (as explained later, the ASSU emulates the behavior of "classical" sticks or yokes, so behaves as a physical, linear second order system mass-spring-damper).

SAFRAN E&D wishes to highlight the performance of its active stick with the demonstration of the capabilities of this new haptic feedback technology. Thus, in the framework of a PhD thesis, the Information and Systems Processing Department (DTIS) of ONERA-Salon de Provence and SAFRAN E&D have started a cooperation to evaluate the interest and the different possibilities offered by the mini active sticks in order to improve the safety and flying qualities of rotary and fixed wing aircraft.

More specifically:

• What are the different fields in which a mini-stick can offer piloting assistance and protection of the flight envelope? • How to define, ab initio, or in an optimal way, the haptic cues and integrate them in the control loop? A state of the art has been done to understand the advances made on this subject, and to evaluate the different possibilities to answer these questions. The main topics which can be found are:

1) ASSU's Parameter setting. Part of the literature dedicated to ASSUs has focused on studying the effects of the characteristics of the stick on the flight qualities of the aircraft. For example, [START_REF] Lusardi | In Flight Evaluation of Active Inceptor Force-Feel[END_REF], [START_REF] Von Grnhagen | Handling Qualities Studies Into the Interaction Between Active Sidestick Parameters and Helicopter Response Types[END_REF], and [START_REF] Greenfield | Side-Stick force-feel parametric study of a cargo-class helicopter[END_REF] studied the variation of the flight quality scales (or HQR, Handling Qualities Rating) as a function of the frequency, damping, inertia or angular displacement of the stick. The aim of these studies was to optimize the ASSU's parameters according to the control law used in the helicopter. 2) Development of Haptic cue functions. The literature reports a large number of articles developing safety flight envelope and navigational aid functions ([4], [START_REF] Einthoven | Tactile Cueing Experiments with a Three Axis Sidestick[END_REF], and [START_REF] Kashawlic | Carefree Maneuvering via Tactile Cueing Flight Demonstrations[END_REF]) These studies often focused on the model to calculate the critical parameter, and not on the type of haptic cue to provide to the pilot. 3) Loop modeling: Pilot + HC + ASSU. Other studies have focused on the modeling of the whole piloting loop, incorporating a helicopter flight mechanics model, a representation of the pilot behavior for a given task, and a model of an ASSU such as [START_REF] Abildgaard | Modeling of a tactile cueing function[END_REF] and [START_REF] Malpica | Handling Qualities analysis of active inceptor force-feel characteristics[END_REF]. This approach allows to study the influence of the stick's parameters on the helicopter's flying qualities. So far, to our knowledge, there is a lack of formal methods for defining haptic cues, other than simulator experimentation with pilots. Thus, we will seek through this thesis to:

1) Develop tools and define criteria that will allow the specification of optimal force feedback laws.

2) The modeling of a complete simulation loop to evaluate the haptic cues defined from these criteria. It is expected that this approach will help the development and testing of haptic cues, by reducing the number of simulator trials. The following sections focus on the modeling of the piloting simulation loop, and in particular, the work done to develop a pilot activity model which takes into account haptic cues. The different criteria to allow the specification of these haptic cues will be addressed in future communications. Results obtained from an experiment, summarizing trends of the haptic detection of pilots are provided, as well as a detection model included in the activity pilot model.

II. CUEING PROTECTION FUNCTION

The first cueing function chosen for the modeling of the simulation loop was the limitation of bank angle to cue the pilot to limit the angle of inclination φ to a maximum of 30˚. The helicopter model is "augmented" with a RCAH (Rate Command Attitude Hold) control law. Thus, the stick positions directly control the angular speeds. Automatic turn coordination (in order to cancel sideslip) has also been added.

III. SIMULATION LOOP

In order to represent the real dynamics of the complete system, each element has been integrated as a module in a SIMULINK model as schematized on Fig. 2. The helicopter model manages four control axes (roll, pitch, yaw, and collective) and several flight aids (ATT, SAS, Vx/Vz hold) can be selected, in addition to the RCAH law.

The model includes, in the feedback loop to the stick, a module (defined in [START_REF] Jeram | Open design for helicopter active control systems[END_REF]) which allows the computation of the parameter to be limited, and its conversion into flight command. Thus, equation (1) provides the maximum pilot command in roll δDDL o before exceeding φ max ,

δDDL o = 2c∆φ φ max ( 1
)
where c is is the SS's return speed to neutral position, and ∆φ the difference between φ max and the current φ.

The Haptic Module provides the force feedback law to be generated on the stick and its definition remains one of the main objectives of this thesis.

A. Pilot Model

A precision pilot model, described by McRuer [START_REF] Mcruer | Mathematical Models of Human Pilot Behavior[END_REF], has been integrated. It provides a list of some aspects of human behavior Y p when controlling an element of an aircraft Y c . This model adds the neuromuscular dynamics of the pilot to the well-known crossover model [START_REF] Von Grnhagen | Handling Qualities Studies Into the Interaction Between Active Sidestick Parameters and Helicopter Response Types[END_REF],

Y p Y c = ω c e -τ s s (2)
where Y p and Y c are the Pilot an Aircraft transfer functions, ω c the crossover frequency, and τ the transport delay time caused by the pilot neuromuscular system.

B. HC Model

Non-linear rotorcraft modelling and simulation is provided by the full non-linear flight mechanics code Flightlab (developed by Advanced Rotorcraft Technology).

C. ASSU Model

The first objective of an active side-stick is to reproduce the behavior of "classical" sticks or yokes. A mass-springdamper system is then generally accepted as a model of the ASSU. Therefore, the active stick can be modeled by a forceinput position-output system of second order for each one of its axes (3):

x F pilot = 1 K ω 2 n s 2 + 2ξω n s + ω n (3)
where ω n /2π is the model frequency, K is the spring's stiffness of the system, and ξ is the damping.

Other elements have been also integrated in the simulation, such as a guidance module, which allows the transmission of pilot instructions in terms of flight parameters (for instance, hold an inclination angle or a forward speed). This module will be improved, enabling more complex piloting tasks.

IV. SENSIBILITY EXPERIMENT

In our study case, the classical approach for building a pilot model could be to measure the error between the current roll angle and the targeted one and to apply corrections (in terms of lateral stick position) to minimize this error. Delays, transfer functions, etc. could be added and adapted to fit to a realistic pilot behavior. This approach would be quite similar to classical and very well-known control methods. The goal being to be able to develop, ab initio, an haptic feedback, our pilot model has to "feel" the variation of the ASSU behavior due to the presence of haptic cues such as force gradients, SS, detent, etc. Moreover, it appeared several questions concerning the cue ability of pilots:

• How do ASSU and cue feedback parameters affect detection? What are their values? • What is the probability that the pilot detect the cue feedback as a function of these parameters? • What are the force thresholds [N] leading to the haptic feedback detection? For this reason, an experiment was designed to analyze the detection of a cue feedback of type SS. It was chosen to test two populations: non-expert subjects and student pilots from the French Air Force Academy at Salon de Provence. The objective was not to compare the results obtained between these populations, but rather to allow the setting up of the experiment with non-pilot subjects and its improvement with more experienced subjects. As the subjects were not helicopter pilots, a simplified linear model was used for the lateral dynamics. The forward speed was kept constant while the vertical motion highly damped. This experiment took place in the prototyping and simulation environment PycsHel at Salon de Provence ONERA center, which enables real-time piloted simulation scenarios. The simulator has two ASSUs equipping the cyclic and collective controls of the pilot.

In the experiment, the subject had to sweep the grip alternatively between right and left positions, and pressed the trigger each time he felt a SS. In addition, a work task had been added forcing the subject to follow a target animated of a random movement along the vertical axis (to be followed with the collective stick). A total of 6 sets of SS combinations were tested in which the frequency or damping of the ASSU model is modified (Table I).

In each set, 3 levels of QF, 4 SS positions and 3 amplitudes were tried (Table II). Finally, we also asked the subjects to change the sweeping velocity between two conditions, namely "slow" and "fast" (both being self-appreciated).

A total of 6 pilots performed all the sets with a previously learning phase. Each one of the sets consists of 144 combinations of SS randomly distributed and repeated 3 times. Each subject started by a different set number. A pause took place at the end of each set, during which the subjects were asked to complete a questionnaire about their sensibility and detection based on the different parameters tested.

V. RESULTS

In order to develop a pilot model based on the detection of effort feedback, statistical analyzes were carried out on the results of the experiment. These analyzes provided:

• Trends about the detection during the modification of the ASSU or haptic feedback parameters, as well as, • A multiple regression model of the force/displacement applied by the pilot on the stick at the time of detection of the SS, as a function of ASSU and haptic feedback parameters.

According to our research, there have been no scientific studies to understand the effects of ASSU's and haptic parameters on human detection. For this reason, an exploratory analysis of variance (ANOVA) of the mean detection per person and combination was carried out.

The analysis results and pilot's comments collected on the questionnaires are presented below:

1) No effect of the ASSU's model frequency was perceived on the SS detection. However, a correct combination of Frequency and Damping is necessary so as not to generate resonance oscillations on the grip. In fact, to the question "Does this combination generates any oscillation at the starting point of the SS?" most of the subjects affirmed to have felt them during the set 1 and especially in the set 2. Some others pointed out that felt vibrations all during the experiment in the rest of the sets, which is certainly due to a slight force ripple/cogging. 2) We found a significant effect of the interaction between the QF and the Amplitudes of the SS. The Fig. 3 shows mean detection for each amplitude of SS and QF tested. We observe how QF augmentation decrease the mean detection in particular for 9 N amplitude (respectively Mn 0.78, Mn 0.66, Mn 0.61), 6 N (respectively Mn. 0.61, Mn 0.55, Mn 0.45) and to a lesser extend 3 N (respectively Mn. 0.26, Mn. 0.22, Mn. 0.20). The subjects have been asked to rank and comment the different values of QF tried in the context of the detection and pilot task. Most of them (4/6) preferred the low-value QF and highlighted the ease to detect at QF 0.5, including the extremely-placed SS (close to the mechanical stop). Only one subject chose to rank first the QF 1.25, and the other admitted that "I feel less [SoftStops] with a QF 1.25. [...] A QF 0.5 makes me to detect less near the neutral point because the grip becomes really sensitive. In the other cases, it helps to improve detection. [...] The best concession is, in my opinion, QF 1".

3) The ANOVA also signaled a significant effect of the interaction Damping-QF. In Fig. 4, we observe again how the increasing value of the QF drops the mean detection for the different damping factors. In the case of the damping, the low-value 0.5 implied better detections, excepting for the QF 0.5.

The student pilots noted their preference concerning the different sets tested. The set 4 (Damping 0.75 Frequency 5) resulted to be the best rated, the set 1, 2 and 3 were all rated equally, followed by the set 5 and finally the number 6. The last two were said to be tiring and difficult to allow SS detection. One of the military students commented: "We tend to feel more SS in sets 1 and 2, because we had to exert less force on the grip and we remain sensitive to the variations of force". He also mentioned what seems to explain why the combination QF 0.5 -Damping 0.5 was less detected/ "However these sets [1 and 2] do not allow, in my opinion, to fly a helicopter because the ASSU is too sensitive. With a QF 0.5, we often arrive at the mechanical stop without any effort". 4) Finally, two significant interactions of the piloting speed of the grip appeared on the ANOVA. The first one, QF-Speed interaction, revealed that the different QF combinations were differently affected by the ASSU's speed, which appears logical. The mean detection decreased of ∆Mn 0.27, 0.22 and 0.18 for the QF values of 0.5, 1 and 1.25, respectively. In the second interaction, Position-Amplitude-Speed, appears that outward SS's positions (i.e. movements to the right) were slightly better detected and less affected by the piloting speed. Piloting speed also affected more to the positions close to the near point, and much less to positions -14/14˚.

B. Regression Model

We performed a multiple linear regression to predict the force applied by the pilot at the moment of the SS detection, as a function of the different haptic and ASSU's parameters.

So far no representative regression model based on the force has been obtained, which could led to the conclusion that the force detection as a function of the different ASSU's parameters do not follow a linear law. On Fig. 5 we can observe how the response of the stick to a command differs from the linear static law. This dynamic behavior being more significant as the force introduced by the pilot is important. In addition, it can be seen that the stick is unable to stop its motion while traversing the SS. This undoubtedly affects haptic detection, since the pilots could detect the force gradient variation at different points, or even "fly over" the SS without noticing it. This behavior is not noticeable with the SAFRAN's ASSU, and future experiments led with this system should provide more repeatable and consistent results.

In parallel with this work, a regression model predicting the overtaking reached from the beginning point of the SS (∆DDL) has been performed. This model explains 43% of the variability of the data obtained from the simulator (Table III 1 ). Since the somatosensory system depends on a large number of factors (human physiology, experience, etc.), a probabilistic model of detection has been added (Table IV 1 ). Despite a poor compliance with the measures, this model allowed us to set up the detection logic in the pilot activity model.

Alternatively, a simple model based on an interpolation of the mean force detection has been performed. This interpolation (Table V) shows some of the trends already mentioned,

• Greater overruns on positions close to the neutral (± 3˚) and a tendency to stabilize the detection force average for the other positions. • Greater efforts with increasing SS amplitude. As shown in Table V, providing mean force detection for a given set of ASSU's parameters (QF, Frequency and Damping), a multi-dimensional interpolation will be used for the other ASSU's parameters, and used as force detection thresholds in the pilot activity model.

VI. RESULTATS ANALYSIS

Based on the results obtained during the experiment, and with the objective of improving the detection of haptic cues, we could consider a nonlinear nominal law. In Fig. 6, the different parameters of the ASSU have been adapted according to the conclusions exposed above. Thus, the force nominal law, as a function of ASSU's displacement, and adapted as to improve pilot detection could be,

• The Nearest SS positions should be protected with larger amplitudes. On inward movements (pronotation), should be even larger since we are stronger. • In the same way, the QF should be stronger on inward movements than on outward. • Decrease of the QF near the SS position. This would reduce arm tension and improve detection. Additionally, to make sure that pilots detect the haptic cue, the damping ratio could be increased when the pilot is operating near or in the SS. Due to the technical characteristics of the sticks used at Onera, this technic (proposed in [START_REF] Einthoven | The HACT Vertical Controller[END_REF]), could not be tested.

VII. CONCLUSIONS

A complete simulation loop is in development, which should allow the definition of haptic feedback laws for specific mission or piloting aid function. The definition of the haptic feedback should rely on criteria which remain to be defined, but which could be a composition of several ones such as the respect of critical parameter limitation, or some flight parameters, or the optimal task completion.

An experiment has been conducted in the PycsHel simulator to initiate the development of a pilot activity model able to detect the different force feedback generated by an ASSU. Statistical methods, adapted to the analysis of the ASSU's parameters, have been developed and could be used for the analyses of future experiments.

An analysis of the results obtained in simulator has highlighted the importance of adapting the different ASSU's parameters to facilitate the detection of haptic cues.

According to these results, it can be concluded that, 1) Pilots are sensitive to force gradient variations, since subjects have shown better detections with little gradients of the nominal law (QF=0.5). Little QF values have also shown to be more sensitive to high speeds, leading to a higher number of overtaking. A good compromise could be a QF value of 1. 2) Pilots are sensitive to ASSU's damping. Sets tested with high damping values have proven less SS's detection because it implies the arm's muscles to be contracted. On the other hand, higher damping values improve precision during helicopter's command and avoid overtaking the SS during high speeds. A good compromise could be a damping value of 0.75.

3) SS's positions influence on their detection. SSs positioned outwards have been more detected than SSs placed inwards. This could be explained physiologically by the well-known fact that we are stronger on inward (pronation) movements, and thereby less sensitive. Subjects also avoided confusing furthest SS with the mechanical stop. At furthest positions, the subject needs to apply a highest force, which implies a tension on the arm's muscles and a worse sensibility to force gradients. Additionally, SSs placed near the trim position presented fewer detections, which can be explained by the second order dynamics of the system. In fact, during rapid force inputs the response of the system deviates from the linear static law. Moreover, a breakout force (to avoid any stick displacement due to small/unintentional applied force) is generally placed at the stick trim position, needing an additional force to initiate the stick motion and which could reduce the SS detection.

VIII. NEXT STEPS The next steps will focus on the improvement of the pilot's detection models by means of the ASSU developed by SAFRAN E&D. Once the model of the activity of the pilot will have a representative response, depending on the haptic parameters and the requirements of the piloting task, the next steps to be done will be:

• The selection of an application case using the lateral cue feedback protection. These tasks should meet certain performance criteria and could be extracted from the ADS-33 standard specification, or from a piloting complex task. • The integration of this task into the haptic loop simulation, and the optimization of the haptic feedback according to piloting and sensitive performance criteria that should be set up. • The implementation of the selected case on the simulator, the analysis of the results and, if needed, the adaptation and improvement of the simulation loop. • Ideally, the implementation of the same procedure for a haptic feedback law on the collective axis.

Fig. 1 .

 1 Fig. 1. Static law showing different type of force feedbacks.

Fig. 2 .

 2 Fig. 2. Haptic feedback loop including an activity pilot and helicopter models.

Fig. 3 .

 3 Fig. 3. Mean detection as a function of the Amplitude of the SoftStop and the gradient of the nominal force-displacement law (QF).

Fig. 4 .

 4 Fig. 4. Mean detection as a function of the ASSU's Damping and the gradient of the nominal force-displacement law (QF).

Fig. 5 .

 5 Fig. 5. Overshoot generated, in comparaison with the linear law, in response to a command with the ASSU.

Fig. 6 .

 6 Fig.6. The linear nominal law (in black) may be not optimally adapted to SS detection. A nominal law adapted (in red) depending on the SS position is here proposed.

TABLE I A

 I SUMMARY TABLE WITH DAMPING AN FREQUENCY COMBINATIONS.

	N˚Set Frequency [Hz] Damping
	1	3	0.5
	2	5	0.5
	3	3	0.75
	4	5	0.75
	5	3	1.25
	6	5	1.25

TABLE II A

 II SUMMARYTABLE WITH VALUES OF QF, SPEED OF THE GREEP, AMPLITUDE AND POSITION OF SS.

	Parameter	Value
	QF	[0.5; 1; 1.25](N. deg -1 )
	Grip speed	[Slow; Fast]
	Amplitude of SS	[3; 6; 9] (N)
	SS Position	[± 3, ± 7; ± 10; ± 14] (deg)

TABLE III A

 III SUMMARYTABLE OF THE REGRESSION MODEL PREDICTING THE OVERTAKING SINCE THE SS'S POSITION .

		∆R 2	B	SE B	β	P
	Model δDDL	0.43				<.001
	Constant		4.12	0.30		<.001
	Damping		0.5	0.11	0.08	<.001
	Amplitude		-0.53	0.04	-0.63 <.001
	SS's Position		0.11	0.01	0.53	<.001
	QF		-1.36	0.31	-0.22 <.001
	Amplitude*QF		0.23	0.05	0.04	<.001

TABLE IV A

 IV SUMMARY TABLE OF THE REGRESSION MODEL PREDICTING THE PROBABILITY OF DETECTING A COMBINAISON.

		∆R 2	B	SE B	β	P
	Model Prob	0.31				<.001
	Constant		0.18	0.04		<.001
	QF		-0.09	0.04	-0.07	0.028
	Amplitude		0.11	0.01	0.63	<.001
	SS's Position		-0.002	0.00	-0.06 <.001
	Speed		-0.26	0.04	-0.31 <.001
	Amplitude*QF		-0.02	0.01	-0.02 <.001
	Speed*QF		0.13	0.03	0.75	<.001
	Amplitude*QF		-0.01	0.01	-0.31 <.001

TABLE V MEAN

 V FORCE TRENDS AS A FUNCTION OF AMPLITUDE AND POSITION OF SOFTSTOP AND FOR DAMPING, FREQUENCY, QF AND SPEED SET.

∆R

explains the proportion of variance explained by the model, B, or b-values, are the estimates coefficients of the linear regression model, SE B are the associated standard error of the b-values, β are the standardized betas and P are the significance associated with the t-test for each b-value of our model.


domain: spi phys math info

text:
Camille Palmier

Karim Dahia

Nicolas Merlinge

Pierre Del Moral

Dann Laneuville

Christian Musso

Pierre Del Moral

Adaptive Approximate Bayesian Computational Particle Filters for Underwater Terrain Aided Navigation

published or not. The documents may come

I. INTRODUCTION

The ability of an underwater vehicle (UV) to accomplish a mission depends on the performance of on-board navigation algorithms. This paper considers Terrain Aided Navigation (TAN) which can provide a drift-free navigation tool for UV, yielding a powerful alternative to current navigation methods which include resurfacing for GPS [START_REF] Melo | Survey on advances on terrain based navigation for autonomous underwater vehicles[END_REF]. TAN generates vehicle position estimates by correlating terrain measurements obtained by a multi-beam sensor together with stored terrain maps, which is related to a filtering problem.

The aim of filtering is to estimate the state of an evolving system, customarily modeled by a stochastic process and called the state process. The state process cannot be measured directly but only via a related process named the observation process. The filtering problem consists in computing the posterior density of the state at the current time given the observation data accumulated up to that time. This state-space problem can be solved analytically when the propagation of the state process and the observation process are both linear and when the noises are independent white Gaussian. In this case, the Kalman Filter (KF) is known to be optimal. The KF was extended to non-linear models, but may become highly unstable in case of non-Gaussian processes or severe non-linearities such as multimodalities, i.e. multiple maxima, also called modes, in the posterior density. This paper considers higly non-linear and non-Gaussian observation processes. Indeed, the terrain profile generates severely non-linear and ambiguous measurements, which yields multimodalities. The Kalman-like filters cannot be used for this class of problems. To tackle non-linearities such as multimodalities, particle filtering methods were introduced [START_REF] Neil J Gordon | Novel approach to nonlinear/non-gaussian bayesian state estimation[END_REF]. Particle filters approximate the posterior density by a mixture of weighted Dirac functions called particles.

The present work addresses the case where the noise of the observation process (also called measurement noise) is unknown. The observation model is the result of a numerical black box model, therefore the likelihood is unknown. Indeed, in the case of underwater TAN using a multi-beam bathymetric sensor, the observation equation is implicit and requires a numerical approximation in practice. The stochastic characteristics of the error between the actual measurements and the numerical approximation of the observation cannot be inferred.

To address cases where the measurement noise distribution is unknown (therefore the likelihood is unknown), Approximate Bayesian Computational (ABC) methods were introduced ([2], [START_REF] Ristic | Bayesian likelihood-free localisation of a biochemical source using multiple dispersion models[END_REF], [START_REF] Marin | Approximate bayesian computational methods[END_REF]). ABC is a class of estimation methods used to overcome the lack of knowledge about the measurement distribution. This method bypasses the evaluation of the likelihood function. An acceptance/rejection test on a uniform law is applied to the samples state. When a candidate sample state is rejected, it is drawn again until it is accepted. This algorithm is computationally demanding since it leads to a rise of the number of samples. To reduce the computational load of ABC methods, a filtering approach was derived -called ABC filter. ABC filter replaces the measurement likelihood by an analytical density functioncalled kernel -whose bandwidth is adaptive [START_REF] Marin | Approximate bayesian computational methods[END_REF]. This scale parameter aims to account for every state density's mode while the measurement ambiguity is not removed. However, it is difficult to tune, which may yield filter unstability.

We propose a new way of choosing this setting, based on the efficiency degeneracy criterion [START_REF] Kong | Sequential imputations and bayesian missing data problems[END_REF]. This criterion is commonly used in particle filters to trigger the resampling step. The proposed method is called Adaptive Approximate Bayesian Computation (A2BC) and can be integrated in any particle-like filter. Two filters have been selected in this work: Regularized Particle Filter (RPF [START_REF] Musso | Improving regularised particle filters[END_REF]) and Rao-Blackwellized Particle Filter (RBPF [START_REF] Schon | Marginalized particle filters for mixed linear/nonlinear state-space models[END_REF]).

The general formulation of the estimation problem, Particle Filter and ABC filter approaches are recalled in Section II. In Section III, we describe the principle of the A2BC method and its implementation within conventional particle filters. In Section IV, we present the problem of submarine TAN navigation. After describing the performance criteria, we numerically compare the performance of two A2BC filters with these of conventional particle filters. Finally, we will conclude in Section V.

II. PROBLEM STATEMENT

A. Filtering problem and particle filter

Consider the following discrete-time state-space model with hidden states {X k } k∈N and observations {Y k } k∈N given by

X k = f k (X k-1 , η k ) (1) Y k = h k (X k , Y k , ν k ) (2)
where f k and h k are possibly non-linear functions and η k and ν k are independent white noises. The underwater application motivates the choice of an implicit observation model [START_REF] Melo | Survey on advances on terrain based navigation for autonomous underwater vehicles[END_REF].

We shift from this implicit formulation (2) of the observation process to an explicit model through numerical approximations: the observation model is the result of a numerical black box code. The function h k is injective: when X k is fixed, only one Y k satisfies equation [START_REF] Del Moral | The monte-carlo method for filtering with discrete-time observations[END_REF]. The implicit problem can therefore be solved by using the fixed point method or sample approximation methods. Thus, we change the observation model [START_REF] Del Moral | The monte-carlo method for filtering with discrete-time observations[END_REF] by

Y k = h ′ k (X k , ν ′ k ) (3)
where h ′ k is still a possibly non-linear function and where ν ′ k is a new measurement noise. Noise ν ′ k is different from ν k . Indeed, the stochastic characteristics of the error between the actual measurements (2) and the numerical approximation (3) of the observation cannot be inferred. ν ′ k is divided into two parts: one part is derived from the numerical approximation and the other is due to an additional measurement noise. The measurement sensor is the output of a numerical code which is not available, thus the likelihood cannot be computed.

We want to estimate the posterior density

p k (x) p(X k = x|Y 1:k ) (4)
where

Y 1:k = [Y 1 , Y 2 , ..., Y k ]
is the vector of all the observation data accumulated up to the time k.

The state estimation consists of two steps: prediction and correction.

• The prediction step determines a predictive distribution p(X k |Y k-1 ) with respect to the dynamical model uncer-tainty p(X k |X k-1 ) and the previous posterior distribution p(X k-1 |Y k-1 ) via the Chapman-Kolmogorov equation:

p(X k |Y k-1 ) = p(X k |X k-1 )p(X k-1 |Y k-1 )dX k-1
(5) • The correction step determines the posterior distribution of the state with respect to the predictive distribution (5) and the likelihood p(Y k |X k ). From Bayes' law, one obtains:

p(X k |Y k ) = p(X k |Y k-1 )p(Y k |X k ) p(X k |Y k-1 )p(Y k |X k )dX k (6)
The particle filter estimates the posterior density by means of a Dirac mixture of N weighted particles (X i k , w i k ):

p k (x) N i=1 w i k δ x=X i k (7)
There are a multitude particle-based algorithms in the literature. We focuse on two of them, which are known for their robustness to non-linearities: the Regularized Particle Filter (RPF) ( [START_REF] Musso | Improving regularised particle filters[END_REF], [START_REF] Oudjane | Stabilité et approximations particulaires en filtrage non linéaire application au pistage[END_REF]) and the Rao-Blackwellized Particle Filter (RBPF) [START_REF] Schon | Marginalized particle filters for mixed linear/nonlinear state-space models[END_REF]. The particularity of RPF is the use of kernel density estimation sampling, in order to smooth the state density estimate. RPF is particularly useful when the process density is narrow. RBPF takes advantage of the structures of the state process by setting aside a linear part that will be solved by a Kalman filter (which is optimal if the processes are Gaussian). Using partially optimal filter allows the filter's variance to be smaller and reducing the particle state vector size limits the computational load. When the likelihood p(Y k |X k ) is unknown, the particle filter correction step cannot be explicitly achieved. To tackle this issue, ABC methods were introduced [START_REF] Del Moral | The monte-carlo method for filtering with discrete-time observations[END_REF].

B. ABC methods

ABC filters are useful when the likelihood is unknown. These methods match observations with simulated pseudoobservations and therefore avoids to assess the likelihood ([2], [START_REF] Ristic | Bayesian likelihood-free localisation of a biochemical source using multiple dispersion models[END_REF], [START_REF] Marin | Approximate bayesian computational methods[END_REF]). The objective is to evaluate the particle weights as a function of the distance between the actual observation and the simulated pseudo-observations. The posterior density p k (x) is derived by marginalizing the joint density p(X 0:k |Y 1:k ). We can get X i k samples obtained from [START_REF] Bar-Shalom | Tracking and data fusion[END_REF]. From these samples, we generate N pseudoobservations,

U i k = M(X i k ) where M is the numerical model. The model M is the h ′ function (3) without the additional measurement noise included in ν ′ k .
The closeness between the pseudo-observations U k and the true observed Y k is determined by a kernel probability density function K ǫ (Y k , U i k ). The posterior density is estimated by:

p(X 0:k |Y 1:k ) = p(X 0 ) p(X 1:k , U 1:k |Y 1:k )dU 1:k (8) p(X 0:k |Y 1:k ) ∝ p(X 0 ) k≥1 K ǫ (Y k , U k )p(U k |Y k , X k )dU k p(X k |X k-1 ) (9)
The ABC importance weights update is as follows:

w i k = w i k-1 K ǫ (Y k , U i k ) (10)
The choice of the kernel scale parameter ǫ has a strong impact on estimation convergence. It needs to be chosen in an adaptive manner, otherwise the filter may abruptly fail if the true observation at one time Y k falls in distribution tails. In the following sub-section, we propose an approach to choose this parameter.

III. ADAPTATIVE APPROXIMATE BAYESIAN COMPUTATIONAL PARTICLE FILTER (A2BC-PF)

A. A2BC

The scale parameter aims to account for every state density's modes while the measurement ambiguity is not removed. The loss of these modes leads to the degeneration of the filter weights. The parameter ǫ will be chosen adaptively in order to robustify ABC methods.

We introduce a method to determine ǫ based on the efficiency degeneracy criterion [START_REF] Kong | Sequential imputations and bayesian missing data problems[END_REF]:

N eff = 1 N i=1 w i k 2 (11)
where w i k are the particle weights and N the number of particles. This criterion is commonly used in particle filters to trigger the resampling step. Resampling is initiated whenever N eff < N th where N th = θN and 0 < θ ≤ 1 is a given threshold.

Since the likelihood is unknown, the ABC filter method is used to compute particle weights. From [START_REF] Oudjane | Stabilité et approximations particulaires en filtrage non linéaire application au pistage[END_REF], the efficiency criterion can be written as:

N eff (ǫ k ) = 1 N i=1 (w i k-1 K ǫ k (Y k , U i k )) 2 (12)
The idea is to choose ǫ ∈ I such that the efficiency criterion ( 12) is greater than the resampling threshold. At each timestep, an optimal ǫk is determined by solving the following problem:

ǫk = argmin ǫ∈I |N eff (ǫ) -N th | (13)
We want to avoid the case where the distribution of weights tends towards a Dirac distribution (i.e. when a single weight tends to unity and all the others tend to zero, see Figure 1). In order to keep the information, the case where epsilon is large is also avoided otherwise neither particle is favored. The scale parameter ǫ must be such that the efficiency criterion is equal or just above the resampling threshold. Since ǫ ∈ I, the solution ǫ to [START_REF] Schon | Marginalized particle filters for mixed linear/nonlinear state-space models[END_REF] is not guaranteed to yield N eff greater than N th . In case of N eff < N th , a conventional resample step is triggered to prevent from the filter degenerency. The kernel bandwidth has a strong impact on the degeneracy of the weights distribution, which may cause some estimation divergence. In practice, finding ǫ (13) can be solved by making numerical approximations (e.g. gradient descent or grid search).

B. A2BC-RPF and A2BC-RBPF

In this section, two improved Particle Filter are derived from the conventional RPF and RBPF, so that their correction steps satisfy [START_REF] Schon | Marginalized particle filters for mixed linear/nonlinear state-space models[END_REF].

◮ RPF ( [START_REF] Oudjane | Stabilité et approximations particulaires en filtrage non linéaire application au pistage[END_REF], [START_REF] Musso | Improving regularised particle filters[END_REF]) is based on the kernel estimation theory [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]. After the resampling step is triggered, the posterior density defined by a mixture of Dirac functions [START_REF] Kathleen | Terrain relative navigation for sensorlimited systems with application to underwater vehicles[END_REF] can be rewritten as a mixture of weighted kernels:

pk (x) = N i=1 w i k K h (X k -X i k ) (14)
When all the particles have the same weight, an optimal kernel K and an optimal bandwidth h opt can be determined by minimizing the Mean Intergrater Square Error criterion [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF], [START_REF] Musso | Improving regularised particle filters[END_REF]. This additional step is called regularization, see Algorithm 1.

◮ RBPF is used when the state vector X k can be decomposed into two sub-vectors (X n k , X l k ) such that the model is linear with respect to X l k conditionally to X n k and Y 0:k :

X n k+1 = f n k (X n k ) + F n k (X n k )X l k + G n k (x n k )η n k (15) X l k+1 = f l k (X n k ) + F l k (X n k )X l k + G l k (x n k )η l k (16) Y k = h k (X n k ) + H k (X n k )X l k + ν k ( 17
)
where

η k = [η n k , η l k ] T is a Gaussian noise with covariance ma- trix Q k = Q n k Q ln k (Q ln k ) T Q l k , v k is Gaussian with covariance
Algorithm 1 A2BC-RPF Algorithm Initialization: For i = 1, ..., N , initialize the particles X i 0 ∼ p(X 0 ) from a prior distribution and set w i -1 = 1/N . For k ≥ 0 do 1) ABC correction:

• For i = 1, ..., N , the pseudo-observations are

U i k = M(X i k ). • Determine ǫk such that ǫk = argmin ǫ k ∈I |N ef f (ǫ k ) -N th | • Update the weights wi k = w i k-1 K ǫk (Y k , U i k ) • Normalization w i k = wi k i wi k 2) Compute the estimate Xk = i w i k X i k 3) Compute the empirical covariance matrix Pk = i w i k (X i k -Xk )(X i k -Xk ) T 4) If N eff < N th do
• Apply one of the resampling procedure (see [START_REF] Li | Resampling methods for particle filtering: classification, implementation, and strategies[END_REF] for a review on resampling methods). Discard/multiply particles X i k according to high/low weights w i k and denote by X i k the selected states. Set [START_REF] Bar-Shalom | Tracking and data fusion[END_REF]. end matrix R k , the nonlinear part follows the known law p(X n 0 ) and the linear part has a Gaussian density N (X l 0 , P 0 ), see Algorithm 2. The posterior density can be written as follows:

w i k = 1/N . • Regularization step: -Compute D k such that D k D T k = Pk -Draw ζ i from a kernel -X i k = X i k + h opt D k ζ i end 5) Prediction step: For i = 1, ..., N , sample X i k+1 ∼ p(X k+1 |X i k ) see
p(X l k , X n 0:k |Y 0:k ) = p(X l k |X n 0:k , Y 0:k )p(X n 0:k |Y 0:k ) (18)
Under the hypothesis that p(X l k |X n 0:k , Y 0:k ) is Gaussian, this density is estimated by a Kalman filter. The non-linear and multimodes density p(X n 0:k |Y 0:k ) is estimated by a particle filter. The linear part follows a Gaussian law, so it is sufficient to calculate at each time step its average and its covariance by using the Kalman filter (see [START_REF] Schon | Marginalized particle filters for mixed linear/nonlinear state-space models[END_REF]). The particle prediction will be given by:

p(X n 0:k |Y 0:k-1 ) = p(X n k |X n 0:k-1 , Y 0:k-1 )p(X n 0:k-1 |Y 0:k-1 ) (19)
This density is Gaussian with mean and covariance matrix given by equation (25) in [START_REF] Schon | Marginalized particle filters for mixed linear/nonlinear state-space models[END_REF]. The resulting approximation is given by:

pk (x) = N i=1 w i k N (x l ; X l,i k , P i k )δ x n =X n,i k ( 20
)
where P i k is the covariance matrix of the linear part X l,i k .

Algorithm 2 A2BC-RBPF Algorithm

Initialization: For i = 1, ..., N , sample X n,i 0 ∼ p(X n 0 ) and set {X l,i 0|-1 , P i 0|-1 } i=1:N = {X l 0 , P l 0 } and w i -1 = 1/N . For k ≥ 0 do 1) ABC correction:

• For i = 1, ..., N , the pseudo-observations are

U i k = M(X n,i k , X l,i k|k-1 ) • Determine ǫk such that ǫk = argmin ǫ k ∈I |N ef f (ǫ k ) -N th | • Update the weights wi k = w i k-1 K ǫk (Y k , U i k ) • Normalization w i k = wi k i wi k 2) Compute the estimate Xn k = i w i k X n,i k 3) Compute the empirical covariance matrix P n k = i w i k (X n,i k -Xn k )(X n,i k -Xn k ) T 4) If N eff (ǫ k ) < N th do • Discard/multiply particles X i k = [X n,i k , X l,i k|k-1
] according to high/low weights w i k and denote by X i k the selected states. Set

w i k = 1/N . • Regularization step: -Compute D k such that D k D T k = Pk where Pk = diag([ P n k , P k|k-1 ]). -Draw ζ i from a kernel -X i k = X i k + h opt D k ζ i end 5) Kalman correction: compute X l,i
k and P i k according to equation ( 22) in [START_REF] Schon | Marginalized particle filters for mixed linear/nonlinear state-space models[END_REF]. 6) Particle prediction: sample X n,i k+1 using (19) (see equation (25) in [START_REF] Schon | Marginalized particle filters for mixed linear/nonlinear state-space models[END_REF]). 7) Kalman prediction: compute X l,i k+1|k and P i k+1|k according to equation (23) in [START_REF] Schon | Marginalized particle filters for mixed linear/nonlinear state-space models[END_REF]. end A regularization step is added in order to increase the robustness of the A2BC-RBPF filter. For a same number of particles, RBPF is more computationally intensive since, it requires one Kalman filter per particle. However, it is possible to obtain a compromise between calculation cost and performance in favour of the RBPF for many applications. For example, in TAN, the Kalman filter covariance matrix is independent of the non-linear part, which reduces the complexity of the RBPF since there is only one covariance matrix to update.

IV. SIMULATION EXAMPLE

To illustrate the behavior of the resulting A2BC filters, an application to underwater Terrain Aided Navigation (TAN) is presented. The observation equation is described and numerical results are presented.

A. Underwater terrain aided navigation

TAN can provide a drift-free navigation tool for underwater vehicles (UV), yielding a powerful alternative to current navigation methods which include resurfacing for GPS. TAN generates vehicle position estimates by correlating terrain measurements obtained by a multi-beam sensor (see Figure 2) together with stored terrain maps. Fig. 2: Multi-beam sensor.

Define X = [P, V ] T the state vector composed of P = [p x , p y , p z ] T the position vector expressed in meters and V = [v x , v y , v z ] T the vector of velocities expressed in meter per second. A simple case of state dynamics is chosen, known as NCV model [START_REF] Bar-Shalom | Tracking and data fusion[END_REF]. The dynamical model is written in the following discrete way:

P k = P k-1 + ∆tV k-1 + η X k V k = V k-1 + η V k (21)
where ∆t is the discretization time-step and η

k = [η X k , η V k ] T is the state noise.
The measurement is made of m beams: Y = [r 1 , ..., r m ] T . Each measurement r i returns the distance between the UV and the seabed (see Meduna [START_REF] Kathleen | Terrain relative navigation for sensorlimited systems with application to underwater vehicles[END_REF] for details):

r i = (p x -p i x ) 2 + (p y -p i y ) 2 + (p z -h(p i x , p i y )) 2 + ν i (22)
where h is the terrain depth. The terrain function is the Digital Elevation Model (DEM). In our case, DEM is defined by a regular grid of elevation values, see Figure 3. The spatial resolution of the chosen map is about 100 meters. The intersection point of beam r i direction vector with the terrain is denoted (p i x , p i y , p i z = h(p i x , p i y )). Since the coordinates of the intersection point depend on the beam vector r i , the measurement model is implicit. The measurement equation ( 22) has to be rewritten as [START_REF] Neil J Gordon | Novel approach to nonlinear/non-gaussian bayesian state estimation[END_REF]. There are several ways to get an explicit formula, e.g. ray tracing, grid search. In this paper, the explicit formula is the output of a numerical model named M, that introduces some unknown sampling noise.

The law of this noise cannot be easily inferred which motivate the use of the method developed in Section III.

B. Performance criteria

In order to compare the algorithms, we use the following criteria, evaluated for N mc ∈ N Monte-Carlo simulations:

• The Root Mean Square Error (RMSE):

RM SE k = Nmc i=1 ( X i k -X k ) T ( X i k -X k ) N mc (23)
where X k is the kth estimated of the state. We will compare these RMSEs with the Posterior Cramér Rao Bound (PCRB) which is calculated according to the Tichavský recursive formula [START_REF] Tichavsky | Posterior cramér-rao bounds for discrete-time nonlinear filtering[END_REF]. The PCRB is approximated over 100 state samples at each time-step.

• The number of non-convergence:

The filter is said to not converge if, at the end of the trajectory, during the last 5 consecutive measurement times, the state estimate X k leaves the confidence ellipsoid Γ k given by the PCRB, such that

Γ k = {X k |(X k -X k ) T P CRB -1 k (X k -X k ) ≤ α 2 } (24)
where the threshold α is equal to the probability P(X 2 (d) ≤ α 2 ) = 0.99 with d the dimension of the state vector.

C. Simulation and results

The kernel in [START_REF] Oudjane | Stabilité et approximations particulaires en filtrage non linéaire application au pistage[END_REF] is chosen as a Cauchy function whose scale parameter ǫ is fixed for conventional filters and determine adaptively with A2BC filters. The linear part of the state in RBPF is the velocity.

For 100 Monte Carlo trials, the number of non-convergence is shown in Table I. The initial number of non-convergence for RBPF was too high (∼ 70%). Thereafter, RBPF will designate an improved version where a regularization step was added. The number of non-convergence is smaller for A2BC filters that are more robust than the benchmark filters. The reference trajectory is located in an ambiguous area of the map. Figure 4 illustrates the trajectories estimated by the two proposed algorithms in the same conditions (the same initial errors and measurements realizations). The filters converge quickly towards the reference trajectory. The criterion N eff of the A2BC filters is plotted in Figure 5 with the resampling threshold N th . The criterion varies over time to avoid the resampling step and therefore maintains the modes of the posterior density as long as possible (i.e. N eff is forced to remain near N th ). It's the reason why in Figure 6, the A2BC-RPF has a higher error than the other filters between time 7 minutes and 17 minutes. Figures 6 and7 show the RMSE results of the four filters tested with similar simulation conditions. Only convergence cases are used to plot the curves on these two figures. Convergence generally occurs around 20 minutes. The curves follow the tendency of the approximate PCRB. The final accuracy of the filters is less than 100 meters in position, which corresponds to the spatial resolution of the map. The velocity accuracy is also satisfactory.

Filters

V. CONCLUSIONS

In some cases, the observation model is a numerical black box, which leads us to consider an unknown likelihood. We introduce a measurement updating procedure based on an adaptive kernel derived from Approximate Bayesian Computational filters. The proposed method is implemented in two wellknown particles filter: Regularized Particle Filter and Rao-Blackwellized Particle Filter. Simulation results demonstrate that the proposed method significantly increases the robustness and the accuracy of particle-like filters while remaining computationally efficient.

  Correction step using a given kernel bandwidth epsilon and the subsequent weight distribution. State space (&) (b) If the kernel bandwidth is too narrow, the degeneracy phenomenon arrises.

Fig. 1 :

 1 Fig. 1: Relation between the scale parameter of the kernel and the weights distribution.

Fig. 3 :

 3 Fig. 3: Map of the contours lines of the California coast (31 ˚88' N, 121 ˚27' W). The colorbar represents the depth levels in meters.

•

  Sampling period: ∆T = 10 seconds • Number of bathymetric measurements: 200 • Time at the end of the trajectory: T = 33 minutes • Number of beams: 3 • Number of particles: N = 5000 • Resampling threshold: N th = 0.75N • ǫ domain (13): I = [0.1, 20] meters • Standard deviation of each beam range: σ R = 10 meters The initial uncertainty for the position P is σ P = diag([1000, 1000, 100]) meters and the initial uncertainty for the velocity V is σ V = diag([0.5, 0.5, 0.5]) meters per second. The initial state is X 0 = [120000, 180000, -100, 5, 5, 0.05] T . The horizontal velocity vector is 7 meters per second. For the regularization step, we chose an Epanechnikov kernel. The filter results are compared for the same number of particles N .

Fig. 4 :

 4 Fig. 4: Comparison between the true trajectory and the tajectories estimated by A2BC-RPF and A2BC-RBPF on the map of the California coast. The colorbar represents the depth levels in meters.

Fig. 5 :

 5 Fig. 5: Values taken by N eff over time and resampling threshold N th .

  (a) RMSE of the position on the axis X.

  RMSE of the position on the axis Y.

Fig. 6 :

 6 Fig. 6: PCRB and RMSE of the horizontal position in meters.

  (a) RMSE of the velocity on the axis X.

  RMSE of the velocity on the axis Y.

Fig. 7 :

 7 Fig. 7: PCRB and RMSE of the horizontal velocity in meters per second. The error at the end of the trajectory of the filters in horizontal position RMSE HP T and horizontal velocity RMSE HV T

Filters

TABLE I :

 I Table of percentage rate of non-convergence for 100 Monte Carlo trials.

		% of non-convergence
	RPF	8
	A2BC-RPF	3
	RBPF	10
	A2BC-RBPF	6

TABLE II :

 II Table of RMSE at the end of the trajectory (in meters for position and in meters per second for velocity).

		RMSE HP T	RMSE pz T	RMSE HV T
	RPF	106.53	1.47	0.195
	A2BC-RPF	104.67	0.83	0.188
	RBPF	105.23	1.46	0.197
	A2BC-RBPF	99.90	0.78	0.177