Keywords: mercury, Antarctica, stable isotopes, food web

can reach the environment through natural and human-related sources, threatening ecosystems all over the planet due to its well-known deleterious effects. Therefore, Antarctic trophic webs, despite being relatively isolated, are not exempt of its influence. To evaluate Hg concentrations in an Antarctic ecosystem, different tissues from 2 species of invertebrates, 2 of fish, 8 of birds, 4 of pinnipeds and at least 5 of vegetation were investigated (n=176). For animals, values ranged from 0.018 to 48.7 μg g -1 dw (whole Antarctic krill and Antarctic Fur Seal liver). They were generally correlated to trophic position (assessed by δ 15 N and δ 13 C) but also to cephalopods and myctophids consumption. For vegetation, values ranged from 0.014 to 0.227 μg g -1 dw (Colobanthus quitensis and an unidentified lichen), with lichens presenting significantly higher values than mosses, likely due to year-round exposure and absorption of animal derived organic matter, as hypothesized by literature.

Mercury (Hg), besides its natural sources, can reach the environment through humanrelated activities, such as mining and the burning of fossil fuels [START_REF] Fitzgerald | The case for atmospheric mercury contamination in remote areas[END_REF][START_REF] Fitzgerald | Marine biogeochemical cycling of mercury[END_REF][START_REF] Selin | Global biogeochemical cycling of mercury: a review[END_REF]. Indeed, non-natural sources play a major role in Hg deposition, estimated to occur in a current rate three to five times higher than its pre-industrial times [START_REF] Selin | Global biogeochemical cycling of mercury: a review[END_REF].

In the open ocean the main Hg source is atmospheric deposition of inorganic Hg (II) by both wet and dry processes (e.g. [START_REF] Fitzgerald | Marine biogeochemical cycling of mercury[END_REF]. Over 80% of the Hg deposited in marine ecosystems is reemitted to the atmosphere (as Hg(0) predominantly, but some as (CH3)2Hg), increasing the residence time of Hg cycling through the reservoirs of the surface biosphere [START_REF] Driscoll | Mercury as a Global Pollutant: Sources, Pathways, and Effects[END_REF]. Once in the water column, inorganic Hg can be methylated by anaerobic microorganisms, such as sulphate-reducing and iron-reducing bacteria (e.g. [START_REF] Hsu-Kim | Mechanisms regulating mercury bioavailability for methylating microorganisms in the aquatic environment: a critical review[END_REF]. The produced Me-Hg (mono-methyl-Hg) is then easily assimilated by the biota and biomagnifies up the food web (e.g. [START_REF] Morel | The chemical cycle and bioaccumulation of mercury[END_REF]. Because microorganisms activity is importantly influenced by temperature, springtime conversion of Hg into Me-Hg may provide an important environmental pathway for its introduction into the biosphere in a time of the year when biota are preparing for peak summertime activity [START_REF] Schroeder | Arctic springtime depletion of mercury[END_REF].

In the Southern Hemisphere, high concentrations of atmospheric Hg have been observed close to the Antarctic continent when compared to lower latitudes [START_REF] Soerensen | Global concentrations of gaseous elemental mercury and reactive gaseous mercury in the marine boundary layer[END_REF].

Antarctica, in spite of being the most isolated continent on Earth, is far from exempt of the input of several sorts of contaminants, including Persistent Organic Pollutants such as PCBs or organochlorine pesticides, as reported in literature (e.g. [START_REF] Bargagli | Environmental biogeochemistry of mercury in Antarctic ecosystems[END_REF][START_REF] Corsolini | Industrial contaminants in Antarctic biota[END_REF]. Antarctica is thus contaminated by anthropogenic Hg which can affect the fitness of organisms and have consequences at the population level as shown on polar skuas [START_REF] Goutte | Demographic responses to mercury exposure in two closely-related Antarctic top predators[END_REF].

Investigating contaminants such as Hg in trophic webs can be greatly improved by using stable isotopes of carbon and nitrogen as an ecological tool to shed light, respectively, the carbon sources exploited by consumers and their trophic position [START_REF] Lesage | Diet-tissue fractionation of stable carbon and nitrogen isotopes in phocid seals[END_REF].

Previous studies on Arctic ecosystems have already shown that Hg analyses coupled to δ 15 N (e.g. [START_REF] Atwell | Biomagnification and bioaccumulation of mercury in an arctic marine food web: insights from stable nitrogen isotope analysis[END_REF] or both δ 13 C and δ 15 N (e.g. [START_REF] Aubail | Investigation of mercury concentrations in fur phocid seals using stable isotopes as tracers of trophic levels and geographical regions[END_REF] allows a better understanding of Hg biomagnification pathways and moreover, the influence of feeding ecology factors such as change in feeding sites, gender or age class trophic segregation, i.e., trophic position as a whole.

Since the late 70s, organic contaminants have received interest in King George Island [START_REF] Lukowski | DDT and its metabolites in Antarctic birds[END_REF] but with a relative lack of data for stable isotopes, with only a small quantity of studies taking SIA into account (e.g. [START_REF] Corbisier | Trophic relationships in the nearshore zone of Martel Inlet (King George Island, Antarctica): δ 13 C stable-isotope analysis[END_REF][START_REF] Majer | Bioaccumulation of potentially toxic trace elements in benthic organisms of Admiralty Bay (King George Island, Antarctica)[END_REF] for benthic invertebrates; Cipro et al., 2011 for vegetation), particularly for upper vertebrates (e.g. [START_REF] Cipro | Persistent organic pollutants and stable isotopes in pinnipeds from King George Island, Antarctica[END_REF]. Also, a lack of Hg data (e.g. dos [START_REF] Dos Santos | Baseline mercury and zinc concentrations in terrestrial and coastal organisms of Admiralty Bay, Antarctica[END_REF] is obvious for the organisms from this area. In this context, this study is the first, to the best of our knowledge, coupling Hg and isotopic data for organisms from King George Island.

Taking all of that into account, the study aims at assessing Hg concentrations in the light of SI ratios; at raising hypotheses on how these results are influenced by one another and also by ecological factors such as diet, distribution, exposure pathways and so on; and finally, at comparing the results with those from literature, investigating the factors that might play important roles in these eventual differences.

Sample collection took place in King George Island (62°05′S, 058°23′ W), the largest one of the South Shetland Islands, separated from the northern portion of the Antarctic Peninsula by the Bransfield Strait. The number of samples of each species, tissue and analysis is given in Table 1. All samples were kept in previously combusted (at least 420 °C for at least 4 h) containers, frozen upon arrival (-20 °C) at the Brazilian Antarctic Station and kept frozen until freeze-drying in the LabQOM(University of São Paulo, Brazil).Hg and stable isotope analyses (SIA) analyses were performed at the University of La Rochelle (France).

Briefly, bird eggs and tissues, fish and invertebrate samples were collected according to [START_REF] Cipro | Occurrence of organochlorine compounds in Euphausia superba and unhatched eggs of Pygoscelis genus penguins from Admiralty Bay (King George Island, Antarctica) and estimation of biomagnification factors[END_REF][START_REF] Cipro | Persistent organic pollutants in bird, fish and invertebrate samples from King George Island, Antarctica[END_REF]. Limpets (Antarctic limpet, Nacella concinna) were manually collected and southern elephant seal, Mirounga leonina, were collected during the austral summers of 2004/05 and 2005/06, in a fully opportunistic manner, i.e., only from animals found already dead, with no signs of degradation. Finally, vegetation samples were collected according to [START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF], which contains also the SIA data: the angiosperm Colobanthus quitensis, the mosses Brachythecium sp., Syntrichia princeps and Sanionia uncinata and the lichens Usnea aurantiaco-atra and Usnea antarctica were collected from early December 2004 to early January 2005.

Prior to stable isotope and Hg analyses, samples were lyophilized and ground to obtain a fine powder and stored in plastic vials.

Mercury analyses were carried out with an Automatic Mercury Analyser spectrophotometer, ALTEC AMA 254, which does not require an acid digestion of the samples. Aliquots ranging from 10 to 50 mg of freeze-dried sample (with no delipidation whatsoever) were directly analysed after being inserted in the oven of the apparatus. After drying, the samples were heated under an oxygen atmosphere for 3 min, and the Hg liberated and subsequently amalgamated on an Au-net. The net was then heated to liberate the collected Hg, which was measured by Atomic Absorption Spectrometry (AAS). Accuracy and reproducibility of the method were tested using dogfish liver (DOLT-2) and muscle (DORM-2) and lobster hepatopancreas (TORT-2) (National Research Council, Canada) reference standards. Standard and blanks were analysed along with each set of samples, and recoveries of the certified values and recoveries of the metals ranged from 93 to 109%. Measurements were also validated by IAEA inter-calibration exercises [START_REF] Coquery | The analytical performance study for the Medpol region: determination of trace elements and methylmercury in estuarine sediment sample[END_REF].

Concentrations are expressed in dry weight in order to compensate eventual moisture loss during freezing and to facilitate comparison between tissues and with other studies. Blanks were analysed at the beginning of each set of samples, and the detection limit of the method was 0.005 μg g -1 dry mass. Lipid content values and obtaining methods have been previously published [START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF][START_REF] Cipro | Persistent organic pollutants and stable isotopes in pinnipeds from King George Island, Antarctica[END_REF][START_REF] Cipro | Persistent organic pollutants in bird, fish and invertebrate samples from King George Island, Antarctica[END_REF].

For SIA, one aliquot of approx. 100mg of sample was placed in a test tube with 4mL of cyclohexane to remove lipids according to Chouvelon et al. (2012a) as lipids are depleted in 13 C. The mixture was shaken for an hour, then centrifuged for separation and the liquid discarded. This procedure was repeated as many times as needed until the liquid phase comes out clear. The delipidated pellet was then dried at 50 °C for 48 h. Aliquots were analysed using a Thermo Scientific Delta V Advantage, ConFlo IV interface (NoBlank and SmartEA) and Thermo Scientific Flash EA1112 Elemental Analyser. Each injection corresponded to 0.4±0.1mg for animal sample or 1.0±0.1mg for vegetation samples, encapsulated in tin cups, and there were no replicates unless C/N ratio were above than 4. In this case, delipidation would be repeated until this condition was met. Pee Dee Belemnite and atmospheric nitrogen were used as standards for calculation of δ 13 C and δ 15 N, respectively. Based on replicate measurements of internal laboratory standards, experimental precision is of ±0.15‰ and ±0.20‰ for δ 13 C and δ 15 N, respectively.

Statistics were performed using Microsoft Excel and Statsoft Statistica 12. Before analyses, data were checked for normality of distribution and homogeneity of variances using Shapiro-Wilk and Brown-Forsythe tests, respectively. Parametrical (Pearson's productmoment correlation, Tukey's HSD/ANOVA) and non-parametrical tests (Spearman's rank correlation, Kruskal-Wallis/ANOVA) followed accordingly. Statistically significant results were set at α =0.05.

Finally, normal distribution and log-link function generalized linear models (GLM) were built as follows: Hg concentrations as the dependent variable, species and tissue as categorical factors and δ 13 C and δ 15 N values as continuous predictors. Biologically relevant models were constructed incorporating the different variables and their interactions.

Continuous variables (δ 13 C and δ 15 N) were not significantly correlated in the whole dataset and could be included in the same models. Model selection was based on Akaïke's Information Criteria adjusted for small sample sizes (AICc). The model with the lowest AICc value was considered to be the most accurate. Models with AICc values differing by less than 2 have a similar level of support in the data, and the model including the least number of parameters can be regarded as the most accurate, according to the principle of parsimony [START_REF] Burnham | Model Selection and Multimodel Inference[END_REF]. Overall model support was assessed using Akaike weights (wi), following [START_REF] Johnson | Model selection in ecology and evolution[END_REF]. Residual (R² adj) analyses should be restricted to description and not be used in model selection (Burnhamand Anderson, 2002). Only models with AICc values differing by less than 10 are presented.

The results for Hg and SIA are presented in Table 1. In regard to intraspecific data for animals, the Antarctic krill, Euphausia superba, presented the lowest Hg concentration of the whole animal dataset (0.018 ±0.005 μg g -1 ), which comes in agreement with its low trophic level. Reported Hg concentrations for this species from different areas were in the same order of magnitude ranging from 0.016 to 0.041 μg g -1 along the Antarctic Peninsula [START_REF] Locarnini | Trace element concentrations in Antarctic krill, Euphausia superba[END_REF], averaged 0.033±0.013 μg g -1 (after dry weight conversion) in several locations ranging from 63 to 68°S and from 126 to 147°W and also from 107 to 157°E [START_REF] Yamamoto | Heavy Metal Accumulation in Antarctic Krill Euphausia superba[END_REF] and, with values averaging 0.035 μg g -1 in the same area as in the present study [START_REF] Dos Santos | Baseline mercury and zinc concentrations in terrestrial and coastal organisms of Admiralty Bay, Antarctica[END_REF]. These authors also bring data for Nacella concinna, averaging 0.026 μg g -1 , also consistent with our findings.

The SIA data for the invertebrates show two clearly distinct δ 13 C signatures: a depleted pelagic/oceanic for E. superba and a much more enriched for the benthic/coastal N. concinna.

These values are similar to the ones reported by a comprehensive study in the same region [START_REF] Corbisier | Trophic relationships in the nearshore zone of Martel Inlet (King George Island, Antarctica): δ 13 C stable-isotope analysis[END_REF] and are key to understand the mixing of benthic/pelagic organic matter sources throughout the local food web. This is due to the fact that, besides the feeble trophic enrichment of δ 13 C (DeNiro and [START_REF] Deniro | Influence of diet on the distribution of carbon isotopes in animals[END_REF] and also possible migration (e.g. [START_REF] Cherel | Stable isotope evidence of diverse species-specific and individual wintering strategies in seabirds[END_REF], the consumption of both sources will result in a proportionally intermediate signature (e.g. [START_REF] Carravieri | Penguins as bioindicators of mercury contamination in the Southern Ocean: birds from the Kerguelen Islands as a case study[END_REF], and therefore, clarify feeding ecology for a given species. Even though both these invertebrates are primary consumers, there is an important difference in their δ 15 N (i.e., 2.7‰, Table 1).However, the signature of the primary producers they consume presents a similar difference: 0.5‰ for phytoplankton and 4‰ for microphytobenthos in this region [START_REF] Dunton | δ 15 N and δ 13 Cmeasurements of Antarctic peninsula fauna: trophic relationships and assimilation of benthic seaweeds[END_REF], i.e., 3.5‰. Thus, the difference found is indeed a matter of baseline nitrogen signature of the organisms on which they prey and not of trophic level itself (Chouvelon et al., 2012b). According to the significant negative correlation between Hg and δ 13 C (ρ = -0.828, paired n = 6) in Nacella concinna, it seems that some δ 13 C depleted food items might represent a Hg source for this mollusc, since it feeds mainly on the relatively δ 13 C-enriched microphytobenthos [START_REF] Corbisier | Trophic relationships in the nearshore zone of Martel Inlet (King George Island, Antarctica): δ 13 C stable-isotope analysis[END_REF], but the diet of this species also comprises microbial films (grazing on the microepiflora), calcareous rhodophytes, seaweed, bryozoans and sessile spirorbid polychaetes [START_REF] Suda | The biology and ecology of the Antarctic limpet Nacella concinna[END_REF].

Moreover, a local abiotic compartment may act like an indirect source, similarly as described by Mão de Ferro et al. ( 2014): limpet samples were collected in the intertidal zone, in most cases next to glaciers, therefore subject to snow melting water, glaciers seasonal melting water (which percolates soil and bedrock before reaching the shore), marine spray, which might all function as Hg sources to the limpets themselves via respiration or to the organisms they feed on. These factors might influence δ 15 N as well, but apparently not in a significant way.

As for the Notothenia spp., muscle Hg concentrations for N. rossii from Kerguelen waters (averaging 0.255 μg g -1 dw) were three times higher than the present results [START_REF] Bustamante | Distribution of trace elements in the tissues of benthic and pelagic fish from the Kerguelen Islands[END_REF]. N. coriiceps from Laurie Island (South Orkney Islands) presented Hg concentration of 0.200 μg g -1 dw (converted from the fresh wt assuming a 75% water content) [START_REF] Moreno | Heavy metals in Antarctic organisms[END_REF]. Surprisingly, this value is closer to [START_REF] Bustamante | Distribution of trace elements in the tissues of benthic and pelagic fish from the Kerguelen Islands[END_REF] than to our findings whereas Laurie Island are situated south from the Polar Front and Kerguelen Islands, north. Indeed, predators north of the Polar Front show higher Hg concentrations than those from the southern area. For example, wandering albatrosses foraging South the Polar front are less exposed to Hg than those feeding north [START_REF] Bustamante | High feather mercury concentrations in the wandering albatross are related to sex, breeding status and trophic ecology with no demographic consequences[END_REF]Carravieri et al., 2014b). Lower values compared to our results were also reported for Notothenia spp. (with data for whole fish averaging 0.016 μg g -1 ; dos [START_REF] Dos Santos | Baseline mercury and zinc concentrations in terrestrial and coastal organisms of Admiralty Bay, Antarctica[END_REF], but the analysed individuals are assumed to be juveniles. In adult N. coriiceps from Adélie land, muscle Hg concentrations were one order of magnitude higher, i.e., 0.221 ± 0.085 μg g -1 dw [START_REF] Goutte | Trace elements in Antarctic fish species and the influence of foraging habitats and dietary habits on mercury levels[END_REF]. The δ 13 C signature found in these fish is intermediate in regard to E. superba and N. concinna which is likely due to the consumption of both oceanic/pelagic and coastal/benthic prey. Moreover, the significant positive correlation between δ 13 C and δ 15 N (ρ=0.662) in N. rossii confirms that this fish displays ontogenic changes with growth: as the fish grows, higher trophic level prey is increasingly preferred, specially small fish [START_REF] Burchett | Food, feeding and behaviour of Notothenia rossii marmorata nearshore at South Georgia[END_REF].

Data for Catharacta sp. eggs appears as the most contaminated for this tissue. This might be due to the high trophic level of this species, as shown by elevated δ 15 N. Skuas are opportunistic feeders [START_REF] Reinhardt | A review of the diets of Southern Hemisphere skuas[END_REF] and prey on other birds: eggs, chicks or the adults themselves. In locations where skuas have developed a specialisation of their diet such as in Mayes Island in the Kerguelen archipelago, chicks are fed mainly with blue petrels [START_REF] Mougeot | Predation on burrowing petrels by the brown skua (Catharacta skua lonnbergi) at Mayes Island, Kerguelen[END_REF] and they consequently present higher feather Hg levels than adults as a reflect of age class feeding habits' differences (Carravieri et al., 2014a). Indeed, the parents feed on low trophic level prey in Antarctic waters [START_REF] Cherel | Diet and feeding ecology of blue petrels Halobaena caerulea at Iles Kerguelen, Southern Indian Ocean[END_REF][START_REF] Cherel | Stable isotope evidence of diverse species-specific and individual wintering strategies in seabirds[END_REF], with a large proportion of mesopelagic fish [START_REF] Cherel | Diet and feeding ecology of blue petrels Halobaena caerulea at Iles Kerguelen, Southern Indian Ocean[END_REF][START_REF] Connan | Interannual dietary changes and demographic consequences in breeding blue petrels from Kerguelen Islands[END_REF]. In turn, breeding skuas in King George Island feed mainly on penguins and their eggs, and fish, varying somewhat between these two prey groups for the two occurring species in this location [START_REF] Peter | Ecological and morphological investigations on South Polar Skuas (Catharacta maccormicki) and Brown Skuas (Catharacta skua lonnbergi) on Fildes Peninsula, King George Island, South Shetland Islands[END_REF]. Nevertheless, considering feather data (Carravieri et al., 2014b), as well as diet [START_REF] Bocher | Influence of the diet on the bioaccumulation of heavy metals in zooplankton-eating petrels at Kerguelen archipelago, Southern Indian Ocean[END_REF], it is reasonable to assume that this diet would lead to high Hg exposure to skua from the present work. In eggs, the present values (2.520 ± 1.709 μg g -1 dw) are higher when compared to the 1.61 ± 1.22 μg g -1 dw for C. maccormicki eggs from a previous study conducted in Terra Nova Bay (Bargagli et al., 1998a). Even though in seabirds showing gender differences in their foraging strategies such as wandering albatrosses whose females had higher Hg concentrations than males despite the excretion through the eggs because they forage in more contaminated waters (north) than males (south) (Carravieri et al., 2014a), in regard to skuas there seems to be no gender segregation in both Hg concentration and SIA, not even during the breeding season (Bearhop et al., 2000a).

SIA data puts skuas in roughly the same ecological niche as previous literature (Bearhop et al., 2000a(Bearhop et al., , 2000b)), particularly when liver is compared to blood, since both these tissues have short and similar turnover rates [START_REF] Hobson | Assessing avian diets using stable isotopes I: turnover of 13 C in tissues[END_REF]. Unpublished data (Cipro et al., personal data) confirms these previous statements, also presenting sub-Antarctic species (from Kerguelen) with higher Hg concentrations than Antarctic ones (from Adélie Land) occupying similar ecological niches, the Cape petrel Daption capense presented remarkably high Hg concentrations despite its δ 15 N closer to penguins which showed much lower Hg concentrations. Such high concentrations might be due to two factors. The first one is that its diet is, in a general way, almost entirely composed by oceanic/pelagic prey [START_REF] Coria | Diet of Cape petrel Daption capense during the posthatching period at Laurie Island, South Orkney Islands, Antarctica[END_REF] whereas pygoscelids penguins might, in a lesser way, consume some coastal/benthic prey (especially Gentoo penguins, P. papua). Pygoscelids also show a reproductive temporal shift, mostly due to their respective distributions, and not necessarily to food availability. These penguins breed in different periods of the austral summer in this area, with Adélies early breeding, followed by Gentoos and then Chinstraps [START_REF] Trivelpiece | Ecological segregation of Adelie, gentoo, and chinstrap penguins at King George Island, Antarctica[END_REF]. The second one is that most of the Cape petrel diet in the region [START_REF] Casaux | Lanternfish (Myctophidae) in the diet of the Cape petrel Daption capense at the South Shetland and South Orkney Islands, Antarctica[END_REF][START_REF] Coria | Diet of Cape petrel Daption capense during the posthatching period at Laurie Island, South Orkney Islands, Antarctica[END_REF] is composed by Myctophidae fish, which have higher Hg concentrations compared to epipelagic ones (e.g., Chouvelon et al., 2012a;[START_REF] Monteiro | Mercury concentrations in prey fish indicate enhanced bioaccumulation in mesopelagic environments[END_REF] and species that prey on these resources are therefore more prone to accumulate Hg [START_REF] Blévin | Wide range of mercury contamination in chicks of southern ocean seabirds[END_REF]Carravieri et al., 2014b;[START_REF] Cipro | Trace elements in tissues of white-chinned petrels (Procellaria aequinoctialis) from Kerguelen waters, Southern Indian Ocean[END_REF].

The giant Petrel Macronectes giganteus showed much higher Hg concentrations in the present study than the reports from the literature. In giant petrels from Bird Island, South Georgia, Hg concentrations in the liver were averaging 2.16 (males) and 4.95 (females) μg g -1 dw [START_REF] González-Solís | Metals and seleniumas bioindicators of geographic and trophic segregation in giant petrels[END_REF] i.e., from 5 to 10 times lower than our findings. Birds sampled in the South Shetland Islands and the Antarctic Peninsula had also Hg concentrations one order of magnitude lower, averaging 1.37 μg g -1 dw [START_REF] Szefer | Mercury and major essential elements in seals, penguins, and other representative fauna of the Antarctic[END_REF]. This marked difference could be due to our small sampling number or to some local effect, since our samples present similar SIA levels to those from Bird Island, South Georgia, the closest site that could be found in the literature for comparison [START_REF] Forero | Stable isotopes reveal trophic segregation by sex and age in the southern giant petrel in two different food webs[END_REF].

Hg concentrations in the eggs of Gentoo penguins P. papua and Adélie penguins P. adeliae were similar to that reported in the literature [START_REF] Moreno | Heavy metals in Antarctic organisms[END_REF], after dry weight conversion assuming 70% water content), i.e., ~0.13 μg g -1 dw. This same work also presents data for P. adeliae liver, ~1.5 μg g -1 dw, also in the same range with our findings. As for liver, literature shows more contaminated liver samples for two pygoscelids in [START_REF] Szefer | Mercury and major essential elements in seals, penguins, and other representative fauna of the Antarctic[END_REF], averaging 2.01 (P. adeliae) and 34.7 μg g -1 dw (P. papua). However, [START_REF] Smichowski | Monitoring trace elements in selected organs of Antarctic penguin (Pygoscelis adeliae) by plasmabased techniques[END_REF] presents data six times lower for P. adeliae chicks' liver collected also in King George Island, which provides evidence of Hg bioaccumulation in this species. Significant correlations for birds are broken down in Fig. 1. Only paired analyses are shown, therefore the number of samples may vary in regard to the total.

Since pygoscelids fast during the egg formation period [START_REF] Astheimer | The timing and energetic consequences of egg formation in the Adélie penguin[END_REF][START_REF] Polito | Stable isotopes reveal regional heterogeneity in the pre-breeding distribution and diets of sympatrically breeding Pygoscelis spp. penguins[END_REF], the influence on egg results of the trophic level of prey during breeding season is not to be considered. The significant negative correlation between Hg and δ 13 C (ρ = -0.627) in P. antarctica eggs could be related to the fact that fasting breeding adult penguins present a decrease in δ 13 C plasma signature during fasting [START_REF] Cherel | Nutrition, physiology, and stable isotopes: new information from fasting and molting penguins[END_REF], whereas Hg levels might increase simply due to the opposite of the growth dilution effect on egg laying females and consequent transfer (e.g. [START_REF] Agusa | Body distribution of trace elements in black-tailed gulls from Rishiri Island, Japan: age-dependent accumulation and transfer to feather and eggs[END_REF]. A significant negative correlation was found between Hg and δ 15 N (ρ= -0.900) in P. papua egg could be related to different proportions of fish and squid in the diet prior to breeding [START_REF] Polito | Stable isotopes reveal regional heterogeneity in the pre-breeding distribution and diets of sympatrically breeding Pygoscelis spp. penguins[END_REF]. Since the Gentoo penguin presents considerable plasticity in its diet when compared to the other two less flexible pygoscelids [START_REF] Miller | Flexible foraging strategies of Gentoo penguins Pygoscelis papua over 5 years in the South Shetland Islands, Antarctica[END_REF], consuming prey from coastal/oceanic and benthic/pelagic origins could result in different δ 15 N baselines and therefore misleading conclusions (see Chouvelon et al., 2012b). This is a similar issue as previously discussed for krill and limpets.

The significant positive correlation between δ 13 C and δ 15 N (ρ = 0.562) in P. antarctica liver is due to the concomitant enrichment of these elements even though it occurs in a much lesser degree for 13 C (DeNiro and Epstein, 1978). Finally, a significant positive correlation was found between Hg and δ 15 N (ρ = 0.615) in P. papua liver likely due to the consumption of higher trophic level prey for adults all year round. Hg concentrations for the eggs of the kelp gull Larus dominicanus fall within the same order of magnitude than the value (i.e., 0.46 μg g -1 ) reported by [START_REF] Moreno | Heavy metals in Antarctic organisms[END_REF] but data for deeper comparison are dramatically lacking for this species. However, Hg concentrations in the liver was the second highest value for seabirds, which is likely due to the influence of high consumption rates of N. concinna during the austral summer [START_REF] Favero | Trophic relationships between the kelp gull and the Antarctic limpet at King George Island (South Shetland Islands, Antarctica) during the breeding season[END_REF]. The influence of this benthic prey can also be seen in the highly enriched δ 13 C signature in kelp gull liver and moreover, in the resulting PCB congener profiles in L. dominicanus, already highlighted in a previous study [START_REF] Cipro | Persistent organic pollutants in bird, fish and invertebrate samples from King George Island, Antarctica[END_REF].

The Antarctic tern Sterna vittata is, in this study, a sui generis case, firstly because only one sample was available and secondly because its SIA results would suggest a coastal/benthic diet whereas the literature [START_REF] Casaux | Aspects of the foraging behaviour of the Antarctic Tern Sterna vittata gaini at Harmony Point, South Shetland Islands[END_REF] reports it as a pelagic predator.

More samples are needed to explore this contradictory observation.

Pinnipeds will be further discussed as a group due to the very small and disperse sampling.

For animals interspecific data, two statistically (Tukey HSD tests) different intrahomogenous groups for Hg levels: one group with the giant petrel Macronectes giganteus liver and one group with all the other samples. This relative homogeneity around the Antarctic Peninsula had previously been reported by [START_REF] Brasso | Penguin eggshell membranes reflect homogeneity of mercury in the marine food web surrounding the Antarctic Peninsula[END_REF].

Particularly for pinnipeds, the data for A. gazella and M. leonina liver are clearly distinguished from the others. A previous analysis of the SIA data [START_REF] Cipro | Persistent organic pollutants and stable isotopes in pinnipeds from King George Island, Antarctica[END_REF] for A. gazella had already shown a clear shift in foraging area (closer to the shore) and towards higher trophic level, probably due to the increase of fish consumption. All of this corroborates this much increased Hg value, adding up to the fact that they stand on a higher trophic level and are more prone to biomagnification. Interpretation for M. leonina is yet more complicated since sex-specific foraging strategies are reported [START_REF] Cherel | Stable isotopes reveal the trophic position and mesopelagic fish diet of female southern elephant seals breeding on the Kerguelen Islands[END_REF][START_REF] Lewis | Sex-specific foraging strategies and resource partitioning in the southern elephant seal (Mirounga leonina)[END_REF]: in a general way females' diets vary less than the ones of males. Females forage mostly on low trophic level prey from deep waters, such as cephalopods or myctophid fish (according to the location), whereas, males appear to be specialist foragers, divided among several different strategies, which might also include pelagic invertebrates, benthic fish and adult pelagic fish.

Such a difference can be seen as a strategy to avoid intra and interspecies competition, especially during the breeding season [START_REF] Cherel | Stable isotopes reveal the trophic position and mesopelagic fish diet of female southern elephant seals breeding on the Kerguelen Islands[END_REF][START_REF] Lewis | Sex-specific foraging strategies and resource partitioning in the southern elephant seal (Mirounga leonina)[END_REF]. Comparison with [START_REF] Szefer | Mercury and major essential elements in seals, penguins, and other representative fauna of the Antarctic[END_REF] shows that Lobodon carcinophagus muscle Hg concentrations are one order of magnitude lower and the ones of Leptonychotes weddellii kidney were two orders of magnitude lower. Since no sex, age or biometrics data could be taken, the interpretation cannot go any further.

Concerning intraspecific data for vegetation, Hg values for several moss species from Edmonson Point, Antarctica ranged from 0.05 to 0.15 μg g -1 dw (Bargagli et al., 1998b), which comes in reasonable agreement with the data from the present study. Another work, in the same area as the present one [START_REF] Dos Santos | Baseline mercury and zinc concentrations in terrestrial and coastal organisms of Admiralty Bay, Antarctica[END_REF] presents values for mosses of 0.023 (Bryum spp.) and 0.040 (Polytrichum spp.) μg g -1 , which are roughly 1.5 to 5 times lower than the ones hereby reported. Also, Hg concentrations of 0.036 μg g -1 for lichens (Usnea spp.) are 4 times lower than the ones from the present work, in average. Significant positive correlation between Hg and δ 13 C (ρ = 0.828) in Sanionia uncinata could be related to a lesser marine influence, as further discussed.

For vegetation, Tukey HSD tests separated the less Hg contaminated mosses Syntrichia princeps and Sanionia uncinata from the more Hg contaminated lichens Usnea spp. in two statistically different intra-homogenous groups. This finding comes in agreement with physiologic differences, mainly water dependence. An analogous separation between lichens and mosses is noticed in [START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF] for POPs. And since moister habitats and marine influence are related to lower δ 13 C and higher δ 15 N [START_REF] Lee | Carbon and nitrogen isotope composition of vegetation on King George Island, maritime Antarctic[END_REF], it would be reasonable to assume that dry deposition phenomena are more significant as Hg source in this case.

However, the clear influence of secondary sources in recent literature [START_REF] Choy | An isotopic investigation of mercury accumulation in terrestrial food webs adjacent to an Arctic seabird colony[END_REF][START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF][START_REF] Roosens | Penguin colonies as secondary sources of contamination with persistent organic pollutants[END_REF] might overcome that conclusion. Therefore our results differ largely from the ones by dos [START_REF] Dos Santos | Baseline mercury and zinc concentrations in terrestrial and coastal organisms of Admiralty Bay, Antarctica[END_REF] not only in an intraspecific way, but also, and mainly, in the clear difference between lichens and mosses Hg levels.

In order to understand the Influence of species, tissue, δ 13 C and δ 15 N on how these variables affect the resulting Hg concentrations in this large array of samples, normal distribution and log-link function generalized linear models (GLM) were built as previously described. The whole dataset was split in animals and vegetation subsets and model results are shown in Table 2. Taking the results into account, the best model to explain Hg levels in the fauna is composed by both δ 15 N and tissue, followed by δ 15 N, species and tissue. Examples of the importance of such variables in the accumulation of trace elements are found throughout the scientific literature: δ 15 N (e.g. [START_REF] Anderson | Element patterns in albatrosses and petrels: influence of trophic position, foraging range, and prey type[END_REF], tissue (e.g. [START_REF] Bustamante | Distribution of trace elements in the tissues of benthic and pelagic fish from the Kerguelen Islands[END_REF] and species (e.g. [START_REF] Borgå | Regional and species specific bioaccumulation of major and trace elements in Arctic seabirds[END_REF] as well as for both δ 15 N and tissue (e.g. [START_REF] Campbell | Mercury and other trace elements in a pelagic Arctic marine food web (Northwater Polynya, Baffin Bay)[END_REF], as in our findings.

The models for vegetation, in turn, are composed by δ 13 C and species followed by δ 13 C, δ 15 N and species. Firstly, it is, in a broad sense, the repetition of Tukey HSD tests previously performed, since mosses and lichens were placed in two statistically different groups and moreover, a clear separation in δ 13 C for lichens and mosses, with no overlapping whatsoever, can be noticed in [START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF]. The influence of δ 15 N in Hg concentrations, which in vegetation cannot be linked to trophic level is, firstly, linked to δ 13 C [START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF][START_REF] Park | Polychlorinated biphenyl congeners in soils and lichens from King George Island, South Shetland Islands, Antarctica[END_REF] and also thought to be linked to the exposure to secondary contaminant sources, particularly the degradation of animal derived organic matter [START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF].

Figure 1. Significant correlations between Hg (μg g -1 dwt) and stable isotopes (‰) for A) Pygocelis antarctica eggs; B) P. antarctica liver and C) Pygocelis papua egg and liver. Empty markers for P. antarctica, filled ones for P. papua. Round markers for eggs, triangular ones for liver.

(

  2004-05 summer) in the intertidal zone. Fish (marbled rockcod, Notothenia rossii and black rockcod N. coriiceps) were collected (2006-07 summer) by mid-water nets or by line and hook. Only unhatched bird eggs (skua, Catharacta sp., identified only until the genus level; kelp gull, Larus dominicanus; and Antarctic tern, Sterna vittata) were collected (2004-05 and 2005-06 summers) so as not to interfere with breeding success. Bird livers (cape petrel, Daption capense and giant petrel, Macronectes giganteus in addition to the previous species) were collected (from 2004-05 up to 2007-08 summers) only from already dead animals with no evident signs of disease, decomposition or emaciation. No further attempt to determine the cause of death was performed. Pinnipeds samples were collected according to Cipro et al. (2012), which already contains the SIA data: samples of the Weddel seal, Leptonychotes weddellii; Antarctic fur seal, Arctocephalus gazella; crabeater seal, Lobodon carcinophagus;

Acknowledgements

The authors acknowledge the financing organisations in the University of La Rochelle: CPER (Contrat de Plan État-Région) and FASEP (Fonds d'étude et d'aide au secteur privé). C.V.Z. Cipro received scholarships from CAPES (Coordination for the Improvement of Higher Education Personnel via the Science without Borders programme, Brazil, grant number 11490/13-6) during the course of this work. The present work is a result of the projects "Environmental Management of Admiralty Bay, King George Island, Antarctica: persistent organic pollutants and sewage" (CNPq No. 55.0348/2002-6), and "Modeling the bioaccumulation of organic pollutants throughout an Antarctic food web" (CNPq 550018/2007-7) and also contributes to the Brazilian National Science and Technology Institute on Antarctic Environmental Research (INCTAPA). Financial support was obtained from the Brazilian Antarctic Program (PROANTAR), Ministry of the Environment (MMA) and Conselho Nacional de Desenvolvimento Cientifico e Tecnológico (CNPq). Antarctic logistics was provided by Secretaria da Comissão Interministerial para os Recursos do Mar. (SECIRM). C. Churlaud and M. Brault-Favrou from the "Plateforme Analyses Elémentaires" of LIENSs are thanked for their assistance during Hg analyses. G. Guillou and B. Lebreton from the "Plateforme Analyses Isotopiques" of LIENSs are thanked for their assistance during SIA analyses. Finally, authors wish many thanks to everyone involved in the fieldwork.


Sylvain Lobry
email: sylvain.lobry@telecom-paristech.fr.

Loïc Denis
email: loic.denis@univ-st-etienne.fr.

Florence Tupin
email: florence.tupin@telecom-paristech.fr

Multi-temporal SAR image decomposition into strong scatterers, background, and speckle

Keywords: SAR, Image decomposition, TV, L0, Change detection

published or not. The documents may come

I. INTRODUCTION

Due to the coherent nature of SAR imaging technique, SAR images suffer from strong fluctuations related to the speckle phenomenon. While containing information about the subresolution texture of the scene, speckle is often regarded as undesirable noise for image interpretation tasks. Speckle is then typically modeled as a multiplicative noise.

While image denoising and restoration are among the oldest challenges addressed in image processing, most algorithms developed in the field of image processing are designed to deal with an additive Gaussian noise and are thus not directly applicable to SAR images.

The simplest way to reduce speckle fluctuations is spatial multi-looking which amounts to averaging pixel values within a sliding window. This speckle variance reduction is obtained at the cost of a resolution loss proportional to the size of the averaging window. While it may produce satisfying results in homogeneous areas, it strongly blurs textured areas, edges between regions and bright scatterers.

Numerous approaches have been proposed to prevent the introduction of blur by mixing values from distinct regions. Lee et al. [START_REF] Lee | Speckle filtering and coherence estimation of polarimetric SAR interferometry data for forest applications[END_REF] locally select the best window among a few oriented windows. The IDAN algorithm [START_REF] Vasile | Intensity-Driven Adaptive-Neighborhood Technique for Polarimetric and Interferometric SAR Parameters Estimation[END_REF] builds an adaptive window by region growing. Several methods have been derived from the non-local means (NL-means) approach by Buades et al. [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]. These methods select similar pixels in an extended window based on patch-similarity [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF]- [START_REF] Deledalle | Nlsar: A unified nonlocal framework for resolution-preserving (pol)(in)sar denoising[END_REF]. Another family of methods reduce speckle by regularization, i.e., by computing the maximum a posteriori estimate under a given prior. Wavelet-based approaches model the distribution of wavelet coefficients [START_REF] Achim | SAR image denoising via Bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF]- [START_REF] Xie | SAR speckle reduction using wavelet denoising and Markov random field modeling[END_REF]. Total variation (TV) regularization penalizes variations between neighboring pixels while preserving sharp edges [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. Total variation has been applied to the regularization of SAR amplitudes [START_REF] Aujol | Image decomposition application to sar images[END_REF]- [START_REF] Palsson | SAR image denoising using total variation based regularization with SUREbased optimization of the regularization parameter[END_REF], SAR intensities [START_REF] Aubert | A variational approach to removing multiplicative noise[END_REF], [START_REF] Steidl | Removing multiplicative noise by douglasrachford splitting methods[END_REF] and log-transformed intensities [START_REF] Bioucas-Dias | Multiplicative noise removal using variable splitting and constrained optimization[END_REF] using different optimization strategies (discrete optimization by graph-cuts, gradient descent, Douglas-Rachford splitting or the alternating directions method of multipliers). We refer the reader to the two recent review papers [START_REF] Argenti | A tutorial on speckle reduction in synthetic aperture radar images[END_REF] and [START_REF] Deledalle | Exploiting Patch Similarity for SAR Image Processing: The nonlocal paradigm[END_REF] for an in-depth analysis of speckle reduction methods.

Man-made structures such as buildings, fences or transmission towers produce very strong back-scattering, with intensities much larger than the surrounding area. Such scatterers are especially numerous in urban areas. Isolated strong scatterers are challenging for speckle reduction methods because they generally break the statistical assumptions made about the radar scene: repetition of similar patches within the search window (patch-based methods), sparse representation in the wavelets domain (wavelets-based methods), or piece-wise constant regions (TV minimization). It is then necessary to identify and process these points separately to prevent from spreading these large values. Strong scatterers can be detected using likelihood ratio tests [START_REF] Lopes | Adaptive processing of multilook complex SAR images[END_REF], [START_REF] Lopes | Structure detection and statistical adaptive speckle filtering in SAR images[END_REF]. These detectors compare the values in the center of a window with the rest of the window, considered as purely background. In dense urban areas, the presence of other point-like scatterers in the vicinity strongly deteriorates the performance. It is thus necessary to perform jointly the detection of strong scatterers with the estimation of the background radiometry. The decomposition method that we propose in the present paper achieves such a joint detection of point-like scatterers and estimation of the background. It applies to multi-temporal stacks of images by detecting scatterers at each date, even if they are not present at other dates due to some change occurring.

This paper extends the decomposition models presented in [START_REF] Denis | Exact discrete minimization for TV+L0 image decomposition models[END_REF], [START_REF] Lobry | Sparse & smooth decomposition models for multi-temporal SAR images[END_REF] in the following ways: We derive the decomposition model from a connection between scatterers detection and maximum a posteriori estimation with a regularization term based on L0 pseudo-norm (section II). The proposed decomposition model can provide a multi-temporal image of detected scatterers and a smoothly varying multi-temporal background image, or a simpler decomposition into multitemporal scatterers image and a time-invariant background. Computing the decomposition requires solving a challenging large-scale combinatorial problem. We describe two algorithmic approaches (section III): an exact discrete optimization based on finding the minimum cut in a large 4-D graph, and a faster but approximate method based on block processing with overlaps. We show on numerical experiments that the choice of the L0 non-convex penalty significantly improves the detection of scatterers compared to the usual L1 norm surrogate (section IV). We illustrate how the detected scatterers can be used in a subsequent processing step for change detection purposes (section V).

II. SPARSE + SMOOTH DECOMPOSITION MODEL FOR MULTI-TEMPORAL DATA

A. From detection to estimation under L0 prior

Let us consider first the problem of detecting strong scatterers in a speckled region of given radiometry u B . At each pixel i, the detection problem amounts to deciding between two hypotheses: the absence (H 0 ) or presence (H 1 ) of a strong scatterer in addition to the other scatterers that form the background:

H 0 : u S i = 0 H 1 : u S i > 0 , (1)
where u S i is the radiometry of the additional scatterer. Under the assumption of uncorrelated speckle and point-like scatterers, this hypothesis test can be performed independently at each pixel by studying the likelihood radio:

log p(v i |u B i + u S i ) p(v i |u B i ) H1 ≷ H0 λ , (2)
with v i the observed amplitude at pixel i (i.e., specklecorrupted) and λ a threshold. This likelihood ratio test can not be readily applied since the radiometry u S i is unknown. Replacing u S i by the maximum likelihood estimate leads to the generalized likelihood ratio test (GLRT) [START_REF] Kay | Fundamentals of statistical signal processing[END_REF]:

arg max u S ≥0 log p(v i |u B i + u S i ) H1 ≷ H0 λ + log(p(v i |u B i )) . (3)
In words, the presence of an additional scatterer is detected at pixel i if the log-likelihood of the observed amplitude v i is improved at least by λ when hypothesizing that a scatterer is present in addition to the background scatterers. The threshold λ sets the false alarm rate: increasing λ reduces the probability of wrongly selecting hypothesis H 1 under H 0 , at the cost of reducing the detection probability. This detection problem can be stated as an estimation problem, i.e., as the estimation of the radiometry of the additional scatterer u S i , where a zero estimated radiometry corresponds to the absence of the additional scatterer at pixel i:

u S i = arg min u S i ≥0 -log p(v i |u B i + u S i ) + λ • (u S i = 0) , (4)
where the notation λ • (u S = 0) refers to value λ whenever u S = 0, and 0 when u S = 0. Equation ( 4) is completely equivalent to equation (3), and can be extended to a whole image. The estimation of the vector of all radiometries u S , for a given vector of corresponding background radiometries u B and under our pixel-independence assumption, is given by:

u S = arg min u S ≥0 - i log p(v i |u B i + u S i ) + λ u S 0 , (5)
where the positivity constraint u S ≥ 0 applies componentwise and the L0 pseudo-norm corresponds to the number of non-zero components of vector u S .

B. Joint estimation of the background and of strong scatterers

Direct applicability of the estimation formulation given in equation ( 5) is limited because it requires the knowledge of the background radiometry u B . Estimation of this background radiometry requires to exclude strong scatterers, i.e., the detection of all scatterers in u S . The two components u B and u S must thus be jointly estimated:

( u S , u B ) = arg min u S ≥0,u B ≥0 - i log p(v i |u B i + u S i ) + λ u S 0 + ψ(u B ) , ( 6
)
where ψ is a smoothness term necessary to enforce the regularity of the background and thus prevent problem degeneracy 1 . Regularization term ψ introduces a coupling between all unknowns, i.e., the minimization problem in ( 6) is not separable. With the L0 term, the optimization problem is combinatorial (non-convex and discontinuous). Given the large scale of the problem (typically millions to billions of unknown pixel radiometries to estimate), two approaches can be considered:

• convex relaxation: the L0 term is generally replaced by the L1 norm to turn the combinatorial minimization problem into a convex minimization problem; • combinatorial optimization: direct minimization of [START_REF] Chen | Nonlocal filtering for polarimetric SAR data: A pretest approach[END_REF] is possible for some specific choices of the regularization term ψ. The first approach suffers from two drawbacks: (i) under speckle noise, the neg-log-likelihood is non-convex, which prevents direct application of standard convex optimization methods; (ii) with the L1 norm, the estimation problem is no longer equivalent to the detection problem. We show in section IV that this leads to a non-constant false alarm rate: scatterers are mainly detected in darker areas.

The second approach consists of applying a combinatorial optimization method to the problem, for a given discretization of the radiometries. It restates the original optimization problem as a minimum cut / maximum flow search on a graph [START_REF] Ishikawa | Exact optimization for Markov random fields with convex priors[END_REF]. It applies to discrete optimization problems involving (possibly non-convex) separable terms and convex pairwise terms, and is well suited (i.e., the graph is sparsest) for anisotropic total variation:

T V (u) = i∼j w i,j |u i -u j | , (7)
where i ∼ j indicates that i and j are indices of two neighboring pixels, and w i,j are weights. In the simplest form, only the 4 nearest neighbors are considered and all weights w i,j are equal to 1. To reduce the metrication effects, farther neighbors can also be considered, with smaller weights, see for example [START_REF] Boykov | Computing geodesics and minimal surfaces via graph cuts[END_REF]. In the following, we will consider a smoothness term ψ based on total variation (7) so that a graph-cut method can be applied to solve the non-convex minimization problem.

C. Multi-temporal decomposition model

We consider the decomposition of a time series of SAR images into 3 components: a spatio-temporal background, additional point-like scatterers, and speckle noise. A multitemporal stack, represented by a vector v, is decomposed by jointly estimating the background u B and the additional scatterers u S . The speckle component n then corresponds to the ratio between observed amplitudes and the estimated radiometries, according to the multiplicative speckle model:

v = (u B + u S ) × n , (8)
where the multiplication is applied component-wise. Under Goodman's fully developed speckle model [START_REF] Goodman | Speckle phenomena in optics: theory and applications[END_REF], the observed amplitude v t,i at pixel i and time t follows a Rayleigh distribution:

p(v t,i |u t,i ) = 2 v t,i u 2 t,i exp -v 2 t,i u 2 t,i , (9)
where u t,i = u B t,i + u S t,i is the radiometry at pixel i and time t. Other statistical distributions may be considered, for example Rice distribution which is more accurate when one strong scatterer is dominant in the resolution cell (u S t,i u B t,i ). Note that with the choice of Rayleigh distribution, the generalized likelihood ratio test (3) takes the form:

v 2 i u B 2 i -log v 2 i u B 2 i H1 ≷ H0 λ + 1 . ( 10
)
Since the test depends only on the ratio v i /u B which is Rayleigh distributed with parameter 1 under H 0 , the threshold λ sets a constant false alarm rate (i.e., λ need not be tuned with the level u B ).

To apply the joint estimation of the background and detection of scatterers described in the preceding sections, the regularization term ψ(u B ) needs to be set. Apart from the presence of many strong-scatterers, SAR images can reasonably well be described by smooth regions (with few variations of the radiometry) separated by sharp edges. In a time series, changes with time are rather smooth (e.g., vegetation changes) and pixels within a region of homogeneous radiometry generally keep comparable radiometries. When an abrupt change occurs, it is important to prevent from smoothing in time its apparition. We consider using a spatio-temporal total variation to impose some spatial and temporal smoothing while preserving sharp edges / changes. The proposed spatio-temporal regularization is:

TV α 3D (u B ) = t i∼j w i,j |u B t,i -u B t,j | + α i t |u B t+1,i -u B t,i | , ( 11
)
where α is the weight along the temporal direction. If α is set to a very large value, no temporal change is possible and a constant (in time) background component is estimated. We denote TV ∞ 3D the regularization in this simplified case of a time-invariant background.

With Rayleigh likelihood and the TV 3D regularization, the joint decomposition amounts to solving the following minimization problem:

arg min u S ≥0,u B ≥0 i,t 2 log(u B t,i + u S t,i ) + v 2 t,i (u B t,i + u S t,i ) 2 + λ u S 0 + β TV α 3D (u B ) . ( 12
)

III. OPTIMIZATION USING GRAPH-CUTS

A. Exact discrete optimization

To apply the graph-cuts combinatorial optimization method to the minimization problem defined in equation [START_REF] Argenti | Multiresolution MAP Despeckling of SAR Images Based on Locally Adaptive Generalized Gaussian pdf Modeling[END_REF], we need that it be formulated as:

arg min x t,i f 0 (x i,t ) + (i,t)∼(i ,t ) f 1 (x i,t , x i ,t ) , (13)
i.e., as the sum of separable terms f 0 and pairwise terms f 1 involving (spatio-temporal) neighbors (i, t) and (i , t ).

For any fixed background u B , the equivalence between the estimation formulation (4) and the detection formulation [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] shows that the computation of the optimal vector of scatterers u S is readily obtained by simple testing against a threshold, for each pixel. It follows from [START_REF] Deledalle | Nlsar: A unified nonlocal framework for resolution-preserving (pol)(in)sar denoising[END_REF] that the optimal u S is given by: By substituting u S with its optimal value u S (u B ) in the minimization problem [START_REF] Argenti | Multiresolution MAP Despeckling of SAR Images Based on Locally Adaptive Generalized Gaussian pdf Modeling[END_REF], we get the sought formulation of equation ( 13), with

u S t,i (u B t,i ) =      v t,i -u B t,i if v t,i > u B t,i and v 2 i u 2 B -log v 2 i u 2 B ≥ λ + 1 , 0 otherwise . (14)
x = u B (15) f 0 (x i,t ) = 2 log(u B t,i + u S t,i (u B t,i )) + v 2 t,i (u B t,i + u S t,i (u B t,i )) 2 + λ • (u S t,i (u B t,i ) = 0) (16) f 1 (x i,t , x j,t ) = β • w i,j • |u B t,i -u B t,j | (17) f 1 (x i,t , x i,t+1 ) = β • α • |u B t+1,i -u B t,i | . ( 18
)
We apply the graph construction method of Ishikawa [START_REF] Ishikawa | Exact optimization for Markov random fields with convex priors[END_REF] to restate the minimization problem as the search of a minimum cut in a graph. We choose a (not necessarily uniform) quantization of values of the background: {q 1 , q 2 , . . . , q max }. We then create a set of nodes representing all quantization levels for each (spatio-temporal) pixel, see Fig. 1(a). Oriented edges are created from the source to the nodes associated with level q max , then within a set of nodes associated to a given pixel, each node is connected to its preceding quantization level. Finally, the nodes representing level q 1 are connected to the sink. The capacity of the edge pointing to the node associated with pixel (i, t) and quantization level q c is set to f 0 (q c ). Upstream edges with infinite capacity are also created to enforce that a minimum cut separates the source and the sink at a unique level for each pixel.

Spatial neighbors are then connected with bi-directional edges: two nodes representing the same quantization level q c are connected and the capacity of the edge is set to β • w i,j • (q c+1 -q c ). Finally, temporal neighbors are connected by creating bi-directional edges between nodes representing the same quantization level q c , with a capacity set to α • β.

In figure 2, we show the first and the last images obtained by applying the proposed decomposition on a SAR time series from Saint-Gervais (France) acquired by TerraSAR-X satellite.

Computational and memory complexity: The graph construction requires a large amount of memory that limits its application to small regions of interest and/or short time series. The number of vertices is proportional to the number of images in the time series, the number of pixels of each image, and the number of quantization levels. The number of edges is proportional to the number of nodes (about 8 times the number of nodes). We used the graph-cuts implementation described in [START_REF] Boykov | An experimental comparison of mincut/max-flow algorithms for energy minimization in vision[END_REF]. In this graph implementation, each vertex requires 48 bytes of storage, and each edge requires 32 bytes. Hence, the memory footprint of the graph representation limits the size of the series of images that can be processed. For example, a series of 20 images of size 300×400 with a quantification into 50 levels requires 33.7GiB of RAM for the graph construction. However, when ψ = TV ∞ 3D , the problem only involves one background. Solving the same problem with this regularization only requires 1.35GiB, simplifying the exact optimization of the proposed model.

Regarding the computational cost, the worst-case complexity of the minimum cut is O(EV 2 |C|) where C is the value of the minimum cut, E is the number of edges and V the number of vertices. In practice, the experimental complexity scales almost linearly with the number of nodes [START_REF] Boykov | An experimental comparison of mincut/max-flow algorithms for energy minimization in vision[END_REF]. On a computer with an Intel R Xeon(R) CPU E5-1620 with 16Gb of RAM, the algorithm takes 52,04s to compute a decomposition on 2 images of 300 × 400 pixels with 50 levels of quantification. Note that we do not fully benefit from the power of the processor as the implementation of the algorithm used is single-core. Algorithms computing the grid-cut in parallel have also been proposed (see [START_REF] Liu | Parallel graph-cuts by adaptive bottom-up merging[END_REF]) and more compact memory representations that exploit the regularity of the graph have been introduced in [START_REF] Jamriska | Cache-efficient graph cuts on structured grids[END_REF].

To reduce the computational and memory complexity, a subset of all quantization levels can be considered at a time to get an approximate solution, see [START_REF] Denis | SAR image regularization with fast approximate discrete minimization[END_REF], [START_REF] Shabou | A markovian approach for insar phase reconstruction with mixed discrete and continuous optimization[END_REF]. We describe in section III-B another approach based on block processing to further reduce memory requirements.

Quantization of the radiometry: The selection of the quantization levels used to represent the background radiometries u B requires to find a trade-off between accuracy (the more quantization levels the better) and the memory complexity (the graph footprint is proportional to the number of quantization levels). We used an adaptive method to set the quantization levels. We manually select the proportion p of the image considered to be in the background (a typical value is around 95%, but may be adapted according to the density of scatterers). The quantization levels are then computed from the quantiles of the p lowest observed values of the first image. Note that only the background values u B need to be quantified, the bright scatterers component u S with the highest dynamic range is not quantified (see equation 14).

B. Memory efficient optimization

To apply the method to large images and/or long time series, it is necessary to develop a method for limited-memory graphcuts optimization.

a) Method: Memory usage of the graph-cut method is proportional to the number of pixels in the series. The required memory can thus be reduced by computing the optimization locally, on spatio-temporal blocks extracted from the time series. Even though our model involves only second order cliques, limiting the direct interactions to the immediate neighborhood, the maximum a posteriori estimate involves longrange correlations, i.e., during the optimization, regularization effects are propagated over long distances. Simple division of the image into smaller blocks thus results into visible block artifacts, as can be observed in figure 3a. These block artifacts are due to the lack of context: a constant area in the global optimum that gets divided into two regions during the blockprocessing is represented by two different (constant) values with an artificial discontinuity between the blocks.

To introduce context in the window F of interest, it is necessary to perform the optimization on a bigger window C containing all the objects partially presents in F . The proposed method works as follows:

1) To process a given window F (named filling window in the following), extract a larger computation window C such that F ⊂ C. 2) Perform the optimization using the graph-cuts method described in section III-A on the spatio-temporal window C. 3) Keep the results of the decomposition only in the (smaller) filling window F . 4) Slide the filling window F and repeat from step 1 until all the image has been covered by the filling window.

b) Results: We show in figure 4 the root mean squared error obtained by the block-processing approach for different sizes of the computation window (the size of the block F is kept constant and equal to 50 × 50 pixels for this single-date image). The required memory grows quadratically with the spatial window size. It can be observed that on this image, the exact solution is obtained for computation windows C larger or equal to 150 × 150 pixels (which requires only 8.6% of the amount of memory needed to process the whole image at once). When the computation window F is strictly smaller than the image size, the solution computed on the block F is not guaranteed to match the solution obtained from the whole image, i.e., it is only approximate.

IV. ADVANTAGE OF THE L0 PSEUDO-NORM OVER THE L1 RELAXATION

The combinatorial optimization method described in section III is applicable to the L0 pseudo-norm and to the L1 norm, and provides in both cases the global optimum. It is thus possible to compare the performance of the two formulations: the L0 penalty derived from the GLRT and the L1 norm widely used to obtain a convex relaxation of the minimization problem, see for example [START_REF] Tropp | Just relax: Convex programming methods for identifying sparse signals in noise[END_REF].

We either consider the original minimization problem as stated in equation ( 12), or a modified minimization problem where the L0 pseudo-norm is replaced by the L1 norm. When considering the modified minimization problem involving the L1 norm, the optimal u S is no longer given by a simple expression as in equation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. Optimal values u S t,i are either zero, or if strictly positive, they must cancel the first partial derivative of the sum of the log-likelihood and the L1 norm, given below:

∂ ∂u S 2 log(u S + u B ) + v 2 (u S + u B ) 2 + λu S = - 2 v 2 (u S + u B ) 3 + 2 u S + u B + λ . ( 19
)
The only positive and real-valued root, as obtained by finding the roots of a third degree polynomial using a computer algebra system, is: with

u S (u B ) = t 1/3 + 4 9 λ 2 t 1/3 - 3 λ u B + 2 3λ , (20)
t = v √ 27λ 2 v 2 -16 3 3/2 λ 2 + 27λ 2 v 2 -8 27λ 3 . ( 21
)
We consider a numerical experiment to compare the detection performance of L0 and L1 formulations, for various contrasts between point-like scatterers and a piece-wise constant background, see figure 5a. The contrast between each scatterer and the background is given in figure 5b. 100 noisy versions are then generated using a Rayleigh multiplicative model. The L0 and L1 models are applied to each of these images with a fixed value of β (set to 0.05) and various λ values.

The receiver operating characteristic (ROC) curves are drawn on figure 6 and compared with the ROC curve obtained using the method based on a local estimation of the background proposed in [START_REF] Lopes | Structure detection and statistical adaptive speckle filtering in SAR images[END_REF]. The L0 model outperforms both the L1 model and the detection based on local background estimation. We show in figure 7 the scatterers detected by each model, for an identical proportion of correct detections. The L0 model gives a uniform performance (correct detections and false detections are well distributed in the image), while the L1 model fails to detect scatterers in the regions with higher background radiometry and gives more false detections in low-radiometry areas. This phenomenon is confirmed by representing the evolution of the probability of false alarm (Pfa) and the probability of detection (Pd) as a function of the background radiometry, figure 8. The L0 term produces constant Pfa and Pd (note that, in the numerical simulation, when the background radiometry changes, the scatterer radiometry is adapted so as to keep a constant contrast: u S = 40% u B ). In contrast, the Pfa and Pd of the L1 model decrease when the background radiometry increases, which is consistent with the non-uniformity of detected scatterers observed in image 7b.

V. AN APPLICATION TO CHANGE DETECTION

A. Method

In this section, we propose a direct application of our decomposition model for change detection, under the hypothesis that changes mostly affect the scatterers. This assumption is likely to be fulfilled in urban or industrial areas.

Given a time-series v, the model described in section II-C provides multi-temporal decomposition {(u B1 , u S1 , . . . , u Bn , u Sn )}, or in the case of a time-invariant background

where u St and u St are the two images extracted from the decomposition of the time series, δ is a shift to iterate over the neighborhood of the pixel i so C t,t (i) is the difference between the number of scatterers present in the window centered in i in one image and the number of scatterers in the other one. This map is then thresholded to locate the changes between the two images.

B. Results

To illustrate the results of our method, we ran this algorithm on a SAR time series from Saint-Gervais (France) acquired by TerraSAR-X satellite. This time series is composed of 26 single-look images. 13 of these images have been sensed in 2009, 13 in 2011. The ground-truth used for the quantitative results have been labeled manually by [START_REF] Su | NORCAMA: Change Analysis in SAR Time Series by Likelihood Ratio Change Matrix Clustering[END_REF]. We compare the results of this algorithm with other change detection methods in figure 9 and show an example of detected changes in figure 10. The ROC curve is obtained by varying the threshold applied on C t,t . Since all changes do not occur as scatterers changes, the ROC curves with our method do not reach high levels of probability of detection Pd. This simple thresholding approach achieves a performance comparable or superior to the methods described in [START_REF] Lombardo | Maximum likelihood approach to the detection of changes between multitemporal SAR images[END_REF] and [START_REF] Krylov | Change detection with Synthetic Aperture Radar images by Wilcoxon statistic likelihood ratio test[END_REF]. Change detection based solely on the detected scatterers however is less accurate than the recent work by Su et al. [START_REF] Su | NORCAMA: Change Analysis in SAR Time Series by Likelihood Ratio Change Matrix Clustering[END_REF]. Note that the simpler decomposition with a time-invariant background (ψ = TV ∞ 3D ) gives worse detection results. A reason for this worse performance may come from the degradation of the detection of scatterers when the background is not correctly estimated, as is the case of regions where the background undergoes some changes.

VI. CONCLUSION

This paper establishes the connection between a formulation of scatterers detection based on detection theory and the maximum a posteriori estimation with L0 pseudo-norm regularization. The problem of jointly detecting scatterers and estimating the background radiometry is then expressed as an image decomposition, stated as a combinatorial minimization problem. To solve this challenging large-scale combinatorial minimization problem, we propose to compute the minimumcut of a 4-D graph. Application to large areas and time series with numerous dates requires a limited memory algorithm based on block-processing with overlap. The comparison of the proposed L0-based decomposition with an L1-based decomposition shows a clear advantage for our formulation, which justifies tackling the combinatorial optimization problem. Our formulation gives constant false alarm rate and constant correct detection rate when the background level varies, in contrast to the L1 model. Finally, a simple application of the spatio-temporal decomposition is shown in the context of change detection.

Other decomposition models may be considered in the future, provided that they fulfill the constraints required to apply the graph-cut optimization method. Future work includes considering models designed for change detection rather than resorting to a post-processing step. The optimization step could be modified to use recent algorithms for min-cut computation (see [START_REF] Liu | Parallel graph-cuts by adaptive bottom-up merging[END_REF], [START_REF] Jamriska | Cache-efficient graph cuts on structured grids[END_REF]). The fusion approach proposed in [START_REF] Denis | Joint regularization of phase and amplitude of insar data: Application to 3-D reconstruction[END_REF] could also be adapted to include information from an optical image in the computation of the multi-temporal background component.

Fig. 1 .

 1 Fig.1. Graph construction used to solve the image decomposition minimization problem. (a) A set of qmax nodes (one for each quantization level) is created for each pixel of the multi-temporal image u B . Edges directed downstream are created, from the source S to the sink P. The edge in the minimum cut, shown with a thick line, indicates that the quantization level right below (colored in black) is the optimal value. (b) Bi-directional edges are added between corresponding quantization levels for all pairs of spatially adjacent pixels (edges from the source and to the sink that are not in the minimum cut are omitted in the drawing). (c) Bi-directional edges are also added between temporal neighbors. Only edges in the minimum cut are shown here, with the corresponding optimal level at each pixel denoted by a black disk.

MultiFig. 2 .

 2 Fig. 2. Decomposition of a time series of TerraSAR-X images of Saint-Gervais, France. Only the first and the last image of the decomposition are shown.

  (a) Solution obtained by local optimization (computation window of size 50 2 leading to a block effect) (b) Exact solution obtained using the presented algorithm for a computation window of size 150 2

Fig. 3 .

 3 Fig. 3. Crops of the background component of the first image of the Saint-Gervais set obtained with the presented method (filling window of size 50 2 ).

Fig. 4 .

 4 Fig. 4. Root mean squared error (RMSE) and memory used as a function of the size of the spatial computation window C, for the first image of the Saint-Gervais dataset.

Fig. 5 .

 5 Fig. 5. Input image used to compare L0 and L1 models.

Fig. 6 .Fig. 8 .

 68 Fig. 6. Receiver operating curve (ROC) between L0 and L1 version of our model and a method based on local window analysis [24].

Fig. 9 .

 9 Fig. 9. ROC curve of various change detection algorithms

  (a) Input image 1. (b) Input image 26. (c) Changes on image 1. (d) Changes on image 26.

Fig. 10 .

 10 Fig.10. Change detection results using the proposed method on images of Saint-Gervais series. Regions with changes that have been detected are circled in white.

in the absence of the smoothness term ψ, the decomposition is trivial:u S = 0 and u B = u B (ML) ≡ arg min u B ≥0i log p(v i |u B i) is the maximum likelihood estimate of the background.

ACKNOWLEDGMENT

The authors would like to thank the German Aerospace Agency (DLR) for the images (projects MTH0232 and LAN1746). This work has been funded by the CNES and the "Future et Ruptures" (Institut Mines-Télécom) program.