Dorra Abid

Ahmed Ktari
email: ahmed.ktari@enis.tn

Dhouha Mellouli

Nedia Gafsi

Nader Haddar

Failure analysis of shot-sleeves used in brass high pressure die-casting process

Keywords:

A failure investigation has been conducted on a two cases of shot-sleeve used in brass-die casting and made up of AISI H10 tool steels to study their failure mechanisms. The chemical composition of the shot-sleeves material and the hardness profiles were evaluated. A preliminary examination of shot-sleeves reveal the presence of cracks network on their inner surfaces which proves that thermal fatigue was probably the main cause of their failure. A Meticulous investigation of these damaged surfaces reveal the presence of an additional small zone cited in the vicinity of the plunger entry side. This zone presents several scratches sign of abrasive wear. Then, the measurement of the cracks length and their linear density along longitudinal and transversal cuts of the damaged shot-sleeve sample are carried out. Results show that the observed cracks network can be divided into two zones. The first one was the most damaged in terms of cracks density and length given that it is the first zone which enter in contact with the molten metal. However, the cracks network examined in the second zone appears to be superficial.

Introduction

The high pressure die-cast process is used to produce parts mostly from aluminium, magnesium, zinc and copper alloys by injection of the molten metal in the mold cavity. Parts produced by this process conform accurately to the die size, have favourable mechanical features, and are low in cost. Hence, this process has a wide range of application fields including the aircraft and automobile industries. During hot forming operations the tools steel used to this process are exposed to a sever thermal and/or mechanical loads which can lead to their damage. The improvement of these structures lifetime was strongly required due to several economic and environmental reasons. One of the major damage mechanisms occurring in hot forming tool steels, under cyclic thermal loads (e.g. die casting and forging [START_REF] Miquel | Heat checking of hot work tool steels[END_REF][START_REF] Srivastava | Computer modeling and prediction of thermal fatigue cracking in die-casting tooling[END_REF][START_REF] Persson | Thermal fatigue cracking of surface engineered hot work tool steels[END_REF][START_REF] Le Roux | Role of heat-flux density and mechanical loading on the microscopic heat-checking of high temperature tool steels under thermal fatigue experiments[END_REF]), is the formation of a network of interconnected cracks, often named Heat-Checking [START_REF] Jean | An investigation on heat checkingof hot work tool steels[END_REF][START_REF] Medjedoub | Endommagement par faïençage des moules de fonderie sous pression en fatigue thermique: origines, mécanismes et approches[END_REF][START_REF] Norström | Thermal fatigue behaviour of hot work tool steels[END_REF][START_REF] Persson | Simulation and evaluation of thermal fatigue cracking of hot work tool steels[END_REF]. The stress cracking, which presents another variant of thermal fatigue cracks, was clearly marked in areas exposed to local stress concentrations and can lead to crack initiation in a die casting mold [START_REF] Danzer | Thermal fatigue cracks in pressure die casting dies[END_REF]. Then, these cracks can grow and became more pronounced driven by several factor including thermal fatigue, erosion, oxidation, soldering of the metal to the die surface, deformation of die contact surfaces and dangerous fracture [START_REF] Kosec | Analysis of casting die failures[END_REF][START_REF] Neu | Thermo-mechanical fatigue, oxidation, and creep. Part 1: damage mechanisms[END_REF][START_REF] Panjan | Improvement of die-casting tools with duplex treatment[END_REF][START_REF] Mellouli | Hardness effect on thermal fatigue damage of hot-working tool steel[END_REF].

In hot metal forming applications like stamping or die casting the service life time of the tools are limited due to their extreme working conditions in term of thermal and/or mechanical loads which result from the close contact between the tools and the hot work piece. In these applications the tools temperature can reach 1200 °C [START_REF] Mellouli | Hardness effect on thermal fatigue damage of hot-working tool steel[END_REF] when the molten metal attain 675 °C and 930 °C for Al-alloy casting and Cu-alloy casting respectively [START_REF] Tidlund | Heat checking in hot work steels, 2nd International Colloquium on Tool Steels for Hot Working[END_REF]. Under such working conditions, the tools are generally damaged through wear and thermal fatigue cracking processes [START_REF] Cser | Tool life and tool quality in bulk metal forming[END_REF][START_REF] Skumavc | Thermal fatigue study of tungsten alloy WNi28Fe15 cladded on AISI H13 hot work tool steel[END_REF][START_REF] Starling | Thermal fatique of hot work tool steel with hard coatings[END_REF][START_REF] Malm | Material-related model for thermal fatigue applied to tool steels in hot work applications[END_REF]. Ktari et al. [START_REF] Ktari | Fatigue fracture expertise of train engine crankshafts[END_REF] have investigated a diesel engine crankshafts fracture used in train. Mellouli et al. [START_REF] Mellouli | Thermal fatigue failure of brass die-casting dies[END_REF] have investigated the thermal fatigue failure of a prematurely failed brass die-casting.

To improve the structures lifetime several researches have been conducted to study the effect of surface coatings on those structure lives. Panjan et al. [START_REF] Panjan | Improvement of die-casting tools with duplex treatment[END_REF] have studied the effect of duplex treatment composed of plasma nitriding and the deposition of 4mm-thick PVD CrN coating on die-casting tools. They have found that wear of duplex treated tool is smaller than wear of the plasma nitrided tool. In addition, Srivastava et al. [START_REF] Srivastava | A multilayer coating architecture to reduce heat checking of die surfaces[END_REF] have studied the effect of three-layer coating architecture on die surface: the outer layer consists of a thermal barrier coating composed of rare-earth oxide, the middle one is a TiAlN diffusion barrier coating and the inner layer is a thin adhesive Ti layer. These authors have found that this coating can significantly improves the thermal fatigue resistance of the substrate. Klobčar et al. [START_REF] Klobčar | Thermo fatigue cracking of die casting dies[END_REF] have studied the effect of cladding with maraging steels; they have found that this treatment provides good thermal fatigue resistance to the treated metal. The functionally graded materials (FGM) were also applied as a solution to improve the structure lifetime. Fazarinc et al. [START_REF] Fazarinc | Thermal fatigue properties of differently constructed functionally graded materials aimed for refurbishing of pressure-die-casting dies[END_REF] have studied several FGM surfaces by changing amounts of alloying elements, precisely Si and Mo, they have found that their resistance to thermal fatigue is estimated at 27 times more resistant than the basic material.

In our case of study, two failed structures named shot-sleeve have been examined. This structure is a part of the die-casting injection press sited on the fixed mold side between the injection plunger and the sprue bush. In service, during the hot forming operations, the shot-sleeve undergo severe cyclic thermal-mechanical loads (i.e. the alternating of compressive and tensile stresses at the surface that arise from differential thermal expansion / contraction during sudden transient temperature changes [START_REF] Cser | Tool life and tool quality in bulk metal forming[END_REF][START_REF] Jean | An investigation on heat checking of hot work tool steels[END_REF]) which can lead to cracks initiation and propagation. The investigation of the thermo-mechanical fatigue shot sleeve was conducted according to the following points:

• Chemical and mechanical characterization of the shot sleeve material.

• Analysis of the failed zones to understand the main causes and mechanisms leading to the shot sleeve structures.

Shot sleeve materials characterization

In the present study, samples of shot-sleeve failed in service, were investigated to understand the mechanisms, which can lead to their damage. The shot-sleeve (Fig. 1), used in brass high pressure die-casting process, undergoes a severe thermo-mechanical loadings during cast cycles. Each cycle include principally four phases as following:

(i) The filling of the shot sleeve with the molten metal by means of a ladle, (ii) The plunger moves through the shot-sleeve and press the casting material into the mold, (iii) The cooling phase begin once the mold has been filled with molten metal and keep on until the casting cold. During the first part of this phase, named compression, a pressure is applied to the cooling metal by means of the plunger until the alloy has solidified. (iv) The mold is then opened and the casting is removed with the ejector pins and the internal surfaces of the mold are sprayed with a lubricant before every cycle.

After the casting has been removed and the lubricant applied to the mold surfaces, the die are clamped together again then the cycle will repeat itself. The high pressure brass die-casting process parameters are summarized Table 1. Table 3 Monotonic mechanical properties of AISI H10 steel grade at room temperature.

Chemical analysis and microstructure

The chemical analysis of a sample of steel taken from damaged shot-sleeve was determined by means of a spectrometer (Jobin Yvon JY 48®). The typically steel grades used for hot work tools and our studied material are summarized in Table 2.

Materials used in high pressure die-cast process equipments have to combine a multitude properties to resist production constraints that they face. According to Reynoldson [START_REF] Reynoldson | Die Casting Tools Manufactured From H13 Hot Work Tool Steels[END_REF], a good material need to have: (i) low chemical reactivity with the surrounding environment (molten metal and oxygen), (ii) low thermal expansion and high thermal conductivity coefficients, (iii) high tempering resistance and (iv) a good toughness.

Both steel grades H10 (32CrMoV12-28) and H13 (X40CrMoV5-1) possess good resistance to abrasive effects, low sensibility to thermal shocks and good dimensional stability. However, H10 presents relatively high wear resistance compared to AISI H13 steel grade. The contents in chromium, molybdenum and vanadium give good hardenability, wear resistance, fracture toughness and hardness properties at high temperatures [START_REF] Bayrakçeken | Failure of a single cylinder engines crank shafts[END_REF]. The chromium increases the hardenability of the steel, decreases the grain growth during austenitization, delay the softening during the tempering and contribute to the reduction of the high-temperature oxidation [START_REF] Lévêque | Données numériques sur les aciers à outils[END_REF]. The molybdenum slows down the steel softening and give it a very good wear resistance at high temperature through the presence of hard metallic carbides M 6 C and M 2 C with a hardness of 1500 and 2000 HV respectively [START_REF] Lévêque | Données numériques sur les aciers à outils[END_REF]. The vanadium allows to generate a very high carbides hardness (MC ≈ 3000 HV). It added generally, with a small amount, with chromium, molybdenum and tungsten to avoid the harmful oxidation resistance beyond 600 °C. The decrease of corrosion resistance is compensated by the presence of chromium [START_REF] Lévêque | Données numériques sur les aciers à outils[END_REF].

Our studied material is H10, a metallographic observations were realized by means of an optical microscope (Leica®) on the same taken sample polished until mirror state then attacked by Nital (5% of nitric acid diluted in methanol). These observations show the presence of a martensitic structure (Fig. 2).

Mechanical properties 2.2.1. Tensile properties

In order to investigate the tensile properties of the studied material AISI H10, tensile tests were carried out on a three standard flat tensile specimens (NF A 03-151) as presented in Fig. 3. The semi-finishing form of each one was obtained from the damaged shot-sleeve case using a wire EDM machine. Then, specimens flat surfaces were mechanically grinding to obtain a good surface finish. These tests were conducted on a universal testing machine with capacity of 50 kN and a crosshead speed of 0.5 mm/s. The obtained results average are summarized in Table 3.

The material toughness

Standard Charpy V-notched specimens were machined from the failed shot-sleeve by a wire EDM machine (Fig. 4). Then, three tests were conducted at room temperatures to determine the material toughness. The average of the measured CVN energy (CVN = 3.4 ± 0.1 J) and the feature of the fractured surface prove the brittle nature of the specimen material (Fig. 5).

Hardness analysis

The hardness profiles were measured, on a shot-sleeve longitudinal section (Fig. 6 (b)) polished until mirror state, from the entry plunger zone to the mold side and from the inner to the outer surface (Fig. 7). These profiles are performed with a Vickers hardness machine using a load of 15 kg. The measured profiles reveal that the hardness appears to be constant for all tested samples.

The original hardness of the shot-sleeve as provided by manufacturer technical data-sheet was about 42-44 HRC. However, the average of measured hardness is about of 44.8 -2.4 +3.8 HRC. This shows that the hardness was remained quasi-constant even after a few thousands of thermal heating and cooling phases (more than 74,000 cycles) and which indicate that the mechanical behaviour of the material is stable (i.e. absence of significant cyclic hardening or softening) under service conditions. Furthermore, micro-hardness measurements were performed with a Wilson hardness machine 402MVD using a load of 1 kg at a distance of approximately 400 μm from the inner surface. The average of measured micro-hardness is about 49.8 -6.4 +5.5 HRC. This result proves that is no significant increase of the hardness from the heart to the inner surface of the shot-sleeve.

Expertise of the damaged shot-sleeves

The metallic structures, which undergo severe cyclic thermo-mechanical loads in service, are typically characterized by the presence of several cracks. These cracks are the most dangerous defects in term of lifetime on which we have to grant a particular attention. In our study, two shot-sleeves damaged in service, supplied by SIAM company, were investigated to determine the The first and the second shot-sleeves have been failed after respectively 46,000 and 74,000 cycles. This corresponds to a quantity of injected molten metal of 18,000 kg and metal 27,000 kg which give an average volume/cycle of 0.39 kg and 0.36 kg, respectively. In addition, only the brass was used as molten metal in these two studied shot-sleeves. According to these previous data, the first shot-sleeve was failed at smaller number of cycles. Which can be explained by the fact that the average of the injected volume per cycle is significantly more important for the first shotsleeve compared to second one (more than 8.3%).

According to the chemical composition, the microstructure and the preliminary damage observations carried on both shot-sleeves it is obvious that both shot-sleeves are made from the same steel grade (AISI H10) and present the same damage morphology; however, the second shot-sleeve present a more pronounced crack network. As a result, the second shot-sleeve was chosen to conduct the subsequent of the investigation.

Observation of the damaged surfaces

The second shot-sleeve was cut with a wire EDM machine according to their longitudinal direction to examine the integrity of the damaged surface (Fig. 8). After that, the cut-out parts are subjected to a surface cleaning operations: A chemical pickling with an alkaline solution using a dilute nitric acid was performed, at room temperature, to eliminate all traces of oxides. Then, the damaged shot-sleeve surfaces were subjected to a preliminary assessment with a naked eye to identify their major damage features. The examined surfaces (Fig. 9) clearly show the presence of a cracks network on the inner shot-sleeve surfaces.

A finer observations of the shot-sleeve inner side were carried out after a longitudinal cut (Fig. 6 (b)). The examination of the damaged surfaces was made using an experimental device equipped with a digital Baumer® camera GIgE (8 bits). This camera is equipped with a CCD sensor allowing to have a resolution of 1624 × 1236 pixels. Several photos are taken and then assembled by means of Photoshop® software to obtain a large image that cover the entire shot-sleeve examined surface (Fig. 9). This surface reveals the presence of three different zones. The first one is located in the plunger entrance side. It presents a number of scratch according the longitudinal direction which corresponds to abrasive wear caused by the friction between the plunger and the shot-sleeve. The second one corresponds to the first zone which enter in contact with the molten metal (T = 950 °C). It undergo cyclic thermal shocks. This zone appears to be the most damaged surface in the shot-sleeve. It presents a dense crack network sign of the thermal fatigue. In addition, the investigation of this zone according to the OθZ plan shows that some cracks, after initiation, were propagated in depth along the radial direction (Detail A). The last zone reveals also the presence of a crack network however it is less dense compared to this observed in the 2nd zone (Fig. 9).

Quantitative study of the cracks network

In this section, a meticulous investigation on the cracks distribution and density was performed by means of longitudinal and transversal observations (Figs. 10 and 11). The cracks depth according to the axial distance was presented in Fig. 10. It is to be noted that the cracks depth trend along the shot-sleeve follows the Gaussian distribution. Fig. 10 shows that the cracks depth is maximal around the zone 2, which corresponds to the casting zone (fall zone). This zone is characterized by a maximum crack length and average values of 6.23mm mm and 1.66mm respectively. In zone 3 the cracks deep is less than those observed in zone 2 nevertheless their density is more pronounced (i.e. 1.59crack/mm in zone 3 against 0.69crack/mm in zone 2). The presence of the shielding effect between neighbouring cracks existing in zone 2 cause the propagation slows down of the smaller cracks and only a few of the biggest ones continue to propagate. Such a phenomenon was observed in many applications in which the thermal fatigue was the main cause of damage [START_REF] Ktari | Fatigue fracture expertise of train engine crankshafts[END_REF][START_REF] Mellouli | Thermal fatigue failure of brass die-casting dies[END_REF][START_REF] Maillot | Thermal fatigue crack networks parameters and stability: an experimental study[END_REF].

The cracks depth investigation from the centre of shot-sleeve (θ = 0°) to the shot-sleeve horizontal plane (θ = 45°) for every shot-sleeve cross-sections A, B and C (Fig. 11). It is obvious that the maximal depths are sited at the level of the zone of casting. This result can be explained by thermal fatigue, oxidation and erosion by the molten brass flow which results in removal of superficial layer. The initial level of the shot-sleeve casting zone (cited in zone 2) was decreased probably due to oxidation and erosion effects as shown in Fig. 12. The cracks situated surroundings the shot-sleeve centre (near θ = 0°) are in reality longer compared to their effective length due to the small decrease of the shot sleeve inner surface level produced after the oxidation / erosion effect. Others micro cracks are totally removed what favoured the propagation of persistent cracks because the motive force was located in them, is confirmed by previous research [START_REF] Miquel | Heat checking of hot work tool steels[END_REF][START_REF] Medjedoub | Effect of local stress on the heat-checking morphology in high temperature tool steels under thermal fatigue: transition from multi-axiality to uniaxiality[END_REF].

In addition, in order to more understand the mechanisms leading to the shot sleeve failure, a sample of the cross-section surface was examined using SEM (Fig. 13) and Energy Dispersive X-ray Spectroscopy (EDS) to identify the chemical composition of the metal inside cracks (points 1 in Fig. 13 (b)) (Table 4). According to Fig. 13 (a), the cracks tend to propagate perpendicular to the damaged surface. The initial growth of thermal fatigue cracks is facilitated by an oxidation attack on the crack surface, which forms the Fe, Cu, Zn, Cr, Mo, Si and O rich layer (Table 4). The presence of brass and oxides, in the cracks increases plastic yield (Rp 0.2 ) and compressive stresses during the heating phase of the thermal cycle [START_REF] Klobčar | Thermal fatigue of materials for die-casting tooling[END_REF]. As a result, the tensile stress is further increased during the cooling phase and can leads to the acceleration of the crack growth rate. In addition, the tensile stresses applied on the shot-sleeve inner surface during the cooling phase combined with the erosion phenomenon caused by the molten brass fall in the second zone leads to a local failure, which promote the opening of thermal cracks [START_REF] Mellouli | Hardness effect on thermal fatigue damage of hot-working tool steel[END_REF]. These cracks act as channels allowing oxygen to penetrate down to the crack tip area and oxidize both shot sleeve and the filling material. Solidified particles down on the crack act as an obstacle for its closer and the size of the crack openings in the filling material indicates the amount of brass that will fill the crack during the next cycle [START_REF] Mellouli | Thermal fatigue failure of brass die-casting dies[END_REF].

Conclusions

The investigation of a failed shot-sleeves used in high pressure brass die-casting was carried out in this study. The failed shotsleeve was at first investigated by the determination of chemical composition, microstructures and hardness. The cracks surface network, observed at the inner surfaces of all the examined shot-sleeves, were shown the presence of three different zones: (i) entered piston, (ii) entrance of the molten metal and (iii) exit of the molten metal. These zones differ by their feature, cracks network density and depths. The longitudinal section shows that the second zone possess the maximum values of cracks depth. In addition, the cross section shows that the maximal cracks length are situated in the vicinity of the shot-sleeve centre (θ ≈ 0°). According to the foregoing discussion it is clear that thermal fatigue is the first cause of shot-sleeve damage. However other phenomena like oxidation and erosion of the shot-sleeve inner surface cause the increases of the damage rate mainly in the casting zone (fall zone).

Fig. 1 .

 1 Fig. 1. Geometry and dimensions of the studied shot-sleeves.

Fig. 2 .

 2 Fig. 2. Microstructure of the shot-sleeve: (a) Light Optical Microscopy (LOM), (b) The SEM microanalysis.

Fig. 3 .

 3 Fig. 3. Dimensions of tensile specimen (in mm).

Fig. 4 .

 4 Fig. 4. Charpy V-notch specimens: (a) Geometry and dimensions (in mm), (b) machined Charpy V-notch specimens.

Fig. 5 .

 5 Fig. 5. Fracture surface of Charpy specimen. (a) macroscopic observation, (b) binocular microscope observation.

Fig. 6 .

 6 Fig. 6. CAO schema of shot-sleeve cutting planes: (a) Initial shot sleeve shape, (b) After the longitudinal EDM cut, (c) Cross-section cut plans.

Fig. 7 .

 7 Fig. 7. Hardness shot-sleeve profiles, (a) from the entry to the mold zones, (b) from inner to the outer surfaces.

Fig. 8 .

 8 Fig. 8. Cracks network observed on the inner surface of the damaged shot-sleeve.

Fig. 9 .

 9 Fig. 9. Shot-sleeve inner surface damage observed with a digital camera.

Fig. 10 .

 10 Fig. 10. Crack lengths distribution along the shot sleeve (longitudinal sectionθ = 0).

Fig. 11 .

 11 Fig. 11. Crack distribution of damaged shot-sleeve in cross-sections A, B and C.

Fig. 12 .

 12 Fig. 12. Oxidation / erosion effect in the fall zone (Zone 2) of the shot-sleeve.

Fig. 13 .

 13 Fig. 13. SEM observation of a damaged shot sleeve cross-section: (a) Short (initiation) and deep crack, (b) Zoom on the examined crack with EDS (point 1).

Table 1

 1 Process parameters in high pressure brass die-casting case of study.

	Shot-sleeve filling	2.5 s
	Injection phase	0.1 s
	Cooling phase	27 s
	Compression	4 s
	Ejection and mold coating	23 s
	Molten metal temperature	950 °C
	Cast ejection temperature	300 °C

Table 2

 2 Chemical

		C	S i	M n	P	S	C r	M o	V
	Studied material	0.304	0.242	0.303	0.025	0.005	2.89	2.71	0.61
	AISI H10	0.28-0.35	0.1-0.4	0.15-0.4	0-0.03	0-0.02	2.7-3.

composition of the studied material and usual steels grades used in hot work tools (% in weight).

Table 4

 4 The chemical composition inside crack area (EDS).

		Fe	Cu	Zn	O	Si	S	Cr	Mo
	Point 1	5.61	57.03	24.09	1.58	0.13	0.32	0.39	0.77

Acknowledgements

The authors would like to thank Mr. H. Ammar the Head manager of SIAM for fruitful discussions about shot-sleeves failure mode, the material support and for chemical analysis of the studied steel samples.


Xin Su
email: xin.su@telecom-paristech.fr

Charles-Alban Deledalle
email: charles-alban.deledalle@math.u-bordeaux1.fr

Member, IEEE Hong Sun
email: hongsun@whu.edu.cn

F Tupin

C.-A Deledalle

Two-Step Multitemporal Nonlocal Means for Synthetic Aperture Radar Images

Keywords: Multitemporal denoising, nonlocal means (NLM), synthetic aperture radar (SAR) images

This paper presents a denoising approach for multitemporal synthetic aperture radar (SAR) images based on the concept of nonlocal means (NLM). It exploits the information redundancy existing in multitemporal images by a two-step strategy. The first step realizes a nonlocal weighted estimation driven by the redundancy in time, whereas the second step makes use of the nonlocal estimation in space. Using patch similarity missregistration estimation, we also adapted this approach to the case of unregistered SAR images. The experiments illustrate the efficiency of the proposed method to denoise multitemporal images while preserving new information.

I. INTRODUCTION

W ITH increasing number of synthetic aperture radar (SAR) systems, many SAR images of the same region are nowadays available. Those time-series or multitemporal SAR images provide users more comprehensive information, including both spatial and temporal domains. There are different applications using multitemporal series: denoising and estimation [START_REF] Bruniquel | Multi-variate optimal speckle reduction in SAR imagery[END_REF]- [START_REF] Sen | A homomorphic system to reduce speckle in videos[END_REF] and change detection [START_REF] Bovolo | A detail-preserving scale-driven approach to change detection in multitemporal SAR images[END_REF], [START_REF] Inglada | A new statistical similarity measure for change detection in multitemporal SAR images and its extension to multiscale change analysis[END_REF] and classification [START_REF] Hachicha | Multi-temporal SAR classification according to change detection operators[END_REF], [START_REF] Bruzzone | An advanced system for the automatic classification of multitemporal SAR images[END_REF]. In this paper, we will focus on multitemporal denoising with the purpose of exploiting available information existing in the temporal series as much as possible.

Beginning from decades ago, SAR image denoising (or despeckling) has been a well-studied problem in the image processing and remote sensing community. State-of-the-art algorithms such as those by Lee [START_REF] Lee | Digital image enhancement and noise filtering by use of local statistics[END_REF], [START_REF] Lee | Refined filtering of image noise using local statistics[END_REF], Kuan et al. [START_REF] Kuan | Adaptive noise smoothing filter for images with signal-dependent noise[END_REF], [START_REF] Kuan | Adaptive restoration of images with speckle[END_REF], and Frost et al. [START_REF] Frost | A model for radar images and its application to adaptive digital filtering of multiplicative noise[END_REF] consider SAR despeckling based on estimation theory. They develop denoising techniques under minimum mean square error or maximum a posteriori estimation theory. Touzi [START_REF] Touzi | A review of speckle filtering in the context of estimation theory[END_REF] proposed a structural multiresolution framework to handle stationary and nonstationary signals and improve the previous approaches by context adaptation.

Another family of approaches based on Wavelet can be readily extended to SAR despeckling after a logarithm transfor-mation. Typical successful examples are [START_REF] Guo | Wavelet based speckle reduction with application to SAR based ATD/R[END_REF] and [START_REF] Gagnon | Speckle filtering of SAR images: A comparative study between complex-wavelet-based and standard filters[END_REF]. Considering the heavy-tailed feature of speckle distribution, a heavytailed model under Bayesian wavelet shrinkage is proposed in [START_REF] Achim | SAR image denoising via Bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF]. Recently, based on spatially adaptive Wavelet analysis, Ranjani and Thiruvengadam [START_REF] Ranjani | Dual-tree complex wavelet transform based SAR despeckling using interscale dependence[END_REF] and Gleich et al. [START_REF] Gleich | Despeckling of TerraSAR-X data using second-generation wavelets[END_REF] used spatial information to better preserve image details, such as edges and textures.

In 2005, Buades et al. [START_REF] Buades | A non-local algorithm for image denoising[END_REF] presented a nonlocal means (NLM) approach, which uses image redundancy to estimate each pixel. The idea is to search similar pixels and average them with different weights expressing the similarity and obtained by patch correlation. Taking inspiration from NLM, Dabov et al. [START_REF] Dabov | Image denoising by sparse 3-D transform-domain collaborative filtering[END_REF] proposed block-matching denoising based on 3-D transform (BM3D). With collaborative filtering as the core idea, it works on nonlocal regions (the group of patches collected from the whole image based on their similarity). Katkovnik et al. [START_REF] Katkovnik | From local kernel to nonlocal multiple-model image denoising[END_REF] proposed an extension of BM3D with pointwise-based shape-adaptive approach generalizing the advantages of both the BM3D filter and the Pointwise filter [START_REF] Foi | Pointwise shape-adaptive DCT for high-quality denoising and deblocking of grayscale and color images[END_REF].

Deledalle et al. [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF] extended the NLM to SAR image denoising. From a theoretical point of view, this probabilistic patch-based (PPB) algorithm [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF] is a more general form of the original NLM in [START_REF] Buades | A non-local algorithm for image denoising[END_REF]. By considering the distribution of SAR data with the Goodman model [START_REF] Goodman | Some fundamental properties of speckle[END_REF], they suggest using a generalized likelihood ratio (GLR) to measure the similarity between SAR noisy patches. Analogously, Parrilli et al. [START_REF] Parrilli | A nonlocal SAR image denoising algorithm based on LLMMSE wavelet shrinkage[END_REF] extended the BM3D for SAR image denoising and proposed a SAR-oriented version of the BM3D filter. It applies the similarity criterion proposed in [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF] to group patches and then uses the local linear minimum-mean-square-error estimation.

Most studies on SAR image denoising are dealing with a single image. However, many temporal denoising approaches have been also explored, the simplest one being multilooking. Trouvé et al. [START_REF] Trouvé | Statistical and operational performance assessment of multitemporal SAR image filtering[END_REF] compared three filters of a multitemporal SAR image: a texture compensation multichannel filter [START_REF] Bruniquel | Multi-variate optimal speckle reduction in SAR imagery[END_REF] based on the Kuan filter, a time-space filter [START_REF] Coltuc | Time-space filtering of multitemporal SAR images[END_REF], and an extension of the maximum homogeneous region approach to multitemporal images [START_REF] Ciuc | Adaptiveneighborhood speckle removal in multitemporal synthetic aperture radar images[END_REF]. Quegan and Yu [START_REF] Quegan | Filtering of multichannel SAR images[END_REF] have a deep insight of the relationship between the desired equivalent number of looks (ENL) and the number of dates in multitemporal images and proposed a multichannel filter based on the intensity covariance matrix. Those denoising approaches of multitemporal images are local methods. Based on the NLM theory, this paper presents a more effective way to exploit the temporal information for multitemporal image denoising. We propose an adaptation of the iterative weighted PPB denoising algorithm [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF]. First, a direct extension of PPB for multitemporal SAR images is analyzed to illustrate the limited efficiency of this direct extension. Based on this analysis, we propose a two-step denoising framework, a temporal averaging step (the first step), and a spatial denoising step (the second step). This denoising framework is also applied to unregistered temporal SAR images by miss-registration estimation.

This paper is organized as follows: Section II briefly introduces the NLM and its extension for SAR images (PPB). A simple comparative experiment is presented in Section III to analyze the problem of the direct extension of PPB. The proposed method is then detailed in Section IV and extended to unregistered images in Section V. This is followed by evaluation in Section VI and conclusion in Section VII.

II. NLM

This section presents two denoising methods based on the NLM approach. They deal with Gaussian and multiplicative speckle noise, respectively.

We first introduce the notations that will be intensively used in the following.

• Bold y: a patch with size of √ K × √ K. • Nonbold y: an image (y t denotes the image at time t).

• Nonbold ỹ: the temporal mean image. • Nonbold y t (i): the ith pixel in image y t (y t (i) is the center of patch y t (i)). • For the sake of simplicity, y 1 and y 2 (with center pixels y 1 and y 2 ) denote two arbitrary patches of an image, whereas y t (i) (with center pixel y t (i)) denotes the ith patch in image y t . y 1,k is the kth pixel in patch y 1 . • u, u and L, L are similarly defined as representing, respectively, the true parameters and the number of looks. û is the estimation of u.

A. NLM for Gaussian Distributed Data

We denote by y the observed image, by y(i) the noisy intensity value at pixel index i, and by û(i) the estimation of the actual pixel value u(i) (the true value that we are looking for). The NLM perform a weighted average (or weighted maximumlikelihood estimation) of the form

û(i) = 1 Z j ∈ Ω i exp (-d NLM (i, j)) • y(j) ( 1
)
where Ω i denotes the candidate pixels used for the estimation, and Z = j∈Ω i exp(-d NLM (i, j)) is a normalization. The weights exp(-d NLM (i, j)) are based on the weighted similarity between patches extracted at positions i and j, which are defined as

d NLM (i, j) = S NLM (y(i), y(j)) h (2)
where h > 0 is a filtering parameter, and y(i) and y(j) are the two K-dimensional vectors containing the values of the |K| × |K| square patches centered, respectively, at pixel positions i and j (patches will be denoted by bold letters in the following). The similarity function S NLM is chosen as the

S NLM (y 1 , y 2 ) = y 1 -y 2 2 = k ∈ K [y 1,k -y 2,k ] 2 . ( 3
)
The exponential function transforming similarity into weights is called the kernel.

The pixel similarity S NLM (y(i), y(j)) is comparing two small square patches of size K surrounding pixels at positions i and j. y 1,k is the value of pixel k in patch y 1 . This definition is based on the assumption that similar patches have similar central pixels [START_REF] Buades | A non-local algorithm for image denoising[END_REF], [START_REF] Duval | On the parameter choice for the Non-Local Means[END_REF]. Fig. 1 shows the 2-D histogram of the Euclidean distance between noise-free pixels |u(i)u(j)| and the Euclidean distance between noise-free patches u(i)u(j) measured from the house image. From this histogram, an approximate proportional relation between |u(i)u(j)| and u(i)u(j) can be found. When Gaussian noise is considered, the similarity between patches is usually defined by the Euclidean distance, as in (3).

B. NLM for Gamma Distributed Data

For the denoising of a single SAR image, Deledalle et al. [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF] extended the NLM with PPB weights. Since this approach was proposed for the denoising of a single image only, this method will be denoted as 1-PPB. In order to deal with the specific nature of speckle noise as described in [START_REF] Goodman | Some fundamental properties of speckle[END_REF], they suggest using a similarity criterion based on the GLR for patches in SAR images (more details in [START_REF] Deledalle | How to compare noisy patches? Patch similarity beyond Gaussian noise[END_REF]), leading to

S GLR (y 1 , y 2 ) = k∈K log y 1,k y 2,k + y 2,k y 1,k (4)
where y 1 and y 2 are two patches of intensity values. In addition, a refining term is added to the weight exp(-d 1-PPB (i, j)), i.e.,

d 1-PPB (i, j) = S GLR (y(i), y(j)) h + S KL (û(i), û(j)) h ( 5
)
where h and h > 0 are two filtering parameters, and û(i) and û(j) are patches that are iteratively estimated. For the first iteration, û(i) = 1, and then, the current estimate is used for û(i). The similarity function S KL used to iteratively refine the weights is based on the symmetric Kullback-Leibler (KL) divergence defined for any pair of patches u 1 and u 2 as

S KL (u 1 , u 2 ) = k ∈ K [u 1,k -u 2,k ] 2 u 1,k u 2,k . ( 6
)
Some papers reference this approach as IT-PPB, but in this paper, 1-PPB stands for this iterative version.

C. Discussion

Compared with the local approaches, NLM has no connection restriction, since those similar pixels may be unconnected to each other. By exploiting patch similarities, the selection of similar pixels is naturally adapted to nonstationarities. For instance, for pixels located on one side of an edge, the weights are all concentrated on the same side of the transition, and the subsequent nonlocal estimation is thus unbiased, while local filters mix samples arising from both sides. A deeper comparison of NLM with local methods, such as intensitydriven adaptive neighborhood, has been done in [START_REF] Deledalle | NL-InSAR: Nonlocal interferogram estimation[END_REF]. In this paper, our contribution mainly focuses on the introduction of temporal information for denoising.

Our aim in this paper is to exploit the efficiency of nonlocal denoising for multitemporal series. We will first consider that images are finely registered and then relax this hypothesis in Section V. Our framework can be applied to interferometric series, but temporal denoising efficiency decreases with the spectrum correlation of the interferometric data.

III. DIRECT EXTENSION AND ANALYSIS

A direct extension of PPB denoising of SAR images is presented in this section. We tested it on a synthetic set of multitemporal SAR images. This comparison experiment will lead us to the approach proposed in the next section.

A. Direct Extension of PPB

The key idea of the PPB denoising in [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF] similar to the NLM is to estimate actual pixel intensity with image redundancy. The way to exploit image redundancy is the search of similar pixels in a search window and averaging those similar pixels with different weights. For multitemporal data, a direct extension of PPB can be the definition of a cube search window, by aggregating all the search windows of the different dates. Considering {y t 1 , y t 2 , . . . , y t N }, the stack of multitemporal images, the temporal NLM filter can be defined on the cube C i = Ω i × {t 1 , t 2 , . . . , t N } of all pixel indexes i. Meanwhile, the direct temporal extension of PPB (T-PPB) to estimate the true value u t (i) is

ût (i) = 1 Z (j,t ) ∈ C i exp (-d T-PPB (i t , j t )) • y t (j) (7)
where

d T-PPB (i t , j t ) = S GLR (y t (i), y t (j)) h T + S KL û t (i), û t (j) h T
and for a stack of images {y t 1 , . . . , y t N }, y t (i) and y t (i) are, respectively, the ith pixel value of y t and its surrounding patch.

B. Comparison

To analyze the interest of this direct extension for the temporal case, we have tested a synthetic set of multitemporal images {y t 1 , y t 2 , y t 3 }. These images are synthesized supposing there is no change in time (stable case), thus giving three realizations of the same scene. They have been denoised using three methods.

• Method i (1-PPB): Denoise the single image y t 1 without using images y t 2 and y t 3 by the PPB approach. • Method ii (T-PPB): Denoise y t 1 using the set of temporal images {y t 1 , y t 2 , y t 3 } by the temporal PPB approach [see [START_REF] Hachicha | Multi-temporal SAR classification according to change detection operators[END_REF]]. • Method iii (M-PPB): First, get the mean values ȳ of the set of multitemporal images {y t 1 , y t 2 , y t 3 } and then denoise this average image ȳ by using the PPB approach. Note that method i uses only one image, whereas method ii uses the whole image set (three images). Although methods ii and iii share the same input noisy images, method iii first temporally averages the three temporal images. Because the denoising approach of the three methods is PPB and the main difference among them is the input, one can easily predict that method i has the poorest results and that methods ii and iii should have comparable performance.

Fig. 2 shows the comparison results. Indeed, by using three 1-look noisy images, the denoising performance of T-PPB in Fig. 2(d) has improved compared with the result of 1-PPB in Fig. 2(b) (seeing the SNR values). However, it is not sufficient, because M-PPB in Fig. 2(f) has significantly outperformed T-PPB, in spite of the same input noisy images.

C. Analysis

In order to analyze these results, we represent the map of weights w(i t , j t ). Fig. 3 shows the weight maps of several interesting pixels for the three methods. In the noise-free image, we can easily find the similar pixels. However, in method ii, more dissimilar pixels have been assigned large weights (bright points in the weight map) than those in method iii. More weights are computed, but they are less accurate. Due to the temporal average before denoising, method iii reduces the risk in searching similar pixels.

A first solution could be the modification of the kernel that transforms similarity in weights. Indeed, it is shown in [START_REF] Duval | A bias-variance approach for the nonlocal means[END_REF] that truncated weights for NLM (instead of exponential function) could improve the selection. However, for a 1-look image, the improvement is not sufficient.

Averaging in the temporal domain (multilooking) is considered as the efficient and best unbiased estimator for the hypotheses of independent and identical distribution. Thus, method iii naturally performs better than methods i and ii. The purpose of this experiment is to illustrate the improvement of the introduction of multilooking into PPB. In this paper, we combine the advantages of both approaches (multilooking and PPB) to achieve a multitemporal denoising method.

IV. TWO-STEP MULTITEMPORAL NLM

Let us go back to the simple comparison experiment in Section III-B. Method iii [in Fig. 2(e) and (f)] first temporally averages noisy images, which reduces the noise level and improves the weights, as shown in Fig. 3. Although the real SAR images or video data have problems such as registration (or motion) and temporal changes, there are still lots of cases that can be temporally combined. Taking inspiration from method iii, we divide the denoising process into two steps, which deal with temporal and spatial information, respectively. The main purpose of multitemporal denoising is to exploit all available information for temporally stable pixels while keeping the new information as much as possible. Note that in this part, all multiimages are well registered, and only the change detection problem is considered in this section.

In case of stable pixels (no change over time and well registered), method iii illustrates the usefulness of temporal averaging. The proposed denoising framework, namely, a twostep (2S-PPB) approach based on NLM and PPB, is to exploit similar pixels in the temporally average image ȳ rather than in the stack {y t 1 , . . . , y t N } (the premise is that there is no change taking place among multiimages). If pixels located at the same position but at different times (such as y t (i) and y t (i)) have not changed (in other words, they share a same true value u t (i) = u t (i)), they can be averaged together to estimate u t (i). From a probability point of view, this equally weighted average can be considered as an estimation using prior information (relative to the estimate using likelihood information in [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF]). This temporal average can be seen as a preliminary multilooking operation. However, the main problem lying in this temporal denoising is the use of only stable pixels. Thus, we have to use change detection, which will be detailed in the following.

A. Change Criteria

Most change detection methods are concerned about changes between different terrestrial objects, such as rivers, buildings, and other artificial objects. However, in image denoising, we theoretically combine the observed pixels from the same true value or reflectivity (coming from the same distribution). Hence, the changes in this section are defined as samples coming from different distributions.

Several similarity criteria for noisy patches are compared in [START_REF] Deledalle | How to compare noisy patches? Patch similarity beyond Gaussian noise[END_REF]. The GLR is proposed for non-Gaussian noise. This criterion can be used to detect changes in times. Another criterion, i.e., the symmetric KL divergence criterion (hereafter referred to as the KL criterion) in denoised images (using 1-PPB on each image), is proposed in [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF] to improve the performance of the GLR criterion computed in noisy images. Both criteria, GLR on noisy images and KL divergence on denoised data, have been investigated for similar pixel selection.

Fig. 4 illustrates the performance of the GLR criterion in the noisy image and the KL criterion in the denoised image (using 1-PPB on each image) to find the similar pixels. The results show that the KL criterion in denoised images is more capable of selecting similar pixels than the GLR criterion in the noisy image. However, this comparison experiment is not fair, because of the different noise levels in Fig. 4(a) and (c). What is more, the KL criterion in denoised images has an inevitable drawback that the used denoising approaches have great influence on the quality of the selections.

To balance the advantages against the disadvantages of both the GLR in noisy images and the KL criterion in denoised images, we suggest employing both of them to detect the changes between temporal images (see Fig. 5). This suggestion has the similar consideration as in [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF], which proposed that the denoising weights are coming from both the noisy images (by the GLR criterion) and the previous denoised images (by KL divergence).

For the sake of simplicity, we use (8) as a binary criterion to define the temporal relation between the pixel values at i t and the pixel value at j t as

ϕ [d 1st (i t , j t )] = 1, if d 1st (i t , j t ) < 2 0, otherwise (8)
d 1st (i t , j t ) = S GLR (y t (i), y t (j)) h 1st + S KL (û t (i), ût (j)) h 1st
where ût (i) and ût (j) are patches extracted in the denoised result, respectively, in û t and û t using 1-PPB. In order to ensure that both GLR and KL criteria have the same contribution on change detection, we normalize the GLR and KL terms in (8) by parameters h 1st = quantile(S GLR , α) and h 1st = quantile(S KL , α), respectively, where for any similarity measure f , quantile(f, α) denotes the α-quantile of the pure distribution of f (i.e., the distribution when patches have the same underlying reflectivity). In practice, we have chosen α = 0.99, such that a binary weight ϕ[d 1st (i t , j t )] = 1 means that pixel y t (i) and y t (i) have a high probability to be realizations coming from the same underlying reflectivity. If α keeps constant, a threshold larger than "2" will add more changed pixels in the temporal step. On the contrary, a smaller threshold may neglect some similar (unchanged) pixels in the temporal step.

B. Two-Step Denoising

The following formula shows the temporal averaging process with the change criterion (binary weights):

ỹt (i) = 1 Z t ∈ [t 1 ,t N ] ϕ [d 1st (i t , i t )] • y t (i) (9)
where

Z = t ∈[t 1 ,t N ] ϕ[d 1st (i t , i t )].
Denoising on the temporal estimate ỹt (i) is comparable to [START_REF] Hachicha | Multi-temporal SAR classification according to change detection operators[END_REF]. However, different pixels ỹt (i) may have different (equivalent) number of looks depending on the number of averaged data. The resulting number of looks is temporally and spatially varying and is given by [START_REF] Lee | Refined filtering of image noise using local statistics[END_REF] in which L t is the original spatially invariant (equivalent) number of looks of the image y t . In this case, the similarity between patches ỹt (i) and ỹt (j) has to be modified to take into account the spatially varying number of looks. For any patches ỹ1 and ỹ2 , the GLR gives (see Appendix A)

Lt (i) = t ∈ [t 1 ,t N ] ϕ [d 1st (i t , i t )] • L t
S GLR (ỹ 1 , ỹ2 ) = k ∈ K ( L1,k + L2,k ) log[ L1,k ỹ1,k + L2,k ỹ2,k ] -( L1,k + L2,k ) log( L1,k + L2,k ) -L1,k log[ỹ 1,k ] -L2,k log[ỹ 2,k ] (11)
where L1 and L2 are, respectively, the (equivalent) number of looks of patches ỹ1 and ỹ2 . Then, the estimate ût (i) using ỹt will be

ût (i) = 1 Z j ∈ Ω i exp (-d 2nd (i t , j t )) • ỹt (j) (12)
where

d 2nd (i t , j t ) = S GLR (ỹ t (i), ỹt (j)) h 2nd + S KL û t (i), û t (j) h 2nd
and S KL plays the same role as in [START_REF] Hachicha | Multi-temporal SAR classification according to change detection operators[END_REF] for iteratively refining the weights by using patches û(i) and û(j) extracted from the previous estimate û t . However, by considering the different number of looks, the similarity function S KL has also to be modified according to the KL divergence, where for any pair of patches u 1 and u 2 (see Appendix B for details)

S KL (ũ 1 , ũ2 ) = k ∈ K L1,k ũ2,k ũ1,k + L2,k ũ1,k ũ2,k -L1,k -L2,k +( L1,k -L2,k ) ψ( L1,k ) -ψ( L2,k ) + ln(ũ 1,k ) -ln(ũ 2,k ) . ( 13
)
Algorithm 1 summarizes the steps to denoise temporal SAR images by the proposed 2S-PPB method.

Algorithm 1

The proposed 2S-PPB algorithm.

Input:

Well

(i) = ỹt 1 (i) + ϕ[d 1st (i t 1 , i t )] • y t (i); [see (9)] 12: Lt 1 (i) = Lt 1 (i) + ϕ[d 1st (i t 1 , i t )] • L t 13:
end for 14:

ỹt

1 (i) = ỹt 1 (i)/ Lt 1 (i); [see (9)

] 15: end for ----------------------------------Step 2 (Spatial step):

-----16: denoise ỹt 1 (i) using ( 11), [START_REF] Kuan | Adaptive restoration of images with speckle[END_REF], and ( 13) in Section IV-B. This denoising step is similar to the PPB approach.

-----------------------------17: Denoised result ût 1 ;

V. TWO-STEP MULTITEMPORAL NLM FOR UNREGISTERED SAR IMAGES

The two-step multitemporal NLM proposed in the previous section deals with multitemporal SAR images that have been well registered. However, in practice, it is not always easy to get well-registered images because of the lack of accuracy of sensor parameters or of terrain deviation. Here, we propose an adaptation of the two-step multitemporal NLM to deal with unregistered temporal SAR images.

A. Miss-Registration Estimation

From the simple comparison experiment in Section III-B, we have seen that the direct extension of 2-D denoising approaches to well-registered temporal SAR images is not optimal. Hence, unlike in [START_REF] Buades | Nonlocal image and movie denoising[END_REF], we try to consider the offset caused by missregistration between temporal SAR images using patch similarity. This offset between remote sensing SAR images caused by miss-registration is much simpler than the complex scene changes in video. Moreover, our aim is not accurate registration but only accurate detection of similar pixels to perform temporal averaging.

Let y t and y t denote two temporal images without registration, and y t (i) and y t (i) are the ith pixel in y t and y t , respectively. Note that y t (i) and y t (i) are both located at i in Fig. 6. Assumption of image registration. In the intensity image registration, all pixels in a small patch share the same offset to their real position.

images, but they do not denote the same position in the geographic coordinate system before image registration. Suppose that y t (i + v tt (i)) is the [i + v tt (i)]th pixel in y t that shares the same geographic location with y t (i) in y t , and v tt (i) denotes the offset between y t (i) and y t (i + v tt (i)). If no change takes place, y t (i + v tt (i)) should be the most similar pixel of y t (i) in image y t . Consequently, the idea to estimate this offset v tt (i) caused by miss-registration is to use the pixel similarity.

Using the conclusion in Section IV-A, we combine the KL divergence and the GLR criterion as the similarity between patches, illustrated in [START_REF] Bruzzone | An advanced system for the automatic classification of multitemporal SAR images[END_REF]. (Note that this KL divergence criterion is from the 1-PPB denoising results on each image.) For each pixel y t (i) in y t , we measure the similarity between y t (i) and the pixels in region A t (i) (an image patch with y t (i) as the center in y t , shown in Fig. 7). The similarity between y t (i) and y t (j) is d 1st (i t , j t ). Each pair of pixel y t (i + n) in y t and region A t (i + n) in y t has a similarity map M tt (i + n) (0 n N A -1, N A is the size of region A), as shown in the top-right image in Fig. 7. Thus M tt (i) = {d 1st (i t , j t )} j t ∈A t (i) .

(

) 14
The offset v tt (i) should be the offset between the center of map M tt (i) and the position of minimum value in M tt (i). However, it leads to a poor performance because of noise. Based on the assumption in Fig. 6, we suggest that the offset v tt (i) caused by miss-registration can be estimated from an average similarity map Mtt (i), which is the mean of similarity maps M tt (i + n) (1 n N A ). The estimation of offset v tt (i) is from the center of map Mtt (i) to the position of minimum value in Mtt (i), i.e.,

Mtt (i) = d1st (i t , j t ) j t ∈A t (i) d1st (i t , j t ) = i t ∈A t (i) d 1st (i t , j t ). (15)
We test the miss-registration estimation on real SAR images. Fig. 14(b) is the temporal average without image registration or miss-registration estimation, which is blurred. After missregistration estimation, the temporal average is illustrated in Fig. 14(c).

B. Two-Step Denoising

The miss-registration estimation entitles the proposed 2S-PPB denoising approach to deal with unregistered multitem- poral images. The denoising process is exactly the same as detailed in Section IV, except that the candidate pixels y t (i) are replaced by y t (i + v tt (i)) found out by the miss-registration estimation (see Fig. 7). Algorithm 2 summarizes the steps to denoise unregistered temporal SAR images by the proposed 2S-PPB method.

Algorithm 2 The proposed 2S-PPB algorithm (unregistered temporal images).

Input:

Unregistered temporal SAR images {y t ----------------------------------Step 1 (Temporal step):-----15: compute ỹt 1 using ( 8) and ( 9) by taking into account the vector field v t 1 t . This denoising step is similar to the one in Algorithm 1 where i t is substituted to

i t + v t 1 t (i).
- ---------------------------- -----Step 2 (Spatial step):-----16: denoise ỹt 1 (i) using ( 11), [START_REF] Kuan | Adaptive restoration of images with speckle[END_REF], and ( 13) in Section IV-B. This denoising step is similar to the PPB approach.

- ---------------------------

A. Test of Pixel Similarity

The value of d 1st given by ( 8) is a key point of our approach, since it is used both for the change detection step and for the offset estimation step. We compare our change detection criterion with the log-ratio, which is a criterion widely used in SAR images. Synthetic images and real SAR images are used in these comparison experiments: 1) synthetic SAR images corrupted by a multiplicative 1-look speckle noise in Fig. 8(a); 2) two real SAR images (TerraSAR images in Saint-Gervaisles-Bains, France) sensed in 2009 and 2011, respectively, in Fig. 8(b). Fig. 8(a) and(b) shows noisy images at different times. The reference maps of changes are shown in Fig. 8(c) (for the real SAR images, we manually label the changes taking place between date 1 and date 2). Fig. 8(e) shows our similaritybased change detection criterion, and Fig. 8(d) shows the logratio criterion. Our similarity-based change detection criterion has higher receiver operating characteristic curves in the falsepositive and true-positive curves in Fig. 8(f).

We also compare our offset estimation with the offset estimation using the intensity tracking method. Intensity tracking with correlation has been widely used in glacier motion estimation tasks on multitemporal SAR images [START_REF] Gray | Velocities and flux of the Filchner ice shelf and its tributaries determined from speckle tracking interferometry[END_REF], [START_REF] Strozzi | Glacier motion estimation using SAR offset-tracking procedures[END_REF]. In this comparison experiment, we simulate a multitemporal image set with an offset v tt = (3, 3) pixels between them, as shown in Fig. 9(a) and(b).

For the offset estimation using intensity tracking, it is selfimplemented in our comparison experiments, and it also has a local averaging process similar to that in the proposed method to reduce the effect of speckle noise. The patch size is 7 × 7 pixels, and the test area is 21 × 21 pixels. To illustrate the performances, the estimation error is calculated using the Euclidean distance between the true offset v tt = (3, 3) and the estimated offset. Fig. 9(c) and(d) shows that the proposed offset estimation has less estimation error than the log-ratio similarity.

B. Denoising of Well-Registered Images

For all the experiments, we use the parameters as suggested by Deledalle et al. in [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF]. The search window Ω and patch size K enlarge with the increase in the number of iterations, Ω ∈ {3 × 3, 7 × 7, 11 × 11, 21 × 21} and K ∈ {1 × 1, 3 × 3, 5 × 5, 7 × 7} for all experiments. h 2nd , controlling the decay of weights, depends on the distribution of similarity between two noisy patches with the same underlying reflectivity and is defined as the α-quantile of this pure distribution (α = 0.92). h 2nd = 0.2|K|.

The experiments are taken under the MATLAB environment on an Intel Core-2 Quad CPU Q9550 at 2.83-GHz 64-bit computer. The time consumption of the proposed method is about 400 s for a three-date 256 × 256 temporal image set (PPB needs about 15 s for a 256 × 256 image).

1) Synthetic Images: We present visual and numerical results obtained on synthetic images corrupted by multiplicative Goodman's speckle noise. The classical noise-free images are used: house, lena, barbara, boat, and peppers. We use the same noise-free image to synthesize a temporal image set, which means that there are no temporal changes. In addition to the proposed 2S-PPB filter, the comparisons that have been tested here are PPB only on single image y t 1 (method i: 1-PPB), 3-D adaptive neighborhood filter (3-D-ANF, it is self-implemented) [START_REF] Ciuc | Adaptiveneighborhood speckle removal in multitemporal synthetic aperture radar images[END_REF] on the multitemporal image set {y t 1 , . . . , y t N }, and temporal PPB on the multitemporal image set {y t 1 , . . . , y t N } (method ii: T-PPB). For the 3-D-ANF filter, we use a 3 × 3 median filter in the first step of 3-D-ANF and 50 × 50 as the maximum size of adaptive neighborhood in the last step of 3D-ANF. Fig. 10 only shows the denoising results of the house images corrupted by L = 1 multiplicative speckle noise. There are three noisy images in the temporal data set. The image obtained by the T-PPB filter is well smoothed compared with the 1-PPB filter. However, the edge and shape preservation has limited improvement. The proposed 2S-PPB filter provides more details of edges as the eaves and windows of house, whereas smooth regions are comparable. 3D-ANF has less loss of structural information in stable cases [shown in the ratio map in Fig. 10(b)], whereas it has poor noise reduction in homogeneous regions.

To quantify the denoising qualities, Table I presents numerical results for images corrupted by multiplicative speckle noise with different ENL L = 1, 3, 5, and 10 and different number of dates (different number of images in the temporal data set) N = 1, 2, 3, and 5. Note that when N = 1, there are no T-PPB and 2S-PPB denoising results. The used performance criterion is the signal-to-noise ratio (SNR), i.e., SNR(û t , u t ) = 10 log 10

Var

[u t ] Mean [(û t -u t ) 2 ] . ( 16
)
We observe that the 2S-PPB filter improves on the T-PPB filter for high-noise-level images (i.e., L = 1), particularly when N is large. However, for low-noise-level images (i.e., L > 5), 2S-PPB has only limited improvement. This is because the similarity of noisy patches in low-noise-level images is efficient enough, and the improvement of patch similarity using temporal average is relatively less important.

2) Realistic SAR Synthetic Images: This part presents denoising results of realistic SAR synthetic images. It is a 100-look SAR acquisition identified as Toulouse of the CNES in Toulouse (France) sensed by RAMSES and provided by the CNES. We corrupt this 100-look Toulouse image by 1-look multiplicative speckle noise to form three temporal images {y t 1 , y t 2 , y t 3 }. In order to simulate the changes in the multitemporal images, a dark line and a bright target are added to y t 1 , labeled by red rectangles (Region #1 and #2) in Fig. 12 Fig. 12 presents the obtained denoising results for the Toulouse images. The results of the proposed 2S-PPB filter have better edge and shape preservation with comparable smoothed regions than 1-PPB and T-PPB filters (ENL in Table II). Structural information of new objects (changes) is also better preserved, as shown in Fig. 12(a) (Region #1 and #2). Stable objects, as the dark lines in Region #3 and #4, have been better restored in the results of 2S-PPB than 1-PPB and T-PPB filters.

3) Real Multitemporal SAR Images: We also test on real multitemporal SAR data as follows:

• San Francisco (well-registered) IGARSS: six single-look TerraSAR images of San Francisco (USA) provided by San Francisco (well-registered). Moreover, more dark and thin streets have been preserved. We also compare the 2S-PPB filter with the classical temporal filtering method, multilooking approach in Saint-Gervais-les-Bains (well-registered). Generally, the 2S-PPB filter gets more smoothed results than multilooking, as Region #1 in Fig. 13(c) and (d). More than that, the changes over time can be restored by 2S-PPB, as the changes in Region #2 in Fig. 13(c) and (d). Multilooking tends to ignore the temporal changes (loss of temporal resolution), whereas the 2S-PPB filter can well preserve both the spatial and temporal resolution and reduce the speckle effect. Table II shows the ENL calculated in homogeneous regions (green rectangles in Fig. 11).

In general, the proposed filter has higher ENL than other filters.

C. Denoising of Unregistered Images

The real temporal SAR images provided by IGARSS (data information in Section VI-B3) are tested here. The difference is that we manually register those temporal images without using the accurate sensor parameters (see Figs. 12 and13). Moreover, the miss-registration between two images is about 4-7 pixels (shown in Fig. 14(b), the temporal mean of the temporal images illustrates the unregistration). These unregistered temporal images are identified as San Francisco (unregistered). In the miss-registration estimation, the search region A t (i) is a 21 × 21 pixel window. Patch size is 7 × 7 pixels, h and h have been chosen identical to the filtering step. After missregistration estimation, the proposed 2S-PPB filter is applied with the same parameters used in Section VI-B.

The temporal mean after miss-registration estimation in Fig. 14(c) shows the performance of the miss-registration estimation. Similarly, for the thin and dark streets in San Francisco (unregistered), the proposed 2S-PPB preserves more details than 1-PPB and T-PPB. However, its performance is not as good as the one on well-registered San Francisco [in Fig. 14(f)], because of the insufficient estimation of miss-registration. In addition, the miss-registration estimation will fail when the offset between temporal images is too large. Indeed, larger offset estimation needs larger search region A t (i), but this increases the risk of finding the similar pixels.

VII. CONCLUSION AND DISCUSSION

In this paper, we have presented the limitation of spatiotemporal similarity. Inspired by this phenomenon, we have proposed a two-step denoising framework based on iterative weighted PPB. With the help of miss-registration estimation, we also adapted it to unregistered images. The effective performance in the experiments shows that the proposed method exploits more available information for stable objects while comparably keeping new objects.

Future work will be focused on finding a more efficient change detection method for temporal average. In addition, the adoption of the proposed 2S-PPB for interferometric SAR, polarimetric SAR, and polarimetric interferometric SAR images will be another subject of future work.

APPENDIX A GLR WITH DIFFERENT NUMBER OF LOOKS

We denote by p(y, u, L) the Gamma distribution, i.e., p(y,u,L

) = 1 Γ(L) L u Ly u L-1 exp - Ly u ( 18
)
where y is the intensity value, u is the noise-free value, and L is the number of looks. The GLR similarity criterion [START_REF] Deledalle | How to compare noisy patches? Patch similarity beyond Gaussian noise[END_REF] between y 1 and y 2 is Note that in [START_REF] Deledalle | How to compare noisy patches? Patch similarity beyond Gaussian noise[END_REF], the number of looks L is the same (L 1 = L 2 = L). Thus, the parameters of distribution θ only consider the noise-free value u. However, in our case, the number of looks is different, which will have influence on the similarity criterion.

S GLR (y 1 , y 2 ) = -log p(y 1 , u 1 = û12 , L = L 1 )p(y 2 , u 2 = û12 , L = L 2 ) p(y 1 , u 1 = û1 , L = L 1 )p(y 2 , u 2 = û2 , L = L 2 ) .
Substitute the Gamma distribution p(y, u, L) in the GLR criterion with û1 = y 1 , û2 = y 2 , and û12 = (L 1 y

1 + L 2 y 2 /L 1 + L 2 ), i.e., S GLR (y 1 , y 2 ) = log y 1 û12 û1 û12 L 1 -1 exp L 1 - L 1 y 1 û12 -1 + log y 2 û12 û2 û12 L 2 -1 exp L 2 - L 2 y 2 û12 -1 = L 1 log L 1 y 1 + L 2 y 2 y 1 (L 1 + L 2 ) + L 2 log L 1 y 1 + L 2 y 2 y 2 (L 1 + L 2 ) . When L 1 = L 2 = L, we have S GLR (y 1 , y 2 )| L 1 =L 2 =L = 2L log y 1 y 2 + y 2 y 1 -2L log 2
which is the same as the weights formula in [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF] (y is the intensity value, and √ y is the amplitude value).

APPENDIX B SYMMETRIC KL DIVERGENCE WITH DIFFERENT NUMBER OF LOOKS

For distributions p(y, u 1 , L 1 ) and p(y, u 2 , L 2 ) of a continuous random variable, KL divergence is defined to be the integral the KL divergence criterion is

S KL (u 1 , u 2 ) = D KL,1 + D KL,2 = L 1 u 2 u 1 + L 2 u 1 u 2 -L 1 -L 2 + (L 1 -L 2 ) [ψ(L 1 ) -ψ(L 2 )] + (L 1 -L 2 ) [ln(u 1 ) -ln(u 2 )] . ( 26
)
When L 1 = L 2 = L, we have

S KL (u 1 , u 2 ) = D KL,1 + D KL,2 = L u 1 u 2 + u 2 u 1 -2 (27)
which is the same as the weights formula in [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF].

Xin Su is currently working toward the Ph.D. degree at Télécom ParisTech, Paris, France. His main research interests include multitemporal synthetic aperture radar image analysis, signal processing, image processing, pattern recognition, and machine learning.

Fig. 1 .

 1 Fig. 1. Two-dimensional histogram of the Euclidean distance between noisefree patches K and noise-free pixels. The ideal case (the assumption that similar patches have similar central pixels [29]) is shown on the white line. (a) Test image house. (b) Two-dimensional histogram of the Euclidean distance. square of the Euclidean distance defined for any pair of patches y 1 and y 2 by

Fig. 2 .

 2 Fig. 2. Denoising results of synthetic images with multiplicative speckle noise (left: input noisy images; right: denoising results). (a) 1-look noisy image. (b) Method i 1-PPB result of (a); SNR: 10.79. (c) Three 1-look noisy multitemporal images. (d) Method ii T-PPB result of (c); SNR: 11.62. (e) Temporal mean of (c) (a 3-look noisy image). (f) Method iii M-PPB result of (e); SNR: 13.81.

Fig. 3 .

 3 Fig. 3. Weight maps. (a) Three 1-look noisy images. (b) Temporal mean of (a). (c), (d), and (e) Weight maps of y(1), y(2), and y(3) in (a). (f), (g), and (h) Weight maps of y(1), y(2), and y(3) in (b).

Fig. 4 .

 4 Fig. 4. Selected similar pixels by the GLR criterion and the KL divergence criterion. (a) The 3-look noisy house image and the pixel that will be used to illustrate the distributions; the red rectangle is the search window. (b) Selections by the GLR criterion. (c) PPB denoising result of (a). (d) Selections by the KL criterion. The KL criterion in denoised images is more capable of selecting similar pixels than the GLR criterion in noisy images.

Fig. 5 .

 5 Fig. 5. Change detection by the KL divergence and the GLR criterion. (a) and (b) Noisy house images with 1-look speckle noise. (c) Change detection results between (a) and (b). (d) and (e) Noisy SAR images; the red rectangles denote the changes between (d) and (e). (f) Change detection results between (d) and (e). Changes in red rectangles have been found out in detection results.

Fig. 7 .

 7 Fig. 7. Sketch map of the miss-registration estimation for unregistered images.

Fig. 8 .

 8 Fig. 8. Change detection criterion comparison. (a) Noisy image at date t. (b) Noisy image at date t . (c) Reference map of changes. (d) Log-ratio change criterion. (e) Similarity map. (f) Change detection results (true-positive and false-positive curves).

- 17 :

 17 return Denoised result ût 1 ; VI. EXPERIMENTS Here, we analyze the performance of the proposed approaches. Before denoising evaluation, the similarity criterion combining GLR and KL criteria used for change detection or offset estimation is studied.

Fig. 9 .

 9 Fig. 9. Offset estimation. (a) Noisy image 1. (b) Noisy image 2 (with an offset v tt = (3, 3) with image 1). (c) Estimation error of intensity tracking [34], [35] (the distance between intensity tracking offset estimation results and the true offset v tt = (3, 3)). (d) Estimation error of the proposed similarity (the distance between the proposed offset estimation results and the true offset v tt = (3, 3)).

  (a). y t 2 and y t 3 are corrupted by different multiplicative speckle noise without the dark line and the bright target [see Fig. 12(e) and (f)].

Fig. 10 .

 10 Fig. 10. Denoising results for the house images corrupted by multiplicative speckle noise with L = 1. From left to right: the denoising results and the ratio of the noisy image to the denoised image. (a) 1-PPB filter. (b) 3D-ANF. (c) T-PPB filter. (d) Proposed 2S-PPB filter. TABLE I SNR VALUE OF ESTIMATED IMAGES USING 1-PPB (BROWN), 3D-ANF (BLUE), T-PPB (GREEN), AND 2S-PPB (RED) FILTER FOR IMAGES CORRUPTED BY MULTIPLICATIVE SPECKLE NOISE WITH DIFFERENT ENL L = 1, 3, 5, 10 AND DIFFERENT NUMBERS OF IMAGES N = 1, 2, 3, 5 IN TEMPORAL DATA SET

Fig. 11 .

 11 Fig. 11. (a) Tozulouse, (b) San Francisco, and (c) Saint-Gervais-les-Bains noisy images (regions in green rectangles are used to calculated the ENL values).

Fig. 12 .

 12 Fig. 12. Denoising results on a zoom of Toulouse DGA ONERA. From left to right: noisy image with 1-look multiplicative speckle noise, results by PPB on a single image (1-PPB), results by 3D-ANF on temporal images, results by T-PPB on temporal images, and results by 2S-PPB on temporal images. From top to bottom: image yt 1 with new objects (a dark line in Region #1 and a bright target in Region #2), ratio maps of noisy images to denoised images, image yt 2 without new objects, and image yt 3 without new objects.

Fig. 13 .

 13 Fig. 13. Denoising results of well-registered San Francisco IGARSS and Saint-Gervais-les-Bains. (a) San Francisco denoising results. From left to right: 1-PPB on a single image, 3D-ANF, T-PPB, and 2S-PPB. (b) Ratio map of San Francisco images. From left to right: 1-PPB on a single image, 3D-ANF, T-PPB, and 2S-PPB. (c) Saint-Gervais-les-Bains images. From left to right: noisy image, denoising result by multilooking, and denoising results by 1-PPB on a single image. (d) Saint-Gervais-les-Bains images. From left to right: 3D-ANF, denoising result by T-PPB on temporal images, and denoising results by 2S-PPB on temporal images. Stable Region #1 in Saint-Gervais-les-Bains and changed Region #2 in Saint-Gervais-les-Bains.

  , u 1 , L 1 ) ln p(y, u 1 , L 1 ) p(y, u 2 , L 2 ) dy. (19)Substitute the Gamma distribution p(y, u, L) in the GLR criterion with y ≤ 0, i.e.,D KL,1 = +∞ 0 p(y, u 1 , L 1 ) ln[A • B • C]dy[START_REF] Buades | A non-local algorithm for image denoising[END_REF]

Fig. 14 .u 1 + ( 1 -u 2 + (L 1 -u 1 u 2 -

 1411212 Fig. 14. Denoising results of unregistered San Francisco IGARSS. (a) y t1 , one of the noisy images (1-look). (b) Temporal average of the multitemporal SAR images, which shows that the temporal images are unregistered. (c) Temporal average of the multitemporal SAR images after miss-registration estimation. (d) Denoising results by T-PPB on unregistered temporal images. (e) Denoising results by 2S-PPB on unregistered temporal images. (f) Denoising results by 2S-PPB on well-registered temporal images.where

  Charles-Alban Deledalle received the B.S. degree in engineering from the École Pour l'Informatique et les Techniques Avancées, Kremlin Bicêtre Cedex, France; the M.S. degree in science and technology from the University Pierre et Marie Curie (Paris 6), Paris, France, in 2008; and the Ph.D. degree at Télécom ParisTech, Paris, in 2011. He is currently a Research Associate with the Institut de Mathématiques de Bordeaux, Centre National de la Recherche Scientifique-Université Bordeaux 1, Talence, France. His main research interests include image denoising, analysis, and risk estimation for inverse problems. Dr. Deledalle received the Best Student Paper Award at the 17th IEEE International Conference on Image Processing (ICIP 2010) and a Ph.D. award in Signal, Image and Vision (Club EEA/GdR ISIS/GRETSI) for his thesis.

  -registered temporal SAR images {y t 1 , y t 2 , . . . , y t N }. A date t 1 of interest. for each y t in {y t 1 , y t 2 , . . . , y t N } do

	Output:
	ût 1 : the denoising result of image y t 1 .
	-----Step 1 (Temporal step):-----
	1: 2:	denoise y t using the 1-PPB approach;
		(see Section II-B)
	3:	obtain predenoised results û t ;
	4:	for each pixel index i do
	5:	compute change criterion ϕ[d 1st (i t 1 , i t )]; [see (8)]
	6:	end for
	7: end for
	8: for each pixel index i do
	9:	initialize ỹt 1 (i) = 0, Lt 1 (i) = 0;
	10:	for each y

t in {y t 1 , y t 2 , . . . , y t N } do 11: ỹt 1

  for each y t in {y t 1 , y t 2 , . . . , y t N } do

	2:	denoise y t using the 1-PPB approach;
		(see Section II-B)
	3:	obtain predenoised results û t ;
	4: end for
	5: for each pixel index i do
	6: 7:	for each image y t in {y t 1 , y t 2 , . . . , y t N } do Mt 1 t (i) = 0;
	8:	for each pixel index i + n in region A t 1 (i) do
	9: 10:	compute similarity map M t 1 t (i + n); [see (15)] Mt 1 t (i) = Mt 1 t (i) + M t 1 t (i + n);
	11:	end for
	12:	find the offset vector v t 1 t (i) in similarity map Mt 1 t (i);
	13:	end for
	14: end for

1 , y t 2 , . . . , y t N }. A date t 1 of interest. Output: ût 1 : the denoising results of image y t 1 .

-----Miss-Registration estimation:-----1:

TABLE II ENL

 II OF NOISY IMAGES AND ESTIMATED IMAGES USING1-PPB (BROWN), 3D-ANF (BLUE), T-PPB (GREEN), AND 2S-PPB (RED) FILTERS FOR REAL SAR IMAGES

SU et al.: TWO-STEP MULTITEMPORAL NONLOCAL MEANS FOR SYNTHETIC APERTURE RADAR IMAGES

ACKNOWLEDGMENT

The authors would like to thank J.-M. Nicolas for his help on image registration and the IGARSS Data Fusion Contest for providing the multitemporal SAR images. They would also like to thank the German Aerospace Agency (DLR) for the Terra-SAR-X images (project MTH0232) and the Centre National d'Etudes Spatiales and the Office National d'Etudes et de Recherches Aérospatiales for the 100-look SAR images.