J Frédéric Bonnans
email: frederic.bonnans@inria.fr.

Pascal Jaisson
email: pascal.jaisson@math.uvsq.fr.

Optimal control of a parabolic equation with time-dependent state constraints

Keywords: Optimal control, parabolic equations, state constraints, alternative optimality system, quadratic growth, second-order optimality conditions, sensitivity, expansion of value function and solutions

In this paper we study the optimal control problem of the heat equation by a distributed control over a subset of the domain, in the presence of a state constraint. The latter is integral over the space and has to be satisfied at each time. Using for the first time the technique of alternative optimality systems in the context of optimal control of partial differential equations, we show that both the control and multiplier are continuous in time. Under some natural geometric hypotheses, we can prove that extended polyhedricity holds, allowing to obtain no-gap second-order optimality conditions, that characterize quadratic growth. An expansion of the value function and of approximate solutions can be computed for a directional perturbation of the r.h.s. of the state equation.

Commande optimale d'une équation parabolique avec contraintes sur l'état dépendant du temps

Résumé : Cet article étudie le problème de la commande optimale de l'équation de la chaleur par un contrôle distribué sur une partie du domaine, en présence d'une contrainte sur l'état. Cette dernière est intégrale en espace et doit être respectée à chaque instant. Utilisant pour la première fois la technique des multiplicateurs alternatifs dans le contexte de la commande optimale des équations aux dérivées partielles, nous montrons la continuité en temps de la commande et du multiplicateur. Sous des hypothèses naturelles, la condition de polyédricité étendue est satisfaite, ce qui permet d'obtenir des conditions d'optimalité du second ordre sans écart, qui caractérisent la croissance quadratique. On en déduit également un développement de la valeur et des solutions approchées pour une perturbation directionnelle du second membre de l'équation d'état.

Mots-clés : Commande optimale, équations paraboliques, contraintes sur l'état, système d'optimalité alternatif, croissance quadratique, conditions d'optimalité du second ordre, sensibilité, développement de la valeur et de la solution.

Introduction

Optimal control problems of PDE (partial differential equations) with state constraints have been intensively studied since the eighties, starting with the work by Bonnans and Casas [START_REF] Bonnans | On the choice of the function spaces for some state-constrained control problems[END_REF][START_REF] Bonnans | Contrôle de systèmes elliptiques semilinéaires comportant des contraintes sur l'état[END_REF][START_REF] Bonnans | Optimal control of semilinear multistate systems with state constraints[END_REF] and then Abergel and Temam [START_REF] Abergel | Optimality conditions for some nonqualified problems of distributed control[END_REF]. More specifically, problems with mixed control and state constraints, or with "pure" state constraints, for parabolic systems have been studied in Casas [START_REF] Casas | Pontryagin's principle for state-constrained boundary control problems of semilinear parabolic equations[END_REF], where a regularity result for the costate equation with r.h.s. measure was obtained, and Casas, Raymond and Zidani [START_REF] Raymond | Pontryagin's principle for state-constrained control problems governed by parabolic equations with unbounded controls[END_REF][START_REF] Casas | Pontryagin's principle for local solutions of control problems with mixed control-state constraints[END_REF]. The Stone-Čech compactification theorem has been introduced for deriving optimality conditionsfor parabolic problems with pure or mixed state constraints by Arada and Raymond [START_REF] Arada | Necessary optimality conditions for control problems and the Stone-Čech compactification[END_REF][START_REF] Arada | Minimax control of parabolic systems with state constraints[END_REF][START_REF] Arada | Optimal control problems with mixed control-state constraints[END_REF][START_REF] Arada | Dirichlet boundary control of semilinear parabolic equations. II. Problems with pointwise state constraints[END_REF]. We note the recent study by Casas, De los Reyes and Tröltzsch [START_REF] Casas | Sufficient second-order optimality conditions for semilinear control problems with pointwise state constraints[END_REF] on sufficient conditions.

In this paper the system is described by a nonlinear parabolic equation whose properties are recalled in section 2. We consider the case of a standard quadratic control and of a finite number of state constraints for each time; our model example is the one of a unique state constraint which is the L 2 norm (integration w.r.t. space only). This allows us in section 3 to fully characterize qualification and to use the alternative formulation of the optimality conditions, introduced in Jacobson, Lele and Speyer [START_REF] Jacobson | New necessary conditions of optimality for control problems with state-variable inequality contraints[END_REF], rigorously obtained in Maurer [START_REF] Maurer | On the minimum principle for optimal control problems with state constraints[END_REF] and generalized in Bonnans and Hermant [START_REF] Bonnans | Second-order analysis for optimal control problems with pure state constraints and mixed control-state constraints[END_REF] to the case of an arbitrary number of state constraints; see also Hager [24] in the case of first-order state constraints. Thanks to this alternative formulation, we are able to prove that both the multiplier associated with the state constraint and the optimal control are timecontinuous, and to obtain the expression of these two "algebraic variables" as fonctions of the state and alternative costate. Then, under some geometric hypotheses (finitely many junction points and growth condition for the state constraint), we obtain in section 4 no-gap second-order optimality conditions. Finally in section 5 we perform a sensitivity analysis w.r.t. the r.h.s. of the state equation. We obtain under weak assumptions the second-order expansion of the value function as well as the first-order expansion of the soltion path.

Framework

Consider the heat equation with distributed control over a subset of the domain:

y t -∆y + γy 3 = i ω u in Q, (1)
y = 0 over Σ, (2)
y(•, 0) = y 0 over Ω, (3)
where γ ∈ R, Ω is an open set of R n , n ∈ {2, 3} with C 2 -smooth boundary ∂Ω, Q := Ω × [0, T ], Σ := ∂Ω × [0, T ], ω is an open subset of Ω, Q ω = ω × [0, T ], u ∈ L 2 (Q ω )
is the control and T > 0 is the horizon. The function i ω is the injection from L 2 (Q ω ) into L 2 (Q) and the given initial state y 0 belongs to H 1 (Ω). Let

H 2,1 (Q) := {y ∈ L 2 (0, T, H 2 (Ω)); y t ∈ L 2 (Q)}, H 2,1 Σ (Q) := {y ∈ H 2,1 (Q); y = 0 over Σ}. (4)
Norms for functions depending on space variables, or depending only on time, will be denoted by simple bars, and the L s norms will be, if there is no ambiguity, denoted by subscript s. We recall the following inclusions (see Lions and Magenes [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF] and Adams [START_REF] Adams | Sobolev spaces[END_REF], resp.)

H 2,1 (Q) ⊂ C([0, T ], H 1 0 (Ω)) with compact inclusion, (5)
H 1 (Ω) ⊂ L 6 (Ω) and |z| 6 ≤ c 2 |∇z| 2 ,
for all z ∈ H 1 0 (Ω), when n ≤ 3. (6) We say that y ∈ H 2,1 (Q) is a state associated with u ∈ L 2 (Q) if (y, u) satisfies the state equation ( 1)-( 3). The following is known (see [START_REF] Lions | Contrôle des systèmes distribués singuliers[END_REF] and for a similar result in the case of boundary control [START_REF] Bonnans | Analysis and control of a nonlinear parabolic unstable system[END_REF]).

Lemma 1. For any given q ∈ L ∞ (0, T, L 3 (Ω)), f ∈ L 2 (Ω), and f 0 ∈ H 1 0 (Ω), the following equation has a unique solution z in H 2,1 Σ (Q):

z t -∆z + qz = f in Q, (7) y(•, 0) = f 0 over Ω, (8)
Proof. Set c 1 := |q| L ∞ (0,T,L 3 ) . Multiplying ( 7) by z, integrating over space, and adopting the convention to eliminate the space argument when no confusion can occur, get

1 2 d dt |z(t)| 2 2 + |∇z(t)| 2 = Ω (f (x, t)z(x, t) -q(x, t)z 2 (x, t))dxdt ≤ (|f (t)| 2 + c 1 |z(t)| 6 )|z(t)| 2 . (9)
By ( 6) and Hölder's inequality, for all ε > 0, we have that

2|z(t)| 6 |z(t)| 2 ≤ ε|z(t)| 2 6 + 1 ε |z(t)| 2 2 ≤ εc 2 2 |∇z(t)| 2 2 + 1 ε |z(t)| 2 2 . ( 10
)
Taking ε such that εc 1 c 2 2 = 1, obtain with (9)

d dt |z(t)| 2 2 + |∇z(t)| 2 ≤ 2|f (t)| 2 |z(t)| 2 + c 1 ε |z(t)| 2 2 ≤ |f (t)| 2 2 + 1 + c 1 ε |z(t)| 2 2 .
(11) By Gronwall's lemma, we have that for some c 3 > 0:

z 2 L ∞ (0,T,L 2 (Ω) + ∇z 2 2 ≤ c 3 |f 0 | 2 2 + f 2 L 2 (Q) . (12)
In particular when f = 0 and f 0 = 0, this proves the uniqueness property. Integrating [START_REF] Bonnans | Contrôle de systèmes elliptiques semilinéaires comportant des contraintes sur l'état[END_REF] over time and taking [START_REF] Bonnans | Contrôle de systèmes elliptiques semilinéaires comportant des contraintes sur l'état[END_REF] into account, we obtain an apriori estimate of z in the space

W (0, T ) := {z ∈ L ∞ (0, T, L 2 (Ω)); ∇z ∈ L 2 (0, T, L 2 (Ω))}, (13)
therefore also with (6) of z in L 2 (0, T, L 6 (Ω)), and finally of qz in L 2 (Ω). By [START_REF] Bebernes | A mathematical analysis of blowup for thermal reactions-the spatially nonhomogeneous case[END_REF] we get an apriori estimate of z in H 2,1 Σ (Q). The construction of the solution is then easily obtained by a Galerkin approximation. Lemma 2. For given u ∈ L 2 (Q ω ), either the state equation has a unique solution, or there exists a maximal time τ ∈ (0, T ] such that the state equation with time restricted to [0, τε] has, for all ε > 0, a unique solution, and |y(t)| 6 is not bounded over [0, τ ).

Proof. a) Existence of a solution of the state equation for small time. Let a 0 := |y 0 | 6 and, for b > a 0 , let ϕ b : R + → R + be defined by

ϕ(s) = 1 when s ≤ b, max(1 + b -s, 0) when s > b. ( 14
)
Consider the perturbed state equation

y t -∆y + γϕ(|y| 6 )y 3 = i ω u in Q, (15)
y(•, 0) = y 0 over Ω, [START_REF] Bonnans | Stability and sensitivity analysis for optimal control problems with a first-order state constraint[END_REF] with possible solution in H 2,1 Σ (Q). Since ϕ has value in [0, 1] and ϕ(a) = 0 when a ≥ b + 1, we have that ϕ(|y| 6 )y 3 ∈ L ∞ (0, T, L 2 (Ω)). It is therefore easy to obtain an a priori estimate of the solution in H 2,1 Σ (Q) and to construct a solution by Galerkin's method; when passing to the limit in the nonlinear term, we use the compact inclusion [START_REF] Arada | Minimax control of parabolic systems with state constraints[END_REF] and the continuity of ϕ. Since |y| 6 is a continuous function of time, we have that |y| 6 ≤ b when t ≤ t 1 , for some t 1 > 0. This proves that the state equation has a solution for t ≤ t 1 . b) Continuation. Let y be a solution of the state equation for time t ≤ t 0 . If the continuous function |y(t)| 6 is bounded, then y has a limit y(t 0 ). Applying the construction of step a with b > max 0≤t≤t0 |y(t)| 6 , we can construct a solution for t > t 0 , close to t 0 . This proves that either the solution is defined over [0, T ], or there exists τ such that the solution is defined over [0, τ ), and lim sup t↑τ |y(t)| 6 = +∞. c) Uniqueness of the solution. If y 1 and y 2 are two solution over (0, t 0 ), setting z := y 2y 1 and q = γ(3y 2 + 3yz + z 2 ), we obtain that q ∈ L ∞ (0, t 0 , L 3 (Ω)). Applying lemma 1 with f = 0 and f 0 = 0, and T = t 0 , obtain z = 0, which is the desired uniqueness property. In particular, let y i be solutions of the state equation for maximal time t ≤ t i , i = 1, 2. Then y 1 and y 2 coincide over (0, min(t 1 , t 2 )). If say t 1 < t 2 if follows that |y 1 (t)| 6 remains bounded for t ≤ t 1 which is impossible. Hence there is a unique solution defined up to its maximal time.

In view of the estimate

1 2 d dt |y u (t)| 2 + |∇y u (t)| 2 + γ|y u (t)| 4 4 = ω y u (x, t)u(x, t) dxdt, (17)
we have, when γ ≥ 0, given u ∈ L 2 (Q), an apriori estimate of y u in W (0, T ); using a Galerkin approximation argument, one easily deduces that when γ ≥ 0 the state equation has a unique solution y u , and that there exists c > 0 not depending on γ such that y H 2,1 (Q) ≤ c u 2 . When γ < 0, the state equation has at most one solution. The strong solution, always defined for small time, can explode in a finite time, see for similar results Bebernes and Kassoy [START_REF] Bebernes | A mathematical analysis of blowup for thermal reactions-the spatially nonhomogeneous case[END_REF] and Tartar [START_REF] Tartar | Topics in nonlinear analysis[END_REF]. In any case we denote by y u the solution. The implicit function theorem can be applied to the state equation (see [START_REF] Bonnans | Contrôle optimal de systèmes instables gouvernés par des équations aux dérivées partielles non linéaires[END_REF][START_REF] Lions | Contrôle des systèmes distribués singuliers[END_REF]), and since the latter is of class C ∞ , we obtain that the mapping u → y u is of class C ∞ . The associated cost function is, for some N > 0:

J(u, y) = 1 2 Q (y(x, t) -y d (x, t)) 2 dxdt + N 2 Qω u 2 (x, t)dxdt. (18)
For some C > 1 2 Ω |y 0 (x)| 2 dx, define g :

L 2 (Ω) → R by g(z) := 1 2 Ω |z(x)| 2 dx - C. Consider the state constraint g(y(•, t)) ≤ 0. ( 19
)
The paper is concerned with the discussion of optimality conditions of the following optimal control problem:

Min (u,y)∈L 2 (Qω)×H 2,1 (Q)
J(u, y) subject to (1)-( 3) and [START_REF] Casas | Pontryagin's principle for local solutions of control problems with mixed control-state constraints[END_REF].

(P )

Remark 3. When γ ≥ 0, problem (P ) is always feasible and has a nonempty set of solutions. Indeed, it follows from (17) that t → g(y u )(t) is nonincreasing when γ ≥ 0 and u = 0. In particular, u = 0 is feasible when γ ≥ 0. The existence of a solution to (P ) can be easily proved by passing to the limit in a minimizing sequence, since when γ ≥ 0, by [START_REF] Bonnans | No gap second order optimality conditions for optimal control problems with a single state constraint and control[END_REF], a minimizing sequence of ocontrols is bounded in L 2 (Q ω ) and hence the corresponding sequence of states is bounded in y in H 2,1 (Q) When γ < 0 the existence of a solution is unclear.

Of course it makes sense to study optimality conditions even if we cannot prove the existence of a solution.

A trajectory of (P ) is an element (u, y) ∈ L 2 (Q ω ) × H 2,1 (Q) satisfying the state equation ( 1)- [START_REF] Arada | Optimal control problems with mixed control-state constraints[END_REF]. A feasible trajectory is one satisfying the state constraint [START_REF] Casas | Pontryagin's principle for local solutions of control problems with mixed control-state constraints[END_REF]. The set of feasible trajectories is denoted F (P ).

We say that (ū, ȳ) is a local solution of (P ) that satisfies the quadratic growth condition with parameter θ ∈ R, if it belongs to F (P ) and there exists ρ > 0 such that

J(ū, ȳ) ≥ J(u, y) + θ|ū -u| 2 L 2 (Qω) if (u, y) ∈ F (P ) and |u -ū| L 2 (Qω ) ≤ ρ. ( 20
)
If this holds for θ = 0, we say that (ū, ȳ) is a local solution of (P ). We say that (ū, ȳ) ∈ F (P ) satisfies the quadratic growth condition if [START_REF] Casas | Pontryagin's principle for state-constrained boundary control problems of semilinear parabolic equations[END_REF] holds for some θ > 0 and ρ > 0. For (u, y) ∈ F (P ), define the contact set by

I(g(y)) = {t ∈ [0, T ]; g(y)(t) = 0} . (21)
For t ∈ [0, T ], if g(y)(t) = 0 (resp. g(y)(t) < 0), we say that the constraint is active (resp. inactive) at time t. Finally we denote by χ ω the restriction L 2 (Ω) → L 2 (ω) (which is nothing else than the transposition of i ω ).

3 First-order analysis

The unqualified optimality system

We first put problem (P ) under a compact format by defining

     G : L 2 (Q ω ) × y ∈ H 2,1 (Q) s.t. y |Σ = 0 → L 2 (Q) × H 1 0 (Ω), G(u, y) := y t -∆y + γy 3 -i ω u y(•, 0) -y 0 . ( 22
)
The linearized state equation at point (u, y) is defined as

z t -∆z + 3γy 2 u = i ω v in Q; z = 0 on Σ, z(•, 0) = 0. ( 23
)
This linearized state equation is well-posed in the sense that, with each v ∈

L 2 (Q ω ) is associated a unique solution z ∈ H 2,1 (Q). Since G is a C ∞ mapping, it follows from the Implicit Function Theorem that the mapping u → y u is of class C ∞ : L 2 (Q ω ) → H 2,1 (Q),
and the directional derivative of the state at point u in direction v is the solution z of the linearized state equation. Next define J :

L 2 (Q ω ) → R and G : L 2 (Q ω ) → C([0, T ]) by J (u) := J(u, y u ); G(u)(t) := g(y u (t)) = 1 2 |y u (t)| 2 -C. (24)
We can rewrite problem (P ) under the form of an "abstract problem":

Min u J (u); G(u) ∈ K, (AP )
where

K = C([0, T ]) -is the cone of continuous nonpositive functions over [0, T ]. We recall that the topological dual of C([0, T ]) is the set M (0, T ) of regular Borel measures over [0, 1]. Therefore the (negative) polar cone of K is K -= M ([0, T ]) + .
We know (see e.g. [18, Example 2.62 and 2.63]) that since K is a cone, its normal cone at a point h ∈ K is N K (h) = K -∩ h ⊥ and that for h ∈ K, denoting by supp(µ) the support of a measure µ:

N K (h) = {µ ∈ M (0, T ) + ; supp(µ) ⊂ h -1 (0)}. ( 25
)
Define (see e.g. [18, Section 3.1.2]) the generalized Lagrangian of problem (AP ) as L :

L 2 (Q ω ) × R × M ([0, T ]) such that L(u, α, µ) := αJ (u) + µ, G(u) , (26)
and the set of generalized Lagrange multipliers associated with u ∈ F (AP ) as

Λ g (u) := {(α, µ) ∈ R + × N K (G(u)); (α, µ) = 0; D u L(u, α, µ) = 0}. ( 27
)
This set is a (possibly empty) convex cone.

Theorem 4. With a local solution (u, y) of (P ) is associated a non empty set of generalized Lagrange multipliers.

Proof. This results from [START_REF] Bonnans | Perturbation analysis of optimization problems[END_REF]Prop. 3.18]. For applying this result it suffices to verify that the cone defined in equation (3.24) of [START_REF] Bonnans | Perturbation analysis of optimization problems[END_REF] has a nonempty relative interior, which trivially holds since K has a nonempty interior.

We next give a more explicit form to the optimality conditions. The condi-

tion µ ∈ N K (G(u) is equivalent in view of (25) to g(y) ≤ 0, µ ≥ 0, T 0 g(y(t))dµ(t) = 0. ( 28
)
By [START_REF] Kawasaki | An envelope-like effect of infinitely many inequality constraints on second order necessary conditions for minimization problems[END_REF], D u L(u, α, µ) = 0 iff, for y = y u :

αN Qω uv + α Q (y -y d )z v + T 0 Ω yz v dµ = 0, for all v ∈ L 2 (Q ω ), ( 29
)
where z v is the solution of the linearized state equation [START_REF] Dunford | Linear operators, Vol I and II[END_REF]. For the sake of simplicity we will assume in the sequel of this paper that the state constraint is not active at the final time. Define the costate p as the solution in L 2 (Q) of the linear equations

Q [α(y -y d )z + p(∆z -3γy 2 z -z t )]dxdt + Q y(x, t)z(x, t)dxdµ(t) = 0, for all z ∈ Z := {H 2,1 Σ (Q); z(•, 0) = 0}. ( 30
) Setting f = z t -∆z + 3γy 2 z, this is equivalent to Q pf = α Q (y -y d )z + Q y(x, t)z(x, t)dxdµ(t). (31)
For given f ∈ L 2 (Q), in view of lemma 1, there exists a unique z ∈ Z such that z t -∆z + 3γy 2 z = f . Therefore the r.h.s. of ( 31) may be interpreted as a continuous linear form over L 2 (Q). By the Riesz theorem there exists a unique p ∈ L 2 (Q) such that [START_REF] Maurer | On the minimum principle for optimal control problems with state constraints[END_REF] holds for all f ∈ L 2 (Q), or equivalently, (30) has a unique solution p ∈ L 2 (Q). Let D(Q) denote the set of C ∞ functions in Q with compact support, and let D ′ (Q) be the associated dual set of distributions. We see that ( 30) is (formally for the two last equations) equivalent to

-p t -∆p + 3γy 2 p = α(y -y d ) + ydµ(t) in D ′ (Q), (32)
p(•, T ) = 0, ( 33
) p = 0 on Σ. ( 34
)
We can justify ( 33) and ( 34) by observing that, since p is solution of a backwards heat equation with r.h.s. measure, by [START_REF] Casas | Pontryagin's principle for state-constrained boundary control problems of semilinear parabolic equations[END_REF], it belongs to the space L r (0, T, W 1,q (Ω)), for all r, q in [1, 2) such that 2/r + n/p > n + 1, and hence, its traces at time T and on Σ are well-defined. Actually we will derive in a direct manner a stronger regularity result in section 3.3. When f = i ω v, for some v ∈ L 2 (Q ω ), substracting (30) from [START_REF] Lions | Contrôle des systèmes distribués singuliers[END_REF], and taking [START_REF] Dunford | Linear operators, Vol I and II[END_REF] into account, obtain Qω (αN uv + pv) = 0, for all v ∈ L 2 (Q ω ), or equivalently αN u + χ ω p = 0 a.e. over Q ω .

Finally we have proved that

Λ g (u) := (α, µ) ∈ R + × N K (G(u)); (α, µ) = 0; χ ω p = -αN u,
where p is solution of ( 30) .

Denote the set of singular multipliers and Lagrange multipliers, resp., by

Λ s (u) := {µ ∈ N K (G(u)); µ = 0; D u L(u, 0, µ) = 0}, (37) Λ ℓ (u) := {µ ∈ N K (G(u)); D u L(u, 1, µ) = 0}. ( 38
)
A nonsingular multiplier is of the form α(1, µ ′ ) where µ ′ := µ/α is a Lagrange multiplier. Denote by supp(µ) ⊂ [0, T ] the support of the measure µ, and let the singular contact set be defined by

I s (u) = {t ∈ I(g(y u )); y u (•, t) = 0 a.e. on ω}. (39)
In view of (17), we see that the singular contact set is the set of times for which the constraint is active and the control has no influence on its time derivative. We next characterize the absence of singular mulitpliers using the following concept. We will say that the state constraint is nondegenerate at point u if I s (u) = ∅; otherwise we say that the state constraint is degenerate.

Theorem 5. Let (u, y) be a feasible point of (P ). Then Λ s (u) is empty iff the state constraint is nondegenerate at point u.

Proof. Let µ ∈ Λ s (u). By [START_REF] Tartar | Topics in nonlinear analysis[END_REF], p = 0 over Q ω . Then by [START_REF] Raymond | Pontryagin's principle for state-constrained control problems governed by parabolic equations with unbounded controls[END_REF], ydµ = 0 over Q ω , which means that supp(µ) ⊂ I s (u). Since µ = 0, the state constraint is degenerate.

Conversely, if the state constraint is degenerate, there exists τ ∈ I s (u). Let µ be the Dirac measure at time τ . Then the costate equation is satisfied, as well as [START_REF] Tartar | Topics in nonlinear analysis[END_REF], by p = 0. Therefore µ is a singular multiplier.

Constraint qualification condition

By a constraint qualification condition one usually understands a condition implying that Lagrange multipliers associated to a local solution exist. Combining theorems 4 and 5 we obtain that nondegeneracy of the constraint is a qualification condition. Yet it is of interest to relate it to the standard Robinson constraint qualification (see [START_REF] Robinson | First order conditions for general nonlinear optimization[END_REF], [START_REF] Robinson | Stability theorems for systems of inequalities, part II: Differentiable nonlinear systems[END_REF]). Since C([0, T ]) -has a non empty interior, the latter writes for a given (u, y) ∈ F (P ): [START_REF] Dunford | Linear operators, Vol I and II[END_REF], such that g ′ (y(t))z(t) = ω y(x, t)z(x, t)dx < 0, for all t ∈ I(g(y)).

There exists (v, z) ∈ L 2 (Q ω ) × H 2,1 (Q), solution of
(40) Remark 6. (i) It is know that the constraint qualification implies that the set of Lagrange multipliers associated with a local solution is non empty and bounded. (ii) In addition, in the case of a constraint to belong to a convex set with a nonempty interior, then by [START_REF] Zowe | Regularity and stability for the mathematical programming problem in Banach spaces[END_REF] (see also [START_REF] Bonnans | Perturbation analysis of optimization problems[END_REF]Prop. 3.17]), constraint qualification holds at a local solution iff the set of Lagrange multipliers is non empty and bounded.

Theorem 7. Let (u, y) ∈ F (P ). Then the state constraint is not degenerate iff the constraint qualification condition (40) holds.

The proof is based on the following lemma, where we show how to control the time derivative of the linearized state constraint. We first need to set (here C is the constant used in the definition of the state constraint):

T (y) := {t ∈ [0, T ]; ω y 2 (x, t)dx ≥ C, (41)
κ(t) := ω y 2 (x, t)dx -1 if t ∈ T (y), 0 otherwise, ( 42
)
and for z ∈ H 2,1 (Q) η(z, t) := -κ(t) ω z(x, t)u(x, t)dt -2 Ω ∇y(x, t)∇z(x, t)dt +4γ Ω y 3 (x, t)z(x, t)dt . ( 43
) Lemma 8. Let (u, y) ∈ L 2 (Q) × H 2,1 (Q). Then (i) For all f ∈ L 2 (Q), the problem z t -∆z + 3γy 2 z = η(z, t)χ ω y(x, t) + f (x, t) (44) z(x, t 1 ) = z 0 , ∀x ∈ Ω, ( 45
)
has a unique solution in H 2,1 Σ (Q), the mapping f → z is continuous, and

d dt g ′ (y)z = Ω f y over T (y). (46)
(ii) In particular, given ν ∈ L 2 (0, T ), when f = ν(t)κ(t)χ ω y(x, t), we have that

f ∈ L 2 (Q) and d dt g ′ (y)z = ν(t) over T (y). ( 47
)
Proof. (i) For a ∈ H 1 (0, T ), nonnegative and nondecreasing, and z in H 2,1 (Q), set w = e -a(t) z. Then w ∈ H 2,1 (Q), and since z t = e a(t) (w t + ȧw), z is solution of (45) iff w is solution of

w t + ȧw -∆w + 3γy 2 w = η(w, t)χ ω y + e -a(t) f, (48)
w(x, 0) = e -a(t) z 0 , ∀x ∈ Ω. ( 49
)
Multiplying both sides of (48) by w and integrating over space, obtain (skipping arguments of functions)

1 2 d dt |w| 2 2 + Ω ȧ(t)w 2 + |∇w| 2 + 3γy 2 w 2 = + κ(t) ω wu -2 Ω ∇y∇w + 2γ Ω y 3 w ω yw + e -a(t) Ω f w. (50)
Using

a ≥ 0, y ∈ H 2,1 (Q) ⊂ C([0, T ], H 1 (Ω)) ⊂ C([0, T ], L 6 (Ω)), κ(t) ≤ 2C,
and the Cauchy-Schwartz inequality, obtain that for some c 1 > 0:

1 2 d dt |w| 2 2 + ȧ(t)|w| 2 2 + |∇w| 2 2 ≤ c 1 |w| 2 2 |u| 2 + α|∇w| 2 2 + 1 α |w| 2 2 + |w| 2 2 + 1 2 |f | 2 2 + |w| 2 2 .
(51)

Setting α := 1 2c1 , obtain, for some c 2 > 0:

d dt |w| 2 2 + 1 2 |∇w| 2 2 ≤ |w| 2 2 c 1 |u| L 2 (ω) + c 2 -ȧ(t) + |f | 2 2 . ( 52
)
Choosing

a(t) = c 1 t t1 |u| 2 + c 2 (t -t 1 )
, for all t ∈ (0, T ), (this is indeed in H 1 (0, T ), nonnegative and nondecreasing), obtain after a time integration

|w(t 2 )| 2 2 + ∇w(t 2 ) 2 L 2 (Ω×(t1,t2)) ≤ |z 0 | 2 2 + f 2 L 2 (Ω×(t1,t2)) . (53)
So we have an a-priori estimate of

w in L 2 ([0, T ], H 1 0 (Ω)) L ∞ ([0, T ], L 2 (Ω)
). Let us show that (48) gives an a-priori estimate of w t -∆w in L 2 (Q). The a-priori estimate of w, and hence, of z, in H 2,1 (Q) will follow as well as the conclusion of the lemma. Indeed, since

η(w, t)χ ω y 2 L 2 (Q) ≤ y 2 L ∞ (0,T,L 2 (Ω)) T 0 η(w, t) 2 dt, (54)
it suffices to prove that η ∈ L 2 (0, T ). And this follows from

η(w, t) 2 ≤ C w 2 L ∞ (0,T,L 2 (Ω)) |u(t) 2 2 + ∇y 2 L ∞ (0,T,L 2 (Ω)) |∇w| 2 2 + y 6 L ∞ (0,T,L 6 (Ω)) |w| 2 2 . . (55)
Finally (46) follows from the fact that d dt g ′ (y)z = Ω (y t z + yz t ), using the expressions of y t and z t , and (ii) is an easy consequence of (i).

Proof of theorem 7. If the state constraint is degenerate, the constraint qualification (40) obviously cannot hold. Assume now that the state constraint is nondegenerate. Set

I ε := {t ∈ [0, T ]; dist(t, I(g(y))) ≤ ε}; J ε := [0, T ] \ I ε . ( 56
)
Since g(y(t)) is continuous over [0, T ], we may take ε > 0 so small that I ε ⊂ T , The set J ε is a finite union of relatively open subsets. Therefore I ε is the union finitely many closed intervals of the form [a iε, b i + ε] ∩ [0, T ], i = 1, . . . , N , and we have that

I(g(y)) ⊂ ∪ N i=1 [a i , b i ].
We next construct the perturbation v as follows. Take v = 0 over J ε . Next, given v and z over [0, (a

i -ε) + ], fix v over [a i -ε, b i + ε] such that , setting ξ(t) := g ′ (y(t))z(t) is affine over [a i -ε, a i ] and has value -1 over [a i , b i + ε].
In view of lemma 8(ii), it suffices to take z solution of ( 44)-( 45), with f = ν(t)κ(t)χ ω y(x, t), and with

ν(t) := -(1 + ξ(a i -ε))/ε, t ∈ (a i -ε, a i ), 0, t ∈ (a i , b i + ε). (57)
The related control is v = η(z, t)χ ω y + f .

In the sequel, we will say that (u, y) is a qualified feasible point of (P ) if it satisfies (40), and a regular extremal if in addition it has an associated Lagrange multiplier.

Alternative formulation of optimality conditions

In this section we will assume that the state constraint is nondegenerate.

We recall the following integration by parts formula, see [23, Vol. I, p. 154]:

Lemma 9. Let a and b be two functions of bounded variations in [0, T ]. Suppose that one is continuous, and the other is right-continuous. Then

a(T -)b(T -) -a(0 + )b(0 + ) = T 0 a(t)db(t) + T 0 b(t)da(t). (58)
We apply this lemma to the expression Q y(x, t)z(x, t)dxdµ(t) appearing in the costate equation [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF]. The function a(t) := Ω y(x, t)z(x, t)dx, with z ∈ H 2,1 (Q), is absolutely continuous; its derivative is

ȧ(t) = Ω [y t (x, t)z(x, t) + y(x, t)z t (x, t)]dx. ( 59
)
Since z(0) = 0 and taking the convention that µ(T ) = 0, we see that the l.h.s. in (58) is zero, and hence,

T 0 Ω y(x, t)z(x, t)dxdµ(t) = - T 0 Ω [y t (x, t)z(x, t)) + y(x, t)z t (x, t)]µ(t)dxdt (60)
We now introduce the (first) alternative costate defined as the sum of the costate and of the product of the measure by the derivative (w.r.t. the state) of the state constraint; in other words

p 1 := p + g ′ (y)µ = p + yµ in L 2 (Q). (61)
Substracting (60) from ( 30), obtain

Q p 1 z t dxdt = Q [(y -y d ) -3γy 2 p -µy t ]zdxdt + Q p∆zdxdt for all z ∈ Z := {z ∈ H 2,1 (Q); z(•, 0) = 0; z(x, t) = 0 on Σ}. ( 62
)
Eliminating p = p 1yµ and using

Q y∆zdxdt = Q z∆ydxdt, we get Q p 1 z t -∆z + 3γy 2 z dxdt = Q [(y -y d ) + µ(3γy 3 -∆y -y t )]zdxdt, for all z ∈ Z := {z ∈ H 2,1 (Q); z(•, 0) = 0; z(x, t) = 0 on Σ}.
(63) This is a classical adjoint equation with L 2 r.h.s.; using the expression of y t , obtain the following result.

Lemma 10. Let (u, y) be a regular extremal of (P ). The alternative multiplier

p 1 = p + yµ is the unique solution in H 2,1 (Q) of the equation -p 1 t -∆p 1 + 3γy 2 p 1 = y -y d -(2∆y -6γy 3 + i ω u)µ in Q, (64) p 1 (•, T ) = 0, (65) p 1 (•, t) = 0 on Σ. (66)
Then, we have the following proposition stating the alternative formulation of the optimality conditions: Proposition 11. Let (u, y) ∈ F (P ) be qualified. Then with (u, y) is associated a Lagrange multiplier iff there exists a pair (p 1 , µ) ∈ H 2,1 (Q) × M + ([0, T ]), satisfying the state equation (1)-( 3), the alternative costate equation (64)-(66), the complementarity conditions (28) for the state constraint, and the relation

N u + χ ω (p 1 -µy) = 0 a.e. on ω × [0, T ]. ( 67
)
Proof. The necessity of the alternative costate equation ( 64)-(66) has been obtained in the above discussion. Eliminating p = p 1yµ from [START_REF] Tartar | Topics in nonlinear analysis[END_REF], we obtain (67). Conversely, substituting p 1 = p + yµ in the alternative optimality system, and using (60), we recover the "classical" adjoint equation ( 32)-( 34) as well as [START_REF] Tartar | Topics in nonlinear analysis[END_REF].

We now prove the following regularity result.

Lemma 12. Let (u, y) be a regular extremal of (P ) and (p, p 1 , µ) the classical and alternative costate and the multiplier associated with the state constraint. Then (i) µ is a continuous function of time and u ∈ C([0, T ];

H 1 (ω)), (ii) if at
Optimal control of a time-dependent state constrained parabolic equation 13 time t the state constraint is active and ω y 2 (x, t) dx = 0, then

0 = d dt (g(y))(t) = -|∇y(t)| 2 + ω u(x, t)y(x, t)dx -γ|y(t)| 4 4 , (68)
µ(t) = N |∇y(t)| 2 + γN |y(t)| 4 4 + ω y(x, t)p 1 (x, t)dx |χ ω y(t)| 2 2 (69) u = 1 N χ ω N |∇y(t)| 2 + γN |y(t)| 4 4 + ω y(x, t)p 1 (x, t)dx |χ ω y(t)| 2 2 y -p 1 . ( 70
)
Proof of lemma 12. (i) By (67), and since both y and p 1 belong to

H 2,1 (Q) ⊂ C([0, T ]; H 1 (Ω))
, it suffices to prove that µ is continuous to obtain the conclusion. Since µ has left and right limits, and p 1 as well as y belong to C(0, T, H 1 0 (Ω)), the same holds for u by (67); the latter implies also, denoting by [•] the jump function (e.g., [u]

(t) := u(t + ) -u(t -)), that N [u] = [µ]χ ω y.
Taking the scalar product of both sides by [u], get

N |[u]| 2 L 2 (ω) = [µ] ω [u]ydx. ( 71
)
If [µ](t) = 0 for some time t, then the constraint is active and attains its maximum at time t. The formal expression of the derivative of g(y) is

d dt (g(y))(t) = Ω y∆y + ω uy -γ Ω y 4 = - Ω |∇y| 2 + ω uy -γ Ω y 4 . ( 72
)
Since y ∈ C(0, T, H 1 (Ω)) the first and the last terms in the r.h.s. are continuous (we recall that H 1 (Ω) ⊂ L 6 (Ω) with continuous embedding in dimension 2 or 3); since u has left and right limits, so has d dt g(y) and they are given by ( 72), where u is the left or right limit, so that the jump is

d dt g(y)(t) = ω [u]y = N [µ] |[u]| 2 L 2 (ω) . ( 73
)
Since g(y) attains a maximum, [ d dt (g(y))(t)] ≤ 0, and hence, [u] = 0, contradicting [µ] = 0. We have proved point (i) as well as (68). (ii) Elimining the control variable from (67) and using (68), get (69). Eliminating then µ thanks to (69) in (67), obtain (70).

Lemma 13. Let (u, y) be a regular extremal of (P ). Then (p, p 1 , µ) the classical and alternative costate and the multiplier associated with the state constraint are unique.

Proof. Let (p 1 , p 1 1 , µ 1 ) and (p 2 , p 1 2 , µ 2 ) two triples of classical and alternative costates and multipliers associated with the state constraint. Set

(p, p1 , μ) = (p 1 -p 2 , p 1 1 -p 1 2 , µ 1 -µ 2 ). ( 74
)
Then p1 is solution of:

-p1 t -∆p 1 + 3γy 2 p1 = -(2∆y -6γy 3 + i ω u)μ in Q, ( 75
) p(•, T ) = 0, (76) p(•, t) = 0 on Σ. ( 77
)
So if μ = 0 it follows that p1 = 0 and these two relations imply p = 0. Therefore it suffices to prove that μ = 0. In view of the abstract form of the optimality system, we have that μ, G ′ (u)v = 0 for all v ∈ L 2 (Q ω ). Take ε > 0 such that the state constraint is not active for time t ≤ ε. The support of μ is concentrated over [ε, T ]. By lemma 8 we can make G ′ (u)v equal to an arbitrary function of H 1 (ε, T )), and since the inclusion of

H 1 (ε, T ) in C([ε, T ]) is dense, μ, G ′ (u)v = 0 for all v ∈ L 2 (Q ω ) implies μ = 0.
We next prove that when state constraint is active over an interval, the Cantor part of the measure is also null. Lemma 14. Let (u, y) be a regular extremal of (P ). Assume that the state constraint is active over an interval [t 1 , t 2 ], where

0 ≤ t 1 < t 2 ≤ T . Then µ is absolutely continuous over [t 1 , t 2 ].
Proof. It suffices to prove that the expression in the r.h.s. of (69) belongs to W 1,1 (t 1 , t 2 ). Since by assumption the denominator is far from zero and belongs to W 1,1 (t 1 , t 2 ), it suffices to check it for each term of the numerator. This is obviously true for the last two terms, and for the first one we observe that for t ∈ (t 1 , t 2 ):

|∇y(t)| 2 = |∇y(t 1 )| 2 - t t1 Ω ∆y(x, s)y t (x, s)dxdt (78)
Indeed this holds when y is smooth enough, and since each side of the equality is a continuous function in H 2,1 (Q), and the latter is the closure (for its norm) of the set of C ∞ functions over Q, (78) holds. The conclusion follows.

Arc based alternative formulation

We next present a variant of the alternative formulation that is useful when formulating shooting algorithms; see Bonnans and Hermant [START_REF] Bonnans | Well-posedness of the shooting algorithm for state constrained optimal control problems with a single constraint and control[END_REF] for related results in a finite dimensional setting. Let us start with some definitions. A boundary arc (resp. interior arc) is a maximal interval of positive measure I such that g(y)(t) = 0 (resp. g(y)(t) < 0), for all t ∈ I. Left and right endpoints of a boundary arc [τ en , τ ex ] are called entry and exit point, respectively. A touch point τ to is an isolated contact point, satisfying g(y)(τ to ) = 0 and g(y)(t) < 0 for t = τ to in the neighborhood of τ to . The endpoints of interior arcs belonging to (0, T ) are called junction points (or times). If the set of junction points of a trajectory is finite, then it is of the form

T = T en ∪ T ex ∪ T to , (79)
with T en , T ex , T to the disjoint (and possible empty) subsets of respectively entry, exit and touch points. We denote by I b the union of boundary arcs, i.e.

I b = N b i=1 [τ en i , τ ex i ] for T en := τ en 1 < • • • < τ en N b
and similar definition of T ex , and we have I(g(y)) = T to ∪ I b . Definition 15. We say that a trajectory (u, y) is solution of the alternative formulation, if there exist p 1 ∈ L 2 (Q), µ 1 ∈ BV (0, T ) and alternative jump in the constraint statement (in the context of our state constrained optimal control problems, the cone C([0, T ]) -of nonpositive continuous functions over [0, T ]). This was set in a general setting by Cominetti [START_REF] Cominetti | Tangent sets to unilateral convex sets[END_REF], and examples of computations of the sigma-term for various function spaces, including those of continuous functions over a compact set, are provided in Cominetti and Penot [START_REF] Cominetti | Tangent sets to unilateral convex sets[END_REF].

The main result of Bonnans and Hermant [START_REF] Bonnans | No gap second order optimality conditions for optimal control problems with a single state constraint and control[END_REF] is that, under weak conditions, for the optimal control of o.d.e.'s, this sigma-term reduces to the contribution of isolated contact points (called touch points), that are equal to the product of the jump of the multiplier times a quadratic form of the critical direction at the touch point. In addition, for first-order state constraints, the contribution of "regular" touch points is also zero. We will obtain a similar result in our context. In fact we improve here the method of [START_REF] Bonnans | No gap second order optimality conditions for optimal control problems with a single state constraint and control[END_REF] by avoiding any hypothesis on the second-derivative of the state constraint at junction points.

Th Lagrangian of problem (P ) (in the formalism where both the state and control variables appear) in qualified form is the function

L : L 2 (Q ω )×H 2,1 Σ (Q)× L 2 (Q) × H 1 0 (Ω) × M ([0, T ]) defined by L(u, y, p, q, µ) := J(u, y) + Q p ∆y -γy 3 + i ω u -y t dxdt + T 0 g(y(t))dµ(t) + Ω q(x)(y(x, 0) -y 0 (x))dx. (97)
Its second-order directional derivative in direction (v, z) is

∆(v, z) := N v 2 2 + Q (1 -6γp(x, t)y(x, t)) z(x, t) 2 dxdt + T 0 |z(t)| 2 2 dµ(t). (98) We say that (v, z) ∈ L 2 (Q) × H 2,1
Σ (Q) is a critical direction if (i) it satisfies the linearized state equation ( 23) and (ii) it is tangent to the state constraints and satisfies the complementarity condition with the multiplier µ, i.e., g ′ (y(t))z(t) ≤ 0 over I(g(y)) ,

g ′ (y(t))z(t) = 0 over supp(µ).

Denote by C(u, y) the set of critical directions associated with (u, y) ∈ F (P ). We say that the contact set has a finite structure if it is a finite union of touch points and boundary arcs. In that case, we will say that the hypothesis of strict complementarity holds if the support of the measure dµ is the union of the boundary arcs. In that case, a direction (v, z) satisfying ( 23) is critical iff it satisfies g ′ (y(t))z(t) = 0 over boundary arcs, g ′ (y(τ ))z(τ )) ≤ 0 for each touch point τ .

(101)

Theorem 17. Let (u, y) be a qualified local solution of (P ), with associated multiplier µ and costate p. If the contact set has a finite structure and the hypothesis of strict complementarity holds, then

∆(v, z) ≥ 0, for all (v, z) ∈ C(u, y). (102)
The proof will use the following corollary of lemma 8.

We next have to distinguish two cases.

(i) The last connected subset of I(g(y)) is the touch point τ . Take

ν(t) := -2ε(c 2 + 1)/(τ -t 0 ) when t < τ , 0 when t ∈ [τ, T ]. (109)
Since (g ′ (y)z) (τ ) = 0, we have then (g ′ (y)z ε ) (τ ) ≤ -ε, so that for small σ > 0, g(y)+ g ′ (y)z ε ≤ 0 over [t 0 , T ], and at the same time we have constructed v ε such that v εv 2 + z εz 2,1 = O(ε). So we have constructed a radial direction (v ε , z ε ) converging to (v, z).

(ii) The last connected subset of I(g(y)) is the boundary arc [τ, τ ′ ]. As a first step, take ν(t) equal to some constant c 3 so that (g ′ (y)z) (τε) = 0. Using the relation similar to (107) but at time τε, get

0 = (g ′ (y)z ε ) (τ -ε) = c 3 (τ -ε-t 0 )+(g ′ (y)(z ε -z)) (t 0 )+(g ′ (y)z) (τ -ε). ( 110
)
Since g ′ (y)z is a continuous fonction, we have that |(g ′ (y)z) (τε)| = o(1), and hence for ε ≤ 1 2 (τt 0 ), it suffices to take c 3 = o(1). The second step is to set (g ′ (y)z) (t) = 0, for all t ∈ [τε, τ ], i.e., ν(t) = -d dt (g ′ (y)z) (t), and hence, ν L 2 (τ -ε,τ ) = o(1). On [τ, τ ′ ] we take of course ν = 0, and if τ ′ < T , we set ν so that (g ′ (y)z) (t) = 0, for all t ∈ [τ ′ , τ ′ + ε], which similarly implies ν L 2 (τ ′ ,τ +ε) = o(1). So again we have constructed a radial critical direction converging to (v, z).

Second-order sufficient optimality condition

We obtain in this section no-gap second-order optimality conditions. Theorem 19. Let (u, y) be a regular extremal of (P ). Then a sufficient condition for the quadratic growth condition [START_REF] Casas | Pontryagin's principle for state-constrained boundary control problems of semilinear parabolic equations[END_REF] is

∆(v, z) > 0, for all (v, z) ∈ C(u, y) \ {0}. (111)
If, in addition, the contact set has a finite structure, then (111) is a necessary condition for quadratic growth.

Proof. a) The sufficiency condition follows from [START_REF] Bonnans | Perturbation analysis of optimization problems[END_REF]Section 3.3]; for the sake of completeness, we give a direct proof. Assume that (111) holds, and that a feasible sequence (u k , y k ) for (P ) satisfies (u k , y k ) = (u, y) for all k, (u k , y k ) → (u, y), and

J (u k ) = J(u k , y k ) ≤ J(u, y) + o( u k -u 2 2 ). ( 112
)
Set

β k := u k -u 2 and v k := (u k -u)/β k . Then u k = u+β k v k . Since v k 2 = 1,
extracting if necessary a subsequence, we may assume the existence of

v ∈ L 2 (Ω) such that v k ⇀ v (weak convergence). Since J (u k ) = J (u)+β k J ′ (u)v k +o(β k ), we deduce from (112) that J ′ (u)v ≤ 0. Similarly, since T K (G(u)) ⊃ K -G(u k ) β k = G ′ (u)v k + o(1), (113)
and since T K (G(u)) is a (weakly) closed convex set, we obtain that G ′ (u)v ∈ T K (G(u)). It follows that v is a critical direction. Using L(u k , y k , p, q, µ) ≤ J(u k , y k ) and L(u, y, p, q, µ) = J(u, y), deduce from (112) that 0 ≥ L(u k , y k , p, q, µ) -L(u, y, p, q, µ)

+ o(β 2 k ) = 1 2 β 2 k ∆(v k , z k ) + o(β 2 k ). ( 114
)
Since ∆(•, •) is w.l.s.c., obtain ∆(v, z) ≤ 0, which by (111), since (v, z) is a critical direction, implies (v, z) = 0, and hence, ∆

(v k , z k ) → ∆(v, z). Since the mapping v → z (solution of the linearized state equation) is compact L 2 (Q ω ) → L 2 (Q), it follows that v k 2 → v 2
, and therefore v k → v (for the strong topology), contradicting v k 2 = 1. b) Let (u, y) ∈ F (P ) be qualified, have a finite structure, and satisfy the quadratic growth condition. Then for ε > 0 small enough (u, y) is a local solution of the same problem with N changed into N -1 2 ε. Since this does not change the expression of the first-order optimality conditions at the point (u, y), the costate p and multiplier µ are the same, while the perturbed Hessian of Lagrangian is ∆

ε (v, z) = ∆(v, z) -ε v 2 2 .
Since ∆ ε (v, z) ≥ 0 by theorem 19, the conclusion follows.

Sensitivity analysis

We consider now a family of optimal control problems parameterized by an additional f ∈ L 2 (Q) on the r.h.s. of the state equation. The perturbed state equation is therefore

y t -∆y + γy 3 = f + i ω u in Q, (115) y = 0 over Σ, (116) y
(•, 0) = y 0 over Ω. (117)
Let (ū, ȳ) be a local solution of (P ) satisfying the quadratic growth condition [START_REF] Casas | Pontryagin's principle for state-constrained boundary control problems of semilinear parabolic equations[END_REF] for some θ > 0 and ρ > 0. Assume that they satisfy the qualification condition, and let (p, μ) denote the associated costate and Lagrange multiplier. Consider the localizing constraint

u -ū 2 ≤ ρ. (118)
The perturbed optimal control problem is Min (u,y)∈L 2 (Qω )×H 2,1 (Q) J(u, y) subject to (115)-( 117) and ( 118). (P f )

Denote by v(f ) the value of problem (P f ). Using the methodology of [START_REF] Bonnans | Perturbed optimization in Banach spaces I: A general theory based on a weak directional constraint qualification, II: A theory based on a strong directional qualification, III: Semi-infinite optimization[END_REF][START_REF] Bonnans | Perturbation analysis of optimization problems[END_REF], we are able to perform a sensitivity analysis along a path of the form

f (σ) := σf 1 + 1 2 σ 2 f 2 + o(σ 2 ), (119)
where f 1 and f 2 are given in L 2 (Q). Consider the following perturbed linearized equation

z t -∆z + 3γ ȳ2 z = f 1 + i ω v in Q; z = 0 on Σ, z(•, 0) = 0. ( 120
)
The related linearized optimization problem is

Min (v,z)∈L 2 (Qω)×H 2,1 (Q) J ′ (ū, ȳ)(v, z); g ′ (ȳ(t))z(t) ≤ 0 over I(g(ȳ)); (120) holds.
(L f1 ) Note that S(L f1 ) coincides with the critical cone if f 1 = 0. We start with a technical lemma concerning the value of the linearized problem (L f1 ).

Lemma 20. We have that val(L f1 ) = Q p(x, t)f 1 (x, t)dxdt. In addition, the set S(L f1 ) is non empty.

Proof. Since (ū, ȳ) satisfies the qualification condition, and in view of the uniqueness of the Lagrange multiplier, we deduce from [13, Prop. 3.1] that (i) (L f1 ) has the same Lagrange multiplier μ as the original problem (P ), (ii) (v, z) ∈ S(L f1 ) iff (v, z) ∈ F (L f1 ) with strict complementarity with μ, i.e. g ′ (ȳ(t))z(t) = 0 over supp(μ)

and (iii) that the formula for val(L f1 ) holds. The nonemptiness of the solution set follows from the fact that, in view of Lemma 8, we can control the time derivative of the state constraint over a neigbourhood of the contact set, and therefore it is possible to set the first-order variation of the state constraint to zero over the contact set.

Next we need an estimate of stability of approximate solutions. The Lagrangian of the family of perturbed problems (by abuse of notation, denoted in the same way) is L(u, y, f, p, q, µ) := J(u, y) + Q p ∆yγy 3 + f + i ω uy t dxdt + T 0 g(y(t))dµ(t) + Ω q(x)(y(x, 0)y 0 )dx.

(122)

It is easily checked that the value of the multiplier for the initial condition is q = -p(0). Lemma 21. Let (ū, ȳ) be a qualified local solution of (P ) satisfying the quadratic growth condition [START_REF] Casas | Pontryagin's principle for state-constrained boundary control problems of semilinear parabolic equations[END_REF]. If ρ is small enough, then (i) when f 2 is small enough, problem (P f ) has a nonempty (necessarily bounded) set of solutions, (ii) if f k → 0 in L 2 (Q), and u k is a sequence of o( f k 2 2 ) solutions of problem (P f k ), we have that u kū 2 = O( f k 2 ).

Proof. Since the Implicit Function Theorem may be applied to the study of the state equation ( 115)-(117), denoting its solution by y f,u , we have that, in the vicinity of (f = 0, ū), the mapping (f, u) → y f,u is well-defined and of class C ∞ . In particular, when ρ and f 2 are small enough, y f,u is uniformly bounded in H 2,1 Σ (Q). Therefore we can pass to the limit in a minimizing sequence of problem (P f ), using standard compactness arguments in order to deal with the nonlinear term of the state equation. Point (i) follows. (ii) If the conclusion does not hold, for arbitrarily small ρ > 0 there exist sequences f k → 0 in L 2 (Q) and u k of o(σ 2 ) solutions of problem (P f k ) such that f k 2 = o( u kū 2 ). We may write σ k := u kū 2 , and u k = ū + σ k v k , with v k 2 = 1. Extracting if necessary a subsequence, assume that v k ⇀ v. From a first-order expansion of the cost function and constraints we easily obtain that v is a critical direction. Now L(u k , y k , f k , p, q, μ) ≤ J(u k , y k ) ≤ val(P

f k ) + o( f k 2 2 ). ( 123
)
From a second-order expansion of the Lagrangian, whose partial derivative w.r.t. (u, y) is zero at (ū, ȳ, f = 0), and since L(ū, ȳ, 0, p, q, μ) = J(ū, ȳ), denoting by z k the solution of (120) with v := v k , get ∆(v k , z k ) ≤ o(1), and hence, ∆(v, z) ≤ lim inf k ∆(v k , z k ) ≤ 0. Since (v, z) is a critical direction, it follows that (v, z) = 0, and hence ∆(v k , z k ) → ∆(v, z), which implies v k → v. This cannot be since v k is of unit norm and v = 0.

Consider the quadratic subproblem Min

(

  v,z)∈S(L f 1 ) ∆(v, z) + Q p(x, t)f 2 (x, t)dxdt (Q)Theorem 22. Let (ū, ȳ) be a qualified local solution of (P ) satisfying the quadratic growth condition[START_REF] Casas | Pontryagin's principle for state-constrained boundary control problems of semilinear parabolic equations[END_REF]. Then a) we have the following expansionv(f (σ)) = val(P ) + σ val(L f1 ) + 1 2 σ 2 val(Q) + o(σ).(124)b) In addition we have that if(u σ , y σ ) is a path of o(σ 2 ) solutions, then u σū 2 = O(σ), each weak limit-point in L 2 (Q ω ) is a strong limit-point, and is solution of problem (Q).If the latter has a unique solution v, then a path u σ of o(σ 2 ) solutions of (P f (σ) ) satisfiesu σ = ū + σv + o(σ).(125)Proof. a) Apply[START_REF] Bonnans | Perturbation analysis of optimization problems[END_REF] Thms 4.94]. Its hypothesis are (i) directional regularity, that follows from the qualification condition, (ii) existence of a o(σ 2 ) solution path for problem (P f (σ) ), which holds in view of lemma 21, (iii) the directional extended polyhedricity condition which is consequence of extended polyhedricity. The conclusion is, with numbering of equations of[START_REF] Bonnans | Perturbation analysis of optimization problems[END_REF], that the expansion (4.235), based on the expressions defined in (4.222)-(4.224), holds. Therefore, since (as easily checked) val(Q) is finite, (our equation) (124) holds. b) This is an obvious consequence of [18, Thm 4.95].

RR n 6784

parameters ν 1

Ten such that, setting Q -T := Ω × ([0, T ]\T ):

y t -∆y + γy 3 = i ω u a.e. (x, t) ∈ Q, (80) y(0, •) = y 0 , (81) y = 0 on Σ (82) -p 1,t -∆p 1 + 3γy 2 p 1 = yy d + (2∆y -6γy 3 + i ω u)µ 1 (t) on Q -T (83)

g(y)(τ en ) = 0, τ en ∈ T en , d dt g(u, y) = 0, on I b (87)

A solution of the alternative formulation is said to satisfy the additional conditions, if the conditions below hold:

Consider the following relations:

We have the following propositions.

Proposition 16. Let (u, y) ∈ F (P ) satisfy the qualification condition. Then (u, y) is a extremal iff it satisfies both the alternative formulation (80)-(90) and the additional conditions (91)-( 93), the relations between the adjoint states and multiplier being given by (94)-(96).

For the proof it suffices to substitute in the alternative formulation the expressions of (p 1 , µ 1 ) in ( 94)-( 96). Details are left to the reader. The derivation is analogous to the one in e.g. [START_REF] Bonnans | Stability and sensitivity analysis for optimal control problems with a first-order state constraint[END_REF]. The above formulation opens the way to the resolution of the optimality system by shooting methods. This extension is not obvious, however, in view of the ill-posedness of the backward heat equation.

Second-order optimality condition 4.1 Second-order necessary optimality condition

Since the seminal work by Kawasaki [START_REF] Kawasaki | An envelope-like effect of infinitely many inequality constraints on second order necessary conditions for minimization problems[END_REF][START_REF] Kawasaki | The upper and lower second order directional derivatives of a sup-type function[END_REF][START_REF] Kawasaki | Second order necessary optimality conditions for minimizing a sup-type function[END_REF], we know that second-order necessary optimality conditions involve the difference of th Hessian of Lagragian with a "sigma-term" taking into account the curvature of the convex set involved Corollary 18. Let (v, z) be solution of the linearized equation. Consider the linearized equation with r.h.s. in feedback form:

, and in addition

Proof. By lemma 8, where here

relation (104) follows from (46). In addition, the mapping ν

as being a composition of continuous mappings), and (v, z) = (v, z) when ν = 0. The conclusion follows.

Proof of theorem 17. We apply the extended polyhedricity theory of [18, Section 3.2.3]). Denote by C R (u, y) the set of radial critical direction, i.e., of critical directions (v, z) such that g(y) + σg ′ (y)z ∈ K for some σ > 0. It is known that ∆(v, z) ≥ 0, for all (v, z) ∈ C R (u, y), and since ∆(v, z) is a continuous function, this still holds over the closure CR (u, y). We will next check the extended polyhedricity assumption CR (u, y) = C(u, y) under which, by the above arguments, the conclusion holds.

Let us check this hypothesis by induction on N tot := N tc + N b , where N tc and N b are the number of touch points and boundary arcs, respectively. The result obviously holds when N tot = 0; assuming that it holds for N tot -1, let t 0 ∈ (0, T ) be such that g(y(t 0 )) < 0 and that I(g(y)) ∩ (t 0 , T ] is the last connected component of I(g(y)).

Give a critical direction (v, z), consider the restriction of (y, u, z, v) to [0, t 0 ]. Then (v, z) is critical over [0, t 0 ] in the sense that it satisfies the linearized state constraint over the active set while being complementary to the multiplier. Since I(g(y)) has N tot -1 connected components on [0, t 0 ], by our induction argument, for any ε > 0, there exists a "radial critical direction over (0, t 0 )", i.e., v ε ∈ L 2 (ω × (0, t 0 )) and z ε ∈ H 2,1 Σ (Ω × (0, t 0 )) such that, for some σ ε > 0, g(y) + σ ε g ′ (y)z ε ≤ 0 over [0, t 0 ], (v, z) satisfy the linearized state equation over (0, t 0 ), and

Let τ be the smallest point of the last connected component of I(g(y)); we can increase t 0 if necessary in order to be sure that κ(t) > 0 over [t 0 , τ ]. Integrating


Member, IEEE Gabriel Vasile
email: gabriel.vasile@univ-savoie.fr

Member, IEEE Emmanuel Trouvé
email: emmanuel.trouve@univ-savoie.fr

Ivan Petillot
email: ivan.petillot@univ-savoie.fr

Member, IEEE Philippe Bolon
email: philippe.bolon@univ-savoie.fr

Jean-Marie Nicolas
email: nicolas@tsi.enst.fr

Michel Gay
email: michel.gay@gipsa-lab.inpg.fr

Senior Member, IEEE Jocelyn Chanussot
email: jocelyn.chanussot@gipsa-lab.inpg.fr

Tania Landes
email: tania.landes@insa-strasbourg.fr

Pierre Grussenmeyer
email: pierre.grussenmeyer@insa-strasbourg.fr

Vasile Buzuloiu
email: buzuloiu@alpha.imag.pub.ro

Jean-Marie Nicolas

Irena Hajnsek
email: irena.hajnsek@dlr.de

Christian Andres
email: christian.andres@dlr.de

Martin Keller
email: martin.keller@dlr.de

Ralf Horn
email: ralf.horn@dlr.de

J.-M Nicolas

High resolution SAR interferometry: estimation of local frequencies in the context of Alpine glaciers

Keywords: Adaptive neighborhood, Experimental synthetic aperture radar (E-SAR), glacier monitoring, local frequencies, SAR interferometry (InSAR), TerraSAR-X

I. INTRODUCTION

W ITH an increased resolution, new spaceborne synthetic aperture radar (SAR) images will be an important source of information for the monitoring of glacier activity, by providing regular measurements such as surface topography, velocity fields or rocks, and crevasse cartography. After successful SAR experiments on Arctic and polar region glaciers [START_REF] Fischer | Observation of recent surges of Vatnajökull, Iceland, by means of ERS SAR interferometry[END_REF], the study of Alpine temperate glaciers proved to be more complex. Some experiments carried out with existing spaceborne interferometric SAR (InSAR) data have shown that only specific data such as European Remote Sensing satellite (ERS) 1/2 tandem couples provide information on welloriented glaciers at certain times of the year [START_REF] Mattar | Validation of alpine glacier velocity measurements using ERS tandem-mission SAR data[END_REF]. However, the characteristics of future spaceborne SAR missions such as TerraSAR-X (11-day repeat cycle, left/right looking, various incidence angles, up to 2-m resolution, dual polarization, etc.) reinforce the potential of interferometric measurements over Alpine glaciers.

The existing standard interferometric processing chains such as the differential interferometric automated process applied to survey of nature (DIAPASON) [START_REF] Massonnet | CNES general-purpose SAR correlator[END_REF] or the repeated orbit interferometry package (ROI-PAC) [START_REF] Rosen | Synthetic aperture radar interferometry[END_REF], [START_REF] Rosen | Updated repeat orbit interferometry package released[END_REF] provide a complete D-InSAR processing chain. These algorithms have already been tested and successfully applied in the context of SAR interferometry to measure accurate displacements (volcanology or seismology). These processors are built from standard wellestablished InSAR algorithms such as coregistration of singlelook complex (SLC) images, phase scaling, and unwrapping. However, this procedure is often error prone when applied to temperate glacier monitoring because of the complex glacier topography (Fig. 1) or the phase noise and discontinuities at glacier boundaries.

With decameter resolution (about 20 m for ERS SAR images), the amplitude is often difficult to use to extract precise information. Different methods have been developed to derive ice velocity measurement, either by "speckle tracking" [START_REF] Strozzi | Glacier motion estimation using SAR offset-tracking procedures[END_REF], [START_REF] Gray | Velocities and ice flux of the Filchner ice shelf and its tributaries determined from speckle tracking interferometry[END_REF] or by feature tracking [START_REF] De Lange | Improvement of satellite radar feature tracking for ice velocity derivation by spatial frequency filtering[END_REF]. These methods should benefit from the new generation of SAR satellites with meter resolution images and fully polarimetric data. This should allow to detect different features, to identify backscattering mechanisms, and to improve the extraction of the geometric deterministic component of the interferometric phase by adaptive techniques [START_REF] Vasile | General adaptive neighborhood technique for improving SAR interferometric coherence and phase estimation[END_REF], [START_REF] Vasile | Intensity-drivenadaptive-neighborhood technique for polarimetric and interferometric SAR parameters estimation[END_REF] or multiresolution approaches [START_REF] Davidson | Multiresolution phase unwrapping for SAR interferometry[END_REF], [START_REF] Lopez-Martinez | Modeling and reduction of SAR interferometric phase noise in the wavelet domain[END_REF].

InSAR measurements include topographic fringes (phase variations due to the relief), which are useful in estimating the surface elevation for deriving glacier volume change [START_REF] Surazakov | Estimating volume change of mountain glaciers using SRTM and map-based topographic data[END_REF], and displacement fringes (induced by the glacier displacement between the two InSAR acquisitions), which are useful in estimating ice stain rates [START_REF] Forster | Measurement of glacier geophysical properties from InSAR wrapped phase[END_REF]. The main limitations are • temporal decorrelation due to glacier surface and subsurface fast changes, particularly during the warmest season [START_REF] Trouvé | Combining airborne photographs and spaceborne SAR data to monitor temperate glaciers: Potentials and limits[END_REF]; • volume decorrelation caused by microwave penetration into uniform ice and associated elevation bias for topographic measurements [START_REF] Dall | InSAR elevation bias caused by penetration into uniform volumes[END_REF]; • the need to separate the topographic fringes from the displacement fringes and possible atmospheric perturbations [START_REF] Massonnet | Radar interferometry and its application to changes in the Earth's surface[END_REF].

To estimate the surface displacement, one of the safest options consists in using an external digital terrain model (DTM) of the target area. However, the required accuracy of the DTM increases with the baseline between the two SAR acquisitions, making interferograms more and more sensitive to the topography. This is particularly crucial in the regions presenting deeply embedded valleys such as the glacial valleys in the Alps. In this context, reliable estimates of the phase gradient within the interferogram become more and more necessary as the resolution of SAR data increases: smaller ground resolution cells will reveal the local topographic variations (microrelief), and larger baselines (smaller altitude of ambiguity) will transform them into local fringe pattern. These fringes are difficult to detect without using locally adaptive techniques. Previous approaches based on the estimation of local frequencies proved to be rather efficient with low-resolution (LR) data [START_REF] Davidson | Multiresolution phase unwrapping for SAR interferometry[END_REF], [START_REF] Trouvé | Fringe detection in noisy complex interferograms[END_REF]. If correctly estimated, local frequencies can be used in the InSAR processing at different stages, such as interferogram filtering [START_REF] Lee | Intensity and phase statistics of multilook polarimetric and interferometric SAR imagery[END_REF] or phase unwrapping [START_REF] Ghiglia | Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods[END_REF]- [START_REF] Trouvé | Improving phase unwrapping techniques by the use of local frequency estimates[END_REF]. They are also used for accurate coherence estimation obtained by compensating the local phase slope to avoid the bias (underestimation) due to the phase rotation within the estimation window. A recent method proposed by Zebker and Chen consists, for instance, in dividing the interferogram into 8 × 8 or 16 × 16 pixel windows and measuring the local dominant fringe patterns in the Fourier domain [START_REF] Zebker | Accurate estimation of correlation in InSAR observations[END_REF].

In this paper, we propose to improve these approaches and to adapt them to the new context of high-resolution (HR) interferograms. An original two-step method for estimating local frequencies is presented. The first step estimates the main fringe pattern on large boxcar windows by using existing methods. It provides a robust estimation of the average local frequency corresponding to the so-called "LR" fringe pattern. The second step consists in removing this main fringe pattern and estimating the remaining local fringe patterns due to local variations, which are missed in the first step for different reasons such as

• difference between the fringe model estimated in the first step and the real fringes; • local aliasing when the LR fringes are already high frequency and local variations increase the fringe rate; • discontinuities of the fringe pattern within the LR estimation boxcar window.

To estimate the local frequency corresponding to the so-called "HR" fringe pattern, we propose to use intensity-driven adaptive neighborhoods (IDANs) that are more likely to follow the local fringe variations and a new method to estimate the 2-D local frequencies on such neighborhoods with variable shape and size. The key issue when using the algorithm proposed in [START_REF] Vasile | Intensity-drivenadaptive-neighborhood technique for polarimetric and interferometric SAR parameters estimation[END_REF] lies in the fact that surface orientation is highly correlated with SAR intensity values. This observation is valid mainly for glacial areas, where the presence of other factors influencing the SAR intensity signal is highly reduced. Results are presented using two different data sets over the Mer-deglace-the Leschaux and the Tacul glaciers, which are located in the Mont-Blanc area: a simulated interferogram obtained from an HR DTM computed from aerial photographs, and a real interferogram from an airborne Experimental SAR (E-SAR) campaign that took place in this area in October 2006. The remainder of this paper is organized as follows. Section II is dedicated to the presentation of the local frequency estimation algorithm. In Section III, the results obtained using the proposed approach are presented and compared to those given by previous methods developed in the framework of ERS interferograms. Both simulated and real interferograms are used. Results of phase unwrapping and interferogram filtering are also given, starting from the obtained frequency estimates. In Section IV, some conclusions and perspectives are presented. Eventually, Appendix briefly describes the principles of the 2-D local frequency estimation technique [START_REF] Trouvé | Fringe detection in noisy complex interferograms[END_REF] used with boxcar windows for the LR fringe pattern analysis.

II. HIGH-RESOLUTION FRINGE PATTERN ANALYSIS

Under the assumption of a large number of scatterers inside the resolution cell, the statistics of the interferometric phase have been completely characterized [START_REF] Lee | Intensity and phase statistics of multilook polarimetric and interferometric SAR imagery[END_REF]. Several interferogram filtering strategies have been reported [START_REF] Monti Guarnieri | SAR interferometry: A 'quick and dirty' coherence estimator for data browsing[END_REF]- [START_REF] Wu | A locally adaptive filter of interferometric phase images[END_REF].

Lee et al. proposed a spatially adaptive filtering method for improving the accuracy of the interferometric phase estimation [START_REF] Lee | A new technique for noise filtering of SAR interferometric phase images[END_REF]. The basic idea of the filter is to preserve edges and to better filter noise in their vicinity. A set of 16 directional subwindows are defined to locate the optimal estimation samples inside the considered neighborhood. The subwindow selection procedure is performed either after locally applying the phase unwrapping procedure or operating in the complex plane. Only the pixels in the selected subwindow are used to compute the filtered value, which is derived from the locally linear minimum meansquared error estimator. In [START_REF] Lopez-Martinez | Modeling and reduction of SAR interferometric phase noise in the wavelet domain[END_REF], a multiscale interferogram analysis is proposed. The phase image is locally analyzed in the Shannon wavelet domain, and the filter is characterized by good spatial resolution maintenance properties.

This paper presents a new method to obtain the local frequencies within an interferogram by fusing multiscale analysis technique [START_REF] Davidson | Multiresolution phase unwrapping for SAR interferometry[END_REF], multivariate adaptive neighborhood (AN) estimation algorithm [START_REF] Vasile | Intensity-drivenadaptive-neighborhood technique for polarimetric and interferometric SAR parameters estimation[END_REF], and robust LR squared fixed-size windowing procedure [START_REF] Trouvé | Fringe detection in noisy complex interferograms[END_REF]. This can only be achieved due to a new adaptive 2-D frequency estimator dedicated to neighborhoods with variable shape and size.

A. Proposed Model

For each pixel of the interferogram, the 2-D phase signal φ is considered as having two deterministic components as follows:

• the LR phase φ LR corresponding to the average fringe pattern created by the LR features; • the HR phase φ HR corresponding to the local difference of the phase signal with the average fringe pattern. This component is induced by HR features.

Furthermore, one random component corresponding to the phase noise φ n is added. Consequently, the analytical phase signal can be written as e jφ = e jφ LR e jφ HR e jφ n .

(

LR and HR components do not correspond to a partition of the frequency domain such as low or high frequencies, but to the spatial resolution of the associated fringe pattern. In the case of topographic fringes for instance, the LR phase corresponds to the mean slope of the area covered by a large estimation window [Fig. 2(a)], whereas the HR phase corresponds to smaller scale additional variations around this average slope due to the microrelief [Fig. 2(b)]. The first step of the proposed algorithm consists in estimating the LR phase and removing the associated fringe pattern. This is performed by estimating the 2-D LR local frequency with a modified version of the MUSIC algorithm [START_REF] Trouvé | Fringe detection in noisy complex interferograms[END_REF]. The method consists in building a 1-D signal vector by appending the lines of a small rectangular sliding subwindow. Then, the covariance matrix of this vector is estimated within a larger window. The estimation of the two frequency components is based on the structure of the vector covariance matrix (VCM) described in the Appendix. The principle of this VCM method requires large rectangular estimation windows, and it cannot be applied to irregular ANs.

B. Spatial Support

In the second step, the local HR frequencies are estimated. The proposed procedure consists in gathering only the pixels that belong to the same HR spatial feature and fulfill the "local stationarity" hypothesis [START_REF] Touzi | A review of speckle filtering in the context of estimation theory[END_REF]. We assume that small homogeneous regions within the two intensity images correspond to ground areas with a homogeneous cover and a constant orientation which should respect the stationarity hypothesis. Although not generally valid for any type of terrain (particularly flat vegetation or urban areas), this assumption is verified for rather steep-sloped temperate glaciers that are encountered over 1500 m above sea level (ASL) in the Alps. As a consequence, the information that is provided by the intensity multivariate vector is suitable to build IDAN [START_REF] Vasile | Intensity-drivenadaptive-neighborhood technique for polarimetric and interferometric SAR parameters estimation[END_REF]. Even if these neighborhoods are not exactly "fringe pattern driven," they are more likely to respect the stationarity hypothesis than blind boxcar windows according to the following considerations.

• Within a full SAR scene, the local variation of the backscattered radar radiation strongly depends on the local surface orientation of the targeted area. In software such as ROI-PAC for instance, the knowledge of the slope orientation provided by a DTM is sufficient to simulate SAR image and to register them on real ERS or ENVISAT images.

If the LR SAR data are mainly sensitive at the mean slope, the new HR SAR systems will more accurately reflect the small slope variations induced by the local microrelief.

• The local variations of the backscattered radar radiation imposed by local changes of the dielectric and magnetic permittivities may not be directly linked to phase nonstationarity. However, they are potential sources of different backscattering mechanism or temporal evolutions yielding to different coherence level in large baseline or repeatpass interferometry. For instance, it allows discriminating between rocks and ice or different snows in the particular context of Alpine glaciers.

The main advantage of the IDAN method is to gather a significant number of samples in an estimation window where stationarity is simultaneously preserved with respect to the two intensity images.

The two master I 1 and slave I 2 intensity images are considered as a bivariate image as follows:

p(m, n) = [I 1 (m, n), I 2 (m, n)] T .
(

) 2
The proposed AN determination relies on the image p(m, n) and involves processing vectorial images mixed with multiplicative noise model. The algorithm consists of the following steps [START_REF] Vasile | Intensity-drivenadaptive-neighborhood technique for polarimetric and interferometric SAR parameters estimation[END_REF]. 1) Rough estimation of the seed value: Compute the marginal median p(m, n) in the 3 × 3 centered neighborhood. 2) Region growing: The eight direct neighbors p(m , n ) of the seed are accepted inside the AN provided they meet the following condition:

p(m , n ) -p(m, n) p(m, n) ≤ T 1 , T 1 = 2 3 σ n µ n . (3)
Threshold T 1 is set according to the speckle mean µ n and standard deviation σ n , which are a priori known constants, as they both depend on the initial number of looks L and are identical for both intensity images. Then, the same procedure is applied for all of the neighbors of the newly included pixels and so on. The region growing is iterated until either the number of pixels already included in the AN exceeds a predefined upper limit N max or none of the new neighbors fulfills the test condition (3). The pixels that have already been tested but not accepted inside the AN (called background pixels in the following) are stored in a separate list. 3) Refined estimation of the seed value: A more reliable estimator of the unspeckled seed value p(m, n) is now obtained by averaging the pixels included in the AN defined in the previous step. 4) Reinspection of the background pixels: The background pixels p(o, r) of the list created in step 2) are tested again and aggregated in the AN provided that

p(o, r) -p(m, n) p(m, n) ≤ T 2 , T 2 = 2T 1 . ( 4
)
The AN is updated accordingly. The test is less restrictive, as the inclusion threshold is twice as large as the one used in the first step of the region growing procedure.

C. Autocorrelation Function

One major advantage of the AN lies in the fact that its shape adapts to the useful signal. This allows the shape of the spatial support to be irregular. As a consequence, conventional estimators of the local frequency, such as VCM [START_REF] Trouvé | Fringe detection in noisy complex interferograms[END_REF] or Fourier transform [START_REF] Zebker | Accurate estimation of correlation in InSAR observations[END_REF], cannot be used with AN.

This section introduces a new adaptive 2-D frequency estimator that can be successfully applied in the specific case of ANs. The main idea is to use the previously obtained AN as spatial support for computing the 2-D autocorrelation function of the phase signal.

Using 2-D notation, the phase signal is modeled by a firstorder approximation corresponding to a 2-D complex sine wave with frequency (f x , f y ) and additive noise φ n as

s(k, l) = e jφ(k,l) = e j[2π(kf x +lf y )+φ n ] .
(

) 5
The autocorrelation coefficients are defined by

γ(p, q) = E s(k, l) • s(k -p, l -q) (6)
where E{• • •} denotes the mathematical expectation and s(k-p, l-q) denotes the complex conjugate of s(k-p, l-q). Under the assumption that the phase noise is an additive independent identically distributed random process, we obtain γ(p, q) = K e j 2π( pf x +qf y ) [START_REF] Gray | Velocities and ice flux of the Filchner ice shelf and its tributaries determined from speckle tracking interferometry[END_REF] where K = | e j φ n (0) | 2 is a real coefficient that does not affect the autocorrelation phase anymore.

As the stationarity and ergodicity conditions are assured, the autocorrelation function is computed by spatial averaging of the products s(k, l) • s(kp, lq) for all the possible delays (p, q) within the AN. One can notice that the reliability of the γ(p, q) estimation strongly depends on the available number N p,q of pairs of pixels delayed with (p, q). Then, the local frequencies can be estimated from two consecutive values of the autocorrelation function (on lines or on rows). The proposed algorithm estimates the components of the 2-D frequency as the arguments of the following sums:

f x = 1 2π • arg   (p,q) N p,q N p+1,q • γ(p + 1, q)γ(p, q)   (8)
f y = 1 2π • arg   (p,q) N p,q N p,q+1 • γ(p, q + 1)γ(p, q)   . ( 9
)
By weighting the estimated value γ(p, q) by N p,q , the more reliable samples become more important. Moreover, for each component f x and f y , a confidence measure C f x (respectively C f y ) is defined by the ratio of the magnitude of the complex sum over the sum of the magnitudes as

C f x = (p,q) N p,q N p+1,q • γ(p + 1, q)γ(p, q) (p,q) N p,q N p+1,q • γ(p + 1, q)γ(p, q) . ( 10
)
When all the different products γ(p + 1, q)γ(p, q) have the same argument, the confidence is equal to 1 and reveals a good agreement between the fringe pattern and the complex sine-wave model. On the contrary, the confidence decreases when the products γ(p + 1, q)γ(p, q) have different arguments, due to unreliable estimation of the autocorrelation function or disagreement with the model. The global confidence

C f = min(C f x , C f y )
is taken as the minimum of the two confidences to reveal an unreliable 2-D frequency estimate as soon as one of the two components is not reliable.

D. Application to SAR Interferograms

The block diagram of the proposed two-step algorithm for local frequency estimation is presented in Fig. 3. First, an LR phase flattening is performed by using the 2-D LR frequency (f LR x , f LR y ) estimated by the VCM method [START_REF] Trouvé | Fringe detection in noisy complex interferograms[END_REF]. In each pixel (m, n) of the interferogram, the global fringe pattern is approximately removed in the surrounding pixels (k, l) by computing the following flattened phase φ fl : Finally, the obtained local frequencies can be directly used as input for either SAR data filtering [START_REF] Vasile | General adaptive neighborhood technique for improving SAR interferometric coherence and phase estimation[END_REF] or phase unwrapping [START_REF] Ghiglia | Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods[END_REF].

φ fl (k, l) = φ(k, l) -2π (k -m)f LR x + (l -n)f LR y . ( 11
)

III. RESULTS AND DISCUSSION

The proposed method has been tested on the "Chamonix-Mont-Blanc" test site which is located in the Alps, near the borders between France, Italy, and Switzerland (45 • 50 N, 6 • 51 E). It includes the Aiguille-Verte (4122 m ASL), the Chamonix valley (1000 m ASL), and several instrumented glaciers. The results presented in this section have been obtained on simulated TerraSAR-X HR data and on airborne E-SAR data over a group of three glaciers: 1) the well-known Mer-de-glace glacier and in its upper part; 2) the Leschaux glacier; and 3) the Tacul glacier.

A. Simulated TerraSAR-X InSAR Data

To illustrate the proposed method on HR data over Alpine glaciers, a 2-m resolution DTM of the Mer-de-glace [START_REF] Trouvé | Combining airborne photographs and spaceborne SAR data to monitor temperate glaciers: Potentials and limits[END_REF] has been used to simulate the SAR amplitude of TerraSAR-X images using the method proposed in [START_REF] Julea | Simulation of SAR images and radar coding of georeferenced information for temperate glacier monitoring[END_REF]. The parameters used in the simulation are presented in Table I.

Fig. 4(a) and (b) presents a real ERS amplitude (∼20-m resolution) and the simulated TerraSAR-X amplitude (∼2-m resolution), respectively. The area corresponds to a part of the Mer-de-glace glacier illustrated in Fig. 5. The altitude information provided by the DTM has been resampled in a TerraSAR-X slant range geometry [Fig. 4(c)] and converted into topographic fringes modulo a tunable altitude of ambiguity [Fig. 4(d)]. A controlled noise level was also introduced on the resulting HR phase and amplitude images. The comparison with LR image, such as the ERS acquisition illustrated in Fig. 4(a), shows that small relief variations become visible in the HR images, and the presence of crevasses can be noticed.

Fig. 6 illustrates the fringe orientation maps computed from the frequency estimates over the Mer-de-glace glacier and the Les-Echelets flank. A colored lookup table, which is provided in Fig. 6 The use of simulated HR interferogram allows objective performance assessments. A comparison between a conventional frequency estimation on a boxcar window and the multiresolution LR + HR proposed method has been performed by measuring the accuracy of the 2-D frequency estimates. The root-mean-square error (RMSE) has been computed using, as reference, the phase gradient measured by Prewitt kernels in the SAR-coded elevation map [Fig. 4(c)]. The results reported in Table II have been obtained by • single-resolution approach using the VCM method on conventional boxcar windows with a 7 × 7 size which is a good compromise between robust estimation and spatial resolution preservation; • LR + HR frequency estimation, LR by VCM method on large boxcar sliding window (11 × 11) and HR by the new algorithm on AN (50 samples).

The RMSE of the proposed method is significantly lower (almost divided by 2) than the RMSE of conventional boxcar estimation.

The importance of the HR frequency component can be assessed by directly injecting the 2-D local frequencies in least-squares phase unwrapping [START_REF] Ghiglia | Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods[END_REF], [START_REF] Trouvé | Improving phase unwrapping techniques by the use of local frequency estimates[END_REF]. The simulated interferogram has been unwrapped from the LR local frequency estimates and from the LR + HR estimates. After multiplying with the corresponding altitude of ambiguity and zero origin correction, the two results can be compared to the real altitude available from the HR DTM, which are resampled in the TerraSAR-X slant range. Fig. 7 presents a 50-m spatial profile along the surface of the Mer-de-glace glacier. The approximate position of the profile is marked in red in Fig. 5. Although the VCM LR unwrapped phase is rather linear, the proposed method follows the fractal shape induced by the local altitude variation.

B. Airborne InSAR Data

The proposed method has been applied on real HR InSAR data acquired by the E-SAR airborne system over the upper part of the Mer-de-glace glacier (the Tacul and Leschaux glaciers). A campaign with SAR image acquisitions and in situ measurements has been organized in October 2006, in collaboration between the German Aerospace Center (DLR) and the French laboratories of the MEGATOR project [START_REF] Trouvé | Combining airborne photographs and spaceborne SAR data to monitor temperate glaciers: Potentials and limits[END_REF]. The parameters of the data set used to illustrate the performances of the proposed method are reported in Table III.

The amplitude, coherence, and phase images of a 10-m baseline L-band interferogram are illustrated in Fig. 8. On the amplitude and the coherence images, one can observe the curved stripes that correspond to the Forbes band phenomena. This periodical feature indicates the one-year displacement of the Tacul glacier, and it is explained by the fact that during the summer season, a large amount of rocks fall on the ice in the upper part of the glacier and almost none in winter when rocks are much more stable. This phenomenon produces the smooth "staircase" shape in the upper part of the Tacul glacier. On the interferogram [Fig. 8(c)], the fringe pattern has two main components: 1) the topography and 2) the range evolution of the path difference (the so-called orbital fringes in the case of spaceborne data).

The results of the 2-D local frequency estimation performed at LR and the final combined LR and HR estimates are illustrated in Fig. 9 by the orientation of the corresponding phase gradient. For comparison, the orientation of the local fringe pattern computed using VCM on nonadaptive 7 × 7 squared centered neighborhood is presented in Fig. 9(b). As expected, the LR result in Fig. 9(a) is smooth and robust but it does not follow the local variations of the fringe pattern neither along the glacier surface nor on the surrounding mountains. This missing information is present in the sum of the LR and the HR components, which provides a more robust and more precise estimate of the local fringe pattern. One can notice a clear resolution improvement over the Forbes bands, where the local small slope variation is completely lost in the LR estimate and is hardly visible in the 7 × 7 VCM estimate. The local slope orientation along the Western side of the Aiguille-du-Tacul mountain (left side of Fig. 9) is more accurately estimated with the proposed algorithm. Fig. 10(b) presents the coherence filtering results obtained with the IDAN filter [START_REF] Vasile | General adaptive neighborhood technique for improving SAR interferometric coherence and phase estimation[END_REF], using the deterministic fringe pattern compensation provided by the HR and LR local frequency estimation [Fig. 9(c)]. The confidence image in Fig. 10(c) was thresholded to avoid using unreliable estimation of the HR fringe pattern. The value empirically chosen for the threshold is 0.4. In addition, the coherence map filtered with the 7 × 7 boxcar filter and the VCM 7 × 7 estimated local frequencies [Fig. 9(b)] are illustrated in Fig. 10(a).

Several different factors influence the coherence estimation, such as the true coherence value itself, the number of sample averaged to reduce estimation bias and variance, and the respect of the local stationarity assumption. To investigate the influence of the frequency estimates and the use of ANs on coherence Histograms in Fig. 11 show the results obtained on the E-SAR interferogram. The use of the AN instead of the boxcar for complex multilooking has a stronger effect on the coherence level than the use of an adaptive multiresolution frequency estimation technique instead of a nonadaptive one. In this particular case, the bulk of the improvement can be achieved by phase slope compensation with boxcar frequency estimation; hence, the coherence levels are affected only in small areas with local strong variations of the fringe pattern. However, it can be noticed that in both cases (boxcar and IDAN), the LR + HR frequencies have slightly improved the results compared to the VCM frequencies. The lower coherence obtained with the ANs is due to the larger number of samples averaged by IDAN with a stop criterion N max = 100 compared to the 49 samples averaged by the 7 × 7 boxcar. Lower coherence can be interpreted as a benefit of the proposed method (a reduction of the coherence overestimation due to a sufficient averaging)

or by a drawback if the lower values are due to the averaging of heterogeneous pixels in nonstationary areas. According to the previous results, regarding IDAN coherence estimation [START_REF] Vasile | Intensity-drivenadaptive-neighborhood technique for polarimetric and interferometric SAR parameters estimation[END_REF] and the visual improvement of the filtered coherence and phase images, it seems reasonable to interpret the coherence decrease by the "positive" explanation rather than by a lack of stationarity over the AN.

In conclusion, the different factors that affect the coherence and the lack of "ground truth" make it difficult to illustrate an improvement of the frequency estimation through the coherence level histograms. The main benefit of the proposed method, regarding coherence estimation, is the use of the same AN for both frequency estimation and coherence filtering to preserve the HR of the data.

IV. CONCLUSION AND PERSPECTIVES

The increased performances of future HR SAR systems (TerraSAR-X, COSMO-SKYMED, or Radarsat-2) together with the developments of appropriate processing algorithms reinforce the potential of SAR interferometry for temperate glacier periodic monitoring. A new method has been proposed to estimate the local frequencies within HR interferograms by a multiresolution algorithm: LR 2-D frequency are estimated on large boxcar windows by previous estimation methods such as VCM [START_REF] Trouvé | Fringe detection in noisy complex interferograms[END_REF] and HR 2-D fringe pattern remaining after LR phase flattening are analyzed by measuring their local frequency on IDANs. This method is adapted to deal with the small ground resolution cells and small altitudes of ambiguity, which makes the fringe pattern sensitive to the local topographic variations.

In the context of temperate glacier activity monitoring, it can be applied to improve the processing of HR interferograms to derive, for instance, very accurate DTM and differential DTM to obtain volume balances. In this case, HR frequency estimation combined with intensity-driven neighborhood techniques provides very good results since most of the sources of phase nonstationarity are revealed by the SAR intensities, which are mostly influenced by the local slope. It will also be useful for differential interferometry to obtain HR velocity field measurements that require efficient phase filtering and unwrapping methods. In this case, one difficulty that has to be addressed in future works is the presence of phase nonstationarity that is more related to the glacier displacement models and independent from the amplitude information. Future work should also investigate different strategies to perform HR fringe pattern analysis in different contexts, where the intensity information would not be dominated by the topography (urban areas, agricultural regions, etc.). A fusion of different information sources such as the coherence or the phase gradient would make the method applicable in more general terrain.

APPENDIX

The VCM method [START_REF] Trouvé | Fringe detection in noisy complex interferograms[END_REF] used to estimate the 2-D local frequencies (f x , f y ) of SAR interferograms on boxcar windows is based on the structure of a "signal vector" resulting from the concatenation of the lines of a 2-D complex sine-wave signal s(m, n) = e iφ(m,n) = e i2π(m•f x +n•f y ) . If this signal is observed on a N s × N s subwindow (square for the sake of simplicity), the size of the "signal vector" v s is N 2 s and the size of its covariance matrix Γ s is N 2 s × N 2 s . With this model, Γ s can be written as

Γ s = v s • v s † =                    1 e  2π•f x e  2π•2f x • • • e  2π(N s -1)f x --- e  2π•f y e  2π•(f x +f y ) • • • e  2π((N s -1)f x +f y ) --- • • •                    [vs † ]
.

Accordingly, this covariance matrix presents a N s × N s block structure: inside the blocks, lines are equal to the previous lines multiplied by a x = e i2π•f x and each block is equal to the previous block multiplied by a y = e i2π•f y . This block structure can be used to estimate the 2-D local frequencies of SAR interferograms by the following three-step algorithm:

1) estimate the N 2 s × N 2 s VCM under the stationarity and ergodicity hypothesis by complex averaging on a N e × N e sliding window, with N e > N s to average (N e -N s + 1) 2 matrix samples; 2) compute by least-squares estimation the coefficients a x and a y ; 3) derive the frequency estimates fx and fy from the arguments of these coefficients.

Fig. 1 .

 1 Fig. 1. May 2004: Surface topography and texture variability on the Mer-deglace glacier (human size marked in red).

Fig. 2 .

 2 Fig. 2. DTM profile on the Mer-de-glace glacier. (a) Using DTED-1 LR DTM (Res = 80 m). (b) Using HR DTM (Res = 2 m).

Fig. 3 .

 3 Fig. 3. Block diagram of the proposed method.

  The resulting HR phase signal e φ fl exhibits the local differences between the 2-D sine-wave model and the real fringe pattern. Second, the 2-D HR frequency (f HR x , f HR y ) of the HR phase signal is estimated from the 2-D autocorrelation function using IDANs as spatial support. According to the Nyquist-Shannon criterion, when the local variation ∆φ of the phase signal belongs to the [-π, π] interval, its local frequency f belongs to [-0.5, 0.5]. The proposed two-step method first estimates the 2-D LR spatial frequency (f LR x , f LR y ) in the domain [-0.5, 0.5] × [-0.5, 0.5]. After the removal of this average fringe pattern, the 2-D HR frequency (f HR x , f HR y ) can also be estimated in the domain [-0.5, 0.5] × [-0.5, 0.5]. In the case of local aliasing due to HR features, which makes the magnitude of a frequency component higher than 0.5 and lower than 1, the LR + HR model is able to recover this component and yields a total of 2-D frequency (f LR x + f HR x , f LR y + f HR y ), which belongs to the domain [-1, 1] × [-1, 1]. In other words, the proposed method is also able to retrieve phase gradients in limited aliased areas where ∆φ belongs to [-2π, 2π].

Fig. 4 .Fig. 5 .

 45 Fig. 4. Mer-de-glace glacier (840 × 340 pixels). (a) ERS-1 five-look amplitude. (b) Simulated one-look TerraSAR-X amplitude. (c) Elevation image in TerraSAR-X slant range geometry. (d) Simulated one-look TerraSAR-X phase.

Fig. 6 .

 6 Fig. 6. Orientation map of the local fringe pattern with TerraSAR-X (840 × 340 pixels). (a) LR with VCM 11 × 11 (local fringe orientation lookup table in the right upper corner). (b) VCM 7 × 7. (c) Sum of the LR and the HR components. (d) IDAN filtered phase with LR + HR phase flattening.

  Fig. 7. Spatial profile (50 m) along the surface of the Mer-de-glace glacier. (In green) Real altitude resampled in the TerraSAR-X slant range. (In red) Unwrapped HR + LR estimates of the local frequencies. (In blue) Unwrapped LR estimates of the local frequencies.

Fig. 8 .

 8 Fig. 8. E-SAR (16-look) interferogram (392 × 709 pixels). (a) Master amplitude. (b) Coherence map. (c) Interferometric phase.

Fig. 9 .

 9 Fig. 9. Orientation map of the local fringe pattern with E-SAR (392 × 709 pixels). (a) LR with VCM 11 × 11 (local fringe orientation lookup table in the right upper corner). (b) VCM 7 × 7. (c) Sum of the LR and the HR components.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. E-SAR L-band SAR data (392 × 709 pixels). (a) Boxcar 7 × 7 filtered coherence with VCM 7 × 7 phase flattening. (b) IDAN filtered coherence with LR + HR phase flattening. (c) Confidence map for HR frequency estimation.

Gabriel

  Vasile (S'06-M'07) received the B.Eng. degree in electrical engineering and computer science and the M.S. degree in image, shapes, and artificial intelligence from the University Politehnica of Bucharest, Bucharest, Romania, in 2003 and 2004, respectively, and the Ph.D. degree in signal and image processing from the Université de Savoie, Annecy-le-vieux Cedex, France. From 2001 to 2003, he was an Engineering Assistant with the Center of Microsystems and Microstructures for Environmental Monitoring and Biomedical Applications, National Institute for Research and Development in Micro and Nanotechnologies (IMT), Bucharest. From 2006 to 2007, he was a Development Engineer with the Laboratoire Grenoble Images Parole Signal Automatique, Centre National de la Recherche Scientifique, Institut National Polytechnique de Grenoble, France. Since 2007, he has held a Centre National d'Etudes Spatiales (French Space Agency) postdoctoral fellowship and works in the French Aerospace Laboratory (ONERA), Palaiseau, France. His research interests include signal and image processing, SAR remote sensing, polarimetry, interferometry, and glacier monitoring. Emmanuel Trouvé (M'99) received the B.Eng. degree in electrical engineering from the Ecole Nationale Superieure de Techniques Avancées, Paris, France, in 1990, the Ph.D. degree in signal and image processing from the Ecole Nationale Supérieure des Télécommunications, Paris, in 1996, and the "Habilitation à Diriger des Recherches" degree from the Université de Savoie, Annecy-levieux, France, in 2006. From 1996 to 1998, he was with Thomson Marconi Sonar in Underwater Acoustic and Signal Processing. Since September 1998, he has been an Associate Professor of signal processing and computer vision with the Université de Savoie, Polytech'Savoie, and works in the Laboratoire d'Informatique, Systèmes, Traitement de l'Information et de la Connaissance (LISTIC). He was the Coordinator of the French multilaboratory MEGATOR project (Monitoring the Evolution of Alpine Glaciers by Optical and Radar Remote Sensing, 2004-2007). His research interests include SAR image processing and data fusion in remote sensing. Dr. Trouvé is a member of the IEEE Geoscience and Remote Sensing Society Data Fusion Technical Committee. Ivan Petillot received the B.Eng. degree in electrical engineering and computer sciences and the M.S. degree in image, information, and hypermedia from the Institut de Recherche en Informatique de Toulouse (IRIT), Toulouse, France, in 2005. He is currently working toward the Ph.D. degree at the Université de Savoie, Annecy-le-vieux, France. His dissertation is on information fusion systems for remote sensing, in the framework of the MEGATOR project. He is with the Laboratoire d'Informatique, Systèmes, Traitement de l'Information et de la Connaissance, Université de Savoie, Polytech'Savoie, Annecy-le-vieux. His research interests include information fusion, risk assessment, and glacier monitoring.

TABLE I PROCESSING

 I PARAMETERS FOR TerraSAR-X IMAGE SIMULATION

TABLE III AIRBORNE

 III SAR DATA ACQUIRED OVER THE TEST AREA

ACKNOWLEDGMENT

The authors would like to thank the European Space Agency for providing ERS SAR data through the Category 1 proposal No. 3525 and the anonymous reviewers for providing useful comments and advice on this paper.

This work was supported by the French national project ACI-MEGATOR.

Philippe Bolon (M'94) was born in 1956. He received the B.Eng. degree in electrical engineering and the Ph.D. degree in signal processing from the National Polytechnic Institute of Grenoble (INPG), Grenoble, France, in 1978 and 1981, respectively.

From 1984 to 1994, he was an Associate Professor with the Université de Savoie, Annecyle-vieux, France. Since September 1994, he has been a Professor with the Université de Savoie, Polytech'Savoie, and since 2002, he has been the Head of the Laboratoire d'Informatique, Systèmes, Traitement de l'Information et de la Connaissance (LISTIC). Since 1996, he has been a member of the steering committee of the French national research group in Signal and Image Processing (GdR ISIS). His current research interests include information fusion, nonlinear image processing, image segmentation, and image analysis. He joined the German Aerospace Center (DLR), Microwaves and Radar Institute, SAR Technology Department, Wessling, in 1983. From 1983 to 1986, he was a Radar Engineer on the development of DLR's airborne experimental SAR system E-SAR, and later on, became a SAR Systems Engineer. In 1994, he was promoted as a Team Leader of the "Airborne SAR" group of the SAR Technology Department. His current responsibilities are mainly focused on the project management for the new DLR airborne SAR system F-SAR and on the management, coordination, and execution of DLR airborne SAR missions in Europe and abroad.