A Laurent
email: laurent@lirmm.fr

P Poncelet
email: pascal.poncelet@ema.fr

M Teisseire
email: teisseire@lirmm.fr

Fuzzy Data Mining for the Semantic Web: Building XML Mediator Schemas

Keywords: Fuzzy Logic, Semi-Structured Data, Mediator Schemas, XML, Semantic Structures, Frequent Patterns

As highlighted by the World Wide Web Consortium, XML has been proposed to deal with huge volumes of electronic documents and is playing an increasing important role in the exchange of a wide variety of data on the Web. However, when dealing with such large and heterogeneous data sources, it is necessary to have an idea on the way these data sources are structured. This information is indeed essential in order to build mediator schemas. These mediator schemas are required to query data in a uniform way. Moreover, this information is interesting since it provides users with a semantic structure of the data they can query. Recently schema mining approaches have been proposed to extract in an efficient way the commonly occurring schemas from a collection. Nevertheless, according to the semantic point of view, such approaches suffer from different drawbacks. In this work, we propose thus a fuzzy approach, showing why and how fuzziness is useful in order to extract frequent approximate schemas.

Introduction

Large amounts of data are available in XML format, and even if large volumes of legacy data are still marked up in HTML, efficient approaches have been proposed to transform HTML documents into XML documents [START_REF] Chung | Reverse engineering for Web data: From visual to semantic structure[END_REF]. These XML documents, stored in many sources distributed over the Web, contain useful information. For instance, many scientists have to manage and share large amounts of data (e.g. biologists have to deal with huge amounts of complex data). In this framework, XML is of great interest in order to represent this complicated knowledge.

Recently, sharing data has become very popular. For instance, Peer-based Data Management [START_REF] Tatarinov | The Piazza peer data management project[END_REF][START_REF] Halevy | Piazza: data management infrastructure for semantic web applications[END_REF] has become a great challenge for scientific data sharing, as Peer-to-Peer (P2P) systems are becoming more and more popular. Nevertheless the information stored in heterogeneous collections of XML documents has to be integrated in order to be processed or to be queried in a uniform way. In all fields like biology and medicine, there will never be a centralized database, which would be impossible to build. Scientists need thus to be provided with tools in order to reach the data they require from multiple databases. Database concepts can help to solve these problems. One result of the research led in data integration and in database federations is mediatorbased information systems [START_REF] Florescu | Database techniques for the world wide web: A survey[END_REF]. For instance, Xyleme is a huge warehouse integrating XML data on the Web [START_REF] Xyleme | A dynamic warehouse for XML data of the web[END_REF]. This integration can be performed by providing mediator schemas. A mediator schema can be considered as a shared structure through which queries can be defined [START_REF] Halevy | Schema mediation in peer data management systems[END_REF], as shown on Figure 1. However, no complete automatic tool is available to extract semantic knowledge from these large amounts of distributed and heterogeneous data. It is indeed still hardly possible to automatically build mediator schemas, which are thus still hand-crafted. In the past few years, data mining has been extensively studied in the framework of the semantic web. In this paper, we consider the problem of mining XML mediator schemas from a database perspective. By considering that a collection of XML documents could be modelled by labelled ordered trees [START_REF] Abiteboul | Data on the Web[END_REF][START_REF] Asai | Online algorithms for mining semi structured data stream[END_REF] (e.g. Figure 2 illustrates such a transformation [START_REF] Leung | On the use of hierarchical information in sequential mining-based XML document similarity computation[END_REF]), we focus on schema mining which provides tools to mine frequent sub-schemas from large databases. These frequent sub-schemas are then merged in order to build mediator schemas. We show that schema mining approaches have some drawbacks and that a fuzzy approach is very useful to mine approximate schemas. Moreover, frequent mediator schemas mined from our approach provide an interesting source of information about the data available on the web since these schemas can be seen as semantic structures describing the information. Our approach is based on the definition of the tree inclusion. We define several ways to introduce fuzzy logic in this problem. The main idea is to propose a definition of soft inclusion, meaning that a tree is no more included or not in another one, but gradually included within it. A degree of inclusion is defined, depending on the way the fuzzy inclusion is considered. Finally, we introduce fuzzy frequent patterns which aim at representing the strength of the links from a frequent tree.

The sequel of the paper is organized as follows. Section 2 goes deeper into presenting the problems of schema mining. In Section 3 we propose our motivations for introducing fuzzy logic in order to provide approximate mediator schemas. Section 4 addresses the Fuzzy Tree Inclusion while in Section 5 we propose to extend the knowledge on the mined semantic structures by providing fuzzy links within the frequent trees. Section 6 concludes the paper.

Schema Mining Principles

This section is devoted to stating the schema mining principles. As proposed approaches are mainly extensions of level-wise approach for candidate and frequent generation, we first propose an overview of the association rules problem. Second we consider the schema mining problem. We broadly summarize the problem statement and give an overview of the core of this process.

A Brief Overview of the Association Rules Problem

The problem of mining association rules, firstly introduced in 1993 in [START_REF] Agrawal | Mining association rules between sets of items in large database[END_REF], has received a great deal of attention. In brief, the problem is the following.

Let I = {x 1 , . . . , x n } be a set of distinct literals called items. A set X ⊆ I with k = |X| is called a k-itemset or simply an itemset. Let a database D be a multi-set of subsets of I. Each T ∈ D is called a transaction. We say that a transaction T ∈ D supports an itemset X ⊆ I if X ⊆ T holds. An association rule is an expression X ⇒ Y, where X, Y are itemsets and X∩Y = ∅ holds. The fraction of transactions T supporting an itemset X with respect to database

D is called the support of X, supp(X) = |{t ∈ D|X ⊆ T }|/|D|. The support of a rule X ⇒ Y is defined as supp(X ⇒ Y ) = supp(X ∪Y ). The confidence of the rule is defined as conf (X ⇒ Y ) = supp(X ∪ Y )/supp(X).
As the number of rules grows exponentially with |I|, the rule sets are typically restricted by minimal thresholds for the measures of support and confidence, minSupp and minConf respectively.

According to this restriction, the problem can thus be decomposed into two subproblems:

• Find all itemsets having support greater than or equal to minSupp where the support for an itemset is defined as the number of transactions that support the itemset. Such itemsets are called frequent itemsets. • Use the frequent itemsets to generate the set of rules according to minConf.

As solving the second subproblem is straightforward we will now have a look on how to obtain efficiently itemsets. In [START_REF] Agrawal | Fast algorithms for mining association rules[END_REF], the Apriori algorithm was defined to mine such associations. It makes multiple scans over the data (level-wise approach). In the first scan, the support of individual items is counted in order to determine which items are frequent. These frequent items are thus called frequent 1-itemsets. From them, new 2-candidates considered as potentially frequent itemsets are generated. A new scan is performed on the database in order to verify the support of each 2-candidate. The process goes on until no new frequent itemset is found.

Schema Mining: Problem Statement

Firstly introduced by [START_REF] Wang | Discovering typical structures of documents: A road map approach[END_REF][START_REF] Wang | Discovering structured association of semistructured data[END_REF], schema mining has attracted a lot of attention [START_REF] Asai | Efficient substructure discovery from large semi-structured data[END_REF][START_REF] Termier | Dryade: A new approach for discovering closed frequent trees in heterogeneous tree databases[END_REF][START_REF] Yan | CloseGraph: Mining closed frequent graph patterns[END_REF][START_REF] Zaki | Efficiently mining frequent embedded unordered trees[END_REF]. The problem is much more complicated than the classical association rule one since complex structures in the form of labelled hierarchical partially ordered data have to be considered.

The problem could be formulated as follows. A rooted labelled tree T = (V, E) is a direct, acyclic, connected graph with V = {0, 1, . . . , n} as the set of vertices (nodes), E = {(x, y)|x, y ∈ V } as the set of edges. We assume that there is a special vertex r ∈ V designated as a root and for all x ∈ V , there is a unique path from r to x. Then if x, y ∈ V and if there is a path from x to y then x is called an ancestor of y (i.e. y is a descendant of x). If the length of the path from two vertices x, y is reduced to one, then the ancestor relationship is considered as a parent relationship. For an internal node x ∈ V , we assume that its children x 1 , x 2 , . . . , x n (n ≥ 0) are ordered from left to right (i.e. there is a sibling relationship between children).

Let denote supp(S) the number of occurrences of the subtree in a tree T . The support of a subtree in a database D is defined as the number of trees in D that contain at least one occurrence of S. A subtree is frequent if its support is greater than or equal to a user-specified minimum support (minSupp) value. The schema mining problem is thus reduced to mine all subtrees in D where the support value holds.

The core of the process for mining frequent subtrees is usually quite similar to the association rule approach and is briefly described in algorithm 1. The Gen Cand function builds candidates. If k = 1 then all nodes are considered as potential subtree candidates. Otherwise, candidates of size k (i.e. having k nodes) are built by considering frequent subtrees of size k -1. At each step a scan is performed over the database in order to verify if the minSupp constraint holds. The process goes on until no new frequent subtree is found. Refer to [START_REF] Asai | Efficient substructure discovery from large semi-structured data[END_REF][START_REF] Zaki | Xrules: An effective structural classifier for XML data[END_REF] for details.

Why considering Fuzzy Approaches

From the problem statement presented so far, one problem remains about the subtree inclusion. In [START_REF] Kilpelinen | Tree matching problems with applications to structured text databases[END_REF], the author defines ten different tree inclusion relations and in [START_REF] Termier | Extraction of Frequent Trees in an Heterogeneous Corpus of Semi Structured data: Application to XML documents mining[END_REF], different definitions are proposed for strict, exact, ordered and weak tree inclusion. Let us have a closer look on these inclusion differences.

• The strict constraint addresses the labels of trees. When two nodes of dif- k++; until F does not grow any more; ferent trees have similar labels then the inclusion is called strict.

Algorithm 1: Mining Frequent SubTrees Data: Tree Database D Result: Frequent Subtrees F F ← ∅; k ← 1; repeat C k ← Gen Cand(k); foreach c ∈ C k do c.cpt ← 0; foreach T ∈ D do if c
• When the parent relationship is preserved, the inclusion is exact and the subtree is called induced. Otherwise, i.e. when the ancestor relationship is preserved, the inclusion is not exact and the subtree is called embedded. • If the children order is preserved the inclusion is ordered.

• When two labels of a tree correspond to a unique label into the other tree, the inclusion is called weak. According to the previous inclusion definitions, lot of algorithms were proposed for mining subtrees. Generally they consider strict and not weak inclusions. For example, in [START_REF] Wang | Discovering typical structures of documents: A road map approach[END_REF][START_REF] Wang | Discovering structured association of semistructured data[END_REF], the authors propose a new approach for discovering frequently occurring subtrees in XML documents. This algorithm is an extension of the level-wise Apriori approach where induced subtrees are considered. The FreqT approach proposed in [START_REF] Asai | Efficient substructure discovery from large semi-structured data[END_REF][START_REF] Asai | Discovering frequent substructures in large unordered trees[END_REF] also considers induced and ordered trees and an Apriori-like algorithm is proposed. In [START_REF] Zaki | Efficiently mining frequent trees in a forest[END_REF], TreeMiner was defined to mine labelled, embedded, and ordered subtrees. This work was extended in order to build a structural classifier for XML data [START_REF] Zaki | Xrules: An effective structural classifier for XML data[END_REF]. In [START_REF] Zaki | Efficiently mining frequent embedded unordered trees[END_REF], the author introduces Sleuth for mining frequent, unordered and embedded subtrees. The TreeFinder algorithm [START_REF] Termier | TreeFinder: a first step towards XML data mining[END_REF] also considers embedded trees and its extension Dryade [START_REF] Termier | Dryade: A new approach for discovering closed frequent trees in heterogeneous tree databases[END_REF] considers closed frequent subtrees, i.e. trees not any more included in another one. Other approaches such as gSpan [START_REF] Yan | gSpan: Graph-based substructure pattern mining[END_REF] and closeGraph [START_REF] Yan | CloseGraph: Mining closed frequent graph patterns[END_REF] generalize the problem by mining subgraphs.

These works mainly focus on the database scan and on data representation in order to propose scalable algorithms. Furthermore all the existing approaches consider crisp inclusion when mining frequent subtrees, i.e. a subtree is included or not into another one, which is too restrictive from a semantic point of view. In order to avoid this crisp inclusion problem, we propose to fuzzify this definition in order to better describe the data available on the Web. It is indeed very interesting to mine approximate schemas in order to have a better idea of the semantic structures we are provided with.

In a general framework, some works have been proposed in order to integrate fuzzy mining methods with the web [START_REF] Arotaritei | Web mining: a survey in the fuzzy framework[END_REF][START_REF] Pal | Web mining in soft computing framework: Relevance, state of the art and future directions[END_REF]. In the more precise framework of tree and semi-structured data mining, some works have been proposed, especially to deal with approximate tree matching. [START_REF] Ceravolo | Knowledge extraction from semi-structured data based on fuzzy techniques[END_REF] introduces a method to compare XML documents based on fuzzy bags. In [START_REF] Ceravolo | Fuzzy techniques for metadata construction[END_REF][START_REF] Ceravolo | Mining class hierarchies from XML data: Representation techniques[END_REF], XML data are mined to build clusters of XML data, thus defining a distance measure between XML documents. [START_REF] Wang | Approximate graph schema extraction for semi-structured data[END_REF] proposes an approximate method for graph schema extraction. Approximation is achieved by summarizing the semi-structured data graph using an incremental clustering method.

However as far as we know, no fuzzy method for mining frequent subtrees from huge databases has been proposed. We propose thus:

• to define the notion of fuzzy tree inclusion within this framework and • to mine fuzzy frequent patterns.

As shown on algorithm 1, a tree candidate has to be fully included within a tree from the database to be considered when computing its support. In this work, we argue that a tree may be considered to a certain extent within this computation even if it not fully included. For this purpose, we define four ways to consider fuzzy inclusion of a tree within another one:

• The first way aims at considering the vertical paths of trees. Contrary to induced one, embedded inclusion allows us to soften the ancestor-descendant relationship. Nevertheless, no methods exist to control this degree of relationship. • In the second way, we address the horizontal paths of trees. While classical approaches consider that sibling nodes are ordered or not, there is no way to consider the proportion of sibling nodes included. • The third way generalizes the previous one by considering the proportion of nodes. • The last one considers similarities between nodes.

These definitions are then integrated within a new algorithm (algorithm 2) which mines frequent subtrees in a fuzzy way. They aim at describing the extent to which a tree is included in another one. While classical approaches deal with crisp inclusion, meaning that a tree is or is not included within another one, we propose to use a degree, defined between 0 and 1. The four ways this degree can be obtained are described below. It corresponds to four ways the function F uzzy Inclusion Degree from algorithm 2 may be defined.

In algorithm 2, the degree cpt of a candidate is updated after each tree scanning by the Compute new cpt function. Several possibilities are given when merging these degrees, depending on the way the counting is performed. These possibilities are described in algorithm 3 by considering several counting methods: sigma-count, thresholded-count, thresholded sigma-count. This method is supposed to be known, as well as the associated threshold value ω if necessary.

For all these four ways fuzzy inclusion may be defined. We consider a node by node computation as done in the classical approach. However, this computation does not aim at matching exactly one node from the candidate tree to another one from the database tree but rather to compute to which extend a node corresponds to another one. This computation node by node leads to a degree deg returned by a F uzzy M atch function. This function takes into account the node n being considered, the tree S from which this node is originated, and the target tree T . The tree S has to be considered in order to keep the tree structure, and not only the nodes without any connections. We have thus to deal with a set of fuzzy degrees, ranging between 0 and 1 to be merged An OWA operator of dimension n is a mapping F : R n → R that has an associated n vector W = (w 1 , w 2 , . . . , w n ) T such that w i ∈ [0, 1] and n i=1 w i = 1. We have F (a 1 , a 2 , . . . , a n ) = n j=1 w j • b j where b j is the j th largest value of the a i .

In this framework, it is thus possible to compute the final degree in multiple ways. In the most pessimistic way, the minimum is computed (logical and), representing the lower boundary. In the very optimistic way, the maximum is computed (logical or), representing the upper boundary. Between these two boundaries, all possibilities are offered by tuning the weights, defining several degrees of andness/orness. The measure of orness is defined as follows:

orness(W ) = 1 n -1 n i=1 ((n -i) • w i )
The measure of andness is defined as andness(W ) = 1orness(W ).

Therefore, the user can choose whether he wants all nodes to be similar, some of them (using any desired degree of andness/orness), or at least one in a very flexible manner. Note that all the OWA operators are not associative. This is the reason why the computation of the degree d in algorithm 4 can only be done after the whole scan of the database. In some cases, this degree is computed using an associative operator (e.g. minimum), which allows for computing it on the fly. This way, it may be possible to cut off the calculus. For instance, when computing the minimum, the process can be stopped as soon as the degree becomes lower than the minimum threshold, if one threshold is considered in algorithm 3.

Fuzzy Vertical Paths

Let us consider fuzzy indirect links within trees. For instance, Figure 4 shows a tree S which is embedded within a tree T . When taken into account, embedded trees are counted within the final support whatever the ancestor-descendant relationship may be. There is no consideration of the number of nodes separating the ancestor node from the descendant one. Even if it is of great interest to consider two nodes as being related even if they are not directly linked, we argue that it becomes irrelevant to consider that two nodes are related if they are separated by a high number of nodes. In our approach, we propose thus to give a scope to this ancestor-descendant relationship. This scope is defined considering the number of nodes between ancestor and descendant nodes. Since it does not make sense to consider crisp boundaries, we propose to consider a fuzzy scope for the ancestor-descendant relationship. We consider fuzzy membership functions describing the ancestordescendant relationship depending on the number of nodes separating the two nodes being considered, as shown on Figure 5 when considering no more than five nodes in a fuzzy way. For this purpose, we use fuzzy quantifiers [START_REF] Ying | Quantifiers, modifiers and qualifiers in fuzzy logic[END_REF][START_REF] Zadeh | A computational approach to fuzzy quantifiers in natural languages[END_REF]. Fuzzy set theory enables to represent fuzzy quantifiers, such as at least 2, most, by membership functions (as illustrated by Figure 5 and Figure 7). In this framework, absolute quantifiers like at least 2 are distinguished from relative ones (most for instance).

5. Fuzzy Ancestor-Descendant Degree. Absolute Fuzzy Quantifier.

The F uzzy Inclusion Degree from algorithm 2 is then computed by merging the degrees to which nodes can be matched. This matching degree results from the extend to which a node may be linked to another one within a user-defined fuzzy ancestor-descendant relationship.

Fuzzy Horizontal Paths

In this section, we now consider fuzzy level inclusion. When considering ordered trees in the crisp approaches, a subtree S is included within a tree T only if all nodes of S can be mapped to nodes of T in the same order. In our approach, we propose to soften this definition by considering the proportion of nodes being included and well-ordered. For instance, Figure 6 shows an ordered tree S which is not embedded within an ordered tree T in the classical approaches since one of the node is misordered. In our approach, we consider S as being included in T with a certain degree since the other nodes satisfy the inclusion. This degree is reported in the F uzzy Inclusion Degree used in the algorithm. It depends on the proportion of nodes that are well-ordered within the nodes of the level. We consider thus a fuzzy membership function defining the degree (ranging from 0 to 1) to which most nodes are well-ordered, as illustrated by Figure 7.

Partial Inclusion

In this section, we consider partial node inclusion. In a general way, all the nodes of a subtree S must be included in a tree T if S is included in T . However, we argue that this is too restrictive when mining data from the real world when imperfections are often present. For instance, Figure 8 shows a tree S having 75% of its nodes included in T . In our approach, we propose thus to define partial inclusion by considering the proportion of nodes of S being present in T . This proportion is then related to a fuzzy quantifier which estimates to which extent the node has to be taken into account when computing the function F uzzy Inclusion Degree.

Note that in the classical case, mining totally included trees allows to cut in the database scan since whenever a node cannot be matched, there is no need looking for the other ones. In our approach, outliers are accepted, which may be considered as a drawback considering scalability. However, it is still possible to cut off the search when the proportion has been overpassed.

Similar Trees

In this section, we consider fuzzy similarities [START_REF] Bouchon | Towards general measures of comparison of objects[END_REF]. These similarities are related to some knowledge about the domain we are dealing with. For instance, when dealing with biological data, it is possible to know to which extent the node labelled by some bacteria is similar to the node labelled by another one. Let us now consider that we are provided with this kind of knowledge on the data. One of drawbacks of schema mining approaches is that the inclusion detection is only performed on nodes having same labels (C.f. Section 3strict inclusion). According to the semantic Web point of view, this restriction is usefulness since two different labels could describe similar concepts.

In our approach, we consider fuzzy approaches to overcome this drawback. For instance, Figure 8 shows a tree S that should be matched with T if we know that the concept D is close to B. For this purpose, we consider fuzzy ontologies as done in [START_REF] Parry | A fuzzy ontology for medical document retrieval[END_REF], and more generally fuzzy relations [START_REF] Klir | Fuzzy Sets, Uncertainty, and Information[END_REF] describing to which extent two nodes are similar.

Let us recall that a fuzzy relation R between two reference sets X and Y is defined as a fuzzy subset from X × Y . In particular, when X and Y are finite, the relation can be described by a matrix M (R) including the values of its membership function.

It is possible to compute the transitive closure of a fuzzy relation. Let G be the graph associated to the fuzzy relation R and let f R be its membership function. The transitive closure of R is the R T such that f R T (x, y) = 0 iff there exists a path in G from the node x to the node y.

Let R T be the max-min transitive closure of R and R • R be the max-min composition of R and R. R T is computed as follows:

(

1) R = R (2) While R changes • R = R • R = R (R • R) (3) Stop. R T = R
In our approach, we consider a fuzzy relation R defined on the universe of node labels N . We have thus R : N × N → [0, 1] describing to which extent a node n ∈ N is similar to a node n ∈ N . By computing the transitive closure of R, all nodes can be compared to every other nodes. Symmetric relations are considered since we deal with similarities. This degree of similarity is taken into account in the F uzzy M atch function when computing the F uzzy Inclusion Degree, as illustrated by algorithm 5. In this algorithm, the tree structure preservation is checked. This verification is required in order to prevent the tree structure from being destroyed. This verification must thus take into account the previous nodes that have been matched together.

n ← f R (n, n ); d ← OW A |T | n =1 deg n ; return d;
Algorithm 5 proposes to compute the final degree as a merged value from all the similarity degrees with nodes from T . In the classical case, the algorithm is more efficient since there is no need to scan all nodes from T . Another possibility is thus to stop the scan of T as soon as a similar node n ∈ T is found, leading to algorithm 6. In this case, we consider a user-defined γ threshold.

(n, n ) ≥ γ then return f R (n, n );
Fuzzy morphisms may also be considered in order to deal with this problem [START_REF] Perchant | Fuzzy morphisms between graphs[END_REF][START_REF] Kohout | Defining homomorphisms and other generalized morphisms of fuzzy relations in monoidal fuzzy logics by means of bk-products[END_REF].

Fuzzy Frequent Subtrees

Classical approaches mine frequent trees that do not provide much information about the occurrences within the database except the fact that these subtrees are embedded in a sufficient number of trees from the database (depending on the user-defined minimum support value). We propose thus to extend the knowledge on the mined semantic structures by providing fuzzy links within the frequent trees. These fuzzy links are more or less thick. They help knowing whether they are very shared, middle shared or a little shared, as illustrated on Figure 9. In this framework, frequent subtrees are fuzzy acyclic graphs. Each link between two nodes is labelled by a degree ranging from 0 to 1. This degree is obtained by considering the proportion of trees from the database that contain the link. This proportion is then mapped to fuzzy quantifiers such as very shared, middle shared or little shared, as shown by Figure 10. It is graphically represented by taking into account the fuzzy quantifier having the greatest membership degree. Each fuzzy quantifier is then represented by a different thickness.

Conclusion

Recently mediator-based integration of heterogeneous data sources over the web has shown that it is a very efficient way for dealing with different types of Web documents such as XML ones. One of the main advantages of such approaches is that an heterogeneous collection of XML documents could be integrated into mediator schemas in order to be processed or queried in a uniform way. However, we have seen that no complete automatic tool is available to extract semantic knowledge from these large amounts of distributed and heterogeneous data. It is indeed still hardly possible to automatically build mediator schemas, which are thus still hand-crafted. Even if schema mining approaches exist to extract in an efficient way the commonly occurring schemas from a collection, they suffer from different drawbacks. The main drawback is that they consider crisp inclusion when mining frequent subtrees, i.e. a subtree is included or not into another one. In the semantic web context this constraint is too restrictive. In this paper we have thus proposed to soften such an inclusion in order to better describe the data available on the Web.

Data mining is of great interest for the semantic Web. In coming years, the Internet is likely to become one of the principal application areas of fuzzy logic. In this framework, mining fuzzy XML mediator schemas is crucial (i) in order to be automatically provided with a mediator schema to query data and (ii) in order to get semantic structures from the huge amounts of distributed and heterogeneous data. The approach we propose here can be used in order to mine frequent web structures from web sites. Moreover, mediator schemas, fuzzy or not, can be used in the framework of fuzzy queries over web data.

Fig. 1 .

 1 Fig. 1. Querying Data Distributed over the Web

Figure 3 Fig. 3 .

 33 Figure3gives examples of a tree S included in a tree T .

  Degree Function. Data: S candidate tree, T target tree Result: Fuzzy Degree d ∈ [0, 1] foreach node n ∈ S do deg n ← F uzzy M atch(n, S, T ); d ← OW A |S| n=1 deg n ; return d;

Fig. 4 .

 4 Fig. 4. Fuzzy Vertical Paths

Fig. 6 .

 6 Fig. 6. Fuzzy Horizontal Paths

Fig. 7 .

 7 Fig. 7. Relative Fuzzy Quantifier: Proportion of Nodes

Fig. 8 .

 8 Fig. 8. Partial Inclusion

Algorithm 5 :

 5 Computing the Fuzzy Node Matching. F uzzy M atch Function. Data: n node, S candidate tree, T target tree Result: Fuzzy Degree d ∈ [0, 1] foreach node n ∈ T do if n matches n' and respects the S tree structure then deg

Fig. 9 .

 9 Fig. 9. Fuzzy Tree

Fig. 10 .

 10 Fig. 10. Fuzzy Quantifiers for Fuzzy Trees

Algorithm 6 :

 6 Fuzzy Node Matching with Stop Criterion. F uzzy M atch Function. Data: n node, S candidate tree, T target tree Result: Fuzzy Degree d ∈ [0, 1] foreach node n ∈ T do /* γ stands for a user-specified threshold */ if n matches n' and respects the S tree structure and f R

In order to provide the user with a large set of solutions, we consider the OWA operators (Ordered Weighted Aggregator) [START_REF] Yager | Families of owa operators[END_REF], as illustrated by algorithm 4.


J Duan

H Huang

B Dong

D Jung

Z Zhang

J C Norman

J E Bowers

F Grillot

Thermally insensitive determination of the chirp parameter of InAs/GaAs quantum dot lasers epitaxially grown onto silicon

Keywords: Quantum dots, Semiconductor lasers, Chirp parameter;

A common way of extracting the chirp parameter (i.e., the α-factor) of semiconductor lasers is usually performed by extracting the net modal gain and the wavelength from the amplified spontaneous emission (ASE) spectrum. Although this method is straightforward, it remains sensitive to the thermal effects hence leading to a clear underestimation of the α-factor. In this work, we investigate the chirp parameter of InAs/GaAs quantum dot (QD) lasers epitaxially grown on silicon with a measurement technique evaluating the gain and wavelength changes of the suppressed side modes by optical injection locking. Given that the method is thermally insensitive, the presented results confirm our initial measurements conducted with the ASE i.e. the α-factor of the QD lasers directly grown on silicon is as low as 0.15 hence resulting from the low threading dislocation density and high material gain of the active region. These conclusions make such lasers very promising for future integrated photonics where narrow linewidth, feedback resistant and low-chirp on-chip transmitters are required.

INTRODUCTION

Silicon photonics is a promising key enabling technology bringing multiple advantages to overcome the data transmission bottleneck of semiconductor microchips in terms of bandwidth limit and power density [START_REF] Sun | Single-chip microprocessor that communicates directly using light[END_REF]. The integration of multiple optical functions on a microelectronics chip brings many innovative perspectives, along with the possibility to enhance the performance of photonics integrated circuits (PICs). Owing to the atom-like discrete energy levels, quantum dots (QDs) allow producing energy-and cost-efficient devices with outstanding properties such as low threshold current, good temperature stability and narrow linewidth [START_REF] Crowley | GaAs based quantum dot lasers[END_REF][START_REF] Duan | Narrow spectral linewidth in InAs/InP quantum dot distributed feedback lasers[END_REF]. In particular, direct epitaxial growth of GaAs layers onto silicon with InAs QDs as gain media has been proved to be a meaningful solution for inexpensive and monolithically integrated silicon light emitters [START_REF] Norman | Perspective: The future of quantum dot photonic integrated circuits[END_REF][START_REF] Nishi | Development of Quantum Dot Lasers for DataCom and Silicon Photonics Applications[END_REF]. In semiconductor lasers, the α-factor is a key parameter resulting from the phase-amplitude coupling effect, and driving the spectral linewidth [START_REF] Duan | Narrow spectral linewidth in InAs/InP quantum dot distributed feedback lasers[END_REF], the sensitivity to optical feedback [START_REF] Otto | Complex dynamics of semiconductor quantum dot lasers subject to delayed optical feedback[END_REF][START_REF] Huang | Analysis of the optical feedback dynamics in InAs/GaAs quantum dot lasers directly grown on silicon[END_REF], the nonlinear dynamics under optical injection [START_REF] Wieczorek | The dynamical complexity of optically injected semiconductor lasers[END_REF] and four-wave mixing generation [START_REF] Huang | Non-degenerate four-wave mixing in an optically injection-locked InAs/InP quantum dot Fabry-Perot laser[END_REF]. The α-factor typically describes the coupling between the carrierinduced variation of real and imaginary parts of susceptibility and is defined as [START_REF] Osinski | Linewidth broadening factor in semiconductor lasers -An overview[END_REF]:

𝛼 = - $% & '(/'* '+/'* (1)
where dn and dg are the small refractive index and modal gain variations accompanied by a carrier density variation dN and 𝜆 is the lasing wavelength.

At the system level, a large α-factor leads to a frequency chirp (∆ν) under direct modulation, which is expressed as [START_REF] Koch | Effect of nonlinear gain reduction on semiconductor laser wavelength chirping[END_REF]:

∆ν = - 2 $% ( ' '3 𝑙𝑛𝑃 + 89:; < =>? @AB ) (2)
where P is the optical power, 𝛤 is the modal confinement factor, 𝜖 characterizes the gain compression factor, 𝑉 FG3 is the active layer volume and 𝜂 is the total differential quantum efficiency. The frequency chirp typically limits the maximum data rate and transmission distance over a dispersive fiber. Our recent work has previously reported α-factors as low as 0.13 in QD lasers directly grown onto silicon as a result of the low threading dislocation density [START_REF] Duan | Semiconductor quantum dot lasers epitaxially grown on silicon with low linewidth enhancement factor[END_REF]. The α-factors is usually measured from the standard amplified spontaneous emission (ASE) which relies on simultaneous observation of the longitudinal modes wavelength shift and the modal gain changes with varying current below lasing threshold [START_REF] Crowley | GaAs based quantum dot lasers[END_REF]. However, due to residual thermal effects, the ASE measurement can possibly underestimate the α-factors. In order to confirm our initial measurements [START_REF] Duan | Semiconductor quantum dot lasers epitaxially grown on silicon with low linewidth enhancement factor[END_REF], this work investigates the α-factor with a thermally insensitive method based on the analysis of residual side-mode dynamics under optical injection. Such method has already been successfully tested in InAs/InP QD lasers [START_REF] Wang | Thermally insensitive determination of the linewidth broadening factor in nanostructured semiconductor lasers using optical injection locking[END_REF]. Finally, results confirm that the α-factors initially measured in [START_REF] Duan | Semiconductor quantum dot lasers epitaxially grown on silicon with low linewidth enhancement factor[END_REF] are low with a minimum value of 0.15 at the gain peak. Given that the standard deviation does not exceeding 15%, we can conclude that both methods are found in a relative good agreement, which is obviously important when small values of the linewidth broadening factor are addressed. Overall, we believe that this work paves the way of future integrated optical transmitters on silicon with increased maximum transmission capabilities.

QD LASER STRUCTURE

The InAs/GaAs Fabry-Perot (FP) laser is 1.35 mm long with a 3.5 µm wide ridge waveguide directly grown onto an onaxis (001) Si wafer in a Veeco Gen-II molecular beam epitaxy chamber. The full laser epitaxial structure and QD growth conditions can be found elsewhere [START_REF] Jung | High efficiency low threshold current 1.3 µm InAs quantum dot lasers on on-axis (001) GaP/Si[END_REF]. A schematic of the QD laser structure is presented in Figure 1. The p-modulation doped QD laser consists of 5 QD layers spaced by 37.5 nm thick GaAs barrier layers, where the first 10 nm GaAs layer was undoped, followed by a 10 nm p-GaAs layer at a target hole concentration of 5 × 10 17 cm -3 using Be, and the final 17.5 nm GaAs layer was undoped again to complete the GaAs barrier. The QD density is 4.9 × 10 10 cm -2 and the photoluminescence spectra of the full laser sample revealed a very small full-width at half-maximum (FWHM) of 29.3 meV from the ground state (GS) peak (not shown), indicating a highly homogeneous InAs QD size throughout the 5 QD stacks. The front facet has a power reflectivity of ~60 % whereas a high reflectivity of ~99 % is used on the rear and the reduced threading dislocation Please verify that (1) all pages are present, (2) all figures are correct, (3) all fonts and special characters are correct, and (4) all text and figures fit within the red margin lines shown on this review document. Complete formatting information is available at http://SPIE.org/manuscripts Return to the Manage Active Submissions page at http://spie.org/submissions/tasks.aspx and approve or disapprove this submission. Your manuscript will not be published without this approval. Please contact author_help@spie.org with any questions or concerns. density (TDD) of 7.3 × 10 6 cm -2 in the active region has been achieved. At room temperature (20 ℃), the laser has a lasing threshold of 26.5 mA, and the GS gain peak was located around 1300 nm at 80 mA. The light-current characteristic curve is shown in Figure 2, where the marker in black corresponds to the bias current of 80 mA used to measure the optical spectrum depicted in the inset.

EXPERIMENTAL SETUP AND RESULTS

In a prior work [START_REF] Duan | Semiconductor quantum dot lasers epitaxially grown on silicon with low linewidth enhancement factor[END_REF], the α-factor was initially measured by the ASE method, which relies on measuring the differential index by tracking the frequency shift of the longitudinal FP modes whereas the differential gain is obtained by measuring the net modal gain from the FP modulation depth. The definition of α-factor in ( 1) is then equivalently expressed by:

𝛼 = - 8% IJ& '&/'K 'L MN? /'K (3)
where L is the cavity length, 𝛿𝜆 is the free spectral range while 𝑑𝜆 and 𝑑𝐺 (R3 are respectively the mode wavelength and net modal gain variations caused by the bias current change dI. However, it is important to stress that the main weakness of the ASE method results from its poor reliability if thermal effects are not properly eliminated. Indeed, thermal effects cause a red-shift of the optical wavelength and thereby interfere with the tracking of the carrier-induced blue-shift effect.

In order to fully eliminate the thermal effects, a current source operated in the pulsed mode is highly preferred although it is also possible to use continuous-waves operation and eliminate the thermal effects by applying a correcting factor. Thus, the wavelength red-shift due to thermal effects, measured by varying the pump current right above threshold, is subtracted from the wavelength blue-shift measured below threshold. Following this protocol, the fitting error can certainly be minimized, however, the thermal correction holds under the assumption that thermal effects in QD lasers are unchanged below and above threshold as well as that the carriers are clamped above threshold. Nevertheless, the latter assumption is somewhat disputable for QD lasers, as the carrier population in resonant states keeps increasing with the bias currents, which will continue to enhance the gain and change the refractive index.

In what follows, the lasers under test are biased with a pulsed current using the shortest pulse width of 0.1 µs. The peak wavelength is then recorded by varying the duty cycle from 0.1% to 10 % with a step of 2 %. For each sub-threshold bias current, the peak wavelength values are plotted as a function of duty cycle and the extrapolation at 0 % duty cycle allows extracting the corresponding values without thermal effects. Figure 3(a) displays one peak wavelength measured for the Please verify that (1) all pages are present, (2) all figures are correct, (3) all fonts and special characters are correct, and ( 4) all text and figures fit within the red margin lines shown on this review document. Complete formatting information is available at http://SPIE.org/manuscripts

Return to the Manage Active Submissions page at http://spie.org/submissions/tasks.aspx and approve or disapprove this submission. Your manuscript will not be published without this approval. Please contact author_help@spie.org with any questions or concerns.

10939 -27 V. as a function of the bias currents for two duty cycles. The wavelength red-shift and blue-shift are clearly evidenced for the cases of 10% and 0% duty cycle. In addition, it is noted that the pulse width of 0.1 µs corresponds to the lower bound of the pulsewidth compared with a few microseconds usually used in many studies [START_REF] Tan | Dependence of Linewidth Enhancement Factor on Duty Cycle in InGaAs-GaAs Quantum-Dot Lasers[END_REF]. Besides, it is also important to stress that going below 0.1 µs would result in an extremely low output power and the impossibility to retrieve the gain and then the α-factor at the price of a poor accuracy. Therefore, a pulse width of 0.1 µs is typically what is recommended to almost eliminate all the thermal effects in the material. In order to investigate the influence of thermal effects on the α-factor, the pulse width dependence of the measured α-factor is also performed on a p-doped QD laser which has 3 QD layers and a 1.5 mm long cavity with a 3 µm wide ridge waveguide. As shown in figure 4(a), this laser has initial α-factor of about 0.3

Please verify that (1) all pages are present, (2) all figures are correct, (3) all fonts and special characters are correct, and (4) all text and figures fit within the red margin lines shown on this review document. Complete formatting information is available at http://SPIE.org/manuscripts Return to the Manage Active Submissions page at http://spie.org/submissions/tasks.aspx and approve or disapprove this submission. Your manuscript will not be published without this approval. Please contact author_help@spie.org with any questions or concerns.

10939 -27 V. which is found a bit larger than that of 5 QD laser due to the smaller material gain. However, when further increasing the pulse width from 0.1 µs to 5 µs and 10 µs, α-factor values extracted at a fixed duty cycle of 0 % are found reduced from 0.3 to negative values down to -1 at 5 µs and to -1.8 at 10 µs. This effect is attributed to the thermal effect induced wavelength red-shift below threshold which results in the negative α-factor values. Figures 4(b) also displays the measured α-factors extracted at a fixed pulse width of 0.1 µs as a function of different duty cycle ranging from 0 % to 10 % with a step of 2 %. The α-factors decrease from 0.3 to -0.35 along with the enhancement of the duty cycle from 0 % to 10 %. Therefore, the α-factors can be strongly underestimated if the thermal effects are not completed eliminated from the ASE method. In addition, the laser biased in the pulsed mode also resulting in a low signal-to-noise ratio and even irregular spectral line shape. Now, the α-factor is investigated with a thermally insensitive method using the residual side-mode dynamics under optical injection locking. Figure 5 illustrates the experimental setup where a tunable external cavity laser was used as the master laser and the injected light was coupled into the QD laser (slave laser) through an optical circulator. A polarization controller is inserted to align the polarization of both lasers. The 90 % of the output power was captured by a high resolution optical spectrum analyzer (OSA) and the rest of 10 % was monitored by a power meter. The temperature of the slave laser was maintained at 20 ℃ throughout the experiment. The optical injection is performed far from the gain peak e.g. to a

Please verify that (1) all pages are present, (2) all figures are correct, (3) all fonts and special characters are correct, and (4) all text and figures fit within the red margin lines shown on this review document. Complete formatting information is available at http://SPIE.org/manuscripts Return to the Manage Active Submissions page at http://spie.org/submissions/tasks.aspx and approve or disapprove this submission. Your manuscript will not be published without this approval. Please contact author_help@spie.org with any questions or concerns.

10939 -27 V. longitudinal mode around 1315 nm, and the wavelength detuning is operated within the stable-locked regime. Owing to the gain reduction, the side modes are deeply suppressed in such way that the residual power of the side modes at the gain peak become governed by the ASE, and therefore the α-factor can be extracted as discussed hereinafter. Figure 6 illustrates the effect of optical injection on the FP cavity modes (blue lines) biased at 27.5 mA right above the threshold. Once the laser is stably-locked, only the mode subject to the optical injection keeps lasing (at 1315 nm) while all side modes are greatly suppressed (red lines). The frequency detuning (F inj ) is defined as the optical frequency difference between the master and slave lasers, while the injection ratio (R inj ) is the ratio of the injected signal power (P inj ) to the free-running slave laser power (P 0 ). Both the enhanced injection ratio and increased frequency detuning reduce the entire net modal gain compared to that of the free-running laser hence resulting in the suppression of the longitudinal side modes. The optical injection induced gain reduction can be understood through the corpuscular equations of the injection-locked oscillators, which under the steady-state conditions leads to the following expression [START_REF] Wang | Thermally insensitive determination of the linewidth broadening factor in nanostructured semiconductor lasers using optical injection locking[END_REF]:

∆𝐺 (R3 = - 8S > B T 𝑅 V(W 𝑃 X /𝑃 Y cos 𝜑 (4)
where the phase 𝜑 of the laser field is given by: 𝜑 = sin `a -8%b cMd S > e cMd : f /: g (ah2 i )

-tan `a𝛼 [START_REF] Nishi | Development of Quantum Dot Lasers for DataCom and Silicon Photonics Applications[END_REF] with 𝑘 G is the coupling coefficient of the master laser to the slave laser, 𝜈 + is the light group velocity and P S is the intracavity power of the slave laser under optical injection. The above equation shows that the F inj reduces the net modal gain through the phase of the laser field. This gain reduction behavior is similar to that induced by the bias current in the ASE method. Therefore, the equation ( 3) can be re-expressed as follows:

𝛼 = - 8% IJ& '&/'& n 'L MN? /'& n (5)
with 𝜆 o being the injection wavelength of the master laser. Consequently, an enhancement of the wavelength detuning shifts the lasing mode towards the longer wavelength side and reduces the gain to a lower level. Figure 7 describes the wavelength dependence of the net modal gain spectra for the p-doped QD laser. With the increase of the wavelength detuning, the whole gain spectra fall down gradually. Since this method is based on analyzing the residual side-mode ASE spectra under injection locking (IL), we named it ASE IL. In ASE IL method, the laser is biased at 27.5 mA just above the threshold current. The injected light is located at 1315 nm, which is 20 nm away from the gain peak of 1295 nm, while the R inj is fixed at 5 dB. It is noted that this method is not sensitive neither to the injection power nor the choice of the injected mode due to the fact that the side modes are non-lasing and operated below threshold [START_REF] Wang | Thermally insensitive determination of the linewidth broadening factor in nanostructured semiconductor lasers using optical injection locking[END_REF]. As shown in the Figure 8, the α-factors measured by ASE IL are found in a relatively good agreement with that measured under ASE and increase from about 0.14 at 1293 nm to 0.19 at 1299 nm, even though they are slightly higher. For example, compared with the α-factor of 0.13 at 1295 nm measured by ASE method, the α-factor measured by ASE IL is found with value of 0.15 with a standard deviation not exceeding 15%. This slight difference is attributed to the fact that ASE IL method relies on a fixed bias current while the ASE method requires a variation of the bias current hence leading to persistent device heating and underestimation of the α-factor. Therefore, the α-factor values measured under ASE IL are well consistent with the ASE values hence showing that this method is reliable and compatible with QD lasers onto silicon.

CONCLUSIONS

To sum, this work investigates the extraction of the α-factors with a thermally insensitive method analyzing the residual side-mode dynamics under optical injection locking of QD lasers epitaxially grown onto silicon. The presented results perfectly confirm our initial measurements conducted with the standard ASE method i.e. the α-factor of the QD laser directly grown on silicon is very small and as low as 0.15 with a standard deviation not exceeding 15% which is very promising for the realization of integrated optical transmitters and increasing the maximum transmission distance in future silicon photonics systems.

Figure 1 .

 1 Figure 1. Cross-sectional schematic of an InAs QD laser epitaxially grown on (001) GaP/Si.

  10939 -27 V. 1 (p.2 of 8) / Color: No / Format: Letter / Date: 2/2/2019 12:54:09 AM SPIE USE: ____ DB Check, ____ Prod Check, Notes:

Figure 2 .

 2 Figure 2. The light-current characteristic of the laser measured at room temperature (20 ℃). The inset shows the optical spectrum at 80 mA (black marker).

  Figure 3. (a) Measured peak wavelength as a function of the duty cycle for different subthreshold bias currents. (b) Measured peak wavelength as a function of the bias current for two values of the duty cycle.

Figure 4 .

 4 Figure 4. (a) Pulse width dependence of the measured α-factor at a fixed duty cycle of 0%. (b) Duty cycle dependence of the measured α-factor at a fixed pulse width of 0.1 µs.

Figure 5 .

 5 Figure 5. Experimental setup for measuring the α-factors.

Figure 6 .

 6 Figure 6. Measured optical spectra for the p-doped QD laser. The blue lines are for the free-running laser (i.e. without optical injection), and the red lines are for the laser is injection-locked.

Figure 7 .

 7 Figure 7. Wavelength variation of the measured net modal gain with detuning of the injected light from shorter to longer wavelength within the stable-locking range.

Figure 8

 8 Figure 8 reveals the comparison of α-factors of the p-doped QD laser measured by ASE and ASE IL methods.In ASE IL method, the laser is biased at 27.5 mA just above the threshold current. The injected light is located at 1315 nm, which is 20 nm away from the gain peak of 1295 nm, while the R inj is fixed at 5 dB. It is noted that this method is not sensitive neither to the injection power nor the choice of the injected mode due to the fact that the side modes are non-lasing and operated below threshold[START_REF] Wang | Thermally insensitive determination of the linewidth broadening factor in nanostructured semiconductor lasers using optical injection locking[END_REF]. As shown in the Figure8, the α-factors measured by ASE IL are found in a relatively good agreement with that measured under ASE and increase from about 0.14 at 1293 nm to 0.19 at 1299 nm, even though they are slightly higher. For example, compared with the α-factor of 0.13 at 1295 nm measured by ASE method, the α-factor measured by ASE IL is found with value of 0.15 with a standard deviation not exceeding 15%. This slight difference is attributed to the fact that ASE IL method relies on a fixed bias current while the ASE method requires a variation of the bias current hence leading to persistent device heating and underestimation of the α-factor. Therefore, the α-factor values measured under ASE IL are well consistent with the ASE values hence showing that this method is reliable and compatible with QD lasers onto silicon.

Figure 8 .

 8 Figure 8. Wavelength dependence of the measured α-factor by ASE method and ASE IL method for p-doped QD laser. The horizontal dotted line indicates the gain peak α-factor value of 0.15, while the vertical dotted line indicates the position of FP gain peak.

  Please verify that[START_REF] Sun | Single-chip microprocessor that communicates directly using light[END_REF] all pages are present,[START_REF] Crowley | GaAs based quantum dot lasers[END_REF] all figures are correct, (3) all fonts and special characters are correct, and (4) all text and figures fit within the red margin lines shown on this review document. Complete formatting information is available at http://SPIE.org/manuscripts Return to the Manage Active Submissions page at http://spie.org/submissions/tasks.aspx and approve or disapprove this submission. Your manuscript will not be published without this approval. Please contact author_help@spie.org with any questions or concerns.

10939 -27 V. 1 (p.6 of 8) / Color: No / Format: Letter / Date: 2/2/2019 12:54:09 AM SPIE USE: ____ DB Check, ____ Prod Check, Notes:

ACKNOWLEDGMENTS

The authors acknowledge the financial support from ARPA-E (DE-AR000067) and the Institut Mines-Télécom.


Raabid Hussain
email: raabid.hussain@u-bourgogne.fr

Alain Lalande

Kibrom Roberto Marroquin

Berihu Girum

Caroline Guigou

Alexis Bozorg Grayeli

Kibrom Berihu Girum

Real-Time Augmented Reality for Ear Surgery

Keywords: Augmented reality, transtympanic procedures, otology, minimally invasive, image-guided surgery

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

During otologic procedures, when the surgeon places the endoscope inside the external auditory canal, the middle ear cleft structures, concealed behind the opaque tympanic membrane (TM), are not directly accessible. Consequently, surgeons can access these structures using the TM flap approach [START_REF] Gurr | Approaches to the middle ear[END_REF] which is both painful and exposes the patient to risk of infection and bleeding [START_REF] Marroquin | Augmented reality of the middle ear combining otoendoscopy and temporal bone computed tomography[END_REF].

Alternatively, transtympanic procedures have been designed which aim at accessing the middle ear cleft structures through a small and deep cavity inside the ear. These techniques have been used in different applications such as ossicular chain repair, drug administration and labyrinthe fistula diagnosis [START_REF] Dean | Eustachian Tube Dilation via a Transtympanic Approach in 6 Cadaver Heads: A Feasibility Study[END_REF][START_REF] Mood | Use of a microendoscope for transtympanic drug delivery to the round window membrane in chinchillas[END_REF].

The procedures offer many advantages: faster procedure, preservation of TM and reduced bleeding. However, limited operative space, field of view and instrument manoeuvring introduce surgical complications.

Our hypothesis claims that augmented reality (AR) would improve the procedure of middle ear surgery by providing instrument pose information and superimposing preoperative computed tomography (CT) image of the middle ear onto the micro/endoscopic video of TM. The key challenge is to enhance ergonomy while operating in a highly undersquared cylinderical workspace achieving sub-millimetric precision. To our knowledge, AR has not been applied to transtympanic and otoendoscopic procedures, thus the global perspective of the work is to affirm our hypothesis.

In computer assisted surgical systems, image registration plays an integral role in the overall performance. Feature extraction methods generally do not perform well due to the presence of highly textured structures and non-linear biasing [START_REF] Viergever | A survey of medical image registration under review[END_REF]. Many algorithms have been proposed specifically for endoscope-CT registration. Combinations of different intensity based schemes such as crosscorrelation, squared intensity difference, mutual information and pattern intensity have shown promising results [START_REF] Hummel | 2D/3D registration of endoscopic ultrasound to CT volume data[END_REF][START_REF] Yim | Registration of 3D CT data to 2D endoscopic image using a gradient mutual information based viewpoint matching for image-guided medialization laryngoplasty[END_REF]. Similarly, feature based schemes involving natural landmarks, contour based feature points, iterative closest point and k-means clustering have also been exploited [START_REF] Jun | Feature points based image registration between endoscope image and the CT image[END_REF][START_REF] Wengert | Markerless endoscopic registration and referencing[END_REF].

Different techniques involving learned instrument models, artificial markers, pre-known kinematic and gradient information using Hough transform have been proposed to identify instruments in video frames [START_REF] Haase | Laparoscopic instrument localization using a 3D time of flight/RGB endoscope[END_REF]. If the target is frequently changing its appearance, gradient based tracking algorithms need continuous template updating to maintain accurate position estimation. Analogously, reliable amount of training data is required for classifier based techniques. Although extensive research has been undertaken for identification of instruments in image plane, limited work has been accomplished to estimate 3D pose. Trained random forest classifier using instrument geometry as a prior and visual servoing techniques employing four marker points have been proposed [START_REF] Allan | Toward detection and localization of instruments in minimally invasive surgery[END_REF][START_REF] Nageotte | Visual servoing based endoscopic path following for robot-assisted laparoscopic surgery[END_REF]. Three point perspective framework involving collinear markers has been also suggested [START_REF] Liu | A portable 3d vision coordinate measurement system using a light pen[END_REF].

Our proposed approach initially registers the CT image with the microscopic video, based on fiducial markers. This is followed by a feature based motion tracking scheme to maintain synchronisation. The surgical instrument is also tracked and its pose estimated using 3-collinear-point framework.

Methodology

The system is composed of three main processes: initial registration, movement tracking and instrument pose estimation. The overall hierarchy is presented in Fig. 1. The proposed system has two main inputs. Firstly, the reconstructed image which is the display of the temporal bone, depicting middle ear cleft structures behind TM, obtained from preoperative CT data through OsiriX 3D endoscopy function (Pixmeo SARL, Switzerland). Secondly, the endoscopic video which is the real-time video acquired from a calibrated endoscope or surgical microscope during a surgical procedure. The camera projection matrix of the input camera was computed using [START_REF] Zhang | A flexible new technique for camera calibration[END_REF]. The calibration parameters are later used for 3D pose estimation.

There is low similarity between the reconstructed and endoscopic images (Fig. 2), thus the performance of intensity and feature based algorithms is limited in this case. Marroquin et al. [START_REF] Marroquin | Augmented reality of the middle ear combining otoendoscopy and temporal bone computed tomography[END_REF] established correspondence by manually identifying points in the endoscopic and CT images. However, accurately identifying natural landmarks is a tedious and time-consuming task. Thus six stainless steel fiducial markers (≈ 1mm in length) were attached around TM (prior to CT acquisition). Fig. 1: Overall workflow of the proposed system. A colour coding scheme, defined in (a), has been used to differentiate between different processes.

Initial Semi-automatic Endosocope-CT Registration

Since the intensity of fiducials is significantly higher than that of anatomical regions on CT images, contrast enhancement and thresholding is used to obtain fiducial regions. The centre of each fiducial is then obtained using blob detection. The user selects the corresponding fiducials in the first frame of the endoscopic video. There are very few common natural landmarks around TM, thus the fiducials ease up the process of establishing correspondence. In order to eliminate human error, similar pixels in a small neighbourhood around the selected points are also taken into account. A RANdom SAmple Consensus (RANSAC) based homography [START_REF] Hartley | Multiple view geometry in computer vision[END_REF] registration matrix H R , which warps the reconstructed image onto the endoscopic video, is computed using these point correspondences.

An ellipse shaped mask is generated using the fiducial points in the endoscopic frame. Since TM does not have a well-defined boundary in the endoscopic video, this mask is used as an approximation of TM. The mask is used in the tracking process to filter out unwanted (non-planar) features.

Endoscope-Target Motion Tracking

Speeded Up Robust Features (SURF) [START_REF] Bay | Surf: Speeded up robust features[END_REF] was employed in our system for tracking the movement between consecutive video frames [START_REF] Marroquin | Augmented reality of the middle ear combining otoendoscopy and temporal bone computed tomography[END_REF]. For an accurate homography, all the feature points should lie on coplanar surfaces [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. However, the extracted features are spread across the entire image plane comprising of the TM and auditory canal. The ellipse generated in previous step is used to filter out features that do not lie on TM (assumed planar). A robust feature matching scheme based on RANSAC and nearest neighbour (FLANN [START_REF] Muja | Fast approximate nearest neighbors with automatic algorithm configuration[END_REF]) frameworks is used to determine the homography transformation H T between consecutive frames. A chained homography framework is then used to warp the registered reconstructed image onto the endoscopic frame:

H i+1 = H T * H i , (1)
where H 0 is set as identity. H i can then be multiplied with H R to transform the original reconstructed image to the current time step. A linear blend operator is used for warping reconstructed image onto the current endoscopic frame.

3D Pose Estimation of Surgical Instrument

In surgical microscopes, small depth of focus leads to a degradation of gradient and colour information. Thus popular approaches for instrument identification do not perform well. Consequently, three collinear ink markers were attached to the instrument (Fig. 3). The three marker regions are extracted using thresholding. Since, discrepancies are present, a pruning step followed by blob detection is carried out to extract centres of the largest three regions. A linear Kalman filter is used to refine the marker centre points to eliminate any residual degradation. Since, the instrument may enter from any direction and protrude indefinitely, geometric priors are not valid. The proposed approach assumes that no instrument is present in the first endoscopic frame. The first frame undergoes a transformation based on H i . Background subtraction followed by pruning (owing to discrepancies in H i ) is used to extract the tool entry point in the frame boundary. The tool entry point is then used to associate the marker centres to marker labels B, C and D. The instrument tip location can then be obtained using a set of perception based equations that lead to:

a = 1 3 (b + c + d + AB CD (c -d) + AC BD (b -d) + AD BC (b -c)) , ( 2
)
where a is projection of the surgical instrument tip A on the 2D image frame, b, c and d are projections of the markers and alphabet pairs represent the physical distance between the markers. A three point perspective framework is then used to estimate the 3D pose of the instrument. Given focal length of the camera, known physical distance between 3 markers and their projected 2D coordinates, the position of the instrument tip can be estimated by fitting the physical geometry of the tool onto the projected lines Ob, Oc and Od (Fig. 3) [START_REF] Liu | A portable 3d vision coordinate measurement system using a light pen[END_REF].

Experimental Setup and Results

The proposed system was initially evaluated on five temporal bone phantoms (corresponding patient ages: 1-55 years). All specimens underwent a preoperative CT scan (Light speed, 64 detector rows, 0.6 x 0.6 x 0.3 mm 3 voxel size, General Electric Medical Systems, France). Six 1 mm fiducial markers, were attached around TM in a non-linear configuration with their combined centre coinciding with the target [START_REF] West | Fiducial point placement and the accuracy of point-based, rigid body registration[END_REF]. Real-time video was acquired using a microscope lens. Small movements were applied to the microscope in order to test the robustness of the system. The experimental setup and the augmented reality output are shown in Fig. 4. Processing speed of 12 frames per second (fps) was realized. The fiducial marker points in the reconstructed image were automatically detected and displayed on the screen and the user selected their corresponding fiducial points on the microscopic frame. Mean fiducial registration error (physical distance between estimated positions in microscopic and transformed reconstructed images) of 0.21 mm was observed.

During surgery, microscope will remain quasi-static. However to validate robustness of the system, combinations of translation, rotation and scaling with a speed of 0-10 mm/s were applied to the microscope. The system, evaluated at 30 second intervals, maintained synchronisation with a slow propagation error of 0.02 mm/min (Fig. 5a). Fiducial markers were used as reference points for evaluation. Template matching was used to automatically detect the fiducial points in the current frame. These were compared with the fiducial points in transformed reconstructed image to compute the tracking error.

The system was also evaluated on pre-chosen surgical target structures (incus and round window niche). TM of four temporal bone cadavers was removed and the above experiments were repeated. Similarly, a mean target registration error (computation similar to fiducial registration error) of 0.20 mm was observed with a slow propagation error of 0.05 mm/min (Fig. 5a).

For instrument pose estimation, pre-known displacements were applied in each axis and a total of 50 samples per displacement were recorded. Mean pose estimation errors of 0.20, 0.18 and 0.60 mm were observed in X, Y and Z axes respectively (Fig. 5b). The pose estimation in X and Y axes was better than in Z axis because any small deviation in instrument identification constitutes a relatively large deviation in the Z pose estimation.

Conclusion

An AR based robotic assistance system for transtympanic procedures was presented. A preoperative CT scan image of the middle ear cleft was combined with the real-time microscope video of TM using 6 fiducial markers as point correspondences in a semi-automatic RANSAC based homography framework.

The system is independent of marker placement technique and is capable of functioning with endoscopes and mono/stereo microscopes. Initial registration is the most crucial stage as any error introduced during this stage will propagate throughout the procedure. Mean registration error of 0.21 mm was observed.

To keep synchronisation, the relative microscope-target movements were then tracked using a SURF based robust feature matching framework. A microscopic propagation error was observed. Simultaneously, 3D pose of a needle instrument, upto 0.33 mm mean precision, was provided for assistance to the surgeon using a monovision based perspective framework. Additional geometric priors can be incorporated to compute pose of angled instruments. Initial experiments have shown promising results, achieving sub-millimetric precision, and opening new perspectives to the application of minmally invasive procedures in otology.

  (a) Global system workflow. (b) Initial registration. (c) Movement tracking. (d) Instrument identification.

  (a) Tympanic membrane. (b) Middle ear cleft. (c) 3D CT endoscopy. (d) Microscopic image. (e) CT MPR image. (f) Reconstruction image.

Fig. 2 :

 2 Fig. 2: Problem definition. Amalgamation of a reconstructed CT image (c) with the endoscopic video (a) may be used to visualize the middle ear cleft structures (b) without undergoing a TM flap procedure. However, similarity between them is low so fiducial markers are introduced which appear (d) grey in the microscopic image, (e) white in the CT MPR image and (f) as protrusions on the CT reconstructed image.

Fig. 3 :

 3 Fig. 3: Three collinear point framework for 3D pose estimation.

  (a) Experimental setup. (b) Augmented Reality Window.

Fig. 4 :

 4 Fig. 4: (a) Augmented reality system. (b) Real-time video from the microscope (left) and augmented reality window (right).

  (a) Registration and tracking accuracy. (b) 3D pose estimation accuracy.

Fig. 5 :

 5 Fig. 5: Experimental results. (a) Registration and tracking accuracy of the AR system evaluated at fiducial and surgical targets. (b) Displacement accuracy assessment of 3D pose estimation process (displayed statistics are for 50 samples).


Aurélien Cabarbaye

Titouan Verdu

Fabien Garcia

Michel Gorraz

Alexandre Bustico

Murat Bronz

Gautier Hattenberger

Design of a high performance MAV for atmospheric research

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Unmanned Aerial Vehicles (UAV) are now commonly used for scientific research and especially in atmospheric research [START_REF] Brandon | Development of an unmanned aerial vehicle for the measurement of turbulence in the atmospheric boundary layer[END_REF][START_REF] Bärfuss | New setup of the uas aladina for measuring boundary layer properties, atmospheric particles and solar radiation[END_REF] with a wide range of scale and mission. They have proved themselves to be cheaper and more agile than manned aircraft or balloons to probe the low-level atmospheric boundary layer. It is also important to be able to operate them in all sorts of weather conditions, including rain and strong winds [START_REF] Allan | Within limits, how are drones being used to improve weather predictions[END_REF]. However, flying a small UAVs at high speed, in harsh conditions and with a flight time as long as possible raises many technical challenges, especially to keep the overall system easy to operate.

One of the key measurement is the turbulence of the atmosphere which usually requires the use of 3D wind probes [START_REF] Calmer | 3D Wind Vector Measurements using a 5-hole Probe with Remotely Piloted Aircraft[END_REF]. These type of sensors are already commercially available, like the µADC from Aeroprobe, but are pretty expensive and requires cumbersome electronic boards. It is also important to increase the sampling frequency with the speed of the aircraft in order to maintain a large usable bandwidth [START_REF] Lawrence | High-resolution atmospheric sensing of multiple atmospheric variables using the datahawk small airborne measurement system[END_REF].

This article presents the design of a high performance small UAV designed for atmospheric turbulence measurements at high speed. After presenting the mission and the operational constraints in section 2, section 3 will describe the aircraft conceptual design and its innovative launching system, the next section will then focus on the on-board system integration and the design of a 3D wind probe. Finally preliminary flight and sensors tests will conclude the paper.

MISSION OVERVIEW AND CONSTRAINTS

The mission considered in this article is the measurement of the turbulence of the atmosphere with a small UAV. This is part of the NEPHELAE project, which involves atmospheric scientists as well as researchers in aerial robotics.

From the long experience of the authors in similar UAV operations [START_REF] Hattenberger | Flight tests for meteorological studies with MAV[END_REF][START_REF] Reymann | Adaptive sampling of cumulus clouds with UAVs[END_REF], it has been decided that a particular focus will be made on the take-off and landing procedures, as well as the capability of the plane to fly in strong wind conditions.

The UAV should be able to fulfil the following requirements:

• Lift a maximum of 800 grams of scientific payload in a roomy enough cargo bay.

• Cruise more than 2 hours to give the entire observation period coverage without having to land for recharging.

• Cruise at an average altitude of 3000 m above sea level where weather phenomena take place.

• Cruise at a speed of around 25 m/s which is fast enough to both counter the violent turbulences accompanying observed weather phenomena and to move from one measurement point to another in an acceptable time.

• The whole, fully equipped, aircraft should not weight more than 2.5kg to be easily transportable even in remote operating areas.

One of the main issue that is arising from these requirements is coming from the rather high cruise flight speed which prevents the use of off-the shelves airframes that would not be able to fly as fast (or as long, if overpowered) as expected, as it will be seen in section 3.

Figure 1 shows a typical location for operating UAVs without any proper ground infrastructure such as a runway. There is not even a flat space for belly landing. The UAV system must therefore presents VSTOL capability, which will be detailed in sections 3.2 and 3.3.

AIRFRAME DESIGN

Aircraft conception is an iterative process that covers the conceptual design, the definition of the launching and recovering systems and the manufacturing choices.

Aircraft conceptual design

The conceptual design consist in optimizing the aircraft characteristics assessing its performances on a typical mission.

The basic architecture of the aircraft is predefined:

• A pusher configuration is adopted in order not to interfere with the sensors accommodated in the nose of the aircraft (cf. section 4).

• An electric power-plant is selected for its simplicity in operation. Sufficient autonomy can be achieved when combined with lithium-ion batteries. The battery is a homemade pack of 18650 type cells of 3500mAh capacity each.

• The aircraft is likely to withstood high vertical wind gradients. Therefore, a flying wing configuration is adopted in order to reduce turbulence by minimizing the longitudinal distance between the aerodynamic center and the lifting surfaces.

• The electric systems in the fuselage are air cooled by a flow entering through a particle separator located in the nose.

One of the design case mission considered is presented in Figure 2 which is a vertical profile up to the top of the cloud formation.

Main optimization parameters are the chord and the span of the wing, the motor and the propeller. The software used (GENCAB) for this study is based on genetic algorithms, differential evolution and non-linear simplex (Nelder Mead algorithm). This hybridisation of global and local techniques makes the convergence of the overall algorithm quicker and also increases the robustness of the tool over a variety of problems, and in particular problems involving alphanumeric parameters. The optimization process leads to a high wing loading to reduce the zero lift drag produced by the wing wetted area. High wing loading presents the other advantage to be less turbulence sensitive. However its drawback is a higher stall speed which has been limited to 13 m/s for take-off and landing operations.

Moreover, the propulsion system maximum power have been chosen much higher than what is required in cruise in order to reduce the climb time to 3000m.

Launching systems

The relatively high stall speed imposes the use of a launching systems since a whole prototype has been lost during a hand launched test.

Several existing options have been considered:

• Cataplut launching seems to be the most popular solution. Tests have been performed with two off the selves systems. It results that either they do not transmit sufficient energy to the aircraft (SKYWALKERRC) or the acceleration is much too high (Gatewing) for the structure, the autopilot and the payload. Figure 3 shows the position of a dummy UAV along the time obtained during a test. The final velocity is about 17.5m.s -1 which Figure 3: catapult launching performance would be sufficient, but the acceleration is around 9.8g which is unacceptable. In addition, the maximum velocity is reached at less than one meter high which is dangerous for operators and does not provide sufficient clearance.

• Winch launching using a bungee and placing the plane on a ramp at the beginning achieves good repeatability. It presents much lower accelerations, however this solution requires planes with rather low stall speed so that the plane doesn't touch the ground after leaving the ramp. If it is not the case, a longer ramp can be used at the expense of a possible uncontrolled yaw moment due to the friction of the wing on the ramp.

• Car (manned or RC) top launching do not present the latter defect but requires an even longer clear flat area to be operated.

• Trebuchet launching or multicopter lifting are more unusual solutions because they require relatively bulky facilities.

• Several UAVs are rocket launched following the idea of the Zero-length launch of manned aircraft developed in the 1950's. The main issues with this method is the fire hazard represented by the solid fuel rocket.

Therefore, a specific takeoff system based on a water rocket has been designed. Its main advantages are: limited acceleration, high height clearance, very high reachable airspeed, very compact facilities, only little hazard compare to the solutions previously presented. In addition, the launch pad consists in a tube inserted in the water rocket that guide it at the very beginning of the launch and reduce the waste of water acting like a piston. The maximal acceleration can be limited by reducing the nozzle diameter, while the internal volume of the water rocket is computed in order to store the required amount of energy. The acceleration has been limited here to 4g. The UAV is released at 15 m/s and 8 m high with an initial pressure of only 9 bars. The acceleration is thus much more gentle than with a catapult as it can be noticed in Figure 4.

However, its integration is more restrictive than for a classical Jet-Assisted Take-Off (JATO) system that is jettisoned in flight like a drop tank. As water rocket has indeed a much lower specific energy than solid fuel rocket, it cannot be fixed underneath the aircraft and must be placed at the rear of the drone like the first stage of a rocket. In this position, the CG is very far back. An additional lifting surface must therefore be added at the very back of the water rocket to restore the natural longitudinal stability. The additional lifting surface area is computed in order to move the aerodynamic center behind the center of gravity of the set with a static margin of 5% of the UAV aerodynamic chord. The UAV and the rocket are connected thanks to a cone shape coupler made of a male plug mounted on the motor shaft and a female plug put on top of the rocket. A parachute is added to the assembly through a cable ensuring the separation of the two components on time.

Recovering systems

Several solutions can be implemented to obtain a short or vertical landing. The most popular solution consists in accommodating a ballistic parachute in the airframe. Despite the fact that very reliable systems are available on the market for the considered aircraft size, the additional weight is too high not to be detrimental.

Arresting hook is another solution but it still requires a piece of runway and constitute another overweight.

Therefore, a vertical arresting net has been selected since it can be used on almost every cleared ground and do not need any additional on-board system.

Manufacturing process

The manufacturing process is custom to deal with the unusual constraints of the chosen UAV operation mode. In addition, the manufacturing must be easy, fast, low cost in order to produce the whole fleet on time, on cost and with the re-quired quality. To do so, several solutions have been adopted as shown in Figure 6. The high aspect ratio combined with the high wing loading and the high expected turbulence impose the use of a spanwise carbon spar. In order to gain weight, the back half of the wing is built from an extruded polystyrene foam block with a hot wire CNC. This material has a high stiffness suitable for the sharp trailing edge. The front half is built with the same tool but from a block of expended polypropylene which offers a very good shape memory that is of great interest when the UAV impacts the arresting net. The joint between the two parts is made of a mixture of resin and glass micro-sphere. The obtained wing is covered with a high-density polyethylene film that protects it against abrasion and improves its polishing.

The planned number of UAVs is too low to justify the purchase of fuselage moulds. Therefore an innovative manufacturing process has been developed. The external skin of the fuselage is 3D printed in PET. The PET is chosen because of its very good inter-layer cohesion and its very low wrapping despite its higher density compared with ABS. However, the fuselage is built with the minimum possible thickness that ensures good printing result. The strength of the fuselage is obtained by laying carbon fibers on the inside surface of the printed part. The resulting fuselage presents the strength of a carbon molded part at the cost of a 3D printed part.

The water rocket is built by splicing two 2L soda bottle back to back. the maximum tolerable pressure is around the expected 9 bars. Therefore, a carbon fiber fabric is wrapped around the bottles assembly. Test have been performed up to 12 bars without showing any weakness. The wing is made of a carbon spar reinforced extruded polypropylene foam shape to withstood hard landings. The general layout of the water rocket booster is shown in Figure 7.

SYSTEM INTEGRATION

Autopilot

The open-source autopilot system used to equipped the UAVs is called Paparazzi [START_REF] Hattenberger | Using the Paparazzi UAV System for Scientific Research[END_REF]. The previous hardware boards used during past project was reaching end-of-life since some components were declared obsolete by the manufactures, so it was decided the design a new version.

The general guidelines were to keep the board small, eventually with the same footprint than the previous generation for easier upgrade, improve the connectivity to support more sensors and use latest available integrated components such as MCU and IMU. The main characteristics of this new board, named Tawaki, are listed in Table 1 The Figure 8 shows the top and bottom 3D view of the final design of the Tawaki autopilot. A first prototype have already been assembled and successfully tested so far.

Communication system

The communication system on this UAV needs to provide at least a 10 km range and accommodate for several aircraft. Table 2: General characteristics for the P2400 modem.

In order to estimate the communication range, a modem has been installed on board a UAV flying in circles at 150 m altitude while a car was driving away with another modem attached to a computer running a test program. To the day of writing this article, tests have been done with a wireless rate of 172.8 kbps resulting in a maximum range of 4.3km. This range is not sufficient for the planned mission, nevertheless the sensitivity of the modem (and hence its maximum communication range) increases when one configures a lower wireless rate allowing for a range increase. Through simple computation with the Friis free-space transmission formula, the feasible range have been evaluated for a 55 kbps wireless rate to be around 7.7 km. Furthermore, we plan to use a directive antenna which would give us a further boost in range up to 14 km which should provide enough range even counting the attenuation due to clouds.

Turbulence probe

In order to measure the turbulence of the atmosphere, a commonly used type of sensors are the multi-hole probes. Their primary usage are to measure the angle of attack and the sideslip angle, which doesn't require a very high sampling frequency, around 10 Hz. However, the dynamic of the turbulence requires a larger bandwidth with a minimum of 100 Hz.

For this project, a new integrated probe using fast differential pressure sensor is under development. The main characteristics are presented in Table 3 Since it integrates a micro-controller and a SD card directly inside the probe, the recording and processing at high frequency can be done with minimum latency and independently of the main autopilot. A serial UART connection allows to forward the filtered and pre-processed data to autopilot, which will eventually forward the messages to the ground for real-time monitoring.

In addition to this, IMU data can be recorded for angle correction, or processed in real-time using attitude filter. The serial connection can also be used to receive external data such as GPS position, speed and heading for a better reconstruction of the local 3D wind field.

The Figure 9 shows the preliminary design of this new sensor. The expected sampling frequency is 500 Hz for the differential pressure and 50 Hz for the absolute pressure.

Other sensors

In addition to the 3D wind probe, other meteorological sensors can be embedded inside the plane:

• Temperature and humidity, mostly used to build vertical profiles of the atmosphere the same way weather balloons would do.

• Cloud sensor [START_REF] Harrison | Volcanic ash detection[END_REF] used to determine some of the characteristics of the water droplets and the liquid water content (LWC) of the atmosphere.

• Radiation sensors used to measure the albedo.

FLIGHT TESTS AND EXPERIMENTAL RESULTS

Flight patterns simulations

In order to efficiently sample the clouds, dedicated flight patterns have been defined based on the feedback from the scientists that will post-process them. The result is a set of predefined adaptive shapes that allows focusing the acquisition on certain parts of the cloud such as the border, the base or the core as shown on Figure 10. The real-time feedback from the sensors reduces the required time to sample a particular region corresponding to the physical process being studied, thus providing a better set of data to reconstruct the complete evolution of the cloud. The other benefit of this expert-based approach is that the computational load is limited compared to previous work on this topic [START_REF] Reymann | Adaptive sampling of cumulus clouds with UAVs[END_REF] and can be performed in real-time by the on-board autopilot without heavy computations on the ground.

These strategy have been tested in simulations with realistic aircraft and cloud models [START_REF] Verdu | Flight patterns for clouds exploration with a fleet of uavs[END_REF]. The Figure 11

Flight tests

Preliminary test flights have been performed with the optimized planes. At the very beginning, the issues during takeoff using conventional techniques, leading to a loss of the plane, conduct to the design of the rocket-based launching procedure.

The second second session of flight tests have revealed an issue linked to the structural vibrations of the aircraft. Those vibrations, which disturb the attitude (AHRS) estimation filter, should be reduced prior to more flights. Two solutions have been implemented, which are a mechanical damping of the autopilot and a better filtering of the accelerometer data.

Five-Hole Probe Calibration (Prototype)

At the time of writing this draft paper, the final turbulence probe was not ready. In preparation for it, the calibration procedure has been done with prototype electronics, and an off-the-shelf 5-hole probe from AeroProbe. The prototype electronics consists of three differential pressure sensors (SDP31), an absolute pressure sensor, and STM32L4 microcontroller.

The calibration process has been done inside a low speed wind-tunnel facility, where the test section dimensions are 50cm by 50cm. The static ports on the probe are located at the very tip of the tube that has a spherical shape.

During the calibration process, it is important to keep the tip of the probe at the same location (fluid stream tube) in order to avoid possible uncertainties that may arise from nonuniform distribution of the speed. Therefore, a custom calibration bench has been designed and build that holds the tip of the probe in its center during the pitch and yaw movements. Figure 12 shows the calibration mechanism inside wind tunnel's test section. The bench has been designed to move from -20 deg, up to +20 deg both in angle of attack and side slip axes. During calibration, first the airspeed is fixed, and then for each angle of attack, side slip angle has been changed. Once all angle of attack and side slip angles are covered, the wind tunnel speed has been changed and the process is repeated for different airspeed (8m/s, 11m/s, and 14m/s).

A linear regression on the calibration data has been ap- plied in order to estimate the angle of attack, side slip angle and the dynamic pressure by using the three on-board differential pressure sensors as

  α β q   = A   P 1 /P 3 P 2 /P 3 P 3   +   c 1 c 2 c 3   ,
where A is a [3 × 3] coefficient matrix (obtained from the linear regression), and c 1-3 are the offset values.

As the angle of attack and side slip angles should be estimated independent of the airspeed, the linear regression features has been selected in a normalized way being P 1 /P 3 and P 2 /P 3 . It can be understood that P 3 is the differential pressure sensor that is connected to the total pressure and static pressure holes.

Figure 13 shows the result of the predicted angle of attack, side slip angle and the airspeed with respect to ground truth that comes from bench mechanism servo angles and windtunnel airspeed setpoint. The actual airspeed has been calculated by predicted dynamic pressure and the measured onboard absolute pressure and temperature values.

CONCLUSION AND ONGOING WORK

This paper has presented the overall design of a small UAV optimized for meteorological research, such as cloud and turbulence analysis. This includes the conception of the plane itself, taking into account the requirements of the atmosphere researchers, but also the adaptation of the manufacturing process to make the structure easily reproducible and adaptable and the launching system adapted to the specific performances of the aircraft. In addition, a special attention to the system integration has led to the conception of an integrated 3D wind probe based on 5-hole measurements, a new version of the main autopilot board and the integration of long range communication devices.

Future work will focus on the plane performances evaluation in real-flight, completing the final design of the 3D wind probe and manufacturing it in order to test it in flight.

Figure 1 :

 1 Figure 1: Base camp for operations of MAV in harsh environment.

Figure 2 :

 2 Figure 2: UAV typical mission.

Figure 4 :

 4 Figure 4: Water rocket performances

Figure 5 :

 5 Figure 5: Water rocket launch phases.

Figure 6 :

 6 Figure 6: UAV construction.

Figure 7 :

 7 Figure 7: Rocket booster construction.

Figure 8 :

 8 Figure 8: 3D view of the final design of the Tawaki autopilot

Figure 9 :

 9 Figure 9: Inner view of the integrated 3D wind probe

Figure 10 :

 10 Figure 10: The different parts of the cloud that will require a specific focus during the mission.

  is showing the flight tracks of three UAVs where two of them are measuring the vertical wind intensity at two different altitudes to estimate the airflow inside the cloud, while the last one is covering the border in 3D to estimate the volume. The turnarounds at the edge of the cloud are triggered based on the liquid water content estimation in real-time.

Figure 11 :

 11 Figure 11: Simulation of 3 UAVs sampling 2 horizontal planes and the last one covering the outer 3D border.

Figure 12 :

 12 Figure 12: Calibration bench located inside the test section of ENAC's low speed wind-tunnel.

Figure 13 :

 13 Figure 13: Comparison of predicted versus ground truth values of angle of attack, side-slip angle and the airspeed.

Table 1 :

 1 . General characteristics for the Tawaki autopilot board.

	MCU	STM32F7
	IMU	ICM20600 (accel, gyro) + LIS3MDL (mag)
	Baro	BMP3
	Serial	3 UARTS, I2C (5V + 3.3V), SPI
	Servo	8 PWM/DShot output (+ ESC telemetry)
	RC	2 inputs: PPM, SBUS, Spektrum
	AUX	8 multi purpose auxiliary pins (ADC, timers, UART, flow control, GPIO, ...)
	Logger	SD card slot
	USB	DFU flash, mass storage, serial over USB
	Power	6V to 17V input (2-4S LiPo) 3.3V and 5V, 4A output
	Weight	12 grams

Table 3 :

 3 . Main features of the 3D wind probe

	MCU	STM32F7
	Differential pressure	3 x SDP31
	Absolute pressure	LPS33HW
	IMU	ICM20600 (accel, gyro)
	Logger	SD card slot
	Data	UART, USB
	Powering	5V
	dimensions	22 mm, L 110 mm min

ACKNOWLEDGMENTS

The NEPHELAE research project is founded by the French National Research Agency (ANR).

The authors would like to thank the researchers from the LAAS-CNRS and Meteo-France involved in this project.