1.

Introduction

The area of Terra Nova Bay (TNB), Antarctica, has been of particular interest for the Italian scientific community in the last 30 years as the Italian Antarctic base Mario Zucchelli Station (MZS) is located in its vicinity. Thanks to the synthetic aperture radar (SAR) imagery of the area available today, a study of the dynamics of ice formation at TNB was planned for summer 2015. The European Space Agency (ESA) Third Party Mission program, which provides Cosmo-SkyMed (CSK) and Radarsat-2 (R-2) images, was used to carry out the study. A total of 14 CSK and 14 R-2 images for the period of February 20 to March 20 were retrieved; in addition, Sentinel-1 (S-1) images were also downloaded from the ESA scientific data hub. Wind data from the automatic weather station (AWS) “Eneide,” located in the vicinity of MZS, were obtained from the Italian Antarctic Research Program¹ as an ancillary data set.

The first survey of the SAR images revealed an unexpected surface feature, i.e., the presence of a prominent ice eddy. Ocean eddies have been found to play a major role in heat and salt transport² and in water mass transport.³ For example, discrete eddies may account for up to 60% of the eddy kinetic energy in strong eddying currents, such as the Antarctic circumpolar current (ACC).⁴ Field studies have also shown that discrete eddies have a major impact on biological productivity.⁵ A recent study⁶ provided an eddy census data set from a global ocean simulation over a period of 7 years and showed that the ACC is the region with the largest number of eddies in the world.

Satellite remote sensing provides a variety of platforms and sensors for monitoring sea surface and sea–atmosphere interactions and is therefore useful for observing mesoscale eddies, which may reach horizontal scales above 100 km and time periods from days to weeks. A marked variation of the sea surface temperature (SST) is one of the clearest indications of an eddy. Exploiting this feature, a neurofuzzy eddy detector was implemented with morphological descriptors, using the SST patterns of optical/infrared Nation Oceanic and Atmospheric Administration Advanced Very High Resolution Radiometer (NOAA AVHRR) images.⁷ In another study, the eddy shape was detected from the specific geometric structures of the rotating SST thermal wind velocity field.⁸

Vertical displacements of the sea surface are sensed by radar altimeters. Using the pattern of sea surface level anomaly (SLA), the eddy core area was detected with a set of specific delineation rules;⁹ in this case, the accuracy of the detection is limited by the heuristic function of the geometric criteria. Another proposed method combines the SLA gradient with a map of spatial markers for solving the Watershed transform.¹⁰ Although it can observe only large-scale singularities, the radar altimeter is widely used in oceanography. Combining SST and altimeter data with the geophysical knowledge of a particular site leads to a marked improvement in eddy detection.

In SAR imaging, the surface variations produced by the eddy’s revolving motion yield a modulation of the backscattered signal. Similarly, with a favorable wind speed, the surface films can enhance the imaged contrast of the eddy filaments. Tracking the signature of sea surface films can be performed by maximum cross correlation or by optical flow techniques;¹¹ in addition, the detection of slick patterns can be complemented by visual criteria.¹² In changing environmental conditions, the capacity of SAR to reveal eddy singularities by sensing the usual manifestations (sea surface roughness or marine surface films) is a challenging task. For this reason, synergic SAR modeling may include complementary data, such as SST imagery, radar altimetry, ground truth measurements, and chlorophyll observations.¹³

European satellites provide SAR images operating in C- and X-band: as these images are free, they were chosen for this study. Recent papers have shown that C- and X-band images in HH polarization are the most effective in sea-ice detection.^14,15 C- and X-band often display very similar distributions; but, when comparing their performances in discriminating the different kinds of sea ice, Dierking¹⁶ found that, only in the case of very thin ice, X-band is superior to C-band. Where the L-band is concerned, it was demonstrated its greater ability in detecting ridges and other sharp ice structures, but its lower sensibility to microscale ice structures.¹⁷ Sea-ice conditions usually have a dynamic evolution, a circumstance that may introduce ambiguities between different SAR frequency bands; in any case, Robinson¹⁸ estimates that no more than a few percent ($<3\%$ to 5%) of all ocean SAR images, captured in different bands, reveal substantial mesoscale features.

The first part of the paper deals with the detection of a time-varying ocean eddy in SAR images, exploiting a specific example of the singularity: rather than using the radar signal, the analysis is carried out by an indirect solution, which extracts the vortex extent outlined by the floating ice fragments of the scene. This is an alternative method compared with those used in previous surveys.

The task to be carried out is the identification of the eddy structure and the extraction of its main parameters: area and lifetime. In our case, from the sequence of SAR images in which the eddy was last seen on 9th March, its lifetime was about 20 days; indeed, from the image of March 15th, the ice eddy appears to be completely disrupted, and a well-defined polynya comes into sight. With regard to eddy area, a specific image processing scheme was developed. The scheme consists of the following stages: (1) nonlinear filtering, (2) segmentation, based on the Markov random field (MRF) theory, using a contextual approach applied both to the original and to the filtered image, and (3) extraction of the eddy area using an active contour detection algorithm, which works in an iterative manner.

The second part of the paper focuses on the role of the katabatic wind of TNB in generating the eddy singularities. It is well known that small-scale eddies, such as those produced at high latitudes, can be generated by atmospheric forcing;¹⁹ in this paper, the correlation between Eneide wind data and eddy extent is quantified as a joint function by means of two different statistical functions. The analysis provides elements for evaluating the temporal matching of the two variables.

2.

Methods

First, the SAR images were preprocessed using TeraScan software:²⁰ all the images were projected onto a common geographical reference with a pixel size of 100 m. Figure 1 shows one example of a georeferenced image. From each georeferenced image, a subset, or imagette, covering the area of TNB, was extracted; a pair of these imagettes are shown in Fig. 2.

Fig. 1

Georeferenced Sentinel-1 image of February 25, 2015, TNB, after being remapped onto an equidistant cylindrical projection.

Fig. 2

CSK imagettes covering the analyzed area of TNB, $800\times 600\text{pixels}$: (a) March 2 and (b) March 5, 2015.

2.2.

Segmentation

In the second stage, the scene is modeled by means of MRF theory. To define the conditional probability in terms of the Bayesian equation, the proposed algorithm analyzes the information derived from the gray-level data. With this scheme, the a priori model implements a homogeneity criterion based on the local label distribution. The segmentation result is obtained by a recursive stochastic minimization scheme. Before seeking the eddy boundary, the scene data require topological conformation, which means we need to have a clean space inside the eddy area, that is, we need to detect and remove the sea-ice objects located inside. This assignment is carried out by a hierarchical labeling process, defined by the topological properties of the ice fragments. Next, the eddy contour extraction task is performed. The geometry of the eddy is not regular and is usually defined by nonclosed contour segments. Classical edge detection algorithms are not very useful because of the complex contour distribution. To overcome this problem, an active contour detection algorithm is used. Moving in an iterative way, the active contour is placed according to both the gradient and the segmented image information. The final result is a binary mask, which provides the detection of the closed form of the eddy. In this way the final size of the eddy area is derived.

Let us consider two random variables, $X$ and $W$. In a maximum a posteriori (MAP) approach, let $X$ be the image to be segmented and $W$ be the segmentation result, $W=\{w_1,w_2,\dots ,w_M\}$, where $M$ is the number of labels of the sample space. The Bayesian rule is applied to define the equation of the segmentation model $P(W|X)=P(X,W)/P(X)$, where the joint probability term is $P(X,W)=P(X|W)P(W)$. As the total probability of $X$ is a common term throughout the process, a simpler expression is obtained by making this term equal to 1. The resulting expression is $P(W=w_i|X)\approx P(X|W=w_i)P(W=w_i)$.

The segmentation problem can be considered as a probabilistic labeling process. According to the Bayes theory, the MAP estimation equation can be written as $\widehat{W}=\arg \underset{x}{\max }[P(W|X)]\approx \arg \underset{x}{\max }[\log P(X|W)+\log P(W)]$, where $\widehat{W}$ is the label field to be estimated.

The Markov hypothesis makes the assumption that both fields $X$ and $W$ are expressions of MRFs, and, in this case, the Gibbs probability distribution²⁷ can be used to represent the terms of the Bayes estimation. The Gibbs distribution follows an exponential law; hence, the resulting Markovian estimation equation is

Eq. (2)

$$\widehat{W}=\arg \underset{x}{\min }[U(W|X)]\approx \arg \underset{x}{\min }[U(X|W)+U(W)],$$

where the probability that a given pixel $X=x$ belongs to a particular $w_i$ label is derived by minimizing the a posteriori energy function $U(W=w_i|X=x)$. In this paper, the conditional energy term $U(X|W)$ is estimated using the pixel gray-level information, and the a priori law $U(W)$ is achieved by a contextual analysis of the local label structure.

To define the information provided to the conditional energy term, the proposed algorithm considers both the original image $X$ and the filtered one $Z$. The morphological filtering performs an important noise smoothing task; but, depending on the topology of the structuring element $b$, thickening or removal of small objects may occur in connected regions, with a shifting of the image features. Some important information still remains in the raw image, and it can be used to complement the segmentation process. An interesting alternative to the combined use of information sources arises when one considers the probability property of independence²⁸ as being conditioned by the field of labels. The probability of the joint gray-level information of pixels $X_i$ and $Z_j$, given the label $W_k$, is obtained by the product rule:

$P(X_i,Z_j{|W}_k)=P(X_i{|W}_k)P(Z_j{|W}_k)$. Including the former property in the Markovian model, Eq. (2) is approximated as

Eq. (3)

$$\widehat{W}=\arg \underset{x}{\min }[U(W|X,Z)]\approx \arg \underset{x}{\min }[\frac{{\lambda }_{1}U(X|W)+{\lambda }_{2}U(Z|W)}{{\lambda }_1+{\lambda }_2}+U(W)],$$

where $U(Z|W)$ is the conditional energy function of the gray-level $Z$ occurring given the $W$ label. This equation involves the conditional random field of the two-level model. In addition, ${\lambda }_1$ and ${\lambda }_2$ are the weight coefficients of the conditional membership functions. With these coefficients, it is now possible to assess the energy contributions of images $X$ and $Z$.

To define the a priori energy function $U(W)$, the Ising model is used; the Hamiltonian is then given by $U(w_i,w_s)=\{\alpha \text{if}{w}_i\ne w_s,-\alpha \text{if}{w}_i=w_s\}$, where $\alpha$ is the strength of interaction of the ferromagnetic model and $U(w_i,w_s)=-\alpha$ when pixels at coordinates $i$ and $s$ of the segmented image $W$ have the same label $w$. To define homogeneous regions, the minus sign supports the minimization process in Eq. (3). Binary segmentation is required to detect the scene elements. The MRF model provides the required binary field by labeling the set of pixels as ice or nonice (open sea). The segmentation result is obtained by a simulated annealing scheme.²⁹ Using ${\lambda }_1={\lambda }_2=0.5$, the convergence is achieved by applying a number of 35 iterations, with a strength parameter $\alpha =0.3$.

As an example of the preprocessing steps, Fig. 3 shows the results obtained with the nonlinear filter (a) and using the binary Markovian segmentation (b). These results were obtained with the CSK image of March 5th [Fig. 2(b)].

Fig. 3

Preprocessing of a CSK imagette: (a) filtered image and (b) Markovian segmentation result.

2.3.

Toward Eddy Region Bounding

The MRF segmentation provides a binary result, where the background of the SAR scene is displayed as black pixels for the open-sea regions and as white pixels of the detected sea-ice objects. Some dynamic characteristics of the wind can be seen in the spatial distribution of the sea-ice elements; this wind signature is easily detected as these elements form a sort of ring. The sea-ice objects do not define a perfect circle (see Figs. 2 and 3) as its contours are nonstationary. To complete the scene description, we must deal with a number of minor sea-ice fragments, which can be seen inside the eddy area; they are not part of the eddy topology and can lead to spurious detections. Before completing the detection and the extraction of the eddy contour, it is important to remove these sea-ice fragments from the scene; to achieve this goal, the ice objects are classified using a labeling method.

Graph theory³⁰ provides useful concepts for describing the objects of a scene. A graph $G$ is defined as the set $G=(V,E)$, where $V=\{v_i,i=1,\cdots ,N\}$ is a set of “vertices,” and $E$ is a set of “edges.” In our case, the set of vertices refers to the lattice structure of the image, and the set of edges is formed by linking pairs of pixels of the set $V$. One way of characterizing the edge features is done by describing the “adjacency” and “connectivity” of set $V$. Two pixels are four-adjacent when they share a side. A set of pixels is defined as connected when it is not divided by a boundary. A region is called four-connected or ${\nu }_4$ if a finite chain links every pixel to a set of consecutive four-neighbor pixels. The ${\nu }_4$ provides an important relationship for labeling segmented scenes. An “ordering list” is arranged for labeling the four-connected regions,³¹ and the number of pixels of each object $S_i$ is obtained. A percentage of the image size is used to fix a threshold $Tr$. When $S_i<Tr$, the related sea-ice object is removed. Figure 4 shows the result of applying this process to the CSK imagette of March 5th. A spotless area is observed inside the eddy [Fig. 4(b)], and the scene is now suitable for deriving the eddy boundary.

Fig. 4

CSK imagette of March 5. The eddy contour detected and outlined: (a) the remapped SAR image and (b) the binary image with eddy area free of sea-ice objects.

3.

Results

The eddy area values are shown in Fig. 5. They were obtained using the processing scheme described above as applied to the sequence of imagettes extracted from CSK, R-2, and S-1 images.

Fig. 5

Eddy area evolution obtained using the processing scheme described above.

The interaction between variables can be computed by measuring the degree of relationship with the assumption of linear or nonlinear dependence. To reveal the coupling between the forcing wind and the induced eddy area, with an alternative perspective, we preferred to analyze the case of a linear relationship by means of two formal statistical correlation functions. The examination is derived from the time series analysis used to provide explicit indicators of the dependence on wind field. To quantify the linear interaction between geophysical processes and meteorological variables, the running correlation coefficient $R_{\mathrm{cc}}$ function³⁹ is simply the adoption of the second-order statistic moment of a bivariate random variable. Using local windows overlapped in time, the $R_{\mathrm{cc}}$ function is the first analysis method.

Let $N$ be the total data length; $X_k$, the wind speed; and $Y_k$, the measured eddy area; the $R_{\mathrm{cc}}$ of the two centered signals, given the $i$’th step with a window length $(-n,n)$, is computed using both first- and second-order moments of the marginal variables

In the second analysis, we propose to apply the Spearman correlation coefficient,⁴⁰ which is a nonparametric rank statistic measure of the connection between variables. One important feature of this function is that it is less affected by single influential observations (outliers). The central information provided by this function is the monotone dependence between variables. Using the same moving window of the $R_{\mathrm{cc}}$ function, a high correlation value means that both variables have corresponding increasing or decreasing distributions. This may seem a simple feature; but, it is another element that demonstrates the association involving the two time series. The “monotonic relationship” of the Spearman correlation complements the “linear dependence” displayed by the $R_{\mathrm{cc}}$ function.

The AWS Eneide, located on TNB near MZS, collects the wind channeled by the mountains of the inner continent (see Fig. 1) and provides real-time and archived meteorological data. Wind speed and direction of Eneide were composed in one single figure defined as “effective wind field,” i.e., the wind component pushing eastward or at 270 deg. The acquired data are transmitted by the Argos satellite every hour and received in Italy. For the 18-day study period, February 20 to March 9, the collected wind field is shown in Fig. 6 (each value is the average over the previous 6 hrs).

Fig. 6

AWS Eneide wind field collected over a time period of 18 days.

The final step is the statistical inference of the dependence between the two variables. The pairwise linear correlation function was analyzed by Eq. (5) using bounded temporary windows. This function makes it possible to explain the interaction of the geophysical variables in discrete time series. Using a time lag of $n=2$, $R_{\mathrm{cc}}$ is computed with windows of 5 units (time segment of 30 hr); the result is shown in Fig. 7(a). The result of the Spearman correlation is shown in Fig. 7(b). Both variables, eddy extent and wind speed, can be considered as random variables, but a joint connection is observed. Figure 7 shows that 57% of the coefficients have a correlation $|R|>0.6$, which, using statistical criteria, confirms both a strong linear relationship and a monotone association between the two variables.

Fig. 7

Correlation functions between eddy extent and wind speed: (a) the running correlation coefficient and (b) the Spearman correlation coefficient.

By statistical criteria, a very strong correlation arises when $|R|>0.8$. Figure 7 shows that only 26% of the paired data are very strongly correlated. The probability theory states that two random variables are independent when the occurrence of one variable does not affect the manifestation and the evolution of the second one, and, consequently, the variables are uncorrelated. Given that strong and very strong correlations are widely distributed, our results show that the analyzed variables are not statistically independent. It is well known that eddies can be induced by other geophysical events, such as the interaction of sea currents with undersea topography, or by convection caused by the sea surface gradient temperature. Related information of such events was not available for this study; nevertheless, the correlation sequence reveals a dependency between wind field and eddy extent variation.

4.

Discussion

Identifying eddy singularities in SAR images is only a partially solved problem given the complex signature of the sea/atmosphere modulation pattern. For this reason, the analysis of the scene features requires a specific processing procedure as suggested by the present paper. Based essentially on the Bayesian estimation and on the MRF theory, an original aspect of the proposed scheme consists of a methodology to parametrize the eddy structure, to detect its temporal evolution in the time window of 18 days. The advantage of our proposal is that the contextual dependence of the segmentation process model is able to deal with the noisy behavior of the images from different satellites (CSK, R-2, and S-1). The resulting scene characterization is provided by merely binary labeled fields, see Fig. 3(b); this seems a simple outcome, but a functional product cannot be obtained by conventional segmentation methods. Two objective correlation measures are used to prove the correspondence between the geophysical variables involved. Both ocean and atmospheric dynamics impact the floating ice distribution, and this introduces discrepancies in the traced eddy contour. The detection of ice fragments floating in the surrounding environment, being configured by a set of segments, is not very precise and yields some divergence in the correlation progression, see Fig. 4. A clear limitation of the correlation measure is that the statistical interpretation is only a “linear” one, whereas a nonlinear analysis is beyond the scope of this paper. The recursive structure of the MRF algorithm requires a significant amount of processing time: using MATLAB^® and computing with an Intel Pentium^® CPU G2020 of 2.9 GHz and 6 GB RAM memory, the elapsed time for a single scene is about 1150.1 s.

5.

Conclusions

Eddies play an important role in the dynamic distribution of heat, mass, and other geophysical variables, which operate at high latitudes. This paper presents a complex processing scheme developed to model small-scale eddies generated by atmospheric forcing. First, a nonlinear filter is applied to a sequence of imagettes extracted from the SAR images; the filter task is the smoothing of the speckle degradation while preserving primary details. The filtering stage is complemented by a segmentation stage based on MRF theory. Because of the speckle effect, the use of both filtering and MRF stages provides a contextual analysis, which improves the segmentation result.

Two regions of the imagettes are considered for the purposes of this paper: the open sea and the sea-ice, both regions being delimited as objects, and with the floating ice providing basic information for tracing the eddy bounds. The set of sea-ice objects of the binary segmentation is labeled to remove drifting sea-ice objects from the internal eddy region. Using this criterion, small objects, both inside and outside the eddy, are rejected, even though the primary interest is to remove those inside the eddy. In this way, a clear definition of the eddy boundary is obtained.

The use of an active contour algorithm enhances the overall quality of eddy extent detection. The processing scheme presented in this work, when applied to the set of SAR images of TNB, produces a database of area variation in the study period. The link between eddy area and wind field was analyzed by means of the running correlation coefficient function ($R_{\mathrm{cc}}$) and the Spearman correlation coefficient, which revealed the soundness of the linear relationship and the monotone dependence between the two variables.