Jason Lambert

Gilles Rougeron

Lionel Lacassagne

Sylvain Chatillon

A Fast Ultrasonic Simulation Tool based on Massively Parallel Implementations

Keywords: Ultrasonic Inspection Simulation, Graphical Processing Units, Pencil Method, Kirchhoff Model

This paper presents a CIVA optimized ultrasonic inspection simulation tool, which takes benefit of the power of massively parallel architectures : graphical processing units (GPU) and multi-core general purpose processors (GPP). This tool is based on the classical approach used in CIVA : the interaction model is based on Kirchoff, and the ultrasonic field around the defect is computed by the pencil method. The model has been adapted and parallelized for both architectures. At this stage, the configurations addressed by the tool are : multi and mono-element probes, planar specimens made of simple isotropic materials, planar rectangular defects or side drilled holes of small diameter. Validations on the model accuracy and performances measurements are presented.

INTRODUCTION

Inspection simulation is used in a lot of non-destructive evaluation (NDE) application : from designing new inspection methods and probes to qualifying methods and demonstrating performances through virtual testing while developing methods. The CIVA software, developed by CEA-LIST and partners is a multi technique platform (Ultrasonic Testing (UT), Computed Tomography and Radiographic Testing) used to both analyze acquisitions and to run simulations validated against international benchmark [START_REF] Mahaut | An overview of Ultrasonic beam propagation and flaw scattering models in the CIVA software[END_REF]. It benefits from semi-analytics models and user friendly GUI to design an inspection case with the help of CAD tools and imaging. However, due to the potential complexity of the configurations, current semi-analytical models based on asymptotic developments takes minutes to hours to run a simulation. Significant efforts are made in the CIVA framework to reduce computation times by working both on models and computational aspects. The works reported here ties in with this general approach.

Nowadays, new massively parallel architectures can empower computational softwares, at the expense of adapting algorithms to the specificities of those architectures to highlight and improve parallel computation steps. The exploitation of both multicore-GPP and GPU capabilities results in a new intensively parallel algorithm of simulation of UT inspection [2] [3]. This new model relies on analytical solution to the beam propagation. Its implementations on both architectures use high performances signal processing libraries (Intel MKL on GPP and NVidia cuFFT on GPU). The new model is, however, limited to canonical configurations due to the strict parallelism requirements and to the lack of genericity of analytical beam propagation.

This paper is organized as follows. The first section presents the new model with the specificities of the anylitical solutions and the resulting limitations to the model. Then, GPP and GPU implementations are detailed. In the next section, the validity of this model is tested against real sets of data from the benchmark held at QNDE2013 conference, conducted by the WFNDEC (World Federation of NDE Centers). Furthermore performances are measured with both implementations on high end hardware to provide some sort of fair comparison. Finally, the last section concludes the paper and discusses future works.

UT INSPECTION SIMULATION

This section is divided in three distinct parts toward presenting the model of the fast and massively parallel ultrasonic simulation. An inspection simulation rely on three steps :

• computation of the field transmitted by the probe on the flaw, • diffraction by the defect,

• computation of the field on the sensor for echo synthesis.

Due to Auld's reciprocity principle, the field diffracted by the defect to te receiving probe can be obtained by coputing the field transmitted by the receiver (used as an emitter) on the defect. This allow to group computations of the field in reception and in emission together. Transmission of the field transmitted by the probe on the flaw, and field/flaw interaction for received echo synthesis.

Field Computation Using Paraxial Ray Tracing Method

The field computation relies on the pencil method, a generic approach for heterogeneous and anisotropic structures. By evaluating the ray path of the beam, from the transducer to the observation point, it is possible to evaluate the time of flight and the amplitude of the contribution of the beam using energy conservation principle on the tube. As seen in Figure 1, the beam may propagate through different materials and cross multiple interfaces : its contributions Ψ is determined by the propagation matrix obtained by multiplying the elementary contributions of each section of the pencil with the initial contribution Ψ 0 , as shown in equation 1 [START_REF] Gengembre | Pencil method for ultrasonic beam computation[END_REF]. In general case, there is no simple solution to determine the ray path from the source to the computation point. However, in isotropic and homogeneous structures, analytical methods can be used to determine ray path. In direct mode, for standard geometrical surfaces, it is possible to determine a polynomial modelizing the path following Snell-Descartes. The roots of this polynomial, whose degree vary from 4th for planar surfaces to 16th for torical surfaces, correspond to the possible solutions for the ray path. They are determined numerically, through Newton's and Laguerre's Method [START_REF] Lambert | Performance evaluation of total focusing method on GPP and GPU[END_REF]. In the case of half-skip mode, ray path can only be determined analytically for planar surfaces and backwalls, through two dimensional Newton's method solving. Those analytical method allows for a greater regularity benefitings to the requierments of massively parallel architecutes.

Ψ = L prop1 .L rot1 .L inter1 .L prop2 .L rot2 .L inter2 ...L propn .L rotn .L intern * Ψ 0 (1)
Those limits results in an analytical model dedicated only to a specific set of geometries, where the propagation matrix can be fully determined preemptively thus avoiding costly matrix multiplication, for direct and half skip mode. To benefit from highly parallel architecture, this model will rely on a regular pencil distribution over the sampled transducer and a numerical resolution of the analytical equations.

Once the pencils computed on the whole transducer surface, their summation with the application of delay laws results in the computation of the impulse response, as seen in Figure 2. As this model addresses defect response simulation, the field is outlined as the maximum of amplitude, the phase, the corresponding time of flight, the direction of propagation and the wave polarisation.

UT Diffraction Simulation using Kirchhoff Coefficient

UT fields are computed, both in emission and reception (using Auld's reciprocity theorem), over the studied defect, on a coarse grid.

To sum up the contribution over the defect surface, a finer computation is required (on a finer grid, one order of magnitude finer). The field data (amplitude, phase, time of flight, polarization and direction of propagation) are interpolated from the grid and are used to compute the UT diffraction. This model uses the Kirchhoff model, whose coefficients are valid near specular reflection [START_REF] Darmon | Modelling of scattering of ultrasounds by flaws for NDT[END_REF]. All the elementary diffraction contributions computed on the defect contribute to the echo impulse response.

The echo signal is the result of the convolution of the echo impulse response signal an the reference signal.

OVERVIEW OF FAST IMPLEMENTATIONS

In this section, the choices of implementation on both architecture will be discussed under the specificities of the hardware.

GPP Characteristics

Modern day GPP are composed of multiple general purpose cores. Those are aimed at executing independent, heavyweight, tasks. GPP disposes of two parallelism levels.

• A fine grained parallelism relying on specific SIMD instructions (Single Instruction Multiple Data) which execute the same operation on short vectors (128 to 512 bits). For example, with 128-bit vectors, a SIMD instruction can perform simultaneously four additions on four single-precision floating point numbers (32-bit). • A coarse grained parallelism relying on multithreading to enable multiple logical tasks to reside simultaneously on the GPP. The OpenMP API is aimed at shared-memory parallelism : it creates a thread per GPP core, each with its own stack where locals variables are located but they can also communicate through some shared variable. Its work distribution rely on a succession of sequential sections (with only one active thread) with parallel sections (with all threads active) assembled in a fork-join fashion.

GPP UT simulation implementation

On the GPP, the computations are regrouped by coarse step in order to even the load on each core of the GPP and in order to maximize the reuse of data describing the simulation.

On each step, the OpenMP API is used to parallelize the computations. Moreover the Intel Math Kernel Library (MKL) is used to benefit from fast FFT implementation on the GPP. This library benefits from SIMD instructions available on the processor and rely on a precomputed "plan" to prepare the coefficients for a specific signal size. Once precomputed, those can be reused as needed for the computations. In order to benefits from the best performances, the signal size has to be determined first and used repeatedly throughout the simulation.

Signal processing is divided in two main steps, each relying on a specific signal size : one concerns the field data computation and the second is responsible for the response signal computation.

Figure 3 represents the algorithm developed for the GPP.

NVIDIA GPU Architecture Considerations

NVIDIA released a new programming architecture in 2006 to ease GPGPU on its GPU. This architecture comes with both a programming paradigm and a C-like interface, called Compute Unified Device Architecture (CUDA), to program GPU. A CUDA task is called "kernel" and it is represented by a C-like function which will be executed once by all threads. With CUDA, GPU are used as coprocessors to their host systems. They dispose of their own memory, off-chip DRAM, called global memory, which is slow to access memory (comparing to GPU frequency), but it is kept alive for the duration of the application allowing communication between kernels without back-and-forth data movement with the host.

A GPU consists of severals Streaming Multiprocessors (SM, 16 on current high end cards). Each one is composed of multiple elements, whose numbers depend on hardware generation, for example :

• CUDA cores for integer and floating-point arithmetic operations;

• one or two schedulers to manage threads. Beside these computing elements, each SM disposes of :

• a rather small (16KB, up to 48KB on recent devices) but fast shared memory on the SM, for exchanges between its cores; • Special Function Units (SFU) for floating-point transcendentals functions; • a limited, per SM, set of registers divided between the threads of a kernel at compile-time.

CUDA multiprocessors use the Single Instruction Multiple Thread (SIMT) architecture : each cycle, one instruction is executed by all the cores of a single SM. Threads are grouped into packs of 32 threads, called warps, to be scheduled by each SM. A SM can host multiple warps up to its physical limits, depending of the required configuration form the kernel.

This model is most efficient when each thread in a warp executes the same instruction. In case of divergence, each branch is consecutively executed by the warp, with threads following the other branch doing nothing. Moreover, memory accesses to the different memory of a GPU benefit from data locality by optimizing for contiguous threads of a warp. Each thread disposes of its own program counter and of its own registers. Thread contexts of warps stay on the multiprocessor for the duration of their execution enabling a cheap and fast hardware scheduling of ready warps in spite of execution time and/or latency.

Implementation on the GPU

Due to the high need for computational regularity on GPU, the control flow of this implementation will aim at executing the same task over multiple data and proceed in steps, each corresponding to dedicated kernel. To benefits from highest GPU performances, kernel were developed to use single-precision floating point opretaions (IEEE-754). To perform signal processing operations, this implementation relies on the optimized cuFFT library which is optimized to perform efficiently on large batches of signals. This library establishes a "plan", consisting of the precomputed coefficients which are applied to the signals. It is most efficient when performed repeatedly on a batch of signals of the same dimensions (number of signals and signal length), to reuse a previously computed plan.

The following Figure 4 illustrate the overall necessary computations steps detailed as follow.

Loop over transducer positions while applying delay laws 3x, on x,y,z

Signals size determination

The first step of the computation determines the temporal width of temporary data using analytical ray path resolution for time of flight computation. Once the temporal width of the signals is determined, this information is sent back to the host for allocating memory and cuFFT plan creation. This size is the same for all the temporary signals needed for field computation to allow the use of a single cuFFT plan to compute FFT on large batches of signals.

Field data computation

Because of GPU memory limitations, with high end consumer grade hardware disposing of 1.5 GB of memory at most, the whole field computation has to be done on a subset of sensor positions due to the size of temporary data. The following steps are executed on the different subsets of positions by a loop controlled by the host.

For each position, elementary contributions of all the pencils in a specific field point need to be summed following the delay law. The kernel assigns a couple of position/field point to a single CUDA block, to first compute pencils data and store this data in their own shared memory. Threads then loop over the different shots of the delay law to apply delays to the pencils, fetched from shared memory, and proceed with the summation to the signals residing in global memory. It results in completely random accesses to the signal, as there is no way to predict if two distinct pulses will collide in memory due to temporal span : this summation is done using atomic operations. The result of this summation is three impulse responses, one for each displacement coordinate, which are then convoluted with the reference signal.

The main direction is extracted by a new CUDA kernel from those convoluted signals, as a reduction in shared memory performed collaboarively by threads from a single block. This direction is saved for echo computations. The scalar projection of the displacement against this direction is computed.

The envelope of this signal is computed through the use of FFT. Its maximum of amplitude is extracted by another CUDA kernel. Its extraction is again done through a contiguous scan of the signal and a collaborative reduction through shared memory. It is used to determine the phase shift of the field propagation and to extract the associated time of flight.

Impulse response computation

For each shot and for each sensor position, a response signal is computed. A kernel computes, for a defined mode, the interaction of the received and the emitted field.

The defect is modelized by a set of fine grained sampling point, on witch field information are bilinearly interpolated from previously computed field data. Threads of a block works collaboratively over the whole set of defect points to compute the elementary interaction and sum this data to the global response signal, for a given sensor position and a given shot. As for summing elementary contribution to construct UT field signal data, due to the random nature of this summation to the impulse response, it is necessary to use atomic operations on the signal residing in global memory.

Once computed, the host calls to cuFFT back and forth to convolute the impulse response by the reference signal.

PERFORMANCES OF THE MASSIVELY PARALLEL IMPLEMENTATIONS

To provide a realistic analysis of the results obtained by both implementations, among the limits of this model, realistic set of simulations are studied from the benchmark session of the QNDE 2013 conference.

Test cases

The selected cases consist of planar parts of homogeneous isotropic steel, inspected at L45°for several focusing depths : FBH Serie1 A series of 10 3mm flat bottom holes; FBH Serie2 A series of 6 3mm flat bottom holes; Notch A back wall breaking notch residing on the planar back wall of the part. The characteristics of those configurations is summed up in Table 1.

Results Validations

To validate the accuracy of the new parallel model, the results of both of its implementations have been studied. As the CIVA 11.0 software is validated against real experiments through benchmarking, especially during the QNDE 2013 conference, its simulations will serve as a point of reference for result validation [START_REF] Toullelan | Results of the 2013 UT Modeling Benchmark Obtained with Models Implemented in CIVA[END_REF].

The dataset presents an important number of data, for each sensor position and each shot from the delay laws, a response signal is obtained. This analysis focuses on some of the most representative results to highlight the model accuracy on each dataset. Figure 5 illustrate the measurements done on the simulation of the FBH1 series, for the shot corresponding to the echo of maximum of amplitude. Those experiments all fall under 1dB from CIVA results, presenting a gap inferior to the requirements for passing benchmarks against experimental results. Thus, it is possible to validate the new model and its implementations.

Performances

Those two implementations have been benchmarked on high end hardware to aim at a fair evaluation of performances.

2× GPP Intel Xeon 5590. The machine contains 2 GPP, consisting each of 6 cores 3.47GHz with Hyperthreading (12 logical cores for each GPP) ; it also disposes of 24GB of memory 1 . This machine is also the one running CIVA 11.0 for reference. NVidia GeForce GTX580. This high end, consumer grade, GPU disposes of 512 CUDA Cores running @1.544Ghz each, and rely on the Fermi architecture. It disposes of 1.5GB of device memory 2 . The benchmarked performances do not include GPU initialization nor the transfer of the result from the device back to the host (at most 200ms on the studied cases). Those two elements are dismissed as this implementation is aimed at a fast simulation, memory transfers to the host will be overlapped with computations and the GPU is only initialized once for the whole set.

Experiments. Impulse response simulation is run on the two architectures and on the three sets of data. To provide a reference for current simulation platform performances, the dataset are also run on the CIVA 11.0 software which rely on a general model. However, one has to notice that the CIVA 11.0 generic model is also benefits from multicore architectures through multithreading.

Table 3 shows details of the performances on both implementations. On real data, both implementation of the parallel model perform much faster than CIVA 11.0 generic model. It is noteworthy that most time is spent in signal processing routines : in delay law application to extract field information and in echo impulse response generation.

On the GPP, the signal processing step is costly due to the sheer quantity of memory operations required for computing FFT. Whereas on the GPU, the benefits from the higher parallelism through the cuFFT library are compromised by the cost of the atomic operation required to perform the summation of multiple contributions on a single signal.

CONCLUSIONS

This paper describes a new UT inspection simulation tool relying on the computational power of massively parallel architectures. This tool results have been validated on industrial use-cases through the QNDE 2013 benchmark by comparison with the current version of the CIVA software. The new model using an analytical beam propagation computation and computational regularity delivers an accurate simulation of the UT diffraction, under 0.8dB of differences on the benchmark. However, it is restricted to simpler configurations : planar surfaces, isotropic and homogeneous structures. By comparing its performances against CIVA11.0 generic model, the new model obtains an accelerations from ×31 to ×70 on GPP whereas on the GPU, a gain of ×41 to ×101 is shown. The performances shown with those first implementations indicate that interactive simulation is not out of range.

The current implementations are still limited to a subset of canonical configurations. Some extension should be realized to address more configurations, for example : side drilled holes ; other diffraction models which can benefits from the fast field and signal processing computations (GTD, SOV...) ; on the GPU, T mode computations.

Moreover, an analysis aimed at reducing the number of signal processing operations will be performed with the help of CIVA knownledge.

Finally, multiple improvements can be sought to speedup each parallel implementation. For example, the use of SIMD instruction on the GPP to benefit from its fine grained parallel abilities. Besides, by reducing the number of atomic operations in the GPU implementation, performances can improve drastically on signal summation.

Those implementation are the first step toward a fast UT inspection integrated within the CIVA 12 software.

FIGURE 1 .

 1 FIGURE 1. Ray tube visualization.

FIGURE 2 .

 2 FIGURE 2. Summation of the impulse response over the sensor.

FIGURE 5 .

 5 FIGURE 5. FBH1 Series -L45°-25mm depth focusing -FBH @ 35,25,15mm (insonated defects for the law.)

TABLE 1 .

 1 Studied configurations.

	Dataset	FBH Series 1 FBH Series 2	Notch
	# of defects	10 FBH	6 FBH	1 notch
	Probe	64 channels ↔ 128 samples
	# of position	361		201
	# of delay laws	7		4
	UT field sampling #	90	54	100
	defect sampling #	4140	2484	12321
	Mode	L direct	L direct	Half skip L

Table 2

 2 presents a summary of the accuracy validation realized on the new model against CIVA 11.0.

TABLE 2 .

 2 Accuracy analysis summary.

		FBH1	FBH2	Notch
	AMax Delay law	L45°-35mm depth focusing L45°-35mm depth focusing L45°-30mm depth focusing
	AMax gap CIVA/new model	0.2dB -0.7dB	0.6dB -0.8dB	0.6dB (GPP) -0.9dB (GPU)
	Time calibration	<1 sample -OK	<1 sample -OK	<1 sample -OK

TABLE 3 .

 3 Performances of full parallel implementations against CIVA 11.0 generic model.

		FBH-1 FBH-2	Notch
	generic model (CIVA 11.0)	241 s	142 s	363 s
		-	-	-
	Fast GPP (CIVA Dev)	7.6 s	2.0 s	7.3 s
	gain over generic model	×31.7	×71.0	×49.7
	Fast GPU (CIVA Dev)	5.8 s	1.7 s	3.6 s
	gain over fast GPP	×1.31	×1.17	×2.02
	gain over generic model	×41.5	×83.5 ×100.8

The compiler used is Intel®C++ Composer XE

Update 5 Integration for Microsoft Visual Studio 2010, Version 13.0.1211.2010[START_REF] Molero | Accelerating numerical modeling of wave propagation through 2-D anisotropic materials using OpenCL[END_REF] The CUDA implementation was compiled and executed with CUDA 5.0 toolkit


Abdellatif Chafik

Vincent Lecuire

Francis Lepage

Wireless Sensor Network Architecture for Monitoring Large Physical System in Cyclic Mobility

Keywords: monitoring, mobile wireless sensor network, routing algorithms, cyber-physical system, measuring transmitters

This paper addresses a sensor system dedicated to monitor a large mechanical machine. This sensor system is based on a wireless sensor network whose special feature is the mobility of one or more nodes following an invariant path traveled repeatedly. In the study described here the path is a circle. The topology chosen is that of locating relays on the mobile system that has appropriate geometric characteristics. A new dynamic routing algorithm is proposed. It uses the location and moving speed of mobile nodes to choose the optimal route of the packets forwarded to the fixed sink. Network performance is evaluated in terms of delivery time and packet loss rate and compared to this obtained by using flooding routing protocol. The results show that the proposed architecture and mechanisms deliver the appropriate quality of service for the monitoring of large physical systems.

INTRODUCTION

Monitoring a large mobile mechanical system from a fixed host is a great challenge at information sensing level. This one must indeed capture data in many places and transmits them in long distance while crossing a link mobile to fixe. Instrumentation can be designed either to monitor continuously the installation to diagnose malfunctions and to ensure its safety or to make measurements during a verification -certification phase. For the second purpose Wireless Sensor Networks (WSNs) and more precisely mobile WSNs are well suited for rapid instrumentation of existing installations with moving parts without large investments. The sensor or sensors board can measure parameters specific to the mobile part (e.g., vibrations, forces acting on a mechanical part) or environmental parameters throughout the trajectory of the mobile (e.g., wind speed, temperature) or a combination of both. Mobile WSNs are a subtype of mobile ad-hoc networks (MANETs). They are distinguished mainly by their sensors and their low power. Among the large mobile mechanical systems, we address those moving in an invariant way over a fixed path travelled cyclically. This specificity is expected to optimize the sensor network topology and two essential functions: localization of the moving sensor node(s) and mainly routing. Examples are transfer lines of products in industry, material ropeways in cement factory or a great ride, a chairlift or gondola in a leisure park or a resort.

Scenario and assumptions

Let us consider an observation wheel like the Singapore Flyer with a diameter of 150 meters and a total height of 165 meters. Passengers take place in 28 capsules moving at the speed of 0.24 m/s. The safety of such outstanding leisure equipment is fully mandatory. Wind speed is measured permanently but at one point only. Although the lateral movement and vibration of capsules have been carefully studied to be reduced at very low level it will be a safety improvement to easily measure them time to time and their variation with wind speed and capsule load. The corresponding system must perform measurements on the mobile part and send them to a fixed station for data analysis and storage. The sensor must be able to send its data immediately after sampling the information independently of its position. This is a requirement for a monitoring application. For a measurement application for a technical certification or for checking by a safety committee, a recording of measurements and a deferred processing is possible. However, a live processing has two major advantages for this application: one is saving time and second is the possibility of changing operating measurement scenario (e.g., other mobile speeds, other stimuli) if measures approaching some thresholds or if artifacts are detected. For multiple reasons, it is necessary to place the control station in a room on the ground: power supply, station environment, operator comfort. A wireless sensor network with mobile nodes seems well suited to carry out such a measurement system. It enables rapid deployment thanks to its two major advantages: lack of infrastructure and self-configuration.

Research objectives

The overall objective of our research is to design the sensor network system capable of capturing information in the right places, at suitable time, with sufficient precision, and of transmitting data to the control station in real time and without loss. For economic reasons and for easy maintenance, a constraint is to use a wireless sensor network available off the shelf. Thus processing resources, characteristics of the radio links and medium access protocol are predefined. WSN have the already mentioned quality of installation facilitated by the lack of infrastructure and by self-configuration. But their performance is limited in terms of bit rate, reliability of communications and message delivery delay. This performance is closely related to the topology of the network and the routing algorithm. This is why these two essential features must be designed making best use of geometry of the mechanical system and its movement with the aim of meeting the required performance. The topology should ensure permanent connectivity between the sensor(s) and the sink. The transfer from sensor to sink is not straightforward in the general case where the two nodes are not always linked by a direct radio channel. So the packet needs to be relayed by intermediate nodes according to the principle of ad-hoc networks. Whatever the route is used it includes a wireless link between the mobile part and the fixed part located at ground. The place of the mobile-ground connection defines two types of topology:

• A topology with mobile-ground connection at level of the sensor, the rest of the topology being fixed ground relays; • A topology with mobile-ground connection at level of the sink, the rest of the topology being relays embbeded in mobile nodes.

The second type of topology is the only one considered in this paper. It is well suited to mechanical systems consisting of spatially distributed subsystems moving together like materials ropeways. Routing determines the path of the data packet from the mobile sensor node to the fixed sink. Because topology is constantly changing, dynamic routing is required. The used routing criterion is the minimum of delay with the constraint of very low packet loss ratio. The proposed routing algorithm takes into account not only the location but the speed of the mobile too.

The rest of the paper is organized as follows. Section 2 describes the related works. The methodology and algorithms of dynamic packets routing in the network are described in Section 3. The performance of the proposed solutions is evaluated by simulations. The results are reported and discussed in Section 4. A conclusion and some perspectives complete this paper.

RELATED WORKS

Monitoring using Wireless Sensor Network

The use of WSN for monitoring is well known in several fields like military forces, environment, civil infrastructures often as an add-on on existing wired sensor networks. For example a WSN for Structural Health Monitoring has been designed, implemented, deployed and tested on the 4200 ft long main span and the south tower of the Golden Gate Bridge [START_REF] Kim | Health monitoring of civil infrastructures using wireless sensor networks[END_REF]]. WSN are used in the same manner in house monitoring and human heath monitoring [START_REF] Akyildiz | Wireless sensor networks: a survey[END_REF]]. These scenarii are very different from this of our research because almost the whole monitored system is static as well as all network nodes. Monitoring scenarii with mobile nodes connected using wireless technology are targeted by Cartel project at MIT [START_REF] Hull | Cartel: a distributed mobile sensor computing system[END_REF]] and MobEyes at UCLA [START_REF] Lee | Mobeyes: Smart mobs for urban monitoring with a vehicular sensor network[END_REF]] using VANET (Vehicular Adhoc Network). But this monitoring is adapted to a Delay Tolerant Network with non permanent link from a node to a sink and the nodes are powerful compared with those of the MWSN used in the architecture proposed in this paper. Machine monitoring using WSN are not frequently investigated in research papers. But some specific WSN dedicated to industrial automation have been developed by industry network manufacturers [START_REF] Christin | Survey on wireless sensor network technologies for industrial automation: The security and quality of service perspectives[END_REF]]. They are based on WSN protocol standards but focus on quality of service and security. Major contributions focused on application and used very simple WSN. Some uses the wireless capability to monitor electrical motors [START_REF] Lu | Online and remote motor energy monitoring and fault diagnostics using wireless sensor networks[END_REF]Gungor (2009), Lima-Filho et al. (2012)] but some address more global systems in order to optimize energy [START_REF] Salvadori | Monitoring in industrial systems using wireless sensor network with dynamic power management[END_REF]] or to monitor machines-tools especially the wear of the tools [START_REF] Tan | Distributed fault detection in industrial system based on sensor wireless network[END_REF]]. Due to the low number of hops in these networks most contributions improving QoS requirements for industrial uses focused on MAC layer [START_REF] Yigitel | Qos-aware mac protocols for wireless sensor networks: A survey[END_REF]]. According to this survey, we observe that majority of the protocols follow a service differentiation approach. Packets are treated according to their class and the corresponding requirements by tuning the associated network parameters at the MAC layer. WSN sensor nodes defaults could occur and so affect the monitoring. The localization of the faulty sensors has been addressed by [START_REF] Chen | Distributed fault detection of wireless sensor networks[END_REF]. Finally no contribution has been found addressing especially big machines monitoring by WSN.

Mobile WSN routing

Any network routing is based on one or several metrics. Position with respect the sink and residual energy are the main metrics used in WSN routing. Energy is not a critical issue in our scenario because the installation is temporary and the nodes are limited in number and are accessible. So batteries can be renewed and recharged. Among the routing using position metric geographic routing is an attractive routing protocol for WSN due to its simplicity and scalability. In geographic routing [START_REF] Mauve | A survey on position-based routing in mobile ad hoc net-works[END_REF]], a node simply forwards a packet to be closer to the sink in terms of geometric distance, so this position based routing named Greedy Routing is simple and scalable since it does not require global topology information. However, due to the lack of global information, greedy routing can fail at voids or obstacle where there is no neighbor node that is closer to destination than this node. This is well known as the local maximum (or minimum) problem. This can be easily avoid when all node locations are known. Mobility of one or more nodes creates a fully dynamic topology and routing strategy must be adapted. Various geographic information based routing could be used. [START_REF] Zhao | VADD: Vehicle-assisted data delivery in vehicular ad hoc networks[END_REF] improve geo-routing using predicable vehicle mobility, which is limited by the traffic pattern and road layout. This concept of using expected next location in routing is interesting for our research.

In the same way [START_REF] Bertaux | Architecture réseau pour véhicule de transport en commun communiquant[END_REF] uses localization of the vehicle-node in order to anticipate the handover and thus avoid disruption of the information flow due to time to set protocols.

Mobile WSN localization

The location of a node is often required by the application in order to associate this location with value(s) measured by the sensor(s). It is also useful in geographic routing when this type of routing used. Localization of a WSN node can be done by various strategies [START_REF] Li | Survey of localization techniques in wireless sensor networks[END_REF]].

These can be divided regarding the information processing into centralized or distributed methods, regarding the type of used signal into range-based or range-free methods and regarding the obtained position into absolute or relative methods. GPS is the most popular tool for outdoor localization but in the case of WSN mobile nodes, an efficient and accurate method has been proposed by [START_REF] Chen | Accurate and efficient node localization for mobile sensor networks[END_REF]]. However a proven manner for mobiles is the use proprioceptive sensors such as low cost accelerometers [START_REF] Hsu | An accelerometer based approach for indoor localization[END_REF]], gyroscopes, gyrometers, magnetometers, inclinometers or Doppler radar. We used one or other localization method up to now.

A NEW GEOGRAPHIC DYNAMIC ROUTING

Problem statement

Let us consider a wireless sensor network consisting of N mobile nodes equally distributed around a circle of radius R, as shown in figure 1. All the nodes move together around the circle so that they always remain at the same distance from each other, whatever the variation of the angular speed with time. A static sink node, which collects sensed data from the mobile nodes for the end user, is located slightly outside of the circle. Sensor nodes use radio-frequency communications to transmit data packets to the sink. The radio transmission range depends on the transmission output power. If the transmission output power is large enough to that sensor nodes can always reach the sink in one hop only, then the routing problem does not occur. However, using a single-hop communication model can be impossible when R is too large compared to the maximum transmission distance, or even undesirable especially for energy conservation purposes when sender nodes are battery-powered (see [START_REF] Monks | A study of the energy saving and capacity improvement potential of power control in multi-hop wireless networks[END_REF]). In the case where a multi-hop communication model is adopted, it suffices that sensor nodes are able to transmit data packets to their two nearest neighbors for that a route to the sink exists. Let d be the distance between two consecutive nodes around the circle. We have:

d = 2R sin π N (1)
Consequently, sensor nodes need to have a radio transmission range slightly higher than d for the wireless sensor network works normally. Under this condition, indeed, sensor nodes which are far away from the sink have two possible routes to reach the sink using multi-hop communications: a forward route following the rotation direction of mobile nodes and a backward route following the reverse direction.

We assume that each mobile sensor node is able to monitor the dynamic variations of its angular position and speed in real time (see section 2.3 which explains how this hypothesis can be tested). Therefore, the sensor-to-sink routing problem can be addressed locally, at the sensor node level. For the source sensor node, the problem is to find the best route to the sink as a function of its current position and speed, or, in other words, how to alternate between the single-hop route, the backward route and the forward route.

Route selection algorithm

Let O be the center of the circle and let S be the location of the static sink. Location of mobile nodes can be conveniently described in the polar coordinate system with pole O and polar axis -→ OS. Consider any mobile node, for example the one with polar angle θ which is labelled A in figure 2. Node A has two neighbor nodes that are candidates for next hop selection, node C in the rotation direction (forward route) and node B in the reverse direction (backward route). Let T be the one-hop transmission delay to send a data packet. If node A at position θ selects the forward route, then node C will receive the packet at position θ + 2π N + θT due to mobility of sensor nodes, where N is the number of sensor nodes and θ is the angular speed. In the same way, if the reverse route is selected, node B will receive the packet at position θ -2π N + θT . Generally speaking, after k hops, the data packet reaches the node located at θ + k( 2π N + θT ) and θ -k( 2π N -θT ) by forward and reverse routes, respectively. Of course, the route having the lowest end-to-end delay, and hence the shortest hop count, is preferable. The choice between the forward route and the reverse route depends on the position θ of the source node.

The number of hops, k f , to reach the mobile node closest to sink by the forward route is given by: lim

k f ∈Z+ [ θ + k( 2π N + θT ) ] = 2π (2)
Similarly for the reverse route, we have:

lim kr∈Z+ [ θ -k( 2π N -θT ) ] = 0 (3)
From ( 2) and ( 3), the equilibrium point, i.e., the position of source node where k f = k r , is given by:

θ 0 = π - N θT 2 (4)
As a result, in the general case, the source node will select the forward route when its current position is higher or equal to θ 0 , and the reverse route otherwise. If no neighbor node is closer to sink than itself, then the source node directly transmits data packet to the sink. This particular case appears when θ > (2N -1)π

N -T θ 2 or when θ < π N -T θ 2 .
Let us note that T , the one-hop transmission delay, depends on the nominal data rate of the radio transceiver and the protocol used at the MAC layer. The value of T needs to be known a priori and registered as an input parameter of our geographic routing algorithm.

PERFORMANCE EVALUATION

Method

The performance of our routing protocol is evaluated by simulations using WSNet network simulator [START_REF] Chelius | WSNet simulator[END_REF]] and considering the scenario of the Singapore Flyer. Consequently, 28 mobile nodes and a static sink are deployed over an area of size 170 × 170. The center of the circle is at coordinates (85, 85) and its radius is 75. The static sink is located at (85, 5). Mobile nodes are equally distributed around the circle: The initial coordinates of node 1 are (85, 10), those of node 2 are (101.69, 11.88), those of node 3 are (117.54, 17.43), and so on. Mobile nodes move around the circle at speed of 0.0032 rad/s. The radio model adopted in simulations complies with the RF CC2420 transceiver which is used on real wireless sensing devices such as the well-known MicaZ, TelosB and Imote2. It operates in the 2.4 GHz frequency band using O-QPSK modulation and provides a nominal data rate of 250 kbps. At the MAC layer, the famous B-MAC protocol [START_REF] Polastre | Versatile low power media access for wireless sensor networks[END_REF]) is selected. B-MAC can work in low-power listening (LPL) mode, i.e., node periodically sleeps, wakes up to check briefly for channel activity, and then returns to sleep if no activity is detected. B-MAC is configurable with some parameters which, last but not least, can be dynamically adjusted. Especially, the check interval (i.e., the period between consecutive wakes-up) is one of the parameters. Moreover, LPL mode can be enabled or disabled. Obviously, the one-hop packet delay (i.e., the input parameter T in our routing algorithm) depends on the configuration parameters of B-MAC. Especially enabling LPL mode has an influence on the packet delay, and packet delay increases as the check interval increases.

Simulation parameters are given table 1. Simulation duration is 4000 seconds, this time being equivalent to two rounds of the wheel. Simulation results are shown for both cases where LPL mode is either enabled and disabled.

When LPL mode is enabled, the check interval is 10 ms. For comparison purposes, the same simulations are done again with a flooding protocol replacing our geographic protocol. Two performance criteria, namely, end-to-end packet delay and packet loss rate, are used for the comparison.

Results

A synthetic view of various simulation results is given in figures 3 and 4 and in tables 2 and 3. As explained before the two main criteria useful for monitoring application are delay and packet loss. Indeed bandwidth could be considered but it is directly obtained by traffic load when packet loss begins to increase.

Figure 3 shows the variation of end-to-end delay when the wheel turns two rounds. This is the delay of a packet sent by the node 1 and received by the sink. Node 1 is the closest to the sink at t = 0. All 28 nodes send packets during simulation time but only those sent by node 1 are tracked. In order to avoid additional delay due to bandwidth bottleneck CBR (Constant Bit rate) application is set with a period of 30s. This figures shows that delay is directly related to the number of hops between source and destination. So it is a minimum when the sink is inside the radio propagation radius that is for t about 0 or 2000 or 4000s. It is a maximum when the node is in farthest location from the sink that means packets make 15 hops. In order to quantify the influence of duty cycle at MAC layer performance measurements are made in two cases: always-on (LPL disabled) and sleep-wakeup (LPL enabled) with a check interval of 10ms. Sleep-wakeup mode increases dramatically the delay as it was expected. Table 2 gives a statistical analysis of this delay with minimum, average and maximum values. Comparisons are not only with both B-MAC configurations but are extended to two routing algorithms: the simple flooding and our geographic routing. In all cases delay is a little lower with flooding than with our protocol. This can be explained by the fact that flooding is simpler than geographic routing and does not require computing before sending any packet. But the difference is not great and cannot be compared with that resulting from the use of the MAC protocol with LPL mode. The maximum remains very low, about 23 ms, with always-on MAC protocol. Packet loss rate as a function of packet sending period is plotted on figure 4 for the two experimented routing protocols. Both protocols are effective when packet interval is greater than 10 s and packet loss rate is virtually nil (less than 1%). But when the period decreases below 10 s packet loss increases sharply in the case of flooding routing. The network works until a period of 2 s in the case of our routing protocol virtually without packet loss. This is an important result that the threshold of the sampling period is significantly lower when using our geographic routing protocol. It could be necessary to monitor some specific parameter with a very short sampling period during a test. Table 3 gives the minimum packet interval or sampling period with the always-on MAC protocol and with 3 values of the check interval with LPL mode. Our protocol is always used. These values are lower limits that guarantee a packet loss rate less than 1%. That means they can be used in monitoring application. It can be deduced from the values in table 3 that LPL mode is not suitable for short sampling period.

CONCLUSION

In this paper we have proposed a wireless sensor network architecture for monitoring large physical system in cyclic mobility. We have especially studied its deployment on a big wheel. The performance evaluation was obtained by simulation. The experimental results shown that the lower period of data acquisition on the 28 sensor nodes with low packet loss ratio is short enough to monitor very dynamics physical phenomena. The proposed georouting protocol shortens the minimum sampling period in comparison with a general use flooding routing protocol.

It could be noticed that the architecture chosen to test our physical system is very slow. In such a case the speed parameter has no influence in the path selection algorithm which then depends only on the location. Other result is the limited packet transfer delay that matches with online monitoring. Therefore this wireless system is not only easy and inexpensive to deploy but also delivers a high quality of service.

Some monitoring applications require that samples must be accurately time-stamped in order to correlate readings from different sensors. Sensor network time synchronisation schemes like this proposed by [START_REF] Ganeriwal | Timing-sync protocol for sensor networks[END_REF] are sufficient to solve this problem. We have to implement and test such synchronization.

More ambitious future works are to study other large physical systems moving on non circular trajectories. We plan to adapt the architecture to systems such as chairlifts or material ropeways which move on two parallel linear paths in opposite directions. We are also looking for partners in order to experiment our architecture on a real working system.

Fig. 1 .

 1 Fig. 1. Topology of the wireless sensor network with N mobile sensor nodes and one static sink.

Fig. 2 .

 2 Fig. 2. Location of mobile sensor nodes in polar coordinates.

Fig. 3 .

 3 Fig. 3. End-to-end packet delay as a function of time for node 1.

Fig. 4 .

 4 Fig. 4. Packet loss rate as a function of packet sending period.

Table 1 .

 1 Simulation parameters

		Parameter	Value
	Radio models: Propagation range	20 meters
		Modulation	O-QPSK
		Data rate	250 kbps
		Interferences	orthogonal
		Antenna	Omnidirectional
	Protocols:	MAC layer	B-MAC
		Routing layer	Our, flooding
		Application layer	CBR
			(20 bytes of payload)

Table 2 .

 2 Delay comparison (in ms)

		Min.	Avg.	Max.
	Our (always-on)	1.536	12.036	23.041
	Flooding (always-on)	1.408	10.027	19.713
	Our (Sleep-Wakeup)	12.288 97.017 186.040
	Flooding (Sleep-Wakeup) 12.160 87.406 170.241

Table 3 .

 3 Minimum packet interval for CBR application

	LPL mode	Our
	disabled (radio always-on)	175 ms
	enabled (check interval of 10 ms)	1.350 s
	enabled (check interval of 20 ms)	2.450 s
	enabled (check interval of 50 ms)	5.850 s


Thomas Capelle
email: thomas.capelle@inria.fr

Peter Sturm
email: peter.sturm@inria.fr

Arthur Vidard
email: arthur.vidard@imag.fr

Formulating LUTI calibration as an optimisation problem: example of Tranus shadow-price estimation

Introduction

Integrated land use and transport modelling has attracted the attention of researchers since 1960 [START_REF] Wegener | Overview of Land-Use Transport Models[END_REF]. Over the years, a large number of models have come into existence. It is well known that integration of land use and transport models creates a complex nonlinear system, which evolves in different scales. Analysing these complex systems is typically a hard problem, especially in the presence of uncertainty, whose effects may be difficult to assess. In such cases, calibration plays a central role, as it helps us determine optimal parameters. Calibrating this type of models is a process that requires many steps. At first, the data needed is not always readily available, but even if so, finding the set of parameters that replicates better the data is not simple. We propose an optimisation framework for the calibration of Tranus, particularly for obtaining a good estimation of the shadow-price variables.

Methodology

Tranus Description. Tranus [3] provides a generic framework to model land use and transportation in an integrated manner, both in urban and regional levels. The region of interest is divided into economic sectors and spatial zones. Then Tranus combines two modules: the land use and activity module which simulates a spatial economic system by assessing the activity locations and economic sector interactions (LCAL); and a transportation module, which estimates the use of the transport network and the associated disutility (LOC).

The land use and activity module estimates the productions and the consumptions for zones at a given period, and the demands of flows that this activity generates. These demands of flows are then fed to the transportation module. In this way the movements of people or freight are explained as the results of the economic and spatial interaction between activities, the transport system and the real estate market. Then from the transportation network, once the corresponding trips are generated, travel flows are allocated to the network according to travel demand. In turn, the accessibility that results from the transport system influences the location and interaction between activities through transport dis-utilities, also affecting land rent.

The two modules in the system use discrete choice logit models [START_REF] Mcfadden | Mixed MNL Models for Discrete Reponse[END_REF], linked together in a consistent way. This includes activity-location, land-choice, and multi-modal path choice and assignment. The modules are then run iteratively, such that production and consumption demands for each zone are met and equilibrium is achieved. A detail description of the equations of Tranus land use and transportation modules can be found in [3].

The problem. In Tranus, the calibration of the land-use and acitivity module is done in a program called LCAL (Land-use CALibration). The objective of the calibration, is to have a model that can replicate the data of the base year. For achieving this, the model has a set of parameters that can be adjusted (elasticities, attractors, discrete choice models constants, etc.). The job of the modeler, is to adjust the input parameters, to make the model fit the base year production data X 0 as close as possible. To be sure that the model can replicate the base year's production, a correction term is added inside LCAL. This variable is added to the utility function and acts as a correcting price, called the shadowprice. One shadow-price is added for each production. What LCAL does in the actual implementation of Tranus, is for any given set of parameters, it will optimise the shadow-prices to replicate the base year X 0 for a fixed set of input parameters. This is done via an iterative method that tries to find a local optimum just exploiting the economic notion of the price and shadow-price variable, there is an explicit cost function. This can lead to local optima, or not even find a solution.

Instead, we propose a framework, with an actual cost function to optimise and a clear and explicit description of the problem.

The problem to solve in LCAL, is to adjust the shadow-prices h ∈ | n × | j (n sectors, j zones) to replicate the base year production. The problem can be written as an optimisation problem:

min h f (h) = X(h) -X 0 2 (1)
The function X(h) doesn't have a closed form, because it is the result of a double fixed point problem:

X(h) = Λ 1 (h, p) + ∆ 1 (h, p)X (2)
p(h) = Λ 2 (h, p) + ∆ 2 (h, p) p (3)
Where equation (3) represents the computation of the intermediate variable p (price), and imposes that the prices are in equilibrium. The problem here is that for computing X(h), we need to solve (2)-(3) iteratively until convergence. Instead of imposing that the prices are in equilibrium, we can add it to the objective function, hence we can rewrite the problem as follows:

min h,p f (h, p) = X(h, p) -X 0 2 + p(h) -p 2 (4)
where X is the solution of the linear system (2) and p is the solution of (3), X 0 is the observed induced production. This appears to be a good approach to attain reasonable shadow prices, but the construction of the ∆ 1 and ∆ 2 matrices and the inversion of linear systems (2)-(3) is expensive. That's why we propose another idea, based on the concept of imposing X 0 as fixed point for the linear system. Instead of solving the linear systems, we evaluate them in X 0 . We exploit the idea that that if we could find a pair (h * , p * ) so that f (h * , p * ) = 0, then X(h * , p * ) = X 0 and p(h * ) = p * . We define the functions:

X(h, p) =Λ 1 (h, p) + ∆ 1 (h, p)X 0 (5) p(h, p) =Λ 2 (h, p) + ∆ 2 (h, p) p (6)
Then, our optimisation problem is as follows, min

h,p f (h, p) = X(h, p) -X 0 2 + p(h, p) -p 2 (7)
Test data. It general terms, is difficult to evaluate a calibration. This difficulty is realated to the fact that we dont know the optimal set of parameters. To be able to test our optimisation scheme, we constructed a set of input data (X 0 , parameters) that has a perfect calibration (up to round-off error).

In other words, for an specified shadow-price h, the calculated X(h) = X 0 . This is a practical for two reasons. First, we are sure to know the optima, so we can make plots of the function around the optimal value. Secondly, asessing the robustness of the optimisation algorithm, testing different starting points and checking if the convergence to the optimal value is attained. To generate the test files, we start with a initial h (ideally a vector of zeros), and we compute iteratively the production (2) and price (3) to attain convergence between successive iterations, then we replace the input data X 0 by the estimated production X(h), thus generating a perfect data set.

Results

In terms of optimisation, the functions ∆ and Λ are differentiable, even class C ∞ , so we develop an optimisation framework using the optimisation algorithm presented in [START_REF] Levenberg | A Method for the Solution of Certain Non-Linear Problemas in Least Squares[END_REF]. The Jacobian of the function is computed analytically.

Example: We tried a small example, based in the Example C from the Tranus website (http://www. tranus.com/tranus-english/download-install). We generated the perfect files, with shadow-prices h n i = 0 for each sector n and zone i. This example, has 3 zones and 5 sectors. We did some cuts near the optimal value, and plotted the function f (h, p). In the Figure 1, the graph set is generated evaluating the function f for shadow-prices h of the first sector and the three different zones. As expected, the cost function is zero at h = 0, and increases its value when we get far away from the optima. The function appears to be locally convex near the optimal value (at least for these 3 parameters) What is really interesting, is to have a plot of the function f near the optimum, for a fixed sector and zone, and see how the function behaves with h and p. If we consider sector 1 and zone 2, we can

2 ) near the optimal value (depicted in red) plot f near the optimal value h 1 1 = 0, p 1 1 = 2.676 as shown in Figure 2. Here we can observe that as the shadow-price gets larger the value of f increases up to a plateau state (X 1 1 (h) → 0 and then f (h, p) becomes constant). In the case of the price p, if we move away from the optimal value p = 2.676, the cost increases quadratically.

We tested the robustness of the optimisation scheme with 1000 random initial values, always converging to the same optimum. The prices are in the interval [0, 4], so considering the shadow prices in [-10, 10] is highly representative.

Conclusions

We have a developed a optimisation methodology that give us a partial calibration of Tranus and a clear representation near the optimal value. Secondly, we have proposed a technical contribution to the way equilibrium equations ( 2)-(3) are handled in Tranus, permiting us to decouple the double point fix problem. Finally, the procedure of generating syntethic data files for testing the calibration methodology is a practical way to test and check the performance of a general calibration scheme.

Figure 1 :Figure 2 :

 12 Figure 1: Graph of slices function f along h 1 i , for each zone i ∈ [1, 2, 3]


Antoine Gourlay

Erwan Bocher

Gwendall Petit

Towards a WPS platform dedicated to an urban knowledge infrastructure

Towards a WPS platform dedicated to an urban knowledge infrastructure 3 IRSTV researchers and engineers -noise mapping, flooding modeling, atmospheric pollutant dispersion, etc. -and by creating a geospatial knowledge repository to study cities. These processes would then be accessible to the federation's partners: local authorities, agencies, citizens, private organizations and companies.

Introduction

The Institute for Research on Urban Sciences and Techniques (IRSTV) is a Research Federation founded by CNRS (French National Center for Scientific Research) and France's Ministry of Research brings together fifteen research laboratories in order to lead an interdisciplinary research in the field of urban environment and especially in the development of knowledge, models, tools and methodologies for the management, the design, the decision-making and the evaluation of city planning policies. With this aim in view, the IRSTV gathers about two hundreds researchers, engineers and administrative staff.

The research federation structures projects within three main fields:

-Interdisciplinary urban observation system (urban remote detection) -Integrated urban environment modeling (urban micro-climatology, environmental noise analysis, GIS and urban data modeling) -Governance, design and sustainable management of cities. This interdisciplinarity created a discrepancy in the various usages of geographic information: data storage and modeling, tools used to manage it and processing workflows used were very heterogeneous [START_REF] Bocher | Projet d'infrastructure de données spatiales[END_REF]- [START_REF] Bocher | GDMS: an abstraction layer to enhance Spatial Data Infrastructures usability[END_REF]. The diversity of GIS software, of data storage formats, of ways to describe the data lead to a compartmentalization of urban geographic knowledge, which was going against the purpose of the IRSTV: to develop an integrated view of all urban data, methods, tools and protagonists that contribute to sustainable city planning and management.

In order to resolve this paradox, and to strengthen the federative attitude within IRSTV, a Geographic Information System (GIS) for urban modeling and city planning was outlined in the regional project MEIGEVille (Integrated environment modeling and sustainable urban planning) [START_REF] Hegron | MEIGEVille, Modélisation environmentale intégrée et gestion durable de la ville[END_REF]. The purpose of this Towards a WPS platform dedicated to an urban knowledge infrastructure 2 GIS was to lay out the bases for an urban knowledge gathering tool, analysis methods and management techniques [START_REF] Hegron | MEIGEVille, Modélisation environmentale intégrée et gestion durable de la ville[END_REF].

A GIS team ("Atelier SIG") is thus created in December 2006 to build and maintain two platforms:

-a Spatial Data Infrastructure (SDI) -a community-based Geographic Information System (GIS)

Using several documents describing good practices for geographic data exchange [START_REF] Clinton | Coordinating geographic data acquisition and access: The National Spatial Data Infrastructure[END_REF][9] [START_REF] Nebert | Developing Spatial Data Infrastructures: the SDI Cookbook (v2.0)[END_REF], the Atelier SIG built an interoperable architecture consisting of (Erreur : source de la référence non trouvée) [START_REF] Bocher | Projet d'infrastructure de données spatiales[END_REF]- [START_REF] Bocher | Urbsat : outil d'extraction d'informations géographiques et de production d'indicateurs[END_REF]:

-a data repository -a server application providing geographic data services that use OGC standards (WMS, WFS, WCS…) -a web portal for displaying, searching and browsing geographic data -a GIS software called "OrbisGIS" for data displaying, querying, processing and sharing

The project presented in this paper aims at enforcing the INSPIRE Directive [START_REF]Directive of the European Parliament and of the Council establishing an Infrastructure for Spatial Information in the European Community (INSPIRE)[END_REF] on a local scale by extending the capabilities of this SDI in order to providing on-demand geospatial web processing using the WPS Standard. This would enhance the current uses of the SDI by pooling all the processes built by

The WPS Standard

The OGC Web Processing Services (WPS) Interface Standard [START_REF][END_REF] provides a standardized mechanism to invoke geospatial processing services as web services. It defines in detail how inputs and outputs are to be handled, as well as how clients can request the execution of a process.

The specification requires that Extensible Markup Language (XML) be used for all communication. XML documents are made up of elements that act as container for data. These elements can be nested and can be given attributes that describe them. The main benefit of XML is that it is readable by humans, which makes designing and debugging much easier.

A WPS service can expose any number of processes of any kind as web services. WPS defines three operations to communicate with a server:

-GetCapabilities return information on the WPS service, its capabilities and its available processes -DescribeProcess returns a full description of a process, including inputs and output -ExecuteProcess invokes a process and returns its output(s)

These actions are accessible through commonly used protocols such as HTTP GET and POST, as well as SOAP and WSDL.

From a user point of view, GetCapabilities gives access to the list of available processes. Then DescribeProcess is called for those processes to retrieve their descriptions. For every process, a data structure is returned containing the process name, an optional description, a version number, the details about its inputs and outputs, and two other parameters:

-If data output(s) from the process can be requested to be stored by the WPS server as web-accessible resources Once a software client has called DescribeProcess and has gathered from the user the necessary inputs to call the process, it can build the necessary parameters for ExecuteProcess. Normally, the response is returned only after process execution is completed. If the process succeeds, the response document contains either the outputs (Erreur : source de la référence non trouvée), or links to the outputs stored in webaccessible locations (Erreur : source de la référence non trouvée). In the special case of a single complex output, the output can be returned directly instead of the XML response document. However, if instructed the server can return a URL from with the response can later be retrieved. This response document can be regularly updated to provide feedback to the caller through a process status (Erreur : source de la référence non trouvée).

The process status can be "ProcessAccepted", indicating that it has be queued and will be processed; "ProcessStarted", indicating that the process is currently running; "ProcessSucceeded", signifying that the processing is over; or "ProcessFailed", meaning that a problem occurred during the execution. In the latter case, an error code is also returned in the response document.

Existing tools supporting WPS

We distinguish two different categories of tools supporting the WPS Standard: web servers for data processing, and GIS desktop software & libraries that are consumers of such processes.

A. Web Servers for data processing

These servers expose processing functionality that is usually found on desktop GIS software -for example data transformation, analysis or re-projection.

One implementation of the WPS Standard is 52°North WPS (52north.org), with native support for GRASS (grass.osgeo.org) and SEXTANTE (sextantegis.com). Another solution is part of degree 3 (wiki.degegree.org), an project for the Java platform. The PyWPS software (pywps.wald.intevation.org) also implements the WPS Standard, with native support for GRASS GIS. Finally, the ZOO Project (zoo-project.org) also provides an implementation.

All software platforms listed above provide a full WPS 1.0.0 compliant implementation, and most of them provide connectors for several other GIS like Figure 8 -uDig in JGrass Towards a WPS platform dedicated to an urban knowledge infrastructure 7 GRASS GIS. See [START_REF] Steiniger | An Overview on Current Free and Open Source Desktop GIS Developments[END_REF] and [START_REF] Steiniger | Free and Open Source GIS Sorfware for Building a Spatial Data Infrastructure[END_REF] for a detailed comparison.

B. GIS Desktop software & libraries

Few Desktop GIS provide native support for WPS, but plugins are available. The GIS gvSIG (gvsig.org), GRASS GIS (grass.osgeo.org), OpenJUMP (openjump.org) and uDig (udig.refractions.net) all support querying WPS processes via plugins (Figures 5678). The libraries SEXTANTE (sextantegis.com) and OpenLayers (openlayers.org) can also use or display a WPS process result. See [START_REF] Steiniger | An Overview on Current Free and Open Source Desktop GIS Developments[END_REF] and [START_REF] Steiniger | Free and Open Source GIS Sorfware for Building a Spatial Data Infrastructure[END_REF] for a detailed comparison.

Integration of the WPS standard into IRSTV's Spatial Data Infrastructure

A. IRSTV's SDI

A Spatial Data Infrastructure (SDI) is a group of tools interconnected using international standards defined by the Open Geospatial Consortium (OGC). The SDI currently available at IRSTV consists of five modules (Erreur : source de la référence non trouvée):

-A spatial relational database management system (RDBMS) housing national geographic data (databases from the French National Mapping Institute (IGN)) as well as data produced by IRSTV (statistics, processed satellite images, simulation results…). Geographic data is described, structured and queried using the Simple Feature SQL norm [START_REF]OpenGIS Implementation Specification for Geographic Information -Simple feature access -Part 2: SQL option, v1.2.1[END_REF].

-A server that exposes geo-services. This application makes geographic data available using the Web Map Services (WMS) [START_REF]OpenGIS Web Map Service (WMS) Implementation Specification, v1.3.0[END_REF] and Web Feature Services (WFS) norms [START_REF]OpenGIS Web Feature Service 2.0 Interface Standard, v2.0[END_REF]. -A geographic data catalog (Geocatalog). It contains geographic metadata using the norm ISO 19115, exposing information such as temporal and spatial extend of the data, its origin, etc. -A cartographic website that enables users to query the geocatalog, as well as remote catalogs as long as they respect the OGC Catalogue Service (CSW) standard. Based on A Geographic Information System (GIS) called "OrbisGIS". Completely developed at IRSTV by the Atelier SIG, OrbisGIS has all the functionalities that make a GIS: data reading, writing and edition, map creation… It contains an abstraction layer called "Gdms" that encapsulate both alphanumeric and geographic data (geometric and raster) into a relational tabular view as would any Relational Database Management System such as PostgreSQL/PostGIS. This layer allows the use of standard SQL, with support for the Simple Feature SQL norm, as a querying language for advanced processing. Thus OrbisGIS is a user interface for this data management layer.

All the modules of the SDI are based on Open Source components and are themselves Open Source. This choice originated in the highly research-specific needs of the federation, the desire to invest into people instead of software and the will to reuse and share with other actors faces with the same challenges. This is also in line with the public utility duties of a state research federation like IRSTV.

B. The integration of WPS

Towards a WPS platform dedicated to an urban knowledge infrastructure 9 WPS aims at providing processes as on-demand web services. The point is to make it become the key to an SDI where:

-Processes can be pooled among a large user base, and their results cached and shared. -The content of a process is exactly the same whether it runs in OrbisGIS or on a server as a WPS process. -The procedure of submitting a script to become a WPS process is simple and does not requires any extra knowledge about WPS.

This would add to IRSTV's Spatial Data Infrastructure a new kind of web services: geoprocessing services. The different research laboratories could then provide their specialized algorithms as a service, as well as the data they produce, within a coherent multidisciplinary urban knowledge platform (Erreur : source de la référence non trouvée).

Gdms -the data management engine used in OrbisGIS -uses standard SQL syntax [START_REF]Information technology -Database languages -SQL[END_REF] for its treatments, with support for OGC Spatial functions [START_REF]OpenGIS Implementation Specification for Geographic Information -Simple feature access -Part 2: SQL option, v1.2.1[END_REF] and other functions provided by researchers for the specific needs of the federation. These additional functions cover fields from hydrography to triangulation. OrbisGIS can be used to both process and display data -inputs and results -at the same time (Erreur : source de la référence non trouvée). Such a script should be able to run as a WPS service as long as some metadata about inputs and outputs is provided, as well as a name and a description.

The objective is to use the specific data management layer Gdms and its specialized functions to build a WPS-accessible processing server. This server would be able to reuse all processes currently working in OrbisGIS that are available within IRSTV, as well as enabling the creation of new scripts that reuse specific functions in fields like hydrology (Erreur : source de la référence non trouvée).

To that end, we designed a SQL descriptive metalanguage that a user can supply to WPS-enable a script without altering its content. The metadata is present at the beginning of the script, enclosed in SQL comments so as not to change the behavior of the script.

A very simple script to filter only valid features from an input data set would look like this:

--@script --@identifier FilterValidData --@title Filter --@abstract Filters the content of a --spatial table to keep only valid geometries --@/abstract --@input TABLE myTable

( the_geom GEOMETRY ) --@output TABLE result --@/script CREATE TABLE result AS SELECT * FROM myTable WHERE ST_IsValid(the_geom)
This script can run in OrbisGIS or in a spatial DBMS as long as the table "myTable" is present, and can run as a WPS service if a spatial table -whatever format it is in -is provided, as long as it has a geometry field, that will be called "the_geom" in the script -the default geometry name in OrbisGIS. The inclusion of WPS in the SDI is implemented as a module-based architecture on an application server (Erreur : source de la référence non trouvée).

Implementation

A. Server architecture

-The deegree 3 WPS frontend is used to handle client request. This particular implementation was chosen over 52° When a GIS calls GetCapabilities, the following happens:

1. Deegree asks OrbisGIS Server the list of all available processes. 2. OrbisGIS Server queries its catalog and returns a list of script with some details. 3. Deegree sends the response, with the list of processes and its own WPS-specific data.

Then when a GIS calls DescribeProcess, the following happens:

1. Deegree asks OrbisGIS Server for a specific process. 2. OrbisGIS queries its catalog and returns the details for the script. 3. Deegree send the response, with the detailed information required by the WPS Standard.

Finally when a GIS calls ExecuteProcess with a specific process, the following happens:

1. Deegree asks OrbisGIS Server for the process. 2. OrbisGIS Server finds the script in its catalog and wraps it in a deegree Process. 3. Deegree starts the process. 4. The process gets all inputs and gives them to Gdms Server to be registered as DataSource inputs for the current execution. 5. Gdms Server parses, optimizes and validates the script against the available DataSource. 6. Gdms Server runs the script.

Towards a WPS platform dedicated to an urban knowledge infrastructure 14 7. OrbisGIS Server saves the DataSource registered as outputs, discard the whole session (all other data or state is lost) and returns to deegree the outputs. 8. Deegree sends the response to the GIS. This whole process relies on the meta-language being added by the user and interpreted by the server.

B. The meta-language

The meta-language is added by the client to the SQL script before sending, but is parsed and used at a different level:

-The SQL script is parsed in Gdms Server and executed by the built-in SQL engine. -The meta-language is parsed in OrbisGIS Server for extracting metadata about the script, and is subsequently completely ignored.

The meta-language was design with a very simple syntax. It requires that:

-All lines have to start with two dashes (--) followed by a space -A key-value pair is given inline as:

--@key some value -Or with this syntax if is spans over multiple lines:

--@key a long value that spans over --multiple lines --@/key -The whole metadata is a key-value pair whose key is named "script":

--@script --@key value --@key2 value2 --...

--@/script

Expected For example, the following defines an input that contains geometry information:

--@input TABLE myTable ( the_geom GEOMETRY )

OrbisGIS Server will check the input against that schema before running the content of the script. -LITERAL means this input is a constant value provided by the caller. The value can be reference in the script directly using its name.

--... --@input LITERAL buffer_size INT --@/script CREATE TABLE result AS SELECT ST_BUFFER(the_geom, buffer_size) FROM myTable;

The "param" after the name is a valid SQL Type. -BOUNDINGBOX means this input is a Bounding Box. It will be converted for use by the script to a geometry value with the shape of the box.

Outputs can be formatted the same way as inputs:

--@output TYPE name params

The first value after "INPUT" is either

C. Input and output formats

The purpose of the meta-language is to abstract the description of the process from the technical details of WPS. Also, the process should only work with data, not with specific data structures. Thus the server provides the necessary abstraction through the "Tables" accessible in the script, so that all data structures supported by the server are supported for any process.

The data structures currently supported by Gdms Server are: Other formats might be added in the future to Gdms, thus making them accessible to be input/output in a WPS process.

D. Pushing scripts to the server

OrbisGIS Server comes with a web frontend. The current prototype only allows the upload of a script through a basic web form. A complete website will be added with support for:

-browsing of the script catalog and catalog management, -authentication and user management, Towards a WPS platform dedicated to an urban knowledge infrastructure 17 -a web service for submitting scripts directly from OrbisGIS.

Applications

As outlined by the WPS Specification, a WPS Service can be used for any kind of geospatial processing, from a simple buffer computation to a specialized urban area density diachronic analysis process.

Three examples are provided below:

-a simple geometry buffer computation -an example from the hydrology toolbox at the IRSTV: drainage network extraction -an example from the field of urban area analysis: urban area density diachronic analysis

A. Simple script: geometry buffer

A common and simple use of a geo-processing service would be to compute a fixed-sized buffer around a set of geometries.

Such an operation would look like that in OrbisGIS as a SQL Script:

CREATE TABLE buffered AS SELECT ST_BUFFER(the_geom, 10) FROM input;
Provided there is a table "input" with geometry data registered in the current set of sources in OrbisGIS, this short script will produce the expected buffered geometries.

WPS-enabling this script is very easy in this case: --@script --@identifier Buffer --@title Geometry buffer --@abstract Computes a buffer --around all given geometries --@/abstract --@input TABLE ( the_geom GEOMETRY ) --@input LITERAL buffer_size INT --@output TABLE buffered --@/script CREATE TABLE buffered AS SELECT ST_BUFFER(the_geom, buffer_size) FROM input;

The metadata describes 2 input values:

-a table with geometries -a literal integer value specifying the buffer size

The metadata also defines a single output table that will contain the buffered geometries. The script is now ready to be uploaded to the webserver.

Once the user has connected to the web site and uploaded the script (or, in the future, deployed the script from OrbisGIS), any WPS client software can Towards a WPS platform dedicated to an urban knowledge infrastructure 18 browse the WPS Service, find the new process and execute it.

B. Hydrology: drainage network extraction

Geographic Information Systems are used in conjunction with hydrology to map physical water behavior over the territory. Often a raster Digital Elevation Model (DEM) is used to produce morphologic and hydrologic indicators. Coupled with data about land usage and city planning, they can then give answers as to the hydrological behavior of a particular area.

In this example a hydrological drainage network and a watershed are extracted from a DEM (raster) called "dem_chezine". The process follows these steps:

-the DEM is filtered to remove any depression (sinks are harmful to the process) [START_REF] Planchon | A fast, simple and versatile algorithm to fill the depressions of digital elevation models[END_REF] --Filter to fill sinks CREATE TABLE filled AS SELECT ST_FillSinks(raster, 0.1) AS raster FROM dem_chezine;

-then a topographic direction grid is produced using the D8 algorithm --Compute the direction grid CREATE table dir AS SELECT ST_D8Direction(raster) AS raster FROM filled;

-then an accumulation grid is extracted as described in [START_REF] Ocallaghan | The extraction of drainage networks from digital elevation data[END_REF] --Compute the accumulation grid CREATE table acc AS SELECT ST_D8Accumulation(raster) AS raster FROM dir;

-a theoretical drainage network is then produced by connecting pixels whose accumulation value is above 1500 and sorting them using Strahler's ordination method.

--Ordination of Strahler, allowing prioritization of the streams CREATE TABLE strahler AS SELECT ST_D8StrahlerStreamOrder(d.raster, a.raster, 1500) FROM dir d, acc a; -finally the obtained rivers are converted to geometries and filtered to keep only streams of order lower than 6 (removing the smaller streams).

--Vectorize the stream This produces the expected drainage network. Another hydrology-specific function can then be used to extract the watershed from the network: the function extracts the outlet of the network from "river_low_6" and climbs back through the accumulation grid in order to produce the watershed. CREATE TABLE watershed_chezine AS SELECT ST_D8Watershed(acc.raster, 'river_low_6') FROM acc;

The whole process can be run in OrbisGIS and produce the expected data as shown on Erreur : source de la référence non trouvée.

This script can be converted and exposed as a WPS process with the following metadata: --@script --@identifier DrainageNetwork --@title Drainage network and watershed --extraction --@/title --@abstract Computes the drainage network --and the associated watershed from --a DEM raster --@/abstract --@input TABLE dem_chezine ( raster RASTER ) --@output TABLE river_low_6 --@output TABLE watershed_chezine --@/script This metadata describes a single raster input (the necessary DEM) and two outputs: the drainage network and the watershed. The whole script can then be uploaded to the WPS Server with OrbisGIS Server and run the same way, to produce exactly the same result (Erreur : source de la référence non trouvée).

C. Urban areas: density diachronic analysis

Urban areas are characterized by an evolution of the land from cultivated parcels to urbanized ones. This third example studies the evolution of urbanized land density between 1999 and 2004 on the same study area.

The script processes a follow:

- To be exposed as a WPS process, the inputs and outputs need to be described. The process has only one output "Evol_MOS_urban", but 3 inputs:

-the data set from 1999, called "ocs99_1_chezine" in the script -the data set from 2004, called "ocs04_1_chezine" in the script -the study area called "study_area", which is used to compute the grid.

The resulting metadata is: --@script --@identifier UrbanDensityComparison --@title Urban density diachronic analysis --@abstract Computes the urban density --difference per 200 meter wide cell between --two data sets --@/abstract - The script can then be turned into a WPS process and its result consulting into any WPS-compatible GIS software (Erreur : source de la référence non trouvée).

-@input TABLE ocs99_1_chezine --( the_geom GEOMETRY, C_NIVEAU INT ) --@/input --@input TABLE ocs04_1_chezine --( the_geom GEOMETRY, C_NIVEAU INT ) --@/input --@input TABLE study_area --( the_geom GEOMETRY ) --@/input --@output TABLE Evol_MOS_urban --@/script

Conclusion & Future work A. Future work on IRSTV's SDI

WPS is currently a young technology with a great future: exposing geospatial processes as web services is becoming more and more essential at every scale, from the European INSPIRE Directive to the efforts of local cities to gather and use their spatial data.

The on-going project within IRSTV's Spatial Data Infrastructure is part of this wider scheme. Although still a prototype, it is being fully developed and will be completely integrated into the SDI in August. A WPS client within OrbisGIS will also be developed to complete the system. This will be an opportunity to gather feedback from potential users, and form a task group dedicated to the improvement of the processing and sharing capabilities of the SDI. This platform will then be used on a daily basis within IRSTV, as well as by the federation's partners.

B. A feedback on the WPS Standard

The WPS Standard was released in 2007 with the objective of initiating the standardization of ondemand processing as web services. This is completely in line with the current trend of Internet-based services, "Cloud Computing", and "Software as a Service" (SaaS). However there are still few implementations of the standard and even fewer use case scenarios clearly defined and displayed on the web. This study is also an opportunity to produce real Towards a WPS platform dedicated to an urban knowledge infrastructure 24 use cases with real data in one of the many fields ondemand processing can be applied to.

This experimentation is also an occasion to test a full scale implementation of an on-demand geospatial processing service with the WPS Standard. It outlined a potential improvement on the WPS Standard in its current form.

The standard defines a very abstract model for a web service, with the possibility to add profiles for specific uses. There is however no standardized WPS profile for dealing with geographic data yet.

A consequence of this is that geospatial processes deal with inputs and outputs in term of GML format, Shapefile format, etc. instead of Geometry objects, i.e. they are not provided with an abstraction over the particular implementation details of the file formats. Moreover, there already is a standard defining an abstract way of describing and accessing geometric data: the Simple Feature Access (SFA) Implementation Standard [START_REF]OpenGIS Implementation Standard for Geographic information -Simple feature access -Part 1: Common architecture[END_REF].

A profile extension could define a way for geospatial processes to be described with SFA's Geometry objects, and for the WPS Service to be responsible for reading/writing a set of supported formats. The WPS Service implementation would thus, when a client calls DescribeProcess, fill out the input and output structure parts about the supported file formats.

This would provide an abstract model of a geospatial processing web service that abstracts both sides of the service -the client-to-WPS interface and the WPS-to-process interface -using existing OGC Standards. Such a profile dedicated to geospatial processing could also provide ways to specify geometric requirements/constraints on input data as specified by the Simple Feature standard, instead of file formats constraints. For example, an input could be required to provide Polygon objects, or only objects with a spatial dimension of 2, whatever the actual format is as long as it is supported by the WPS Service.

This abstraction is the ambitious idea behind the use of WPS described in this paper and being implemented within the IRSTV.

C. Personal feedback

As a general engineering student with a specific interest in computer science, this one-year long project was an occasion to learn a lot in several different fields:

-Working for the IRSTV gave me a tour into the wide domains of urban modeling, city planning, hydrology, etc. This was a unique opportunity to use my technical knowledge to Towards a WPS platform dedicated to an urban knowledge infrastructure 25 contribute to much applied sciences, as well as grasping scientific knowledge, both interesting and that might prove useful in the future.

-Studying the WPS Standard was an occasion to dive into the fields of geomatics, geographic information science and geodesy. -Using and creating Open-Source software, and in this case an Open-Source computing platform that would benefit many actors, was a perfect example of building, sharing and contributing to a community-based collaborative platform. Being part of a young and dynamic team with people with different backgrounds also contributed to this collaborative approach to innovation and science.

-Working on this implementation was a wonderful opportunity to do heavy software engineering using state-of-the-art technologies like Scala and widely recognized industry standards like Java or SQL.

Figure 1 -

 1 Figure 1 -IRSTV Spatial Data Infrastructure (SDI) and GIS

Figure 2 -

 2 Figure 2 -ExecuteRequest with ouputs in the response

Figure 3 -

 3 Figure 3 -ExecuteResponse with outputs on a web-accessible server

Figure 4 -

 4 Figure 4 -ExecuteResponse with support for status report during execution

Figure 7 -

 7 Figure 7 -OpenJUMP

Figure 9 -

 9 Figure 9 -IRSTV's current Spatial Data Infrastructure

Figure 10 -

 10 Figure 10 -IRSTV's SDI with geoprocessing services

Figure 11 -

 11 Figure 11 -OrbisGIS with data and a SQL Script

Figure 12 -

 12 Figure 12 -Functions from OrbisGIS available for use in scripts on the WPS Server

Figure 13 -

 13 Figure 13 -Architecture of the proposed implementation

-

  AutoCAD Drawing Exchange Format (DXF)input only -Character Separated Values (CSV) -Esri Shapefile (SHP, SHX, DBF; zipped) -dBase file (DBF) -Gdms native format (GDMS) -MapInfo Interchange Format (MIF/MID) -input only -Solene geometry & attributes (CIR, VAL)input only -Virtual Reality Modeling Language (VRML)input only -TIFF, JPG, or PNG with ArcView world file -Esri ASCII Grid (ASC)

Figure 14 -

 14 Figure 14 -drainage network and watershed extraction

Figure 15 -

 15 Figure 15 -outputs of the WPS process (stacked)

Figure 16 -

 16 Figure 16 -output grid with high densities in red and low ones in yellow

  Filter out streams of order > 6. CREATE TABLE rivers_low_6 AS SELECT * FROM allrivers WHERE gid < 6;

	Towards a WPS platform dedicated to an urban
	knowledge infrastructure 19
	CREATE TABLE allrivers AS SELECT ST_VectorizeLine() FROM strahler ;
	--

This project is supported by three scientific research projects: the SOGVille project, and the two ANR projects VegDUD and EvalPDU.