Trung-Dung Vu
email: trung-dung.vu@inrialpes.fr

Julien Burlet
email: julien.burlet@gmail.com

Olivier Aycard
email: olivier.aycard@inrialpes.fr

Grid-based Localization and Local Mapping with Moving Object Detection and Tracking

Keywords: occupancy grid, simultaneous localization and mapping, moving object detection, multiple object tracking, interacting multiple model, laser radar data fusion

We present a real-time algorithm for simultaneous localization and local mapping (local SLAM) with detection and tracking of moving objects (DATMO) in dynamic outdoor environments from a moving vehicle equipped with a laser scanner, short range radars and odometry. To correct the vehicle odometry we introduce a new fast implementation of incremental scan matching method that can work reliably in dynamic outdoor environments. After obtaining a good vehicle localization, the map surrounding of the vehicle is updated incrementally and moving objects are detected without a priori knowledge of the targets. Detected moving objects are finally tracked by a Multiple Hypothesis Tracker (MHT) coupled with an adaptive Interacting Multiple Model (IMM) filter. The experimental results on datasets collected from different scenarios such as: urban streets, country roads and highways demonstrate the efficiency of the proposed algorithm.

INTRODUCTION

Perceiving or understanding the environment surrounding of a vehicle is a very important step in driving assistant systems or autonomous vehicles. The task involves both simultaneous localization and mapping (SLAM) and detection and tracking of moving objects (DATMO). While SLAM provides the vehicle with a map of static parts of the environment as well as its location in the map, DATMO allows the vehicle being aware of dynamic entities around, tracking them and predicting their future behaviors. It is believed that if we are able to accomplish both SLAM and DATMO reliably in real time, we can detect critical situations to warn the driver in advance and this will certainly improve driving safety and prevent traffic accidents. In the literature, SLAM and DATMO have been attracted considerable research works [START_REF] Prassler | Navigating a robotic wheelchair in a railway station during rush hour[END_REF] [2] [START_REF] Wang | Simultaneous localization, mapping and moving object tracking[END_REF] and they also are essential parts of the perception modules in driverless cars [START_REF] Montemerlo | Junior: The Stanford entry in the Urban Challenge[END_REF] [START_REF] Urmson | Autonomous driving in urban environments: Boss and the Urban Challenge[END_REF] winning the recent series of DARPA Grand Challenge competitions 1 . However, for highly dynamic outdoor scenarios like in crowded urban streets, there still remains many open questions. These include, how to represent the vehicle environment, how to obtain a precise location of the vehicle in presence of dynamic entities, and how to differentiate moving objects and stationary objects as well as how to track moving objects reliably over time.

In this context, we designed and developed a generic perception architecture addressing these problems focusing on outdoor dynamic environments [START_REF] Aycard | PUVAME -new French approach for vulnerable road users safety[END_REF]. The architecture (Fig. 1) is comprised of two main parts: the first part where the map of vehicle environment is constructed and dynamic objects are identified; the second part where detected moving objects are verified and tracked.

In the first part of the architecture, to model the environment surrounding the vehicle, we employ the occupancy grid framework proposed by Elfes [START_REF] Elfes | Occupancy grids: A probabilistic framework for robot percpetion and navigation[END_REF]. In order to perform mapping or modeling the environment from a moving vehicle, generally a precise vehicle localization is necessary. To correct vehicle locations from odometry, we introduce a new fast laser-based incremental localization method that can work reliably in dynamic environments. When good vehicle locations are estimated, by integrating laser measurements we are able to build a consistent grid map surrounding of the vehicle. And when new laser measurements are coming, dynamic objects can be then detected based on their discrepancies with the constructed grid map. Related results have been presented in our previous publication [START_REF] Vu | Online localization and mapping with moving object tracking[END_REF] and in this paper we employ the radar data combined with object detection results from laser data in order to obtain a more robust performance.

In the second part, previously detected moving objects in the vehicle environment are passed to the tracking process. Since some objects may be occluded or not detected, some are false alarms, tracking helps to identify occluded objects, recognize false alarms and reduce missed detections. In general, the multiple object tracking problem is complex: it involves the definition of filtering methods as well as the data association methods and maintenance of the list of objects currently present in the environment [START_REF] Bar-Shalom | Tracking and Data Association[END_REF]. Regarding the filtering techniques, Kalman filters [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF] and particle filters [START_REF] Arulampalam | A tutorial on particle filters for online non-linear/non-Gaussian Bayesian tracking[END_REF] are mostly used. These filters require the definition of a specific dynamic model of tracked objects. However, defining a suitable motion model is not trivial and Interacting Multiple Models (IMM) [START_REF] Mazor | Interacting multiple model methods in target tracking: a survey[END_REF] have been successfully applied in practice. In our previous works [START_REF] Burlet | Adaptive interactive multiple models applied on pedestrian tracking in car parks[END_REF], we have developed a fast method to adapt on-line IMM according to trajectories of detected objects and so that we obtain a suitable and robust tracker. To deal with the data association and track maintenance problem, we extend our approach to multiple objects tracking using the Multiple Hypothesis Tracking (MHT) approach [14][15]. The proposed algorithm for solving SLAM and DATMO is tested on data collected from the Daimler demonstrator car equipped with a camera, two short range radars and a laser scanner (Fig. 2). The laser scanner can detect obstacles at a range of 70m under a field of view of 160 • . It provides raw data as a list of impacts with an angular resolution of 1 • . The radars detect targets up to 30m within a field of view of 80 • and return pre-filtered data as a list of "dot" objects with their estimated positions and Doppler velocities (Fig. 2 right). In addition, vehicle odometry information such as velocity and yaw rate are provided. The measurement cycle time of the sensor system is 40ms.

Experimental platform

In our implementation, laser data is used to perform mapping as well as detection and tracking of moving objects. Radar data is then fused with laser data to confirm the results obtained by laser data in order to give a more reliable results on detection and tracking objects in the radar field of view. Images from camera are only for visualization purpose.

RELATED WORK

Before discussing in detail our approach to problems of SLAM and DATMO, it is interesting to recall some notable works in the domain.

One of the first works on SLAM with DATMO was that of Prassler's group [START_REF] Prassler | Navigating a robotic wheelchair in a railway station during rush hour[END_REF]. They described a first system on automated wheelchairs for static and dynamic object detection, moving object tracking and obstacle avoidance. The environment is represented by a time stamp grid map that provide a interesting way to detect and track moving objects. However, this method rely completely on odometry information with suppose that the odometry is ideal and it cannot detect objects moving slowly. Although the proposed solution is not really complete, it identified the need of both SLAM and DATMO for automated mobile systems.

Haehnel et al. in [START_REF] Hähnel | Map building with mobile robots in dynamic environments[END_REF] used a feature-based approach to identify pedestrians from laser range scans and use Joint Probabilistic Data Association particles filters [START_REF] Schulz | Tracking multiple moving targets with a mobile robot using particle filters and statistical data association[END_REF] to track moving pedestrians indoor. The corresponding measurements are then filtered out and classical scan registration and mapping techniques in static environment are used. However, this approach is not able to work in outdoor environment where various dynamic objects can not be described by simple features.

Wang [START_REF] Wang | Simultaneous localization, mapping and moving object tracking[END_REF] developed the first outdoor real-time system solving both SLAM and DATMO simultaneously for urban environments from a ground vehicle. To correct the vehicle odometry he used an ICP-based matching scan method and moving objects are detected based on a simple geometric analysis. He also presented a mathematical framework integrating both SLAM and DATMO and showed that they can be mutually beneficial from each other. The idea is that the results of SLAM will be more accurate if moving objects can be filtered out and thanks to a more accurate pose estimation and a better map from SLAM, DATMO can detect and track moving objects more reliably.

Recently, in the DARPA Urban Challenge competition, we have been witnessed significant advances in efforts of building autonomous vehicles. It is shown that driverless cars, for instance: Boss [START_REF] Urmson | Autonomous driving in urban environments: Boss and the Urban Challenge[END_REF] and Junior [START_REF] Montemerlo | Junior: The Stanford entry in the Urban Challenge[END_REF], are capable of operating autonomously and safely through urban-alike environments. However, testing scenarios for the competition contains only vehicles as moving objects which limits their approaches to be only able to detect and track vehicles. In addition, to obtain a good performance, participant vehicles are equipped with so many precise and expensive sensors, such as 3D laser scanners, 2D laser scanners, radars, vision, precise inertial sensors... Inspired by the pioneer work of Wang on SLAM and DATMO, our objective here is trying to put forward the state-of-the-art solutions to these tasks in order to build a general and reliable vehicle perception system with affordable sensors (e.g. 2D laser scanner, short-range radars). To this end, we introduce a new fast grid-based laser scan matching method to correct vehicle odometry that works extremely well even in the presence of dynamic entities. It will be shown later that this is an important step to build an accurate map of the environment and help to detect moving objects reliably. We also present a new approach of multiple object tracking capable of online adapting movements of moving objects which results in a more robust tracker. Parts of this work have been published separately in [START_REF] Vu | Online localization and mapping with moving object tracking[END_REF] [START_REF] Burlet | Adaptive interactive multiple models applied on pedestrian tracking in car parks[END_REF].

In the following section, we describe in detail our approach to vehicle localization and environment mapping. Algorithm for detecting moving objects is presented in Section 4. Multi objects tracking approach is detailed in Section 5. Experimental results are reported in Section 6 and finally conclusions and future works are given in Section 7.

LOCAL SLAM

To model the environment surrounding of vehicle, we employ the occupancy grid framework proposed by Elfes [START_REF] Elfes | Occupancy grids: A probabilistic framework for robot percpetion and navigation[END_REF]. Compared with feature-based approaches [START_REF] Leonard | Simultaneous map building and localization for an autonomous mobile robot[END_REF], grid maps can represent any environment and are especially suitable for noisy sensors in outdoor environments where features are hard to define and extract. Gridbased approaches also provide an interesting mechanism to integrate different kinds of sensors in the same framework taking the inherent uncertainty of each sensor reading into account.

To perform mapping, only laser data is used. For our purpose of safety vehicle navigation, a good global map is not necessary, so that the problem of revisiting or loop closing in SLAM is not considered in this work. For this reason, we propose an incremental mapping approach based on a fast laser scan matching algorithm in order to build a consistent local vehicle map. The map is updated incrementally when new data measurements arrive along with good estimates of vehicle locations obtained from the scan matching algorithm. The advantages of our incremental approach are that the computation can be carried out very quickly and the whole process is able to run online.

Notation

Before describing our approach in detail, we introduce some notations used. We denote the discrete time index by the variable t, the laser observation from vehicle at time t by the variable z t = {z 1 t , ..., z K t } including K individual measurements corresponding to K laser beams, the vector describing an odometry measurement from time t -1 to time t by the variable u t , the state vector describing the true location of the vehicle at time t by the variable x t .

Occupancy Grid Map

In this representation, the vehicle environment is divided into a two-dimensional lattice M of rectangular cells and each cell is associated with a measure taking a real value in [0, 1] indicating the probability that the cell is occupied by an obstacle. A high value of occupancy grid indicates the cell is occupied and a low value means the cell is free. Assuming that occupancy states of individual grid cells are independent, the objective of a mapping algorithm is to estimate the posterior probability of occupancy P(m | x 1:t , z 1:t ) for each cell m of the grid, given observations z 1:t = {z 1 , ..., z t } at corresponding known poses x 1:t = {x 1 , ..., x t }.

Using Bayes theorem, this probability is determined by:

P(m | x 1:t , z 1:t ) = P(z t | x 1:t , z 1:t-1 , m) . P(m | x 1:t , z 1:t-1 ) P(z t | x 1:t , z 1:t-1 ) (1)
If we assume that current measurement z t is independent from x 1:t-1 and z

Equation [START_REF] Hähnel | Map building with mobile robots in dynamic environments[END_REF] gives the probability for an occupied cell. By analogy, equation [START_REF] Wang | Simultaneous localization, mapping and moving object tracking[END_REF] gives the probability for a free cell:

P(m | x 1:t , z 1:t ) = P(m | x t , z t ) . P(z t | x t ) . P(m | x 1:t , z 1:t-1 ) P(m) . P(z t | x 1:t , z 1:t-1 ) (3)
By dividing equation ( 2) by (3), we obtain:

P(m | x 1:t , z 1:t ) P(m | x 1:t , z 1:t ) = P(m | x t , z t ) P(m | x t , z t ) . P(m) P(m) . P(m | x 1:t-1 , z 1:t-1 ) P(m | x 1:t-1 , z 1:t-1 ) (4)
If we define Odds(x) = P(x) P(x) = P(x) 1-P(x) , equation (4) turns into:

Odds(m | x 1:t , z 1:t ) = Odds(m | x t , z t ).Odds(m) -1 .Odds(m | x 1:t-1 , z 1:t-1 ) (5)
The corresponding log Odds representation of equation ( 5) is:

log Odds(m | x 1:t , z 1:t ) = log Odds(m | z t , x t ) -log Odds(m) + log Odds(m | x 1:t-1 , z 1:t-1 ) (6)
In ( 6), what we need to know are two probability densities, P(m | x t , z t ) and P(m). P(m) is the prior occupancy probability of the map cell which is set to 0.5 representing an unknown state, that makes this component disappear. The remaining probability P(m | x t , z t ), is called the inverse sensor model. It specifies the probability that a grid cell m is occupied based on a single sensor measurement z t at location x t . Fig. 3 shows the function we use to compute the occupancy probability of grid cells along a laser beam measuring a distance of d. From the log Odds representation, the desired probability of occupancy P(m | x 1:t , z 1:t ) can be easily recovered. And since the updating algorithm is recursive, it allows for the map updated incrementally when new sensor data arrives.

The second image in Fig 6 shows an example of an occupancy grid map constructed from laser measurements during the vehicle's movement. The color of grid map cell indicates the probability that corresponding space being occupied: gray=unknown, white=free, black=occupied.

Localization in Occupancy Grid Map

In order to build a consistent map of the environment, a good vehicle localization is required. Because of the inherent error, using only odometry often results in an unsatisfying map (see Fig. 4 left). When features can not be defined and extracted, direct scan matching techniques like ICP [START_REF] Besl | A method for registration of 3D shapes, Pattern Analysis and Machine Intelligence[END_REF] can help to correct the odometry error. The problem is that sparse data in outdoor environments and dynamic entities make correspondence finding difficult. One important disadvantage of the direct scan matching methods is that they do not consider the dynamics of the vehicle. Indeed we have implemented several ICP variants [START_REF] Rusinkiewicz | Efficient variants of the ICP algorithm[END_REF] and found out that scan matching results are unsatisfactory and often lead to unexpected trajectories of vehicle. This is because matching only two consecutive scans may be very hard, ambiguous or weakly constrained, especially in outdoor environment and when the vehicle moves at high speeds. An alternative approach that can overcome these limitations consists in setting up the matching problem as a maximum likelihood estimation. In this approach, given an underlying vehicle dynamics constraint, the current scan position is corrected by comparing with the local grid map constructed from all observations in the past instead of only with one previous scan. By this way, we can reduce the ambiguity and weak constraint especially in outdoor environment and when the vehicle moves at high speed. Mathematically, we calculate a sequence of poses x1 , x2 , ... and sequentially updated maps M 1 , M 2 , ... by maximizing the marginal likelihood of the t-th pose and map relative to the (t -1)-th pose and map:

xt = argmax x t {P(z t | x t , M t-1 ) . P(x t | xt-1 , u t )} (7)
In the equation ( 7), the term P(z t | x t , M t-1 ) is the measurement model which is the probability of the most recent measurement z t given the pose x t and the map M t-1 constructed so far from observations z 1:t-1 at corresponding poses x1:t-1 that were already estimated in the past. The term P(x t | xt-1 , u t ) represents the motion model which is the probability that the vehicle is at location x t given that the vehicle was previously at position xt-1 and executed an action u t . The resulting pose xt is then used to generate a new map M t according to (6): Now the question is how to solve the equation ( 7), but let us first describe the motion model and the measurement model used.

M t = M t-1 ∪ { xt , z t } (8)
For the motion model, we adopt the probabilistic velocity motion model similar to that of [START_REF] Thrun | Probabilistic Robotics[END_REF]. The vehicle motion u t is comprised of two components, the translational velocity v t and the yaw rate ω t . Fig. 5 depicts the probability of a vehicle being at location x t given its previous location x t-1 and control u t . This distribution is obtained from the kinematic equations, assuming that vehicle motion is noisy along its rotational and translational components.

For the measurement model P(z t | x t , M t-1 ), mixture beam-based model is widely used in the literature [21][22]. However, the model come at the expense of high computation since it requires ray casting operation for each beam. This can be a limitation for real time application if we want to estimate a large amount of measurements at the same time. To avoid ray casting, we propose an alternative model that only considers end-points of the beams. Because it is likely that a beam hits an obstacle at its end-point, we focus only on occupied cells in the grid map. A voting scheme is used to compute the probability of a scan measurement z t given the vehicle pose x t and the map M t-1 constructed so far. First, from the vehicle location x t , individual measurement z k t is projected into the coordinate space of the map. Call hit k t the grid cell corresponding to the projected end-point of each beam z k t . If this cell is occupied, a sum proportional to the occupancy value of the cell will be voted. Then the final voted score represents the likelihood of the measurement. Let P(M i t ) denote the posterior probability of occupancy of the grid cell M i estimated at time t (following ( 6)), we can write the measurement model under the sum following:

P(z t | x t , M t-1 ) ∝ K ∑ k=1 { P(M hit k t t-1 ) so that M hit k t t-1 is occupied } (9)
The proposed method is just an approximation to the measurement model because it does not take into account visibility constraints, but experimental evidences show that it works well in practice. Furthermore, with a complexity of O(K), the computation can be done rapidly. It remains to describe how we maximize [START_REF] Elfes | Occupancy grids: A probabilistic framework for robot percpetion and navigation[END_REF] to find the correct pose xt . Hill climbing strategy in [START_REF] Thrun | A real-time algorithm for mobile robot mapping with applications to multi-robot and 3D mapping[END_REF][22] can be used but may suffer from a local maximum. Exploiting the fact that the measurement model can be computed very quickly, we perform an extensive search over vehicle pose space. A sampling version of the motion model (Fig. 5 right) is used to generate all possible poses x t given the previous pose x t-1 and the control u t . The resulting pose will be the pose at which the measurement probability achieves a maximum value. Because of the inherent discretization of the grid, the sampling approach turns out to work very well. In practice, with a grid map resolution of 20 cm, it is enough to generate about four or five hundreds of pose samples to obtain a good estimate of the vehicle pose with the measurement likelihood that is nearly unimproved even with more samples. The total computational time needed for such a single scan matching is about 10 ms on a low-end PC. An example of scan matching result is shown in Fig. 6. The most likely vehicle pose is obtained when the laser scan is aligned with the occupied parts of the map and at the same time the vehicle dynamics constraint is satisfied.

Besides the computational effectiveness, one attraction of our algorithm is that it is not affected by dynamic entities in the environment (see Fig. 4 right). Since we only consider occupied cells, spurious regions in the occupancy grid map with low occupancy probability that might belong to dynamic objects do not contribute to the sum [START_REF] Bar-Shalom | Tracking and Data Association[END_REF]. The voting scheme ensures that measurement likelihood reach a maximum only when the laser scan is aligned with the static parts of the environment. To some meaning, measurements from dynamic entities can be considered as outliers of the alignment process. This property is very useful for moving object detection process that will be described in the next section.

Local mapping

Because we do not need to build a global map nor deal with loop closing problem, only one online map is maintained representing the local environment surrounding of the vehicle. The size of the local map is chosen so that it should not contain loops and the resolution is maintained at a reasonable level. Every time the vehicle arrives near the map boundary, a new grid map is reinitialized. The pose of the new map is computed according to the vehicle global pose and cells inside the intersection area are copied from the old map.

MOVING OBJECTS DETECTION

In the previous section, we represent how to obtain precise vehicle localization and how to build local vehicle grid map from laser data. In this section we will describe how to identify moving objects by using the previously constructed grid map. Detected objects are then confirmed using radar data and their velocities are estimated.

Using Occupancy Grid to detect Moving Objects

After a consistent local grid map of the vehicle is constructed, moving objects can be detected when new laser measurements arrive by comparing with the previously constructed grid map. The principal idea is based on the inconsistencies between observed free space and occupied space in the local map. If an object is detected on a location previously seen as free space, then it is a moving object. If an object is observed on a location previously occupied then it probably is static. If an object appears in a previously not observed location, then it can be static or dynamic and we set the unknown status for the object in this case.

Another important clue which can help to decide whether an object is dynamic or not is evidence about moving objects detected in the past. For example, if there are many moving objects passing through an area then any object that appears in that area should be recognized as a potential moving object. For this reason, in addition to the local static map M constructed as described in the previous section, a local dynamic grid map D is created to store information about previously detected moving objects. The pose, size and resolution of the dynamic map is the same as those of the static map. Each dynamic grid cell D i store a value α i indicating the number of observations that a moving object has been passed by that cell. The bigger value of α i , the more probability that object seen at D i is moving.

From these remarks, our moving object detection process is carried out in two steps as follows. The first step is to detect measurements that might belong to dynamic objects. Here for simplicity, we will temporarily omit the time index. Given a new laser scan z, the corrected vehicle location and the local static map M and the dynamic map D containing information about previously detected moving objects, state of a single measurement z k is classified into one of three types following:

state(z k ) =          static if M hit k = occupied dynamic if M hit k = f ree or D hit k > α undecided if M hit k = unknown
where M hit k and D hit k are the corresponding cells of the static and dynamic map respectively at the end-point hit k of the beam z k , α is a pre-defined threshold.

The second step when measurements that might belong to dynamic objects are determined, moving objects are then identified by clustering end-points of these beams into separate groups, each group represents a single object. Two points are considered as belonging to the same object if the distance between them is less than a threshold of 0.2 m that is chosen empirically . Fig. 7 illustrates the described steps in detecting moving objects. The leftmost image depicts the situation where the vehicle is moving along a street seeing a car moving ahead and a motorbike moving in the opposite direction. The middle image shows the local static map and the vehicle location with the current laser scan drawn in red. Measurements which fall into free region in the static map are detected as dynamic and are displayed in the rightmost image. After the clustering step, two moving objects are identified (in green boxes) and correctly corresponds to the car and the motorbike.

Fusion with radar

The general purpose of the data fusion is to provide a more reliable and more accurate model than a single data would provide. After moving objects are identified from laser data, we confirm the object detection results by fusioning with radar data and estimate velocities of the detected objects.

With the radar sensors being used, a built-in preprocessing of the radar measurements takes place, wherein reflections with a similar distance, relative velocity, and amplitude are grouped together. The radar sensors return pre-filtered data as lists of potential moving objects. The object lists of the two radars are independent from each other. Each object is provided with information about the location and the Doppler velocity. For each moving object detected from laser data as described in the previous section, a rectangular bounding box is calculated and the radar measurements which lie within the box region are then assigned to corresponding object. The velocity of the detected moving object is estimated as the average of these corresponding radar measurements.

Fig. 8 shows an example of how the fusion process takes place. Moving objects detected by laser data are displayed in red with green bounding boxes. The targets detected by two radar sensors are represented as small circles in different colors along with corresponding velocities. We can see in the radar field of view, two objects detected by laser data are also seen by two radars so that they are confirmed and their velocities are estimated. Radar measurements that do not correspond any dynamic object and fall into other region of the grid are not considered. Since the radar is setup with the field of view much smaller than the laser field of view (Fig. 2), the fusion process indeed did not help much to improve the overall detection results but we can see how detection results could be improved if more sensors available.

Fig. 8. Moving objects detected from laser data are confirmed by radar data.

MULTIPLE HYPOTHESIS TRACKING USING ADAPTIVE IMM

The aim of multi-object tracking is to estimate the number and the states of real objects evolving in the environment by generating and maintaining during time a set of tracks according to detected (observed) objects2 obtained at each step. For convenience we call track a tracked object that is composed by a list of the same object detected over time. This involves a choice of filtering methods, but also Fig. 9. Architecture of multi-object tracking system data association methods and a maintenance of the list of objects currently present in the environment. The most known techniques are the the Global Nearest Neighborhood (GNN) combined with filtering, Joint Probabilistic Data Association Filter (JPDAF) and the Multiple Hypothesis Tracking (MHT) [START_REF] Bar-Shalom | Tracking and Data Association[END_REF] [START_REF] Thrun | Probabilistic Robotics[END_REF]. In the conventional GNN only the most likely assignment is considered at each step, allowing only to associate at most one detected object to one track. The JPDAF method permits to assign several detected objects to one track by a weighted probabilistic sum. Nevertheless, it works with a fixed number of tracks and increase the track state uncertainty since several objects with different positions can update on unique track. In MHT alternative association hypotheses are built over time. In conflict situations, instead of taking an immediate decision (GNN) or combining hypotheses (JPDAF), hypotheses are propagated into the future in anticipation that it will resolve the possible association ambiguities.

The basic principle of MHT is to generate and update a set of association hypotheses during tracking process. An hypothesis corresponds to a specific probable assignment of detected objects with tracks. By maintaining and updating several hypotheses, none irreversible association decisions are made and ambiguous cases are solved in further steps. Reid introduces first a complete algorithm given a systematic way in which multiple data association hypotheses can be formed and evaluated for the problem of multiple target tracking [START_REF] Reid | An algorithm for tracking multiple targets[END_REF]. It permits to systematically generate and evaluate hypotheses by building track trees. For these reasons, we based our MHT on this efficient algorithm.

Regarding tracking techniques, Kalman filters [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF] or particle filters [START_REF] Arulampalam | A tutorial on particle filters for online non-linear/non-Gaussian Bayesian tracking[END_REF] are generally used. These filters require the definition of a specific dynamic model of tracked objects (ie, a motion model). However, defining a suitable motion model is a real difficulty. To deal with this problem, Interacting Multiple Models [START_REF] Mazor | Interacting multiple model methods in target tracking: a survey[END_REF] [START_REF] Li | A survey of maneuvering target tracking -Part V: Multiplemodel methods[END_REF] have been successfully applied in several applications.

The IMM approach overcomes the difficulty due to motion uncertainty by using more than one motion model. The principle is to assume a set of models as possible candidates of the true displacement model of the object at one time. To do so, a bank of elemental filters is ran at each time, each corresponding to a specific motion model, and the final state estimation is obtained by merging the results of all elemental filters according to the distribution probability over the set of motion models. By this way different motion models are taken into account during filtering process. In previous works [START_REF] Burlet | Pedestrian tracking in car parks: an adaptive interacting multiple model based filtering method[END_REF] [13], we have developed a fast method to adapt on-line IMM according to trajectories of detected pedestrian and so we obtain a suitable and robust tracker. In this work we extent this method in order to track dynamic objects in the vehicle environment.

As shown in Fig. 9, our multi-object tracking method is composed of four different parts:

• The first one is the gating. In this part, taking as input predictions from previous computed tracks, we compute the set of new detected objects which can be associated to each track. • In the second part, using the result of the gating, we perform object to tracks association and generate association hypotheses, each track corresponding to a previously known moving object. Output is compoed of the computed set of association hypotheses. • In the third part called track management, tracks are confirmed, deleted or created according to the association results and final track trees are output. • In the last part corresponding to the filtering step, estimates are computed for 'surviving' tracks and predictions are performed to be used the next step of the algorithm. In this part we use an adaptive method based on Interacting Multiple Models (IMM).

Gating

In this part, taking as input predictions from previous computed tracks and newly detected objects, a gating is performed. It consists in, according to an arbitrary distance function, determine the detected objects which can be associated with tracks. Also during this stage, clustering is performed in order to reduce the number of association hypotheses. It consists in making clusters of tracks which share at least one detected object. In the next stage, association can be performed independently for each cluster decomposing a large problem in smaller problems which induce generation of less hypotheses.

If we take as an example the situation depict by the Fig.

Association

In this part, taking as input clusters of tracks and detected objects validated by the gating stage, association hypotheses are evaluated. By considering likelihood of objects with tracks, new track apparition probability and non-detection probability, an association matrix is formed.

Let be L(o i ,t j ) the function giving the likelihood of object i with track j, P NT the new track apparition probability and P ND the non detection probability. Alway taking as an example the situation in the Fig. 10, the association matrix is written:

        L(o 1 ,t 1 ) -∞ P NT L(o 2 ,t 1 ) L(o 2 ,t 2 ) P NT -∞ L(o 3 ,t 2 ) P NT P ND P ND -∞        
Thus a possible association hypothesis corresponds to a valid assignation in the matrix of detected objects with tracks i.e one unique element in each row and each column is chosen to compose the assignation. In order to reduce the number of hypothesis, only the m-best association hypotheses are considered as done in Cox work [START_REF] Cox | An efficient implementation of Reid's multiple hypothesis tracking algorithm and its evaluation for the purpose of visual tracking[END_REF] using this matrix. This m-best implementation of the Reid's algorithm permits to reduce the number of hypotheses and thus to control the trees' growth in width. So for each cluster (each set of tracks sharing at least one detected objects) the m-best assignment in the association matrix are computed using the Murty method [START_REF] Murty | An algorithm for ranking all the assignments in order of increasing costs[END_REF] which computes the m-best assignations in the matrix and by this way be obtain the m-best Hypotheses.

Track management

In this third stage, using the m-Best Hypotheses resulting of the association stage, the set of track trees, is maintained i.e tracks are confirmed, deleted or created. The track management consists in only kept the branches with leafs attached to the m-best hypothesis and prune all other branches. New tracks are created if a new track creation hypothesis appears in the m-best hypotheses. A new created track is confirmed if it is updated by detected objects after a fixed number of algorithm steps (three in our implementation). Thus spurious measurement which can be detected as objects in the first step of our method are never confirmed.

If a non-detection hypothesis appears and so to deal with non-detection cases (which can appear for instance when an object is occulted by an other one, tracks without associated detected objects are updated according to their last associated objects and next filtering stage becomes a simple prediction. But if a track is not updated by a detected object for a given number of steps, it is deleted.

Furthermore, to reduce the continuously tracks' growth, an other pruning is performed. Typically trees' growth is controlled in length by the so called N-Scans pruning technique which consists in only kept the N last scans in the trees. By this way, the maximum length of tracks trees is N and it permits to apply the MHT algorithm on realistic problems. As the quality of gating relies directly on the quality of filtering and especially the prediction step, we have chosen Interacting Multiple Models (IMM) [START_REF] Mazor | Interacting multiple model methods in target tracking: a survey[END_REF][25] to deal with motion uncertainties in this filtering part.

Adaptive Filtering using Interacting Multiple Models

Besides, we developed an efficient method in which critical parameter of the IMM is on-line adapted [START_REF] Burlet | Pedestrian tracking in car parks: an adaptive interacting multiple model based filtering method[END_REF][13] according to the most probable trajectories formed by tracks. Thus as Fig. 11 shows our filtering stage is composed of three parts : an IMM filtering part, a part in which most probable trajectories are computed and a last part in which we adapt the IMM filter.

Principle:

As explained, the IMM approach overcomes the difficulty due to motion uncertainty by using a set of M elemental filters at each time, each corresponding to a specific motion model, and the final state estimation is obtained by merging the results of all filters according to the distribution probability P(µ) over the set of motion models. Also, the probability the object changes of displacement model is encoded in a transition probability matrix(TPM) which gives the distribution P( µ t | µ t-1 ), i.e the transition between models which is assumed Markovian. One cycle of an IMM is composed of tree steps (Fig. 12): A step in which filter execution is done and P(µ t ) is updated, a fusion step allowing to compute estimate fusion and a reinitialization step.

An unique filter give us the distribution at time t over object state x t knowing the current detected object z t and previous estimation P(x t |z t ). Also, P(µ t ) is updated according to P(µ t |µ t-1 ) corresponding to the TPM, it gives the transition probability between modes and so is defined as a matrix, P(µ t-1 ) is the previous distribution over models and the likelihood of the current detected object with each filter.

Thus as we use a bank of filters and we want to obtain an estimate fusion P(X t ), according to all filters outputs. The estimate fusion is obtained by:

P(X t ) = M ∑ m=1 P( [µ t = m]) P m (x t |z t ) (10)
Also during the computation process, the new distribution probability over models P(µ t ) is computed and store for each hypothesis.

To obtain new predictions, filters are reinitialized according filters outputs and in each filter the corresponding dynamic model is applied. By this way, we obtain M predictions per leaf which will be use in the gating stage. Nevertheless, to apply IMM on real applications a number of critical parameters have to be defined, for instance the set of motion models and the transition probability matrix(TPM). To cope with this design step which cannot match the reality, we propose an efficient method in which the TPM is adapted online.

Definition of our IMM:

The first step to apply our method is to define an appropriate IMM and, in particular, models which compose it. In specific applications, different objects such as cars or motorcycles can move in any directions and can often change theirs motions. Thus in our aim we choose various IMM's models to cover the set of possible directions and velocities. In previous work for one pedestrian tracking we focused only on directions but here we focus on a range of velocities while keeping a set of directions to cope with directions' transition in vehicles' behavior. As each filter corresponds to a specific motion model, we have to define each motion model. So, assuming we have different possible velocities defined according to the vehicle velocity and eight directions in the set of possible directions an object can follow, we obtain sixteen motion models (Fig. 13).

Hence, according to the definition of these sixteen motion models, our IMM is composed of sixteen filters. Kalman filters are chosen for implementation as they allow fast computation.

We must usually also define the TPM. As we develop a method which computes the TPM online, we do not need specific informations concerning the TPM and no modeling are needed. So the TPM is initially chosen to be uniform. As eight modes are defined, the TPM is an uniform square 16 × 16 matrix. In the next part of the text, we will see how the TPM is on-line adapted.

Computation of the most probable trajectories

Once estimates are performed in all track trees leafs, the most probable trajectory is computed for each track. Basically, it consists in taking the branch having the maximum probability (computed during filtering) to obtain one unique hypothesis for one given track tree. This step permits to give users more readability on what is happening during tracking process and also permits us to adapt on-line the IMM parameter according to these trajectories.

Adaptation of the IMM

To adapt the TPM in our specific situation i.e tracking detected objects, most probable trajectories are considered. Taking as input the set of trajectories computed during filtering process, we will adapt one-line the TPM of the IMM filter in order to obtain a better transition between motion models close to the real behavior of tracked objects.

The principle is the following. For a given number N of trajectories we build sequences of associated models probabilities.And then, using this models probabilities, the TPM is adapted and reused in the IMM filters for the next estimations.In more details, algorithm 1, given in pseudo-code, is the algorithm defined to compute one adaptation of the TPM.

An adaptation of the TPM is done after a given number N of trajectories obtained from tracks, to update TPM using a window on trajectories (cf. loop line 3-19 of algorithm 1). Moreover trajectories are processed one by one in three steps:

1-Models' probabilities are collected by travel through the computed most probable sequence 2-Most probable models' sequence is computed 3-Most probable models' transitions are quantified Collection of models' probabilities: For each part of a given most probable trajectory computed in last stages of the filtering process, we collect the distribution over models(lines 7). Thus a model probabilities' sequence S n obtained in such a way and is stored to be processed (line 8).

Computation of the most probable model sequence:

In a next step, the most probable models' sequence of S n is computed (line 11). More precisely, considering the actual TPM and a set S n = µ 0 ...µ K of model probabilities through time 0 to K, we aim to obtain the most probable models' sequence knowing the estimates computed by the IMM: /* Store µ k ,...µ k ′ from T n the most probable n th trajectory */ 6:

Max P(µ 0 µ 1 ...µ k | x 0 x 1 ... x K ) (11
By this way an on-line adaptation of the TPM is obtained. Thus, the effectiveness of filtering part of our MHT is improved since the prediction quality is enhanced by our method. And so, the quality of the whole MHT is improved. Following the numeration of the different motion models defined in Fig. 13, the 16 × 16 frequencies matrix are shown on Fig. 14, Fig. 15 and Fig. 16 at three different steps of the execution process. We can see that after five trajectories some transitions appear to be more frequent than other (Fig. 14). Also, after twenty five trajectories (Fig. 15) the continuous adaptation makes appear clearly different behaviors, especially transitions between models oriented to the front and the back of the vehicle (models number from two to eight and from nine to fifteen) 3 . After a number of trajectories, an efficient model of the real objects' behaviors is obtained.

Example of adaptation result:

Without our automatic and one-line adaptation it would be difficult to model such behaviors a priori and impossible to continuously model the real behavior of objects during one or several processes. Furthermore, obtain a TPM which model the real objects' motion improve the quality of the IMM filtering and thus the quality of the whole filtering process.

EXPERIMENTAL RESULTS

Our proposed algorithms for objects detection and tracking is tested on datasets collected with the DaimlerChrysler demonstrator car. The vehicle was driven through different kinds of scenarios such as city streets, country roads and highways with a maximum speed of 120 kph. In our implementation, the width and height of local grid map are set to 160 m and 200 m respectively, and the grid resolution is set to 20 cm. Every time the vehicle arrives at 40 m from the grid border, a new grid map is created. The object detection is run for every new laser scan and the tracking process is updated according to the detection results.

Fig. 17 shows some snapshots of the moving object detection and tracking process in action. The images in the first row represent online maps and objects moving in the vicinity of the vehicle are detected and tracked. The current vehicle location is represented by blue box along with its purple trajectories after corrected from the odometry. The red points are current laser measurements that are identified as belonging to dynamic objects. Green boxes indicate detected and tracked moving objects with corresponding tracks displayed in different colors. Information on estimated velocities is displayed next to detected objects. The second row are images for visual references to corresponding situations.

In the figure, the leftmost column depicts a scenario where the demonstrator car is moving at a very high speed of about 100 kph while a car moving in the same direction in front of it is detected and tracked. On the rightmost is a situation where the demonstrator car is moving at 50 kph on a country road. A car moving ahead and two other cars in the opposite direction are all recognized. Note that the two cars on the left lane are only observed during a very short period of time but both are detected and tracked successfully. The third situation in the middle, the demonstrator is moving quite slowly at about 20 kph in a crowded city street. Our system is able to detect and track both the other vehicles and the motorbike surrounding. In all three cases, precise trajectories of the demonstrator are achieved and local maps around the vehicle are constructed consistently. In our implementation, the computational time required to perform both SLAM and DATMO for each scan is about 20 -30 ms on a 1.86GHz, 1Gb RAM laptop running Linux. This confirms that our algorithm is absolutely able to run synchronously with data cycle in real time. More The table above shows quantitative results obtained using our method on three sequences of different types of environments. The first column are different types of scenario. The second column shows the numbers of real objects which entered the vehicle's sensors range which is manually counted. The third number corresponds to the numbers of steps in our algorithm in which one object is not detected but always tracked (non-detection cases). The fourth column are the numbers of false alarms i.e when our detector (in some cases because of vehicle sensors noise) detected moving objects but our tracker recognized these detection. In the last column are the total numbers of tracks computed during the given sequence.

Quantitative results

Data

The results show that during three sequences, most part of object appearances are tracked. We can note that the number of tracks remains more important than the number of real objects. It is due to objects which moves across or close to the sensors' range boundary. Indeed, close to the sensors' range boundary, laser sensor loose precision and so the detection stage become less efficient. Then if an object reappears in the sensor range it is so considered as a new one by our tracker. Also, even if an important number of non-detections and false alarms appears, the tracking part permits to cope with such problems especially since the quality of prediction step is greatly improved by our adaptive IMM. Our two stage approach permits to cope with sensors noise where an efficient detection is reinforced by a robust tracking of objects.

CONCLUSIONS AND FUTURE WORKS

We have presented an approach to accomplish local mapping with detection and tracking moving objects. Experimental results have shown that our system can successfully perform a real time mapping and moving object tracking from a vehicle at high speeds in different dynamic outdoor scenarios. This is done based on a fast scan matching algorithm that allows estimating precise vehicle locations and building a consistent map surrounding of the vehicle. After a consistent local vehicle map is build, moving objects are detected and are tracked reliably using an adaptive Interacting Multiple Models filter coupled with an Multiple Hypothesis tracker.

Future works include incorporating object models of several predefined classes with specific shapes and sizes that give a more meaningful representation of detected objects instead of only sets of contour points as in our current work. In addition, algorithms of road detection and road type classification based on constructed grid map are being considered. The motivation is that road detection can help object detection step to filtering out noisy and irrelevant data off-the-road and focus more on road-likely regions. In all, the fusion of a vehicle map, road detection, moving object classification and tracking modules certainly will enable a better interpretation of driving situations.

Fig. 1 .

 1 Fig. 1. Architecture of the perception system

Fig. 2 .

 2 Fig. 2. Left: the Daimler demonstrator car. Right: an example of sensor data, laser measurements are displayed in small red dots and radar measurements displayed as bigger dots.

Fig. 3 .

 3 Fig. 3. Profile of an inverse sensor model illustrates the occupancy probability along a laser beam measuring a distance of d.

Fig. 4 .

 4 Fig.4. Hit maps build directly from raw laser data collected from a vehicle moving along a straight street: with vehicle localization using odometry (left); and using results of scan matching (right). Note that the scan matching results are not affected by moving objects in the street.

Fig. 5 .

 5 Fig. 5. The probabilistic velocity motion model P(x t | x t-1 , u t ) of the vehicle (left) and its sampling version (right).

Fig. 6 .

 6 Fig.6. An example of scan matching. From left to right: reference image; map constructed so far M t-1 with previous vehicle location x t-1 ; new laser measurement z t ; and matching result is obtained by trading off the consistency of the measurement with the map and the previous vehicle pose.

Fig. 7 .

 7 Fig. 7. Moving object detection example. See text for more details.

Fig. 11 .

 11 Fig. 11. Principle of our adaptive filtering program

Fig. 12 .

 12 Fig. 12. Principle of IMM

Fig. 13 .

 13 Fig. 13. The sixteen chosen motion models in the vehicle's frame

) Algorithm 1 :

 1 Adaptive IMM Algorithm 1Adaptation of TPM(T 0 , ..., T N )

Fig. 14 .Fig. 15 .Fig. 16 .

 141516 Fig. 14. Frequencies matrix obtained after five trajectories

Fig. 17 .

 17 Fig. 17. Experimental results show that our algorithm can successfully perform both SLAM and DATMO in real time for different environments results and videos can be found at http://emotion.inrialpes.fr/ tdvu/videos/.

  10, in this stage one set is computed as T 1 and T 2 share object O 2 . Also according to gates, objects O 1 and O 2 can be assigned to T 1 and objects O 2 and O 3 to T 3 .

	T1	T1's gate	T2's gate	T2
			O3
		O1
		O2
		Fig. 10. Example of association problem

  Type Real Objects Non-detections False Alarms Total Tracks

	City	57	7%	3%	88
	Road	74	11%	3%	109
	Highway	5	7%	1.5%	47

www.darpa.mil/grandchallenge/index.asp

usually the term observation is used in such a case but as in our work raw sensor observations are treated to obtain detections, the term detected object will be used for more clarity

According to nonholonomic constraints we cannot obtain direct transitions from the front model to the back model for instance but as shown in the results transitions between adjacent models occur

http://www.prevent-ip.org/profusion

ACKNOWLEDGMENTS

The work is supported by the European project PReVENT-ProFusion 4 and partially by the French Délégation Générale pour L'Armement (DGA).

for all Ob ject pose x k in T n do 7: {µ k } ← T n (k) We just need to obtain the maximum of the distribution P(µ 1 µ 2 ...µ K | x 0 x 1 ... x K ), thus the inference is made using the Viterbi Data Algorithm [START_REF] Forney | The Viterbi algorithm[END_REF]. As complexity of this algorithm is in O(KM 2 ), we efficiently obtain the most probable models' sequence.

Quantification of most probable model transitions: Using this most probable models' sequence, the number of transitions from one model to an other is quantified (lines 13 to 17). To do so a frequencies matrix is considered. This matrix models the number of transitions which have occurred from one model to an other. We note F this matrix and so F i j gives the number of transitions which has occurred from model i to j. Using the most probable models' sequence corresponding to a specific trajectory and computed by the Viterbi algorithm, the update of F is directly obtained by counting transitions in this sequence. Furthermore, F is kept in memory to be used in next adaptation and before the first update all its elements are set to 1.

Finally, when N trajectories have been treated, the new TPM is obtained by normalization of the frequencies matrix F. Thus the TPM is re-estimated using all model sequences S 1 ...S N and is reused in the IMM for next executions (lines 21 and 22). In practice, before the first run, the TPM is chosen uniform (according to F initialization) as we do not want to introduce a priori data.


Damien J Allain

Symmetric Infinite Impulse Response filtering: methods and application on ocean inertial waves

Keywords:

Symmetric Infinite Impulse Response (IIR) filters have been designed to filter geophysical time series, in particular to remove the signature of Inertial Waves (IWs) as it is often so strong it disturbs the observations and analyses of other phenomena. As the symmetric filters are based on IIR digital filters, their computational costs are reduced to a minimum and their parameters can be adapted to all latitudes when used on IWs. The Symmetric Band-Cut (SBC) filter has been tested. The results show that the SBC filter can remove a lot of the signal from IWs but that a good understanding of the filter is necessary to interpret its output, as it shows transition and boundary effects. For that reason, a synthetic case has been presented to enable users to get the best understanding of the filter.

Introduction

Inertial Waves (IWs) are described by e.g. Millot (1990, sec.6).

IWs have a frequency spectrum that is spread around the Coriolis Frequency (CF). In this article, the CF at will be noted with the latitude as subscript: at a latitude of φ the CF is ω φ .

Taking Ω the spinning frequency of the Earth gives:

ω φ = 2Ω sin φ (1a)
2014/07/08 17:46 1 Symmetric IIR filtering of inertial waves

In an even simpler way, taking ω K2 the average frequency of the K2 tidal wave gives:

ω φ = ω K2 sin φ (1b)
This means:

0 ≤ ω φ ≤ ω K2 (2)
IWs can have a strong signature on analysed parameters, e.g. with water velocity at or near the surface of the deap ocean. This signature is often so strong it disturbs the observations and analyses of other phenomena. For an example, see [START_REF] Petrenko | Barotropic eastward currents in the western Gulf of Lion, north-western Mediterranean Sea, during stratified conditions[END_REF].

Filtering out the signature of IWs is very different from removing the tidal signature. Tides can also have a strong signature on analysed parameters. This signature is also often so strong it disturbs the observations and analyses of other phenomena. Thankfully tides have an almost discrete frequency spectrum. So a frequency analysis will give the amplitudes and phases of the frequency constituents. From these a prediction can be computed. Subtracting this prediction removes the tidal signature. This is called detiding.

Unlike tides, the frequency spectrum of IWs is continuous: it has an infinite number of constituents. This means the only way to remove or attenuate the signature of IWs is filtering out their frequency domain from the signal.

In sec. 2 a base type of digital filter will be chosen and its basic characteristics will be studied. In sec. 3, its main limitation for its use for, at least, geophysical applications will be addressed giving a final type of filter. In sec. 4 the results of a test of the filter on real data will be presented.

Filter choices

Type

It is currently common to take a convolution kernel to carry out filtering in oceanography. An example is the [START_REF] Demerliac | Calcul du niveau moyen journalier[END_REF] filter, whose convolution kernel spans 72 hours. With ever-increasing resolutions, we are currently reaching a point at which simulations on a grid of 40•10 6 vertexes (for example on a 10 3 × 10 3 × 40 grid) or more are common place. A 72 steps filter on such a grid would require 72 × 40•10 6 3•10 9 multiplications per time step and, coding numbers with 8 B double precision type, would take more than 23 GB of memory, which is currently feasible, but still quite a challenge. It will also be shown, e.g. in sec. 2.3.3, that an appropriate filter for IWs can take several hundreds of elements. So, because of its need in memory and its cost in computation, a convolution kernel is not practical for filtering the signature of IWs.

But digital filtering is not limited to convolution. In particular, an Infinite Impulse Response (IIR) filter can be chosen because it can have a very large convolution kernel with a very limited Number of Delay Registers (NDR). For simplicity, Direct-Form (DF) IIR filters are chosen. There are two type of DF IIR filters. DF-I IIR filters do not have an internal register that is vulnerable to overflow when using fixed-point calculation [START_REF] Smith | Introduction to Digital Filters with Audio Applications[END_REF], but have a NDR that is twice the order. As floating-point calculation is used here, and in order to limit further the NDR, a DF-II IIR filter, whose NDR is the same as the order of the filter, is chosen.

As geophysicists are not all familiar with IIR filters, mainly used in electronics, a lot of relevant details will be given. In sec. 2.2 the set-up of the DF-II IIR filter will be detailed. Because the Band-Pass (BP) filter is so complementary to the Band-Cut (BC) filter that understanding the characteristics of the BP filter helps understanding the characteristics of the BC filter, both BP and BC filters will be presented. In sec. 2.3 and 2.4 the parameters of, respectively, the BP and BC filters will be calculated. In sec. 2.5 the phase shift and its importance will be exposed.

Order

Taking z = e jωs -1 the transfer function of a 1 step advance with s the sampling frequency, we have the diagram of second-order DF-II IIR filter on fig. 1. The input I is fed to an adder. Smith (2007, fig.9.2) The output B of this adder is connected to the input of the first delay register. The filter has 2 cascaded delay registers. The outputs of the delay registers are not only multiplied by b 1 and b 2 and added to B, itself multiplied by b 0 , to give the output O but also multiplied by a 1 and a 2 and added to I to give B. This gives, for every element n of the input I, the intermediate value B and the output O:

z -1 a 1 B O z -1 a 2 I b 1 b 2 b 0 Figure 1: Second-order DF-II IIR fil- ter, from
B n = I n + a 1 B n-1 + a 2 B n-2 (3a) O n = b 0 B n + b 1 B n-1 + b 2 B n-2 (3b) = b 0 (I n + a 1 B n-1 + a 2 B n-2 ) + b 1 B n-1 + b 2 B n-2 (3c)
This feed-back makes the response of the filter infinite. Eq. (3) gives the frequency transforms of I, B and O:

F (B) = F (I) + F (B) a 1 z -1 + a 2 z -2 (4a) ⇔ F (B) = F (I) 1 1 -a 1 z -1 -a 2 z -2 (4b) F (O) = F (B) b 0 + b 1 z -1 + b 2 z -2 (4c)
which gives the transfer function H(ω):

H(ω) = b 0 + b 1 z -1 + b 2 z -2 1 -a 1 z -1 -a 2 z -2 (4d)
We are now left with the choice of constants to tune the filter to our needs.

Band-Pass

Constants

For a BP transfer function H p , we need, with z 0 = e jω 0 s -1 and ω s = ω 0 s -1 :

   H p (0) = 0 H p (πs) = 0 H p (ω 0 ) = 1 ⇔    b 0 +b 1 +b 2 = 0 b 0 -b 1 +b 2 = 0 b 0 +b 1 z -1 0 +b 2 z -2 0 = 1 -a 1 z -1 0 -a 2 z -2 0 (5a)
2014/07/08 17:46 3 Symmetric IIR filtering of inertial waves

⇔        b 1 +b 2 = -b 0 -b 1 +b 2 = -b 0 b 1 cos ω s +b 2 cos 2ω s +a 1 cos ω s +a 2 cos 2ω s = 1 -b 0 -b 1 sin ω s -b 2 sin 2ω s -a 1 sin ω s -a 2 sin 2ω s = 0 (5b) ⇔        b 1 = 0 b 2 = -b 0 -b 0 cos 2ω s +a 1 cos ω s +a 2 cos 2ω s = 1 -b 0 b 0 sin 2ω s -a 1 sin ω s -a 2 sin 2ω s = 0 (5c)
which, when sin ω s cos 2ω s -cos ω s sin 2ω s = sin(ω s -2ω s ) = -sin ω s = 0, gives:

           b 0 = b p b 1 = 0 b 2 = -b p a 1 = 2 (1 -b p ) cos ω s a2 = 2b p -1 (5d)
or:

H p (ω) = b p -b p z -2 1 -a 1 z -1 -a 2 z -2 = b p (1 -z -2 ) 1 -a 1 z -1 -a 2 z -2 = b p (z -z -1 ) z -(1 -b p ) z 0 + z -1 0 -(2b p -1) z -1 (5e)

Responses

The convolution kernel is obtained by computing the impulse response. As both the convolution kernel and the impulse response are equivalent, both terms will be used when describing this feature of the filters presented in this article.

The amplitude and phase diagrams and the convolution kernel are shown on fig. 2 for different values of b p . The amplitude and phase diagrams show the specifications of the BP filter, given in eq. (9a), are respected. The amplitude diagram shows that, the lower b p is, the sharper the filter is. The phase diagram shows that when ω → 0 the phase tends to π/2 and that when ω → 2πs the phase tends to -π/2. It also shows that the lower b p is, when ω decreases or increases from ω 0 , the faster the phase will go towards π/2 or -π/2. The impulse response diagram shows slowly decaying reponses. It also shows that the lower b p is, the slower the decay is. For higher values of b p , the impulse reponses show oscillations.

Quality factor

Fig. 2 shows the convolution kernel can be very wide. Estimating the width of the convolution kernel is important in order to take properly into account the contribution of the last elements of the input. For this, both definitions of the quality factor Q are used. We have indeed, with ∆ω the bandwidth:

Q = ∆ω ω 0 (6a)
and also, with E the energy contribution of an element of the input and t the amount of time since that element occurred, derived from Smith (2007, Quality Factor):

Q = ω 0 E ∂t ∂E (
E ∝ exp -t ω 0 Q (6c)
If we want 99.9% of the contribution of the last element, we must take an extra l steps:

l = -ln 10 -3 Q ω 0 s 7 Q ω 0 s = 7 ∆ω ω 2 0 s (6d)
We can approximate ∆ω by deriving the transfer function around the peak, taking = js -1 ∂ω:

z + ∂z = z + js -1 z∂ω = z(1 + ) z -1 + ∂z -1 = z -1 -js -1 z -1 ∂ω = z -1 (1 -) H(ω 0 + ∂ω) = b p z 0 (1 + ) -z -1 0 (1 -) z 0 (1 + ) -(1 -b p ) z 0 + z -1 0 -(2b p -1) z -1 0 (1 -) = b p z 0 -z -1 0 + (z 0 + z -1 0 ) z 0 + z 0 -z 0 -z -1 0 + b p z 0 + z -1 0 -2b p z -1 0 + 2b p z -1 0 + z -1 0 -z -1 0 = b p z 0 -z -1 0 + (z 0 + z -1 0 ) z 0 + b p z 0 -z -1 0 + 2b p z -1 0 -z -1 0 = b p z 0 -z -1 0 + (z 0 + z -1 0 ) z 0 -z -1 0 (b p + ) + 2b p z -1 0 = z 0 -z -1 0 + z 0 + z -1 0 z 0 -z -1 0 + b -1 p z 0 -z -1 0 + 2 z -1 0 = 2j sin ω s + 2 cos ω s 2j sin ω s + 2j b -1 p sin ω s + 2 z -1 0 = sin ω s -j cos ω s sin ω s + b -1 p sin ω s -j z -1 0 (7a)
2014/07/08 17:46 5 Symmetric IIR filtering of inertial waves

H(ω 0 + ∂ω) = sin ω s + s -1 cos ω s ∂ω sin ω s + s -1 jb -1 p sin ω s + z -1 0 ∂ω = sin ω s + s -1 cos ω s ∂ω sin ω s + s -1 jb -1 p sin ω s + cos ω s -j sin ω s ∂ω (7b) |H(ω 0 + ∂ω)| = [sin ω s + s -1 cos ω s ∂ω] 2 [sin ω s + s -1 cos ω s ∂ω] 2 + s -1 sin ω s b -1 p -1 ∂ω 2 (7c)
as ∂ω → 0:

|H(ω 0 + ∂ω)| → 1 - s -1 sin ω s b -1 p -1 ∂ω 2 [sin ω s ] 2 → 1 -2 -1 s -1 b -1 p -1 ∂ω 2 (7d)
when b p → 0:

|H(ω 0 + ∂ω)| → 1 - ∂ω 4b p s 2 (7e)
from which we can approximate the quality factor:

Q w 0 2b p s (8a)
which is given in most electronics text books but I did not find its derivation, as in ( 7), in the many ones I searched in. It gives the following simple relation for 99.9% contribution:

l 7b -1 p (8b)
It will be shown in sec. 4 that, for applications on IWs, we can have b p 0.07, making l 100, which is particularly large for a convolution kernel.

Band-Cut

Constants

For a BC transfer function, we need, with z 0 = e jω 0 s -1 and ω s = ω 0 s -1 :

   H c (ω 0 ) = 0 H c (0) = 1 H c (πs) = 1 ⇔        b 1 cos ω s +b 2 cos 2ω s = -b 0 -b 1 sin ω s -b 2 sin 2ω s = 0 b 1 +b 2 +a 1 +a 2 = 1 -b 0 -b 1 +b 2 -a 1 +a 2 = 1 -b 0 (9a)
which, when sin ω s cos 2ω s -cos ω s sin 2ω s = sin(ω s -2ω s ) = -sin ω s = 0, gives:

           b 0 = b c b 1 = -2b c cos ω s b 2 = b c a 1 = -b 1 = 2b c cos ω s a 2 = 1 -2b c (9b)
or:

H c (ω) = b c + b 1 z -1 + b c z -2 1 + b 1 z -1 -a 2 z -2 = b c z + b 1 + b c z -1 z + b 1 -a 2 z -1 = b c z + z -1 -z 0 -z -1 0 z -b c z 0 + z -1 0 -(1 -2b c ) z -1 (
H c (ω) = (1 -b p ) z + z -1 -z 0 -z -1 0 z -(1 -b p ) z 0 + z -1 0 -(2b p -1) z -1 = 1 - z -(1 -b p ) z 0 + z -1 0 -(2b p -1) z -1 -(1 -b p ) z + z -1 -z 0 -z -1 0 z -(1 -b p ) z 0 + z -1 0 -(2b p -1) z -1 = 1 - b p (z -z -1 ) z -(1 -b p ) z 0 + z -1 0 -(2b p -1) z -1 = 1 -H p (ω) (9d)
which shows that the BC filter is the complement of the BP filter. 9a), are respected. The amplitude diagram shows that, the lower b p is, the sharper the filter is. It must also be noted that the derivative around the minimum is not 0. This means the filter will have an important leak around its cut frequency. The phase diagram shows that when ω → ω 0 the phase tends to -π/2 when ω < ω 0 and to π/2 when ω > ω 0 . It also shows that the lower b p is, when ω decreases or increases from ω 0 , the faster the phase will go towards 0. The impulse response diagram shows the response is the exact difference between 2014/07/08 17:46 the impulse itself and the impulse response of the BP filter, again showing that the BC filter is the complement of the BP filter. So, for the same value of b p , the width of the convolution kernel of the BC filter is the same as that of the BP filter.

Responses

Phase shifts

It must be highlighted that the filters described in the previous sections have a phase shift. Indeed their impulse responses are not symmetric. This is not relevant in the case of a spectral analysis, where the output will only be the amount of energy in a frequency band. Such analyses are of interest, e.g. in the search of a cause of an error known to be around a certain frequency.

Symmetric filtering

However filtering without phase shift is necessary [START_REF] Demerliac | Calcul du niveau moyen journalier[END_REF]. This means the transfer function only has a real part. This can be done by multiplying it with its complex conjugate:

H(ω) = b 0 + b 1 z + b 2 z 2 1 -a 1 z -a 2 z 2 (10)
Eq. ( 10) shows this can be done by running the filter again with, as its input, the reversed output of the previous filter run. This can also be found intuitively: a reverse filter run will shift whatever has been shifted one way back the other way. It is therefore only possible during post-processing. Users who wish to remove the signature of a signal with a certain frequency spectrum on sub-sampled outputs must be reminded to filter the aliased spectrum, and so to also filter whatever else was at that aliased spectrum before aliasing.

The frequency and impulse responses diagrams of the SBP and the SBC filters are shown on fig. 4 for different values of b p . The frequency response of a symmetric filter is just the square of the amplitude response of its base filter. So, again, the frequency response diagrams show that, the lower b p is, the sharper the filters are. The impulse response diagram show the response of the SBC is the exact difference between the impulse itself and the impulse response of the SBP filter. The frequency response diagrams also show the frequency response of the SBC is the exact difference between 1 and the frequency response of the SBP filter. So both types of diagrams show that the SBC filter is the complement of the SBP filter. It must also be noted that the derivative around the minimum is 0 for both filters. This means that, unlike the BP filter, the SBP filter will not have an important leak around its cut frequency making it more efficient at removing signals with a continous spectrum like IWs. The impulse responses diagrams show that this time the impulses reponses are symmetric around 0. Again the impulse response diagrams show slowly decaying reponses, but this time the decay are on both sides of the center. They also again show that the lower b p is, the slower the decays are. The impulse reponses show again oscillations for higher values of b p , actually for the very same values of b p for which the base filters show oscillations themselves.

Test results

Model output

In order to illustrate the behaviour of these filters on real data, an example was taken from a model run in the Bay of Biscay. The model used was SYMPHONIE version S-2010.17 Klaus storm on 23 to 24 Jan 2009. Its grid has 362×370 tracer points and 43 levels, so more than 3.9•10 6 values per parameter per time frame. The saved output has an hourly sampling and spans from 1 Dec 2008 00:00 to 28 Feb 2009 00:00, so also covers the days of the Quinten storm on 9 Feb 2009. The first 15 days were remove from this analysis as spin-up time.

The most relevant data point was searched by running a spectral analysis of the sea water velocity components by taking the Root Mean Square (RMS) of the output of a BP filter taking, with a latitude between 43 • and 51 • , between 17.5h or 20.5 • h -1 and 15.4h or 23.4 • h -1 therefore centred around 16.5h or 21.8 • h -1 with a quality factor of 16. The treatment took around 6min for each components.

A set of frames, covering the period of IWs, of the filtered output are shown on fig. 5. Comparing the directions between the first frame (1440) and the other frames gives the period of the oscillations. On the continental shelf, i.e. north of a line going through 44N 2W, 46N 4W and 48N 8W, the arrows show the directions have a period around 12.5h (frames 1452 and 1453). This shows the dominant signal in this area, even after filtering, remains the diurnal tides. In the deep ocean the directions have the inertial period around 16.5h (frames 1452 and 1453). All the frames also show that the norm of the velocity remains the same from frame to frame. This shows the particles of water are just moving in circles. The maximum speed being around 0.15m s -1 , the circles are up to around 1.4km in diameter.

The results of the analysis are shown on fig. 6. Data on the continental shelf is discarded from this study as signal from the strong tide leaked through the filter. So only the deep 2014/07/08 17:46

.1m s -1 / √ • h -1 .
The frequency analysis of the model output shows a strong large peak between 20 and 22 • h -1 with a maximum of 0.043m s -1 / √ • h -1 at 21.6 • h -1 which is very close to the the inertial frequency of 21.41 • h -1 . This peak is even larger than the peaks of the diurnal tides. There are 3 sharp peaks at the frequencies of semi-diurnal tides between 26.4 and 32 • h -1 . These peaks all reach between 0.014 and 0.016m

s -1 / √ • h -1 for N2 (28.4 • h -1 ), M2 (29 • h -1
) and S2 (30 • h -1 ) tides. Otherwise, this analysis shows what can be interpreted as noise with most amplitude values below 0.006m s -1 / √ • h -1 between 10 • h -1 and 32 • h -1 and all amplitude values below 0.002m s -1 / √

• h -1 above 32 • h -1 .
In particular, it shows no particular peak at the frequencies of diurnal tides between 11.7 and 16.7 • h -1 . The frequency analysis of the detided model output has only few differences with that of the model output. Only the peaks of the diurnal tides have been removed, although their remains a sharp peak just below the N2 frequency (29 • h -1 ) that reaches less than 0.011m s -1 / √ • h -1 . In particular, the IWs peak is preserved by the detiding as the frequencies given to the detidor were far from its frequency domain. The frequency analysis of the filtered detided model output has only few differences with that of the detided model output. Only the IWs peak has been removed, reducing the values in its frequency domain to levels very similar to the noise around.

What model output shows low frequency variations and some strong high frequency oscillations especially during the second half of the presented series. This particular part, from day 38, i.e. 22 Jan, starts with 6 days of oscillations with an amplitude around 0.15m s -1 followed by 8 days around 0.1m s -1 then, from 5 feb, by 5 days of oscillations with an amplitude around 0.15m s -1 then by 6 days of oscillations with an amplitude around 0.2m s -1 then finally decaying to less than 0.1m s -1 in the last 12 days. The detided model output shows only small differences with the model output, in particular preserving the high frequency oscillations. The filtered detided model output removes most of the strong oscillations while preserving the low frequency variations. It does particularly well during the last quarter of the presented series, see the peak of day 56, i.e. 9 Feb, the day of the Quinten storm, but not so well during the first 10 days, when it only halfes the oscillations. This is because the filter is symmetric and therefore half of its convolution kernel sees almost no oscillation.

These results show that the filter can remove a lot of the signal from IWs but that a good understanding of the filter is necessary to interpret its output.

Synthetic case

In order to improve the understanding of the reader, a synthetic case is also presented fig. 8. The synthetic case also lasts 75 days. It starts at a value of 0.2m s -1 with a smooth transition to 0m s -1 over less than 25 days. From day 30 to the end, the case oscillates at 25 • h -1 . The case oscillates between 0m s -1 and 0.4m s -1 decaying to 0.25m s -1 at day 44. On day 45, the sign 2014/07/08 17:46 13 Symmetric IIR filtering of inertial waves is inverted. From day 45 to the end, the case oscillates between 0m s -1 and -0.25m s -1 decaying to -0.2m s -1 at day 75. It must be noted that the case shows oscillations of almost 0.01m s -1 of amplitude on day 1, decaying to almost 0 on day 10. The filtered case shows oscillations from around day 22, increasing from 0m s -1 to 0.1m s -1 of amplitude on day 30, then decreasing to 0.004m s -1 on day 37 then increasing to 0.22m s -1 at the sign invertion on day 45 then decreasing to 0m s -1 on day 56. The filtered case also shows oscillations from around day 65, increasing from 0m s -1 to 0.05m s -1 on day 75. This illustrates transition and boundary effects of the filter. There are small oscillations on step responses. The filter only halves input oscillations when they start and finish, because the filter is symmetric and therefore half of its convolution kernel sees no oscillation. The filter does not attenuate any oscillation around a sign invertion, because the filter is symmetric and therefore half of its convolution kernel sees oscillations one way and the other half sees oscillations the other way, but sign invertions will be most unlikely on real cases.

Conclusion

Symmetric filters have been designed to filter geophysical time series. As the symmetric filters are based on IIR digital filters, their computational costs are reduced to a minimum and their parameters can be adapted to all latitudes when used on IWs. The SBC filter has been tested. The results show that the SBC filter can remove a lot of the signal from IWs but that a good understanding of the filter is necessary to interpret its output, as it shows transition and boundary effects. For that reason, a synthetic case has been presented to enable users to get the best understanding of the filter.

Figure 2 :

 2 Figure 2: Amplitude and phase diagrams and convolution kernels of the DF-II IIR BP filter for different values of b p

  9c) 2014/07/08 17:46 6 Symmetric IIR filtering of inertial waves Taking b c = 1 -b p gives:

Figure 3 :

 3 Figure 3: Amplitude and phase diagrams and convolution kernels of the DF-II IIR BC filter for different values of b p

Figure 4 :

 4 Figure 4: Frequency (left) and impulse (right) responses diagrams of the DF-II IIR SBP (top) and SBC (bottom) filters for different values of b p

Figure 5 :Figure 6 :

 56 Figure 5: Response of the sea water velocity filtered around the inertial frequency at the surface on 13 Feb 2009 from 00:00 (frame 1440) to 17:00 (frame 1457). The background color shows the norm of the velocity and the arrows show the directions. Values range from 0m s -1 (light blue) to 0.16m s -1 (dark purple) and isolines are every 0.02m s -1 .

Figure 7 :

 7 Figure 7: Eastward sea water velocity at 45.3811N 6.01434W, surface. The arrows show the frequencies of the diurnal and semi-diurnal tides (respectively DT and SDT) and the inertial frequency at 45.3811 • of latitude. The filter impulse response and the amplitude of its transfer function have been scaled by, respectively, 10m s -1 and 0.1m s -1 / √ • h -1 .

Figure 8 :

 8 Figure 8: Test case and frequency transform. The filter impulse response and the amplitude of its transfer function have been scaled by, respectively, 10m s -1 and 0.1m s -1 / √ • h -1 .

  08 17:46 11 Symmetric IIR filtering of inertial waves model output, with N2, M2 and S2 removed, and the filtered detided model output at this location are shown on fig. 7. Only frequencies above 10 • h -1 have been shown as meteorological variations below this frequency are not relevant to this study.

		0.10 0.15 0.20 0.25	eastward sea water velocity at 45.3811N 6.01434W, surface	model detided filtered filter
		0.05
		0.00
	(m/s)	0.05
		0.10
		0.15
		0.20
		0.25
		0.30
		days since 15 Dec 2008 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 72 74 0.35
		0.12
		0.10
		0.08
		0.06
	(m s -1
		0.04
		0.02
		( • h -1 ) 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 0.00

Symmetric IIR filtering of inertial waves

Acknowledgements

I am grateful to Gaëlle Herbert for her model outputs. This work was funded by SHOM.