SAR image observations of the A-68 iceberg drift

A methodology for examining a temporal sequence of Synthetic Aperture Radar (SAR) images as applied to the detection of the A-68 iceberg and its drifting trajectory, is presented. Using an improved image processing scheme, the analysis covers a period of eighteen months and makes use of a set of Sentinel-1 images. A-68 iceberg calved from the Larsen C ice shelf in July 2017 and is one of the largest icebergs observed by remote sensing on record. After the calving, there was only a modest decrease in the area (about 1%) in the first six months. It has been drifting along the east coast of the Antarctic Peninsula and 5 it is expected to continue its path for more than a decade. It is important to track the huge A-68 iceberg to retrieve information on the physics of iceberg dynamics and for maritime security reasons. Two relevant problems are addressed by the image processing scheme presented here: (a) How to achieve quasi-automatic analysis using a fuzzy logic approach to image contrast enhancement, and (b) Adoption of ferromagnetic concepts to define a stochastic segmentation. The Ising equation is used to model the energy function of the process, and the segmentation is the result of a stochastic minimization. 10

present, it is the largest iceberg in the world. Because of its size, an iceberg like A-68A can have a life of several years. Iceberg drifting patterns constitute a risk for navigation and shipping routes. Satellite remote sensing imagery can provide the tool for 25 mapping iceberg trajectory progression.
In iceberg monitoring by remote sensing there are two basic objectives: iceberg detection and iceberg drifting forecast. For iceberg detection, a hierarchical object-based segmentation is applied to a set of geometrical parameters of ENVISAT/ASAR images (Mazur et al., 2017). The radar altimeter is an alternative instrument and in (Tournadre et al., 2016) the signatures of icebergs in waveform space are analysed by threshold criteria to parametrize iceberg distribution. In (Scheick et al., 2019), a 30 machine learning technique is applied to mask clouds in multispectral Landsat images. Then, iceberg detection is performed by threshold criteria being careful to notice the radiometric contrast between icebergs and the surrounding open sea. Using Sentinel-1 SAR and CryoSat-2 SIRAL data, Han et al. (Han et al., 2019) describe the topological evolution of iceberg A68 and investigate the effects of environmental forces over a period of 18 months. A review of the remote sensing of the cryosphere and processing techniques for sea ice can be found in (Meier and Markus, 2014;Zakhvatkina et al., 2019). For more than 45 35 years passive microwave images have been used for monitoring polar regions. The natural molecular interaction of the scene elements produces an electromagnetic radiation which can be used by passive microwave (PMW) sensors to discriminate the electromagnetic signatures of water, snow cover and ice extent (Thomas, 1986). Main applications of PMW images are sea-ice concentration analysis (Ivanova et al., 2014;Kern et al., 2019), thin ice studies (Mäkynen and Similä, 2015), sea-ice extent and ice-edge location (Meier and Stroeve, 2008), and sea ice production (Preußer et al., 2019). However, the emitted radiation is 40 very low, and consequently the passive energy must be compiled over large regions. With a swath width of 1450 km, the spatial resolution of the Advanced Microwave Scanning Radiometer 2 (AMSR2) ranges from 5×3 km to 62×35 km. Thus, PMW instruments are appropriate to observe large regions and can give only a rough estimation when applied to iceberg monitoring.
For iceberg and ice tracking forecasting, an unmanned aerial vehicle platform was used to analyse thermal video (Leira et al., 2017). A set of dynamic forecasts was obtained using GPS trackers positioned on icebergs in (Yulmetov et al., 2016). Based on 45 Sentinel-1 images in (Demchev et al., 2017) non-linear diffusion filtering reduces the speckle noise, and features are detected in a non-linear multiscale space representation; nearest-neighbour matching reveals the connections between the extracted features, these being the basis for sea ice drift tracking. In another study (Muckenhuber et al., 2016), the Sentinel-1 SAR image resolution is reduced by a spatial average operation to decrease speckle influence. Then, sea ice tracking is performed using a scale-invariant feature transform algorithm. In a set of ENVISAT/ASAR images, after morphological characterization by pixel-50 based segmentation, tracking is performed using ocean current data (Collares et al., 2018). In (Wesche and Dierking, 2016), a drift model makes use of wind predictions for estimating positions and trajectories of icebergs observed in ENVISAT/ASAR images. More complete delineations, such as the statistical, kinematic and dynamic models, require hydro-meteorological data and both atmosphere and ocean circulation models (Diansky et al., 2018). In general, modelling the interacting forces is a very complex task (Andersson et al., 2016;Bigg et al., 2018).

55
With regard to the image processing domain, SAR reconnaissance capabilities are limited by the peculiar behaviour of radar imaging; indeed, basic problems, such as the irregular image contrast and the multiplicative degradation by speckle noise are still a challenge. Pixel-based techniques, such as K-means, Fuzzy C-means, minimum distance criteria and normalized multi-2 https://doi.org/10.5194/tc-2020-180 Preprint. Discussion started: 21 July 2020 c Author(s) 2020. CC BY 4.0 License. band indexes are well suited for optical and multi-band images, but their algorithmic performance is limited by the random nature of the SAR data (Maître, 2010).

60
For this reason, in this paper, the stochastic process theory is taken into account. For modelling the spatial interaction of pixel data, a model based on concepts of statistical ferromagnetism appears promising. Two relevant problems are addressed by our image processing technique: (a) Low-level fuzzy logic image contrast enhancement, which was derived from medical image analysis, and (b) A segmentation algorithm which considers the random behaviour of the SAR imagery for merging contextual data. A processing scheme was then implemented which consists of the following steps: (1) Contrast enhancement;

Material
This study is based on a set of twelve Sentinel-1 Extra Wide Swath Ground Range Detected (S1 EW GRD, 400 km swath, 20 x 40 m spatial resolution) images at Level-1 in HH polarization; the images were acquired from 22 July 2017 to 26 January 2019 and their geographical coordinates range from latitude 66 • S to 69 • S and from longitude 57 • W to 63 • W. After retrieval 70 from the ESA Scientific Data Hub, the images were remapped onto a regular grid in stereo-polar projection with a pixel size of 200×200 m. The scene size is 400×400 km. Figure 1 shows the image corresponding to 22 July 2017, just a few days after the calving event. In Antarctica, most icebergs are created by the calving of ice shelves and glacier tongues. The flat plateau top appearance is a characteristic feature of the tabular icebergs produced in this region.

75
Electromagnetic variables of radar may introduce undesirable effects in the radiometric quantization so that the grey-level distribution displays a histogram with saturation in local ranges. The subsequent effect is poor image contrast which reduces perception capabilities. Some images in the analysed data set display this characteristic, and, for this reason, intensity transformation was included in the analysis.

SAR scattering 80
Remote sensing by SAR systems is the result of a complex electromagnetic phenomenon and the radargrammetry technique must consider adverse variables, which may affect the function of the imaging system (Leberl, 1990). The physical manifestation of radar reflectivity is the scattering phenomenon. Diffuse and specular reflections are due to the geometric irregularities of the surface. Other electromagnetic properties, such as the dielectric constant, permeability and conductivity complement the scattering models. These properties modify the rate of the incident and reflected energy. Therefore, the backscattered signal 85 determines the radiometric signature of the scene elements.
At high latitudes, the properties of the scene elements change with time. Geophysical and climatological variables, such as the temperature of the medium, wind speed, rain, salinity and humidity introduce dynamic fluctuations in the scattering phenomenon (Xu et al., 2018). For example, new ice produces specular reflections like a thin film with a smooth appearance.
Dry snow produces a very weak reflection. Ice sheets and dry soil have a similar dielectric property (Richards, 2009). However, 90 moisture in snow, salt and air in the ice layers all increase diffuse scattering. The variations of the dielectric property can cause two main problems in SAR image analysis: an irregular contrast of the scene elements and an overlap between the modes of the intensity histogram. These problems must be solved in the feature extraction/segmentation process and a solution to both is proposed in this section.

95
Several contrast enhancement methods are currently applied in image remote sensing analysis. The main options are based on linear and non-linear formulations such as equalization, normalization, matching, logarithmic and exponential transformation functions (Solomon and Breckon, 2011). In the spatial domain, one important parameter is the dynamic range, which is defined by the smallest and the largest grey-level value of the image under analysis. To obtain an improved mapping of the grey-scale distribution, the basic approach is to transform the dynamic range, a task which can be accomplished by fuzzy set theory 100 (Nachtegael et al., 2003). In the case of the fuzzy histogram equalization (Kamel and Campilho, 2005), a membership function µ(g) is defined for each pixel grey-level value g nm at the spatial coordinates (n,m) and this is expressed by: The terms g max and g min are the maximum and the minimum values of the grey-level domain. The parameter used is the dynamic range, which is the normalization term of Eq.
(1). The function µ(g) is interpreted as a homogeneity operator of the 105 input image luminance, or as an adapted measure of the biological perception of contrast (Pratt, 2001). Both the input image X and µ(g) have the same matrix rank.
In the original formulation of the fuzzy sets (known as type-1 fuzzy sets or T1 FSs), the inferred information is structured by membership functions, which are a representation of the probability density function. A limitation of the fuzzification process is that the membership functions are deterministic for a given random variable but usually the histogram of an image exhibits 110 mixed random variables, and this uncertainty requires additional abstraction. The next step of our scheme makes use of type-2 fuzzy sets (T2 FSs): as their membership functions become fuzzy (Mendel, 2007), we obtain a better representation of the uncertainty and the information ambiguity of the inferred probability density function.
In this paper, the T2 FSs were the choice for implementing a contrast enhancement algorithm. Equation (1) where α is a fuzzifier parameter with 0 < α < 1, and µ up and µ low are the upper and lower bounds of the T2 FSs membership function (Tizhoosh, 2005). The µ up and µ low functions and the input image X have the same matrix size. The fuzzy function 120 maps the input image into a grey-level transformation, and this implies multi-criteria decision making. One suitable option for global decision making is the t-conorm operator, and in this paper, the adopted algebraic operator was derived from medical image processing literature (Chaira, 2019): whereX is the expected value of the input image X, and µ up and µ low are the T2 FSs membership functions obtained by 125 Eq.
(2). The process begins with the ingestion of the input image X into Eq. (1), afterwards Eqs. (2) and (3) are computed. The membership functionμ(g) maps the contrast enhancement operation.

Stochastic segmentation approach
SAR images are affected by multiplicative speckle degradation; therefore, even binary segmentation is not a simple task. In an elementary polar environment description, the analysed scene is a binary field composed of open sea and ice sheet objects. In the framework of Bayes' theory, the implementation proposal infers the relevant information from both pixel-based and locally connected pixels. The input image X is a pixel lattice S of N×M, where the pixel coordinates (i, j) are structured by a neighbourhood system η. According to the Euclidean distance, the first and second-order system, η 1 and η 2 , correspond to the 4-connected and to the 8-connected systems, respectively. A clique is a subset C ⊂ S and it represents the primitive image structure of connected pixels or sites. For a system η 2 , the associated set of cliques C 1 and C 2 are, respectively, the central 135 pixel X ij and the set of pixel pairs. The spatial feature field defines a set of n mutually exclusive labels L = {l 1 , l2, · · · , ln}.
The output of the segmentation process is the variable Y = { y ij ∈ L}.
With the term P (X|Y ), the Bayesian theory takes into account the probability distribution of the pixel grey-level, given the label field Y and also the "a priori" information of the labelling process, i.e. the term P (Y ). A Bayesian maximum-a-posteriori (MAP) estimator is: where Y l ∈ L, and P (X) is the probability of realization of the input random variable. Using a Markov random field (MRF) model, the probability terms of the MAP equation can be adapted to introduce contextual information. Once the random variable X is assumed as a MRF realization, a Gibbs function models the region process of Y . Thanks to this concept, the terms of Eq. (4) are approached by the sum of energy functions U ≈ U (X|Y ) + U (Y ).

145
The term U (X|Y ) is considered a realization of the label set in the grey-level range, and, in this paper, the conditional modes are expressed by Gaussian functions U (X|Y ) = ln( √ 2πσ i ) + (X −Ȳ i ) 2 (2σ 2 i ) −1 , whereȲ i is the mean, and σ 2 i is the variance of the label Y i . The MRF theory is based on statistical physics (SIGELLE and RONFARD, 1992), and in this paper for introducing the function U (Y ), the Ising model (Ibe, 2013) was implemented following a ferromagnetic interpretation of the random process. The cardinality of the sites σ i is specified through the local label arrangement of Y : In a ferromagnetic reading, α is a characteristic of the involved element, M is a supplementary magnetic field, and β is the magnetic condition of the material. The effect of M is to induce alignment of the ferromagnetic elements in the direction of the field of M . The β parameter indicates the interactive magnetic forces of adjacent sites. The magnetic attractive case occurs when β > 0. The joint effect of M and β is to produce states of low energy, and in the case of a segmentation process, to 155 generate homogeneous label configurations. Thus, the resulting U function is driven by: To find the optimal estimate of the label field L, a numerical minimization of U is needed. As Eq. (6) is a non-convex function displaying different local minimal energy states (zero slop intervals), in order to induce progressive low energy configurations, a simulated annealing scheme was implemented. To obtain further adjustments in the local energy array, thus allowing to reach 160 a global minimum state, the Gibbs sampler criterion (Chatelain et al., 2011) was applied. A homogeneous grouping of the pixels is obtained at the end of the recursion.

Fuzzy Contrast Enhancement
In order to evaluate its performance, the applied fuzzy algorithm was compared with alternative contrast solutions: (a) the  Fig. 3; Fig. 3(a) shows the iceberg detection of the image of 13 December 2017 while Fig. 3(b) the result corresponding to 18 January 2018. The detected iceberg shape is displayed in white and, for a better visualization, the whole SAR scenes are used.  Figure 4 shows the variation of the area on the left of y-axis and the perimeter on the right of the y-axis. To appreciate the long-term tendency, a curve fitting function was used: for both parameters a 5th-degree polynomial curve fits the series of data points. A decay tendency is observed in both area (blue curve) and perimeter (red curve) parameters.
Two complementary parameters are the major axis length and the rotation angle. The binary pattern of the detection is the basis for expressing the iceberg shape as a polygon with specific properties. The geometric centroid (centre of mass) of the 200 connected iceberg pixels, i.e. concurrency point of both the major and minor axes, was derived. The two axes are marked as the segments AB and CD in Fig. 1. Taking as reference the horizontal axis and the segment defined by the centroid and point A, the rotation angle is computed in a counter-clockwise sense. Figure 5 displays the rotation angles derived from the first and  the last analysed images. Figure 6 shows the time evolution of both the rotation angle and the major axis length parameters.

Discussion
The multiplicative nature of speckle degradation introduces spurious pixel grey level values, and this statistical confusion is a basic difficulty for SAR image segmentation. To address the random nature of the SAR data, two probability abstractions provide the required information: a contextual second-order neighbourhood model and a pixel-based analysis. Therefore, the segmented field is the result of a double segmentation model, and it was implemented using numerical optimization. Two train-210 ing windows are manually fixed to derive the mean values of the ice sheet and non-ice sheet objects. This is done for each image  Fig. 6, point 10), and the speed increased to 16.8 km/month. By January 2019, the angle of the major axis was about 250 • . Using the centroid data, the total displacement distance was 220.6 km. In the analysed period, a slight reduction in the planar shape parameters was observed. The visible iceberg area reduced by about 3.7% and the major axis 225 length by 3.9%. Melting, breakup and forced motion are consequences of the iceberg-environment interaction; main driving force arises from the surrounding ocean with some atmospheric contributions. Large icebergs last for several years and the gravitational force may introduce topological changes. The gravitational force pushes outward the iceberg mass, and, over the years, the cumulated effect produces a decrease in thickness and an increase in iceberg length (Bigg, 2015). The influence of these elements is out of the scope of this paper. In the last analysed image, the A-68A iceberg was approaching the marginal 230 zone of the Antarctic Circle. At this point, the coastal current is expected to be the driving force of its displacement. Moving in the direction of the Scotia Sea, the iceberg must still travel about 400 km to reach the northernmost point of the Antarctic Peninsula.

Conclusions
A methodology is proposed for the analysis of a temporal sequence of SAR images. Two fundamental problems in the remote 235 sensing domain are the irregular image contrast and the mixed multimodal class distribution. This paper takes into account the image uncertainty for proposing the combined use of fuzzy logic and of a ferromagnetic approach which models overlapping class intervals. A pre-processing stage implements a fuzzy contrast enhancement in the spatial domain. In the fuzzification process, a set of image features define the membership functions whose domain and range are a rough fit to the image feature histogram. The concepts of ferromagnetic theory were chosen to define a stochastic segmentation method. In ferromagnetic 240 theory, the effect of an external magnetic field is to induce alignment of the ferromagnetic elements; this, in the case of a segmentation process, simulates the magnetic attractive force by generating local homogeneous pixel configurations. The Ising model and the Bayes equation were the basis for implementing the spatial pixel interaction. The derived binary field is the result of a stochastic minimization process. Because of the scene size and of the recursive nature of the optimization algorithm, the computational requirements of the MRF segmentation are computation-intensive. The final analysis result shows the movement 245 of the A-68A iceberg over a time period. Due to its colossal size, small variations in area, perimeter and major axis length parameters were observed. Up to 26 January 2019, the detected area was 96% of its original size. The surrounding ice in the winter season, wind patterns and sea-floor elevation layers cause irregular displacements and variant iceberg velocities, but the dominant direction seems to be towards the Eastern Weddell Sea. The main contribution of this paper is in the image processing domain with an application to the tracking of the A-68A iceberg. Ancillary information such as meteorological data, 250 ocean currents, wind speed, temperature and geomorphology of the seabed was not available for this study, but the proposed methodology can be integrated to perform dynamic modelling. Competing interests. The authors declare that they have no conflict of interest