Benoit Michel

Nathalie Mazet
email: mazet@univ-perp.fr

Pierre Neveu

Experimental investigation of an innovative thermochemical process operating with a hydrate salt and moist air for thermal storage of solar energy: global performance

Keywords: Thermochemical process, open sorption process, solid-gas reaction, hydration, thermal storage, seasonal storage, solar energy

This paper investigates an innovative open thermochemical system dedicated to high density and long term (seasonal) storage purposes. It involves a hydrate/water reactive pair and operates with moist air. This work focuses on the design of and experimentation with a large scale prototype using SrBr 2 /H 2 O as a reactive pair (400 kg of hydrated salt, 105 kWh of storage capacity and a reactor energy density of 203 kWh/m 3 ). Promising conclusions have been obtained regarding the feasibility and performance of such a storage process. Hydration specific powers from 0.75 to 2 W/kg have been reached for a bed salt energy density of 388 kWh/m 3 . Moreover, two important parameters that control the storage system have been identified and investigated: the equilibrium drop and the mass flow rate of moist air. Both have a strong influence on the reaction kinetics and therefore on the reactor's thermal power.

Nomenclature:

A affinity, J•mol

Introduction

The energy demand for buildings accounts for 25% of total energy consumption in the world and 40% in Europe. Space heating represents the main part of this energy: 53% in the world and about 80% in Europe [1]. The increasing scarcity and cost of fossil fuels and incentives to reduce greenhouse gas emissions have led to a growing interest in solar energy. Solar energy is widely affordable and has the capability to meet household demand over the year. Unfortunately, its intermittency and variability with weather conditions, time, and seasons lead to a mismatch between heating demand and solar energy availability. Thus, a key issue to maximize the use of solar energy for space heating is to valorize the excess solar energy in summer by implementing long-term storage (3-6 months). Such a seasonal storage system must have the lowest heat losses between summer and winter, and the highest energy density to reduce its size and cost.

Several studies have dealt with seasonal storage for residential applications [START_REF] Pinel | A review of available methods for seasonal storage of solar thermal energy in residential applications[END_REF][START_REF] N'tsoukpoe | A review on long-term sorption solar energy storage[END_REF]. Among available storage systems, a sorption process takes advantage of a high storage density (about 100 to 500 kWh•m -3 of storage material), and weak heat losses between the storage and heating periods because the energy is stored as chemical potential and the sensible heat is weak. Therefore, such a system is relevant for seasonal storage for space heating for households. In comparison with competitive systems, the energy density is about 90 kWh•m -3 for latent storage and about 54 kWh•m -3 for water (sensible heat over ΔT=70 °C and 25% heat losses) [START_REF] Hadorn | Advanced storage concepts for active solar energy-IEA SHC Task 32 2003-2007[END_REF].

Sorption processes are increasingly being investigated for storage purposes, either regarding reactive materials or the workings of the storage system [START_REF] N'tsoukpoe | A review on long-term sorption solar energy storage[END_REF]. Several prototypes of sorption storage have been experimented with [START_REF] Stitou | Performance of a high temperature hydrate solid/gas sorption heat pump used as topping cycle for cascaded sorption chillers[END_REF][START_REF] Stitou | Experimental investigation of a solid/gas thermochemical storage process for solar air-conditioning[END_REF][START_REF] Mauran | Solar heating and cooling by a thermochemical process. First experiments of a prototype storing 60 kW h by a solid/gas reaction[END_REF][START_REF] Zondag | Prototype thermochemical heat storage with open reactor system[END_REF][START_REF] Weber | New Solar Thermal Storage Concept -Combined Hot Water and Sorption Store[END_REF][START_REF] Bales | TCA Evaluation -Lab Measurments, Modelling and System Simulations[END_REF][START_REF] Pierrès | CaCl2/H2O absorption seasonal storage of solar heat[END_REF][START_REF] Jaehning | Thermo-chemical storage for solar space heating in single-family house[END_REF][START_REF] Gartler | Development of a high energy density sorption storage system[END_REF][START_REF] Bales | Chemical and sorption heat storage[END_REF][START_REF] Marias | Thermochemical storage: first results of pilot storage with moist air[END_REF][START_REF] Iammak | Investigation of modular chemical energy storage performance[END_REF][START_REF] Boer | Solid sorption cooling with integrated storage: the SWEAT prototype. 3rd international heat powered cycles conference-HPC[END_REF][START_REF] Hongois | Development and characterisation of a new MgSO4-zeolite composite for long-term thermal energy storage[END_REF][START_REF] Bales | Laboratory Tests of Chemical Reactions and Prototype Sorption Storage Units. A Report of IEA Solar Heating and Cooling[END_REF] at a small scale; the largest is the SOLUX system, which stores 60 kWh [START_REF] Mauran | Solar heating and cooling by a thermochemical process. First experiments of a prototype storing 60 kW h by a solid/gas reaction[END_REF]. Nevertheless, no large-scale seasonal project has yet been completed [START_REF] Xu | A review of available technologies for seasonal thermal energy storage[END_REF].

Let's recall that different energy densities can be defined for these systems, referring to volumes of active material only, or to the whole reactor implementing the active material and including dead volumes and volumes dedicated to heat and mass transfer within the reactor and the coolant loops (heat exchanger, gas diffuser, collector…). Figure 1 compares these two energy densities for most of the prototypes, referring to the reactive material (x-axis) and the whole reactor prototype volume (y-axis). The energy density of the prototypes is significantly less than the reactive material energy density, i.e. only about half. Among all these experiments, the highest values of the material and prototype energy densities are, respectively, 280 kWh•m -3 for the reactive material and 140 kWh•m -3 for the storage system. The work described in this paper aims at developing a solid/gas sorption storage system beyond this range and experimenting with it on a significant scale (about 1/10 the capacity of a long term storage system for house heating).

This thermochemical storage process involves a reversible chemical reaction between a solid and a gas (equation 1). The reaction is a hydration/dehydration and a well-known hydrate/water pair has been selected (equation 2). The synthesis (or hydration) of the solid is exothermic (destorage step), while its decomposition (or dehydration) requires heat input (storage step). A thermochemical process based on such a hydrate/water pair can operate according to two different modes: the salt can react with pure water vapor at low pressure or with a moist air flow at atmospheric pressure. The first working mode is standard and simple, and most thermochemical systems in the literature operate with pure vapor [START_REF] Stitou | Experimental investigation of a solid/gas thermochemical storage process for solar air-conditioning[END_REF][START_REF] Mauran | Solar heating and cooling by a thermochemical process. First experiments of a prototype storing 60 kW h by a solid/gas reaction[END_REF][START_REF] Jaehning | Thermo-chemical storage for solar space heating in single-family house[END_REF][START_REF] Iammak | Investigation of modular chemical energy storage performance[END_REF][START_REF] Boer | Solid sorption cooling with integrated storage: the SWEAT prototype. 3rd international heat powered cycles conference-HPC[END_REF]. Nevertheless, this low pressure generates strong technological constraints for the reactor. On the other hand, the innovative operating mode with moist air allows a less expensive reactor conception [START_REF] Yu | Sorption thermal storage for solar energy[END_REF] but its management is more complex. Only a few authors are currently investigating the feasibility and performances of systems running with moist air [START_REF] Zondag | Prototype thermochemical heat storage with open reactor system[END_REF][START_REF] Marias | Thermochemical storage: first results of pilot storage with moist air[END_REF][START_REF] Michel | Thermochemical process for seasonal storage of solar energy: Characterization and modeling of a high density reactive bed[END_REF][START_REF] Hauer | Open absorption systems for air conditionning and themal energy storage, Thermal energy storage for sustainable energy consumption[END_REF][START_REF] Bertsch | Low temperature chemical heat storage-an investigation of hydration reactions[END_REF][START_REF] Abedin | Assessment of a closed thermochemical energy storage using energy and exergy methods[END_REF]. So far, the successful operation of a seasonal storage system based on a thermochemical process operating with moist air has never been demonstrated at the pilotscale. This challenge is the reason for this study.

Therefore, a prototype thermochemical storage reactor operating with moist air has been designed and experimented. Several previous works and literature reports outlined the background of this study. There are summarized in the following:

Beside the high storage density criteria, a seasonal storage system also has to fulfill requirements about thermal power production in order to meet user demand during the house heating (destorage) step. A previous simulation work using TRNSYS software and a large set of simulations and sensitivity studies reported and detailed in [START_REF] Tanguy | Seasonal storage coupled to solar combisystem: dynamic simulations for process dimensioning[END_REF] defined this demand, its requirements and target values for the storage system. The main results of this study were: a thermochemical seasonal storage system involving 5 to 10 tons of SrBr 2 ,6H 2 O (1300-2600 kWh of storage capacity) could allow a fractional energy savings (i.e. rate of energy consumption reduction compared to conventional non-solar installation) of about 65% to be reached for typical houses [START_REF] Heimrath | Project Report A2 of Subtask A, the Reference Heating System, the Template Solar System[END_REF] and French climates. Under such conditions this storage system should provide a thermal power from 3 to 4 kW. Thus, the required specific power at the reactor output was deduced: 0.3 to 0.8 W.kg -1 . This specific power range is the target of the prototype developed in this study.

The strontium bromide/H 2 O reactive pair has already been successfully experimented with in previous works [START_REF] Mauran | Solar heating and cooling by a thermochemical process. First experiments of a prototype storing 60 kW h by a solid/gas reaction[END_REF][START_REF] Mauran | 3rd Intermediate Report[END_REF]. Its ideal energy storage density is very high: 629 kWh•m -3 (referring to the bulk density, 2390 kg•m -3 , and the molar mass of non-porous hydrated salt, 0.3555 kg•mol -1 ) [START_REF] Lide | CRC Handbook of Chemistry and Physics[END_REF]. Regarding the energy density related to the reactive bed volume, a previous experiment on a small scale [START_REF] Michel | Thermochemical process for seasonal storage of solar energy: Characterization and modeling of a high density reactive bed[END_REF] demonstrated that promising values, higher than 400 kWh/m 3 , can be reached, highlighting its suitability as a long term thermal storage system.

In an open thermochemical system operating with moist air, the total gas flow rate passing through the bed is high and mass transfers have to be carefully investigated. The permeability of the strontium bromide porous bed was measured in a dedicated experimental set-up over a large range of packing densities, including the bed density of the current prototype [START_REF] Michel | Thermochemical process for seasonal storage of solar energy: Characterization and modeling of a high density reactive bed[END_REF].

Then, a 2D model of such a thermochemical system was developed [START_REF] Michel | Procédé thermochimique pour le stockage intersaisonnier de l'énergie solaire: modélisation multi-échelles et expérimentation d'un prototype sous air humide[END_REF]. Simulations demonstrated that in the range of packing densities corresponding to seasonal storage, the two operating modes (with pure vapor or moist air) lead to closely similar values of thermal power at the reactor output. For both modes, the mean specific power value over the reaction is higher than the target value (0.3 to 0.8 W/kg) as long as the bed energy density is less than 410 kWh/m 3 .

Thus, an experimental study of a thermochemical storage system operating with moist air was performed. It aims at:

-demonstrating the feasibility of such a storage system at a pilot-scale, -analyzing its capability to fulfill both energy density and specific output power criteria, and -determining the parameters allowing control of the output performances of the storage system.

Thermochemical system for thermal storage

2.1. Solid/gas thermochemical process

A thermochemical system is based on the thermal effect of a monovariant reversible reaction between a solid and a gas:

(1) For the prototype under study, the strontium bromide/water reactive pair was used. The reaction is:

<SrBr 2 , 1H 2 O> + 5(H 2 O) ↔ <SrBr 2 , 6H 2 O> + Δh 0 r (2)
<SrBr 2 ,1H 2 O> and <SrBr 2 ,6H 2 O> are respectively the dehydrated (S 0 ) and hydrated (S 1 ) salts, the reactive gas (G) is water and the stoichiometric coefficient ( ) is 5 mol G /mol s . The equilibrium conditions (p eqSG , T eqSG ) of such solid/gas reactions follow the Clausius-Clapeyron relation:

(3) p 0 is the reference pressure (1 bar). The standard enthalpy and entropy of reaction have been measured [START_REF] Mauran | Solar heating and cooling by a thermochemical process. First experiments of a prototype storing 60 kW h by a solid/gas reaction[END_REF] : =337000 J/mol s and =875 J/K/mol s (per mol of salt). Thus, the thermodynamic equilibrium conditions are determined by only one intensive variable: the gas pressure p eqSG or the temperature of the solid T eqSG . For the specific case of a thermochemical system operating with moist air at atmospheric pressure, the thermodynamic equilibrium conditions (Equation 2) are related to the partial pressure of water in moist air.

Working mode of a thermochemical storage system with moist air

The seasonal thermochemical storage process for space heating works as follows:

-Initially, at the end of summer, the reactive bed is ideally fully dehydrated. The salt is SrBr 2 ,1H 2 O -During the heating period (winter): a flow of moist air passes through the reactor, and the dehydrated salt reacts with water vapor. Exothermic synthesis (hydration) occurs and the storage system releases heat to the user ( per mol of salt). At the end of this step, the salt bed is fully hydrated (SrBr 2 ,6H 2 O).

-During the following storage period (summer): moist air is heated thanks to solar collectors, then this hot air flow passes through the porous bed of SrBr 2 ,6H 2 O. That leads to an endothermic decomposition reaction (dehydration) of the SrBr 2 ,6H 2 O salt to SrBr 2 ,1H 2 O.

-Between the storage and heating (destorage) periods, the reactor is closed and disconnected from the reactive gas flow, allowing the reaction heat to be stored over a long time. Heat losses are weak because they are linked to the sensible heat of the materials, which is significantly less than the heat of reaction. Thus, a thermochemical process is particularly suitable for seasonal heat storage.

Prototype design and experimental set up

Reactor design

The thermochemical reactor prototype was designed and sized in order to achieve both criteria that define a seasonal storage system as described above: a high energy density and a specific power suitable for space heating for residential housing.

Various solid/gas reactor configurations have been proposed in the literature [START_REF] Zondag | Comparison of reactor concepts for thermochemical storage of solar heat[END_REF][START_REF] Zondag | Engineering assessment of reactor designs for thermochemical storage of solar heat[END_REF] : screw reactor, fluidized reactor, fixed bed reactor…. Nevertheless, according to PROMES laboratory knowledge [START_REF] Stitou | Performance of a high temperature hydrate solid/gas sorption heat pump used as topping cycle for cascaded sorption chillers[END_REF][START_REF] Stitou | Experimental investigation of a solid/gas thermochemical storage process for solar air-conditioning[END_REF][START_REF] Mauran | Solar heating and cooling by a thermochemical process. First experiments of a prototype storing 60 kW h by a solid/gas reaction[END_REF][START_REF] Pierrès | Experimental results of a solar powered cooling system at low temperature[END_REF][START_REF] Mauran | Heat and mass transfer in consolidated reacting beds for thermochemical systems[END_REF], a simple and cheap reactor configuration involving a porous fixed bed of reactive solid was chosen for this experimentation.

The prototype was designed in a modular way by stacking eight rectangular modules filled with a fixed bed of reactive salt (Figure 2 and Figure 3). A ninth module was planned but finally not included due to manufacturing difficulties.

Both the top and bottom of each module are in contact with the moist air flow thanks to a 1.5 cm thick air blade which is used as a mass diffuser. Obviously, the diffuser volume diminishes the energy density of the reactor. In order to maximize the reactor energy density, the number of diffusers should be minimized by implementing thicker reactive beds. However, a high thickness can limit mass transfer and reduce the specific power delivered by the thermochemical reactor (the power decreases with the square of the bed thickness [START_REF] Michel | Thermochemical process for seasonal storage of solar energy: Characterization and modeling of a high density reactive bed[END_REF]). Thus, the thickness of the reactive bed is a key parameter in designing a thermochemical reactor and achieving a tradeoff between the two performance criteria.

A previous experimental study at small scale demonstrated that a reactive SrBr 2 bed featuring about 400 kWh/m 3 of energy density and a thickness of 5-10 cm reached the target value of specific thermal power in the winter phase [START_REF] Michel | Thermochemical process for seasonal storage of solar energy: Characterization and modeling of a high density reactive bed[END_REF]. Thus, the thermochemical reactor prototype was designed with a reactive bed 7.5 cm thick and 388 kWh/m 3 as the energy density. The other module dimensions (width: 65.1 cm, and depth: 69.4 cm, see Figure 2), were chosen in order to simplify the handling of the prototype. Each module contains 50 kg of hydrated salt; the whole prototype contains 400 kg of hydrated salt and can store 105 kWh.

The performance of this prototype has been assessed thanks to the simplified 1D model (sharp front model) developed and validated previously for similar reactive beds [START_REF] Michel | Thermochemical process for seasonal storage of solar energy: Characterization and modeling of a high density reactive bed[END_REF][START_REF] Lu | Comparison of a general model with a simplified approach for the transformation of solid-gas media used in chemical heat transformers[END_REF] and the permeabilities of the hydrated and dehydrated salt bed measured (k 1 and k 0 , respectively, see Table 1). The moist air flow rate through the bed is fixed at a constant value (300 m 3 /h). This set of operating and design conditions allows the storage reactor to supply a specific power matching the target range (see Table 1). Moreover, the chosen bed thickness leads to a pressure drop across the bed that can be technically managed easily (from 60 to 560 Pa, over the reaction). The whole design of the prototype is described in Figure 2 and Figure 3. On the bottom side of each module, a perforated plate and a metallic mesh avoid any salt losses, while the top side is free (Figure 2).

Operating and design conditions (input data)

On one side of each reactive bed, a gas diffuser (3 cm width) connected to the reactor air inlet supplies moist air to the reactive salt. This air flows axially through the reactive bed, then is collected on the opposite side by an outlet air collector (Figure 3). Finally, the prototype is insulated by a double layer of insulation: 5 cm of "Pyrogel XT 10", =0.022 W•m -1 .K -1 , and 3.2 cm of "Armaflex AC", =0.039 W•m -1 .K -1 . 2. Excluding feet, insulation and the ninth module, the reactor energy density is 203 kWh/m 3 . Note that this value can easily be enhanced in an industrial manufacturing process by optimizing the moist air flow through the whole prototype and reducing the volumes devoted to air diffusers and collectors. The prototype channels are oversized, as standard components were used for the sake of simplicity.

The air input and output of this reactor are connected to an aeraulic bench test providing a controlled moist air flow (temperature, moisture, flow rate).

Instrumentation and measurement uncertainties

The prototype was thoroughly instrumented at the following points in order to record the numerous variables and monitor the progress of the thermochemical reaction:

At the inlet and outlet air pipe: flow rate ( and , ±0.8 %), relative humidity (HR i and HR j , ±0.8 %), and total pressure (p ti and p tj , ±1 Pa).

Note that a chilled mirror hygrometer is used at the reactor outlet. It allows the moist air specific humidity (w j , ±0,20 K at the dew point measurement) to be measured with high accuracy, even when the value of the relative humidity is low, in the dehydration step.

At the reactor boundaries: the total pressure difference between the input and output of the reactor Δp ±20 Pa (differential pressure sensor), input/output moist air temperature T i , T j (PT100, ±0.2 K).

These measurements allow estimation of the reactor equivalent permeability, and the enthalpy balance of moist air flow between the reactor air input and output.

Within the reactor modules: 24 thermocouples are distributed in the salt beds (in modules 1, 3, 4, 5, 6 and 8), measuring axial and longitudinal temperature gradients across the salt layer. Several thermocouples located at the same position in all the modules permit temperature differences along the reactor axis to be evaluated. The temperature measurement uncertainties are ±0.3 K.

At ambiance: ambient relative humidity: HR amb ±0.8 %, ambient temperature: T amb ±0.3 K, ambient total pressure: p amb ±1 Pa.

Additionally, the reactor was placed on a large scale and weighted continuously to measure the reaction advancement defined by equation 3 during the reaction.

(4)

N s1,X the number of moles of hydrated salt in the reactor at the global advancement X, N s,t the total number of moles of salt. m s0 , m s1 and m sX refer to the total mass of the reactor, when the salt is fully hydrated, dehydrated and at advancement X, respectively. The weight measure uncertainties are ±0.2 kg. X is 0 for the dehydrated solid and X=1 at the end of the hydration reaction. From these measurements and using correlations based on the French standard NFX15-110 and the Ashrae Handbook Fundamentals [START_REF]Ashrae Handbook: Fundamentals[END_REF], it is possible to calculate all the characteristics of the moist air flow (temperature, humidity, enthalpy, density, ...). Measurements are recorded every 30 seconds. The results presented in the following sections are averaged values over 30 min of recorded data.

Experimental operating conditions

The rate of reversible solid gas reaction strongly depends on the equilibrium drop, i.e. the difference between the operating temperature and pressure and equilibrium thermodynamic values defined by equation 2 [START_REF] Lu | Modelling of gas-solid reaction-Coupling of heat and mass transfer with chemical reaction[END_REF][START_REF] Mazet | Improving the performance of the reactor under unfavourable operating conditions of low pressure[END_REF]. Thus, for the hydration and dehydration steps the prototype was tested at two equilibrium drops, a high and a low value. The two sets of operating conditions are presented in a Clausius-Clapeyron diagram in Figure 4 and in Table 3 (temperature, water partial pressure, the corresponding specific humidity and affinity). These equilibrium drops were chosen to be in the range of moist air conditions that can be encountered by a seasonal storage process: the high equilibrium drop corresponds to standard operation of the reactor, while the low equilibrium drop corresponds to more restrictive operating conditions. These equilibrium drops can be related to the affinity A by: [START_REF] Stitou | Performance of a high temperature hydrate solid/gas sorption heat pump used as topping cycle for cascaded sorption chillers[END_REF] with p v the partial pressure of water, and p eqSG the solid/gas equilibrium pressure at T. This affinity is negative for dehydration conditions and positive for hydration, and increases with the equilibrium drop.

Experimental results

The thermochemical storage prototype was tested for five months. Seven hydration/dehydration cycles were performed and various operating parameters were applied. Figure 5 shows the progress of the reaction (X) as a function of time for these cycles. As the study focuses on a seasonal storage system, the cycle times are rather long (15 to 20 days for the high equilibrium drop), so only a few cycles could be performed during this five-month experiment. Moreover, the cycles were not completed ( X<1) because the kinetics decrease drastically at the very end of the reaction, so completing these reactions fully would lead to a too long experiment time.

The experimentations were mostly performed with a high equilibrium drop and a constant moist air flow rate (313 m 3 /h during dehydration, 290 m 3 /h during hydration). In addition, the influence of these two parameters was investigated by experimenting with a variable flow rate (from 313 to 150 m 3 /h, cycles 6 and 7) and a variable equilibrium drop (alternately high and low value, cycle n°5). Intermittent operation, including several days stopped, was also experimented with (cycles 3 and 4) in order to analyze the off-season operation of a seasonal storage system. The following sections analyze these experimental results in order to better understand the operation of the moist air thermochemical storage system and assess its global performance. First, the performance of cycles operating under fixed conditions (a high equilibrium drop and constant flow rate) were analyzed ( § 4.1). Second, cycles run with variable operating conditions (flow rate, air temperature and humidity,...) were investigated in order to understand the dynamic influence of these parameters and ways to control the storage power and air temperature at the reactor outlet ( §4.2).

Global performance of the prototype

Reaction rate

Advancement as a function of time is shown in Figure 6a) and Figure 6b) for, respectively, hydration and dehydration, working at a high equilibrium drop. The operating conditions are listed in Table 4. First, Figure 6 shows that the kinetics is not constant over the reaction but decreases as the reaction progresses. This phenomenon is well-known for such solid/gas reactions. Second, the performance changes significantly between each of the first four cycles. Reaction times at mid-reaction (X = 0.5) range from approximately 60 h for the first hydration to about 100 h for the fourth hydration. Similarly, during dehydration, the reaction time at X=0.5 changes from 35 h for the first dehydration to 55 h for the third. Thus, the mid-reaction time increases 1.7 fold for hydration and 1.6 for dehydration compared to the first cycle. However, from the sixth cycle on, the kinetics remained within this range: the mid-reaction time is 70 h and 40 h for, respectively, hydration and dehydration. Therefore, despite the reaction rate evolution over the first cycles, the reaction kinetics seems to stabilize from the sixth cycle onward. This kinetics diminution followed by stabilization after a few cycles is a well-known effect for such solid/gas reactions [START_REF] Rambaud | Problématique des transferts en milieu poreux réactif déformable pour procédés de rafraîchissement solaire[END_REF]. It is usually explained by a change in the porous bed texture during the reaction due to the significant increase in the salt grain volume from the dehydrated to hydrated state, which causes strong mechanical stresses. Further experiments are in progress to measure this kinetics over a long period of time cycling at the prototype scale. This section deals with the thermal power released or consumed by the reaction (respectively, in hydration or dehydration). Thanks to the numerous sensors implemented in the prototype, it can be estimated in two ways. The first is based on the mass balance of the reactor and the second on the enthalpy balance of the moist air flow. Both methods are redundant as they lead to results within the same uncertainty range [START_REF] Michel | Procédé thermochimique pour le stockage intersaisonnier de l'énergie solaire: modélisation multi-échelles et expérimentation d'un prototype sous air humide[END_REF]. For the sake of simplicity, only the mass balance method is presented in this paper: The variation in the water contained in the salt over the reaction time is deduced from monitoring the mass of the reactor. The thermal power resulting from the reaction and the specific power can be expressed as: (6) [START_REF] Mauran | Solar heating and cooling by a thermochemical process. First experiments of a prototype storing 60 kW h by a solid/gas reaction[END_REF] With the salt mass variation during time, m s1 the mass of hydrated salt and M v the water molar mass.

Thermal power of the prototype

Figure 7 presents the thermal power due to the reaction, , as a function of the advancement for several cycles operating at a high equilibrium drop and constant flow rate (see Table 4). As discussed above, performances change significantly over the first four cycles. The reaction thermal power at X=0.5 is higher than 800 W for the first hydration and about 300 W for the fourth hydration. In dehydration, the reaction thermal power is higher than in hydration phase because of the larger equilibrium drop (Table 3). But the same power decrease occurs over the first cycles. At the sixth cycle, the thermal power increases and seems to become quite stable for hydration. At X=0.5 of the sixth cycle, the reaction power is -770 W in dehydration and 485 W in hydration. Further experimentations are in progress in order to validate the reaction rate stabilization over numerous cycles.

Regarding the specific power, the four first hydrations lead to values from 2 W/kg to 0.75 W/kg of hydrated salt (at X=0.5), and the sixth to 1.21 W/kg of hydrated salt. For dehydration, the specific power is higher: 1.9 W/kg of hydrated salt at X=0.5 for the sixth cycle.

Thus, these various experimentations lead to promising results as the specific power of the prototype is significantly higher than the target value to W/kg, see §1) required to fulfill the space heating demands of classic houses.

Effect of intermittent operation of the thermochemical reactor

Two steps of long-time stop (five and ten days) were introduced in the third hydration experiment. For both stop steps, the reactor temperature dropped naturally down to ambient temperature (20°C). These stops are similar to the operation of a storage system during offseason periods (fall and spring) when storage can be started and stopped several times for a few days each time.

Figure 8 presents local temperatures and thermal power evolutions. Figure 8a shows that, after each stop, the reaction power rises very quickly up to the same value as before the stop. The temperature behavior gauged by four thermocouples located in the middle of module 5 at different heights in the salt bed was also similar (Figure 8b). The temperature profile in the bed does not change before and after each stop step. Note that the decrease in the salt temperature and reaction power during hydration is not caused by the stop steps, but is due to the classic decrease in such solid/gas reaction kinetics that occurs with reaction progress (see §4.1.1). Moreover, Figure 7 shows that the reaction power of this third discontinuous hydration was similar to that in other continuous hydrations.

Therefore, long-term stop steps don't have any significant effect on the global performance of the thermochemical reactor.

System control

So far, this paper has fully demonstrated the feasibility of a thermochemical storage system and analyzed its working mode. In addition, it is important to control the power and temperature output supplied by the storage in order to meet user demand. Thus, this portion investigates the influence of two parameters that could allow control of such a thermochemical storage system: the moist air flow rate and the equilibrium drop.

Influence of the moist air flow rate

This experiment (cycle 7) involved four steps carried out at different mass flow rates of moist air: 313, 250, 200 and 150 m 3 /h, and a constant high equilibrium drop. Each step lasted seven hours, and the whole experiment ran from advancement X=0.55 to 0.27 in the dehydration phase, and from X=0.3 to 0.43 in the hydration phase.

Figure 9a and Figure 9b, present the reaction thermal power as a function of the moist air mass flow rate for hydration and dehydration. Each point is an averaged value calculated over the seven hour step at a fixed flow rate. These figures demonstrate the quite linear relationship between the thermal power of the reaction and the moist air flow rate for both reactions. As discussed in [START_REF] Michel | Procédé thermochimique pour le stockage intersaisonnier de l'énergie solaire: modélisation multi-échelles et expérimentation d'un prototype sous air humide[END_REF], this linear correlation also indicates that mass transfer in the porous reactive bed is the limiting phenomenon of the reaction in the prototype under study.

Thus, in this case, the thermal power consumed/supplied by the thermochemical reactor can be easily managed by controlling the flow of moist air through the reactive bed.

Influence of the equilibrium drop

The fifth cycle (hydration/dehydration) was carried out by alternating periods with high and low equilibrium drops (see Table 3). All other parameters were identical to previous cycles. Each period was long enough to reach quasi-stationary state. Each one lasted at least 20 h, and the whole hydration / dehydration cycle lasted respectively 160 h and 350 h.

Figure 10 presents the thermal power as a function of advancement over this unique cycle. This figure exhibits a strong change in the power according to the equilibrium drop. In the case of dehydration, at mid-hydration (X=0.5) the thermal power changed from -600 W to -180 W for an equilibrium drop which corresponded, respectively, to air input temperatures of 80°C and 60°C. The lower temperature leads to a very small reaction affinity (Table 3), which explains the weak power released by the reactor in this case.

During hydration, at X=0.5, the power changed from 500 to 250 W going from a high to a low equilibrium drop.

Notice that in the case of hydration the experimental protocol is slightly different as the high and low equilibrium drops correspond to two different water partial pressures, respectively 1000 Pa and 700 Pa (see Figure 4). Consequently, the input mass flow rate of water vapor also diminishes and results in a decrease in the reaction power. The effects of the equilibrium drop and water partial pressure overlap.

However, it is possible to quantify the portion resulting from each effect. Let's recall that in the previous paragraph Figure 9 correlates the reaction power and moist air mass flow rate, and allows the correlation between this power and the water vapor mass flow rate (linear relationship) to be deduced. Therefore, the power change at X=0.5 (figure 10) due to the decrease in the water mass flow rate can be evaluated as 148 W (at X=0.5) and, consequently, the change in the equilibrium deviation leads to a decrease of 102 W. Thus, the equilibrium drop has a strong influence: a rather slight variation in the operating conditions (of 300 Pa in hydration, of 20°C in dehydration, see Table 3) can change the reactor's thermal power a significant extent. It is therefore essential to control this parameter by acting on the humidity or temperature of the inlet moist air (see Figure 4). Selecting a high humidity moist air source for the hydration step and using a solar collector able to efficiently heat this moist air flow in the dehydration step are both key ways to control the input/output thermal power of the storage reactor. Moreover, this strong sensitivity of reactor performance to the equilibrium drop confirms the need for accurate models, as developed in a previous paper [START_REF] Michel | Procédé thermochimique pour le stockage intersaisonnier de l'énergie solaire: modélisation multi-échelles et expérimentation d'un prototype sous air humide[END_REF], describing the particular workings of the thermochemical reactor.

Finally, this study demonstrates that managing the equilibrium drop and the moist air flow is an efficient way to control the thermal power supplied by the thermochemical system.

Output air temperature

Classical thermochemical systems operating with pure vapor are generally limited by conductive heat transfer within the porous reactive bed [START_REF] Mauran | Heat and mass transfer in consolidated reacting beds for thermochemical systems[END_REF][START_REF] Lu | Modelling of gas-solid reaction-Coupling of heat and mass transfer with chemical reaction[END_REF][START_REF] Michel | Procédé thermochimique pour le stockage intersaisonnier de l'énergie solaire: modélisation multi-échelles et expérimentation d'un prototype sous air humide[END_REF]. Thus, they usually operate in an isobaric manner and their working conditions (temperature, T, and vapor pressure p v ) can be plotted on a horizontal isobaric line in a Clausius-Clapeyron diagram. Conversely, for a system working with moist air, the conditions at the input (T i , p vi ) and output (T j , p vj ) of the reactive bed are not located on an isobaric line (i.e. the water pressure evolves as the moist air flows across the reactive bed). Assuming stationary operation and neglecting heat losses, a simplified linear relationship between the input-output temperature lift T j -T i , and the outlet partial pressure of moist air, p vj , (or the outlet specific humidity w j ) has been developed by several authors [START_REF] Michel | Procédé thermochimique pour le stockage intersaisonnier de l'énergie solaire: modélisation multi-échelles et expérimentation d'un prototype sous air humide[END_REF][START_REF] Marias | Experimental results of thermochemical storage pilot reactor with moist air[END_REF][START_REF] Hauer | Development of a commercial thermal energy storage -From thermodynamics to the final product[END_REF]. This linear relationship depends only on the input humidity conditions (w i or p vi ) and characteristics of the solid/gas reactive pair.

T' is the temperature in C, and c ma and c mv are, respectively, the heat capacity of the dry air and the vapor. An example of this linear relationship is plotted as p vj T j in Figure 11 for the seventh hydration and dehydration, and for input operating conditions T i , p vi given in Table 4 (high equilibrium drop). Thus, for these input conditions (T i and p vi or w i ), the output moist air conditions T j and p vj are located on this 'temperature lift' line defined by Equation 8. Moreover, the reactive bed conditions T j , p vj are obviously limited by the thermodynamic equilibrium conditions defined by Equation 2. Therefore, during hydration, the maximum output temperature T j is the crosspoint of this temperature lift line (equation 8) and the equilibrium line (equation 2), both plotted in Figure 11a. Symmetrically, this crosspoint is the minimum output temperature during dehydration (Figure 11b).

Thus, for the 7th cycle, Figure 11a shows that the maximal moist air temperature at the reactor output during hydration is 35.4°C. For input conditions corresponding to the low equilibrium drop (not plotted), this maximum output temperature would be 31.8 °C.

Figure 11 also presents the experimental output temperature and pressure for each portion of the seventh cycle carried out at different mass flow rates during hydration and dehydration as described in part 4.2.2. These figures clearly show that there is no significant variation in the output moist conditions over the experimental range of flow rates (from 150 to 3 /h).

The slight displacement of the output point is due to the reaction's progress and the corresponding decrease in the reaction rate, as presented for hydration 3 in Figure 7. At the end of the reaction this output point will move quickly up to the input point. Thus, this experiment and analysis clearly demonstrate that controlling the air flow rate is not a relevant way to control the temperature at the reactor output.

Conclusion

This work investigates a thermochemical system dedicated to high density and long term (seasonal) storage purposes. This paper focuses on an innovative working mode with moist air. It leads to a simpler reactor design and a cost than the usual pure vapor working mode but the management of this original system has to be thoroughly investigated. Experiments were carried out on a prototype at a significant scale. They allow promising conclusions to be drawn regarding the feasibility, working mode and performance of such a storage process.

Thanks to previous works about reactive materials and solid/gas reactor modeling, a large prototype (400 kg of hydrated salt, 105 kWh storage capacity, and Der = 203 kWh/m 3 ) was designed, built and experimented with for several hydration/dehydration cycles. The experiments showed a decrease in the reaction rate for the very first cycles, then the reactor kinetics seemed to stabilize from the sixth cycle on. The specific thermal power released by the reactor during hydration was higher than the target range defined for space heating in typical houses in the French climate. An intermittent operating mode combining reactive and non-reactive steps, such as can occur during offseason periods, did not affect the temperature profile within the reactive bed nor the thermal power output released by the storage system.

Furthermore, ways to control this storage system in order to meet user demand were investigated. This control involves, on the one hand, the thermal power released by the reactor in the destorage step, and on the other hand the moist air temperature at the reactor outlet. Two important parameters were investigated: the equilibrium drop (i.e. the difference between moist air input conditions and thermodynamic equilibrium of the reaction) and the mass flow rate of the moist air. Both have a strong influence on the reactor kinetics and therefore on the reactor thermal power. A direct and quite linear relationship between the moist air mass flow rate and the reactor thermal power output has been highlighted. Conversely, the output moist air temperature is not changed by the air flow rate (in the range tested) but it depends on the input conditions of the moist air flow (temperature and humidity).

Therefore, the mass flow rate is a simple control parameter but it can only regulate the reaction power supplied by the thermochemical reactor. On the other hand, controlling the moist air conditions at the reactor inlet is relevant to regulating both reactor features: its thermal power and the moist air outlet temperature.

Thus, this study clearly demonstrates the feasibility of a long-term storage system involving a thermochemical process operating with moist air, and has enabled the identification of efficient and simple ways to control the thermal power output in order to meet user demand.

Figure 1 :

 1 Figure 1: Energy densities of various sorption storage prototypes as function of the energy density of the reactive materials implemented in the prototype: TCA (LiBr/H2O, [10]), HYDES and MODSTORE (Silica gel/H2O, [12,13]), SPF (Zeolithe 13X/H2O, [19]), MONOSORP (Zeolithe 4A/H2O, [14]), CETHIL (MgSO4+Zeolithe 13X/H2O, [18]), SOLUX (SrBr2/H2O, [7]), ECN (MgCl2/H2O, [8]), and ESSI: the prototype reported in this paper (SrBr2/H2O).

Figure 2 :Figure 3 :

 23 Figure 2: Module of the prototype a) Schematic overview of a module b) Module before filling with salt, with two thermal probes. c) Module filled with salt

Table 2 :

 2 *, h x L x l (cm) 99.3x77.5x72* Reactor external dimensions, h x L x l (cm) 130x106x90 Module dimensions, h' x L' x l' (cm) 8x69.4x65.1 Air collector dimensions, h x L x l (cm) 99.3x3x72 Mass diffuser dimensions, h x L x l (cm) 1.5x69.4x65.1 Bed thickness, Zs (cm) 7.5 Mass of hydrated salt, m 1 (kg) 400 Salt bed energy density, Dec (kWh•m -3 ) 388 Reactor energy density*, Der (kWh•m -3 ) 203* Geometric characteristics of the prototype. (* excluding feet, insulation and ninth module)

Figure 4 :

 4 Figure 4 : Experimental operating conditions for the reactor prototype based on the SrBr2,1-6H2O reaction, in a Clausius-Clapeyron diagram.

Figure 5 :

 5 Figure 5 : Advancement vs. time for seven cycles performed by the thermochemical storage reactor prototype.

  /dehydration cycle T i (°C) p vi (Pa) Air flow rate Advancement X

Figure 6 :

 6 Figure 6: Experimental advancement vs. time for reactions at a high equilibrium drop and constant air flow rate: a) hydrations; b) dehydrations. Most of the reactions were not complete and begin at different advancements. Thus, the initial time of each reaction has been shifted in order to compare the set of hydrations and dehydrations to the first cycle in the same advancement range.

Figure 7 :

 7 Figure 7 : Thermal power ( ) released or consumed by the reaction in the prototype as a function of advancement for different cycles operating at a high equilibrium drop. The target values of specific power during hydration (0.2 to 0.8 W/kg of hydrated salt, see §1) correspond to a reactor power ranging from 80 W to 320 W (green area).

Figure 8 :

 8 Figure 8 : Intermittent operation of the prototype including two long time stops (third hydration) a) Temperature of the salt bed (module 5), at different heights (1.8, 3.3, 4.9 and 6.7 cm), b) Thermal power released by this hydration.

Figure 9 :

 9 Figure 9 : Reaction thermal power as a function of the mass flow rate of moist air. a) hydration, b) dehydration (Cycle 7).

Figure 10 :

 10 Figure 10: Thermal power ( ) released or consumed by the reaction as a function of the advancement for a cycle alternating high and low equilibrium drops.

Figure 11 :

 11 Figure 11 : Average input and output air conditions (pv and T) for the 7th cycle. a) Hydration b) Dehydration Symbols refer to different moist air flow rates. Blue line: temperature lift line (equation 8). Red dotted line: thermodynamic equilibrium (equation 2).

Table 1 : Simulation of the working prototype using the sharp front model under winter (destorage) conditions: input data (reactive bed energy density, Dec, dehydrated and hydrated permeabilities, k0 and k1, inlet water partial pressure, pvi, moist air flow rate, , and reaction temperature, Tc) and resulting specific power ( ) and pressure drop across the bed (Δp).

 1

	Simulated performances

Table 3 : Sets of thermodynamic conditions imposed on the moist input air for prototype experimentation.

 3

			Hydration	Dehydration
			High	Low	High	Low
			equilibrium drop	equilibrium drop	equilibrium drop	equilibrium drop
	T i (°C)		25	25	82	62
	p vi (Pa)		1 000	700	2 500	2 500
	ω i (kg v .kg -1	a )	0,0062	0,0043	0,0157	0,0157
	Affinity (J•mol -1 s )	19 040	14 620	-28 200	-7 650

Table 4 : Operating conditions of each reaction carried out at a high equilibrium drop and constant air flow rate. The difference between final and initial advancement of successive reaction corresponds to experimentations carried out with different operating conditions.

 4

Acknowledgement

We thank the French National Research Agency (Agence Nationale de la Recherche) for their financial support for the ESSI project (ANR-08-STOCK-E-04) and the project partner INES developing the aeraulic bench.


Xin Su

Charles-Alban Deledalle

Florence Tupin

Hong Sun

NORCAMA: Change Analysis in SAR Time Series by Likelihood Ratio Change Matrix Clustering

Keywords: NORmalized Cut on chAnge criterion MAtrix (NORCAMA), Change Criterion Matrix (CCM), Change detection, Change classification, Likelihood Ratio Test (LRT), SAR time series

This paper presents a likelihood ratio test based method of change detection and classification for synthetic aperture radar (SAR) time series, namely NORmalized Cut on chAnge criterion MAtrix (NORCAMA). This method involves three steps: 1) multitemporal pre-denoising step over the whole image series to reduce the effect of the speckle noise; 2) likelihood ratio test based change criteria between two images using both the original noisy images and the denoised images; 3) change classification by a normalized cut based clustering-and-recognizing method on change criterion matrix (CCM). The experiments on both synthetic and real SAR image series show the effective performance of the proposed framework.

Introduction

Change analysis in remote sensing images is the process of analyzing differences (including identifying, recognizing and so on) in regions of interest by observing them at different dates [START_REF] Singh | Digital change detection techniques using remotely-sensed data[END_REF]. Many applications of remote sensing images involve change analysis, such as rapid mapping of disaster, landuse and land-cover monitoring and so on. [START_REF] Lu | Object-oriented change detection for landslide rapid mapping[END_REF] used change detection method to detect and locate the landslides for rapid mapping of landslides. Similarly, a multi-sensor change detection method between optical and synthetic aperture radar (SAR) imagery is proposed in [START_REF] Brunner | Earthquake damage assessment of buildings using VHR optical and SAR imagery[END_REF] for earthquake damage assessment of buildings. For urbanization monitoring, post-classification change detection methods are proposed (Taubenböck et al., 2012;Yin et al., 2011). From a methodological point of view, change analysis methods can be classified into two classes, binary-temporal change analysis and multi-temporal change analysis according to the number of images.

In the binary-temporal change analysis of two optical images, the most widely used operator is difference operator [START_REF] Singh | Digital change detection techniques using remotely-sensed data[END_REF]. For multi-spectral images, change vector analysis [START_REF] Bruzzone | Automatic analysis of the difference image for unsupervised change detection[END_REF] is proposed. People also perform the analysis on the transformed data instead of the spectral data directly, such as Tasseled Cap transformation [START_REF] Fung | An assessment of TM imagery for land-cover change detection[END_REF][START_REF] Huang | Derivation of a tasselled cap transformation based on Landsat 7 at-satellite reflectance[END_REF], principal component analysis [START_REF] Fung | Application of principal components analysis to change detection[END_REF][START_REF] Deng | PCA-based land-use change detection and analysis using multitemporal and multisensor satellite data[END_REF] and independent component analysis [START_REF] Marchesi | ICA and kernel ICA for change detection in multispectral remote sensing images[END_REF]. Beyond change detection, [START_REF] Bruzzone | An iterative technique for the detection of land-cover transitions in multitemporal remote-sensing images[END_REF] explicitly identified land-cover transitions (changes among Bare soil, Corn, Soybean, Sugar beet, Wheat) in multi-temporal remote-sensing images based on supervised classification. Given SAR images, two main approaches have been developed in the literature for change analysis: coherent change detection and incoherent change detection. The former uses the phase information in the SAR images through the study of the coherence map [START_REF] Preiss | Coherent change detection: theoretical description and experimental results[END_REF]. In incoherent change detection, the log-ratio [START_REF] Rignot | Change detection techniques for ERS-1 SAR data[END_REF] is the most common operator. Improvements have been proposed thanks to automatic thresholding methods [START_REF] Bazi | An unsupervised approach based on the generalized Gaussian model to automatic change detection in multitemporal SAR images[END_REF] or multi-scale analysis to preserve details [START_REF] Bovolo | A detail-preserving scale-driven approach to change detection in multitemporal SAR images[END_REF]. [START_REF] Lombardo | Maximum likelihood approach to the detection of changes between multitemporal SAR images[END_REF] proposed a generalized likelihood ratio test given by the ratio between geometric and arithmetic means for SAR images. [START_REF] Quin | MIMOSA: An Automatic Change Detection Method for SAR Time Series[END_REF] extended this ratio of different means to a more general way with an adaptive and nonlinear threshold, which can be applied to not only SAR image pairs but also SAR time series.

Beyond change analysis between two dates, multi-temporal change analysis (more than 2 dates) mainly focuses on the longterm change information. SAR image features consisting of long-term coherence and temporal backscattering is proposed for a classification purpose [START_REF] Bruzzone | An advanced system for the automatic classification of multitemporal SAR images[END_REF]. [START_REF] Julea | Unsupervised spatiotemporal mining of satellite image time series using grouped frequent sequential patterns[END_REF][START_REF] Julea | Efficient spatio-temporal mining of satellite image time series for agricultural monitoring[END_REF] propose a crop monitoring using satellite image time series by a frequent sequential pattern (a group of pixels sharing common temporal patterns and satisfying a minimum spatial connectivity). A generic change detection approach is proposed in (Verbesselt et al., 2010a) for multi-temporal images by detecting and characterizing breaks for additive seasonal and trend changes. It integrates the decomposition of time series into trend, seasonal, and remaining components within a longterm time series. An improved harmonic seasonal model which requires fewer observations has been presented in (Verbesselt et al., 2010b). Transform tools [START_REF] Jong | Analysis of monotonic greening and browning trends from global NDVI time-series[END_REF][START_REF] Martínez | Vegetation dynamics from NDVI time series analysis using the wavelet transform[END_REF] have also been used for analysis of the normalized difference vegetation index time series.

In this paper, we address the problem of change classification of multi-temporal SAR series. We will focus on samesensor same-incidence case and particularly on urban areas. A global processing chain in 3 steps, namely NORmalized Cut on chAnge criterion MAtrix (NORCAMA), is defined, as shown in Fig. 1. In the pre-processing step, a multi-temporal SAR image filter is used to reduce speckle phenomenon. After that, two change criteria based on likelihood ratio test combining noisy and denoised data are developed and compared. The final step is a clustering-and-recognizing classification based on the change criterion matrix, in which changes are classified into different types (including step change, impulse change, cycle change and complex change). The last 2 steps as well as the global frame work are the main contribution of this paper. Evaluation on synthetic and real images show the good performance of the proposed approach.

This paper is organized as follows. Section 2 briefly recalls the multi-temporal denoising method for SAR images. The proposed approximate likelihood ratio test and generalized likelihood ratio test change criteria are presented in Section 3. The proposed clustering-and-recognizing change classification method is then detailed in Section 4. This is followed by evaluation (Section 5) and conclusion (Section 6).

Pre-Processing: Multi-Temporal Denoising

When dealing with multi-temporal images, lots of information is available and useful for estimation in the time series. Based on this motive, we proposed a two-step probabilistic patch based (2SPPB) denoising method [START_REF] Su | Two steps multi-temporal non-local means for SAR images[END_REF] based on non local means [START_REF] Buades | A non-local algorithm for image denoising[END_REF] and probabilistic patch based weights (PPB) [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF] adapted to multitemporal SAR images. To allow a self-content reading of this paper and to introduce the useful notations, we briefly summarize the main steps of this approach. A complete description can be found in [START_REF] Su | Two steps multi-temporal non-local means for SAR images[END_REF]. It consists of a temporal averaging step and a spatial denoising step (summarized in Algorithm 1). Firstly, an average image is created by combining stable pixels while keeping unchanged the pixels not in accordance with the other dates (temporal averaging step). Then, on this improved image, a spatial denoising step is applied. A key point in both the temporal and spatial averaging is the weights based on the similarity between pixels, which are measured by similarity between patches. This section presents a brief summary of 2SPPB which will be useful for the following steps.

Pixel similarity

We denote by y t the observed SAR image, by y t (i) the noisy intensity value at pixel index i at time t, and by ût (i) the estimation of the actual pixel value u t (i) (the true value that we are looking for). Considering {y t 1 , y t 2 , . . . , y t N } as the stack of multitemporal images, the similarity S (i t , i t ′ , h, h ′ ) between pixels y t ′ (i) and y t (i) is defined through their patches y t (i) and y t ′ (i):

S i t , i t ′ , h, h ′ = S GLR (y t (i), y t ′ (i)) h + S KL ( ût (i), ût ′ (i)) h ′ (1)
S (i t , i t ′ , h, h ′ ) consists of a generalized likelihood ratio S GLR (GLR) from noisy images y t and y t ′ and a Kullback-Leibler divergence S KL (KL) from currently denoised images ût and ût ′ . S GLR and S KL are normalized by parameters h and h ′ . For any pair of pixels y 1 and y 2 , the GLR criterion S GLR is comparing two small square patches y 1 and y 2 of size K surrounding pixels at the positions of y 1 and y 2 .

S GLR (y 1 , y 2 ) = k∈K L 1,k + L 2,k log L 1,k y 1,k + L 2,k y 2,k -L 1,k + L 2,k log L 1,k + L 2,k -L 1,k log y 1,k -L 2,k log y 2,k (2)
where, y 1,k is pixel k in patch y 1 and L 1,k is the (equivalent) number of looks of y 1,k (idem for y 2,k and L 2,k ). The KL criterion S KL is computed iteratively from the denoised results ût and ût ′ (corresponding to the noisy images y t and y t ′ ). For any pair of denoised pixels û1 and û2 , S KL is also defined on two patches û1 and û2 .

S KL ( û1 , û2 ) = k∈K L 1,k û2,k û1,k + L 2,k û1,k û2,k -L 1,k -L 2,k +L 1,k ψ(L 1,k ) -ψ(L 2,k ) + ln( û1,k ) + ln( û2,k ) -L 2,k ψ(L 1,k ) -ψ(L 2,k ) + ln( û1,k ) + ln( û2,k ) (3)
where ψ(L) is the digamma function.

Two steps denoising

The first step of 2SPPB is to average the temporal pixels with binary weights:

y 1st t (i) = 1 Z t ′ ∈[t 1 ,t N ] ϕ S i t , i t ′ , h 1st , h ′ 1st • y t ′ (i) (4) ϕ S i t , i t ′ , h 1st , h ′ 1st = 1, if S i t , i t ′ , h 1st , h ′ 1st < 1 0, otherwise (5)
where, S i t , i t ′ , h 1st , h ′ 1st is computed from noisy images y t and y t ′ , and denoised images ûPPB t and ûPPB t ′ using PPB. The second step of 2SPPB approach is to exploit similar pixels in the temporally average image y 1st t rather than in the stack {y t 1 , . . . , y t N }. The estimation at time t is thus given by ût

(i) = 1 Z j∈Ω i exp(-S i t , j t , h 2nd , h ′ 2nd ) • y 1st t ( j) (6)
where the GLR criterion in S i t , j t , h 2nd , h ′ 2nd is computed using the temporal step (the first step) result y 1st t . To improve the ût estimation, the KL criterion is iteratively refined by using the previous estimate ût . The final result of this temporal denoising step will be denoted by ût for date t (corresponding number of looks Lt ). Algorithm 1 summarizes the multi-temporal denoising processing.

Algorithm 1 The multi-temporal denoising (2S-PPB) algorithm.

Input:

Registered temporal SAR images {y t 1 , y t 2 , . . . , y t N }.

A date t 1 of interest. Output:

ût 1 : the denoising result of image y t 1 .

----Step 1 (Temporal step):----1: for each y t in {y t 1 , y t 2 , . . . , y t N } do 2: denoise y t using PPB filter [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF]; 3: obtain pre-denoised results ûPPB t ; 4: end for 5: for each pixel index i do

6: Calculate y 1st t 1 (i); (Eq. 4) 7: end for ------------------ ----
Step 2 (Spatial step):-----8: Set ût 1 = 1; 9: for Iteration from 1 to 10 do 10:

for each pixel index i do 11:

Calculate S i t , j t , h 2nd , h ′ 2nd using y 1st t 1 (i) and ût 1 (i);

12:

Update ût 1 (i) using Eq.6;

13:

end for 14: end for ------------------15: return Denoised result ût 1 ;

Change Detection

In statistics, change detection problem can be considered as a comparison of two hypotheses H 0 and H 1 [START_REF] Radke | Image change detection algorithms: a systematic survey[END_REF]:

H 0 : u 1 = u 2 = u 12 (null hypothesis) H 1 : u 1 u 2 (alternative hypothesis) (7)
where H 0 is unchanged case and hypothesis H 1 is changed case. Likelihood ratio test is a state-of-the-art test which determines the decision by a likelihood ratio [START_REF] Kay | Fundamentals of statistical signal processing, volume II: Detection theory[END_REF]:

R = p(Y|u 12 , H 0 ) p(Y|u 1 , u 2 , H 1 ) ( 8
)
where Y is the observed data. In this section, we propose to develop the likelihood ratio test using the multi-temporal denoising results. Contrary to most likelihood ratio tests which only use noisy data, both the denoised data and the noisy data are involved in the proposed criteria.

Change Criterion by Approximate Likelihood Ratio Test

Recall that by y we denote a pixel intensity value. Under the speckle noise model described in [START_REF] Goodman | Some fundamental properties of speckle[END_REF], y is a realization of a random number characterized by the Gamma probability density function (pdf) p(y|u). u is the parameter of this Gamma pdf or the noise-free pixel value. Under the test in Eq.8, the change criterion between y 1 and y 2 using likelihood ratio test can be defined as:

R ALRT (y 1 , y 2 ) = p (y 1 , y 2 |u 12 , H 0 ) p (y 1 , y 2 |u 1 , u 2 , H 1 ) (9)
The criterion in Eq. 9 is a composite hypothesis problem because it requires the noise-free value u 1 and u 2 which are unavailable in practice. Instead of using usual generalized likelihood ratio extensions, we propose to use the estimated results ût to replace the noise-free u t (u 1 = û1 , u 2 = û2 and u 12 = û1 +û 2 2 ). Combined with the Gamma probability density function, the likelihood ratio in Eq. 8 becomes:

R ALRT (y 1 , y 2 ) = p (y 1 , y 2 |u 12 , H 0 ) p (y 1 , y 2 |u 1 , u 2 , H 1 ) = 1 4 û2 û1 + û1 û2 + 2 -L exp L y 1 û1 + y 2 û2 - 2y 1 + 2y 2 û1 + û2 ( 10
)
where L is the original spatially-invariant (equivalent) number of looks. Note that the approximate likelihood ratio R ALRT highly depends on the denoised values û1 and û2 , since R ALRT (y 1 , y 2 ) ≡ 1 when û1 = û2 .

Change Criterion by Generalized Likelihood Ratio Test

In a more general way, we can take into account the denoised values and consider the likelihood probability of H 0 and H 1 as p(y 1 , y 2 , û1 , û2 |H 0 ) and p(y 1 , y 2 , û1 , û2 |H 1 ). To simplify this likelihood probability, we can assume that {y 1 , û1 } and {y 2 , û2 } are independent, although this assumption is not well justified. Thus,

p(y 1 , y 2 , û1 , û2 |H 0 ) =p(y 1 , û1 |u 12 , H 0 )p(y 2 , û2 |u 12 , H 0 ) p(y 1 , y 2 , û1 , û2 |H 1 ) =p(y 1 , û1 |u 1 , H 1 )p(y 2 , û2 |u 2 , H 1 ) .
The likelihood ratio test is given by:

R GLRT (y 1 , y 2 ) = p(y 1 , û1 | u 12 , H 0 )p(y 2 , û2 | u 12 , H 0 ) p(y 1 , û1 | u 1 , H 1 )p(y 2 , û2 | u 2 , H 1 ) . (11)
Since u 12 , u 1 and u 2 are not available, they can be replaced by their maximum likelihood (ML) estimation:

u 1 = Ly 1 + L1 û1 L + L1 u 2 = Ly 2 + L2 û2 L + L2 ( 12
)
u 12 = Ly 1 + Ly 2 + L1 û1 + L2 û2 2L + L1 + L2 .
L, L1 and L2 are the number of looks of y, û1 and û2 respectively. Note that this is very similar to [START_REF] Lombardo | Maximum likelihood approach to the detection of changes between multitemporal SAR images[END_REF]. Nevertheless, the multi-temporal denoised values used in the proposed approach can provide more accurate estimation without loss of spatial resolution. It was not the case in [START_REF] Lombardo | Maximum likelihood approach to the detection of changes between multitemporal SAR images[END_REF] where spatial partitioning and averaging are introduced. In case of Gamma distributions with different number of looks, each probability term p(y, û|u) in Eq.11 can be approximated under conditional independence assumption by:

p(y, û|u) = p(y|u)p( û|u) = y -1 û-1 Γ(L)Γ( L) (Ly) L ( Lû) L u L+ L exp - Ly + Lû u . (13)
Finally, the change criterion by generalized likelihood ratio test is given by:

R GLRT (y 1 , y 2 ) = Ly 1 + L1 û1 L + L1 L+ L1 Ly 2 + L2 û2 L + L2 L+ L2 2L + L1 + L2 Ly 1 + L1 û1 + Ly 2 + L2 û2 2L+ L1 + L2 (14)
When û1 = û2 , the generalized likelihood ratio R GLRT still depends on the noisy values y 1 and y 2 which is different from the approximate likelihood ratio R ALRT .

Thresholds for Change Detection

In [START_REF] Kervrann | Optimal spatial adaptation for patch-based image denoising[END_REF], the authors proposed to define the parameters according to the quantiles of the similarity criterion when it is subject to identical and independent distributed random variables. Pursuing this idea, we propose to choose the thresholds according to the quantiles of R ALRT and R GLRT . The change detection threshold can be set by τ ALRT = quantile(R ALRT , α = 0.01) (and τ GLRT = quantile(R GLRT , α = 0.01)), which means the false alarm rate is 1%. However, it is not easy to obtain the distribution of R ALRT and R GLRT since they depend on the number of looks of noisy images, the number of images used in the denoising process and all the parameters of multi-temporal filter (such as the h, h ′ , search window size, patch size and so on). Thus, R ALRT and R GLRT distribu- tions are simulated using synthetic SAR images to choose the thresholds.

For our purpose of parameter setting, any picture can be used to generate multiple speckle images, such as the famous Lena image. Note that all the synthetic multi-temporal noisy images use the same true image, which guarantees no changes among them. The same number of images and the same number of looks as the real SAR images to be processed have to be used. Then, the multi-temporal denoising process in section 2 is performed. The approximate and generalized likelihood ratio test change criteria calculated from these images are considered as pure distributions of R ALRT and R GLRT . As shown in Fig. 2, the histograms of R ALRT and R GLRT are truncated by the thresholds (red lines) with false alarm 1%. The parts on the right of the thresholds are considered as unchanged case, the left part is changed case. Fig. 3 shows the R ALRT and R GLRT histograms of changed and unchanged pixels in real SAR data Paris (image information detailed in section 5). Those changed pixels are labeled manually. The robustness of the proposed change criteria R ALRT and R GLRT can be epitomized by the overlap of R ALRT and R GLRT histograms.

Change Classification

Change analysis between 2 dates aims at detecting a binary pattern (change or no-change). When dealing with a multitemporal data set (more than 2 dates), the analysis among them is much more complex. For instance, the temporal behaviors of a new building usually can be considered as a step change, which means that comparing the oldest date with other dates, it was unchanged at the beginning but it changed since a certain date (shown in Fig. 4.a). Similarly, we can define the boats in rivers or cars on the roads as impulse changes (Fig. 4.b). These change information can be used in the multi-temporal image interpretation tasks. Therefore, a clustering-and-recognizing method is proposed to classify changes into different types. This method consists of two steps, clustering using normalized cut on a change criterion matrix (to assign a same label to similar or unchanged temporal pixels) and classification according to their temporal behaviors. The following subsections detail the proposed approach.

Change Criterion Matrix (CCM)

At position i of a multi-temporal SAR series {y t 1 , ..., y t N }, we have the two pixel series {y t 1 (i), ..., y t N (i)} (original noisy data), {û t 1 (i), ..., ût N (i)} (denoised data by multi-temporal filter of section 2) and associated equivalent number of looks { Lt 1 (i), ..., Lt N (i)}. The change criterion matrix (CCM) at position i is defined as:

M(i) =               R(y 1 , y 1 ) R(y 1 , y 2 ) ... R(y 1 , y N ) R(y 2 , y 1 ) R(y 2 , y 2 ) ... R(y 2 , y N ) ... ... ... ... R(y N , y 1 ) R(y N , y 2 ) ... R(y N , y N )               (15)
where R(y n , y m ) shorts for R y t n (i), y t m (i) denotes the change criterion (R ALRT or R GLRT ) between pixel y t n (i) and y t m (i). Note that R(y n , y m ) = 1 when n = m. Contrary to the multi-date divergence matrix in [START_REF] Atto | Multidate Divergence Matrices for the Analysis of SAR Image Time Series[END_REF] performing at the image or sub-image level, the CCM presents the change information at pixel level. Each CCM M(i) denotes the temporal behavior of the pixel series at position i.

Clustering by Normalized Cut

Spectral clustering techniques make use of the similarity matrix of the data to perform clustering. Since the CCM can be considered as a similarity matrix of the time series, spectral clustering method has been applied on CCM to cluster the temporal pixels. In this case, no more similarity measurements is needed compared with other clustering methods (like K-Means algorithm for which new similarity to cluster center has to be computed).

Normalized spectral clustering proposed by [START_REF] Shi | Normalized cuts and image segmentation[END_REF] is employed in this work, which can be summarized in Algorithm. 2. In this algorithm, the unnormalized Laplacian matrix M u (i) is computed by:

M u (i) = M s (i) -M(i) (16) M s (i) =               R(y 1 , y n ) 0 ... 0 0 R(y 2 , y n ) ... 0 ... ... ... ... 0 0 ... R(y N , y n )               R(y m , y n ) = n=1,...,N R(y m , y n )
The only parameter in Algorithm. 2 is the number of clusters p.

Choosing the number of clusters p is a general problem for all clustering algorithms, and a variety of successful methods have been devised (more details in (Von Luxburg, 2007)). Eigengap heuristic is one of them and particularly designed for spectral clustering. The main idea is to choose the number p such that all eigenvalues λ 1 , ..., λ p are very small, but λ p+1 is relatively larger. However, this heuristic fails when the clusters of the data are overlapping (because of noise). To solve this problem, we binarize the CCM M(i) using the change detection threshold.

M b (i) =               R b (y 1 , y 1 ) ... R b (y 1 , y N ) R b (y 2 , y 1 ) ... R b (y 2 , y N ) ... ... ... R b (y N , y 1 ) ... R b (y N , y N )               (17) R b (y m , y n ) =        0 if R b (y m , y n ) < τ 1 if R b (y m , y n ) τ
The Eigengap heuristic performed on the binary change criterion matrix M b (i) can easily estimate the number of clusters p (see the example shown in Fig. 5). The toy model in Fig. 5.a has 2 clusters. Using CCM M(i) in 5.b, the difference between eigenvalues λ 2 and λ 3 is not large enough compared with the one between λ 1 and λ 2 . It is very easy to find the best estimation of p using the binary CCM M b (i) in 5.c because of the large gap between λ 2 and λ 3 . It is obvious that this estimation of p highly depends on the choice of the thresholds. However, the robustness of the proposed change criteria (especially R GLRT shown in Fig. 3) can guaranty the estimation accuracy of p.

Recognizing

After clustering, each pixel series {y t 1 (i), ..., y t N (i)} has a cluster label series {l t 1 (i), ..., l t N (i)}, in which l t n ∈ {1, ..., p}. We can identify different types of change according to the transformation in the cluster label series {l t 1 (i), ..., l t N (i)}. For example, if p = 1, there is no change among this pixel series. If p = 2 with cluster label series {1, 1, ..., 1, 2, 2, ..., 2}, it is a step change.

Algorithm 2 Clustering of the pixel series (Normalized spectral clustering [START_REF] Shi | Normalized cuts and image segmentation[END_REF]) Input:

A change criterion matrix M(i) of pixel series {y t 1 (i), ..., y t N (i)}, number p of clusters to construct.

Output:

The clustering labels {l t 1 (i), ..., l t N (i)} for pixel series {y t 1 (i), ..., y t N (i)} 1: Compute the unnormalized Laplacian matrix M u (i) using Eq.16. 2: Compute the first p generalized eigenvectors v 1 , ..., Impulse change usually has p = 2 and cluster label series is {1, 1, ..., 1, 2, 2, ..., 2, 1, 1, ..., 1}. When p 3, the transformation is complex and changes are defined as complex case. According to these identifications (details in Table 4.3), changes can be classified into several classes.

v p of M u (i) (M u (i)v = λIv). 3: Let v ′ be

Experiments

The proposed methods are evaluated on both synthetic images and real multi-temporal SAR images.

Experiments of Change Detection

Data Set

Synthetic images: Fig. 6.a show the noisy synthetic images y t and y t ′ corrupted by single-look multiplicative speckle noise respectively (with the ground truth of changes between them in Fig. 6.a). The four squares are 32×32 pixels with true value 128. The darker frame is 8 pixels width with 32 as true value and the true value of background is 64.

Realistic SAR Synthetic Images: A denoised image of 21 single-look TerraSAR X-band images in Paris (France) sensed in 2011 is considered as the noise-free image (multi-temporal denoising approach in section 2), as shown in Fig. 6.b. Two single-look noisy images y t and y t ′ are generated with changes added in y t ′ . These changed regions are about 15-25 pixels width and length, for instance a 20×20 pixels patch of vegetation is replaced by a same size patch of building and so on. Fig. 6.b shows the ground truth of changes.

Real SAR images 1: 26 single-look TerraSAR images in Saint-Gervais-les-Bains (France) (13 images are sensed in 2009 and the other 13 images in 2011) are shown in Fig. 7.a-c, identified as Saint-Gervais-les-Bains. Reference Ground truth of changes is labeled manually in Fig. 7.c.

Real SAR images 2: Experiment in Fig. 7.d-f uses 21 singlelook TerraSAR X-band images identified as Paris in Paris (France) sensed in 2011. We label the ground truth of changes manually, as shown in Fig. 7.f.

Real SAR images 3: Experiment in Fig. 12.a-d uses 24 CARABAS-II magnitude images acquired in Vidsel, Sweden 2002, identified as CARABAS ( Sensor Data Management System (SDMS) Public web site, 2008). We only detect the changes between image v02 2 1 1 and image v02 4 1 1, while all the 24 images are used in the multi-temporal denoising process.

Real SAR images 4: Experiment in Fig. 12.e-h uses 9 single-look TerraSAR X-band images identified as Sendai in Sendai Harbor (Japan) sensed in 2011. Fig. 12.e and f show the images acquired respectively on 06/05/2011 and 08/06/2011. All the 9 images are used in the multi-temporal denoising step.

Change detection methods

The proposed change criteria approximate likelihood ratio test R ALRT and generalized likelihood ratio test R GLRT are compared with some state-of-the-arts methods, such as Log-Ratio operator [START_REF] Rignot | Change detection techniques for ERS-1 SAR data[END_REF], the generalized likelihood ratio test (GLRT) proposed in [START_REF] Lombardo | Maximum likelihood approach to the detection of changes between multitemporal SAR images[END_REF], Wilcoxon Test based change criterion [START_REF] Krylov | Change detection with synthetic aperture radar images by Wilcoxon statistic likelihood ratio test[END_REF] and Method for generalIzed Means Ordered Series Analysis [START_REF] Quin | MIMOSA: An Automatic Change Detection Method for SAR Time Series[END_REF], summarized in Tab.5.1.2.

Results

The change detection results are assessed by the True-Positive versus False-Positive curves using the reference map of changes. Fig. 6.c and Fig. 7.c show the the True-Positive versus False-Positive curves. The proposed methods R ALRT and R GLRT can generally obtain higher ROC curves than others. The experiments on CARABAS and Sendai data in Fig. 9 and 10 show that the proposed GLRT change detection has comparable performance with MIMOSA [START_REF] Quin | MIMOSA: An Automatic Change Detection Method for SAR Time Series[END_REF]. R ALRT and R GLRT outperform other change criteria, but the latter is more reliable than the former (the ROC curves of R GLRT are higher than R ALRT in Fig. 6 and 11).

Experiments of Change Classification 5.2.1. Test on realistic SAR synthetic images

This experiment use one denoised image of 21 single-look TerraSAR X-band images of Paris (France) sensed in 2011 as the noise-free image (multi-temporal denoising approach in section 2), as shown in Fig. 6

Confusion

Classification results Matrix (%) Unch.

Step

Test on real SAR images

We have 21 single-look TerraSAR X-band images identified as Paris in Paris (France) sensed in 2011 and 6 single-look Ter-raSAR X-band images identified as San-Francisco sensed in San-Francisco, U.S.A. 2007 and 2011. These images have been accurately registered using the sensor parameters.

Conclusion

In this work, a global framework NORmalized Cut on chAnge criterion MAtrix (NORCAMA) for change classification of multi-temporal SAR time series has been presented. To reduce the effect of speckle, a multi-temporal denoising approach is applied in the pre-processing step. Using both noisy data and denoised data, the approximate likelihood ratio and the generalized likelihood ratio are computed as change criteria. The change classification is performed by clustering on change criterion matrix and classifying of label transformation. Different types of change, like step changes, impulse changes and cycle changes, have been defined by the proposed method, which can be used for multi-temporal SAR image interpretation. Note that the methods in the 3 steps of the proposed framework can be replaced by alternative methods, for instance other denoising approaches for the pre-processing step and other change criteria for the change detection step.

The future work will be focused in introducing spatial information into change analysis such as combining types of change and shapes to identify objects and developing a temporal-spatial SAR image analysis framework. The spatial information should allow us to extend this approach to multi-incidence images.

Figure 1 :

 1 Figure 1: The global diagram of the proposed framework NORmalized Cut on chAnge criterion MAtrix (NORCAMA). It consists of 3 steps, from left to right: 1) the pre-processing step using a multi-temporal SAR filter (2S-PPB) (Su et al., 2014) to denoise; 2) change detection step using the proposed change criteria based on likelihood ratio test; 3) change classification by the proposed clustering-and-recognizing method.

Figure 2 :

 2 Figure 2: The simulated histograms of R ALRT and R GLRT using synthetic images. The red lines are thresholds τ ALRT and τ GLRT with false alarm 1%. The blue lines are the histograms of unchanged R ALRT and R GLRT .

Figure 4 :

 4 Figure 4: Examples of step changes and impulse changes. From left to right: original multi-temporal SAR images at time t 1 , t 2 , t 3 , t 4 , t 5 , change criterion matrix of a pixel in the red rectangle (cold color: unchanged; warm color: changed).

Figure 3 :

 3 Figure 3: An example of the normalized histograms of R ALRT and R GLRT using real SAR images Paris (detailed in section 5). The red lines are thresholds τ ALRT and τ GLRT with false alarm 1%. The blue lines are the histograms of unchanged R ALRT and R GLRT . The green lines are the histograms of changed R ALRT and R GLRT .

  A pixel series {y t 1 (i), ..., y t 12 (i)} which should be clustered into 2 groups (red and blue).

  ,...,12) (b) CCM M(i) and its eigenvalues.

  Binary CCM M b (i) and its eigenvalues.

Figure 5 :

 5 Figure 5: Estimation of number of clusters. (a). a pixel series {y t 1 (i), ..., y t 12 (i)} which should be clustered into 2 groups; (b). Estimation of the number of clusters using CCM M(i), the gap between λ 2 and λ 3 is not obvious; (b). Estimation of the number of clusters using binary CCM M b (i), the gap between λ 2 and λ 3 is larger.

  .a. 6 single-look images are generated with different changes added in them. As shown in Fig.6.c, different kinds of changes have been introduced, such as step change (in red), impulse change (in green) and cycle change (in blue). Fig.6.d shows the change classification result by R GLRT with confusion matrix shown in Tab.5.2.1.

  (a) Saint-Gervais-les-Bains data set. From left to right: noisy image y 1 , noisy image y 26 and the reference map of changes. (b) Paris data set. From left to right: noisy image y 1 , noisy image y 21 and the reference map of changes.

  positive alarm vs true positive curves of Saint-Gervais-les-Bains and Paris data set.

Figure 7 :

 7 Figure 7: Change detection results for real SAR images Saint-Gervais-les-Bains and Paris data set.

Figure 8 :

 8 Figure 8: Change detection results.

Figure 9 :

 9 Figure 9: Change detection results of real SAR images CARABAS ( Sensor Data Management System (SDMS) Public web site, 2008).

Figure 10 :

 10 Figure 10: Change detection results of real SAR images Sendai.

  Fig.12.a and b only show the first noisy image and its denoising result. Fig.12.c shows the results of the change classification approach by R GLRT , in which red regions denote step changes, green are impulse changes and blue are cycle changes. We can observe that many boats in river have been classified as impulse change. To analyze the performance of change classification, some interesting regions (examples of step change, impulse change, cycle change and complex change) of all the dates have been shown in Fig.13.

  (a) Example of step change (some sports facilities in AT&T ballpark) in San-Francisco. From left to right, the change classification results, the noisy image from date t 1 to t 6 . (b) Example of impulse change (two ships at the pier 48) in San-Francisco. From left to right, the change classification results, the noisy image from date t 1 to t 6 . (c) Example of cycle change (the bank of the river) in San-Francisco. From left to right, the change classification results, the noisy image from date t 1 to t 6 . (d) Example of complex change (the boats at the pier) in San-Francisco. From left to right, the change classification results, the noisy image from date t 1 to t 6 .

Figure 13 :

 13 Figure 13: Details of change classification results of Fig.9.b. From top to bottom: example of step change, impulse change; cycle change and complex change. (black: no change, red: step change, green: impulse change, blue: cycle change and yellow: complex change)

Table 1 :

 1 the matrix containing the vectors v 1 , ..., v p as columns. 4: Consider each row of v ′ as a sample, v ′ = {v ′ 1 , ..., v ′ N }. 5: Cluster the samples v ′ 1 , ..., v ′ N with the k-means algorithm into clusters with p as the number of clusters. The cluster labels of v ′ 1 , ..., v ′ N are l 1 , ...l n , ..., l N (l n ∈ {1, ..., p}). 6: return Cluster labels l 1 , ...l n , ..., l N The identifications of different types of change.

	Types	p	Label series {l t 1 (i), ..., l t N (i)}
	Unchanged	1	{1, 1, ...1}
	Step	2	{1, 1, ...1, 2, 2, ...2}
	Impulse	2	{1, 1

, ...1, 2, 2, ...2, 1, 1, ...1} Cycle 2 {1, ...1, 2, ...2, 1, ...1, 2, ...2, ...} Complex 3 {1, 1, ..., 2, 2..., 3, 3...4, 4...}

Table 3 :

 3 Confusion matrix of change classification results. Unch.: unchanged, Step: step change, Impl.:impulse change, Cyc.: cycle change and Comp.: complex change.

Table 2 :

 2 The change detection methods used in the comparison experiments.

Acknowledgment

The authors would like to thank Jean-Marie Nicolas for his help on image registration, the IGARSS Data Fusion Contest for providing the multi-temporal SAR images and the Department of Radar Systems, Swedish Defense Research Agency (FOI), for the CARABAS-II data distribution. They also wish to acknowledge the German Aerospace Agency (DLR) for the Terra-SAR-X images (project MTH0232 and LAN1746).