Damien Magoni

Network Topology Analysis and Internet Modelling with nem

Keywords:

The design of network protocols is greatly accelerated by the use of simulators particularly when studying a protocol's behavior in a large internetwork topology. However, the accuracy of the simulation results is heavily affected by the input network topology. As taking a real map as an input is not always feasible, artificially created topologies are often used. There exist many network topology generators available but recent discoveries on the Internet topology have made most of them obsolete when it comes to modelling a typical part of the Internet. This paper presents a free open source software designed for network topology analysis and modelling called network manipulator (nem). It is capable of creating realistic Internet-like topologies and it can check proof them on the fly by a thorough topology analysis. The paper especially focuses on the architecture and the capabilities of the software.

I. INTRODUCTION

The use of simulation for protocol design has increased the need for accurate network topologies. These topologies used as input for simulations must be realistic with respect to the OSI level of the protocol being evaluated. In a simulation, a layer 2 protocol for example should be evaluated on a topology similar to its real topology (e.g. a bus, a ring, etc.) 1 . Thus, as many of the layer 3 protocols are created in order to be deployed in the Internet (e.g. protocols for routing, multicasting, mobility, etc.), their simulation should be performed upon Internet-like topologies.

Early network topology generators were rather basic and not accurate in internetwork topology modelling [START_REF] Waxman | Routing of multipoint connections[END_REF], [START_REF] Zegura | A quantitative comparison of graph-based models for internetworks[END_REF]. However, since the discovery of the power laws in the collected maps of the Internet, many topology generators have been created to fulfill this task. They use different method to reach their goals and the graphs that they can generate should still be checked against the highest number of criteria to ensure their accuracy in modelling the topological properties found in the Internet maps. We introduce a free open source software that not only generates Internet-like graphs (by several methods) but also analyzes their topological properties for an immediate control of their realism. The paper contains five sections. We briefly review other similar tools in section 2. The capabilities of our software are detailed in section 3. The architecture of nem is presented in section 4. Section 5 describes some important parts of the nem implementation. Finally, section 6 contains results showing the accuracy of the nem graphs vs the collected Internet maps and a performance comparison of nem vs other power law topology generators.

II. PREVIOUS WORK

Among the topology generators, one of the earliest and most famous model was designed by Waxman in 1988 [START_REF] Waxman | Routing of multipoint connections[END_REF]. This kind of model is usually called flat topology model. The nodes are randomly placed on an Euclidean plane irrespective of any hierarchy order among them. This model was later replaced by hierarchical topology models such as the Tiers [START_REF] Doar | A better model for generating test networks[END_REF] and the Transit-Stub [START_REF] Zegura | A quantitative comparison of graph-based models for internetworks[END_REF] models. These models try to enforce a multi-level hierarchy that can be found in the Internet (e.g. host-router-AS). The discovery of power laws in the Internet [START_REF] Faloutsos | On power-law relationships of the internet topology[END_REF]- [START_REF] Magoni | Analysis of the autonomous system network topology[END_REF] has brought the creation of a new kind of topology model. We call it the power law topology model because, as the name suggests, this model tries to reproduce power law distributions in order to generate Internet-like graphs. nem, like all the tested generators presented in section 6, belongs to this topology model category and it can produce router level graphs like all the other generators or AS level graphs like the Inet2 generator [START_REF] Jin | Inet: Internet topology generator[END_REF]. A recent study by Radoslavov et al. [START_REF] Radoslavov | On characterizing network topologies and analyzing their impact on protocol design[END_REF] provides a detailed analysis of topology characterization of real, generated and canonical topologies. Generated topologies include Tiers [START_REF] Doar | A better model for generating test networks[END_REF], Transit-Stub [START_REF] Zegura | A quantitative comparison of graph-based models for internetworks[END_REF] and Waxman [START_REF] Waxman | Routing of multipoint connections[END_REF] graphs. Several multicast scenarios are then played on all these topologies to examine the influence of the network topology on protocol simulation. Unfortunately no power law topology model generators are evaluated in this work. We partly fill this gap by looking at the values of some topology properties measured in graphs produced by recent power law topology models such as Model A [START_REF] Aiello | A random graph model for massive graphs[END_REF], BRITE [START_REF] Medina | On the origin of power laws in internet topologies[END_REF], Extended Scale Free (ESF) [START_REF] Albert | Topology of evolving networks: local events and universality[END_REF], Inet2 [START_REF] Jin | Inet: Internet topology generator[END_REF], PLOD [START_REF] Palmer | Generating network topologies that obey power laws[END_REF] and nem although we do not run multicast scenarios on the produced topologies. This is left for future work. Concerning the analysis of the topological properties of graphs, few network generation tools provide an in-depth topology analysis of the produced graphs and the researcher must usually use another software than the one which generated the graphs in order to obtain detailed topological information. As the aim of the topology models and their corresponding generators is the production of realistic network topologies (especially Internet-like) for feeding network simulators such as ns-2 [START_REF] Project | network simulator (ns-2)[END_REF], OPNET [START_REF]network simulator[END_REF] and GloMoSim [START_REF] Glomosim | Parallel Computing Laboratory[END_REF], it is interesting for the researcher to simultaneously get the generated graph and its topology analysis. The latter should be immediately comparable to an Internet map topology analysis in order to check proof the generated graph realism. nem provides this functionality thus making the validation of a generated graph easier. Internet reference maps include maps created by Govindan et al. at the Information Sciences Institute by the Mercator software [START_REF] Govindan | Heuristics for internet map discovery[END_REF] and by Burch et al. at Lucent labs [START_REF] Burch | Mapping the internet[END_REF].

III. CAPABILITIES

In this section we present the capabilities of our software. nem can do three fundamental tasks:

1) Convert network files from one format to another.

2) Analyze the topology of networks.

3) Generate random networks. nem has no user interface but works with an ASCII batch file called a process file. Each line contains one or more instructions to be executed on a network. The detailed syntax of a process file is described in the nem manual.

A. Conversion and Interaction with other Software

nem can handle a variety of network file formats generated by other software. nem is able to load measured Internet maps such as the ones produced by the famous inter-domain router route-views (Autonomous System level maps) or by the mapping software Mercator (router level maps). These maps can be converted into the format of the ns-2 or OPNET simulators to serve as input topologies although they are usually too big. Our map sampling algorithm enables the generation of graphs from these maps that keep their fundamental properties while having a much smaller size (see section 6). The graphs can also be viewed in 3 dimensions by choosing the HypViewer output format and by using the HypViewer software.

B. Topology Analysis

Running an analysis on a network produces a file that has the same prefix than the network filename and has the extension * .analysis. The analysis concerns only the topology of the network (seen as a graph). Among the most important topology properties whose values are calculated by nem, we have the degree, the tree, the mesh, the distance, the number of shortest paths and the bicomponent related properties.

Depending on the characteristics of the graph (directed or not, connected or not, etc.) some properties are not calculated (see section 5 for details on the course of an analysis). Also depending on the kind of graph (random graph, Internet-like graph, Internet map, etc.) some properties are not significant (e.g. power law related properties only make sense in Internetlike graphs and maps). However, many calculated properties are valid and interesting whatever the kind of topology of the analyzed network.

Furthermore, the Dijkstra's outer loop can be parallelized in order to analyze huge networks. In the analysis of a graph, this algorithm is particularly time consuming. The all-pairs Dijkstra algorithm calculates the shortest path from any one node to another. It is the Dijkstra algorithm [START_REF] Aho | Data Structures and Algorithms[END_REF] implemented with a binary heap which runs in O(nlog 2 n) and which is run for every node in the graph. Thus the whole algorithm runs in O(n 2 log 2 n). In nem, the outer loop of the algorithm can be cut into slices to enable parallelization. At the end of the slice analysis, a partial analysis dump file of the slice is produced. The latter use of the merge operator on all the partial analysis dump files produces a file which contains the final analysis of the properties depending on the Dijkstra's algorithm (distance properties).

C. Graph Generation

The generation of a network topology is done through the use of a specification file. It is an ASCII file containing the values of the parameters needed to create a network. This file must have the extension * .specif. Common parameters include the number of nodes, the number of edges, and the name of the network output file. The user has the choice between several generation methods given in table I. Only the map sampling method has been created by us and we describe it below. The other methods available in the nem generation module have been defined in research papers. However we found convenient to implement them in nem (i.e. because they were not implemented elsewhere or their implementations were not available) by carefully following the algorithms described in the corresponding papers. Some parameters are common to all generation methods while others only make sense and work for a given generation method.

Detailed explanations for the methods, their parameters and their effects on network generation are provided in their respective papers whose references are given in table I. It's worth noticing that the user can generate a random graph if needed instead of an Internet-like graph. For doing this we provide the modified Waxman method. It works like the regular Waxman [START_REF] Waxman | Routing of multipoint connections[END_REF] algorithm but connectivity is ensured by repeatedly linking connected components together until there is only one left. Another possibility is to use our map sampling algorithm on a topology that is not a specific Internet-like topology. The resulting graph will not a priori have an Internetlike layout.

We describe here our own graph generation method. It makes use of an algorithm that randomly extract a subgraph of a real Internet map. As we focus on the router level of the Internet, we use router level Internet maps as input graphs (but it works in the same way for the AS level). By extracting nodes in a mostly random fashion, we ensure the creation of a wide diversity of subgraphs. Fig. 1 shows how the tree creation algorithm works. Fig. 1 (a) shows a small part of the map being sampled. In the first step (b), node 1 is selected randomly among all the nodes of the map. Nodes 2 and 5 are then added to a candidate list (see below) with candidate link a and b respectively. In step (c), node 5 is selected randomly among the nodes of the candidate list. Having candidate link b, it is connected to the tree (composed only of node 1 for the moment) by this link. Nodes 4, 6 and 7 are added to the candidate list with candidate link e, f and g respectively. Node 2 is already in the candidate list so it is not added to it. As node 2 is connected to 5 by link d, a random roll is performed to check if its candidate link changes to d (it is currently a). The roll fails and the candidate link is not changed. In step (d), node 2 is selected randomly among the nodes of the candidate list. Having candidate link a, it is connected to the tree (composed of nodes 1 and 5 for the moment) by this link. Node 3 is added to the candidate list with candidate link c. As node 5 is already selected, link d is added to the redundancy link list for later use (after the tree construction). In step (e), node 6 is selected randomly among the nodes of the candidate list. Having candidate link f, it is connected to the tree (composed of nodes 1, 2 and 5 for the moment) by this link. Node 8 is added to the candidate list with candidate link i. Node 7 is already in the candidate list so it is not added to it. As node 6 is connected to 7 by Fig. 1. tree creation process Fig. 2. nem architecture link h, a random roll is performed to check if its candidate link changes to h (it is currently g). The roll succeeds and the candidate link is changed to h. In the last step shown (f), node 7 is selected randomly among the nodes of the candidate list. Having candidate link h, it is connected to the tree (composed of nodes 1, 2, 5 and 6) by this link. As node 5 is already selected, link g is added to the redundancy link list for use after the tree construction. This algorithm enables us to extract a tree having a number of nodes equal to the one of our future graph and to create a supply of redundancy links. The complete and detailed pseudo code of our map sampling method can be found in [START_REF] Magoni | Internet topology modeler based on map sampling[END_REF]. The amount of time needed to generate nem graphs with our map sampling method is in the order of seconds.

IV. ARCHITECTURE

The architecture of nem is shown in fig. 2. It is a schematic view of the modules constituting nem and how they interact with each other. This picture does not represent how nem is implemented.

The manipulator module takes as argument a process file. It is a batch file that contains a description of the actions to be carried out by nem. An action consists of one input part and one or more output parts. In the input part, nem can either load a network file (already produced by nem or another Fig. 3. UML activity diagram of the network topology analysis generator) or generate a network from a set of parameters given in a specification file. The input converter module is controlled by the network filename extension (if known by nem). If the extension is .specif, the generator module is called, otherwise the loader module is called by the converter if the network file format is recognized (see table 1 for a list of the supported input file formats). In the output part, nem can analyze the graph and/or convert it in another format available for output.

V. IMPLEMENTATION

nem has been implemented in C++ and contains around 9000 lines of code. The C++ language has been chosen for the speed of the binary code and the efficiency of the source code management due to the object oriented nature of this programming language.

The function calls in the analysis module of nem is depicted as an UML activity diagram in fig. 3. The degrees' phase calculates the degree distribution and the related power law values. The connectivity phase determines if the graph is directed or not. This depends on the network input file (most of them are undirected though). If the graph is undirected, the Baase algorithm [START_REF] Baase | Computer Algorithms[END_REF] is used to determine if the graph is connected otherwise the Tarjan algorithm [START_REF] Aho | Data Structures and Algorithms[END_REF] is used (although it also works for undirected graphs). If the graph is not connected, the number and the distribution size of the connected components are calculated and the analysis ends. Otherwise, if the graph is directed, the Dijkstra algorithm [START_REF] Aho | Data Structures and Algorithms[END_REF] is run upon it in order to determine the distance property values such as average path length, diameter, etc. If it is undirected and not a unique tree, the number of cut points is determined and the Hopcroft algorithm [START_REF] Aho | Data Structures and Algorithms[END_REF] is performed in order to obtain information on the biconnected components of the graph (i.e. number, distribution size, number of vertices on a cycle, number of vertices on a bridge, etc.).

The map sampling algorithm of the generator module of nem is described in the previous section and the explanations concerning its implementation can be found in our previous work [START_REF] Magoni | Internet topology modeler based on map sampling[END_REF].

VI. PERFORMANCE EVALUATION

In this section we present the performances of the graph generation module of nem. When we speak of nem graphs in this section, we speak about graphs generated by using our map sampling method (see table 1). It is important to explain this because nem can generate graphs by using the Extended Scale Free (ESF) method for example and in this case the graphs would be called ESF graphs (even though nem has produced them). We first look at the accuracy of the nem graphs in comparison to Internet maps and then we compare graphs generated by using several power law generation methods including ours.

A. Accuracy with Respect to Internet Maps

In order to evaluate the accuracy of the nem graphs in comparison to Internet maps, we have generated 20 graphs of sizes ranging from 125 nodes to 4000 nodes using the Mercator Internet map of 1999. The average analysis of the 20 graphs has been produced by nem just after their creation by the generation module. The figures contain the plots of these values as well as the minimum and maximum values observed in the 20 graphs (by looking in the individual analysis files). Most properties present in the Internet map are accurately reproduced in the nem graphs. This was the objective of our sampling method (i.e. the graphs are extracted from the map). In particular, we have looked at the average path length distributions of a sample of the extracted graphs and they all have a Gaussian shape like the one found in the Internet map (overlooking the small oscillations due to the sampling scale). This means that we can consider the average values of these distributions as being interesting distance indicators. Fig. 4 clearly shows that it is possible to extract graphs that have the same average path length than a 200 times bigger map (shown by the dotted line). All graphs of size equal or above 1000 nodes up to 4000 nodes have average distances between 8 and 11 which are typical values measured in Internet router level maps [START_REF] Govindan | Heuristics for internet map discovery[END_REF], [START_REF] Hoerdt | Completeness of the internet core topology collected by a fast mapping software[END_REF]. For graphs of less than 1000 nodes it is hardly possible to obtain satisfying results.

It is worth noticing that our filtering function presented in [START_REF] Magoni | Internet topology modeler based on map sampling[END_REF] and which avoids extracting very high degree nodes, is necessary to achieve these results. The basic extraction is not enough to obtain accurate average path lengths at these size levels (a few thousand nodes). As we observed similar results for the node eccentricity distributions of the nem graphs, we can claim that the distance properties of the graphs extracted by our map sampling algorithm are satisfactory for graphs of sizes 1000 up to 4000, but we also found that these graphs do comply with the power laws found in the Internet map. As an example, fig. 5 shows the absolute correlation coefficients (ACC) for the rank and degree power laws [START_REF] Faloutsos | On power-law relationships of the internet topology[END_REF] of the graphs extracted by the map sampling algorithm. All extracted graphs (from size 125 to 4000) comply with these laws for any of the generated sizes with values above the 0.95 commonly used threshold [START_REF] Faloutsos | On power-law relationships of the internet topology[END_REF]. Also we found that all extracted graphs comply Fig. 4. average path length of the graphs Fig. 5. rank and degree ACC of the graphs with the tree size and tree rank power laws described in [START_REF] Magoni | Analysis of the autonomous system network topology[END_REF] when the graph size is equal or above 1000 nodes and equal or below 4000 nodes as shown in fig. 6.

To conclude, we showed that the graphs generated by nem comply with the rank and degree power laws found by Faloutsos et al. [START_REF] Faloutsos | On power-law relationships of the internet topology[END_REF] for all evaluated sizes and comply with the tree power laws found in our prior work [START_REF] Magoni | Analysis of the autonomous system network topology[END_REF] for sizes comprised between 1000 and 4000 nodes. Concerning the distance properties (average path length and average eccentricity), we saw that accurate results can be achieved for evaluated sizes ranging from 1000 to 4000 nodes. These graph sizes should be appropriate for most protocol simulations (i.e. neither too small nor too big).

B. Comparison Between nem and Other Generators

In order to evaluate the accuracy of the nem graphs in comparison to other power law topology generators, we have generated graphs of sizes ranging from 500 nodes to 8000 nodes using the ESF, Model A, Inet2, PLOD, BRITE and nem. Inet2 uses an Autonomous System map of 3025 nodes collected in 1999 as a seed for its graph generation method and thus is not able to produce graphs containing less than 3025 nodes. That is why there are only 2 points for Inet2 plots in the figures. As the process of generating graphs involves random selection (and thus random rolls), we have used a sequential scenario of simulation [START_REF] Law | Simulation Modelling and Analysis[END_REF] to produce the results shown in Fig. 6. tree rank and tree size ACC of the graphs this section. We have used the Mersenne Twister code [START_REF] Matsumoto | Mersenne twister: a 623dimensionally equidistributed uniform pseudo-random number generator[END_REF] for producing the random numbers needed in our simulation. As the random rolls are the only source of randomness in our simulation, we can reasonably assume that the simulation output data obey the central limit theorem. We have performed a terminating simulation where each run consists in generating a graph and determining its topological properties (i.e. one run is the time horizon). We have performed a sequential checkpoint each time after having created one graph. Convergence is reached between 8 to 21 graphs depending on the generation method and the topological property. All the simulation results have been obtained assuming a confidence level of 0.95 with a relative statistical error threshold of 5% for all measured metrics. We show the plots of the resulting values for the main topological properties concerning degree, distance, mesh and trees.

1) Degree Properties: Fig. 7 shows the degree absolute correlation coefficient (ACC) of the graphs. This indicator measures the compliance of the graph degree distribution to the degree power law (Faloutsos et al. [START_REF] Faloutsos | On power-law relationships of the internet topology[END_REF] power law number 2). Only BRITE does not comply with this power law. BRITE uses a parameter m that indicates how new nodes connect to the graph during its construction. We have examined some degree distribution samples of BRITE graphs with m = 2. They contain a few degree 1 nodes that cause these graphs to not comply with power law 2. These few nodes of degree 1 are encountered in all the BRITE graphs with m = 2 whatever their size. We have verified that BRITE, unlike the other generators, generates graphs having a degree distribution that starts at m instead of one (when m > 1). A few outliers of degree less than m cause BRITE graphs to fail compliance with power law 2.

Fig. 8 shows the degree exponent of the graphs. This indicator measures the slope of the line given by the least squares fitting of the log-log plot of the node degree distribution. Like the average concisely describes a Gaussian distribution, the exponent concisely describes a power law distribution. As BRITE graphs do not comply with the degree power law, no exponent can be calculated and thus there is no plot for BRITE graphs in fig. 8. The values of most graphs are around the threshold but only ESF and nem graphs reach the target 2) Distance Properties: Fig. 9 shows the average path length of the graphs. Only nem graphs get close to the average path length measured in Internet maps. PLOD graphs also manage to get above the minimum threshold. All other graphs including ESF graphs clearly do not accurately model the average path length. It's worth noticing that this property is very important because the hop count is one of the only few values that can be measured and exploited by IP routing and transport protocols. Fig. 10 shows the average eccentricity of the graphs. All graphs have difficulty reaching the target value of the indicator. However PLOD, nem and Model A manage to be around or above the minimum threshold. The other graphs are quite below it and thus their nodes are not spread as they could be in a real Internet map.

3) Mesh Properties: The definition given in [START_REF] Magoni | Analysis of the autonomous system network topology[END_REF] says that a node belonging to a cycle (i.e. a closed path of disjoint nodes) or lying on a path connecting two cycles is called an in-mesh node. The mesh is the set of all the in-mesh nodes of a graph. The mesh is important because this is the area where can be found all the redundant edges of the graph. Accurately modelling the mesh has a crucial impact on multiple shortest path or alternate path simulations. Fig. 11 shows the mesh degree exponent of the graphs. PLOD and Model A are out of the threshold area and thus do not match this indicator. Fig. 12 shows the mesh average path length of the graphs. nem graphs are close to the target value especially when the graph size increases. All other graphs are below the minimum threshold of this indicator (except for the 8000-node ESF graphs) and thus do not comply with this indicator.

4) Tree Properties:

We have also studied the trees in the graphs. The trees in an Internet map follow the tree size power law described in [START_REF] Magoni | Analysis of the autonomous system network topology[END_REF]. In the previous section we defined the mesh. The forest is simply the set of nodes of the graph that do not belong to the mesh. These nodes are located in trees and the union of these trees form the forest. The trees are connected to the mesh by special nodes called roots. We consider the roots as nodes belonging to the mesh. As trees are sensitive areas for link failure, it is important to model them with accuracy. We have not shown the tree rank ACC of the graphs because all graphs have values above the 0.95 threshold except BRITE graphs. Actually BRITE graphs do not comply with tree rank and tree size power laws because nearly 100% of their nodes are in the mesh. Thus there is not enough trees Fig. 11. mesh degree exponent of the graphs Fig. 12. mesh average path length of the graphs in order to compute tree rank and tree size least squares fitting. Fig. 13 shows the tree size ACC of the graphs (except BRITE graphs for reasons previously explained). Model A graphs do not match the threshold of this indicator and PLOD is just close enough. Other generators have good results.

The careful study of all these indicators have shown us that only the topological properties of nem graphs do comply with all of them (or at least reach their thresholds) as shown in table II. A 'yes' entry means that the graphs produced by the generator comply with the indicator, a 'thre' entry means that they are close from the indicator's threshold and a 'no' entry means that they do not comply with the indicator.

VII. CONCLUSION

Accurate network topology generation is an important research field for protocol simulation. In this paper we have presented a software designed for network topology analysis and modelling. It is particularly suited for Internet-like topology modelling. We have detailed its capabilities, architecture and implementation and we have evaluated its performances. We have shown that nem is appropriate for analyzing graphs and maps obtained from various sources and that nem is accurate when modelling Internet-like graphs. We hope that Fig. 13. tree size ACC of the graphs this software will be a helpful tool for the network research community. nem is open source and available on the Web at http://www-r2.u-strasbg.fr/nem. For future work, we are planning to develop a graphical user interface for nem and to increase the number of topological parameters analyzed in the graphs (such as the clustering coefficient) We have also already added several features such as the ability to produce an overlay map from one router level map and one Autonomous System level map and the ability to load maps produced by our network cartographer software.

Fig. 7 .

 7 Fig. 7. degree ACC of the graphs

Fig. 9 .

 9 Fig. 9. average path length of the graphs

This is an extended version of an article presented at the 10th IEEE/ACM International Symposium on Modeling, Analysis and Simulation of Computer and Telecommunication Systems in


RADARS ON THE AIM MISSION A Herique
email: alain.herique@obs.ujf-grenoble.fr

V Ciarletti

Aim Team

A DIRECT OBSERVATION OF THE ASTEROID'S STRUCTURE FROM DEEP INTERIOR TO REGOLITH: TWO

Introduction:

The internal structure of asteroids is still poorly known and has never been measured directly. Our knowledge is relying entirely on inferences from remote sensing observations of the surface and theoretical modeling. Is the body a monolithic piece of rock or a rubble-pile, an aggregate of boulders held together by gravity and how much porosity it contains, both in the form of micro-scale or macro-scale porosity? What is the typical size of the constituent blocs? Are these blocs homogeneous or heterogeneous? The body is covered by a regolith whose properties remain largely unknown in term of depth, size distribution and spatial variation. Is it resulting from fine particles reaccretion or from thermal fracturing? What are its coherent forces? How to model its thermal conductivity, while this parameter is so important to estimate Yarkowsky and Yorp effects?

After several asteroid orbiting missions, theses crucial and yet basic questions remain open. Direct measurements of asteroid deep interior and regolith structure are needed to better understand the asteroid accretion and dynamical evolution and to provide answers that will directly improve our ability to understand and model the mechanisms driving Near Earth Asteroids (NEA) deflection and other risk mitigation techniques. There is no way to determine this from ground-based observation. Radar operating from a spacecraft is the only technique capable of achieving this science objective of characterizing the internal structure and heterogeneity from submetric to global scale for the benefit of science as well as for planetary defence or exploration [1].

Low Frequency Radar: The deep interior structure tomography requires low-frequency radar to penetrate throughout the complete body. The radar wave propagation delay and the received power are related to the complex dielectric permittivity (i.e to the composition and microporosity) and the small scale heterogeneities (scattering losses) while the spatial variation of the signal and the multiple paths provide information on the presence of heterogeneities (variations in composition or porosity), layers, ice lens. A partial coverage will provide "cuts" of the body when a dense coverage will allow a complete tomography. Two instruments concepts can be considered: a monostatic radar like Marsis/Mars Express (ESA) that will analyze radar waves transmitted by the orbiter and received after reflection by the asteroid, its surface and its internal structures; a bistatic radar like Consert/Rosetta (ESA)

[2],[3] that will analyze radar waves transmitted by a lander, propagated through the body and received by the orbiter (Figure 1). Monostatic radar requires very low frequencies in the range 10 to 20 MHz necessitating the use of large antennas and is more consuming in term of mission resources (mass, data flow), driving all the mission specification. On the other hand, bistatic radar can use slightly higher frequencies in the range 60-90 MHz, simplifying the accommodation on mission carrying a surface package. This concept is fully compliant with medium class planetary missions.

2096.pdf 47th Lunar and Planetary Science Conference (2016)

AIM will carry Mascot2, a lander inheriting from Mascot/Hayabusa2 [6] to land on Didymoon. On Mas-cot2 and AIM, the bistatic will porbe the Didymoon's internal structure, with a typical resolution of 30 meters to characterize the structural homogeneity in order to discriminate monolithic structure vs. building blocks, to derive the possible presence of various constituting blocks and to derive an estimate of the average complex dielectric permittivity, which relates to the mineralogy and porosity of the constituting material. Assuming a full 3D coverage of the body, the radar will determine Didymoon's 3D structure: deep layering, spatial variability of the density, of the block size distribution, of the average permittivity.

When AIM is combined with DART, the bistatic experiment will characterize possible structural modification induced by DART impact. It will also support mass determination and orbit characterization with range measurements during and after descent. Finally, it will contribute to the characterization of the primary (called "Didymain").

On AIM mothership, the shallow subsurface radar will determine the structure and layering of Didymoon and Didymain shallow sub-surface down to a few meters with a metric resolution and will map spatial variation of the regolith texture which is related to the size and mineralogy of the constituting grains and macroporosity and spatial distribution of geomorphological elements (rocks, boulders, etc) that are embedded in the subsurface.

With DART, it is a key instrument to assess the regolith tomography before and after impact in order to characterize the crater topography, the internal structure modifications and the mass loss. Then it is also a monitoring of the impact ejecta generated by the collision of the DART spacecraft in the vicinity of the secondary asteroid in order to estimate size distribution, speed, and total mass.

It will also contribute to shape modeling, mass determination and orbital characterization with altimeter mode. And finally, more prospective objectives will be considered, such as the support to ground-based radar measurements like Arecibo or Goldstone: orbital radar measurement is indeed a unique opportunity to crossvalidate ground-based NEA characterization with radar signal in the same frequency range and with better resolution, better SNR and more favorable geometry.

[1] Herique A. et al. (2015)

Figure 1 :

 1 Figure 1: Illustrations of Bistatic tomographic Hight Frequency Radar: Imaging the first tens meters of the subsurface with a metric resolution to identify layering and to reconnect surface measurements to internal structure can be achieved with a higher frequency radar on Orbiter only with a 300MHz -800MHz frequency range typically. An enlarged frequency range up to 3GHz, like the WISDOM radar designed and developed for the ESA ExoMars 2018 Rover mission [4] adds valuable science return contributing into the shape modeling, mass estimation, orbital parameters determination or close asteroid navigation with an altimeter mode. AIM Mission: Bistatic tomography radar and high frequency radar are presently under phase A/B1 study in the frame of the ESA's Asteroid Impact Monitoring mission. AIM as a stand-alone mission or constituting the Asteroid Imact & Deflection Assessment (AIDA) with the Double Asteroid Redirection Test (DART) mission under study by APL is a mission to characterize "Didymoon", the secondary of the binary NEA Didymos and to contribute to the evaluation of impact mitigation strategies [5],[6].

Lunar and Planetary Science Conference (2016)

  , A direct observation of asteroid's structure fromm deep interior to regolith: why and how? 4th IAA Planetary Defense Conference. -PDC. [2] Kofman, W. et al. (2007) The Comet Nucleus Sounding Experiment by Radiowave Transmis-sion (CONSERT). A short description of the instrument and of the commissioning stages. SSR, Vol 128. [3] Kofman,W., et al. (2015) Internal properties of the nucleus of Comet 67P/Churyumov-Gerasimenko as measured by the CONSERT instrument,

	47th		2096.pdf
			Science 349
	(6247). [4] Ciarletti V. et al. (2011), WISDOM a GPR
	designed for shallow and high resolution sounding of
	the Martian subsurface, Proc. of the IEEE. [5] Michel
	P. et al. (2016), Asteroid Impact Mission (AIM): the
	european component of the AIDA space project. LPSC
	2016 [6] Cheng A. F. et al. (2016) Asteroid Impact &
	Deflection Assessment Mission: Kinetic Impactor, PSS
	in press. [7] Ulamec, S., et al. (2014), Landing on
	Small Bodies: From the Rosetta Lander to MASCOT
	and beyond; Acta Astronautica 93, pp. 460-466.
		HFR	LFR
	Target	Regolith	Deep interior
	Radar	Monostatic	Bistatic with Mascot
	Signal	Step Frequency	BPSK
	Heritage	WISDOM / Exomars	CONSERT / Rosetta
	Wave propaga-tion	Reflec-tion/Scattering	Transmission
	Antenna	Vivaldi	Dipoles
	Tx Polarization	circular	Linear
	Rx Polarization	2 linear	Circular
	Penetration (m)	10 to 20	170
	Resolution (m)	~2	15 -30
	Nom. BW (MHz)	300-800	50-70
	Ext.. BW (MHz)	300-2500	45-75
	PRF (Hz)	0.5-5	0.1-0.3
	Step durat. (µs)	60	0.05
	Integration	3000	261120
	Tx Power (W)	20	16
	Primary pow (W)	88	10 (each unit)
	data rate (kbit/s)	300	5


Nahed Jalloul
email: nahed.jalloul@univ-rennes1.fr

Fabienne Porée

Geoffrey Viardot

Philippe L'hostis

Guy Carrault

Activity Recognition Using Multiple Inertial Measurement Units

Keywords: Bioinformatics, Biomedical sensors, Signal Processing

Objectives: This paper addresses the design of an ambulatory monitoring system based on a set of wearable, wireless inertial measurement units able to perform activity recognition for healthy individuals and Parkinson's disease patients, as well as analyze and assess the severity of levodopa induced dyskinesia.

Material and methods: The monitoring system is composed of six Shimmer3 modules placed at different positions of the individual's body. Both healthy individuals and one patient performed a protocol of simple daily life activities while wearing the Shimmer3 modules. As an initial step, validity of the monitoring system in identifying healthy individuals' activities is presented. Data corresponding to the activities was separated and features in both the time and frequency domains were extracted. Multiple factor analysis was used to evaluate and infer the relationships between the different module positions. A method of feature selection was implemented to determine the most important features, positions and sensors included in the different modules. The classification of activities was done using a KNN classifier. 1 Results: Promising results were obtained in classifying the activities of healthy individuals, with a global accuracy of 77.6%, while certain adaptation is required for the application on Parkinson's disease patients.

Conclusion:

While activity recognition for healthy individuals using this system was successful, further evaluation of the contribution of each module needs to be done in order to determine optimal module positions. To validate the obtained results on Parkinson's disease patients, more acquisitions on patients are required.

Introduction

In the past few decades, activity recognition through wearable sensing devices has become a widely researched area finding applications in various domains such as medicine, security and entertainment. In particular, several medical applications have been concerned with utilising these non-invasive, wearable devices for detecting motor symptoms of patients.

Parkinson's disease is the second most common neurodegenerative disorder, following Alzheimer's [START_REF] Alves | Epidemiology of Parkinson'â Ćňâ Ďćs disease[END_REF]. This disease is characterised by several motor symptoms such as hypokinesia, bradykinesia, freezing of gate, and tremor as well as motor complications of treatment methods such as dyskinesia. The assessment of the motor symptoms today is based on the clinician rated scales and patient diaries. Although some of these scales, such as the unified Parkinson's disease rating scale, have proven to have reliable test-retest performance, they remain to be subjective, short-based assessments that lack the ability to quantify the disease's motor complications. Therefore, research has been aimed towards finding a more objective, longterm automated method for the detection and quantification of Parkinson's disease motor symptoms.

Recent studies have been conducted that aim to detect certain motor symptoms such as tremor, bradykinesia and hypokinesia using kinematic sensors. However, these studies are usually concerned with the assessment of one or two motor symptoms and are based on short time recordings of patients in somewhat restricted environments (see [START_REF] Hobert | Technical and clinical view on ambulatory assessment in Parkinson's disease[END_REF] for review). The enhancement of therapeutic measures lies in understanding the complexity and evolution of motor complications that manifest both due to disease progression as well as the implementation of the therapeutic plan itself. In order to do that, there is a need for a long-term and reliable monitoring device that is able to assess the patient's general motor state, detect and quantify motor symptoms, providing an objective means for determining disease progression [START_REF] Zwartjes | Ambulatory monitoring of activities and motor symptoms in Parkinson's disease[END_REF]. Because of that need, there exists a link between the field of ambulatory monitoring of Parkinson's disease motor symptoms and the field of activity recognition.

Although there have been many advances in the field of activity recognition, there are several issues that still remain urging further exploration of the currently available technologies in enhancing the accuracy of classifying daily activities [START_REF] Lara | A survey on human activity recognition using wearable sensors[END_REF]. Some of the most commonly used types of sensors in activity recognition today are inertial sensors such as accelerometers and gyroscopes, magnetic sensors, and physiological signal measurement devices. Studies have reported varying results based on these different types of sensors, mostly due to differences in the types of activities to be recognised and depending on the application at hand. The objective of this paper is to propose the design of a new ambulatory monitoring system based on a set of wearable, wireless inertial measurement units able to perform daily life activity recognition for healthy individuals and Parkinson's disease patients.

Section II introduces the system of data acquisition, along with the placement of modules, the incorporated sensors and the activity protocol performed. Following that, the signal processing method is explained and the extracted features are described, including a brief summary of multiple factor analysis, the feature selection method and K Nearest Neighbour (KNN) classification. Section III presents the obtained results for both multiple factor analysis and KNN classification, while Section IV shows the results of the primary application of the developed algorithm on a Parkinson's disease patient. In Section V, a brief conclusion is presented.

Materials and Methods

Data acquisition

The entire system of acquisition is composed of six Shimmer3 modules [START_REF] Burns | SHIMMER T M -A wireless sensor platform for noninvasive biomedical research[END_REF] placed at different segments of the subject's body (Figure 1). The modules were attached to the subject's body by elastic straps and did not affect the subject's free movement. Each of the sensing modules contains the five following sensors: triaxial low noise accelerometer, triaxial wide range accelerometer, triaxial gyroscope, triaxial magnetometer, temperature sensor, and altimeter. The sensors all operated under the same sampling rate of 512 Hz and data was collected on each module's internal memory card. For the analysis of movement, only the signals collected from both accelerometers, gyroscope and magnetometer were considered. Twelve different signals were therefore collected from each shimmer module, and a total of seventy-two signals were analysed for all module positions. Nine healthy subjects (3 males and 6 females) participated in the data acquisition protocol which consisted of a series of N = 7 simple daily life activities (see Figure 2). The activities were performed during the session without interruption and data acquisition was done continuously. Each subject performed the same protocol twice, on two different days. During the acquisition sessions, subjects were instructed to act freely and perform each activity at their own pace. The period of each activity was timed for an average of two minutes and subjects were then instructed to move on to the next activity. The sessions were video-taped, and the exact start and end times of each activity were annotated based on the collected video recordings.

Feature extraction

Signals were then processed according to the schema in Figure 4. Each activity's raw signal data was first extracted from the total acquired data based on that activity's annotated start and end times. Then, the data corresponding to each raw signal was segmented into M = 7 second time windows.

Then, for each time window, x[n], n = 1,..., Mx512, features in both the time and frequency domain were extracted.

• DC and Std are the mean and standard deviation of each sensor data taken over the specified time window. Prior to calculating the frequency domain parameters, the zero-frequency component of the collected signals was removed.

• Energy was calculated as the sum of squared discrete Fourier transform component magnitudes for each time window [START_REF] Bao | Activity recognition from user-annotated acceleration data[END_REF].

• In order to obtain the mean and median frequencies, denoted as MNF and MDF respectively, the power spectral density of each window was first estimated using the Burg method [START_REF] Huang | Study of Autoregressive (AR) Spectrum Estimation Algorithm for Vibration Signals of Industrial Steam Turbines[END_REF], which fits an autoregressive model to the signal by minimising the forward and backward prediction errors.

The mean frequency is defined as the average of the obtained power spectrum, while the median frequency is defined as the value dividing the power spectrum into two equal areas. The order of the autoregressive model was selected after an exhaustive search on all the database using the Akaike information criterion.

The total number of features calculated for each module is therefore 96 (8 features x 12 signals).

Multiple Factor Analysis

Multiple factor analysis is a geometrical approach to analyse data comprised of a number of observations, in which each observation is described by variables of quantitative and categorical nature simultaneously [START_REF] Escofier | Multiple factor analysis (AFMULT package)[END_REF]. Different sets of variables can be organised into groups, given that each group contains variables of the same nature. This method is based on performing factor analysis on the entire set of variables while being able to apply weighting to the observations or different groups [START_REF] Hair | Multivariate data analysis[END_REF]. This leads to a representation of the observations and variables, resulting in a categorisation of the observations, a categorisation of the variables, and the association of the observations to each group of variables. The factor analysis performed is either equivalent to principle component analysis, in the case of quantitative variables, or multiple correspondence analysis, in the case of categorical variables.

Feature Selection and KNN classification

In certain cases, a large number of descriptive variables can have a negative effect on the performance of classification systems. Feature selection is the process of selecting a subset of variables that best distinguish between the different classes of a set of observations [START_REF] Jain | Feature selection: Evaluation, application, and small sample performance[END_REF]. This method provides the ability to reduce the number of variables and allows only the most pertinent variables to be used for classification. It is possible to achieve better classification results using this process since the noise caused by unimportant or redundant variables is eliminated. The Relief Method, is a feature selection tool that works by estimating how well each of the provided variables can separate between the given classes. According to their weights, the variables are then re-arranged in a decreasing order of importance [START_REF] Kira | The feature selection problem: Traditional methods and a new algorithm[END_REF]. In order to test the effect of feature selection on the classification accuracy, several datasets will be considered. The first, consisting of all calculated variables, and then datasets that consist of a varying number of selected variables by the relief method to find the optimal number of variables.

KNN classification [START_REF] Cover | Nearest neighbor pattern classification[END_REF] was chosen due to its simplicity and the fact that it is one of the most fundamental classification methods when there is little or no prior knowledge about the distribution of data. This method is based on the Euclidean distance measure between a test sample and the specified training samples.

The predicted class of any test sample is set equal to the most frequent true class of the k nearest surrounding training samples.

Results

Multiple Factor Analysis

The previously mentioned variables were computed for the subjects activity protocol data performed on the first day. Each time window corresponding to the different performed activities is considered as an observation, and the variables belonging to the same module position are considered as a single group.

The total matrix is composed of 2266 observations x 576 variables. Initially, the groups of variables are equally weighted and a representation of the relationship between the different classes of observations and the groups of variables, or in other words the module positions, is obtained. The factor analysis is applied to the data including all groups of variables. In Figure 5, the mean of observations belonging to each class is projected on the common factors. Using this projection, a clear separation between the classes of activities is observed. On the one hand, when projecting the mean values on the first factor, dynamic activities, such as walking, and static activities, such as standing still or lying down, show a large difference. On the other hand, a moderate separation can be observed between the static activities when projected onto the second factor.

Feature Selection and KNN Classification

Classification of the activities was based on the extracted features using a KNN classifier. The classifier was trained and tested using two different approaches.

In the first approach, no feature selection was implemented and the computed features corresponding to all module sensors and positions were included. The classifier was trained on all subjects activity protocol data performed on the first day, and tested using the subjects activity protocol data taken from the second day.

This resulted in a global accuracy of 53% for the activity classification, with the best performance obtained for the activity of walking which had sensitivity and specificity values of 100% and 87.1% respectively. The poorest performance obtained corresponded to low-energy activities such as reading and standing still.

In the second approach, feature selection was applied to the subjects activity protocol data performed on the first day. The classifier was then trained by the selected features and tested using the same features computed for the subjects activity protocol data of the second day. The obtained rankings of the variables based on the given weight by the relief method is shown in Figure 6. Each time, a different number of ranked variables by the relief method was chosen, and the classification process was repeated. The resulting classification accuracy of all the activities for a different number of selected features is represented in Figure 7.

The best overall classification accuracy of 77.6% was obtained when classification was done using the first 20 ranked variables by the relief method. Table 1 shows the sensitivity and specificity obtained for each activity classification and the contingency matrix.

Application on Parkinson's Disease Patient

A primary evaluation of the developed algorithm was done using data from a Parkinson's disease patient who performed the same protocol of activities while wearing the sensing modules. During the recording session, the patient suffered from varying intensities of Levodopa induced dyskinesia [START_REF] Fabbrini | Levodopa-induced dyskinesias[END_REF]. In certain cases, the dyskinesia was at a debilitating severity, preventing the patient from performing the activity. In such cases, the protocol was paused for a few minutes giving the patient enough time to recuperate before moving onto the next activity. The patient displayed dyskinesia primarily in the left side of the body. It was also evident, during the performance of the protocol, that certain activities were difficult to perform while experiencing dyskinesia, while others were relatively unaffected by it.

Similarly to the analysis performed on the data collected from healthy subjects, MFA was performed on the data collected from the Parkinson's disease patient. Dyskinesia was prominent during lying down, standing and sitting still for this particular patient. As a result, there were noticeably significant movements during these activities that distinguished between their nature when performed by the patient, and the static nature in which they were performed by the healthy subjects. The mean of the observations belonging to each class corresponding to the patient's data is projected onto the first and second factors (Figure 8).

The same process of feature selection using the relief method and KNN classification was also implemented on the data collected from the patient. However, a leave one out training and testing method was applied here due to the lack of available data. Without feature selection, and using the features extracted from all module sensors and positions, an overall classification accuracy of 43.8% was obtained. The activity walking displayed the highest sensitivity and specificity values of 77.8% and 83.96% respectively, while the activity of lying showed the poorest classification sensitivity and specificity, in the majority of cases being misclassified as walking. These results must be considered as preliminary due to the small size of observations compared to the number of features. Nevertheless, for our final goal which is the detection and classification of dyskinesia and akinesia, these results are promising, demonstrating that during a dyskinetic episode, we were able to separate between dynamic and static activities.

Discussion-Conclusion

Preliminary evaluation of activity classification using the Shimmer3 inertial measurement units shows promising results. The sensing modules small size, portability and simple operating manner make them ideal for the design of an ambulatory system for long-term monitoring of activities and the motor states of patients.

In the case of healthy subjects, the sensing modules were able to successfully classify the different activities with a global accuracy of 77.6%. Further evaluation of the contribution of each module position with respect to the performed activity needs to be done in order to determine the optimal placement of the modules. Following this step, further features will be explored and more complex classification methods will be tested, such as Support Vector Machine, to enhance the accuracy of identifying different activities in healthy subjects.

Optimizing the activity classification process is an important step before testing the system's ability to detect and estimate Parkinson's disease motor symptoms. Disease progression and the prevalence of motor symptoms, specifically Levodopa induced dyskinesia, will affect the patient's general motor condition and their ability to perform simple daily activities. This change will be significant in estimating the severity of motor symptoms and objectively monitoring the patient's disease progression.

Figure 1 :

 1 Figure 1: Shimmer3 module positions.

Figure 2 :

 2 Figure 2: Description of the protocol including six module positions, six sensors (two not used) and seven activities.

Figure 3

 3 Figure 3 shows examples of collected signals from sensors of a single module (Ankle) during two different activities (walking and lying down).

Figure 3 :

 3 Figure 3: Example of collected signals from a single module's sensors during walking and lying down.

Figure 4 :

 4 Figure 4: Data processing and feature extraction.

Figure 5 :

 5 Figure 5: (A) Projection of the mean of observations of each class on the first and second axes obtained after performing MFA on the data obtained from healthy subjects. (B) Total variance for each type of activity for the first (blue) and second (red) axes.

Figure 6 :

 6 Figure 6: Variable rankings by decreasing order of weight.

Figure 7 :

 7 Figure 7: Classification accuracy obtained from different numbers of selected variables. The best performance was found for R = 20 variables.

Table 1 :

 1 Classification results obtained using relief feature selection with R = 20 features and KNN classification and their corresponding sensitivity and specificity values.

	Activity	Walking Standing	Lying	Sitting	Writing	Reading	Eating	Total
	Walking	234	32	0	5	4	1	0	276
	Standing	75	176	0	3	0	0	0	254
	Lying	0	2	510	4	0	4	0	520
	Sitting	5	53	0	442	3	1	37	541
	Writing	0	0	0	1	252	0	9	262
	Reading	0	10	0	35	58	56	0	159
	Eating	1	0	0	88	33	44	88	254
	Total	315	273	510	578	350	106	134	2266
	Sensitivity	84.78	69.29	98.08	81.7	96.18	35.22	34.65
	Specificity	95.93	95.18	100	92.12	95.11	97.63	97.71

Acknowledgment

The authors gratefully acknowledge volunteers and technicians who participated in the acquisition process. This work was funded by Biotrial.

Conflicts of interest

No conflict of interest.