text:
Xiaojing Fan

Deqiang Han
email: deqhan@xjtu.edu.cn

Jean Dezert
email: jean.dezert@onera.fr

Yi Yang

User-specified Optimization Based Transformation of Fuzzy Membership into Basic Belief Assignment

Keywords: evidence theory, basic belief assignment, fuzzy membership function, optimization, transformation

To combine different types of uncertain information from different sources under different frameworks, we need transformations between different frameworks. For the transformation of a fuzzy membership function (FMF) into a basic belief assignment (BBA), several approaches have been proposed. Among these approaches, the uncertainty optimization based transformations can provide BBAs without predefining focal elements. However, these two transformations, which respectively use the uncertainty maximization and minimization criteria, emphasize the extreme cases of uncertainty. We expect to obtain a BBA, which is the trade-off between the two BBAs obtained by solving the uncertainty maximization and minimization, to avoid extreme attitudinal bias. In this paper, we propose two transformations of an FMF into a BBA by using a user-specified weighting factor to obtain such a trade-off (or balanced) BBA. Some examples and related analyses are provided to show the rationality and effectiveness of the proposed transformations.

I. INTRODUCTION

In the information fusion, we need to deal with a large amount of uncertain information. Various types of uncertainty theories have been proposed to deal with different types of uncertainty, e.g., the probability theory, fuzzy set theory [START_REF] Zadeh | Fuzzy Sets[END_REF], possibility theory [START_REF] Zadeh | Fuzzy sets as a basis for a theory of possibility[END_REF], rough set theory [START_REF] Pawlak | Rough sets[END_REF] and Dempster-Shafer evidence theory (DST) [START_REF] Shafer | A Mathematical Theory of Evidence[END_REF] etc. When we fuse the information from different sources under different theoretical frameworks, we need the transformation between different frameworks.

For the information represented by the FMF and BBA, we can transform an FMF into a BBA. Then, we can combine the BBAs to implement the information fusion. There have been proposed many transformations of an FMF into a BBA [START_REF] Bi | Combining Multiple Classifiers Using Dempsters Rule of Combination for Text Categorization[END_REF]- [START_REF] Ma | An Improved Method to Transform Triangular Fuzzy Number Into Basic Belief Assignment in Evidence Theory[END_REF]. In [START_REF] Bi | Combining Multiple Classifiers Using Dempsters Rule of Combination for Text Categorization[END_REF], Bi et al. proposed a transformation that normalizes a given FMF to generate a BBA with singleton focal elements only. By using the α-cut approach, Florea et al. [START_REF] Florea | Approximation techniques for the transformation of fuzzy sets into random sets[END_REF] transformed an FMF into a BBA with focal elements nested in order. However, these two approaches above have to predefine the focal elements, which lack of intuitiveness and objectiveness. Han et al. [START_REF] Han | Novel Approaches for the Transformation of Fuzzy Membership Function into Basic Probability Assignment based on Uncertainty Optimization[END_REF] proposed two approaches without predefining focal elements. These two approaches can provide BBAs by solving constrained uncertainty maximization and minimization.

For the two transformations of Han et al. [START_REF] Han | Novel Approaches for the Transformation of Fuzzy Membership Function into Basic Probability Assignment based on Uncertainty Optimization[END_REF], both two objective functions are the ambiguity measure (AM ) and their constraints are mainly constructed based on the given FMF. Their rationality and effectiveness are both justified in [START_REF] Han | Novel Approaches for the Transformation of Fuzzy Membership Function into Basic Probability Assignment based on Uncertainty Optimization[END_REF]. During the process of solving optimization problems, these two transformations emphasize on the minimum and maximal uncertainties of the BBA, respectively. We think that the BBA being the trade-off (or balanced) between the two BBAs obtained by solving the uncertainty maximization and minimization is more preferred, which might avoid being "onesided" on the uncertainty degree. In this paper, we propose two approaches by using a user-specified weighting factor to determine BBAs. One transformation is the weighted average by using the user-specified weighting factor with the two BBAs obtained by optimization based transformations [START_REF] Han | Novel Approaches for the Transformation of Fuzzy Membership Function into Basic Probability Assignment based on Uncertainty Optimization[END_REF]. The other transformation brings out a trade-off BBA by solving a constrained minimization problem. The objective function is based on the user-specified weighting factor, the distance of evidence and the two BBAs obtained by uncertainty optimization. The constraints are mainly based on the given FMF. That is, each of our proposed transformations can transform an FMF into a BBA, which can be considered as the trade-off between the two BBAs obtained with uncertainty optimization. Some examples and related analyses are provided to justify the proposed transformations.

II. PRELIMINARY

A. Basics of the Theory of Belief Functions

The theory of belief functions [START_REF] Shafer | A Mathematical Theory of Evidence[END_REF], introduced historically by Shafer in DST, is a powerful framework for uncertainty modeling and reasoning. Let Θ = {θ 1 , θ 2 , ..., θ n } be the frame of discernment (FOD). Under the closed world assumption, the FOD is mutually exclusive and exhaustive. The BBA (also called a mass function) is defined on the power set of Θ, which can be denoted by a function m :

2 Θ → [0, 1] satisfying A⊆Θ m(A) = 1, m(∅) = 0 (1)
where ∅ denotes the empty set. ∀A, if m(A) > 0, then A is called a focal element. m(A) denotes the evidence support to the proposition A.

The belief function Bel for all A ⊆ Θ, as:

Bel(A) = B⊆A m(B) (2)
The plausibility function P l for all A ⊆ Θ, as:

P l(A) = A∩B =∅ m(B) (3)
Suppose there are two independent BBAs m 1 and m 2 on the same FOD. Historically Shafer proposed Dempster's rule to combine two (or more) BBAs. Dempster's rule of combination is

m(A) =      0, A = ∅ B∩C=A m 1 (B)m 2 (C) 1 -K , A = ∅ (4)
where K = B∩C=∅ m 1 (B)m 2 (C) represents the conflict coefficient between two BBAs. There exist other alternative combination rules [START_REF] Ilin | Information Fusion with Belief Functions: a Comparison of Proportional Conflict Redistribution PCR5 and PCR6 Rules for Networked Sensors[END_REF], [START_REF] Smarandache | On the consistency of PCR6 with the averaging rule and its application to probability estimation[END_REF].

B. Uncertainty Measure of a BBA

The uncertainty of a BBA includes two types: the discord and the non-specificity. Different measures of uncertainty [START_REF] Höhle | Entropy with respect to plausibility measures[END_REF]- [START_REF] Harmanec | Measuring total uncertainty in Dempster-Shafer theory: a novel approach[END_REF] have been proposed, e.g., the non-specificity measure [START_REF] Dubois | A note on measures of specificity for fuzzy sets[END_REF], the ambiguity measure (AM ) [START_REF] Jousselme | Measuring ambiguity in the evidence theory[END_REF] and the aggregated uncertainty (AU ) [START_REF] Harmanec | Measuring total uncertainty in Dempster-Shafer theory: a novel approach[END_REF]. The definition of AM is as follows:

AM (m) = - θ∈Θ BetP m (θ) log 2 (BetP m (θ)) (5)
where BetP m (θ) = θ∈A⊆Θ m(A)/ |A| is the pignistic probability [START_REF] Smets | The transferable belief model[END_REF]. |A| denotes the cardinality of the set A.

III. TRANSFORMATION OF FMF INTO BBA

A. Concept of Fuzzy Set

Fuzzy sets [START_REF] Zadeh | Fuzzy Sets[END_REF] were proposed by Zadeh to describe the concepts without precise definitions. Let Θ be the universe of discourse (equivalent to FOD in the belief functions). A fuzzy membership function is denoted by u = µ(θ), θ ∈ Θ.

For µ : Θ → [0, 1], µ(θ) ∈ [0, 1] is called the degree of membership for θ.
B. Traditional Transformations of FMF into BBA a) Transformations with the predefinition of focal elements: For a given FMF, two available types of transformations below can provide a BBA, which have to predefine the focal elements. Suppose that the FOD is Θ = {θ 1 , θ 2 , ..., θ n } and the given FMF is µ = [µ(θ 1 ), µ(θ 2 ), ..., µ(θ n )]. The obtained BBA is represented by m.

In the work of Bi et al. [START_REF] Bi | Combining Multiple Classifiers Using Dempsters Rule of Combination for Text Categorization[END_REF], the BBA is determined as follows:

m({θ i }) = µ(θ i ) n j=1 µ(θ j ) (6)
This approach predefines all focal elements as singletons, and it is the result of normalization for the given FMF.

Another transformation with the predefinition of focal elements is the work of Florea et al. [START_REF] Florea | Approximation techniques for the transformation of fuzzy sets into random sets[END_REF] by using the α-cut approach. Suppose that µ(θ 1 ), µ(θ 2 ), ..., µ(θ n ) are sorted into ascending order as 0 = α 0 < α 1 < α 2 < ... < α M ≤ 1, where M ≤ |Θ|. The BBA is determined by using the transformation [START_REF] Florea | Approximation techniques for the transformation of fuzzy sets into random sets[END_REF] as follows:

m(A j ) = α j -α j-1 α M (7)
where A j = {θ i ∈ Θ|µ(θ i ) ≥ α j }, i = 1, 2, ..., n, j = 1, 2, ..., M . This transformation predefines the focal elements nested in order for the given FMF. Both two approaches can transform an FMF into a BBA. However, the transformations with the predefinition of focal elements lack of intuitiveness and objectiveness. For a given FMF, the optimization based transformations can obtain a BBA without predefining the focal elements.

b) Transformations based on the uncertainty optimization: In the work of Han et al. [START_REF] Han | Novel Approaches for the Transformation of Fuzzy Membership Function into Basic Probability Assignment based on Uncertainty Optimization[END_REF], the two transformations that have no predefinition of focal elements are obtained by solving the uncertainty maximization and minimization. Suppose that the FOD is Θ = {θ 1 , θ 2 , ..., θ n } and the given FMF is µ = [µ(θ 1 ), µ(θ 2 ), ..., µ(θ n )]. The obtained BBA is represented by m.

There exists a relationship [START_REF] Kampe De Feriet | Interpretation of fuzzy membership functions of fuzzy sets in terms of plausibility and belief[END_REF] between the FMF and the belief function or plausibility function. When n i=1 µ(θ i ) ≥ 1, the FMF is equivalent to a singleton plausibility function, which is denoted by

P l({θ i }) = {θi}∩A =∅ m(A) = µ(θ i ), ∀ {θ i } ⊆ Θ (8)
It is the necessary and sufficient condition for the FMF to be a singleton plausibility function. When n i=1 µ(θ i ) ≤ 1, the FMF is equivalent to a singleton belief function, which is denoted by

Bel({θ i }) = A⊆{θi} m(A) = µ(θ i ), ∀ {θ i } ⊆ Θ (9)
Similarly, it is the necessary and sufficient condition for the FMF to be a singleton belief function.

The detailed proof of the above relationships are given in [START_REF] Kampe De Feriet | Interpretation of fuzzy membership functions of fuzzy sets in terms of plausibility and belief[END_REF].

There is a BBA transformed from a given FMF, and the FMF and BBA satisfy Eq. ( 8) or Eq. ( 9). Then, n linear equations for the corresponding relations can be obtained. In addition, one has A⊆Θ m(A) = 1. There exist n + 1 linear equations. However, except for m(∅) = 0, in the worst case there are 2 n -1 focal elements to assign the belief. The n + 1 linear equations with respect to 2 n -1 undetermined variables, which is an under-determined problem, i.e., it usually has multiple solutions.

Therefore, to obtain a unique BBA, Han et al. [START_REF] Han | Novel Approaches for the Transformation of Fuzzy Membership Function into Basic Probability Assignment based on Uncertainty Optimization[END_REF] established two uncertainty optimization based transformations for the given FMF. Both two objective functions are AM and the constraints are mainly based on the given FMF.

The objective function of the uncertainty maximization problem and the corresponding constraints are as follows: When

n i=1 µ(θ i ) ≥ 1, max m - n i=1 [BetP m (θ i ) log 2 (BetP m (θ i ))] s.t.        {θi}∩A =∅ m(A) = µ(θ i ), ∀ {θ i } ⊆ Θ A⊆Θ m(A) = 1 0 ≤ m(A) ≤ 1 (10) When n i=1 µ(θ i ) ≤ 1, max m - n i=1 [BetP m (θ i ) log 2 (BetP m (θ i ))] s.t.        A⊆{θi} m(A) = µ(θ i ), ∀ {θ i } ⊆ Θ A⊆Θ m(A) = 1 0 ≤ m(A) ≤ 1 (11)
In the sequel, this transformation is represented by "T max " for convenience.

The objective function of the uncertainty minimization problem and the corresponding constraints are as follows: When

n i=1 µ(θ i ) ≥ 1, min m - n i=1 [BetP m (θ i ) log 2 (BetP m (θ i ))] s.t.        {θi}∩A =∅ m(A) = µ(θ i ), ∀ {θ i } ⊆ Θ A⊆Θ m(A) = 1 0 ≤ m(A) ≤ 1 (12) When n i=1 µ(θ i ) ≤ 1, min m - n i=1 [BetP m (θ i ) log 2 (BetP m (θ i ))] s.t.        A⊆{θi} m(A) = µ(θ i ), ∀ {θ i } ⊆ Θ A⊆Θ m(A) = 1 0 ≤ m(A) ≤ 1 (13)
In the sequel, this transformation is represented by "T min " for convenience.

The unique BBA can be determined without predefining focal elements by using "T max " or "T min ". The obtained BBA is the optimal solution of the uncertainty maximization or minimization. During the process of transforming an FMF into a BBA, "T max " and "T min " emphasize on the maximal and minimum uncertainty cases of the obtained BBA, respectively. We think that the BBA being the trade-off between the two BBAs obtained by using "T max " and "T min " is more preferred, which might avoid bias in terms of uncertainty degree.

IV. TRANSFORMATIONS WITH USER-SPECIFIED WEIGHTING FACTOR

As aforementioned, we can obtain two BBAs by using "T max " and "T min ", respectively. Based on these two BBAs, we aim to construct a transformation to determine a tradeoff BBA. The trade-off BBA which satisfies the relationship between the FMF and the singleton plausibility or singleton belief. We use a user-specified weighting factor to influence how close the trade-off BBA is to each of the two BBAs above. Suppose that the user-specified weighting factor is represented by α and 0 ≤ α ≤ 1. When α → 0, the trade-off BBA is close to the BBA obtained by using "T min ". When α → 1, the trade-off BBA is close to the BBA obtained by using "T max ". To meet the requirements above, we propose two different approaches to determine the trade-off BBAs.

A. Weighted Average based Transformation

Let Θ = {θ 1 , θ 2 , ..., θ n } be the FOD. The given FMF is represented by µ = [µ(θ 1 ), µ(θ 2 ), ..., µ(θ n )]. Suppose that the BBA obtained by using "T min " is denoted by m min , and the BBA obtained by using "T max " is denoted by m max . The user-specified weighting factor is denoted by α (0 ≤ α ≤ 1). The trade-off BBA is denoted by m. The Weighted Average of m min and m max can bring out a trade-off BBA as follows:

m(A) = (1 -α) • m min (A) + α • m max (A) (14)
where A ⊆ Θ. In the sequel, the transformation based on the weighted average (WA) is denoted by "T WA " for convenience. The BBA obtained in ( 14) is an admissible BBA and it satisfies the constraints established based on FMF.

According to Eq. ( 14), the following conditions can be satisfied:

A⊆Θ m(A) = 1, 0 ≤ m(A) ≤ 1 (15)
For the transformation of an FMF into a BBA, it is necessary that the obtained BBA satisfies the relationship between the FMF and the singleton plausibility or singleton belief. Although "T WA " is a simple and direct transformation of an FMF into a trade-off BBA, it also satisfies the relationship. The proof is provided below. When n i=1 µ(θ i ) ≥ 1, m min and m max satisfy Eq. ( 6), respectively, i.e., P l min ({θ i }) = P l max ({θ i }) = µ(θ i ), i = 1, 2, ..., n. According to Eq. ( 14),

P l({θ i }) = {θi}∩A =∅ m(A) = {θi}∩A =∅ [(1-α)•m min (A)+α•m max (A)] = (1-α)• {θi}∩A =∅ m min (A)+α• {θi}∩A =∅ m max (A) = (1 -α)•P l min ({θ i }) + α•P l max ({θ i }) = (1 -α) • µ(θ i ) + α • µ(θ i ) = µ(θ i ) (16)
Similarly, when n i=1 µ(θ i ) ≤ 1, m min and m max satisfy Eq. ( 7), respectively, i.e., Bel min ({θ i }) = Bel max ({θ i }) = µ(θ i ), i = 1, 2, ..., n. According to Eq. ( 14),

Bel ({θ i }) = m(A) = (1 -α) • m min (A) + α • m max (A) = (1 -α)•Bel min ({θ i })+α•Bel max ({θ i }) = (1 -α) • µ(θ i ) + α • µ(θ i ) = µ(θ i ) (17)
For the trade-off BBA, when n i=1 µ(θ i ) ≥ 1, the given FMF is equivalent to the corresponding singleton plausibility. When n i=1 µ(θ i ) ≤ 1, the given FMF is equivalent to the corresponding singleton belief. That is, "T WA " can transform the given FMF into the trade-off BBA, which satisfies the relationship between the FMF and the BBA.

B. User-specified Optimization based Transformation

Let the FOD be Θ = {θ 1 , θ 2 , ..., θ n }. The given FMF is denoted by µ = [µ(θ 1 ), µ(θ 2 ), ..., µ(θ n )]. Suppose that m min and m max denote the BBAs obtained by "T min " and "T max ", respectively. The user-specified weighting factor is denoted by α (0 ≤ α ≤ 1). The trade-off BBA is represented by m.

The user-specified weighting factor is used to influence the similarity between the trade-off BBA and m min (or m max ). The degree of similarity between two BBAs is represented by the distance of evidence. We can use the Jousselme's distance [START_REF] Jousselme | A new distance between two bodies of evidence[END_REF], which is a strict metric defined as

d J (m a , m b ) = 1 2 (m a -m b )D(m a -m b ) (18)
where

D(A, B) = |A ∩ B| / |A ∪ B|, A ⊆ Θ, B ⊆ Θ.
According to Eq. ( 18), the obtained BBA is more similar to m min , if d J (m, m min ) is smaller. If d J (m, m max ) is smaller, the obtained BBA is more similar to m max .

To obtain the trade-off BBA between m min and m max , a relationship between the user-specified weighting factor and the distance of evidence can be constructed. When α is given from 0 to 1, with the decreasing of d J (m, m min ), the value of d J (m, m max ) is increasing. Then we can establish the following equation:

d J (m, m min ) d J (m, m max ) = α 1 -α (19)
The BBA satisfies Eq. ( 19) may not always exist. If the following function (equivalent to Eq. ( 19)) achieves the minimum value, then the trade-off BBA is obtained.

obj(m) = [(1 -α)•d J (m, m min ) -α•d J (m, m max )] 2 (20)
When α is given, we can establish a constrained minimization problem to transform an FMF into a BBA. The objective function is Eq. ( 20) and the constraints are mainly based on Eq. ( 6) or Eq. [START_REF] Han | Novel Approaches for the Transformation of Fuzzy Membership Function into Basic Probability Assignment based on Uncertainty Optimization[END_REF]. The transformation of an FMF into a tradeoff BBA is obtained by solving the user-specified optimization problem as follows: When

n i=1 µ(θ i ) ≥ 1, min m [(1 -α)•d J (m, m min ) -α•d J (m, m max )] 2 s.t.        {θi}∩A =∅ m(A) = µ(θ i ), ∀ {θ i } ⊆ Θ A⊆Θ m(A) = 1 0 ≤ m(A) ≤ 1 (21) When n i=1 µ(θ i ) ≤ 1, min m [(1 -α)•d J (m, m min ) -α•d J (m, m max )] 2 s.t.        A⊆{θi} m(A) = µ(θ i ), ∀ {θ i } ⊆ Θ A⊆Θ m(A) = 1 0 ≤ m(A) ≤ 1 (22)
In the sequel, the transformation based on the user-specified optimization (USO) is denoted by "T USO " for convenience.

For a given FMF, the trade-off BBA can be obtained by using "T USO ", which is a user-specified optimization based transformation.

V. EXPERIMENTS

In this section, we provide some examples to illustrate how to transform an FMF into a trade-off BBA using our approaches. Here, we use the optimization toolbox in the Matlab TM to solve the optimization problems under constraints.

A. Example 1

Let the FOD be Θ = {θ 1 , θ 2 , θ 3 }. The given FMF is µ(θ 1 ) = 0.9, µ(θ 2 ) = 0.7, µ(θ 3 ) = 0.3. Suppose that m min and m max are the BBAs obtained by using "T min " and "T max ", respectively. We just list the corresponding BBAs for α = 0, 0.3, 0.7 and 1.

This FMF satisfies 3 i=1 µ(θ i ) = 1.9 > 1. Therefore, the given FMF is equivalent to the singleton plausibility. The BBAs obtained by using "T WA " and "T USO " are listed in the Table I and Table II, respectively. By using "T WA " and "T USO ", when α = 0, the obtained BBAs are identical to m min , and the values of AM are the minimum uncertainty. When α → 0, the obtained BBA is similar to m min , and its uncertainty is close to the minimum uncertainty.

Similarly, when α = 1, the obtained BBAs are identical to m max , and the values of AM are the maximal uncertainty. When α → 1, the obtained BBA is similar to m max , and its uncertainty is close to the maximal uncertainty.

B. Example 2

Let the FOD be Θ = {θ 1 , θ 2 , θ 3 , θ 4 }. The given FMF is µ(θ 1 ) = 1, µ(θ 2 ) = 0.2, µ(θ 3 ) = 0.3, µ(θ 4 ) = 0.3. Suppose that m min and m max are the BBAs obtained by using "T min " and "T max ", respectively. We just list the corresponding BBAs for α = 0, 0.2, 0.8 and 1.

According to

4 i=1 µ(θ i ) = 1.8 > 1,
the FMF is equivalent to the singleton plausibility. In the Table III

α BBA AM α = 0 m({θ 1 }) = 0.3, m({θ 1 , θ 2 }) = 0.3 1.3367 m({θ 2 }) = 0.1, m(Θ) = 0.3 α = 0.3 m({θ 1 }) = 0.21, m({θ 1 , θ 2 }) = 0.42 1.4003 m({θ 2 }) = 0.07, m({θ 1 , θ 3 }) = 0.06 m({θ 3 }) = 0.03, m(Θ) = 0.21 α = 0.7 m({θ 1 }) = 0.09, m({θ 1 , θ 2 }) = 0.58 1.4730 m({θ 2 }) = 0.03, m({θ 1 , θ 3 }) = 0.14 m({θ 3 }) = 0.07, m(Θ) = 0.09 α = 1 m({θ 1 , θ 2 }) = 0.7, m({θ 3 }) = 0.1 1.5129 m({θ 1 , θ 3 }) = 0.2 TABLE II USING "T USO " TO OBTAIN BBAS IN EXAMPLE 1. α BBA AM α = 0 m({θ 1 }) = 0.3, m({θ 1 , θ 2 }) = 0.3 1.3367 m({θ 2 }) = 0.1, m(Θ) = 0.3 α = 0.3 m({θ 1 }) = 0.2027, m({θ 1 , θ 2 }) = 0.4145 1.3969 m({θ 2 }) = 0.0828, m({θ 1 , θ 3 }) = 0.0802 m({θ 3 }) = 0.0172, m(Θ) = 0.2026 α = 0.7 m({θ 1 }) = 0.1037, m({θ 1 , θ 2 }) = 0.5962 1.4828 m({θ 2 }) = 0.0001, m({θ 1 , θ 3 }) = 0.0965 m({θ 3 }) = 0.0999, m(Θ) = 0.1036 α = 1 m({θ 1 , θ 2 }) = 0.7, m({θ 3 }) = 0.1 1.5129 m({θ 1 , θ 3 }) = 0.2
the BBAs obtained by using "T WA " and "T USO " are listed, respectively.

When α = 0, the obtained BBAs are identical to m min , and the values of AM are the minimum uncertainty. When α = 1, the obtained BBAs are identical to m max , and the values of AM are the maximal uncertainty.

In the Table III and Table IV, when α → 0, the obtained BBA is similar to m min , and its uncertainty is close to the minimum uncertainty. When α → 1, the obtained BBA is similar to m max , and its uncertainty is close to the maximal uncertainty.

α BBA AM α = 0 m({θ 1 }) = 0.7, m(Θ) = 0.2 1.0896 m({θ 1 , θ 3 , θ 4 }) = 0.1 α = 0.2 m({θ 1 }) = 0.60, m({θ 1 , θ 2 }) = 0.04 1.2099 m({θ 1 , θ 3 }) = 0.06, m({θ 1 , θ 4 }) = 0.06 m({θ 1 , θ 3 , θ 4 }) = 0.08, m(Θ) = 0.16 α = 0.8 m({θ 1 }) = 0.3, m({θ 1 , θ 2 }) = 0.16 1.5122 m({θ 1 , θ 3 }) = 0.24, m({θ 1 , θ 4 }) = 0.24 m({θ 1 , θ 3 , θ 4 }) = 0.02, m(Θ) = 0.04 α = 1 m({θ 1 }) = 0.2, m({θ 1 , θ 2 }) = 0.2 1.5955 m({θ 1 , θ 3 }) = 0.3, m({θ 1 , θ 4 }) = 0.3 TABLE IV USING "T USO " TO OBTAIN BBAS IN EXAMPLE 2. α BBA AM α = 0 m({θ 1 }) = 0.7, m(Θ) = 0.2 1.0896 m({θ 1 , θ 3 , θ 4 }) = 0.1 α = 0.2 m({θ 1 }) = 0.6008, m({θ 1 , θ 2 }) = 0.0471 1.2133 m({θ 1 , θ 3 }) = 0.0484 m({θ 1 , θ 2 , θ 3 }) = 0.0037 m({θ 1 , θ 4 }) = 0.0515 m({θ 1 , θ 2 , θ 4 }) = 0.0006 m({θ 1 , θ 3 , θ 4 }) = 0.0993, m(Θ) = 0.1486 α = 0.8 m({θ 1 }) = 0.2990, m({θ 1 , θ 2 }) = 0.1665 1.5227 m({θ 1 , θ 3 }) = 0.2345, m({θ 1 , θ 4 }) = 0.2010 m({θ 1 , θ 2 , θ 4 }) = 0.0335 m({θ 1 , θ 3 , θ 4 }) = 0.0655 α = 1 m({θ 1 }) = 0.2, m({θ 1 , θ 2 }) = 0.2 1.5955 m({θ 1 , θ 3 }) = 0.3, m({θ 1 , θ 4 }) = 0.3 C. Example 3 Let the FOD be Θ = {θ 1 , θ 2 , θ 3 , θ 4 }. The given FMF is µ(θ 1 ) = 0.6, µ(θ 2 ) = 0.1, µ(θ 3 ) = 0.2, µ(θ 4 ) = 0.1.
We just list the corresponding BBAs for α = 0, 0.4, 0.9 and 1.

This FMF satisfies

4 i=1 µ(θ i ) = 1.
The FMF is equivalent to the singleton plausibility or singleton belief. In the Table V and Table VI, the BBAs obtained by using "T WA " and "T USO " are listed, respectively.

In the Table V and Table VI, the BBAs obtained when α = 0 are identical to the BBAs obtained when α = 1, i.e., the two BBAs obtained by using "T min " and "T max " are the same. Therefore, ∀α ∈ [0, 1], the obtained BBAs are without the influence of α. When α is given from 0 to 1, all the obtained BBAs are Bayesian belief functions and are identical.

α BBA AM α = 0 m({θ 1 }) = 0.6, m({θ 2 }) = 0.1 1.5710 m({θ 3 }) = 0.2, m({θ 4 }) = 0.1 α = 0.4 m({θ 1 }) = 0.6, m({θ 2 }) = 0.1 1.5710 m({θ 3 }) = 0.2, m({θ 4 }) = 0.1 α = 0.9 m({θ 1 }) = 0.6, m({θ 2 }) = 0.1 1.5710 m({θ 3 }) = 0.2, m({θ 4 }) = 0.1 α = 1 m({θ 1 }) = 0.6, m({θ 2 }) = 0.1 1.5710 m({θ 3 }) = 0.2, m({θ 4 }) = 0.1

D. Example 4

Let the FOD be Θ = {θ 1 , θ 2 , θ 3 }. The given FMF is µ(θ 1 ) = 0.6, µ(θ 2 ) = 0.2, µ(θ 3 ) = 0.1. Suppose that m min and m max are the BBAs obtained by using "T min " and "T max ", respectively. We just list the corresponding BBAs for α = 0, 0.3, 0.8 and 1.

This FMF satisfies VIII, when α = 0, the obtained BBAs are identical to m min , and the values of AM are the minimum uncertainty. When α = 1, the obtained BBAs are identical to m max , and the values of AM are the maximal uncertainty. When α → 0, the obtained BBA is similar to m min , and its uncertainty is close to the minimum uncertainty. When α → 1, the obtained BBA is similar to m max , and its uncertainty is close to the maximal uncertainty.

In the Table VII and Table

α BBA AM α = 0 m({θ 1 }) = 0.6, m({θ 3 }) = 0.1 1.2362 m({θ 2 }) = 0.2, m({θ 1 , θ 2 }) = 0.1 α = 0.3 m({θ 1 }) = 0.6, m({θ 1 , θ 2 }) = 0.
α BBA AM α = 0 m({θ 1 }) = 0.6, m({θ 3 }) = 0.1 1.2362 m({θ 2 }) = 0.2, m({θ 1 , θ 2 }) = 0.1 α = 0.3 m({θ 1 }) = 0.6, m({θ 1 , θ 2 }) = 0.0697 1.2748 m({θ 2 }) = 0.2, m({θ 1 , θ 3 }) = 0.0006 m({θ 3 }) = 0.1, m({θ 2 , θ 3 }) = 0.0297 α = 0.8 m({θ 1 }) = 0.6, m({θ 1 , θ 2 }) = 0.02 1.3321 m({θ 2 }) = 0.2, m({θ 2 , θ 3 }) = 0.08 m({θ 3 }) = 0.1 α = 1 m({θ 1 }) = 0.6, m({θ 2 }) = 0.2 1.3527 m({θ 3 }) = 0.1, m({θ 2 , θ 3 }) = 0.1 E. Example 5
Suppose that a system of classification with three sensors, including displacement sensor S 1 , pressure sensor S 2 and image sensor S 3 . Let the FOD be Θ = {θ 1 , θ 2 , θ 3 }. Three sensors are used for measuring the size, weight and state of the sample, respectively. The measurements of sensors are used to obtain two FMFs and a BBA. According to the parameters and the measurements of the sensor, the FMF is defined as

µ(θ i ) =          x -min i ave i -min i , x ∈ [min i , ave i ] x -max i ave i -max i , x ∈ [ave i , max i ] 0, others (23)
where i = 1, 2. min i and max i are the minimum and maximal values of the class θ i , respectively. ave i is the average value of the class θ i . In the Table IX, the parameters of S 1 and S 2 and the measurements of a sample are listed. The class of this sample is θ 2 . According to (23), two FMFs are as follows: S 1 : µ(θ 1 ) = 0.3158, µ(θ 2 ) = 0.5930, µ(θ 3 ) = 0.4049; S 2 : µ(θ 1 ) = 0.6, µ(θ 2 ) = 0.6667, µ(θ 3 ) = 0.3333.

According to the image of S 3 , the expert determined the BBA directly as follows: The BBAs obtained by using "T WA " and "T USO " are listed in the Table X XI. With the increasing value of α, m S1 and m S2 are more close to m max (i.e., the BBA obtained by using "T max " or the BBA obtained when α = 1), which is the reason of the decreasing value of BetP ({θ 2 }) in the Table XIII.

m S3 ({θ 3 })=0.51, m S3 ({θ 2 , θ 3 })=0.38, m S3 (Θ)=0.11.
α = 0 m({θ 1 }) = 0.0276 m({θ 1 }) = 0.0276 m({θ 2 }) = 0.5731 m({θ 2 }) = 0.5731 m({θ 3 }) = 0.3652 m({θ 3 }) = 0.3652 m({θ 1 , θ 3 }) = 0.0340 m({θ 1 , θ 3 }) = 0.0340 α = 0.3 m({θ 1 }) = 0.0429 m({θ 1 }) = 0.0408 m({θ 2 }) = 0.5531 m({θ 2 }) = 0.5547 m({θ 1 , θ 2 }) = 0.0024 m({θ 1 , θ 2 }) = 0.0069 m({θ 3 }) = 0.3637 m({θ 3 }) = 0.375 m({θ 1 , θ 3 }) = 0.0168 m({θ 1 , θ 3 }) = 0.0135 m({θ 2 , θ 3 }) = 0.0212 m({θ 2 , θ 3 }) = 0.0089 m({Θ}) = 0.0002 α = 0.8 m({θ 1 }) = 0.0513 m({θ 1 }) = 0.0418 m({θ 2 }) = 0.5409 m({θ 2 }) = 0.5341 m({θ 1 , θ 2 }) = 0.0089 m({θ 1 , θ 2 }) = 0.0199 m({θ 3 }) = 0.3747 m({θ 3 }) = 0.3780 m({θ 1 , θ 3 }) = 0.0014 m({θ 2 , θ 3 }) = 0.0255 m({θ 2 , θ 3 }) = 0.0228 m({Θ}) = 0.0007 α = 1 m({θ 1 }) = 0.0487 m({θ 1 }) = 0.0487 m({θ 2 }) = 0.5435 m({θ 2 }) = 0.5435 m({θ 1 , θ 2 }) = 0.0124 m({θ 1 , θ 2 }) = 0.0124 m({θ 3 }) = 0.3837 m({θ 3 }) = 0.3837 m({θ 2 , θ 3 }) = 0.0118 m({θ 2 , θ 3 }) = 0.0118

VI. CONCLUSIONS

In this paper, we have proposed two approaches with a user-specified weighting factor to transform a given FMF into a trade-off BBA. These two approaches are both effective approaches for obtaining a trade-off BBA. The users can transform a given FMF into a BBA by their preferred approach. With the cardinality of FOD increasing, the computational complexity of the optimization will become exponential growth. The reason for this is the structure of the belief functions. By using the user-specified weighting factor to influence how close the trade-off BBA is to each of the two BBAs obtained by solving the uncertainty maximization and minimization. The example of using our transformations in the practical application is provided. The numerical examples indicate that the uncertainty of the obtained BBA is between the minimum and maximal uncertainties. In a future work, we will try to use and compare different types of the distance of evidence as objective function to expect a better trade-off BBA.

  This work was supported by the National Natural Science Foundation (Nos. 61573275, 61671370), Postdoctoral Science Foundation of China (No. 2016M592790), Postdoctoral Science Research Foundation of Shaanxi Province (No. 2016BSHEDZZ46) and Fundamental Research Funds for the Central Universities (No. xjj2016066).

  1 }) = 0.6, m({θ 2 }) = 0.1 1.5710 m({θ 3 }) = 0.2, m({θ 4 }) = 0.1 α = 0.4 m({θ 1 }) = 0.6, m({θ 2 }) = 0.1 1.5710 m({θ 3 }) = 0.2, m({θ 4 }) = 0.1 α = 0.9 m({θ 1 }) = 0.6, m({θ 2 }) = 0.1 1.5710 m({θ 3 }) = 0.2, m({θ 4 }) = 0.1 α = 1 m({θ 1 }) = 0.6, m({θ 2 }) = 0.1 1.5710 m({θ 3 }) = 0.2, m({θ 4 }) = 0.1

  07 1.2749 m({θ 2 }) = 0.2, m({θ 2 , θ 3 }) = 0.03 m({θ 3 }) = 0.1 α = 0.8 m({θ 1 }) = 0.6, m({θ 1 , θ 2 }) = 0.02 1.3321 m({θ 2 }) = 0.2, m({θ 2 , θ 3 }) = 0.08 m({θ 3 }) = 0.1 α = 1 m({θ 1 }) = 0.6, m({θ 2 }) = 0.2 1.3527 m({θ 3 }) = 0.1, m({θ 2 , θ 3 }) = 0.1

1 BetP

 1 BetP ({θ 1 }) = 0.0446 BetP ({θ 1 }) = 0.0446 BetP ({θ 2 }) = 0.5731 BetP ({θ 2 }) = 0.5731 BetP ({θ 3 }) = 0.3822 BetP ({θ 3 }) = 0.3822 α = 0.3 BetP ({θ 1 }) = 0.0525 BetP ({θ 1 }) = 0.0511 BetP ({θ 2 }) = 0.5648 BetP ({θ 2 }) = 0.5627 BetP ({θ 3 }) = 0.3827 BetP ({θ 3 }) = 0.3863 α = 0.8 BetP ({θ 1 }) = 0.0565 BetP ({θ 1 }) = 0.0520 BetP ({θ 2 }) = 0.5567 BetP ({θ 2 }) = 0.5570 BetP ({θ 3 }) = 0.3868 BetP ({θ 3 }) = 0.3910 α = ({θ 1 }) = 0.0549 BetP ({θ 1 }) = 0.0549 BetP ({θ 2 }) = 0.5556 BetP ({θ 2 }) = 0.5556 BetP ({θ 3 }) = 0.3895 BetP ({θ 3 }) = 0.3895

TABLE I USING " T

 I" WA " TO OBTAIN BBAS IN EXAMPLE 1.

TABLE III USING " T

 III" WA " TO OBTAIN BBAS IN EXAMPLE 2.

TABLE V USING " T

 V" WA " TO OBTAIN BBAS IN EXAMPLE 3.

TABLE VI USING " T

 VI" USO " TO OBTAIN BBAS IN EXAMPLE 3.

TABLE VII USING " T

 VII" WA " TO OBTAIN BBAS IN EXAMPLE 4.

TABLE VIII USING " T

 VIII" USO " TO OBTAIN BBAS INEXAMPLE 4.

TABLE IX THE

 IX PARAMETERS AND THE MEASUREMENTS OF SENSORS.

	Class	min i	S 1 max i	ave i	min i	S 2 max i	ave i
	θ 1	43.3	58.4	50.8	2.9	4.1	3.4
	θ 2	50.9	70.1	59.5	2.0	3.4	2.8
	θ 3	49.4	79.3	65.7	2.2	3.8	2.9
	Sample		56			3.2

TABLE X USING " T

 X" WA " TO OBTAIN BBAS IN EXAMPLE 5.

	α	m S 1	m S 2

TABLE XIII THE

 XIII PIGNISTIC PROBABILITIES IN EXAMPLE 5.

  and Table XI, respectively. Then, we combine these three BBAs (i.e., m S1 , m S2 and m S3 ). The combined BBA is represented by m. The combined BBAs and the pignistic probabilities are listed in the Table XII and TableXIII, respectively. We just list the corresponding BBAs for α = 0, 0.3, 0.8 and 1.In the Table XIII, all the classification results are θ 2 and are correct. When α is given from 0 to 1, m S1 ({θ 2 }) and m S2 ({θ 2 }) are decreasing in the TableX and Table

i=1 µ(θ i ) = 0.9 < 1. The FMF is equivalent to the singleton belief. The BBAs obtained by using "T WA " and "T USO " are listed in the Table VII and Table VIII, respectively.


domain: spi phys math info

text:
C Musso

B Sacleux

A Bresson

J M Allard

K Dahia

Y Bidel

N Zahzam

C Palmier

Terrain-aided navigation with an atomic gravimeter

Cold atom interferometer is a promising technology to obtain a highly sensitive and accurate absolute gravimeter. With the help of an anomalies gravity map, local measurements of gravity allow a terrain-based navigation. This paper follows the one we published in Fusion 2017. Based on an atomic gravimeter we present a method to map the gravity anomaly. We propose a modification of the Laplace-based particle filter so as to make it more robust. Some simulation results demonstrate a better robustness of the proposed filter.

I. INTRODUCTION

This paper follows the one published in Fusion 2017 [START_REF] Musso | Absolute gravimeter for terrain-aided navigation[END_REF]. The purpose of the paper is to provide a solution for surface or sub-surface navigation by Terrain Matching using an absolute gravimeter. Application of atom interferometry to the inertial navigation sensors has been proposed in [START_REF] Hugh | Submarine navigation applications of atom interferometry[END_REF]. ONERA has developed a compact and robust atom gravimeter, with the support of the French Defense Agency (DGA). Dedicated to measure gravity aboard a boat or a plane, the GIRAFE2 atom gravimeter considered is similar to the one described and used in reference [START_REF] Bidel | Compact cold atom gravimeter for field applications[END_REF].

In section II, we present the improved absolute atom gravimeter GIRAFE2. It is based on the acceleration measurement of a free falling gas of ultra-cold atoms thanks to atom interferometry. In section III, we propose a method for mapping the gravity anomaly. Based on this atomic gravimeter and with the aid of an external sensor delivering the height of the carrier, we develop a method to extract the gravity anomaly. In section IV, we present an improvement of the Laplace-based particle filter for terrain-based navigation [START_REF] Musso | Absolute gravimeter for terrain-aided navigation[END_REF]. It is based on a multimodal importance function with a heavy tail. Finally, in section V, some simulation results demonstrate, in case of high terrain ambiguities, the robustness of the proposed filter compared with the previous one.

II. ABSOLUTE GRAVIMETER BASED ON ATOM

INTERFEROMETRY

The inertial sensor used in numerical simulations provides an absolute measurement of the acceleration (its continuous value and all temporal variations), with the help of a gyrostabilized platform which maintains an atom accelerometer aligned with the gravity direction despite angular movements of the carrier (Fig. 1). The GIRAFE2 atom gravimeter considered is similar to the one described and used in reference [START_REF] Bidel | Compact cold atom gravimeter for field applications[END_REF]. It is based on the acceleration measurement of a free falling gas of ultra-cold atoms thanks to atom interferometry [START_REF] Ch | Atomic interferometry with internal state labelling[END_REF], [START_REF] John | Ultra-high sensitivity accelerometers and gyroscopes using neutral atom matter-wave interferometry[END_REF]. Precisely, this absolute gravimeter is measuring vertical proper acceleration. It is absolute in the sense that it provides absolute values of acceleration of free-falling rubidium atoms with respect to the baseplate of the instrument.

The test mass is produced from a magneto-optical trap loaded from a background vapor. After the trap loading, a stage of optical molasses and a microwave selection, we obtain a cloud of few millions of atoms at a temperature about few µK. The acceleration of the free falling cloud of ultra-cold atoms is then measured by light pulse atom interferometry (Fig. 2). After this sequence, the interference signal is then obtained by measuring the population of atoms in the two hyperfine states corresponding to the two output ports of the interferometer. This measurement is obtained by a standard fluorescence method. The proportion P of atoms in a specific state can be written as

P = P m - C 2 cos 4πT 2 λ Γ (1)
where P m is a constant called offset of the fringe, C is an experimental constant called contrast, λ is the laser wavelength, T the time separation between Raman laser pulses and Γ is the acceleration of the atoms along the direction of the laser beam. This simple relation ( 1) is well verified with no drift or unknown bias at sub mGal level [START_REF] Gillot | Stability comparison of two absolute gravimeters: optical versus atomic interferometers[END_REF] (1mGal = 10 -6 g ≈ 10 -5 m/s 2 ). This new technology has already been used to build absolute static gravimeter and demonstrates better or equal performances than classical absolute gravimeter. Typical stability demonstrated in this kind of instrument is largely below 100 µ Gal/Hz 1/2 in static mode. Recent tests have shown the possibility to operate this next generation of inertial sensors on a moving platform such as aircrafts or marine ship [START_REF] Geiger | Detecting inertial effects with airborne matter-wave interferometry[END_REF], [START_REF] Bidel | Absolute marine gravimetry with matter-wave interferometry[END_REF].

For instance, considering marine ship, due to higher vibration level aboard, absolute measurement at a level of 1mGal / Hz 1/2 could be awaited under calm sea level.

III. A METHOD TO MAP THE GRAVITY ANOMALY WITH THE ATOMIC GRAVIMETER

The absolute gravimeter described in section II provides local gravity measurements. In the earth-centered, earth-fixed coordinate system (ECEF), the highly nonlinear vertical acceleration model is expressed as follows

Γ = -g(φ, λ) + z + γ c (λ) + 2 V × Ω(λ), u z γe→ Eötvös effect (2)
where γ c is the centripetal acceleration

γ c = V 2 E R(λ) -z + V 2 N R(λ)(1-e 2 ) (1-e 2 sin 2 (λ)) -z (3)
We assume an ellipsoidal earth modeling. ., . denotes the inner product and × the cross product. φ is the longitude and λ the latitude. g(φ, λ) is the modulus of the local gravity acceleration vector, V = [V N , V E , V Z = ż] T is the carrier speed and u z is a unit down vertical vector provided by the gyro-stabilized platform (Section II). z is the vertical acceleration due to the carrier movement. Ω(λ) = ω [cos(λ), 0, -sin(λ)]

T is the earth's axis of rotation expressed in the ECEF frame, ω being the earth rotation speed. R(λ) is Earth's radius at latitude λ and e is the earth eccentricity.

The local gravity acceleration g(φ, λ) is expressed as the sum of the standard gravity g 0 and the gravity anomaly g a (φ, λ) (Fig. 3).

       g 0 (λ) =
a g e cos 2 (λ) + b g p sin 2 (λ)

a 2 cos 2 (λ) + b 2 sin 2 (λ) g(φ, λ) = g 0 (λ) + g a (φ, λ) (4)
where a and b are the equatorial and polar semi-axes, and where g e , g p are the gravity accelerations at the equator and poles respectively. For the purpose of navigating, developing gravity anomaly maps is required. The use of a carrier instrumentation to obtain coverage over an extended geographical area is a convenient approach in terms of fast mapping. The aim of this section is to propose a method to estimate the gravity anomaly g a . The difficulty of the experimentation aboard a carrier lies in the spectrum of the acceleration z(t) (2) to which the atomic gravimeter (AG) is subjected. A suitable location of the AG suspension should be not far from the carrier's center of mass. To pick up g a from the measure Γ(t) (2) it is required to estimate z(t) accurately. We propose a method to determine the vertical acceleration z(t) by means of an external measurement of z(t).

A. Mapping the gravity anomaly

The approach described here responds to the scenario where the knowledge of the displacement z according to the local vertical direction is available and accurately known. The vertical acceleration z(t) signal is first determined before extracting g a from the gravimeter's measure Γ(t) (2). We model z(t) as a polynomial of degree n ≥ 2. We estimate the coefficient a = (a 1 , ..., a n+1 ) T of this polynomial by means of the weighted least square method (WLS). Let us introduce the N-vector of the vertical measurements z m (t) up to time t

Z m (t) = [z m (t -(N -1)∆t), z m (t -(N -2)∆t), ..., z m (t)] T (5)
The WLS leads to minimize the following criterion w.r.t a

min a (Z m (t) -Ha) T R (Z m (t) -Ha) (6)
where the matrix R is a positive semi-definite matrix such that the weights of the measures z m (t) are greater around the center of Z m (t) and where

H = [H i,j
] for i = 1, ..., N and j = 1, ..., n + 1

H i,j = (i -1) j-1 ∆t j-1 (7)
The solution of ( 6) is the following

â = K Z m K = (H T RH) -1 H T R (8)
The acceleration zc is estimated at the midpoint

t c = t-∆t c = t -(N +1)∆t 2 as follows      ẑc = G â G = [G 1 , ..., G n+1 ] with G 1 = G 2 = 0 and G i = i(i -1)∆t i-2 c
for i ≥ 2 (9) Indeed, the acceleration ẑc is determined by the second derivative of the fitted polynomial with coefficient â at time t c . The following equation allows us to extract g a g a (t) = zc (t) -g 0 (t) -Γ(t) -γ e (t) + ν(t) [START_REF] Murangira | Proposal distribution for particle filtering applied to terrain navigation[END_REF] where zc (t) is estimated by ( 9), g 0 is the standard gravity (4), γ e is the Eötvös effect (2) and where Γ is provided by the atomic gravimeter (2). ν(t) is a white gaussian noise. The direct extraction of g a , based on this equation, gives an inaccurate estimation. The signal evolution of the gravity anomaly has a very low frequency w.r.t the error signal affecting z. In this context, a low band Kalman filter based on the following 2nd order dynamic is proposed for accurately estimating g a (t).

ga + 2 ξω ġa + ω 2 g a = η (11)
η is a white gaussian noise, ω is the pulsation depending on the carrier velocity and ξ is a damping coefficient. In order to establish the map (such as Fig. 3), it is necessary to know the position and the velocity of the carrier, the standard gravity g 0 (4), and the Eötvös effect γ e (2).

B. Simulation results

As an example, a typical z profile during 900s of a carrier at a speed of 5m/s is illustrated on Fig. 4. The carrier is subjected to an anomaly generated by a second order stochastic process [START_REF] Sanjeev | A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking[END_REF] having a pulsation ω close to 5.10 -5 rd/s. The power spectral density of the noise η is such that the output of this process gives a 30 mGal s.t.d which is typical for anomaly signal variation (Fig. 5). The discretization of the process is done with a sample frequency 1/∆t = 10 Hz (5). The degree of the polynomial fitting z is n = 8. The size of the window ( 5) is 2∆T c = 10.2s (N = 101). The s.t.d. of the noises affecting the measurement of z(t) and the gravimeter outputs Γ(t) are equal respectively to 0.01m and 1 mGal.

Simulations show that the s.t.d of z(t) -ẑ(t) is 383 mGal (Fig. 4). The error of the estimated anomaly obtained by Kalman filtering [START_REF] Sanjeev | A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking[END_REF] shows a s.t.d equal to 0.62 mGal (Fig. 6) with a negligible bias. The estimation error of g a is a colored noise, mainly due to the smoothing method. It is clear that the performance crucially depends on the accuracy of the measurement of z(t).

IV. GRAVITY ANOMALIES FOR TERRAIN-AIDED

NAVIGATION

The measurement equation is based on the model (2). The carrier is here a boat: the vertical acceleration zk is considered negligible when the integration time producing the gravity measurement is long enough. The term γ c will also be ignored as it is negligible in comparison with the anomalies map errors. This leads to the simplified measurement equation where g(φ k , λ k ) is defined in [START_REF] Ch | Atomic interferometry with internal state labelling[END_REF]. k is a Gaussian with s.t.d σ g . It is assumed that the carrier follows a constant heading course (rhumb line) with a constant speed at zero altitude. The state is as follows

Y k = -g(φ k , λ k ) + 2 V × Ω(λ k ), u z + k (12)
X k = [φ k , λ k , V N , V E ] T (13)
where φ k and λ k are the longitude and the latitude at time (k) and where V N and V E are the north and east velocities respectively. The vertical velocity is supposed to be null. The dynamics is non-linear and is given by

   λ k = λ k-1 + ∆T V R cos(K) φ k = φ k-1 + [ϕ(λ k ) -ϕ(λ k-1 )] tan(K) (14)
where ϕ(λ) = log 1+sin(λ) cos(λ)

, where ∆T is the sampling period and K is the heading relative to the north.

In the case of a large uncertainty in the initial state, it is necessary to use a nonlinear filter since the measure (12) is highly nonlinear and multimodal. The particle filter is very convenient and was introduced for the terrain-aided navigation in [START_REF] Bergman | Recursive bayesian estimation[END_REF] and was improved later for example in [START_REF] Murangira | Proposal distribution for particle filtering applied to terrain navigation[END_REF]. A large variety of particle filters (PF) has been proposed (see tutorial [START_REF] Sanjeev | A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking[END_REF]). We recall the basics of the PF and we propose an improvement of the Multimodal Laplace Particle Filter presented in [START_REF] Musso | Absolute gravimeter for terrain-aided navigation[END_REF]. This filter has proved to be effective in a multimodal context for the TBD (Track Before Detect) [START_REF] Bui | Approximation particulaire et méthode de laplace pour le filtrage bayésien[END_REF], [START_REF] Murangira | Nouvelles approches en filtrage particulaire. application au recalage de la navigation inertielle[END_REF], [START_REF] Musso | Improvement of the laplace-based particle filter for track-before-detect[END_REF].

A. Generic particle filter

Consider the following dynamical and measurement model.

X k = f k (X k-1 ) + η k Y k = h k (X k ) + k ( 15
)
where f k and h k being possibly non-linear functions and where k , η k are zero mean Gaussian noises. The aim of the PF is to estimate the posterior density

p k (x) P (X k = x|Y 1:k ) (16)
The PF estimates the state given all the measurements until time (k) 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 pk (x) (17)
The particles X i k evolve according the dynamical model ( 15) and are corrected through the likelihood

Ψ k (x) P (Y k |X k = x) (18)
The following steps summarize the Sampling Importance Resampling algorithm (SIR) with prior proposal [START_REF] Neil J Gordon | Novel approach to nonlinear/non-gaussian bayesian state estimation[END_REF]:

1) Initialisation: sample X i 0 ∼ p(x 0 ) for i = 1, . . . , N and set

ω i 0 = 1 N 2) Prediction: sample X i k ∼ p X k = x|X i k-1
for i = 1, . . . , N 3) Measurement update:

• compute the weights ωi k = ω i k-1 Ψ k (X i k ) • set ω i k = ωi k N j=1 ωj k • compute Xk = N i=1 ω i k X i k 4) Resampling step: discard/multiply particles X i k ac- cording to high/low weights ω i k
This PF algorithm can be significantly improved when considering importance functions (IF) which take account the current measurement such as the Laplace IF. We recall briefly this algorithm, details are provided in [START_REF] Musso | A laplacebased particle filter for track-before-detect[END_REF], [START_REF] Musso | Absolute gravimeter for terrain-aided navigation[END_REF].

B. Laplace-based importance function

We ignore for a while the time index (k). Consider an unknown d-dimensional state x distributed according to a prior q and observed through a measurement

y = h(x) + (19)
h being a non-linear function and a zero mean Gaussian noise. The prior q is not expressed in closed form, only a sample from q is available (it is the case in particle filtering). We suppose in this section that the posterior has a predominant mode. We recall the Laplace method [START_REF] Bui | Multidimensional laplace formulas for nonlinear bayesian estimation[END_REF] which is useful to design an efficient proposal function. The likelihood p(y|x) is denoted by Ψ(x), the posterior is written as

P (X = x|Y = y) p(x|y) ∝ Ψ(x) q(x) (20)
We wish to compute an estimate of the posterior. The importance sampling (IS) estimator p of p(x|y) is obtained by drawing N samples X i from a proposal distribution q so that p(x|y)

≈ N i=1 w i δ x=X i p(x) (21)
where

w i = wi N i=1
wi with the following importance weights

wi = Ψ(X i )q(X i ) q(X i ) (22)
The choice of the proposal distribution is crucial for controlling the Monte Carlo error which can be large if the prior and the likelihood have little overlap (Fig. 7). This is the case for example if the gravity anomalies measurements are very accurate (which is the case with the atom gravimeter).

In this case, the information contained in the measurement function is very high and most samples fall in a region where the likelihood is low. It is well-known that the optimal importance function is the posterior qopt (x) = p(x|y). For this purpose, we choose a proposal, for instance a Gaussian, which has moments nearly equal to those of the posterior. The posterior expectation E[X|Y ] and the posterior covariance matrix V[X|Y ] are well approximated by the Laplace formula if the posterior has a predominant mode. The following approximations are in general sufficient [START_REF] Musso | Improvement of the laplace-based particle filter for track-before-detect[END_REF] Fig. 7. Small overlap between the prior (red) and the narrow likelihood (blue). Posterior (green)

E[X|Y ] ≈ x V[X|Y ] ≈ Ĵ-1 ( 23
)
where x is the maximum a posteriori (MAP)

x = arg max

x∈R d {Ψ(x) q(x)} (24)
Assuming that the prior q(x) is Gaussian the posterior information matrix J can be approximated by

Ĵ ≈ -(log Ψ) (x) + P -1 ( 25
)
where P is the covariance matrix of the prior q. The proposal density q is chosen so that it has x as mean and Ĵ-1 as covariance matrix

q ∼ q (x, Ĵ-1 ) (26)
Of course, q must be chosen among densities that are easy to sample from. A good criterion for evaluating the accuracy of the approximation (21) is the asymptotic variance of the unnormalized weights wi [18]

V(q) ≈ 1 N   R d Ψ(x) 2 q(x) 2 q(x) dx ( R d Ψ(x) q(x)dx) 2 -1   ( 27
)
It is suitable for robustness reasons to chose q close to the posterior with heavy tails. Indeed, if, for example, q and q are Gaussian in (27), q 2 /q can be unbounded if the s.t.d of q is too small. Therefore, in the sequel, we take q as a multivariate Student distribution.

t d (x; ν, µ, Σ) = Γ( ν+d 2 )(πν) -ν 2 |Σ| -1/2 Γ( ν 2 ) [1 + (x -µ) T Σ -1 (x -µ)/ν]
ν+d 2

(28) The smaller the value of ν, the heavier tail is.

C. Approximation of the MAP adapted to the gravity measurement model

In most filtering problems, only a small part of the state vector X k = [x 1,k , x 2,k ] T contributes non linearly in the measurement equation:

h k (X k ) = h k (x 1,k , x 2,k ) = h 1 (x 1,k ) + M (x 1,k ) x 2,k (29)
This is the case in our context where the latitude and the longitude (φ k , λ k ) are the nonlinear contribution in the measurement equation [START_REF] Bui | Approximation particulaire et méthode de laplace pour le filtrage bayésien[END_REF]. The velocities (V N , V E ) contribute linearly in the measurement equation. Splitting the state as

X k = [ φ k , λ k x 1,k , V N , V E x 2,k ] T
the measurement equation ( 12) can be expressed as follows:

Y k = -g(φ k , λ k ) + 2 V × Ω(λ k ), u z + k = h 1 (x 1,k ) + M (x 1,k ) x 2,k + k = h k (X k ) + k where M (x 1,k ) = -2 ω u T z   0 sin(λ k ) -sin(λ k ) 0 0 cos(λ k )   ( 30
)
We aim to approximate the MAP xk defined as:

xk = arg max x∈R 4 {Ψ k (x) q k (x)} (31)
where Ψ k the likelihood [START_REF] Del | Genealogical and interacting particle systems with applications, probability and its applications[END_REF] and the q k is the prior [START_REF] David T Sandwell | New global marine gravity model from cryosat-2 and jason-1 reveals buried tectonic structure[END_REF].

The nonlinear contribution can be handled separately in the maximization (31) under the Gaussian assumption on the prior q k . This leads to the following algorithm [START_REF] Musso | Improvement of the laplace-based particle filter for track-before-detect[END_REF] which approximates the MAP.

1) Minimize the conditional criterion with respect to x

1 x1 = arg min x1∈R 2 G(x 1 ) 2) Evaluate x2 x2 = γ(x 1 ) (32)
with x 1 ∈ R 2 and where G(x 1 ), described in [START_REF] Musso | Absolute gravimeter for terrain-aided navigation[END_REF], depends entirely on x 1 and γ is an explicit function of x 1 . The prior q k (x) = q k (x 1 , x 2 ) is assumed to be Gaussian. The MAP (31) is simply estimated by xk = [x 1 , x2 ] T . Finally, the initial 4dimensional maximization (24) to get the MAP boils down to a 2-dimensional minimization.

D. Improvement of the MM Laplace Particle Filter

In order to address the problem of multimodality due to ambiguities of the gravity anomalies map, we can use a multimodal importance function. At time (k) the posterior has the following form

p(X k = x k |Y k ) = p(x 1,k , x 2,k |Y k ) ∝ Ψ k (x 1,k , x 2,k ) q k (x 1,k , x 2,k ) (33)
This posterior can be multimodal, so we need to compute m local MAP xj k for j = 1, • • • , m of the posterior (33). This is done in 2 steps as described in (32). We maximize first the 2dimensional criterion G(x 1,k ) = G(φ k , λ k ) over the latitudes and the longitudes to get m local MAP xj 1,k . The domain of this maximization is restricted to the confidence ellipse (uniform mesh) given by the predicted covariance matrix P k . P k is computed empirically with the predicted particles. Then, we compute (32) xj 2,k = γ(x j 1,k ) (which corresponds to the north and east velocities) to get the full components of the local MAP

xj k = (x j 1,k , xj 2,k ).

The proposed importance function q2

k consists in a mixture of m Student distributions (28) centered on the local MAP

           q1 k (x) = m j=1 ρ j k φ x, xj k , [ Ĵj k ] -1 q2 k (x) = 1 m m j=1 t d x, xj k , Σ k (34)
where Σ k = ν-2 ν P k such that the covariance of t d is equal to P k and where q1 k is defined in [START_REF] Musso | Absolute gravimeter for terrain-aided navigation[END_REF] with φ being a Gaussian. Laplace-based resampling, i.e generating particles Xi k according to qj k (x) (34), is performed only in case of degeneracy of the weights w i k . It replaces the traditional (multinomial) resampling. This degeneracy is detected by the effective sample size of particles

N ef f = 1 N i=1 [w i k ] 2
when it is less than a threshold N th [START_REF] Kong | Sequential imputations and bayesian missing data problems[END_REF]. The evaluation of q k ( Xi k ) needed for computing the weights ( 22) is done by taking the prior equal to a Student distribution (28) centered on the mean of the predicted particles µ k with Σ k = ν-2 ν P k . The algorithm of the proposed PF is described below. 1) Initialisation (k=1) For i = 1, . . . , N . Generate the particles

X i k-1 = (x i 1,k-1 , x i 2,k-1
) (positions and velocities) according to the prior. with w i k-1 ≡ 1/N

2) Prediction For i = 1, . . . , N . Propagate the particles by applying the carrier dynamics ( 14)

X i k|k-1 = f k (X i k-1 ) + η i k
3) Correction For i = 1, . . . , N . Compute the likelihood (12) Ψ k (X i k|k-1 ) = P (Y k |X i k|k-1 ) and the weights

w i k ∝ Ψ k (X i k|k-1 ) w i k-1 such that N i=1 w i k = 1. Compute N ef f = 1 N i=1 [w i k ] 2 * If N ef f ≥ N th . The corrected particles are X i k , w i k = X i k|k-1 , w i k * If N ef f < N th then
∝ Ψ k ( Xi k ) q k ( Xi k ) q2 k ( Xi k ) such that N i=1 wi k = 1 with the prior q k ∼ t d (., µ k , Σ k ) The corrected particles are X i k , w i k = Xi k , wi k 4) State estimation The state is estimated by Xk = N i=1 w i k X i k Go to the prediction step k → k + 1 V. SIMULATION RESULTS
We evaluate the performances of the proposed algorithm compared to the original one [START_REF] Musso | Absolute gravimeter for terrain-aided navigation[END_REF]. For that purpose, we simulate 2 trajectories generating anomaly gravity signals which can be problematic for a particle filter (Fig. 10). The anomaly gravity map (Fig. 8) covers the zone [-9

• -4 • ]x[46 • 49 • ]
in the longitude-latitude plan. This map is extracted from the map presented in [START_REF] David T Sandwell | New global marine gravity model from cryosat-2 and jason-1 reveals buried tectonic structure[END_REF] and has spatial resolution of 12 arc second (which corresponds to 400m at this latitude). Below, we present the characteristics of the 2 scenarios and of the parameters of the 2 filters.

Scenario 1

The zone where the trajectory 1 starts presents severe ambiguities: many carrier trajectories starting from the uncertainty zone collect a similar anomaly measurements history (taking account the measurement noise) for the first 20 iterations (Fig. 9). At the beginning of this trajectory, there is a strong variation of the gravity anomaly (Fig. 10) which provides a great (local) information to the filter. But, in this context where the ambiguity is strong and where the initial uncertainty zone is large, we expect the filter to have a poor behavior. Indeed, we have observed 25 divergences for algorithm 1 and 11 divergences for algorithm 2. The latter shows a better robustness in this difficult quasi unobservable context. The improved filter can estimate temporarily a wrong trajectory but, due to the heavy tail of the IF, it can join the true one. We can see (Fig. 9) that some particles generated by the IF act as a kind of "trailblazer". The support (99.9 confidence ellipse) of the mesh is enlarged. The RMSE for the 4 components of the state vector is computed for the 2 algorithms only on the convergent trials (Fig. 11 &12). During the first 30 iterations the anomaly gravity signal varies widely (Fig. 10) which leads to a fast convergence of the filters. Then the signal varies slowly and the estimation does not improve. We observe that the latitude is better estimated. This is due to the fact that the Eötvös effect, which depends on the latitude provides acceleration information (2) [START_REF] Musso | Absolute gravimeter for terrain-aided navigation[END_REF]. The Laplace-based resampling rate for algorithm 1 is about 21% and 15% for algorithm 2. For the 60 iterations, the computational cost is 13 s for the algorithm 1 and 9 s for the algorithm 2 on a 2.5 GHz Intel Core I5. The trajectory 2 (Fig. 8) shows a slowly varying anomaly gravity signal (Fig. 10). However, the zone is less ambiguous that the one of the scenario 1. For the two algorithms we have observed 2% of divergences. The RMSE of these algorithms are comparable (Fig. 13 & 14) for the 4 components. The Laplace-based resampling (section ??) rate for algorithm 1 is about 25% and 10% for algorithm 2. For the 60 iterations, the computational cost is 11 s for algorithm 1 and 6 s for algorithm 2.

The algorithm 2 offers a better robustness with a slightly better RMSE. Moreover, as the Laplace-redistributions are less frequent for this algorithm, the computational cost is reduced.

Fig. 1 .

 1 Fig. 1. Cold atom gravimeter prototype (GIRAFE2), aboard Norlandair Twin-Otter during an airborne gravimetric survey.

Fig. 2 .

 2 Fig. 2. Spatio-temporal diagram of the atomic free-fall. Output signal is simply proportional to the sinus of the absolute acceleration with a welldefined period.

Fig. 3 .

 3 Fig. 3. Gravity anomalies. North atlantic map. Bureau Gravimétrique International (BGI). View of the mid-Atlantic ridge

Fig. 4 .Fig. 5 .

 45 Fig. 4. Time history of zt

Fig. 6 .

 6 Fig. 6. Time history of the gravity anomaly error after Kalman filtering

Fig. 8 .•

 8 Fig. 8. Anomalies gravity map with the carrier's trajectories

Fig. 9 .

 9 Fig. 9. illustration of the behavior of the heavy-tailed importance function (iterations 7 and 8)

Fig. 10 .

 10 Fig. 10. Time history of the gravity anomaly for the 2 scenarios

Fig. 11 .

 11 Fig. 11. Scenario 1. Position RMSE comparison of algorithm 1 (left figure) with algorithm 2 (right figure).

Fig. 12 .

 12 Fig. 12. Scenario 1. Velocity RMSE comparison of algorithm 1 (left figure) with algorithm 2 (right figure).

Fig. 13 .

 13 Fig. 13. Scenario 2. Position RMSE comparison of algorithm 1 (left figure) with algorithm 2 (right figure).

Fig. 14 .

 14 Fig. 14. Scenario 2. Velocity RMSE comparison of algorithm 1 (left figure) with algorithm 2 (right figure).VI. CONCLUSIONWe present the new version of the absolute gravimeter based on atom interferometry developed by ONERA. Based on it and on external measurements of the height of the carrier, we propose a method to map the gravity anomaly. The modified Laplace-based particle filter proposed for terrainbased navigation shows a better robustness in difficult contexts.

  perform Laplace-based resampling (34). -Compute the local modes of the posterior (32) : xj k = [x j 1,k , xj 2,k ] T for j = 1, . . . , m -Generate samples Xi k from the proposal q2 k (34) and compute the importance weights wi k