Claire Paya

Euriell Le Corronc

Yannick Pencolé

Philippe Vialletelle

Observer-based detection and localization of time shift failures in (max,+)-linear systems*

This paper addresses a problem of failure detection and localization in production lines modeled, through Timed Event Graphs (TEG), as (max,+) linear systems with disturbances, over which observers can be developed. The state of the observed system is estimated and an indicator that returns true if a time shift failure is detected is defined. The localization step is proposed for elementary structures of TEG through the results of the detection.

I. INTRODUCTION

In the industry, Discrete Event Systems (DES) can be used to model and diagnose faults, malfunctions on production lines. The objective is to detect, locate and identify failures as soon as possible to avoid further equipment unavailability. At STMicroelectronics Crolles 300 plant, we investigate this problem in order to detect production drifts, especially failures that generate time shifts in the production lines. STMicroelectronics is among the world's largest semiconductor companies, serving all electronics segments. Semiconductor manufacturing is a complex industries and one of its most important challenges is to succeed in detecting production drifts before they have real impact on production plans.

Fault diagnosis involving time in DES has been introduced in [START_REF] Tripakis | Fault diagnosis for timed automata[END_REF] using timed automata to refine diagnosis decisions based on timed observations. In [START_REF] Ghazel | State observer for des under partial observation with time petri nets[END_REF], time Petri nets allow to easily model concurrency. This paper represents production lines as Timed Event Graphs (TEG) which are a subclass of Petri nets where places are associated with a duration. A TEG can be modeled by (max,+) algebra as introduced in [BCOQ92], [Max91], [START_REF] Komenda | Max-plus algebra in the history of discrete event systems[END_REF]. For example [START_REF] Kim | Feedback control of cluster tools for regulating wafer delays[END_REF] uses (max,+) algebra techniques to control wafer delays in cluster tools for semiconductor production.

In this paper, we propose an observer-based (max,+) algebra method that detects and then localizes the potential sources of time shifts in a TEG. Our proposal extends the following results. [LCPSP21] presents a detection method of time shift failures by simply comparing simulated and real outputs without estimating the internal state of the system. Once the detection made, the localization proposes a -possibly large -set of admissible sources of the failure using signature matrices and characteristic signatures. To improve the precision of the detection, [START_REF] Paya | Observer-based detection of time shift failures in (max,+)linear systems[END_REF] uses a (max,+) observer [START_REF] Hardouin | Observer design for (max,+) linear systems[END_REF] to actually compare the estimated state with a fault-free simulated state. Our proposal goes one step further by providing suspected sources of the time shifts in a more precise way, by analysing elementary structures of TEG (tandem, synchronization and parallelism). Using (max,+) algebra makes this work of polynomial complexity.

The paper is organized as follows. Section II recalls the necessary mathematical background. Section III refines the detection method detailed in [START_REF] Paya | Observer-based detection of time shift failures in (max,+)linear systems[END_REF] that will be used to setup the rules for localizing the source of time shift failure in TEG defined in Section IV.

II. MATHEMATICAL BACKGROUND

This section recalls the mathematical background used for describing (max,+)-linear systems [START_REF] Baccelli | Synchronization and linearity: an algebra for discrete event systems[END_REF], [Max91].

A. Dioid theory

Dioid theory is the mathematical framework for modeling Timed Event Graphs (TEG) as (max,+) algebra. A dioid D is a set composed of two internal operations ⊕ and ⊗. The addition ⊕ is associative, commutative, idempotent (i.e. ∀a ∈ D, a⊕a = a) and has a neutral element ε. The multiplication ⊗ is associative, distributive on the right and the left over the addition ⊕ and has a neutral element e. Element ε is absorbing by ⊗. When there is no ambiguity, the symbol ⊗ is omitted. A dioid is complete if it is closed for infinite sums and if ⊗ is distributive over infinite sums.

For instance, the set Z max = (Z ∪ -∞), endowed with the max operation as addition ⊕ and the addition as multiplication ⊗ with neutral element ε = -∞ and e = 0 is a dioid. By adding +∞ to the dioid Z max , we get the complete dioid Z max where (-∞) + (+∞) = (-∞).

Definition 1: For a dioid D, denotes the order relation

such that ∀a, b ∈ D, a b ⇔ a ⊕ b = b.
Theorem 1 ([BCOQ92]): Let D be a complete dioid, x = a * b is the solution of x = ax ⊕ b where a * = i≥0 a i is the Kleene star operator with a 0 = e and a i+1 = a ⊗ a i .

To model TEG as (max,+)-linear systems, a specific dioid has to be defined: the dioid M ax in [[γ, δ]]. It is a particular dioid of formal series with two commutative variables γ and δ that will represent event and time shifts of a TEG as explained in the next subsection. This dioid comes from the following

B[[γ, δ]] dioid. The complete dioid B[[γ, δ]]
is the set of formal series with two commutative variables γ and δ with Boolean coefficients in {ε, e} and exponents in Z. A series s ∈ B[[γ, δ]] is written s = (n,t)∈Z 2 s(n, t)γ n δ t where s(n, t) = e or ε (respectively representing the presence or the absence of the monomial). The neutral elements are ε(γ, δ) = (n,t)∈Z 2 εγ n δ t and e(γ, δ) = γ 0 δ 0 . Graphically, a series of B[[γ, δ]] is described by a collection of point of coordinates (n, t) in Z 2 with γ as horizontal axis and δ as vertical axis. For instance, Figure 1 shows series u 1 = u 2 = γ 0 δ1 ⊕ γ 1 δ 2 ⊕ γ 2 δ 3 ⊕ γ 3 δ 4 ⊕ γ 4 δ 5 ⊕ γ 5 δ 6 ⊕ γ 6 δ 7 ⊕ γ 7 δ +∞ (monomials with e as Boolean coefficient). The complete dioid

M ax in [[γ, δ]] is the quotient of B[[γ, δ]] modulo γ * (δ -1 ) * where ∀a, b ∈ M ax in [[γ, δ]]: a = b ⇔ aγ * (δ -1 ) * = bγ * (δ -
[[γ, δ]] is s = K k=0 γ n k δ t k with K ∈ N∪{+∞} and n 0 < n 1 < . . . t 0 < t 1 < . . .

B. Models of (max,+)-linear systems

Figure 2 presents a TEG with two inputs u 1 and u 2 and one output y 1 . Durations associated to places are expressed here by letters a, b, c, d, f and transition x 3 is a synchronization between paths from u 1 and u 2 . The entire structure of such a TEG can be modeled by equations in M ax in [[γ, δ]] through a set of matrices A, B and C. Then, relations between input u, state x and output y transitions are expressed by its state representation.

Let u ∈ M ax in [[γ, δ]] p×1 , x ∈ M ax in [[γ, δ]] n×1 and y ∈ M ax in [[γ, δ]] q×1
, the state representation of a TEG is:

x = Ax ⊕ Bu, y = Cx, where A ∈ M ax in [[γ, δ]] n×n , B ∈ M ax in [[γ, δ]] n×p and C ∈ M ax in [[γ, δ
]] q×n . Equality x = Ax ⊕ Bu can be transformed to x = A * Bu thanks to Theorem 1 so we have y = CA * Bu.

Matrix H = CA * B represents the transfer function that is the dynamic of the system between inputs and outputs.

Example 1: In the TEG of Figure 2, suppose that the durations are a = b = d = 1 for places p 1 , p 3 , p 4 , c = 2 for p 2 and f = 0 for p 5 . The state representation's matrices are: The indices of γ represent the backward event shift between transitions and the indices of δ represent the backward time shift. For instance in B(1, 1) = γ 0 δ 1 , the n + 1 th firing of x 1 depends on the n th firing of u 1 (no event shift); the firing date of x 1 is exactly shifted of duration a = 1 from the firing date of u 1 . A trigger of an input transition u i represents the occurrence of an event. Then, input flow events of TEG are represented by series of M ax in [[γ, δ]] as series u1 = u2 illustrated in Figure 1. In u 1 and u 2 , index 0 of γ is the first event occurrence of this input (index 1 is the second, etc.); index 1 of δ is the date of this first occurrence. The absence of an 8 th event is indicated by +∞ in monomial γ 7 δ +∞ .

A =   . . . . . . γ 0 δ 1 γ 0 δ 1 .   , B =   γ 0 δ 1 . . γ 0 δ 2 . .   , C = . . γ 0 δ 0 .

C. Time comparison of series

In the rest of the paper, time comparison of series are made thanks to a particular element called residual.

Definition 3: Let Π : D → C be an isotone 1 mapping, where D and C are complete dioids. The largest solution of Π(x) = b, if it exists, is called the residual of Π and is denoted Π . When Π is residuated, Π is the unique isotone mapping such that Π • Π I dC and Π • Π I dD where I dC and I dD are resp. the identity mappings on C and D.

The right product R a : x → x ⊗ a defined over a complete dioid D is residuated. Its residual is R a (x) = x• /a and corresponds to the pseudo-inverse of the product. Now, time comparison between series can then be defined based on residuals. It requires the use of dater functions to get times of

M ax in [[γ, δ]] series. Definition 4: Let s ∈ M ax in [[γ, δ]] be a series, its dater function is the non-decreasing function D s (n) from Z → Z such that s = n∈Z γ n δ Ds(n) .
Example 2: Series u 1 has for dater function

D u1 (0) = 1, D u1 (1) = 2, D u1 (2) = 3, D u1 (3) = 4, D u1 (4) = 5, D u1 (5) =
∀n ∈ Z, D b• /a (0) ≤ T a,b (n) ≤ -D a• /b (0), where D b• /a (0) is obtained from monomial γ 0 δ D b• /a (0) of series b• /a and D a• /b (0) is obtained from γ 0 δ D a• /b (0) of a• /b. Definition 6: Let a, b ∈ M ax in [[γ, δ]],
the time shift interval between series a and b is

∆(a, b) = [D b• /a (0); -D a• /b (0)], (1)
where

γ 0 δ D b• /a (0) ∈ b• /a and γ 0 δ D a• /b (0) ∈ a• /b.
Example 3: Generally speaking, let us consider two different series a = γ 0 δ 12 ⊕ γ 1 δ 15 ⊕ γ 2 δ 18 ⊕ γ 3 δ 21 ⊕ γ 4 δ +∞ and b = γ 0 δ 12 ⊕ γ 1 δ 15 ⊕ γ 2 δ 19 ⊕ γ 3 δ 23 ⊕ γ 4 δ +∞ . The minimal time shift between a and b is

D b• /a (0) = 0 from b• /a = γ 0 δ 0 ⊕ γ 1 δ 3 ⊕ γ 2 δ 7 ⊕ γ 3 δ 11 ⊕ γ 4 δ +∞ . The maximal time shift is -D a• /b (0) = 2 from a• /b = γ 0 δ -2 ⊕ γ 1 δ 2 ⊕ γ 2 δ 6 ⊕ γ 3 δ 9 ⊕ γ 4 δ +∞ . The time shift interval is ∆(a, b) = [0; 2]
meaning that a is equal or faster than b by a maximum of 2 time units.

III. DETECTION OF TIME SHIFT FAILURES OF (MAX,+)-LINEAR SYSTEMS WITH DISTURBANCES

This section presents an observer-based indicator that will be used for detecting and then localizing time shift failures. Like the indicator proposed in [START_REF] Paya | Observer-based detection of time shift failures in (max,+)linear systems[END_REF], it uses observers of (max,+)-linear systems as defined in [START_REF] Hardouin | Observer design for (max,+) linear systems[END_REF] (see Section III-A) and relies on the fact that time shift failures can be characterized by input disturbances (Section III-B). This new indicator is also based on a refinement of the estimated states (Section III-C) that is used to improve the accuracy of the localization information.

A. Observer of a disturbed (max,+) linear system An observer of a Timed Event Graph, as defined by [START_REF] Hardouin | Observer design for (max,+) linear systems[END_REF] and shown in Figure 3, aims at estimating the internal states of a system based on the measurement of its inputs u and outputs y o in the presence of unobservable disturbances characterized by specific inputs w. Then, the state representation of the system is:

x = Ax ⊕ Bu ⊕ Rw, y = Cx.
(

) 2
where w ∈ M ax in [[γ, δ]] l×1 with l the number of disturbed internal transitions and matrix R ∈ M ax in [[γ, δ]] n×l is filled with γ 0 δ 0 monomials where internal transitions are disturbed, with ε otherwise. The observer's equations are:

     x e = Ax e ⊕ Bu ⊕ L(y r ⊕ y o ) = (A ⊕ LC) * Bu ⊕ (A ⊕ LC) * LCA * Rw, y e = Cx e , (3)
in which x e and y e are the estimated state and the estimated output of the system and L ∈ M ax in [[γ, δ]] n×q is a new matrix called the observer matrix. Now, to obtain x e as close as possible to real state x, the observer relies on the largest matrix L ∈ M ax in [[γ, δ]] n×q such that x e x meaning:

(A ⊕ LC) * Bu ⊕ (A ⊕ LC) * LCA * Rw A * Bu ⊕ A * Rw which is given by L = (A * B• /CA * B) ∧ (A * R• /CA * R).
Matrix L actually represents the connections between the Example 4: In the TEG of Figure 2 with data of Example 1, the observer matrix is

L =   . . γ 0 δ 0 . γ 0 δ 0 . γ 0 δ 0 γ 0 δ 1 γ 0 δ 1   .
With inputs u 1 and u 2 as given in Figure 1, the estimated state

x e = [x e1 , . . . , x e3 ] T is   γ 0 δ 2 ⊕ γ 1 δ 3 ⊕ γ 2 δ 4 ⊕ γ 3 δ 5 ⊕ γ 4 δ 6 ⊕ γ 5 δ 7 ⊕ γ 6 δ 8 ⊕ γ 7 δ +∞ γ 0 δ 4 ⊕ γ 1 δ 5 ⊕ γ 2 δ 6 ⊕ γ 3 δ 7 ⊕ γ 4 δ 8 ⊕ γ 5 δ 9 ⊕ γ 6 δ 10 ⊕ γ 7 δ +∞ γ 0 δ 5 ⊕ γ 1 δ 6 ⊕ γ 2 δ 7 ⊕ γ 3 δ 8 ⊕ γ 4 δ 9 ⊕ γ 5 δ 10 ⊕ γ 6 δ 11 ⊕ γ 7 δ +∞  

B. Time shift failure as input disturbance

To take advantage of the state estimation of an observer for the detection and the localization of time shift failures, we need to characterize a time shift failure as an input disturbance. We recall here the characterization of [START_REF] Paya | Observer-based detection of time shift failures in (max,+)linear systems[END_REF]. A time shift failure is considered as a permanent phenomena and is formally defined by an unknown delay θ in a place p of a TEG. The place p, with an upstream transition x i-1 and a downstream transition x i , is characterized by a number of o tokens and a duration t (see Figure 4). Upstream transition is described by x i-1 = K n=0 γ sn δ hn (see Definition 2), where s n is the transition firing number, h n is the firing date and K the number of firing events. The downstream transition is x i = K n=0 γ sn+o δ hn+t . If a time shift failure θ > 0 holds in a place, the downstream transition then becomes: x i = K n=0 γ sn+o δ hn+t+θ . The firing dates of x i will be slowed down by θ comparing to the firing dates of x i during a nominal behavior. The same time shift failure over a place p can then be characterized by an input disturbance w i as shown in Figure 5. If this disturbance is w i = K n=0 γ sn+o δ hn+t+θ , with s n and h n defined in series x i-1 , it has the same effect on transition x i . In other words, having an unknown time shift failure θ > 0 in place p is equivalent to say that such an unknown disturbance w i exists.

Example 5: By taking back Example 4, to characterize a time shift failure in p 2 by an offset θ = 1, an input disturbance w 2 is added to the transition x 2 as illustrated in Figure 5. The transition that takes into account this offset is γ 4 δ 5+1+1 ⊕ γ 5 δ 6+2+1 ⊕ γ 6 δ 7+2+1 ⊕ γ 7 δ +∞ . Then, to copy this time shift failure effect, the disturbance has to be w 2 = γ 0 δ 4 ⊕ γ 1 δ 5 ⊕ γ 2 δ 6 ⊕ γ 3 δ 7 ⊕ γ 4 δ 8 ⊕ γ 5 δ 9 ⊕ γ 6 δ 10 ⊕ γ 7 δ +∞ .

x 2 = γ 0 δ 1+2+1 ⊕ γ 1 δ 2+2+1 ⊕ γ 2 δ 3+2+1 ⊕ γ 3 δ 4+2+1 ⊕ t p x i-1 x i o

C. Refinement of the estimated state

The indicator proposed in [START_REF] Paya | Observer-based detection of time shift failures in (max,+)linear systems[END_REF] is based on the computation of the time shift ∆(x e , x s ) (see Definition 6) where x e is the state estimated by the observer and x s results from the simulation of a fault-free model of the system based on the real observable input u. As the effect of a time shift failure might not be immediate, the first monomials of x e representing the first event occurrences might not be affected by the failure and then are considered as normal occurrences in the indicator of [START_REF] Paya | Observer-based detection of time shift failures in (max,+)linear systems[END_REF] while the failure is present. The consequence is that the indicator in [START_REF] Paya | Observer-based detection of time shift failures in (max,+)linear systems[END_REF] may generate ambiguous information and we propose here to filter out of x e these unaffected monomials and only keep as the estimated state the sensitive-to-disturbance (STD for short) part of x e . The way to extract this STD state is through a deeper analysis of the structure of the observer. It depends on values of matrix L that establishes connections between observed output y o of the system and the internal transitions x e of the observer. In the following, for a given matrix M , M ij denotes the element of M at line i, column j; M i• the i th line of M and M • i its i th column. A row L i• of matrix L represents an estimated state x ei of the system.

Proposition 1: Suppose that

L i• = ε . . . ε , then x ei = A * i• B •i u i , i
.e the estimated state x ei does not depend on disturbances.

Proof sketch. From (3), estimated state is

x ei = (A i• ⊕ L i• C • i ) * (B • i u ⊕ L i• C • i A * i• R • i w i• ). Using algebraic oper- ations of M ax in [[γ, δ]] and recalling that ε is absorbing by ⊗, it follows that x ei = A * i• B • i u i does not depend on w. Let the canonical form L ij = K k=0 γ n k δ t k be an entry of L, in the following (n 0 ) ij denotes the smallest γ index n 0 in series L ij . It is possible that (n 0 ) ij differs from 0 if tokens are already in the place represented by γ n0 δ t0 in L ij .
Let N i = min 1 j q {(n 0 ) ij } where q be the number of outputs of the system. The following proposition gives the event occurrence of the estimated state trajectory from which the effect of the disturbance can be observed.

Proposition 2: Let x ei = K κ=0 γ νκ δ τκ be the estimated trajectory of x i in its canonical form. The first event occurrence of x ei that takes into account the unknown disturbances from w i• is the event occurrence of index ν d such that

ν d = N i + ν b 0 + ν wi • 0
where ν b 0 and ν

wi • 0 are the smallest γ indices of series b = C • i A * i• R • i and w i• . Proof sketch. From (3), estimated state is x ei = (A i• ⊕ L i• C • i ) * (B • i u⊕L i• bw i• ) where b = C • i A * i• R • i
is the only part depending on disturbances w and ν d is defined as the smallest γ index of series L i• bw i• . Now, the STD state x edi is based only on γ indices from the state x ei that takes into account disturbances, which are the ones from ν d .

Definition 7: The sensitive-to-disturbance (STD) estimated state x edi is defined from x ei by

x edi = K k=d γ ν k δ τ k
where ν d is defined by Proposition 2.

Example 6: Generally speaking, if the complete estimate of a state x 1 is x e1 = γ 0 δ 4 ⊕ γ 1 δ 5 ⊕ γ 2 δ 6 ⊕ γ 3 δ 7 ⊕ . . . and if ν d = 1, the event occurrence of index ν d from which the disturbance is taken into account is the event numbered by 1. The STD estimated state is then

x ed1 = γ 1 δ 5 ⊕ γ 2 δ 6 ⊕ γ 3 δ 7 ⊕ . . .

D. Indicator of time shift failure

Based on dater and time shift functions introduced in Section 2, we propose a new time shift indicator that compares the fault-free model with the observed system through u and y o , by only using information of the STD estimated state x e and the fault-free model state x si . To perform a consistent comparison between x e and x si , event occurrences from x si before ν d must also be filtered out.

Definition 8: Let x si = A * Bu = K k=0 γ ν k δ τ k be the complete fault-free model state. Series x , si denotes the faultfree model state comparable with x edi :

x , si = K k=d γ ν k δ τ k .
Example 7: From Example 6, if x s1 = γ 0 δ 3 ⊕ γ 1 δ 4 ⊕ γ 2 δ 5 ⊕ γ 3 δ 6 ⊕ . . . , then x , s1 = γ 1 δ 4 ⊕ γ 2 δ 5 ⊕ γ 3 δ 6 ⊕ . . .. Definition 9: Let x edi be the STD estimated state of x i and x , si its fault-free model state. Indicator I xi (u, y o ) of state x i is defined as the Boolean function: Example 8: Back to Examples 1-4-5 where a time shift failure of 1 is in place p 2 . The state x , s comparable with the state x ed is the vector

I xi (u, y o ) = false if ∆(x edi , x , si ) = [0; 0], true otherwise, with ∆(x edi , x , si ) = [D x edi • /x , si (0)
[x , s1 , . . . , x , s3 ] T =   γ 0 δ 2 ⊕ γ 1 δ 3 ⊕ γ 2 δ 4 ⊕ γ 3 δ 5 ⊕ γ 4 δ 6 ⊕ γ 5 δ 7 ⊕ γ 6 δ 8 ⊕ γ 7 δ +∞ γ 0 δ 3 ⊕ γ 1 δ 4 ⊕ γ 2 δ 5 ⊕ γ 3 δ 6 ⊕ γ 4 δ 7 ⊕ γ 5 δ 8 ⊕ γ 6 δ 9 ⊕ γ 7 δ +∞ γ 0 δ 4 ⊕ γ 1 δ 5 ⊕ γ 2 δ 6 ⊕ γ 3 δ 7 ⊕ γ 4 δ 8 ⊕ γ 5 δ 9 ⊕ γ 6 δ 10 ⊕ γ 7 δ +∞   .
A failure is then detected by the indicators whose computed intervals are:

     ∆(x ed1 , x , s1 ) = [D x ed1 • /x , s1 (0) , -D x , s1 • /x ed1 (0)] = [0, 0], ∆(x ed2 , x , s2 ) = [D x ed2 • /x , s2 (0) , -D x , s2 • /x ed2 (0)] = [1, 1], ∆(x ed3 , x , s3 ) = [D x ed3 • /x , s3 (0) , -D x , s3 • /x ed3 (0)] = [1, 1].
IV. TIME SHIFT FAILURE LOCALIZATION IN (MAX,+)-LINEAR SYSTEMS WITH OBSERVER This section details how to extract information from the computed intervals of the indicators to determine the source of the failure in the TEG. All along this section, we assume that only one permanent time shift failure is present in the TEG, in a place p with a delay θ. We propose in this section to analyse how to retrieve the localization for three types of elementary structures that are contained in any TEG: tandem (see Figure 6), parallelism (see Figure 7) and synchronization (see Figure 2). Computing the localization of the failure in a specific TEG then consists in assembling these elementary analyses to get the global conclusion.

A. Tandem and parallel structures

Table I summarizes the localization results that can be expected when a tandem structure is involved given all the possible configurations. The first set of configurations (Rows 1-4) is when x 1 and x 2 can be associated with an indicator (i.e. the observer is able to provide a STD estimated state for both x 1 and x 2 ). The second set of configurations (Rows 5-7) is when only the state x 2 of the tandem has an indicator. Each row corresponds to a possible configuration of the indicator intervals and the conclusion about where the failure is in the analysed structure. For instance, in Row 3, ∆(x ed1 , x , s1 ) = [θ, θ] and ∆(x ed2 , x , s2 ) = [θ, θ] the time shift failure θ is located either in place p 1 or in a place in the upstream of the tandem.2 Indeed, x 1 = K k=0 γ n k δ t k +a and the next transition is x 2 = K k=0 γ n k δ t k +a+b . If the time shift failure θ is in place p 1 (or in a place in its upstream), x ed1 = K k=d1 γ n k δ t k +a+θ and x ed2 = K k=d2 γ n k δ t k +a+θ+b , hence the results in the indicators. Similarly, Table III summarizes the results for the parallel structure (see Figure 7). Parallelism is an association of multiple tandems which shares a common transition. The results for parallelism are a direct consequence of the ones for the tandem structure. To perform the exhaustive analysis, we suppose without loss of generality that the synchronization structure is governed by two inputs u 1 = n∈Z γ n δ Du 1 (n) and u 2 = n∈Z γ n δ Du 2 (n) where D ui (n) is a dater function (see Definition 4). Table II details all the results into 3 blocks defining the three estimation configurations: in the first one, all the x i 's are estimated; in the second one, x 3 and only one among x 1 and x 2 are estimated; and finally in the last one, only x 3 is estimated. Consider for instance Row 5. First, as ∆(x ed3 , x , s3 ) = [θ, θ], the effect of the time shift failure is always present at the synchronisation time. There are two possible situations. The first one is when ∀n ∈ Z, t 1 = D u1 (n), t 2 = D u2 (n) such that a+θ+b+t 1 > c+d+t 2 (i.e. the time shift failure is along the path p 1 , p 3 and its effect makes this path always slower than the path p 2 , p 4 ). As ∆(x ed1 , x , s1 ) = [0, 0], the time shift failure θ is necessarily in place p 3 . Indeed, x 1 = K k=0 γ n k δ Du 1 (n k )+a , x 2 = K k=0 γ n k δ Du 2 (n k )+c and the next transition is x 3 = K k=0 γ n k δ Du 1 (n k )+a+b ⊕ K k=0 γ n k δ Du 2 (n k )+c+d . If a time shift failure θ is in place p 3 the transition is x ed1 = K k=d1 γ n k δ Du 1 (n k )+a , x ed2 = K k=d2 γ n k δ Du 2 (n k )+c and the next transition is x ed3 = K k=d3 γ n k δ Du 1 (n k )+a+b+θ ⊕ K k=d3 γ n k δ Du 2 (n k )+c+d . Therefore x ed3 = K k=d3 γ n k δ Du 1 (n k )+a+b+θ . The second situation is the symmetrical one that is if a + b + t 1 < c + d + θ + t 2 , and for the same reasons, the time shift failure is in p 4 . Row 6 describes the case where the effect of the time shift failure is intermittent which leads to ∆(x ed3 , x , s3 ) = [0, θ] instead of ∆(x ed3 , x , s3 ) = [θ, θ]. In this case p 3 is suspected under the weaker condition that ∃n ∈ Z, t 1 = D u1 (n), t 2 = D u2 (n) such that a + θ + b + t 1 > c + d + t 2 . Note that here the indicator takes advantage of the underlying STD estimated state as the one of [START_REF] Paya | Observer-based detection of time shift failures in (max,+)linear systems[END_REF] would generally return [0, θ] in any of these cases.

∆(x ed1 , x , s1 ) ∆(x ed2 , x , s2 ) Localisation 2 [0,0] [0,0] ∅ 3 [θ,θ] [θ,θ] (p 1 ) 4 [0,0] [θ,θ] p 2 ∆(x ed2 , x , s2 ) Localisation 6 [0,0] ∅ 7 [θ,θ] ( p

V. CONCLUSIONS

This paper defines a method for detecting and localizing time shift failures in systems modeled as TEG using observer

Fig. 1 .

 1 Fig. 1. Representation of series u 1 = u 2

Fig. 2 .

 2 Fig. 2. A Timed Event Graph (TEG)

  6 and D u1 (6) = 7 that lists all the dates of the event occurrences. As u 1 = u 2 , u 2 has obviously the same dater function. Definition 5: Let a, b ∈ M ax in [[γ, δ]] and their respective dater functions D a and D b . The time shift function representing the time shift between a and b for each n ∈ Z is defined by T a,b (n) = D a -D b . Theorem 2 ([Max91]): Let a, b ∈ M ax in [[γ, δ]], the time shift function T a,b (n) can be bounded by:

Fig. 3 .

 3 Fig. 3. Observer structure with disturbance

FigFig. 5 .

 5 Fig. 4. Representation of a place

  ; -D x , si • /x edi (0)]. Proposition 3: The indicator returns true only if a time shift failure has occurred in the system. Proof sketch. By construction this indicator has the same detection capabilities as the indicator of [PLCPV20], hence the result.

FigFigure 2

 2 Fig. 6. Tandem Fig. 7. Parallel

  1 ) * . Internal operations and neutral elements are identical to those of B[[γ, δ]]. All series can be expressed by the following canonical form. Definition 2: The canonical form of s ∈ M ax in

STMicroelectronics, Crolles,France email:{claire.paya,philippe.vialletelle}@st.com

2 LAAS-CNRS, CNRS, Université de Toulouse, Toulouse, France {euriell.le.corronc,yannick.pencolé}@laas.fr

Π isotone ≡ ∀s, s ∈ D s s ⇒ Π(s) Π(s ).

In every Table, notation p i means that the localization is p i while the notation (p i ) means the localization is p i or a place in the upstream of p i in the TEG.

Note that transitions u 1 and u 2 are actually not part of the synchronization structure that we describe here.

that estimates their sensitive-to-disturbance states. An indicator is proposed using these STD estimated states whose interval values help to localize the failures in three types of elementary TEG structures.

A first direct perspective is to improve the localization for an other TEG elementary structure called loop, meaning any elementary structure looped between two transitions. Then, all the conclusions made separately have to be unified in a general conclusion about where the failure is localized on the TEG. Finally, we will look for exploiting more deeply the interval values of the indicator to characterize the nature of the time shift failure.


K Ogawa

A N French

T J Schmugge

Gerald Thomas

J C Ritchie

A Hsu

F Jacob

Detecting Landcover Change at Jornada, New Mexico with ASTER Emissivities

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

Introduction

Monitoring the spatial and temporal changes in land cover for semi-arid and arid land regions is required for hydrologists, ecologists, and agronomists. Land cover information is used by hydrologists to update surface conditions affecting stream flow, infiltration, and evapotranspiration [START_REF] Menenti | Energy and water flow through the soil-vegetation-atmosphere system: the fiction of measurements and the reality of models[END_REF][START_REF] Su | Remote sensing of land use and vegetation for mesoscale hydrological studies[END_REF]; by agronomists for acreage and yield prediction [START_REF] Fang | Rice crop area estimation of an administrative division in China using remote sensing data[END_REF]; and by ecologists for assessing the relationships between land degradation, human activities, and global climate change [START_REF] Chehbouni | A preliminary synthesis of major scientific results during the SALSA program[END_REF]. The need for this monitoring is especially important in dry regions, since many of these are productive agricultural lands under pressure from extreme drought and population increase [START_REF] Falkenmark | Meeting water requirements of an expanding world population[END_REF].

Land cover change has been documented for decades over critical semi-arid regions, such as West Africa [START_REF] Diouf | Monitoring land-cover changes in semi-arid regions: remote sensing data and field observations in the Ferlo, Senegal[END_REF], the Mediterranean Basin [START_REF] Sobrino | WATERMED-Water use efficiency in natural vegetation and agricultural areas by remote sensing in the Mediterranean Basin, Arid regions monitored by satellites from observing to modelling for sustainable management[END_REF], and for more than 100 years, southern New Mexico, USA [START_REF] Buffington | Vegetational changes on a semidesert grassland range from 1853 to 1963[END_REF]. In the latter instance, observations have revealed large and deleterious changes where rangeland previously dominated by grass has become dominated by mesquite and creosotebush [START_REF] Gibbens | Vegetation changes in the Jornada Basin from 1858 to 1998[END_REF], Havstad et al. 2000).

The significance of these long-term observations, though previously recognized, has only begun to be fully appreciated by the use of remotely sensed data spanning multiple decades [START_REF] Rango | Using historic data to assess effectiveness of shrub removal in southern New Mexico[END_REF]. Using frequent synoptic remote sensing observations of arid lands, especially if they were available at resolutions better than 100 m, would greatly improve our ability to monitor, analyze, and understand the implications of rapid land cover changes. objective of long term monitoring is to collect well-calibrated vegetation index (VI) data [e.g., the normalized difference vegetation index (NDVI) and its variants]. These observations are now routinely collected at daily to bi-weekly temporal sampling and at spatial resolutions between 1 m and 1 km from sensors such as Ikonos [START_REF] Dial | IKONOS satellite, imagery, and products[END_REF], Advanced Spaceborne Thermal and Reflection radiometer (ASTER) [START_REF] Yamaguchi | Overview of Advanced Space-borne Thermal Emission and Reflection Radiometer (ASTER)[END_REF], Landsat [START_REF] Goward | The Landsat 7 mission-terrestrial research and applications for the 21st century[END_REF], the Vegetation instruments aboard Satellites Pour l'Observation de la Terre (SPOT4 and SPOT5), and MOderate-resolution Imaging Spectrometer (MODIS) (Justice andTownshend 2002, Justice et al. 1998), each of which provides useful estimates of living, green vegetation. The relationship between vegetation indices and fractional cover and leaf area indices has been widely investigated (e.g., [START_REF] Gutman | The derivation of green vegetation fraction from NOAA/AVHRR data for use in numerical weather prediction models[END_REF][START_REF] Choudhury | Relationships between vegetation indices, radiation absorption, and net photosynthesis evaluated by a sensitivity analysis, Remote Sens[END_REF][START_REF] Baret | The robustness of canopy gap fraction estimates from red and near-infrared reflectances: a comparison of approaches[END_REF][START_REF] Carlson | On the relation between NDVI, fractional vegetation cover, and leaf area index[END_REF][START_REF] Jiang | Analysis of NDVI and scaled difference vegetation index retrievals of vegetation fraction[END_REF], meaning that accurate and reliable estimates of green biomass are usually obtainable from remote sensing platforms, provided ancillary data about observational conditions and expected plant cover types are also available. Based on these established relationships, time series of VI data have been used to detect both seasonal and yearly land cover changes [START_REF] Wardlow | Analysis of timeseries MODIS 250 m vegetation index data for crop classification in the u.s. central great plains[END_REF][START_REF] Telesca | Quantifying intra-annual persistent behaviour in SPOT-VEGETATION NDVI data for Mediterranean ecosystems of southern Italy, Remote Sens[END_REF][START_REF] Anyamba | Interannual variability of NDVI over Africa and its relation to El Niño/Southern Oscillation[END_REF][START_REF] Justice | Monitoring East African vegetation using AVHRR data[END_REF][START_REF] Byrne | Monitoring land-cover change by principal component analysis of multitemporal Landsat data[END_REF].

Although the vegetation index approach has demonstrated its value for land cover change monitoring, it has a frequently unmentioned shortcoming-during plant dormancy, VI values are similar to and possibly indistinguishable from bare soils. Even for indices specifically designed to minimize soil background effects, such as the Modified Soil Adjusted Vegetation Index (MSAVI) [START_REF] Qi | A modified soil adjusted vegetation index[END_REF], the spectral inputs from red and near infrared wavelengths do not readily distinguish soils from non-green vegetation. In instances where the landscape is bare soil, the VI maps are likely representative of true conditions. But for other instances where the above-ground biomass remains during winter months and is senescent or inactive, the resultant VI maps are not representative. For much of the year in arid lands, living plants are non-green; and discriminating them from soil, using indices such as NDVI, is difficult.

When assessing vegetation cover changes over a period of years, the difficulty is particularly significant because within-season VI variability is just as large, or larger than, interannual VI variability. The difficulty is exacerbated by an inability to obtain frequent remote sensing images because of cloud cover (e.g., from Landsat) or because the spatial resolution is too coarse for areas of interest (e.g., from MODIS). Unless one can be confident that VI extremes within a plant growth cycle are captured, postulated VI trends will be biased by insufficient temporal sampling. Hence, use of VI methods to infer rates of land cover change over sparsely vegetated landscapes can be misleading.

More recent work indicates VI can be improved by incorporation of hyperspectral near infrared (NIR) data, such as in the 2.2 µm region, where cellulose and lignin absorption features appear [START_REF] Bannari | Estimating and mapping crop residues cover on agricultural lands using hyperspectral and IKONOS data[END_REF][START_REF] Daughtry | Remote sensing the spatial distribution of crop residues[END_REF]. In these ongoing studies, good discrimination between soil and vegetation appear possible using end-member analyses.

When evaluating land cover over months to years, an alternative monitoring approach does exist and is based upon spatial estimation of thermal infrared emissivities. Emissivity ( λ ), a measure of thermal radiation efficiency, is defined as the ratio between actual emitted radiation (L λ ) and emitted radiation from a blackbody (L λ,BB ) at the same temperature:

λ = L λ L λ,BB (1)
Emissivity is the proportionality factor that distinguishes brightness temperatures from true radiometric temperatures [START_REF] Norman | Terminology in thermal infrared remote sensing of natural surfaces[END_REF] and is important for modeling the earth's surface energy balance. The energy balance at long wavelengths also plays a role in the ability to observe land surface emissivities:

L sensor = L surface + (1 -)L ↓ (2)
where thermal infrared radiation (TIR) observed by a sensor (L sensor ) just above the surface is the weighted sum of band surface emitted radiation (L surface ) and band downwelling sky radiation (L ↓). Band-averaged emissivity ( ) and its complement are the weighting factors, which means that estimation of surface emissivity is only possible when contrast exists between L surface and L ↓. Emissivity is independent of temperature itself and varies spectrally according to surface composition and geometry. These latter properties make emissivity observations potentially useful for land cover characterization. Measurements at laboratory scales (Salisbury andD'Aria 1992, Elvidge 1988), for example, show that emissivities of soil and vegetation are commonly distinct and do not rely upon plant chlorophyll content. This makes emissivity a potential tool for discriminating soil and vegetation emissivities. At field scales, soil and vegetation can still be distinguished (Humes et al. 1994), as illustrated by observations of senescent vegetation during a 1997 study over Oklahoma (French et al. 2000). In that instance, emissivities at wavelengths between 8 and 9.5 µm were close to 1.0 for vegetation, while those for soils were ∼0.91.

Considering the capabilities of VI and emissivity data together, a three-way characterization of land cover is possible by distinguishing bare soils, green vegetation, and non-green vegetation.

Since 2000 an important remote sensing capability became available through the launch of the ASTER instrument in December 1999 onboard the Terra satellite. ASTER uniquely offers multispectral thermal infrared (TIR) remote sensing at moderate resolution (90 m), data unavailable elsewhere. Using ASTER's five thermal bands, land surface emissivities could be estimated in ways similar to that used by French et al. (2000) for the Oklahoma study, allowing more extensive evaluation of the land cover discrimination approach.

One such evaluation arose fortuitously from ongoing studies of a semi-arid rangeland in southern New Mexico, USA, at a research site utilizing remote sensing image data to monitor landscape structural evolution and land surfaceatmosphere exchanges. The study site, known as the Jornada Experimental Range (Jornada), has been the focus of local and regional remote sensing studies since 1995. The primary theme of these studies is assessment of rangeland health, including that of vegetation land cover changes.

Vegetation cover changes at Jornada have been dramatic and occur over a wide range of time scales. Using surveys from 1858, [START_REF] Gibbens | Vegetation changes in the Jornada Basin from 1858 to 1998[END_REF] Analysis and comparisons proceed in six main steps. First, emissivity retrieval approaches from the multispectral thermal infrared sensor ASTER are considered. Three different techniques are described and implemented. This is done because consensus is lacking for optimal temperature-emissivity separation and because of concerns that inferred emissivity changes could be data processing artifacts. Second, information about the Jornada setting is discussed.

In the third step (the study focus) emissivity changes observed over Jornada are described, with particular attention paid to anomalously large changes indicative of vegetation land cover degradation. In the fourth step the significance of the changes is considered. Here tests are employed to determine the importance of interfering effects upon emissivity retrieval precision and bias. Fifth, an important related issue is addressed, namely the relationship between vegetation cover and thermal emissivities. Frequently emissivity data are estimated from VI's, an approach that effectively questions the utility of emissivity observations. Here the non-equivalence of remotely sensed VI's and emissivities are shown. Sixth, the overall significance of the ASTER-Jornada study is discussed, and ASTER emissivity patterns are interpreted.

A major benefit of multispectral thermal infrared remote sensing is its ability to estimate land surface temperatures, spectral emissivities, and broadband emissivities within the TIR window. These properties are especially valuable when observed simultaneously because they help achieve land surface temperature (LST) accuracies better than 1 • C and help distinguish between different land surfaces according to soil, rock and vegetation cover [START_REF] Tonooka | An atmospheric correction algorithm for thermal infrared multispectral data over land-a water-vapor scaling method[END_REF][START_REF] Gillespie | A temperature and emissivity separation algorithm for advanced spaceborne thermal emission and reflection radiometer (ASTER) images[END_REF], Hook et al. 1992[START_REF] Tonooka | An atmospheric correction algorithm for thermal infrared multispectral data over land-a water-vapor scaling method[END_REF]. By contrast, spectral emissivity characterization is not feasible with single-window TIR, nor with most dual-window TIR instruments. Although it is true that multiple bands also cannot resolve the under-determined TIR remote sensing problem (i.e., that for N bands, there exist N+1 unknowns, namely N emissivities plus one surface temperature), they do allow objective characterization of the emissivity spectrum, an accomplishment otherwise infeasible. This makes it possible to relate spectral variations to surface composition and surface cover, as well as help improve atmospheric corrections [START_REF] Tonooka | An atmospheric correction algorithm for thermal infrared multispectral data over land-a water-vapor scaling method[END_REF].

Various approaches to estimating emissivity from multi-band TIR data exist, examples of which are discussed in [START_REF] Jiménez-Muñoz | Improved land surface emissivities over agricultural areas using ASTER NDVI[END_REF]), Jacob et al. (2007), [START_REF] Dash | Land surface temperature and emissivity estimation from passive sensor data: Theory and practice-current trends[END_REF][START_REF] Li | Evaluation of six methods for extracting relative emissivity spectra from thermal infrared images[END_REF]. Analyses considered in this study were the Temperature-Emissivity Separation (TES, [START_REF] Gillespie | A temperature and emissivity separation algorithm for advanced spaceborne thermal emission and reflection radiometer (ASTER) images[END_REF]) approach, the Temperature-Independent Spectral Indices (TISI, [START_REF] Becker | Temperature-Independent Spectral Indices in thermal infrared bands[END_REF]) approach, and the Normalized Emissivity Method (NEM, [START_REF] Gillespie | Lithologic mapping of silicate rocks using TIMS, The TIMS data users' workshop[END_REF][START_REF] Kahle | Separation of temperature and emittance in remotely sensed radiance measurements[END_REF][START_REF] Kahle | Middle infrared multispectral aircraft scanner data: analysis for geologic applications[END_REF]) approach. For reasons discussed below, the primary approach used for the Jornada study was NEM. However, TES and TISI were also employed to demonstrate that results from the land cover change assessment were not contingent upon the choice of algorithm.

One of the better-known temperature and emissivity algorithms for ASTER images is TES, a procedure that resolves the under-determinacy 1 by employing an empirical functional relationship between emissivity spectral contrast and minimum emissivity [START_REF] Matsunaga | A temperature-emissivity separation method using an empirical relationship between the mean, the maximum, and the minimum of the thermal emissivity spectrum[END_REF]. With TES the minimum emissivity ( min ) of a sample can be closely approximated by estimating the range of either absolute or relative emissivities (∆ ) via the formula:

min = b 0 + b 1 ∆ b 2 (3)
where parameters b 0 , b 1 , and b 2 are derived from laboratory measurements of emissivities. For standard ASTER data products, the parameters are based on over 200 laboratory samples. For this study, the parameters used were 0.994, -0.687, and 0.737, respectively. By iteratively combining Eq. 3 with atmospherically corrected surface emitted radiance for each thermal band, the radiometric surface temperature can be separated from emissivities for each TIR band.

Based on simulation studies, [START_REF] Gillespie | A temperature and emissivity separation algorithm for advanced spaceborne thermal emission and reflection radiometer (ASTER) images[END_REF] show that accuracies of ±1.5 • C and ±0.015 in emissivities are possible with TES. Recent experience with ASTER and TES generally confirms this assessment, particularly for geological applications, where emissivity contrasts are large (Tonooka 2001, 1 The under-determination can also be resolved with multiple observations, as is done in the day-night method [START_REF] Watson | Two-temperature method for measuring emissivity[END_REF]. That method however is highly sensitive to measurement errors, atmospheric correction errors, and co-registration of images [START_REF] Mushkin | Extending surface temperature and emissivity retrieval to the mid-infrared (3-5 µm) using the Multispectral Thermal Imager (MTI)[END_REF].

Due to the method's complexity, readers interested in TES implementation details should consult [START_REF] Gillespie | A temperature and emissivity separation algorithm for advanced spaceborne thermal emission and reflection radiometer (ASTER) images[END_REF] and the ASTER Theoretical Basis Document for the standard data product AST05.2

For agricultural applications the TES procedure is less satisfactory because the land surfaces are dominated by high emissivities. In these cases, TES's differencing approach tends to under-estimate emissivities and consequently overestimate land surface temperatures, particularly over gray body targets such as vegetation (Jacob et al. 2004). Potential emissivity errors can exceed 2.0%, resulting in errors of 2-3 • C . This outcome is a consequence of the method's inability to distinguish between true and apparent spectral contrast,

where the latter can be due to signal noise or inaccurate atmospheric compensation. Improving temperature and emissivity retrievals for studies therefore requires alternatives to ensure that high emissivity targets remain so after the temperature-emissivity separation.

One way to achieve such a result is to utilize the normalized emissivity method (NEM, Gillespie 1985, Kahle and[START_REF] Kahle | Separation of temperature and emittance in remotely sensed radiance measurements[END_REF], an initialization procedure for TES. NEM is straightforward: a maximum spectral emissivity max appropriate for the scene (0.98 is used for this Jornada study) is specified for whichever band has the maximum surface brightness temperature. Commonly this maximum temperature appears in bands sampling wavelengths between 10 and 13 µm. Once identified, the radiometric surface temperature at that band is computed, then applied to the remaining TIR bands to extract spectral emis-sivities.

The NEM approach of pre-specifying max is not an arbitrary one and requires some judgment based on landscape context and targets of interest [START_REF] Kahle | Separation of temperature and emittance in remotely sensed radiance measurements[END_REF]. Pre-specification has a direct effect upon emissivity retrieval, where consistent selection of max allows meaningful relative comparisons of emissivities between different scenes but may cause substantial differences with respect to other temperature-emissivity separation techniques. Since max is close to 1.0, sensitivity to its pre-specification is small when accounting for downwelling sky radiation (last term in Eq. 2). [START_REF] Nerry | Bidirectional reflectivity in AVHRR channel 3: application to a region in northern Africa[END_REF] show that the impact of assuming an initial emissivity while computing a downwelling radiation correction term was small, ≤0.14%. vegetation is viewed as a canopy rather than as individual leaves and multiscattering effects would increase maximum emissivities to over 0.99 and reduce the variability of maximum emissivities [START_REF] Fuchs | Infrared thermometry of vegetation[END_REF][START_REF] Fuchs | Surface temperature measurements of bare soils[END_REF][START_REF] Sutherland | Significance of vegetation in interpreting thermal radiation from a terrestrial surface[END_REF], Norman et al. 1995[START_REF] Palluconi | The spectral emissivity of prairie and pasture grasses at Konza prairie, Kansas[END_REF][START_REF] Zhang | The spectral variation of surface emissivity within the 8-12µm "window[END_REF]. Another method that can overcome the low max problem is TISI [START_REF] Becker | Temperature-Independent Spectral Indices in thermal infrared bands[END_REF]), a technique that computes relative emissivities (T ISI i,j ) from power-scaled brightness temperatures:

T ISI i,j = T i T j n i (4)
where T i and T j are atmospherically corrected brightness temperatures for spectral channels i and j. The exponent n i can be determined from a least squares analysis:

L i ∼ α i T n i i (5)
where channel spectral radiance is L i (mW m -2 sr -1 µ m -1 ) and radiometric temperature for a channel is T i (K).

Table 1 lists the best-fit terms for Eq. 5, considering ASTER channels 10-14.

TISI values are nearly independent of land surface temperature, which means that their uncertainties can potentially be low, given accurate atmospheric corrections. However, TISI values are non-unique, dependent upon the chosen bands, and indicate emissivity ratios, as discussed by [START_REF] Becker | Temperature-Independent Spectral Indices in thermal infrared bands[END_REF].

To retrieve emissivities themselves, further observations or additional assumptions are needed. In the former case, nighttime shortwave infrared data are suitable [START_REF] Petitcolin | Land surface reflectance, emissivity and temperature from MODIS middle and thermal infrared data[END_REF], but in this study, such data at 90 m spatial scales are not available. In the latter case, a reference channel emissivity can be used in the same way as for the NEM approach .

In summary, remote sensing estimation of emissivity at Jornada can be achieved in several ways, three of which are considered here. None can be accomplished in the current context without underlying assumptions. These include assumptions about spectral emissivity contrast and the existence of an accurate emissivity reference channel.

The Jornada Experiment

The Jornada Experimental Range (Jornada) research area is a semi-arid rangeland in southern New Mexico, 30 km northeast of the city of Las Cruces and 40 km west of White Sands (Fig. 1). From 47 scene acquisitions, 27 mostly cloud-free images acquired between February 2001 and October 2003 were selected for processing (Table 2). To ensure the best possible analyses, all visible/near-infrared (VNIR) images were georegistered to within 15 m accuracies and thermal infrared (TIR) images to within 90 m accuracies using common ground control points. Estimates of spectral emissivities were done by procedures previously described in section 2, with primary emphasis upon the NEM approach.

Since the objective of this study was to investigate the ability of emissivities to detect land cover change, analysis focused on shorter TIR wavelengths where most of the emissivity variability in soils occurs (Fig. 2). At wavelengths 10.5-12 µm, emissivity variability is small for most land cover conditions, whereas for 8-10 µm interval variations potentially range from less than 0.8 to 0.98.

Not all soils exhibit this high spectral variability. Examples of low spectral variability include fine-grained, moist, and quartz-poor soils. In these circumstances, distinction between bare soils and soils covered with vegetation could be difficult. At Jornada soil emissivity variability is high, as exemplified by the dashed line in Fig. 2. Three ASTER bands sample the high-variability interval: 10, 11, and 12. Band 11 (8.475-8.825 µm) was used for this study because it is less sensitive to atmospheric correction uncertainties than Band 10 (8.125-8.475 µm) and because it did not have the anomalous gain changes of Band 12 (8.925-9.275 µm) discussed by [START_REF] Tonooka | ASTER/TIR calibration status and user-based recalibration[END_REF][START_REF] Tonooka | Vicarious calibration of ASTER thermal infrared bands[END_REF]. Furthermore, Band 11 is nearly identical to MODIS band 29 (8.4-8.7 µm), allowing spectral comparisons with this coarser resolution sensor.

To ensure the best possible TIR retrievals, ASTER thermal data were calibrated using both reference radiometric database coefficients and temporally based interpolation coefficients (needed to account for detector gain changes between routine updates). These interpolated values were particularly important for calibration of ASTER band 12.

To estimate vegetative cover, both remote sensing and ground measurements To help verify ASTER emissivities at Jornada, ground collections of multiband TIR radiometric data from four sites were used (Table 3). In situ samples were acquired with a tripod-mounted 5 band CIMEL-312T instrument, while loosely-bagged soil samples were measured (hemispheric reflectance) in the laboratory using a Beckman spectrophotometer at the Jet Propulsion Laboratory, Pasadena, California (http://speclib.jpl.nasa.gov). Four of the resulting emissivity spectra were collected at the transect sites (1 at the Grass site, 1 at the Transition site, and 2 at the Mesquite site) and corresponded to quartz-bearing soils with low emissivities between 8 and 12 µm. The fifth spectrum corresponded to gypsum from the western edge of White Sands (32.8222 • N, 106.4125 • W). In contrast to quartz-rich soils, gypsum has generally high emissivities (>0.96) except for 8.5-8.9 µm wavelengths, where it Although accurate retrieval of land surface emissivities is important, a greater concern for land cover change assessment is consistent emissivity retrieval, since it is the temporal change of emissivity that will determine the effectiveness of the analysis. When considering interfering factors for change assessment, four are critical: instrumental noise, georegistration errors, anomalously wet soils, and inaccurate atmospheric corrections. Fortunately, most of these errors for the Jornada study were small and unimportant. In particular, instrumental errors in the TIR bands were <0.4 • C [START_REF] Tonooka | Vicarious calibration of ASTER thermal infrared bands[END_REF]. Georegistration errors were also small, <90 m, because all scenes were registered with ground control points. Considering the variability scale of vegetation at Jornada ranges between 3 and 8 m [START_REF] Pelgrum | Length-scale analysis of surface albedo, temperature, and normalized difference vegetation index in desert grassland[END_REF], in conjunction with a 90 m pixel size, apparent temporal changes in emissivities due to misregistrations were not significant. Wet soils too were rarely a problem because rainfall shortly before selected overpasses did not occur. The one exception occurred in September 2002, when rainfall was heavy shortly before overpass time. The effect of this event is discussed further in Section 5.

Since three of four critical interfering factors were not significant concerns for the Jornada study, the greatest concern for change analyses was ensuring adequate removal of atmospheric effects. Inaccurate corrections could overwhelm small emissivity changes and invalidate results. Thus, quantifying the potential effects of atmospheric correction uncertainty became important, especially when using atmospheric profiles that were not acquired at overpass times, nor co-located over Jornada.

To assess correction uncertainties, a sensitivity analysis of atmospheric profile correction was conducted using simulations. By assuming temporal variability equivalence between the profile over the radiosonde site and profiles above Jornada, and by developing relationships between the correction terms and columnar water vapor (the dominant factor affecting TIR corrections), the expected effect of uncertain TIR corrections upon emissivity could be estimated.

The In step one, columnar water vapor for 2159 EPZ radiosonde profiles was computed and plotted by time (Fig. 3). Water vapor amounts ranged ∼0.2-3.4 cm with strong seasonal increases associated with summer monsoons. Day-today water vapor also varied seasonally, ranging ∼25-45% of columnar water amounts. To include a wide range of realistic simulated atmospheric conditions throughout the year, the upper limit value (45%) was used for subsequent simulations. In step two, MODTRAN runs were performed using a subset (359) of the 2001-2003 atmospheric profile database. Outputs were filtered with ASTER TIR spectral response functions, modeled with power functions, and checked against the remainder of the database. The estimation functions were:

τ i = b 0i + exp [b 1i w p i ] (6)
for atmospheric transmissivity, and:

L i = b 0i + b 1i w p i (7)
for path radiance (upwelling and downwelling). The terms b 0,i , b 1,i , and p i correspond to the parameters in Tab. 4. w is columnar water vapor (cm). Eqs. Atmospheric corrections were done according to the atmospheric radiative transfer equation:

L surface,i = L sensor,i -L ↑ i τ i -(1 -i )L ↓ i (8)
where L surface,i is radiance emitted by the surface for channel i and L sensor,i is radiance observed by the sensor. For each scene, the closest-in-time radiosonde profile was used to derive the mean value of each of the atmospheric correction terms in Eq. 8. Apparent instead of observed columnar water vapor values for each term were then computed using the inverse of Eqs. 6 and 7. This approach accounts for differences between predicted and observed atmospheric transmissivities and path radiances. Using the columnar water variability estimated in step 1 as the standard error of mean water vapor, random values of apparent water vapor were generated to allow forward solutions of Eqs. 6.

Having obtained simulated the atmospheric correction terms τ , L ↑ , and L ↓ , surface emissivities could be estimated from the NEM algorithm. To estimate emissivity uncertainties caused by random errors in atmospheric corrections, linear trend lines were fit to data in each site. De-trending data for the 'Reference' site was statistically insignificant (R 2 =0.008, p<2.2E-16).

Residual error for the 'Degraded' site was 0.014, and for the 'Reference' site was 0.016, meaning that random atmospheric perturbations to NEM-derived Band 11 emissivities were ≤1.6%.

Observed Temporal Changes

Temporal change in land cover between 2001 and 2003 over Jornada was assessed in two ways: one in terms of change in land surface emissivities, and another in terms of change in NDVI. Although seasonal variations were likely embedded within the data series, only 3-year linear trends were considered be-variations at higher orders would have placed greater demands on accurate temporal sampling and would not allow suppression of erroneous short term (weekly to monthly) changes.

ASTER Band 11 Emissivity Change

Considering the georegistered stack of 27 ASTER images, linear fits were computed for each pixel location by regressing band 11 emissivities against cumulative days since 1 January 2001. The resulting emissivity slope values (Fig. 6) ranged between -3% year -1 to +1% year -1 . For most of the Jornada region band 11 emissivity trends were small (<1%), suggesting no significant linear change in emissivity patterns occurred for this three year period. However, the west and along rangeland boundaries to the north. These terminations could not have been caused by instrumental, processing, or atmospheric correction errors, meaning that the observed decreased emissivities were not data artifacts but were indicative of physical land cover changes corresponding to rangeland management practices.

These emissivity changes are consistent with retrieval algorithms other than NEM. To show this, the previously discussed 'Degraded' and 'Reference' sites were selected. Comparison of trends (Fig. 8) obtained from three different emissivity retrieval algorithms-NEM, TES, and TISI-showed similar linear trends for both sites. Shown in black are emissivities for the 'Degraded' site, and shown in gray are emissivities for the 'Reference' site. Absolute emissivity values from NEM, TES, and TISI significantly differed, despite use of identical remote sensing data and initial emissivities. However, the trend line slopes were similar, -2.09 to -2.70% year -1 , with R 2 values ranging 0.60 to 0.66 (Table 5). Slopes for the NEM and TISI emissivity trends were statistically identical and only slightly different from the TES trend. Standard errors for methods were also similar (∼0.04 %/year -1 ).

The significance of the coherent regions can also be tested by comparison against regions with no land cover change. One such region is White Sands, a ∼20 km x 25 km area covered by generally dry, unvegetated gypsum sand.

Gypsum is an especially good stability test for ASTER band 11 due to its large potential sensitivity, ranging between 0.72 for a dry surface and 0.98 for a wet surface. Selecting a 12 x 20 pixel patch well within the White Sands area (1080 m x 1800 m, 32.8348 N, 106.2925 W)3 showed no significant emissivity change occurred over the 2001-2003 period (Fig. 9). Linear trend analysis for 9

scenes returned a slope of 0.3%/year -1 , with R 2 =0.02 and p value of 9.5E -15 .

MODIS Band 29 Emissivity Change

The spatially coherent decreasing emissivity zones identified by ASTER are sufficiently large to compare with 1-km Terra MODIS observations. If sensors such as MODIS corroborate ASTER-based results, there will be not only the possibility of greatly expanded regional studies at 1 km resolution but also much more confidence that the changes are real and not instrumentation artifacts. MODIS has only three TIR bands comparable to ASTER's five (Fig.

2

) but can still be analyzed with the NEM, TES or TISI approaches. For compatibility, the NEM approach was used for MODIS. The spectral band equivalent to ASTER band 11 is MODIS band 29, with the MODIS response weighted for slightly shorter wavelengths relative to ASTER band 11.

MODIS band 29 emissivity changes were estimated using the same 27 overpass days identified in Table 2. The resulting patterns (Fig. 10) are remarkable for their similarity to changes observed by ASTER (Fig. 6). Despite the coarser spatial resolution, band 29 emissivities from MODIS show the same 3% year -1 decrease within the 'Degraded' site (larger circle). Just as for the To the extent that NDVI observations were correlative to vegetation cover, these results can be qualitatively checked by comparison against LAI transects collected at the Grass, Transition, and Mesquite sites. The transects were important independent observations of vegetation cover, but unfortunately because of their locations and sampling frequencies, they could not be used in a rigorous way to confirm vegetation conditions in the 'Degraded' region. Explanation for the relatively poor discrimination ability of NDVI data may be provided by comparison against LAI transect data (Fig. 13). The relationship between LAI and NDVI agrees with the previous time series (Fig. 11) that showed trends with little distinction between sites. Statistical results (Table 7) show moderate to good explained variance (0.50-0.81).

LAI/NDVI slope relationships were similar and agreed within standard errors.

The relationship between LAI and emissivity, on the other hand, suggests a stronger correlation exists for the sparse canopies at Jornada. Explained variance was good to excellent (0.81-0.99) for the three sites with significant dis-crimination between vegetation-emissivity slope relationships at the Mesquite 6 Emissivity-NDVI relationship

An issue not discussed so far is the relationship between emissivity and vegetation indices such as NDVI over arid landscapes. Some suggest a meaningful relationship between them exists (e.g., [START_REF] Van De Griend | On the relationship between thermal emissivity and the normalized difference vegetation index for natural surfaces[END_REF][START_REF] Valor | Mapping land surface emissivity from NDVI: application to European, African, and South American areas, Remote Sens[END_REF][START_REF] Bolle | Meditteranean Land-Surface Processes Assessed from Space[END_REF]), meaning that difficult-to-obtain emissivity data could be replaced with much more available NDVI data. This substitution would be especially valuable because VNIR data also typically have greater spatial resolutions than possible with TIR data. The basis for the relationship is that surface emissivities are low over bare soils and high over vegetation canopies. Soil, when dry, often exhibits low emissivities (0.9 or less), while vegetation canopies are efficient scatterers and effectively blackbodies. Hence, as long as vegetation canopies are also photosynthetically active, vegetation indices should be able to track emissivity changes indirectly. [START_REF] Van De Griend | On the relationship between thermal emissivity and the normalized difference vegetation index for natural surfaces[END_REF], for example, proposed a logarithmic formula relating broadband emissivity to ground-based NDVI observations based on calibrations over sites in Botswana:

= 1.0094 + 0.047 ln (NDVI) (9) [START_REF] Van De Griend | On the relationship between thermal emissivity and the normalized difference vegetation index for natural surfaces[END_REF] suggest scale independence for the -NDVI relationship and demonstrate its use at satellite remote sensing scales. For longer wavelengths in the TIR window, 10.5-12.5 µm, [START_REF] Valor | Mapping land surface emissivity from NDVI: application to European, African, and South American areas, Remote Sens[END_REF] show that by considering vegetation structure, modification of Eq. 9 allows its applicability to regions outside of Botswana. The question remains, however, how well models such as Eq. 9 could predict emissivities for land covered by sparse, dormant vegetation.

Using to NDVI, and there is a non-unique relationship for NDVI values between 0.1 and 0.2. For NDVI values less than 0.2, the slope of the scatter plots is essentially vertical, meaning there is no meaningful relationship between NDVI and spectral emissivities. For most dates displayed in Fig. 14, the maximum emissivity was nearly reached by NDVI ∼0.2. The black curve (Eq. 9), even if translated vertically, did not accurately represent observations, with emissivity discrepancies sometimes exceeding 3%.

Performing these analyses for broadband emissivities, instead of narrow band ones, likely would not substantively change these results [START_REF] Gieske | Temperature-emissivity separation with ASTER and LANDSAT 7 imagery validation on the fring of the Okavango Delta, Botswana[END_REF]).

Using a wider spectral interval would reduce the large emissivity variability seen at low NDVI while maintaining the flat response for higher NDVI. This lower variability may explain some of the success reported by [START_REF] Momeni | Evaluating NDVI-based emissivities of MODIS bands 31 and 32 using emissivities derived by Day/Night LST algorithm[END_REF], who evaluated emissivities at 10.5-11.5 µm wavelengths at 1 km spatial resolution. In contrast, Jornada results show no reliable relationship between emissivity and NDVI exists at 90-m scales. While there could be ways to infer emissivities, using NDVI improved with geometrical parameters [START_REF] Jiménez-Muñoz | Improved land surface emissivities over agricultural areas using ASTER NDVI[END_REF], the likelihood of its success in Jornada-like environments is poor because emissivity variability is not well-correlated to plant greenness.

Discussion

Retrieval and analysis of ASTER TIR data over southern New Mexico range- We interpret these patches as areas of decreased vegetation densities where sparsely vegetated land were further degraded, exposing greater amounts of bare soil over a three year period. What is notable for the emissivity data is not the short term, seasonal changes in vegetation-these might be observed in a better way using VI data -but rather the inter-annual changes regardless of seasonality. The physical reason for this difference is explained by the dynamics of rangeland vegetation and the fact that VI and emissivity observations respond to different biochemical and structural properties of vegetation.

VI responds to plant chlorophyll densities, while emissivity responds to plant canopy geometry and patterns between plant canopies. This difference in response was described in Section 6. Since VI changes are seasonal according to plant growth cycles, their variability within a year is just as large as between years, which means that it can be difficult to distinguish long-term change from short-term changes. This difficulty has been noted by others (e.g., [START_REF] Li | Study on land cover change detection method based on NDVI time series datasets: change detection indexes design[END_REF].

Emissivity changes, on the other hand, are less likely to show such seasonal variability over persistent vegetation because the multi-scattering effects result from plant matter regardless of chlorophyll content. At time scales greater than a year, plant distributional patterns do change, and in these cases corresponding emissivities changes can be detected. Hence, emissivity variability within a year (absent surface wetting events) is small, while variability be- Alternative interpretations of the Jornada emissivity patterns were also considered, such as the possibility of data collection and processing artifacts. These included TIR calibration errors, detector degradation, and inaccurate atmo-TIR calibration and detector problems were unlikely because the ASTER gains were episodically monitored and updated to reflect correct gains and offsets.

Poor atmospheric corrections, though certainly a concern, were demonstrated to potentially cause emissivity errors no greater than 1.6%. Furthermore, the correction errors from radiosonde profiles would not cause systematic trends, nor would such errors induce the discrete and the well-delineated patches evident in Fig. 6.

Another interpretation is that the emissivity changes were due to changes in soil emissivity. These could be induced by various factors such as change in surface organic matter, surface texture, and surface moisture. Soil emissivities in particular have some dependency upon surface grain size (e.g., [START_REF] Salisbury | The effect of particle size and porosity on spectral contrast in the mid-infrared[END_REF]. It is plausible-though in this instance unverifiable-that the observed emissivity changes were related to grazing activity which may have disturbed the surface soil texture. Soil emissivity changes can also be seen with changing surface moisture, but in these cases the changes occur over short time spans, such as observed for the September 17, 2002 overpass.

Considering that these alternatives were either unlikely or conjectural, the more probable explanation for the three-year emissivity patterns, whether decreasing or increasing, is due to change in vegetation canopy densities. For sparsely vegetated terrain, the interaction between high-emissivity vegetation and low-emissivity soils explains in a physically meaningful way the spatial and temporal distribution of observed emissivity patches. As vegetation densities increase, low emissivity surfaces are increasingly masked by higher emissivity surfaces. There is also evidence from ground transect data that the temporal emissivity decrease is associated with LAI decreases at three peripheral sites (Fig. 12 and Table 7). With two exceptions, canopy LAI values at the Grass, Transition, and Mesquite sites dropped from ∼1.4 to ∼0.6, a change of more than 50%. The environmental causes for the vegetation density changes are unknown and apparently are not linked to short term precipitation patterns (Fig. 12).

Conclusions

Detection of land cover change is an important role for remote sensing whereby long-term observations can be used to monitor spatial and temporal seasonal to multi-year patterns in vegetation cover. The usual approach to this task is to construct VI maps, derived from normalized red and near infrared data.

Generally these maps have good quality and are extremely valuable for tracking seasonal vegetation changes. Nevertheless, VI techniques distinguish nongreen vegetation from background soils with difficulty, which means that signals from long term changes can be overwhelmed by seasonal effects. Such circumstances commonly arise for both cultivated and uncultivated lands.

In this study we have shown a complementary approach that may improve abilities to discriminate vegetation cover regardless of plant color. Using thermal infrared observations from a one-of-a-kind remote sensing instrument, ASTER, seasonal land cover changes can be detected by using spectral emis- These results highlight the importance of multispectral thermal infrared data that includes observations at wavelengths within 8-9.5 µm. The value of TIR data extends beyond land surface temperature retrieval to include the information rich portion of the TIR window. Commonly, TIR detectors sample TIR data at wavelengths between 10-13.5 µm for split window analyses. For land surface applications, variability of surface emissivities at these wavelengths is typically small and difficult to use for change assessment in the manner described in this study. A TIR remote sensing strategy that is better for both temperature and emissivity estimation is to include observations within the 8-9.5 µm interval, where emissivity variations due to soils and vegetation are frequently large. These data, in combination with longer wavelength TIR data, will provide valuable land cover information unavailable using other remote sensing bands. Future work will expand the emissivity change assessment to more sites throughout the U.S. Southwest.

  [START_REF] Dash | Separating surface emissivity and temperature using two-channel spectral indices and emissivity composites and comparison with a vegetation fraction method[END_REF] found yet smaller uncertainties. For combined soil and vegetated terrain, band averaged maximum emissivities are commonly within 1.5% of 0.98. Considering soil emissivity spectra from the Johns Hopkins collection within the ASTER spectral library(Salisbury and D'Aria 1992, Jet Propulsion Laboratory 2001), 39 of 41 samples (95%) had max >=0.965 somewhere in the TIR window. For this library set, 38 of the 41 had the maximum emissivity located in either ASTER band 13 or 14 (10.25-11.65 µm), while 3 were associated with the shorter wavelength. No maximum values were found in ASTER bands 11 or 12(respectively, 8.475-8.825 and 8.925- 9.275 µm). At more typical remote sensing scales (e.g., 90 m for ASTER),

  were used. For the remote sensing data, NDVI values were derived from ASTER bands 2 and 3N reflectances previously aggregated to 90 m. The reflectances were inferred from at-sensor radiances by using the 6S atmospheric radiative transfer model[START_REF] Vermote | Second simulation of the satellite signal in the solar spectrum, 6S: and overview[END_REF]) and NOAA radiosonde profile data (http://raob.fsl.noaa.gov) from Santa Teresa, NM (EPZ, 31.90 • N, 106.70 • W, ∼ 75 km south of Jornada). For the ground measurements, the three sites previously mentioned (Grass, Transition, Mesquite) were used for semi-annual leaf area index (LAI) transect data. Each transect was 150 m long with LICOR LAI-2000 measurements collected at 1-meter intervals for three 30-m sections through vegetation including grass and shrubs. The measurements approximate total canopy cover, whether green or not.

  tests were performed in four steps. The first step quantified variability of columnar water vapor over the Santa Teresa radiosonde site (EPZ) by considering the atmospheric database from 2001 to 2003. The second step established empirical relationships between columnar water vapor and atmospheric correction terms for each of the five ASTER TIR bands. This was done by performing simulations using MODTRAN radiative transfer code and a subset of the EPZ database. The third step generated apparent emissivities using the relationships from step 2 and the variability data obtained from step 1. Specifically, MODTRAN simulations were performed and applied to the NEM algorithm for each of the 27 ASTER scenes using radiosonde profiles and normally distributed perturbations to the correction terms. The fourth step extracted the statistical outcomes for ASTER band 11 emissivities from each site of interest.

6 and 7

 7 are valid for columnar water vapor amounts between 0.25 cm and 2.5 cm. Graphical examples of these functions compared with MODTRAN output values are shown in Fig. 4 for ASTER band 11. In step three, 1001 MODTRAN simulations were run using randomly gener-ated, Gaussian distributed, τ (transmissivity), L ↑ (upwelling path radiance), and L ↓ (hemispherically integrated downwelling path radiance) correction terms for each of the 27 ASTER scenes. The chosen number of simulations was based on a compromise between the need for large sample sizes and time needed for computations. Distribution plots from test runs indicated 101 simulations were insufficient, whereas output from 10001 runs were only marginally different from 1001 runs (odd numbers were used to simplify median sampling).

  Lastly, in step four ASTER band 11 estimates for each pixel (1.05E9 of these) were consolidated for statistical analysis. Aggregation by target site resulted in data such as shown in Fig.5, where ASTER band 11 emissivity quantiles are plotted for two sites spanning 2001-2003. Fig.5illustrates outcomes for two sites denoted as 'Degraded' and 'Reference'. These sites did not have ground observations and were selected for comparative purposes. The 'Degraded' site is the location of apparent land cover degradation and will be further discussed below. The 'Reference' site was chosen because it had no significant land cover change, as seen from remote sensing data. The site coincides with a fenced pasture having restricted access. In Fig5a, the temporal progression of ASTER band 11 emissivities are displayed for the 'Degraded' site. Fig.5bdisplays the temporal progression for the nearby 'Reference' site, where emissivities showed less change. The symbols are conventional box and whisker displays of quantiles (1st quartile, median, 3rd quartile), except for the whiskers, which represent the full emissivity range. All 27 selected overpass times could be used for the 'Reference' site, but two fewer were used for the 'Degraded' site due to limited ASTER coverage for overpasses on 8 February and 7 November 2002.

Fig. 6

 6 Fig. 6 also revealed spatially coherent regions-one of them ∼5 km x 10 kmnear the southern and western boundaries of the USDA Jornada Range, where emissivities decreased by 2-3% year -1 . These decreases were greater than the 1.6% residual uncertainty estimated from the atmospheric simulation tests. A portion of this emissivity decrease region is delimited by the large black circle (repeated in Figs. 7 and 10). By comparison, a trapezoidal-shaped area (surrounding the small black circle indicating the 'Reference' site) showed negligible emissivity change. The explained variations of these emissivity changes are indicated in Fig. 7, which shows that the estimated linear trends have R 2 values up to 0.85 within the larger circled area. The trends were highly significant, with p values <2E10 -16 . The extents of both the trend and the R 2 values, furthermore, were strongly correlated to land use patterns. Note for example the abrupt spatial termination of decreased emissivities along a roadway to

Fig. 12

 12 Fig. 12 summarizes the ground observations for 2001-2003. The three monthly rainfall bars bars at the base of the plot show that cumulative rainfall did not correspond to monsoon periods.

  Plotted in the top figure are LAI vs. ASTER band 11 emissivities, and plotted in the bottom figure are LAI vs. NDVI. For conditions ranging between bare soil and full cover, an linear-exponential relationship between LAI, emissivity, and NDVI could be expected, where asymptotes were reached for LAI values greater than ∼3.0. For Jornada land cover conditions linear fits (indicated as solid lines) were sufficient because LAI values were less than 2.0. Six LAI transect observations were made, but one set made on 17 September 2001 was excluded (open symbols) because of very wet surface conditions. Such conditions affected all sites by increasing emissivities.

  site and the Transition and Grass sites(58.02 LAI/emissivity vs. 22.41 and 21.38). This discrimination is possibly related to spatial clumping of vegetation. At the Mesquite site, shrubs are strongly clumped with large bare soil expanses between dunes but at the Transition and Grass sites the vegetation is more uniformly distributed. A consequence of this clumping is that increases of LAI at the Mesquite site obscure significantly less soil than do increases at the other two sites. NDVI appears insensitive to clumping because the bulk of the canopies are non-green.

  ASTER TIR and VNIR data from 15 different overpasses during 2001-2003 at Jornada (Fig. 14), the question can be partly answered. For this semiarid landscape, dependence of ASTER Band 11 (8.475-8.825 µm) emissivity upon NDVI is in most instances questionable. With possible exceptions for scenes acquired after 15 May 2002, Band 11 emissivities are poorly correlated

  land shows strong, consistent, and coherent regions of decreasing band 11 emissivities over the three-year period, indicating that land cover change for this environment can be monitored from space. ASTER thermal infrared observations over the Jornada Experimental Range and the adjacent New Mexico State University Ranch have revealed patches of land, ∼5 km x 10 km, where emissivities at 8.6 µm decreased on the order of 3% year -1 , with early 2001 NEM-based emissivities of ∼0.93 decreasing to ∼0.87 by late 2003.

  tween years could be large. To illustrate how the land cover change could be estimated, the regression results from comparing LAI transect data to band 11 emissivities (Tab. 7) were applied to the entire collection of 27 ASTER scenes and plotted as dashed gray lines in Fig 12. The tentative trends (i.e., regression results were based on too few observations to warrant greater certainty) show how LAI values at the Mesquite site may have decreased from ∼2.0 to 1.3 between 2001 and 2003. For the Transition and Grass sites, the decrease appears to be from 1.5 to 0.8.

  Figures_1

Table 1

 1 Power function parameters for ASTER TIR bands according to: L = αT n , with spectral radiance L (mW m -2 sr -1 µ m -1 ) and temperature (K). λ is the band's

	central wavelength (µm).
	Band	λ	α	n
	10	8.2819 1.7825739e-10 5.5393314
	11	8.6313 6.2117374e-10 5.3234390
	12	9.0757	2.6602331e-9 5.0742792
	13	10.650	1.5377459e-7 4.3605276
	14	11.2812 5.4625771e-7 4.1325621

Table 3

 3 Field, laboratory, and ASTER band 11 emissivities at the Grass, Transition, Mesquite, and White Sands sites. Two soil sample types were acquired at Mesquite.

	Data sources were in situ Cimel-312T radiometric observations and bagged surface
	samples for laboratory spectrophotometry.
	Site	Source B11 Emissivity
	Grass	Cimel	0.893
		Lab	0.882
		ASTER	0.886
	Transition	Cimel	0.800
		Lab	0.803
		ASTER	0.872
	Mesquite	Cimel	0.718
	(Bright Sand)	Lab	0.703
		ASTER	-
	Mesquite	Cimel	0.847
	(Dark Sand)	Lab	0.880
		ASTER	0.853
	White Sands	Cimel	0.740
	(Gypsum)	Lab	0.724
		ASTER	0.725
	is strongly reflective ( ∼0.72-0.74). Agreement between field, laboratory, and
	ASTER data are generally very good. Emissivity differences were less than

0.016 in most instances. Exceptions occurred where samples could not be collected intact (surface crusts at Mesquite) and where ASTER could not resolve spatial heterogeneities (at Transition).

Table 5

 5 Emissivity trend estimation over the 'Degraded' site with NEM, TES, and TISI.

	Method Slope (%/year) Standard Error (%/year) R 2
	NEM	-2.09	0.036	0.63
	TES	-2.70	0.043	0.66
	TISI	-2.09	0.039	0.60

Table 6

 6 NDVI change statistics for five Jornada sites. Slope and standard error of slope shown as NDVI % year -1 .

	Site	Slope Standard Error Residual (%) R 2
	Degraded	2.3	0.2	6.1	0.06
	Reference	2.3	0.2	5.8	0.07
	Grass	2.5	0.2	6.0	0.08
	Transition	2.3	0.2	6.0	0.07
	Mesquite	2.6	0.2	6.1	0.08

Available at: http://eospso.gsfc.nasa.gov/eos homepage/for scientists/atbd/docs/-ASTER/atbd-ast-05-08.pdf

ASTER resolution did not allow selecting the same site (Tab. 3) as used previously for ground samples because of mixed-pixel problems.

This work was supported in part by NASA EOS Grant 03-OES-02 and by the USDA/ARS Jornada Experimental Station. The laboratory soil sample measurements were made by Cindy Grove, Jet Propulsion Laboratory, California Institute of Technology, Pasadena, California. Data analysis and presentation were greatly facilitated by using the R package (R Development Core Team 2006). Precipitation data were provided by the Jornada Long-Term Ecological Research (LTER) project and were funded by the U.S. National Science Foundation (Grant DEB-92-40261).

site (smaller circle) showed no significant change over the 2001-2003 period. Five km eastward of the main region with decreasing emissivities lies a smaller region (with light gray tones) where emissivities increased by 1%. The same patch of increasing emissivity was observed by ASTER (Fig. 6).

ASTER NDVI Change

Viewing Jornada from an NDVI perspective showed similar but less significant temporal changes than seen with emissivity (Fig. 11). Shown are trends for five sites: the 'Degraded' site, the 'Reference' site, and the three transect sites: Grass, Transition, and Mesquite. Plotted are NDVI quantiles for each site by elapsed days since 1 January 2001. The trendlines are least squares fits to the NDVI data over each site. The gray boxes approximate the extent of the monsoonal season, the most likely period for rainfall.

Visual inspection of trends suggested that NDVI decreased over time in the same way as it did for emissivity, including a weak suggestion of seasonally periodic changes. However, NDVI values at Jornada were usually low-ranging between 0.0 and 0.2-and there was no apparent distinction in vegetation changes between the 'Degraded' and 'Reference' sites. Statistical analyses confirmed this observation. Unlike emissivity trends (Tab. 5), the NDVI linear trends were not significant (Tab. 6), with p values < 2.2E-16. Considering standard errors, NDVI slopes for all five sites were indistinguishable. R 2 values were all negligible, meaning that NDVI changes over the three years could not be explained by linear trends. Jacob, F., Schmugge, T., Olioso, A., Courault, D., French, A., Ogawa, K., Petitcolin, F., Chehbouni, G., Pinheiro, A. and Privette, J.: 2007