text:
Yilin Dong

Xinde Li
email: xindeli@seu.edu.cn

Jean Dezert
email: jean.dezert@onera.fr

Combination of Sources of Evidence with Distinct Frames of Discernment

Keywords: Information fusion, Distinct frames of discernment, Belief functions, Combination Rule, Granular Structure, Rough set

Multi-source information fusion strategies in target recognition have been widely applied. Generally, each source is defined and modelled over a common frame composed of the hypotheses to discern. However, in practice, the independent sources of evidence can refer to distinct frames of discernment in terms of the hypotheses they consider. Under this condition, the classical combination process cannot be applied directly. Working with distinct frames of discernment for information fusion is a problem often encountered in the development of recognition systems which requires a particular attention. In order to combine such sources, this paper presents a new combination method which splits the process of fusion into two steps: construction of granular structure, calculation of belief mass, followed by the fusion process. Our simulations results show that the proposed method can effectively solve the problem of fusion of sources defined on distinct frames.

I. INTRODUCTION

The aim of this paper is to propose a new approach for combining sources of evidence referring to distinct (uncommon) frame of discernment, that can be eventually (partially) incompatible.

Motivation

In Dempster-Shafer Theory (DST) also known as the theory of evidence [START_REF] Shafer | A mathematical theory of evidence[END_REF], each Source of Evidence (SoE) refers to a (common) Frame of Discernment (FoD) with exhaustive and exclusive elements. Each element of the FoD corresponds to a choice for the solution of the problem under concern. In some applications, the FoD related with a source is not always exhaustive and exclusive and that is why Dezert-Smarandache Theory (DSmT) [START_REF] Smarandache | Advances and applications of DSmT for information fusion[END_REF] has been developed. Moreover, in more realistic situations, the FoDs related with each source can be distinct so that no global common FoD is directly available in order to apply classical fusion techniques developed in DST or DSmT. The fusion of sources of evidence related to distinct FoDs encountered in real applications motivates the works presented in this paper.

State of the art and challenges

In order to fuse sources of evidence referring to distinct FoDs, some scholars have proposed different methods. The most representative ones are presented in [START_REF] Janez | Theory of evidence and non-exhaustive frames of discernment: Plausibilities correction methods[END_REF], [START_REF] Kulasekere | Conditioning and updating evidence[END_REF], [START_REF] Premaratne | Evidence combination in an environment with heterogeneous sources[END_REF] and [START_REF] Wickramarathne | Belief theoretic methods for soft and hard data fusion[END_REF]. These methods can be classified into three categories as follows.

1) Expanded methods: In these methods, for example in [START_REF] Janez | Theory of evidence and non-exhaustive frames of discernment: Plausibilities correction methods[END_REF], a common (expanded) FoD is built and each SoE assigns a zero mass value to elements of this common FoD that are not included in its own specific FoD. This kind of method is simple and effective up to some degree. The implicit assumption done in these methods is that each source of evidence can identify the expanded common FoD as being the disjunction of all own specific FoDs of the sources. This implicit assumption is however unreasonable because it overlooks the main fact that the elements of the expanded common FoD are not within the scope of expertise of each source under concern.

2) Deconditioning methods: To overcome the problem of expanded methods, Janez [START_REF] Janez | Theory of evidence and non-exhaustive frames of discernment: Plausibilities correction methods[END_REF] adopted a new standpoint, and he developed a deconditioning fusion method. This method can be interpreted as the opposite (dual) of the extension method. The so-called discriminable domain is divided according to the FoD of each source of evidence. For example, we consider two SoE S 1 and S 2 with FoDs Θ 1 and Θ 2 and let1 Θ c Θ 1 ∩ Θ 2 be the common part of Θ 1 and Θ 2 (assumed to be non empty). The elements of Θ 1 ∩ Θc will be discriminated only by source S 1 , whereas the elements of Θ 2 ∩ Θc will be discriminated only by the S 2 . The elements of Θ c will be discriminated by both sources S 1 and S 2 . Such kind of separation is based on the range of each evidence, which may yield to some ambiguities of the extension methods. This separation principle splits the original fusion process into three sub-steps, in which one needs to make the fusion. Hence, with these methods one can degrade the computational efficiency of the fusion strategy, especially if large databases have to be browsed.

3) Conditioning updating methods: Unlike the previous methods, the more recent conditional method proposed by Premaratne et al. in [START_REF] Kulasekere | Conditioning and updating evidence[END_REF], [START_REF] Premaratne | Evidence combination in an environment with heterogeneous sources[END_REF] and [START_REF] Wickramarathne | Belief theoretic methods for soft and hard data fusion[END_REF], is a different approach, which does not need to expand the FoDs of each source into a common FoD, nor to separate and identify the common elements of FoDs. In this new (conditioning and updating) method each source updates its knowledge base by selecting the incoming evidence that it can discern (or it is interested in). The performance of such method is however very dependent on the ad-hoc choice of two weighting tuning parameters.

Contribution

As aforementioned discussed, the existing methods have all their own limitations and drawbacks. In this paper, we propose to combine the granular computing developed by Pawlak [START_REF] Pawlak | Rough sets: Theoretical aspects of reasoning about data[END_REF] in his Rough Set Theory (RST) and belief mass allocation in Belief Functions (BFs). Specifically, we split the process of fusion into two steps: 1) the construction of granular structure (GS) based on the defined relations, and 2) the calculation of belief mass of target focal element (TFE). The innovation of GS construction is that at first we determine the indistinguishability relation according to TFE. Then we construct the relevant GS about the indistinguishability relation in the framework of the available FoDs. The calculation of the belief mass of TFE is the summation of the related masses of TFE in each GS with conflict redistribution, followed by a normalization step. This paper is organized as follows. Section 2 reviews some basic concepts of belief functions. Section 3 discusses fusion strategy based on RST for SoE defined on distinct FoDs. Section 4 gives the summary of the new proposed combination method. In section 5, a target identification case is chosen to show the performance of this new approach. Section 6 concludes the paper and proposes some perspectives.

II. BASICS OF BELIEF FUNCTIONS

In DST framework [START_REF] Shafer | A mathematical theory of evidence[END_REF], the FoD Θ {θ 1 , . . . , θ n } (n ≥ 2) is a set of exhaustive and exclusive elements. This is called Shafer's model of FoD. The set of all subsets of Θ is called the powerset of Θ and it is denoted 2 Θ . A basic belief assignment (BBA) m(•) is defined as a mapping m(•) : 2 Θ → [0, 1], verifying m(∅) = 0 and A∈2 Θ m(A) = 1. In DSmT [START_REF] Smarandache | Advances and applications of DSmT for information fusion[END_REF], one can abandon Shafer's model2 and refute the principle of the third excluded middle. Instead of defining the BBAs on the power set 2 Θ (Θ, ∪), the BBAs are defined on the so-called hyper-power set (or Dedekind's lattice) denoted D Θ (Θ, ∪, ∩) whose cardinalities follows Dedekind's numbers sequence [START_REF] Smarandache | Advances and applications of DSmT for information fusion[END_REF], Vol.1. A (generalized) BBA is defined as the mapping m(•) :

D Θ → [0, 1], verifying m(∅) = 0 and A∈D Θ m(A) = 1. The DSmT framework encompasses DST framework because 2 Θ ⊂ D Θ . If m(A) > 0, A is called focal element of the BBA m(•).
A BBA is called a Bayesian BBA if all of its focal elements are singletons and Shafer's model is assumed, otherwise it is called non-Bayesian. A full ignorance source is represented by the vacuous BBA m v (Θ) = 1. In DST, the combination of several distinct SoE is done with the controversial Dempster's rule of combination which is the normalized version of conjunctive rule [START_REF] Shafer | A mathematical theory of evidence[END_REF]. To palliate the drawbacks of Dempster's rule, the Proportional Conflict Redistribution Rule # 6 (PCR6) has been proposed in DSmT. Because we do not use these rules in this work, we do not need to present them in details here.

III. FUSION OF EVIDENCES BASED ON RST A. Construction of GS of TFE

Consider a FoD Θ = {θ 1 , θ 2 , • • • , θ n }, (n ≥ 2).
We denote X a specific TFE under concern. We denote U = {u 1 , u 2 , • • • , u i } a chosen subset of D Θ corresponding to the set of focal elements of the BBA m(•) of the SoE. According to Pawlak's RST [START_REF] Pawlak | Rough sets: Theoretical aspects of reasoning about data[END_REF], any set U can be partitioned from some chosen rules which help to determine if two elements of U can be considered as equivalent or relevant, or not. For example, if the chosen rule is the belief mass, we will say that u i ∈ U and u j ∈ U are equivalent if and only if m(u i ) = m(u j ). Using RST notations and formalism [START_REF] Pawlak | Rough sets: Theoretical aspects of reasoning about data[END_REF], the relation IN D(R) defined by ( 1) is called a R-indiscernibility relation.

IN D(R) {(u i , u j ) ∈ U × U | ∀r ∈ R, r(u i ) = r(u j )} (1) If (u i , u j ) ∈ IN D(R), then u i and u j are indiscernible (or indistinguishable) from R. We denote U (R) {[u i ] R |u i ∈ U } the partition of U based on R. [u i ] R denotes the equivalent class of u i with respect to R, i.e., [u i ] R = {u j ∈ U |r(u j ) = r(u i )}. The partition U (R) is interpreted as a GS built from equivalence relation R.
In fact, R is a partitioning method which can discriminate elements of U into separate subsets. With different equivalence relations, we can construct various meaningful partitions (GSs) of U .

R : u i ≡ u j if and only if m(u i ) = m(u j )
we see that U can be partitioned as follows:

U (R) = {{u 1 , u 3 }, {u 2 , u 4 , u 5 }, {u 6 , u 7 , u 8 }, {u 9 , u 10 }}
Of course, if a different rule is adopted, an other GS of U is generated. In fusion problems, one has to consider several initial FoDs which are directly provided by the SoE. Each FoD associated with a source is related to the FoD corresponding to the scope of expertise of the SoE.

Example 2: Consider three distinct FoDs Θ 1 = {x 1 , x 2 , x 3 }, Θ 2 = {y 1 , y 2 } and Θ 3 = {z 1 , z 2 }
representing radar cross section (RCS) attributes, shape attributes and speed attributes, respectively. Suppose that one observes a target T among a family of potential targets Θ T {T 1 , T 2 , T 3 , T 4 , T 5 , T 6 } and each target is characterized by its 3-tuple (x i , y j , z k ). For example, T 1 = (x 1 , y 1 , z 1 ), T 2 = (x 1 , y 2 , z 1 ), T 3 = (x 2 , y 2 , z 2 ), T 4 = (x 3 , y 2 , z 1 ), T 5 = (x 3 , y 1 , z 1 ), and T 6 = (x 2 , y 1 , z 2 ). Suppose that three smart sensors (i.e. sources of evidence) provide three reports expressed by the following BBAs

m 1 (x 1 ) = 0.1, m 1 (x 3 ) = 0.2, m 1 (x 1 ∪ x 2 ) = 0.3, m 1 (x 1 ∪ x 2 ∪ x 3 ) = 0.4 m 2 (y 1 ) = 0.2, m 2 (y 2 ) = 0.3, m 2 (y 1 ∪ y 2 ) = 0.5 m 3 (z 1 ) = 0.8, m 3 (z 2 ) = 0, m 3 (z 1 ∪ z 2 ) = 0.2.
From these three SoE, one gets the three FoDs (i.e. sets of focal elements) U 1 , U 2 and U 3

U 1 = {x 1 , x 3 , x 1 ∪ x 2 , x 1 ∪ x 2 ∪ x 3 }; U 2 = {y 1 , y 2 , y 1 ∪ y 2 }; U 3 = {z 1 , z 2 , z 1 ∪ z 2 }.
These three sets of focal elements can be understood as three different GSs obtained by three independent experts using their prior knowledge in order to partition the set

U 123 U 1 × U 2 × U 3 .
In order to compute the global (fusioned) belief mass of any TFE X from the available FoDs U 1 , U 2 and U 3 , we propose two works with two distinct domains called: the compatible domain and the incompatible domain of X, whose definitions are given below.

Definition 1 (Compatible domain of X): Consider q ≥ 2 focal sets U 1 , U 2 , . . ., U q provided by q SoE. The compatible domain of a TFE X = (x i1 , . . . , x iq ), denoted U C 1,...,q (X), corresponds to the set of all q-tuples of U 1 ×. . .×U q and each componentwise intersection of q-tuple of U 1 ×U 2 ×U q with X must be equal to the component of X which means that the qtuples of U 1 × . . . × U q is compatible with X. Mathematically, the compatible-domain of X is defined by

U C 1,...,q (X) {(u i1 , . . . , u iq ) ∈ U 1 × . . . × U q |∀u i k ∈ U k , u i k ∩ x i k = x i k , x i k ∈ X, k = 1, . . . , q}. (2)
For instance, taking back the example 2 and considering the TFE T 1 = (x 1 , y 1 , z 1 ), all 3-tuples that satisfy definition 1 are listed in Table I.

Interpreting the conflict as a situation which cannot occur (i.e. an impossible case), we define the incompatible domain of a given TFE X as the set of all elements of the product space U 1 × U 2 × . . . × U q which have at least one component incompatible (or in conflict) with X. This yields to the following definition: Definition 2 (Incompatible domain with X): Consider a TFE X = (x i1 , . . . , x iq ), and q ≥ 2 focal sets U 1 , U 2 , . . ., U q provided by q SoE. The incompatible domain with X (or conflict domain) denoted U I(X) 1,...,q (∅) corresponds to the set of all q-tuples of U 1 × . . . × U q having at least a componentwise empty intersection with X, that is

U I(X) 1,...,q (∅) {(u i1 , . . . , u iq ) ∈ U 1 × . . . × U q | ∃u i k ∈ U k , u i k ∩ x i k = ∅}. ( 3
)
Example 2 (revisited):

We show how to build these two domains for Example 2. As mentioned before the set of focal elements available from the three SoE are

U 1 = {x 1 , x 3 , x 1 ∪ x 2 , x 1 ∪ x 2 ∪ x 3 }; U 2 = {y 1 , y 2 , y 1 ∪ y 2 }; U 3 = {z 1 , z 2 , z 1 ∪ z 2 }.
Because U 1 has 4 elements, U 2 has 3 elements and U 3 has 3 elements, the Cartesian product U 1 ×U 2 ×U 3 has 4×3×3 = 36 elements. We suppose that the TFE under concern is X = T 1 = (x 1 , y 1 , z 1 ). Based on the definitions 1 and 2, the self domain and the conflict domain of T 1 , are given as follows:

• Compatible domain of T 1 : All the elements of U C 1,2,3 (T 1 ) are listed in Table I. • Incompatible domain of T 1 : The set U I(T1) 1,2,3 (∅) is the complement of U C 1,2,3 (T 1 ) in U 1 × U 2 × U 3 , that is U I(T1) 1,2,3 (∅) = U 1 × U 2 × U 3 -U C 1,2,3 (T 1 )
, where the minus operator symbol corresponds to the set difference operator. U I(T1) 1,2,3 (∅) has therefore 24 = 36 -12 elements (not listed here due to space limitation constraint). The product set U 1 × U 2 × U 3 has been partitioned into two subsets, the compatible domain of T 1 including the twelve elements of U 1 × U 2 × U 3 compatible with T 1 , and the incompatible domain of T 1 including the 24 elements incompatible with T 1 .

Link with RST: From the standpoint of RST, we can define some equivalence rules based on definitions 1 and 2, denoted R = {r 1 = rule based on definition 1, r 2 = rule based on definition 2}. Under a chosen rule r ∈ R, some particular elements (u i1 , . . . , u iq ) of U 1 × . . . × U q can be selected to characterize some special relationship with TFE (based on equivalence relation using definition 1 and conflict relation using definition 2). According to the principles of RST: GSs3 can be induced by various relations depending on the specificity of the problem, such as a neighborhood relation or a tolerance relation [START_REF] Pawlak | Rough sets: Theoretical aspects of reasoning about data[END_REF]. In the following example, we will show how a tolerance relation is used to solve the case where we have incomplete information of targets.

Dealing with incomplete information: In many practical problems, the available information can be incomplete. For example, it could happen that an object under interest has a (or more) missing (non observed) attribute. In such situation, this requires to deal with incomplete information with a special notation using the question mark symbol ? to represents the missing information. For instance, if the third attribute for the target T 1 of example 2 is missing (unknown), we will denote it as (x 1 , y 1 , ?) and question mark symbol ? represents the third missing attribute. Because of the possibility of missing attribute, the definition 1 needs to be revised as follows 4 :

TABLE I: Elements of U C 1,2,3 (T 1 ) for Example 2. (x 1 , y 1 , z 1 ) (x 1 , y 1 ∪ y 2 , z 1 ) (x 1 , y 1 , z 1 ∪ z 2 ) (x 1 , y 1 ∪ y 2 , z 1 ∪ z 2 ) (x 1 ∪ x 2 , y 1 , z 1 ) (x 1 ∪ x 2 , y 1 ∪ y 2 , z 1 ) (x 1 ∪ x 2 , y 1 , z 1 ∪ z 2 ) (x 1 ∪ x 2 , y 1 ∪ y 2 , z 1 ∪ z 2 ) (x 1 ∪ x 2 ∪ x 3 , y 1 , z 1 ) (x 1 ∪ x 2 ∪ x 3 , y 1 ∪ y 2 , z 1 ) (x 1 ∪ x 2 ∪ x 3 , y 1 , z 1 ∪ z 2 ) (x 1 ∪ x 2 ∪ x 3 , y 1 ∪ y 2 , z 1 ∪ z 2 )
Definition 3 (Compatible domain with incomplete information): Consider q ≥ 2 focal sets U 1 , U 2 , . . ., U q provided by q SoE. The compatible domain of a TFE X with missing attributes, denoted U C 1,...,q (X), corresponds to the set of all qtuples of U 1 × . . . × U q and each componentwise intersection of q-tuple of U 1 × U 2 × U q with X must be equal to the component of X, that is

U C 1,...,q (X) {(u i1 , . . . , u iq ) ∈ U 1 × . . . × U q |∀u i k ∈ U k , u i k ∩ x i k = x i k or u i k →?, k = 1, . . . , q}. (4)
Here, u i k →? means that for all u i k in U k replacing the missing attribute in k-th dimension of X, respectively. And

u i k →? actually means u i k ∩ ? ≡ u i k , k = 1, . . . , q.
Considering back the example 2, if X = (x 1 , y 1 , ?), we will get U C 1,...,3 (X) as follows:

U C 1,...,3 (X) = {(x 1 , y 1 , z 1 ), (x 1 , y 1 ∪ y 2 , z 1 ), (x 1 , y 1 , z 1 ∪ z 2 ), (x 1 , y 1 ∪ y 2 , z 1 ∪ z 2 ), (x 1 ∪ x 2 , y 1 , z 1 ), (x 1 ∪ x 2 , y 1 ∪ y 2 , z 1 ), (x 1 ∪ x 2 , y 1 , z 1 ∪ z 2 ), (x 1 ∪ x 2 , y 1 ∪ y 2 , z 1 ∪ z 2 ), (x 1 ∪ x 2 ∪ x 3 , y 1 , z 1 ), (x 1 ∪ x 2 ∪ x 3 , y 1 ∪ y 2 , z 1 ), (x 1 ∪ x 2 ∪ x 3 , y 1 , z 1 ∪ z 2 ), (x 1 ∪ x 2 ∪ x 3 , y 1 ∪ y 2 , z 1 ∪ z 2 ), (x 1 , y 1 , z 2 ), (x 1 ∪ x 2 , y 1 , z 2 ), (x 1 ∪ x 2 ∪ x 3 , y 1 , z 2 ), (x 1 , y 1 ∪ y 2 , z 2 ), (x 1 ∪ x 2 , y 1 ∪ y 2 , z 2 ), (x 1 ∪ x 2 ∪ x 3 , y 1 ∪ y 2 , z 2 )}
The elements of U C 1,...,3 (X) are the triplets of U 1 × U 2 × U 3 compatible with the two first components of X for which the missing attribute of X is replaced either by z 1 , z 2 , or z 1 ∪ z 2 .

The definition 3 can also be regarded as a kind of equivalence relation 5 which is applied to find the relevant sets of TFE in U 1 × . . . × U q .

The main contribution of this paper is to propose some rules to handle the different domains of TFE based on the specificity of the problem. Based on these rules, we can identify the compatible-domain and the incompatible domain for a chosen TFE, and then define belief allocation strategies and combination in order to get a final quantitative result as proposed in the sequel.

B. Calculation of the belief mass of TFE

Once the compatible domain and the incompatible domain of TFE are identified, we use the elements in these domains to calculate the belief mass of a chosen TFE. In fact, in RST, the granular space obtained by each partitioning rule will be used to approximate the target concepts, namely the upper approximation set and the lower approximation set, as well as the positive, negative and uncertainty domains about the target concept. Such results of approximation are to some extent qualitative expressions. Once equivalent sets of TFE in each domain are found, we calculate the belief mass of TFE in each domain space in order to get a quantitative result. In a certain sense, rough set solves the qualitative solution of target concept and the multi-granular rough set fusion is a kind of qualitative fusion method, whereas the evidence fusion provides a quantitative fusion solution. Although the two approaches are focused on different aspects, from a fusion point of view, the rough set approach can be seen as complementary to the fusion rules of belief functions. In what follows, we propose a method to combine BBAs to get the belief mass of TFE based on the construction of own domain and conflict domain.

Fusion based on partially uniform redistribution rule

Consider q ≥ 2 BBAs m 1 (•),. . . , m q (•) with q sets of focal elements U 1 ,. . . , U q , and a subset X of

U 1 × U 2 × • • • × U q .
The (conjunctive) fusion for the TFE X is calculated by m C 1,...,q (X) =

(ui 1 ,...,ui q )∈U C 1,...,q (X)

m 1 (u i1 ) • . . . • m q (u iq ).
(5) To reduce the complexity of the fusion combination rule, we propose to average the conflicting masses allocated to elements of U I(X) 1...q (∅) and to distribute it back to X, which is expressed as

m 1,...,q (X) = m C 1,...,q (X)+ 1 |U I(X) 1...q (∅)| (ui 1 ,...,ui q )∈U I(X) 1...q (∅) m 1 (u i ) × . . . × m q (u iq ). (6)
where |U

I(X) 1...q (∅)| is the cardinality of U I(X) 1...q (∅).
With such basic redistribution principle, after having computed the BBA m 1...q (X) of all possible TFE X, we don't get a normalized BBA. That is why a normalization step must be applied to get final fusion result, Example 2 (revisited): Let's consider X = T 1 = (x 1 , y 1 , z 1 ), then based on Table I one has

U C 1,2,3 (T 1 ) = {(x 1 , y 1 , z 1 ), (x 1 , y 1 ∪ y 2 , z 1 ), (x 1 , y 1 , z 1 ∪ z 2 ), (x 1 , y 1 ∪ y 2 , z 1 ∪ z 2 ), (x 1 ∪ x 2 , y 1 , z 1 ), (x 1 ∪ x 2 , y 1 ∪ y 2 , z 1 ), (x 1 ∪ x 2 , y 1 , z 1 ∪ z 2 ), (x 1 ∪ x 2 , y 1 ∪ y 2 , z 1 ∪ z 2 ), (x 1 ∪ x 2 ∪ x 3 , y 1 , z 1 ), (x 1 ∪ x 2 ∪ x 3 , y 1 ∪ y 2 , z 1 ), (x 1 ∪ x 2 ∪ x 3 , y 1 , z 1 ∪ z 2 ), (x 1 ∪ x 2 ∪ x 3 , y 1 ∪ y 2 , z 1 ∪ z 2 )}.
The (conjunctive) mass of belief of T 1 is given by

m C 1,...,3 (T 1 ) = (ui 1 ,ui 2 ,ui 3 )∈U C 1,2,3 (X) m 1 (u i1 )m 2 (u i2 )m 3 (u i3 ) = m 1 (x 1 )m 2 (y 1 )m 3 (z 1 ) + m 1 (x 1 )m 2 (y 1 )m 3 (z 1 ∪ z 2 ) + m 1 (x 1 ∪ x 2 )m 2 (y 1 )m 3 (z 1 ) + m 1 (x 1 ∪ x 2 )m 2 (y 1 )m 3 (z 1 ∪ z 2 ) + m 1 (x 1 ∪ x 2 ∪ x 3 )m 2 (y 1 )m 3 (z 1 ) + m 1 (x 1 ∪ x 2 ∪ x 3 )m 2 (y 1 )m 3 (z 1 ∪ z 2 ) + m 1 (x 1 )m 2 (y 1 ∪ y 2 )m 3 (z 1 ) + m 1 (x 1 )m 2 (y 1 ∪ y 2 )m 3 (z 1 ∪ z 2 ) + m 1 (x 1 ∪ x 2 )m 2 (y 1 ∪ y 2 )m 3 (z 1 ) + m 1 (x 1 ∪ x 2 )m 2 (y 1 ∪ y 2 )m 3 (z 1 ∪ z 2 ) + m 1 (x 1 ∪ x 2 ∪ x 3 )m 2 (y 1 ∪ y 2 )m 3 (z 1 ) + m 1 (x 1 ∪ x 2 ∪ x 3 )m 2 (y 1 ∪ y 2 )m 3 (z 1 )
= 0.56.

The incompatible domain U I(X)

1...q (∅) has 24 elements. They correspond altogether to a conflicting belief mass value given by the total mass value (equal to one) minus the mass of belief of the 12 q-tuples compatible with T 1 . That is Example 3 (revisited): We consider the following TFE X = T * = (x 1 , y 1 , ?) with a missing attribute (the third component is missing). Then, one has

U C 1,2,3 (T * ) = {(x 1 , y 1 , z 1 ), (x 1 , y 1 ∪ y 2 , z 1 ), (x 1 , y 1 , z 1 ∪ z 2 ), (x 1 , y 1 ∪ y 2 , z 1 ∪ z 2 ), (x 1 ∪ x 2 , y 1 , z 1 ), (x 1 ∪ x 2 , y 1 ∪ y 2 , z 1 ), (x 1 ∪ x 2 , y 1 , z 1 ∪ z 2 ), (x 1 ∪ x 2 , y 1 ∪ y 2 , z 1 ∪ z 2 ), (x 1 ∪ x 2 ∪ x 3 , y 1 , z 1 ), (x 1 ∪ x 2 ∪ x 3 , y 1 ∪ y 2 , z 1 ), (x 1 ∪ x 2 ∪ x 3 , y 1 , z 1 ∪ z 2 ), (x 1 ∪ x 2 ∪ x 3 , y 1 ∪ y 2 , z 1 ∪ z 2 ), (x 1 , y 1 , z 2 ), (x 1 ∪ x 2 , y 1 , z 2 ), (x 1 ∪ x 2 ∪ x 3 , y 1 , z 2 ), (x 1 , y 1 ∪ y 2 , z 2 ), (x 1 ∪ x 2 , y 1 ∪ y 2 , z 2 ), (x 1 ∪ x 2 ∪ x 3 , y 1 ∪ y 2 , z 2 )} Compared to the U C 1,2,3 (T 1 ), U C 1,2,3 (T * ) contains six extra elements. The (conjunctive) mass of belief of T * is given by m C 1,2,3 (T * ) = 0.56 + m 1 (x 1 )m 2 (y 1 )m 3 (z 2 ) + m 1 (x 1 ∪ x 2 )m 2 (y 1 )m 3 (z 2 ) + m 1 (x 1 ∪ x 2 ∪ x 3 )m 2 (y 1 )m 3 (z 2 ) + m 1 (x 1 )m 2 (y 1 ∪ y 2 )m 3 (z 2 ) + m 1 (x 1 ∪ x 2 )m 2 (y 1 ∪ y 2 )m 3 (z 2 ) + m 1 (x 1 ∪ x 2 ∪ x 3 )m 2 (y 1 ∪ y 2 )m 3 (z 2 ) = 0.56 + 0 = 0.56.
The conflicting mass of T * is given by the total mass value (equal to one) minus the mass of belief of all the 18 triplets compatible with T * , that is

IV. SUMMARY OF FUSION BASED ON RST

Figure 1 shows the whole diagram of our proposed combination method and the most important steps are the construction of GSs and calculation of belief mass of TFE. In the following example, we choose conflict allocation in new combination method to fuse sources with distinct FoDs.

V. APPLICATION EXAMPLE

Here, we will use the target identification example available in the literature [START_REF] Ristic | Target identification using belief functions and implication rules[END_REF] to illustrate our new method for combining SoE with distinct FoDs. This example is about an air surveillance system, which determines the targets' ID by combining the outputs of multiple sensors. Suppose that the FoD for the target allegiance is Θ A , and the FoD for the target classes

Example 1 :

 1 Consider a set U = {u 1 , . . . , u 10 } of focal elements of a normalized BBA m(•) whose values are m(u 1 ) = 0.0156 m(u 2 ) = 0.0874 m(u 3 ) = 0.0156; m(u 4 ) = 0.0874 m(u 5 ) = 0.0874 m(u 6 ) = 0.2354 m(u 7 ) = 0.2354 m(u 8 ) = 0.2354 m(u 9 ) = 0.0002; m(u 10 ) = 0.0002 If we consider the simple equivalence rule R defined by

1 |U

 1 .3 (∅) = 1.0 -0.56 = 0.44 Therefore the average conflicting mass value is. .3 (∅) = 1 24 • 0.44 ≈ 0.0183. Thus the total belief mass committed to T 1 (before normalization) will be m 1,...,3 (T 1 ) = m C 1,...,3 (T 1 ) + We apply the same principle of derivation for all targets of Θ T . After the normalization step, we finally get m 123 (T 1 ) ≈ 0.2274 m 123 (T 2 ) ≈ 0.2485 m 123 (T 3 ) ≈ 0.0558 m 123 (T 4 ) ≈ 0.2174 m 123 (T 5 ) ≈ 0.1998 m 123 (T 6 ) ≈ 0.0511.

  [START_REF] Janez | Theory of evidence and non-exhaustive frames of discernment: Plausibilities correction methods[END_REF] (∅) has 36 -18 = 18 elements, one getsm 1,2,3 (T * ) =m C 1,2,3 (T * )

s × b 23 {f × b 23 {

 2323 B 3 ),(Θ A , Θ ref B )} {(s, b 23 ),(Θ A , b 23 ),(f, b 22 ), (f, b 23 ),(f, B 2 ),(f, B 2 B 3 ), (s, b 23 ), (s, B 2 ),(s, B 2 B 3 ), (s, Θ ref B ),(Θ A , b 23),(Θ A , B 2 ), (Θ A , B 2 B 3 ),(Θ A , Θ ref B )} {(s, b 22 ),(Θ A , b 22 ),(f, b 22 ), (f, b 23 ),(f, B 2 ),(f, B 2 B 3 ), (f, Θ ref B )} 0.0984 s × B 2 {(s, B 2 ),(s, B 2 B 3 ),(s, Θ ref B ), (Θ A , B 2 ),(Θ A , B 2 B 3 ),(Θ A , Θ ref B )} {(f, b 22 ),(f, b 23 ),(f, B 2 ), (f, B 2 B 3 ),(f, Θ ref B )} 0.0504 s × B 2 B 3 {(s, B 2 B 3 ),(s, Θ ref B ), (Θ A , B 2 B 3 ),(Θ A , Θ ref B )} {(f, b 22 ),(f, b 23 ),(f, B 2 ), (f, B 2 B 3 ),(f, Θ ref B )} 0.0399 s × Θ ref B {(s, Θ ref B ),(Θ A , Θ ref B )} {(f, b 22 ),(f, b 23 ),(f, B 2 ), (f, B 2 B 3 ),(f, Θ ref B )} 0.0369 f × b 22 {(f, b 22 ), (f, B 2 ),(f, B 2 B 3 ), (f, Θ ref B ),(Θ A , b 22 ),(Θ A , B 2 ), (Θ A , B 2 B 3 ),(Θ A , Θ ref B )} {(f, b 23 ),(Θ A , b 23 ),(s, b 22 ), (s, b 23 ),(s, B 2 ),(s, B 2 B 3 ), (f, b 23 ), (f, B 2 ),(f, B 2 B 3 ), (f, Θ ref B ),(Θ A , b 23 ),(Θ A , B 2 ), (Θ A , B 2 B 3 ),(Θ A , Θ ref B )} {(f, b 22 ),(Θ A , b 22 ),(s, b 22 ), (s, b 23 ),(s, B 2 ),(s, B 2 B 3 ), (s, Θ ref B )} 0.1321f × B 2 {(f, B 2 ),(f, B 2 B 3 ),(f, Θ ref B ), (Θ A , B 2 ),(Θ A , B 2 B 3 ),(Θ A , Θ ref B )} {(s, b 22 ),(s, b 23 ),(s, B 2 ), (s, B 2 B 3 ),(s, Θ ref B )} 0.0406 f × B 2 B 3 {(f, B 2 B 3 ),(f, Θ ref B ), (Θ A , B 2 B 3 ),(Θ A , Θ ref B )} {(s, b 22 ),(s, b 23 ),(s, B 2 ), (s, B 2 B 3 ),(s, Θ ref B )} 0.0232 f × Θ ref B {(f, Θ ref B ),(Θ A , Θ ref B )} {(s, b 22 ),(s, b 23 ),(s, B 2 ), (s, B 2 B 3 ),(s, Θ ref B )} 0.0182 Θ A × b 22 {(Θ A , b 22 ),(Θ A , B 2 ), (Θ A , B 2 B 3 ),(Θ A , Θ ref B )} {(s, b 23 ),(f, b 23 )} 0.0887 Θ A × b 23 {(Θ A , b 23 ),(Θ A , B 2 ), (Θ A , B 2 B 3 ),(Θ A , Θ ref B )} {(s, b 22 ),(f, b 22 )} 0.1054 Θ A × B 2 {(Θ A , B 2 ),(Θ A , B 2 B 3 ),(Θ A , Θ ref B )} {∅} 0.0072 Θ A × B 2 B 3 {(Θ A , B 2 B 3 ),(Θ A , Θ

TABLE II :

 II Available Sensor Reports.

TABLE V :

 V Fused BoE generated from four sensor reports. A , b 22 ),(Θ A , B 2 ), (Θ A , B 2

	Target X	Domain	Belief Mass
		Compatible domain	Incompatible Domain
		{(s, b 22 ), (s, B 2 ),(s, B 2 B 3 ),
	s × b 22	(s, Θ ref B ),(Θ

Here, we use the symbol to mean equal by definition.

If Shafer's model doesn't fit with the problem.

In this paper, U C 1,...,q (X) and U I(X) 1,...,q (∅) are regarded as a variety of GSs in RST.

Modification of definition 2 is similar.

This equivalence relation is also called tolerance relation in RST.

VII. ACKNOWLEDGMENT This work was supported in part by the National Natural Science Foundation of China under Grant 61573097 and 91748106, in part by Key Laboratory of Integrated Automation of Process Industry (PAL-N201704), the Fundamental Research Funds for the Central Universities (3208008401), the Qing Lan Project and Six Major Top-talent Plan, and in part by the Priority Academic Program Development of Jiangsu Higher Education Institutions. The authors thank the reviewers and editors for giving valuable comments, which are very helpful for improving this manuscript.

Report ID

Corresponding evidence BoE

Report (1) m

(1) Θ B (B 2 ) = 0.7 m Hence, the FoD II, and fusion center combines all the sensor reports and generates a BBA for decisionmaking. Suppose that the pairs {Report (1) , Report (2) } and {Report (3) , Report (4) } are fused separately. The resulting BBAs are vacuously extended and fused by the new combination method in this paper so as to obtain the fused mass m Θ (•) as

Based on the original data given by four sensor reports in Table II, we can obtain four initial GSs:

Based on definitions of self-domain, conflict-domain, we can obtain some new and interesting results which are shown in Table III, IV, V: In Table III, since focal elements in U 1 and U 2 exist inclusion relations shown in Figure 2, for example,

thus, the intersections between each element in U 1 and U 2 are not ∅, and according to definitions 3, the conflict domain is empty set, and the elements of its self domain are shown in Table III. Furthermore, according to definition 4, we can calculate the belief mass of B 2 : m

Similarly, in Table IV, the conflict between U 3 and U 4 is caused by s and f , thus, conflict domain includes these two elements. Then, based on Eq.( 5) and Eq.( 6), we can get the corresponding masses of s, f and A, which are given in Table IV

Then, in Table V, we can get the self domains and conflict domains of all TFEs. Finally, we can get the belief masses of all elements in Θ A × Θ B . Based on the maximum of belief mass principle, they appear to favor a target ID of f × b 22 (i.e., a friendly b 22 ), which is consistent with the conclusions in [START_REF] Ristic | Target identification using belief functions and implication rules[END_REF] and [START_REF] Wickramarathne | An analytical framework for soft and hard data fusion: A Dempster-Shafer belief theoretic approach[END_REF]. Remark 3: This example shows two advantages of our new combination method in this paper: [START_REF] Shafer | A mathematical theory of evidence[END_REF] we do not need to expand FoDs of four sensor reports to ensure that they are under the unique FoD. For example, if we use the classical TABLE III: Domains and belief masses of TFEs in {Report (1) , Report (2) }.

TABLE IV: Domains and belief masses of TFEs in {Report (3) , Report (4) }.

A )} {(s, f )} 0.24

A , Θ

A )} {∅} 0.22 ), Θ B }, and then we update the BBAs by assigning zero masses to the elements that are not in their specific FoDs. This expanded way will obviously increase the computational complexity; (2) On the other hand, the new combination method can effectively get the belief mass of the special single target (for example in this case, hostile Airbus-320) when there are no need to pay more attentions to other potential targets. However, all belief masses of targets will be calculated in a single calculation step in traditional combination ways.

VI. CONCLUSION

In this paper we did propose a new combination method to solve the fusion problem of sources with distinct FoDs. Specifically, the new fusion process is divided into two steps: 1) according to GSs provided by original sources, we construct the compatible domain and the incompatible domain of TFE based on several rules; and 2) we calculate the belief mass of TFE with very basic belief allocation strategy. Experiments show that the approach presented in this paper can provide a solution to the fusion of sources with distinct FoDs compatible with more sophisticate methods. In future works, we will continue to explore more rules to get better relevant sets of TFE and to improve the conflict redistribution strategies in order to combine more efficiently (hopefully) different sources defined on distinct FoDs.


domain: spi phys math info

text:
Pierre Godet
email: pierre.godet@onera.fr

Frédéric Champagnat

Guy Le Besnerais

Aurélien Plyer

Learning fluid trajectory models for time-resolved PIV

Introduction

Compared to classical particle image velocimetry, time-resolved PIV (TR-PIV) allows to characterize the time evolution of unsteady phenomena instead of their instantaneous and mean properties. However, as sensors and laser sources work at higher rates, it leads to a reduced signal-to-noise ratio. To deal with this issue, one can exploit the temporal coherence in a multiframe estimation approach. Furthermore, multiframe motion estimation allows computing short trajectories, instead of two-frame displacements. Thus, other quantities than instantaneous velocity can be evaluated such as fluid acceleration.

FTC of [START_REF] Lynch | A high-order time-accurate interrogation method for time-resolved PIV[END_REF], FTEE of [START_REF] Jeon | Fluid trajectory evaluation based on an ensemble-averaged crosscorrelation in time-resolved piv[END_REF], and the LKFT approach, Yegavian et al. (2016), perform multiframe estimation for PIV assuming a polynomial time dependency. Thanks to this polynomial model, multiframe estimation gains robustness against the noise. It can be seen in the figure 1 which compares the errors on the estimated displacement fields with LKFT or with the two-frame algorithm FOLKI, [START_REF] Champagnat | Fast and accurate PIV computation using highly parallel iterative correlation maximization[END_REF]. Nevertheless, LKFT presents localized errors which, as we will see later, appear in regions where the polynomial model is not perfectly adapted to the movement. Here we propose to investigate the use of more precise motion models derived from the data.

For instance, in the context of computer vision, [START_REF] Garg | Dense multi-frame optic flow for non-rigid objects using subspace constraints[END_REF] proposed to use principal component analysis (PCA, known as POD for Proper Orthogonal Decomposition in the fluid mechanics' community) to derive trajectory models adapted to a given video sequence. More generally, the current trend of data-based approaches motivates the development of generic multiframe motion estimation methods able to exploit a learned basis of trajectories. In this purpose, we elaborate on ideas of [START_REF] Champagnat | Fast and accurate PIV computation using highly parallel iterative correlation maximization[END_REF]Yegavian et al. (2016) and derive an efficient a lgorithm t o e stimate t he c oefficients of th e de composition of fluid trajectories on an arbitrary basis by window-based cross-correlation analysis.

We present a first application of this algorithm on a PCA-learned basis of trajectory m odels. Compared to [START_REF] Garg | Dense multi-frame optic flow for non-rigid objects using subspace constraints[END_REF], our estimation is based on a Lucas-Kanade algorithm instead of a variational approach. The former is more adapted to the noise level encountered in PIV application, and also leads to faster implementation.

Previous work

FOLKI is a Lucas-Kanade algorithm first b een p ublished f or c omputer v ision a pplications, L e [START_REF] Besnerais | Dense optical flow by iterative local window registration[END_REF], and applied later to the context of PIV, [START_REF] Champagnat | Fast and accurate PIV computation using highly parallel iterative correlation maximization[END_REF]. The Lucas-Kanade Fluid Trajectories (LKFT) approach, Yegavian et al. (2016), extended some ideas of FOLKI to the multiframe motion estimation for time-resolved PIV. LKFT assumes a polynomial formulation of the motion and directly estimates the coefficients by cross-correlation analysis on an interrogation w indow. Yegavian et al. (2016) showed in particular that temporal modelling improves the robustness to noise of LKFT. Such gain had been observed in previous references. FTEE [START_REF] Jeon | Fluid trajectory evaluation based on an ensemble-averaged crosscorrelation in time-resolved piv[END_REF] also assumes a polynomial time 2016), both with an interrogation window of 33 × 33 pixels. Estimation realized on the simulated sequence presented below in section 4.1, degraded by an additive Gaussian noise with a level of 4 %. LKFT was runned with a polynomial degree of 3 on a sequence of 8 images. dependency but uses a different method to estimate the coefficients. Fluid Trajectory Correlation (FTC) of [START_REF] Lynch | A high-order time-accurate interrogation method for time-resolved PIV[END_REF] tracks independent fluid elements over multiple image-pairs and then applies a polynomial fit to each estimated trajectory.

In the computer vision community, it has been shown by Tomasi and Kanade (1992) that the image motion of a rigid scene could be described in a very low-dimensional linear subspace. [START_REF] Irani | Multi-frame correspondence estimation using subspace constraints[END_REF] used this idea to compute a motion basis to constrain the multiframe optical flow estimation for a rigid motion, and [START_REF] Garg | Dense multi-frame optic flow for non-rigid objects using subspace constraints[END_REF] extended this principle to the non-rigid case using PCA. PCA has also been used to learn a motion basis in [START_REF] Black | Learning parameterized models of image motion[END_REF] In fluid mechanics, PCA (or POD) is used under to estimate spatial motion models. POD on spatiotemporal patches is studied in [START_REF] Stapf | A learning-based approach for highly accurate measurements of turbulent fluid flows[END_REF]. In contrast, we focus here on learning temporal motion models.

Approach 3.1 A generic multiframe motion algorithm

In [START_REF] Champagnat | Fast and accurate PIV computation using highly parallel iterative correlation maximization[END_REF], the displacement field between two frames is estimated, at each pixel location k, by minimizing the sum over an interrogation window W of squared differences between image intensities. In a multiframe context, Yegavian et al. ( 2016) chose a reference frame I re f among the N considered consecutive images, and defined a multiframe criterion as the sum of the two-frame criteria for every image pair involving the reference image:

N ∑ n=1 ∑ k W (k -k ) (I re f (k ) -I n (k + u(k,t n ))) 2 (1)
LKFT reduces the number of indeterminates by approximating the time dependency of motion as a polynomial model. Here we consider a generalization of this approach where the motion is projected on some basis of temporal models. To be more specific: each component (horizontal u and vertical v) of the displacement is written as a linear combination of temporal models, e i and h i , which are vectors of R N , N being the length of the sequence:

u(k,t n ) = ∑ d i=1 θ (u) k,i e i [n] v(k,t n ) = ∑ d i=1 θ (v) k,i h i [n]
(2) the temporal models are polynomial (ie. LKFT (Yegavian et al., 2016), or are learned on motion trajectory sample in our new method. As a result, for each pixel k of the reference image, instead of estimating the displacement at every time-step, a small set of 2d coefficients θ k,i describes the whole trajectory passing over the pixel's location at the reference time instant.

e i [n] = h i [n] = t i n ) in
The 2d coefficients are estimated by iterative minimization of the multiframe cost function (1). At each iteration, (2) is approximated by a linear least-squares criterion and its optimization is reduced to the resolution of a 2d × 2d linear system for each pixel. In our approach, the first order expansion of (1) follows the inverse form of FOLKI described in [START_REF] Champagnat | Fast and accurate PIV computation using highly parallel iterative correlation maximization[END_REF]. As shown in the appendix below (in the case e i = h i ), we take advantage of a space-time separability property to propose a very efficient two-step resolution: a d × d-matrix inversion made once at the beginning (for all pixels, time-steps and iterations) and a 2 × 2-matrix inversion at each pixel and each iteration. Note that, although this has not been demonstrated in the appendix because of the lengthy mathematical developments, this new algorithm can handle two different bases for the horizontal and vertical components of the motion field.

Our formulation includes the LKFT paradigm. Indeed, using a polynomial model leads to results identical to the initial implementation of LKFT as described in Yegavian et al. ( 2016), but with a drastic reduction of the computational time and memory usage. However, the genericity of this algorithm, as it can work with any trajectory basis, makes it possible to use temporal models learned from the data.

Learning a motion basis from the data

The proposed algorithm offers the possibility to regularize motion estimation with a data-driven model: a new specialized basis of temporal models can be used for every new image sequence.

To build a basis adapted to the data, we use PCA on a set of trajectories representative of the studied sequence. Representative trajectories could be estimated from the image sequence in a preprocessing step by particle tracking velocimetry methods, or be extracted from a simulation corresponding to the experimental scenario at hand. We gather such trajectories in a large data matrix and apply a singular-value decomposition. The selected motion basis is made of the vectors related to the highest singular values. The rank of the basis can be set manually or chosen automatically by thresholding the cumulated energy of singular values.

Two different learning processes are proposed here. In the first one, a PCA is applied to temporal samples of the horizontal and vertical displacements taken indifferently together. This approach leads to a single basis for both components -e i = h i in Eq. (2) -and is called "isotropic PCA" in the following. In the second one, horizontal and vertical samples are separated, and two PCAs are done independently on each set. Two different bases are obtained, each specialized for one component. This second process is called "asymmetric PCA" in the following.

Experiments

Temporal modelling allows reducing the number of estimated coefficients with respect to the classical PIV processing of each pair of consecutive images. This not only reduces computational time but also increases robustness to noise thanks to temporal redundancy. However, the estimation quality depends on the choice of the model: a trajectory stepping out of the model cannot be well estimated. The aim is to provide new trade-offs between noise and model errors.

Through different simulated experiments, we have compared the estimations given by our new generic algorithm with different basis choices:

• a polynomial basis, leading to the same result than LKFT • a basis shared by both horizontal and vertical components estimated by the isotropic PCA process.

• two different bases, one for each component, estimated by the asymmetric PCA process.

To focus on the ability of these models to adapt to a complex motion sequence, we compute the PCA on the best data-suited trajectories we can hope for by using a subset of the simulated trajectories. The estimated multiframe motion fields are 3D-tensors which provide for each spatial location, given by the pixel grid of the reference image, the trajectory on the whole given sequence. In all the following experiments the reference image is at the center of the sequence (for instance, for a sequence of 5 or 6 images, the reference frame is the third one). On the one hand, the result can be evaluated in the whole spatial field at one specific time-step: for the first time-step following the reference image this corresponds to an instantaneous velocity, it can be compared with other PIV results. In Sec. 4.2 we study instantaneous spatial fields of multiframe methods, a comparison with the two-frame FOLKI algorithm of [START_REF] Champagnat | Fast and accurate PIV computation using highly parallel iterative correlation maximization[END_REF] using the same interrogation window size is displayed in Fig. 1. On the other hand, by selecting a particular spatial location, the whole estimated trajectory can be studied, as proposed in Sec. 4.4.

Testing sequence

In the following experiments, we use the velocity field in the wake of a flapping foil at Reynolds number 250, obtained by direct numerical simulation, [START_REF] Jallas | Linear and nonlinear perturbation analysis of the symmetry breaking in time-periodic propulsive wakes[END_REF]. The streamwise and cross-stream components of the instantaneous displacement at the central time-step are shown on figure 2. Using these simulated velocity fields, we compute trajectories and generate a particle image sequence on which the PIV methods are tested. The particles are Gaussian-shaped with a standard deviation of 0.4 pixels. To evaluate robustness to noise, an additive Gaussian noise is added to the sequence for some of the following experiments. To represent a time-resolved PIV context, this sequence is temporally well sampled, leading to interframe displacements of the order of one pixel. Working on such simulated data makes it possible to derive quantitative comparison metrics with regard to the ground truth motion. We use the "endpoint error" defined as the norm of the difference between the estimated vector and the ground truth vector for each pixel location.

Instantaneous motion field: spatial analysis

In this section, we study the errors on the instantaneous displacement field at the time-step following the reference frame. We first compare the potential of the PCA decomposition to the polynomial model, by computing the projection error of the simulated known motion field of 2 on both bases. The projection error is the endpoint error between the ground truth field and its projection on the considered basis (polynomial or obtained by PCA). By doing so, it can be determined to what extent the different models are adapted to the studied motion. All trajectories are projected on the PCA basis on the one hand, and on the polynomial basis on the other. We consider 8 time-steps, an order 3 polynomial model and a PCA basis made of 3 vectors. According to the projection error displayed in figure 3 (upper line), the PCA leads to a better approximation in regions with high motion gradients.

We also present results of multiframe PIV processing of the generated sequence of images introduced before. Error fields on the cross-stream displacement are shown on the lower line of Figure 3. LKFT with an order 3 polynomial is compared to the proposed PCA-based multiframe estimation with 3 vectors; both methods use a 33 × 33 interrogation windows and on the same sequence of 8 images. These parameter choices are justified in the parametric study in section 4.3. We use a noise level of 4 % of the maximal intensity. The localized model errors seen on the projected field (upper line) match the ones found on the estimated field, and are also slightly higher with the polynomial basis. Thus, the basis learned by PCA appears better suited to the studied motion, and its use in the whole dense motion estimation process leads to better results than the polynomial basis. We did not present the results for the "asymmetric PCA", because they are very similar to those obtained with the "isotropic PCA".

The multiframe paradigm can be compared to classical two-frame PIV on figure 1 and4. Figure 1 shows that unlike the two-frame estimation which is globally noisy, the multiframe one is well regularized except in some areas where the error is higher because the local motion steps out of the chosen model. Figure 4 enlights that when the noise level increases, the error of the two-frame method increases faster than the error of the multiframe algorithms. Multiframe estimation is more robust to noise. However, for low noise levels, the two-frame estimation error is lower than the multiframe one in the areas of high bias.

To estimate an instantaneous motion field, multiframe approaches appear more robust to noise but also exhibit localized errors in the high-gradient regions, where the model is not fully suited to the motion. In these regions, the learned PCA basis allows a reduction of the bias compared to the polynomial approach.

Parametric study

For LKFT as for our new algorithm, some parameters can be tuned :

• The radius of the interrogation window (a radius of 16 means an interrogation window of size 33×33).

• The number N of consecutive images considered for the estimation.

• The number d of vectors (temporal models) in the basis. For LKFT it corresponds to the polynomial degree.

For the problem to be well-defined, N has to be greater than d. By plotting the averaged endpoint error in the instantaneous field at the center of the sequence, as a function of these parameters, we determine their influence on the estimation quality.

The influence of the size of the interrogation window is shown on figure 5. As the window gets bigger, the error of every tested method decreases sharply, reaches an optimum and then increases slowly. If the window is too small the problem is not well-conditionned, if it is too big the estimation fails to capture fine details. The optimum of the multiframe methods varies with the dimension d of the basis: a small number of temporal models means a stronger temporal regularisation which can compensate a lack of spatial regularisation due to a smaller window. The two-frame method error is always higher, and the gap with the multiframe approaches is larger for small windows. In our case, a 33 × 33 window is a good choice for both d = 2 and d = 3.

Figure 6 shows the influence of the length of the sequence. The two-frame PIV error level is given for comparison but, of course, does not depend on the sequence length. By increasing the number of frames, one can hope a growing robustness to noise and a lower estimation error. This is true if the number d of temporal models in the basis is high enough to capture the complexity of the motion. As this complexity generally grows when raising the sequence length, using many frames implies to increase the rank of the basis to limit the bias error. This can be seen on figure 6: with two temporal models, the multiframe methods are better than the two-frame's only for short sequences (4 or 5 frames), and then the error radically increases; longer sequences can be used with 3 models, until 8 or 10 frames before the error goes up.

Finally, note that in almost all these experiments, PCA models lead to lower or equivalent errors than the polynomial model, with a slight advantage for the asymmetric PCA model.

Trajectories: temporal analysis

Multiframe paradigm also gives the possibility to consider the trajectories over the whole sequence, some are presented on figure 7. For this evaluation, we consider a longer sequence, with 25 frames, 3 temporal models and no added noise on the image sequence. We focus on two different areas of the field (underlined by a white square in Fig. 3): the first one exhibit high motion gradients, the second one has smoother Second row: error on the displacement fields estimated by multiframe estimations with LKFT's polynomial approach (left) and the proposed algorithm using an "isotropic PCA" basis (right), both with an interrogation window of 33 × 33 pixels. A sequence of 8 images was used, LKFT was tuned with a polynomial degree of 3, and the PCA basis was composed of 3 vectors. A Gaussian noise with a level of 4 % was added to the images. motions. We compare the estimated trajectories by our algorithm for each basis choice to the ground truth trajectories. In the high motion gradient case (figure 7, first row), compared to the polynomial model, the "isotropic PCA" significantly improves the trajectory estimation. This is particularly noticeable on the left trajectories of the high-gradient area. The "asymmetric PCA" provides some further minor improvements. In the smoother case (figure 7, second row), the polynomial model is globally better. The PCA approaches make some errors, which are a bit higher with the "isotropic PCA".

The PCA bases represent well the non-polynomial motions; however, by construction, such bases focus on the statistically dominant trajectories. As a consequence, this method is unlikely to give a reasonable estimate of a motion which only appears in a small part of the field, and which is too different from the rest of the field.

Conclusion

We have proposed a novel and efficient algorithm for multiframe PIV using an arbitrary temporal basis. The genericity of the proposed algorithm opens the way to better noise/model error trade-off compared to polynomial approaches. As a first application, we proposed to use a basis learned on the studied sequence by PCA. We have evaluated this approach on a simulated unsteady flow, and obtained better results than with the polynomial model.

This work can be considered as a proof of concept on learned temporal models for multiframe motion estimation. The next step will be to apply it on real time-resolved PIV data. This may be challenging because of the following issues. First of all, image degradation encountered in real datasets are much more complicated than the Gaussian noise used in our simulation study. Moreover, a strategy for learning basis in a real experiment has to be designed. It could rely on a preprocessing of the data, for instance exploiting trajectories obtained by particle tracking, or use trajectories sampled from a representative numerical simulation.

Finally, in the proposed method, the basis used is the same for every spatial location. Further improvement could be made by using different bases for different regions of the field. Methods combining motion estimation and segmentation could be explored here. A Algorithm details: the case of a single basis for both coordinates

To minimize the criterion of equation 1 we adopt an inverse Lucas-Kanade iteration, [START_REF] Champagnat | Fast and accurate PIV computation using highly parallel iterative correlation maximization[END_REF], for better computational performances. The aim is to find, at every pixel location k in the reference frame, the trajectory (u(k,t n )) n∈[1.

.N] minimizing the linearized criterion:

∑ n=1..N ∑ k ∈W (k) ε(k , n) -(∇I re f ) T (k ) • u(k,t n ) 2 with ε(k , n) = I re f (k ) -Ĩn (k + u (0) (k ,t n )) + (∇I re f ) T (k ) • u (0) (k ,t n )
The number of indeterminates is reduced thanks to the decomposition:

u(k,t n ) = u(k,t n ) v(k,t n ) = T n θ k (3)
where θ k stands for the 2d coefficients:

θ k = [θ (u)
k,1 , . . . θ can be written as:

u(k,t n ) = d ∑ i=1 θ (u) k,i e i [t n ] v(k,t n ) = d ∑ i=1 θ (v) k,i e i [t n ]
The criterion to minimize becomes:

J k (θ) = ∑ n=1..N ∑ k ∈W (k) ε(k , n) -(∇I re f ) T (k )T n θ 2 (4) = θ T A k θ -2b T k θ + constant (5)
its solution is:

θ k = A -1 k b k , (6)
with A k a 2d × 2d-matrix given by:

We then rewrite the matrix A k of Eq. ( 7) by defining:

α k β k β k γ k ∆ = ∑ k ∇I re f (k )(∇I re f ) T (k )
and the d × d-matrix T:

T ∆ = ∑ n E T n E n
so that we get:

A k = α k T β k T β k T γ k T ,
From this, we have the following expression of the inverse:

A -1 k = 1 α k γ k -β 2 k γ k T -1 -β k T -1 -β k T -1 α k T -1 , (9)
The main operation to compute the inverse of A k is then the computation of the inverse of T, which can be done only once, for all frames and iterations. At each iteration, it will only remain the computations of α k , β k and γ k .

The time and space separability and the use of a shared temporal basis (making T n diagonal) reduce computational costs. In the case of two different bases, the separability still exists but is less straightforward and needs additional computations such as Cholesky factorizations and eigenvalues decompositions.

Figure 1 :

 1 Figure 1: Error on the cross-stream displacement fields estimated by two-frame PIV, Champagnat et al. (2011) and multiframe estimation with LKFT's polynomial approach, Yegavian et al. (2016), both with an interrogation window of 33 × 33 pixels. Estimation realized on the simulated sequence presented below in section 4.1, degraded by an additive Gaussian noise with a level of 4 %. LKFT was runned with a polynomial degree of 3 on a sequence of 8 images.

  ; Wulff and Black (2015) but their basis vectors represent spatial dependency whereas Garg et al. (2010); Yegavian et al. (2016) model time dependency.

Figure 2 :

 2 Figure 2: Streamwise and cross-stream components of the ground truth instantaneous displacement at the central time-step.

Figure 3 :

 3 Figure 3: Error (in pixels) on the cross-stream component of the instantaneous displacement at the central time-step. First row: projection error on the polynomial (left) and PCA (right) bases.Second row: error on the displacement fields estimated by multiframe estimations with LKFT's polynomial approach (left) and the proposed algorithm using an "isotropic PCA" basis (right), both with an interrogation window of 33 × 33 pixels. A sequence of 8 images was used, LKFT was tuned with a polynomial degree of 3, and the PCA basis was composed of 3 vectors. A Gaussian noise with a level of 4 % was added to the images.

Figure 4 :

 4 Figure 4: Endpoint error [px] on the instantaneous displacement, spatial average and maximum (for the maximum we ignore the 50 pixels nearest to each boundaries), as a function of the additive Gaussian noise level. Sequence of 5 images, basis of 2 vectors for multiframe methods, interrogation window of 33 × 33 pixels. A noise level of ten corresponds to 4% of the maximum intensity.

Figure 5 :

 5 Figure 5: Averaged endpoint error [px] on the instantaneous displacement, as a function of the size of the interrogation window. Sequence of 5 images, basis of d vectors for multiframe methods. A Gaussian noise with a level of 4 % was added to the images.

Figure 6 :

 6 Figure 6: Averaged endpoint error [px] on the instantaneous displacement, as a function of the number of images in the sequence, with an interrogation window of 33 × 33 pixels. A Gaussian noise with a level of 4 % was added to the images. The two-frame PIV error level does not depend on the sequence length but is given for comparaison.

Figure 7 :

 7 Figure 7: Estimated (yellow) and ground truth (blue) trajectories for the different methods, with 3 temporal models and a 25-image sequence. Two areas are represented : where motion gradients are strong (first row) or weak (second row). The corresponding areas are squared in white in other preceding figures.

  1 [t n ] • • • e d [t n ]] Twith e 1 , • • • , e d any basis of temporal models. To be more specific, the two components ofu(k,t n ) = u(k,t n ) v(k,t n )

T

  ∇I re f (k )(∇I re f ) T (k )T n ∇I re f (k )(∇I re f ) T (k ) T n (7)This comes from the separability in space and time, equation 3. The 2d-vector b k is given by n ε(k , n)T ∇I re f (k )

Acknowledgements

The authors are grateful to Olivier Marquet and Benjamin Leclaire (ONERA/DAAA) who provided the test data used in this work. Benjamin Leclaire also brought a precious help through several discussions. This


Deviance of Poisson regression

Enwu Liu, Mary MacKillop Institute for Health Reaserch, ACU October 2019

In Poisson regression there are two Deviances. The Null Deviance 1 shows how well the response variable is predicted by a model that includes only the intercept (grand mean).

And the Residual Deviance 2 is -2 times the difference between the log likelihood evaluated at the maximum likelihood estimate (MLE) and the log-likelihood for a "saturated model" (a theoretical model with a separate parameter for each observation and thus a perfect fit). Now let us write down those likelihood functions. Suppose Y has a Poisson distribution whose mean depends on vector x, for simplicity, we will suppose x only has one predictor variable. We write

E(Y |x) = λ(x)
For Poisson regression we can choose a log or an identity link function, we choose a log link here.

Log(λ(x)) = β 0 + β 1 x β 0 is the intercept.
The Likelihood function with the parameter β 0 and β 1 is

L(β 0 , β 1 ; y i ) = n i=1 e -λ(xi) [λ(x i )] yi y i ! = n i=1 e e -(β 0 +β 1 x i ) e (β0+β1xi) yi y i !
The log likelihood function is:

l(β 0 , β 1 ; y i ) = - n i=1 e (β0+β1xi) + n i=1 y i (β 0 + β 1 x i ) - n i=1 log(y i !) (1)
When we calculate null deviance we will plug in β 0 into (1). β 0 will be calculated by a intercept only regression, β 1 will be set to zero. We write

l(β 0 ; y i ) = - n i=1 e β0 + n i=1 y i β 0 - n i-1 log(y i !) (2)

Figure 1 :

 1 Figure 1: Poisson regression results

Next we need to calculate the log likelihood for "saturated model" (a theoretical model with a separate parameter for each observation), therefore, we have µ 1 , µ 2 , ..., µ n parameters here.

(Note, in (1), we only have two parameters, i.e. as long as the subject have same value for the predictor variables we think they are the same).

The log likelihood function for "saturated model" is

Then it can be written as:

(Note, y i ≥ 0, when y i = 0, y i logµ i = 0 and log(y i !) = 0, this will be useful later, not now)

Now put μi into (3), since when y i = 0 we can see that μi will be zero (this is the reason why we need (3). Now for the likelihood function (3)of the "saturated model" we can only care y i > 0, we write

From ( 4) you can see why we need (3) since logy i will be undefied when y i = 0

Now the let us calcualte the deviances.

The

The Null Deviance=

Ok, next let us calculate the two Deviances by R then by "hand".

x< -c(2,15,19,14,16,15,9,17,10,23,14,14,9,5,17,16,13,6,16,19,24,9,12,7,9,7,15,21,20,20) y<-c(0,6,4,1,5,2,2,10,3,10,2,6,5,2,2,7,6,2,5,5,6,2,5,1,3,3,3,4,6,9) p_data<-data.frame(y,x) p_glm<-glm(y~x, family=poisson, data=p_data) summary(p_glm)