Zhenhua Wang
email: zhenhua.wang@hit.edu.cn

Wentao Tang

Qinghua Zhang
email: qinghua.zhang@inria.fr

Vicenç Puig
email: vpuig@iri.upc.edu

Yi Shen

Zonotopic state estimation and fault detection for systems with time-invariant uncertainties ⋆

Keywords: Set-membership, time-invariant uncertainty, zonotopes, H ∞ design

This paper proposes a robust guaranteed state estimation method with application to fault detection by combining H ∞ observer design with zonotopic analysis for discretetime systems with both time-variant and time-invariant uncertainties. In order to improve the estimation accuracy, based on the H ∞ technique, the observer design is achieved by solving a linear matrix inequality. The main contribution of this method lies in that the time invariance of some uncertainties is considered to reduce the conservatism of interval estimation. Finally, a numerical example is provided to demonstrate the effectiveness of the proposed method.

INTRODUCTION

State estimation is widely studied and practiced in control theory and applications. Due to ubiquitous uncertainties in system models and measurements, state estimation does not converge to true states. To handle uncertainties, many methods have been proposed, such as the Kalman filter addressing stochastic noises [START_REF] Kalman | New methods in Wiener filtering theory[END_REF], the Kalman filter incorporating unknown inputs [START_REF] Keller | Two-stage Kalman estimator with unknown exogenous inputs[END_REF] and H ∞ observer designs [START_REF] Wang | Robust H ∞ observer design of linear time-delay systems with parametric uncertainty[END_REF][START_REF] Darouach | H ∞ observers design for a class of nonlinear singular systems[END_REF]. Recently, state estimation with guaranteed error bounds, based on set-membership computations, has received a lot of attention, since error bounds are useful in some applications, such as fault diagnosis [START_REF] Xu | Actuator-fault detection and isolation based on set-theoretic approaches[END_REF][START_REF] Xu | Set-theoretic methods in robust detection and isolation of sensor faults[END_REF] and reachability analysis [START_REF] Althoff | Computing reachable sets of hybrid systems using a combination of zonotopes and polytopes[END_REF]. There exist several set-membership methods in the literature, e.g. the zonotope-based method [START_REF] Puig | Worst-case state estimation and simulation of uncertain discretetime systems using zonotopes[END_REF][START_REF] Combastel | A state bounding observer based on zonotopes[END_REF][START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF], ellipsoid-based method [START_REF] Maksarov | State bounding with ellipsoidal set description of the uncertainty[END_REF][START_REF] Kurzhanski | Ellipsoidal techniques for reachability analysis: internal approximation[END_REF][START_REF] Chernousko | Ellipsoidal state estimation for dynamical systems[END_REF] and interval observer [START_REF] Bernard | Closed loop observers bundle for uncertain biotechnological models[END_REF][START_REF] Efimov | On interval observer design for time-invariant discrete-time systems[END_REF]). Among these existing set-membership methods, the zonotopebased method has recently received much attention, and ⋆ This work was supported by National Natural Science Foundation of China (Grant No. 61403104, 61773145), State Key Laboratory of Robotics (Grant No. 2014-O04), and the Fundamental Research Funds for the Central Universities (HIT. KLOF.2015.076). some criteria such as P-radius [START_REF] Le | Zonotopic guaranteed state estimation for uncertain systems[END_REF] and F-radius [START_REF] Combastel | Zonotopes and Kalman observers: Gain optimality under distinct uncertainty paradigms and robust convergence[END_REF] have been proposed to improve the estimation accuracy. Although the zonotope-based method is effective in set-membership estimation, like the other aforementioned methods, it has a certain conservatism due to the wrapping effect. More accurate estimation can be obtained by using higher dimensional zonotopes, but at the cost of increasing computation burden. Therefore, it remains an important issue to establish a trade-off between the estimation accuracy and the computational burden in set-membership methods.

This paper addresses the particular case where the considered systems are subject to time-invariant uncertainties, in addition to time varying uncertainties as usually considered in set-membership based state estimators. Typical examples of time-invariant uncertainties are biases in a system model, such as actuator bias and sensor bias. Most existing methods based on set-membership computations do not consider time-invariant uncertainties. Of course, it is possible to treat time-invariant uncertainties as time varying uncertainties in order to apply existing methods, but ignoring the time-invariant nature of the uncertainties increases the conservatism in set-membership computations.

Few studies have been reported on state estimation for systems subject to both time-invariant and time varying uncertainties. In [START_REF] Combastel | On computing envelopes for discrete-time linear systems with affine parametric uncertainties and bounded inputs[END_REF], discrete-time linear systems with both time-varying bounded inputs and timeinvariant bounded parametric uncertainties are considered by using the so-called parameterized families of zonotopes. However, the proposed method in [START_REF] Combastel | On computing envelopes for discrete-time linear systems with affine parametric uncertainties and bounded inputs[END_REF] can only obtain the out-approximation of the reachable state set at a finite given time. [START_REF] Pourasghar | Reduced-order interval-observer design for dynamic systems with time-invariant uncertainty[END_REF] considers a large class of linear dynamic systems subject to time-invariant uncertainties. The proposed solution in this reported work requires a monotonicity property of the designed interval observer. But the conditions for achieving such a design are not always easy to satisfy.

In this paper, we propose a guaranteed state estimation method for discrete-time state-space systems with both time-varying and time-invariant uncertainties. The proposed method consists of two steps: First, an optimized Luenbeger observer is designed for a single state trajectory estimation, then its state estimation error is bounded by set-membership computations based on zonotopes. The Luenbeger observer and the zonotope based error bound computation then constitute a set-membership estimation estimator. For the class of considered systems, sufficient conditions for achieving the design of such an estimator are explicitly formulated. Next, the proposed method is applied to fault detection for systems with bounded disturbances and time-invariant uncertainties. The main contribution of this paper is an algorithm for error bound computation taking into account time invariant uncertainties. By incorporating the information that some uncertainties are time invariant, the conservatism of zonotopic state estimation is significantly reduced.

The structure of the rest part of the paper is as follows: In Section 2, background on zonotopic analysis is introduced. Section 3 presents the problem formulation. The robust zonotopic estimation is presented in Section 4. Section 5 introduces the application of the proposed method in fault detection. Simulation results are provided in Section 6 to demonstrate the effectiveness of the proposed method. Finally, Section 7 presents the conclusion of this paper.

PRELIMINARIES

Definition 1. An s-order zonotope Z is the affine transformation of a hypercube

B s = [-1, +1] s , i.e. Z = p ⊕ HB s = {p + Hz, z ∈ B s }
where p ∈ R n is the center of Z and H ∈ R n×s , which defines the shape of Z, is called the generator matrix of Z. Herein and throughout this paper, we use ⊕ to denote the Minkowski sum.

For the sake of simplicity and readability, we also use Z = ⟨p, H⟩ to denote a zonotope.

In the sequel, the following properties of zonotopes will be used:

⟨p 1 , H 1 ⟩ ⊕ ⟨p 2 , H 2 ⟩ = ⟨p 1 + p 2 , [H 1 H 2 ]⟩ (1a) L ⊙ ⟨p, H⟩ = ⟨Lp, LH⟩ (1b) ⟨p, H⟩ ⊆ ⟨p, H⟩ (1c)
where p, p 1 , p 2 ∈ R n and H, H 1 , H 2 ∈ R n×s denote the centers and generator matrices of zonotopes, respectively, L ∈ R l×n is a matrix with appropriate dimensions, ⊙ denotes linear map, and H ∈ R n×n is a diagonal matrix as follows:

H =          s ∑ j=1 |H 1,j | • • • 0 . . . . . . . . . 0 • • • s ∑ j=1 |H n,j |          ,
and ⟨p, H⟩ is also a box.

In the applications of zonotopes, the reduction operator will be used to bound a high dimensional zonotope by a lower one, which can be described by the following property [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF] and is originally proposed by [START_REF] Combastel | A state bounding observer based on zonotopes[END_REF].

Property 1. Given the zonotope Z = p ⊕ HB s ⊂ R n , reorder the columns of matrix H in decreasing Euclidean norm and denote the obtained matrix as Ĥ. Then Z ⊆ p⊕ R q (H)B q (n ≤ q ≤ s), where R q (•) is the reduction operator that can reduce the order of Z from s to q that is a chosen reduced dimension. And

R q (H) = [H T Q],
where H T consists of the first q -n columns of Ĥ and Q is a diagonal matrix with

Q i,i = s ∑ j=q-n+1 | Ĥi,j |, i = 1, . . .

, n

For a symmetric matrix P ∈ R n×n , P ≻ 0(P ≺ 0) denotes that P is positive(negative) definite. In the following, the comparison operator ≤ and the absolute value operator | • | are understood elementwise. In a symmetric block matrix, we use * to represent a term that can be induced by symmetry.

PROBLEM FORMULATION

In this paper, we consider the following discrete-time system

{ x k+1 = Ax k + Bu k + Gθ + Ew k y k = Cx k + F v k (2)
where We assume that x 0 , w k , v k and θ are unknown but bounded, i.e.

x k ∈ R nx , u k ∈ R nu ,
|x 0 -p 0 | ≤ x, |w k | ≤ w, |v k | ≤ v, |θ| ≤ θ (3) where p 0 ∈ R nx , x ∈ R nx , w ∈ R nw , v ∈ R nv and θ ∈ R n θ are known vectors.
According to the definition of zonotopes, (3) can be reformulated as

x 0 ∈ ⟨p 0 , H 0 ⟩, w k ∈ ⟨0, W ⟩, v k ∈ ⟨0, V ⟩, θ ∈ ⟨0, Θ⟩ (4) where H 0 = diag(x), W = diag(w), V = diag(v) and Θ = diag(θ).
The aim of this paper is to design a robust guaranteed state estimator for system (2). In this paper, we integrate robust observer design with zonotopic analysis to achieve guaranteed state estimation. The estimator to be designed is based on the observer for system (2) with the following form:

xk+1 = Ax k + Bu k + L(y k -C xk )
(5) where xk ∈ R nx is the state estimation and L ∈ R nx×ny is the gain matrix to be designed. In order to increase the estimation accuracy, we first use H ∞ technique to design L such that the estimation error is robust to uncertainties. Then, based on the obtained observer, guaranteed state estimation is achieved via zonotopic analysis. Moreover, the proposed method is applied to fault detection for systems with bounded disturbances and time-invariant uncertainties.

ROBUST ZONOTOPIC STATE ESTIMATION

In this section, we aim to obtain robust state estimation described by zonotope sets, which are consistent with the initial state set, measurements and uncertainties.

Robust observer design based on H ∞ technique

In order to attenuate the effects of uncertainties and thus to increase the estimation accuracy, we use the well-known H ∞ technique to design the observer gain L for observer (5). To this end, we first analyze the error dynamics of observer (5).

Define the estimation error as e k = x k -xk . By subtracting (5) from (2), we obtain the following error dynamics: e k+1 = (A -LC)e k + Gθ + Ew k -LF v k (6) that can be rewritten as e k+1 = (A -LC)e k + Dd k (7) where

D = [G E -LF ] , d k = [ θ w k v k ] .
We denote the transfer function from d k to e k as G ed (z). The H ∞ norm of G ed (z), i.e. its maximum singular value, is denoted as ∥G ed (z)∥ ∞ .

Based on (7), we propose the following theorem to design a robust observer, which is based on H ∞ technique. Theorem 1. Given a scalar γ > 0, the error system ( 7) is stable and satisfies the following performance index

∥G ed (z)∥ ∞ < γ (8) if there exist a symmetric positive definite matrix P ∈ R nx×nx (P ≻ 0) and a matrix W ∈ R nx×ny such that      0 0 0 -γ 2 I ⋆ P A -W C P G P E -W F -P      ≺ 0, (9)
Moreover, if the LMI in ( 9) is solvable, the gain matrix L can be determined by L = P -1 W .

Proof. Applying the well-known bounded real lemma to error dynamics (7), we know that ( 7) is stable and satisfies the H ∞ performance index in (8), if and only if there exists a symmetric positive definite matrix P ≻ 0 and a matrix L such that

[ (A -LC) T P (A -LC) -P + I ⋆ D T P (A -LC) D T P D -γ 2 I ] ≺ 0.
(10)

Note that ( 10) is not a standard LMI, we use Schur complement lemma to convert (10) as

  -P + I ⋆ ⋆ 0 -γ 2 I ⋆ P (A -LC) P D -P   ≺ 0. ( 11
) By substituting D = [G E -LF ] into (11), we obtain      0 0 0 -γ 2 I ⋆ P (A -LC) P G P E -P LF -P      ≺ 0. ( 12
)
Letting W = P L, (12) becomes (9). 2 Remark 1. To optimize the estimation accuracy, the minimal γ can be obtained by solving the following optimization problem: min γ 2 , s.t. (9). ( 13)

and the corresponding observer gain can be obtained by L = P -1 W .

Guaranteed state estimation based on zonotopes

Based on the observer (5) designed via H ∞ technique, we can obtain guaranteed state estimation described by zonotopes sets, which can enclose all the admissible values of state. We note that the Zonotopic Kalman Filter(ZKF) [START_REF] Combastel | Zonotopes and Kalman observers: Gain optimality under distinct uncertainty paradigms and robust convergence[END_REF] is based on a similar observer structure with that of (5). The ZKF for system (2) can be described by the following lemma. Lemma 1. For system (2), if

x k ∈ ⟨c k , H k ⟩, then x k+1 ∈ ⟨c k+1 , H k+1 ⟩ with c k+1 = Ac k + Bu k + L k (y k -Cc k ) H k+1 = [ (A -LC) Hk , G, E, -LF ] , Hk = R q (H k ). (14) Moreover, L k can be obtained by L k = A Hk HT k C T (C Hk HT k C T + F F T ) -1 (15)
Proof. Lemma 1 can be easily proved by applying Proposition 1 and Theorem 5 in [START_REF] Combastel | Zonotopes and Kalman observers: Gain optimality under distinct uncertainty paradigms and robust convergence[END_REF] to system (2). Hence, readers can refer to this reference for more details about the proof. 2

We note that the time invariance of some uncertainties is not considered in ( 14), which will introduce conservatism. Moreover, the observer gain L k is computed by ( 15) on line, which may cause heavy computational burden since the dimensions of Hk may be considerably large. In order to reduce conservatism, we consider the time invariance of θ via the reachability analysis of the estimation error and then obtain more accurate guaranteed state estimation than ZKF. To this end, we first introduce the following useful lemma. Lemma 2. For a discrete-time linear system as follows

x k+1 = Ax k + Ew k (16) where x 0 ∈ X 0 = ⟨p 0 , H 0 ⟩ and w k ∈ W = ⟨0, W ⟩, the state x k belongs to zonotope Xk = ⟨p k , Ĥk ⟩ with p k+1 = Ap k , (17) Ĥk+1 = [ AR q ( Ĥk ) EW ] (18)
Proof. Assuming x k ∈ X k , based on ( 16), X k+1 is updated as follows:

X k+1 = A ⊙ X k ⊕ E ⊙ W (19) Denote X k = ⟨p k , H k ⟩.
Using the properties of zonotope in ( 1), ( 19) can be rewritten as:

p k+1 = Ap k , H k+1 = [AH k EW ] .
According to Property 1, we have ⟨p k , H k ⟩ ⊆ ⟨p k , R q ( Ĥk )⟩ and it follows that ⟨p k+1 , H k+1 ⟩ ⊆ ⟨p k+1 , Ĥk+1 ⟩. Therefore, we have X k ⊆ Xk . 2

Based on Lemma 2, we propose the following theorem to achieve guaranteed state estimation for system (2), where the time invariance of θ is considered. Theorem 2. For system (2), under the condition (4), x k can be bounded by a zonotope Xk as follows:

Xk = ⟨x k , Ĥk ⟩ ( 20
)
where xk is obtained by ( 5) with x0 = p 0 and Ĥk has the following form:

Ĥk = [ Hk Hk ] . ( 21
)
with Hk and Hk obtained from

Hk = k-1 ∑ i=0 (A -LC) i GΘ (22) and Hk+1 = [ (A -LC)R q ( Hk ) EW -LF V ] (23)
and H0 = H 0 .

Proof. Since e k = x k -xk , it follows that x k = xk + e k Since xk can be obtained from (5), the guaranteed state estimation can be achieved via estimating the reachable set of e k . In the following, we first obtain the guaranteed estimation of e k , and then give that of x k .

Motivated by [START_REF] Wang | H -/H ∞ fault detection observer in finite frequency domain for linear parameter-varying descriptor systems[END_REF], we split the error dynamics in (6) into two subsystems as follows:

ēk+1 = (A -LC)ē k + Gθ (24) and ẽk+1 = (A -LC)ẽ k + Ew k -LF v k (25
) where e k = ēk + ẽk with ē0 = 0 and ẽ0 = e 0 . It is obvious that ēk is only affected by θ and ẽk is driven by w k and v k . By using (24) iteratively, we obtain

ēk = k-1 ∑ i=0 (A -LC) i Gθ (26)
Note that θ ∈ ⟨0, Θ⟩, (26) follows that ēk ∈ Ēk = ⟨0, Hk ⟩ (27) where Hk is given by ( 22).

On the other hand, applying Lemma 2 to (25) yields ẽk ∈ Ẽk = ⟨0, Hk ⟩ (

where Hk is given by ( 23).

Since e k = ēk + ẽk , we obtain

e k ∈ Êk = Ēk ⊕ Ẽk = ⟨0, Ĥk ⟩ (29)
where Ĥk = [ Hk Hk ] .

Finally, according to x k = xk + e k , we have

x k ∈ Xk = xk ⊕ Êk = ⟨x k , Ĥk ⟩, ( 30
) which ends the proof. 2 Remark 2. From (1c), we know that a zonotope can be bounded by a box. Then, we can obtain the interval estimation of the state including the upper bound x k and lower bound x k , which can be obtained from

           x k (i) = xk (i) + s ∑ j=1 | Ĥk (i, j)|, i = 1, . . . , n x k (i) = xk (i) - s ∑ j=1 | Ĥk (i, j)|, i = 1, . . . , n (31)
where s is the column number of Ĥk .

ZONOTOPIC FAULT DETECTION

When system (2) is affected by additive actuator faults, it can be formulated as follows

{ x k+1 = Ax k + Bu k + Gθ + Ew k + Bf k y k = Cx k + F v k (32)
where f k ∈ R nu represents the actuator fault vector.

The state estimations obtained by ZKF and the proposed method are both based on (2) without consideration of actuator fault.

When no fault occurs, the state estimation based on zonotopes can guarantee to enclose the admissible state values, i.e. x k ∈ Xk . Since

y k = Cx k + F v k , it follows that y k ∈ C Xk ⊕ F ⟨0, V ⟩ (33)
According to the properties of zonotopes, (1), we have

y k ∈ Ŷk = ⟨p y k , H y k ⟩ (34) where { p y k = Cp k H y k = [ C Ĥk F ] (35)
However, when a fault occurs and the magnitude of it is big enough to be detected, y k will fall outside Ŷk . Therefore, we propose a zonotopic fault detection strategy:

If y k / ∈ Ŷk , it can be inferred that a fault has occurred.

An interval estimation enclosing Ŷk can be obtained similar to (31) for the simplicity of fault detection. Remark 3. The zonotopic fault detection can also be achieved based on the state estimation obtained by ZK-F. However, the fault detection based on the proposed method can achieve higher fault detection rate (FDR) than that based on zonotopes, since the time-invariance of θ is considered in the proposed method.

6. SIMULATIONS

State estimation results

In this section, an example of the decoupled linearized longitudinal motion dynamics of the F-18 aircraft adapted from [START_REF] Yang | Adaptive fault-tolerant H ∞ control via state feedback for linear systems against actuator faults[END_REF] is used to demonstrate the effectiveness of the proposed method. The state space equation of the consider system is as follows

[ α q ] = A m 7h14 long [ α q ] + B m 7h14 long [ δ E δ P T V ] + w, ( 36
)
where α is the angle of attack, q is the pitch rate, δ E si the symmetric elevator position, δ P T V is the symmetric pitch thrust velocity nozzle position and w is the external disturbance. A m 7h14 long and B m 7h14 long denote the system matrices at Mach 0.7 and 14-kf t altitude. q is assumed to be measurable. By using the Euler method to discretize the continuous-time model ( 36) with sample time T s = 0.1s, we obtain a discrete-time system in the form of ( 2

k | ≤ [1 1] T , |v k | ≤ 1.
The system suffers an unknown constant actuator bias, which can be described by G = B and θ ∈ ⟨0, 0.03I 3 ⟩.

The initial sate is

x 0 = [0.2 -0.2] T . The uncertainties are set as p 0 = [0 0] T , H 0 = 0.5I 2 , W = I 2 , V = 1
and let x0 = p 0 .

By solving the optimization problem (13), we obtain the observer gain as

L = [ -0.0051 1.0117 ]
with the minimal γ = 0.0646.

In the simulation study, the proposed method is compared with ZKF. The simulation results are shown in Figure 1-2. They show that the guaranteed state estimations obtained by the proposed method are more accurate than those by ZKF, which is due to that the proposed method has considered the time-invariance of θ.

Fault detection results

In order to demonstrate the effectiveness of the zonotopic fault detection based on the proposed method, we consider a case where an actuator fault occurs, which can be simulated as follows

f k = { [0, 0] T , k < 50; [0.15, 0] T , k ≥ 50.
The simulation results are shown in Figure 3. It shows that y k quickly exceeds the interval estimation obtained based on the proposed method after the occurrence of fault, which indicates that a fault has been detected. On the other hand, it also shows that the interval estimations based on ZKF has lower fault detection rate because of the extra conservatism introduced by neglecting the time invariant nature of the uncertainties.

CONCLUSION

In this paper, a novel guaranteed state estimation method is proposed for the discrete-time linear systems affected by both disturbances and time-invariant uncertainties.

In order to increase the estimation accuracy, the proposed method uses the H ∞ technique to design robust observer. Based on the designed observer, guaranteed state estimation is obtained via zonotopic analysis. The proposed method has considered the time dependency of uncertainties to reduce conservatism. Moreover, the proposed method is applied to actuator fault detection to increase FDR. The simulation results show that the proposed method can obtain more accurate estimation than ZKF, which does not consider explicitly such timeinvariance. The fault detection approach based on the proposed method is shown able to obtain higher FDR than that based on ZKF.

  01I 2 , C = [0 1] , F = 0.01. The disturbance w k and measurement noise v k are bounded by |w

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. x k (1) and its interval estimations by ZKF and the proposed method.

-P + I ⋆ ⋆ ⋆ ⋆ 0 -γ 2 I ⋆ ⋆ ⋆ 0 0 -γ 2 I ⋆ ⋆


Guillaume Becq

Stéphane Bonnet

Lorella Minotti

Michel Antonakios

Régis Guillemaud

Philippe Kahane

Collection and Exploratory Analysis of Attitude Sensor Data in an Epilepsy Monitoring Unit

Keywords: epilepsy, movement, seizures, attitude sensors, feature extraction, principal component analysis, monitoring, accelerometer, magnetometer

The aim of this paper is to present the collection of attitude sensor data from an epilepsy monitoring unit and the results of standard exploration using principal component analysis. The collection of data from attitude sensors positioned on three limbs of epileptic patients at their bedside is described. The analysis of the data focuses, on one hand, on motor features extraction from attitude sensor data and on the other hand, on visual segmentation of seizures into events corresponding to motor manifestations classes by an expert. Principal component analysis is then realized over these features and groups of data are localized according to the expert classification. This exploration indicates a possible discrimination between these motor manifestation classes.

I. INTRODUCTION

E PILEPSY is one of the most important neurologi- cal disorders affecting both adults and children [START_REF] Thomas | Epilepsies[END_REF]. Epilepsy is characterized by the spontaneous repetition of epileptic seizures which denote brain dysfunction. During an epileptic seizure, an excessive electrical discharge, arising from one brain area, may propagate into other brain areas and induce different clinical manifestations. The chronological organization of these ictal clinical manifestations is actually reflecting the electrical pathway during the discharge. For instance, frontal epileptic seizures are often characterized by the occurrence of one-sided or two-sided motor manifestations of limbs, head or body axis.

The semiological seizure analysis is usually performed in an Epilepsy Monitoring Unit (EMU) by recordings and visualizing seizures with video-EEG recordings. This is essential for the identification of epileptogenic areas in the brain and the comprehension of different clinical and EEG patterns.

Although EEG recording is the gold standard to analyze epileptic seizure, video gives also many clues for the diagnostic. However, video essentially presents qualitative information, which is often hard and long to process. Moreover data visual inspection by EEG technicians is rather subjective and does not allow comparing between different motor patterns. Several solutions have been proposed recently to Asterisk indicates corresponding author. This work was supported by the ITS 2006. G. becq, L. Minotti and P. Kahane are with Inserm U704, Dynamic of neural networks, in the CHU Grenoble, BP217, 38043 Grenoble cedex 9, France guillaume.becq@laposte.net) S. Bonnet, M. Antonakios and R. Guillemaud are with the CEA/ Leti/DTBS/STD/Lisa, 17 rue des martyrs, 38054 Grenoble cedex, France (stephane.bonnet@cea.fr) quantify movement in an automated, unsupervised and realtime manner.

Different authors have proposed to detect and diagnose epileptic seizures using video analysis with or without markers fixed on human body [START_REF] Karayiannis | Quantifying motion in video recordings of neonatal seizures by robust motion trackers based on block motion models[END_REF]- [START_REF] Li | Movement quantification in epileptic seizures: a new approach to video-EEG analysis[END_REF]. This approach is attractive because it can be based on the routine video-EEG procedure performed in every EMU but a major limitation is that the person (or the markers) should remain in the field-of-view of the camera.

A much powerful approach is based on body-mounted sensors like accelerometers or rate gyroscopes [START_REF] Frost | Triaxial accelerometry: a method for quantifying tremor and ataxia[END_REF]- [START_REF] Salarian | An ambulatory system to quantify bradykinesia and tremor in parkinson's disease, 4th ann[END_REF]. MEMS technology is well-suited to monitor epilepsy since it allows quantifying intensity, frequency and duration of movements. Furthermore inertial sensors are well adapted for applications of ambulatory measurement of human movements: high integration, small weight and low cost. As a consequence they can be used outside the clinical setting.

In order to obtain a full characterization of 3D epileptic movements, we propose to combine magnetic sensors with accelerometers into the same attitude sensor module [START_REF] Bonnet | A magnetometer-based approach for studying human movements[END_REF]. Magnetometers are primarily intended to characterize movement in the horizontal plane as it can be observed during tonic episodes. Furthermore, these sensors are also sensitive to tremors while being insensitive to accelerations.

This paper aims at showing that quantitative analysis of wrist and head movements, by means of attitude sensors, provides objective measures for the characterization of epileptic seizures and reveals that an automated classification of motor manifestations could be possible. This is done by using principal component analysis (PCA).

In the first section, we describe the clinical environment of this study and then we present the different features that were extracted to characterize epileptic motor manifestations occurring during seizures. Preliminary results are proposed in the last section.

II. MATERIALS

The EMU in CHU Grenoble is composed of two bedrooms, each one equipped with Brain Quick video-EEG acquisition system from Micromed [START_REF] Micromed S.R.L | via Giotto[END_REF]. It is composed of an amplifier accepting up to 128 leads with selectable sampling rate from 256 to 1024 Hz (16 bits resolution). An optical link insures data transmission to a computer in the monitoring room. Audios-videos and EEG recordings are synchronized with Micromed System Plus patient data management system software. In addition to the daily routine, a protocol was elaborated with the medical staff to collect motion information during seizure from body-mounted attitude sensors. To capture upper limb activities, three sensor modules, each containing a tri-axial accelerometer and a tri-axial magnetometer, were used. These attitude sensors (dimension: 50 mm x 30 mm x 15 mm) were mounted on both wrist and head segments using straps (Fig. 1). The fabric is washable and well supported by patients.

Every module was linked by a 4 m wire to a data acquisition system, a so-called Trident system, positioned at the patients bedside. All signals were multiplexed and sampled at 46 Hz (16 bits resolution) using the Trident system. The energy was supplied by an AC-DC voltage transformer with medical safety approval. Data (18 digital signals) were redirected to a computer in the monitoring room via RS-232 link for recording. An opto-coupling device was introduced to isolate the patient from the computer electrical source coming from the hospital network.

In order to synchronize the video-EEG system and the Trident system, a synchronization was designed. This unit collected a TTL signal (0-5 V , 1 Hz) from the Micromed system as it is currently done currently to drive intermittent photic stimulation (IPS). The synchronization unit provided then a dry contact to the Trident system and this binary information was stored altogether with data.

In practice, the technical staff should activate the IPS at the beginning and end of the recording session. The diagram of the data acquisition process is illustrated in Fig. 2. Finally, video-EEG data and motion data were merged in an offline process into a single file by dealing with time synchronization and resampling.

Data were recorded on 17 patients, 12 women, 5 men, aged from 8 to 45 years during 211 half days of 3 h recordings session. Within this period, only 9 patients generated motor seizures with a total of 29 seizures. This corresponds to a ratio of approximately 1.4 episodes per week of working days.

III. METHODS

A. Motor manifestation classification

Seizures were visually detected by the medical staff using video EEG and temporally defined by an electrical beginning t 0 , a clinical beginning t 1 and an electrical end t 2 as in Fig. 3. Each seizure was visually analyzed by an expert in order to segment each movement from each limb into an event associated to a class of motor manifestation. For each event, we used the classification proposed in [START_REF] Lüders | Semiological Seizure Classification[END_REF] as a reference and simplified it into: No movements (NOMVT) no movements observed; Automatisms (AUTO), repetitive automatic movements of limbs and trunks; Clonic manifestations (CLONIC), successive contractions of agonist and antagonist muscles (including myoclonies and clonies); Tonic manifestations (TONIC), sustained contraction of muscle resulting or not in a change of position; Tonico-clonic manifestations (TC), movements combining both clonic and tonic manifestations; Hypermotor manifestations (HYPER), complex movements more or less rapid, involving the proximal segment of limbs or trunk; Versive manifestations (VERSIVE), deviation of the head with or without eyes deviation; Others manifestations (OTHERS), unclassified manifestations (e.g. grinning, paining).

For this study VERSIVE was considered as TONIC, TC was consider as CLONIC and data of the OTHERS class were discarded. A number of artifacts generated by the medical staff for clinical evaluation or patient's care were removed for data cleaning.

B. Feature extraction 1) General concept: AS data coming from the different events were collected into a database. They follow a feature extraction process in order to obtain for each event a vector containing a set of angular, dynamic, spectral and wavelet features supposed to be relevant for each class an possibly useful for classification. A block diagram of the feature extraction process is presented Fig. 4.

The principal asset of our method based on AS data is the evaluation of the principal Euler angle and features extraction from it.

2) Principal Euler angle computation: Using both vector observations Acc /b and M ag /b and combining them into AS /b , it is possible to estimate the sensors attitude with Fig. 3. Data from AS and EEG are merged together. EEG enables an expert to localize the seizure defining: t 0 the electrical beginning of the seizure, t 1 the clinical beginning of the seizure, and t 2 the end of the seizure. Visual inspection with video enables an expert to segment the seizure into different events. Data presented here come from the left wrist. Three events corresponding to tonic, automatism and others manifestation are succeeding one after another. respect to a frame of reference [START_REF] Bonnet | A magnetometer-based approach for studying human movements[END_REF] 1 . In this study, we chose to use the Euler angle/axis representation: a rotation from the reference frame to the current frame is performed about the unit axis u through angle θ. A sample to sample or instantaneous rotation was determined using:

θ(t i ) = f Θ ( AS /b (t i ), AS /b (t i-1
)) with t i = t s + i/F s , t e = t s + (N -1)/F s , t s start time of the event, t e end time of the event, and f Θ an estimation function based on least min squared minimization [START_REF] Bonnet | A magnetometer-based approach for studying human movements[END_REF].

3) Angular features: Angular features like in [START_REF] Salarian | An ambulatory system to quantify bradykinesia and tremor in parkinson's disease, 4th ann[END_REF] using gyroscopes for measuring rotation rates were computed: the root mean square of the signal, referred as the mobility by Salarian, M obility =

1 (N -1)
N -1 i=1 θ 2 (t i ); from the cumulated rotation Θ(t j ) = j i=1 θ s (t i ) by removing instantaneous rotations θ s with angular velocities θ(t) * F s fewer than 10 o / s, the median value med Θ and the range of values Range.

4) Dynamical features: We used the norm of acceleration Acc (t) to compute dynamical features as a good candidate to discriminate movements versus no movements Acc (t) = Acc x (t) 2 + Acc y (t) 2 + Acc z (t) 2 . The entropy of the norm of acceleration, En Acc = i h i log(h i ), was also retained with a quantification of 16 bits between 0 and 3 G to compute h i = n i /N , with n i 1 The measurement model for the accelerometer is given by: Acc = gγ + noise. The components of this vector quantity are resolved in the bodyfixed frame: Acc /b = (Accx, Accy, Accz) t . The measured acceleration signal is a combination of a low-frequency component due to the sensor orientation with respect to gravity ( g) and a high-frequency component due to the sensor movement itself (γ). For slow movements, the accelerometer is principally acting as an inclinometer and the norm of the accelerometer is close to 1 g. The separation between the two components is approximately achieved by low-pass filtering (Butterworth, 4th order, cut-off frequency: 0.5 Hz). The measurement model for the magnetometer is simplified into: M ag = M ag earth + noise. Where M ag earth denotes the Earth magnetic field vector. The components of this vector quantity are resolved in the bodyfixed frame: M ag /b = (M agx, M agy, M agz) t .

} AS t = { Acct   x y z  , Mag t   x y z } Acc t   x y z  Mag t   x y z  Acc LF  t  Norm Norm Acc t 

Low pass

For each event For all samples at time t From to t s t e F1,,F d  t

Vector of features associated to this event Fig. 4. Feature extraction block diagram. Data from accelerometer Acc(t) and from magnetometer M ag(t) are grouped together to form an attitude sensor vector AS(t). Filtering, Euler angles computation and norm computation leads to new variables for features extraction. This is done for data of each event, resulting in one vector of features associated to each event.

the number of samples at level q i (n i = #( Acc (t) = q i )), N the number of samples in the event and h i log (h i ) = 0 for h i = 0.

5) Spectral features: Spectral features were evaluated using the spectral edge frequency (SEF), indicating the frequency value for which α % of the spectral power is obtained for 10 and 95 % (respectively SEF 10 and SEF 95 ). Features for Θ and medio-lateral axis component for Acc and M ag (Acc ml , M ag ml ) were computed.

6) Wavelet features: Wavelet features were evaluated using a ratio of wavelet coefficients in different scales bands. Coefficients came from a continuous wavelet transform, C x (a, b) = ∞ -∞ x(t) ψ(a, b, t)dt, a being a scale parameter and b a lag parameter, using a Morlet atom wavelet ψ(t) ∼ exp (-x 2 /2) cos(5x), with the Matlab wavelet toolbox [12]. The different scales retained where chosen to correspond to the pseudofrequencies grouped into three different bands δ 3 = {23, 20, 17, 14, 12, 10, 9} Hz, δ 2 = {8, 7, 6, 5, 4, 3, 2} Hz, δ 1 = {1.25, 1, 0.75, 0.5, 0.25, 0.1, 0.05} Hz. Features for Acc ml , M ag ml and Θ were computed.

C. Principal component analysis

Features data were standardized using the transformation Y j = Xj -µj σj with X j vectors made from the samples of the jth feature, µ j and σ j the mean and standard deviation of these samples. 30 uniform random samples from each sample were introduced into the database for data regularization. A principal component analysis was realized in order to find the best linear combinations of all features retaining the best singular values of the correlation matrix. This often results in better separation of data in the new basis and indicates possible discrimination. V. DISCUSSION Using AS data and applying PCA indicates that an automatic classification of features extracted from this data into motor manifestations relative to epilepsy can be possible.

IV. RESULTS

manifestations events

The discrimination between classes suffers principally from the proximal nature of some manifestations and on the weak discriminative power of some retained features that must be investigated.

We expect to learn classifiers with the best features. This will allow us to obtain automatic succession of motor manifestations classes and generate graphs as a tool for neurologists to characterize the stereotaxic evolution of motor manifestations during seizures.

Fig. 1 .

 1 Fig. 1. Attitude sensors (AS) set on patients in an epilepsy unit. On the left, a silicon rectangular parallelepiped was first set on a wrist with adhesive tape. For patient and staff convenience, special dedicated self gripping belts were designed and set on the left wrist (middle), right wrist (not shown) and on the head (right).

Fig. 2 .

 2 Fig.2. AS data coming from the different locations are recorded on an isolated PC thanks to the Trident System conversion. In order to synchronize EEG -video data from the Micromed system and AS data, an intermittent photic stimulation (IPS) signal is routed to be used as a marker present in both recordings. Data synchronization is made offline with own dedicated softwares.

  {AS t -1|AccLF  , AS t|Acc LF } Compute Euler angles Angular Spectral Dynamic Wavelet {Euler Principal angle , Euler Principal axis

  were scored, representing in each class: 42 NOMVT, 25 AUTO, 17 CLONIC, 70 TONIC, 21 TC, 20 HYPER, 11 VERSIVE, 20 OTHERS.

Fig. 5 .

 5 Fig. 5. Data projection on the four first principal components. In the diagonal, figure indicates histograms of data for this component. Each data corresponding to one event is represented by its new coordinates (scores) in the principal components basis. A color indicates the motor manifestation class of the event. 30 random samples from each original sample were generated to simulate densities. 3D scatters in a 3 first components basis are represented on the right top with different fields of view.

ACKNOWLEDGMENT

The authors would like to thank the medical staff, Patricia Boschetti, Véronique Dorlin, Martine Juillard and the patients that participate to the study.