ESTIMATING SEA ICE PARAMETERS FROM MULTI-LOOK SAR IMAGES USING FIRST-AND SECOND-ORDER VARIOGRAMS

The spatial structures revealed in SAR intensity imagery provide essential information characterizing the natural variation processes of sea ice. This paper proposes a new method to extract the spatial structures of sea ice based on two spatial stochastic models. One is a multi-Gamma model, which characterizes continuous variations corresponding to ice-free area or the background. The other is a Poisson line mosaic model, which characterizes the regional variations of sea ice with different types. The linear combination of the two models builds the mixture model to represent spatial structures of sea ice within SAR intensity imagery. To estimate different sea ice parameters, such as its concentration, scale etc. We define two kinds of geostatistic metrics, theoretical firstand secondorder variograms. Their experimental alternatives can be calculated from the SAR intensity imagery directly, then the parameters of the mixture model are estimated through fitting the theoretical variograms to the experimental ones, and by comparing the estimated parameters to the egg code, it is verified that the estimated parameters can indicate sea ice structure information showing in the egg code. The proposed method is applied to simulated images and Radarsat-1 images. The results of the experiments show that the proposed method can estimate the sea ice concentration and floe size accurately and stably within SAR testing images. * Corresponding author


INTRODUCTION
Sea ice plays a very important role on hydrology, atmospheric thermal cycle, ocean currents and ecosystem in high latitudes, especially the Polar Regions (Zhu et al., 2012;Emilio et al., 2010).Sea ice limits shipping in high latitudes and requires special attentions when planning and running offshore activities in the Arctic Ocean and its peripheral seas (Johannessen et al., 2007).Consequently, timely obtaining the information of sea ice in a wide range is becoming more and more urgent for many governments and organizations.With the development of aerospace and aviation technology, remote sensing has proved to be a powerful tool for monitoring sea ice.Compared with the optical remote sensing, synthetic aperture radar (SAR) systems have been shown to be very useful due to its capability of acquiring data under all weather condition, during day and night, and are widely applied to marine , forestry, environment, agriculture fields (Kerman, 1999).
Estimating parameters of sea ice such as its concentration, size, age from SAR imagery is always a difficult task.Currently, many methods have been proposed to extract sea ice information.Nystuen and Garcia (1992) used both texture and standard statistics to classify sea ice types, where both of them are derived from a grey-level co-occurrence metrics (GLCM).
In their experiments, they found that in some special ice zone, because of the movement of sea ice, a kind of new ice, called "odden" and multiyear ice, are of a SAR backscatter property and difficult to be distinguished.Soh and Tsatsoulis (1999) used GLCM to quantitatively evaluate textural parameters and determined the optimal parameters and representations for mapping sea ice texture.They investigated the quantization, displacement and orientation factors of GLCM on SAR sea ice imagery.They concluded that matrices with a range of displacement values are more adequately representative of a texture than a single χ 2-optimal matrix.Because of the movement of ocean currents and winds, different types of sea ice may have the same backscatter and vice verse.These phenomena will lead to incorrect classification of sea ice.Some researchers (Clausi & Deng, 2004;Ochilov & Clausi, 2012) investigated operational segmentation and classification methods by comparing three kinds of texture feature extraction methods and four types of segmentation models.After that, different extraction methods are combined with different segmentation models to do experiments.The gamma mixture model followed by a MRF label model is able to demarcate ice from open water region.All above studies used the traditional statistical metrics and considered the autocorrelation of sea ice only.Although the ice age is derived only from texture structure, it is not so accurate, and there is so much unknown information about sea ice.This paper considers both the autocorrelation of sea ice and the correlation between sea ice and the open water.Spatial statistical metrics are used to estimate the basic parameters of sea ice, such as total concentration and floe size of sea ice.
According to the properties of SAR imagery and sea ice spatial structure, we build two stochastic models, one is multi-Gamma model and the other is Poisson mosaic model.Then two kinds of geostatistic metrics are used to measure the spatial structure of sea ice.Theoretical variograms are built based on the two stochastic models, and the experimental variograms are calculated from the images.Then fit the theoretical variograms to the experimental variograms and the parameters are estimated in the theoretical variograms.
The rest of the paper is organized as follows.Section 2 describes the spatial stochastic models used to model the spatial structures.In Section 3, the proposed method is validated on simulated images and then applied to the actual RADARSAT-1 images of the Ungava Bay of Canada.The paper concludes with a summary of findings and suggestions for future enhancements in Section 4.

Random function and variogram
From spatial statistics point of view, the spatial structure and variability revealed in SAR data can be modeled by a regionalized variable z(x), which is simply considered as the realization of a random function Z(x) constructed at all points x of a given region D ⊂ R 2 , while Z={Z(x) : x ∈ D ⊂ R 2 } defines a random field on D.
The SAR data have been translated into general datum for analyzing.A stationarity condition should be set up first.Let Z(x) be the second-order stationarity random function, then it is of the following properties, where m is the expectation value of Z(x) and C(h) is the covariance function for all pairs of pixels x and x+h.When h = 0, he covariance function divided by the variance is called the correlation function Under second-order stationarity assumption, the second-order variogram γ 2 (h) of Z(x), which describes the variability between the intensity values of two pixels separated with h, is From Eqs. (1) and (2) (Schabenberger & Gotway, 2005), In isotropic stochastic processes, the covariance is a function of the Euclidian distance h=||h|| only, that is, C(h) = C(h).σ 2 is the variance of Z(x).As a result, Eq. ( 5) can be simplified as follows, Fig. 1 shows the relationship between C(h) and γ 2 (h).R is the largest distance between two relative positions; Sill is the limitation value of γ 2 (h), its equal to C(0).When h exceeds the range, C(h) tends to zero and γ 2 (h) is stabilized.It means that data separated by a distance larger than the R are uncorrelated.The R is an important parameter related to the spatial scale of the data (Feng, 2008).It can be used to reflect the floe size of sea ice, so how to estimate R from data (images) accurately is a significant part of this method.
For a given random function Z(x), under second-order stationarity assumption and in isotropic stochastic processes, its first-order variogram γ 1 (h) is defined as follow (Chiles & Delfiner, 1999) where γ 1 (h) is seen as first-order moment.When h approaches its range R, the value of γ 1 (h) tends to be constant.The first-and second-order variograms will both be used to characterize the spatial structures and their validities and accuracies will be compared in the following sections.

Stochastic models
In order to characterize the sea ice spatial structures, SAR intensity images are considered as a combination of two stochastic second-order stationary models.According to the SAR data properties (Keith et al. 2006) and the distribution properties of sea ice, we build multi-Gamma model Z g (x) and Poisson mosaic model Z m (x).The former one is used to characterize continuous variations corresponding to the background of sea ice, the other is a tessellation model in which the image domain is randomly separated into non-overlapping cells to simulate the real sea ice areas.(Chilès & Delfiner, 1999).

Bivariate Gamma model Z g (x):
A random function Z g (x) is said to be a bivariate Gamma random function if the vector Z = (Z g (x), Z g (x + h)) for x, x + h ∈ D is distributed according to a bivariate Gamma distribution.The bivariate Gamma distribution on R 2 has been defined in several forms.Most of them exploit various properties of the univariate Gamma distribution to construct bivariate families (Kotz et al., 2000).
In this paper, bivariate Gamma distribution of random vector Z = (Z g (x), Z g (x + h)) in R 2 is defined by its moment generating function or Laplace transform, which is defined with an affine polynomial (Bernardoff, 2006).We use an affine polynomial where the parameters satisfy the conditions: β > 0 and 1 > ρ > 0. Then the moment of the generating function of Z can be defined as ) From the defined generating function, it is obvious that Z g (x) is a distribution according to a univariate Gamma distribution with shape parameter α and scale parameter β, that is, Z g (x) ∼ G(x; α, β) (Zhang et al., 2010).Its probability density function can be expressed as follows where Γ(α) is the Gamma function.
The moment of the bivariate Gamma distribution can be obtained by differentiating Eq. ( 10).The mean and variance can be obtained as In isotropic stochastic processes, the covariance Cov (Z g (x), Z g (x+h)) of the bivariate Gamma function can be expressed as an exponential form 2 ( ), ( ) In this paper, an exponential family of covariance function is chosen for the bivariate Gamma model.The second-order variogram of Z g (x), denoted γ 2,g (h), is thus an exponential variogram defined by (Garrigues et al., 2007 where r g is the R of the variogram, it will be used to indicate the global texture structure of water in the real sea ice SAR images.αβ 2 is the sill of the variogram. Let X and Y be Gamma random variables with shape parameter α 1 and α 2 , scale parameter β 1 and β 2 , respectively.A random variable Z = X − Y has the following distribution function (Zhang, 1983).
For a stationary Gamma function Z g (x), the difference Z g (x) − Z g (x + h) is a random variable with mean and variance equal to zero and 2γ g (h).After the calculation of the distribution function of Z g (x) − Z g (x + h) with Eq. ( 14), the expectation of |Z g (x) − Z g (x + h)| can be figured out.Furthermore, from Eq. ( 7) γ 1,g (h) can be obtained (Sheldon, 2004).According to different values of α, γ 1,g (h) have different expressions.Taking α = 2 (corresponding to 2-look SAR imagery) as an example, γ 1,g (h) can be expressed as follows The combination of Poisson random lines and independent Gamma random variables within each cell defines our mosaic random function Z m (x).Lantuéjoul (2002) shows that the covariance function C m (h) of the mosaic model has an exponential form

Mosaic model Z m (x):
The second-order variogram of the mosaic model, denoted γ 2,m (h), is thus an exponential function where r m is the R of γ 2,m (h), it can be used to estimate the average size of sea ice floe in the real sea ice SAR images.For α = 2 the first-order variogram of mosaic model is defined as follows (more details will be discussed in Subsection 2.3.2)

Mixture model
To characterize different spatial scale we make a linear combination of Z m (x) and Z g (x), where σ 2 is the variance of Z c (x), ω is the weighting factor, it will be used to indicate the total concentration of sea ice.
Based on Eq. ( 17), the theoretical first-and second-order variogram of the mixture model will be calculated as follows.

Second-order variogram:
Since Z c (x) is a weighted sum of the independent random functions Z m (x) and Z g (x), its second-order variogram γ 2,c (h) is also the weighted sum of the second-order variograms γ 2,g (h) and γ 2,m (h) ( ) (1 ) ( )  13) and ( 17), we can write Eq. ( 20) as follow

First-order variogram:
As described above, the tessellation process of Poisson lines partitions the domain into non-overlapping cells within the image.The cell values are realizations of independent Gamma random variables.Two pixels x and x + h separated by a distance ||h|| belong to the same cell, denoted event A. Its probability is related to the covariance function of Z m (x) and thus to the second-order variogram of Z m (x) (Garrigues et al., 2007) The probability of A is equal to Conditioning on the event A and A , the theoretical first- order variogram can be decomposed as follow, where the first expectation value can be calculated from Eq. ( 19) In the case of event A, the locations x and x + h are in the same cell.As a result, Z m (x + h) -Z m (x) = 0, Eq. ( 25) reduce as follows, from Eq. ( 14) the expectation value of |Z g (x) -Z g (x + h)| can be figured out.By the same token, the second one will be In the case of the event A , the pixels x and x + h do not belong to the same cell.Z m (x + h), Z m (x), Z g (x) and Z g (x + h) are independent random variables with the same Gamma distribution.Using Eq. ( 14) the probability density function of the difference Z g (x + h) -Z g (x), as same as Z m (x + h) -Z m (x), can be figured out.ω and ( ) ] can be represented as a function of γ 2, g (h).Here also gives the expression of γ 1 (h) when α = 2 as an example.
( ) From Eq. ( 28), for the multi-Gamma and mosaic models, the first-and second-order variograms are related quadratically and linearly, respectively.

Parameter estimation by least-squares criterion
In this paper, least-squares criterion is employed to estimate the mixture model parameters, including first-and second-order variograms, by fitting theoretical variograms to experimental ones (Schabenberger & Gotway, 2005).Let H = [h 1 , ..., h n ] and Θ = (α, β, ω, r g , r m ).In real SAR images, α is a constant equal to the number of looks of SAR sensors.These parameters are used to characterize the spatial structure in a SAR intensity image.Theoretical variogram family γ(H, Θ) can be calculated from Eqs. ( 20) and ( 24). Let where N(h) is the number of the pairs of pixels separated by the distance h.For fitting γ(H, Θ) to where ε(H) is the n × 1 error vector with zero mean and variance-covariance matrix By weighted least squares, the parameter vector Θ can be estimated so as to minimize the weighted sum of squares

Comparison of theoretical variograms with different values of parameters
To test the feasibility of the mixture model, we compare the theoretical variograms with different values of parameters.
First is for the second-order variogram (Eq.( 21)).We set σ 2 = 1, ω 2 from 0 to 1 with interval 0.125, that is, ω 0 2 = 0, ω 1 2 = 0.125, ... In (a) it can be seen that the theoretical second-order variograms are close to each other all range of h, it is difficult to distinguish them.In (b) and (c), the theoretical second-order variograms for Z 5 to Z 7 can be distinguished when 3 pixels < h < 30 pixels, out of the range they are difficult to be discriminated.Z 1 to Z 4 are hard to be discriminated in all range of h.
We do the same comparison for the first-order variogram (Eq.( 28)).The given parameters' values are the same as the secondorder variogram comparison.Fig. 4 is the comparison results.The first-order variogram of Z 7 is very different from the other six ones, so it is not shown in the figure .From the figure we can find that the three results are similar.When 3 < h < 20 (pixel), the first-order variogram for Z 1 to Z 3 can be distinguished.Z 2 and Z 6 are mixed with each other and Z 3 to Z 5 are mixed with each other.After h = 20 pixels, Z 2 and Z 6 can be discriminated, but Z 3 to Z 5 are still difficult to discriminated.
Through the comparison above we can know that when r g ≠ r m and the weight of Poisson mosaic model is more than 50%, the mixture model with different values of ω 2 can be discriminated by the second-order variogram.When the weight of Poisson mosaic model is less than 50%, the mixture model with different values of ω 2 can be discriminated by the firstorder variogram.

Experiment on simulated images
To prove the correctness and accuracy of the first-and second-order variograms, the experiment is performed on simulated images first.
Two basic simulated images are generated from multi-Gamma model and Poisson mosaic model, respectively.Then combine them with different weight parameters ω i 2 , i = 1, ..., 7, to generate mixture simulated images.  2 = 0.875 to ω 1 2 = 0.125, respectively.Fig. 7 shows the results of the estimation of the first-and second-order variograms from simulated images.For ω 7 2 = 0.875 to ω 1 2 = 0.125, the theoretical variograms fits well the experimental variograms with accuracy over 90%.Because of the randomness of Gamma model, when h exceeds the R the values of variograms still have some fluctuations.When the weight of Gamma model is over 80%, the fluctuations become obvious.
Table 1 lists the results of estimated parameter ω 2 of mixture model Z(x).The values of ω 2 are the real weights of the two models, and the values of ω est 2 are the estimated weights.From the results, we can find that from Z 1 (x) to Z 3 (x) the estimation results from first-order variogram are more accurate than the second-order variogram's.But from Z 5 (x) to Z 7 (x) the secondorder variogram has better estimation results.) of the mixture model Z(x) estimated from each simulated images.

Experiment on real sea ice images
The proposed algorithm is tested by using SAR intensity images to identify sea ice structures.).Experiments will be performed for entire polygons, namely X, W, H, HH (Fig. 9), and squareshaped sub-regions of them (Fig. 10).
We choose some polygons which have more representatives for different types of sea ice to do the experiments.Fig. 9 is the Fig. 7 Result of the estimations of the first-order and second-order variograms from the simulated images.From (a) to (g) are corresponding to Z1(x) to Z7(x).From (a1) to (g1) are first-order variograms and from (a2) to (g2) are second-order variograms.
polygons of sea ice cut from Fig. 8. Their shapes are irregular and the spatial correlations of them are complicated.
According to the information given in the egg code we can know that even in one type of sea ice, such as X, there are two or three types of sea ice with different floe size.So in each polygon a square area is selected as a sample used for further experiments.Fig. 10 shows those samples with 100×100 pixels from corresponding polygons in Fig. 9.The samples have regular shape and in their areas the concentrations of sea ice are homogeneous.Figs.11 and 12 display the experimental first-and secondorder variograms computed over the polygon images and sample images along with the first-and second-order variograms computed with the estimated parameters of the mixture model.
Tables 2 and 3 list the information (floe size and total concentration) of sea ice derived from egg code and the estimated parameters of the mixture model.The ω est 2 of the mixture model is the proportion of the Poisson mosaic structure, which corresponds to the total concentration in the egg code.The R of the mosaic model r m is related to the mean size of the sea ice, which corresponds to the floe size in the egg code.The R of the multi-Gamma model r g is related to the global texture structure of the sea.Comparing the estimated results with the information derived from egg code we can find that the parameters in the mixture model can be used to indicate some sea ice parameters.Fig. 11 Result of the estimation of the first-order and second-order variograms from polygons (Fig. 9) Table 3. Parameter estimation results of polygons in Fig. 9

DISCUSSION AND CONCLUSION
This paper provides a new method to characterize the spatial structures of sea ice in SAR intensity imagery.The sea ice spatial structures are modeled as a weighted linear combination of two stochastic models: a Poisson line mosaic model and a multi-Gamma model.Then two kinds of geostatistic metrics are defined to describe the image spatial structures.By the experiment, we first compared the theoretical variograms with different values of parameters and found that each of the two different variograms have their own advantages in discriminating different spatial structures.Then the proposed method was applied to the simulated images and RARDARSAT-1 images.Through comparing the parameters estimated from the mixture model with the sea ice information obtaining from egg code, we found that two parameters (concentration and floe size) of sea ice can be retrieved from the mixture model.
Several notes can be remarked from this research.How to estimate other information in the egg code, such as stage of development, is another significant parameter of sea ice.The Voronoi mosaic model is another tessellation model, maybe its way of dividing is more suitable for describing sea ice structures.It will be used to replace Poisson mosaic model in the further study.In the real SAR images, the polygon images are irregular and the sea ice structures are not homogenous, the isotropy assumption is difficult to be satisfied.So the experimental variograms calculated from polygons are not as accurate as the types of sea ice r g (pixel) r m (pixel)  samples'.We should build new models under anisotropy assumption and research them.With the development of the data of remote sensing, the spatial resolution will be higher and higher.The Gamma distribution for modeling the intensity of SAR image is not applicable.Instead of Gamma distribution, some elaborate statistic distributions, for example: K distribution (Keith et al. 2006), are needed for the purpose.

Fig. 1
Fig.1 Diagram of C(h) and γ 2 (h).h represents the distance between two positions.C(h) is the covariance function and γ 2 (h) is the second-order variogram.
A mosaic model on a domain D can be defined by partitioning the domain D into a tessellation and assigning each cell of the tessellation a value independently drawn from a distribution.Poisson mosaic model is built based on the Poisson line tessellation, in which the image domain is partitioned by Poisson lines (see Fig. 2 (b)).Each line can be characterized by two parameters: its distance to the origin, denoted d with d > 0, and a random orientation θ with [0, 2π] (Kotz et al., 2000) (see Fig. 2 (a)).

Fig. 2
Fig.2 Plane linear expression (a), Poisson Division (b) of the two random variables, then the probability density functions of ωZ and variable Z c , Eq. (27) becomes σE[Z c ]. Using the probability density function of Z, the expectation of Z c can be calculated.Through β the E[Z c , where the first-and secondorder experimental variograms are computed according to , ω 8 2 = 1.Then we set three examples: (a) r g = r m = 30 pixels.(b) r g = 10 pixels, r m = 50 pixels.(c) r g = 50 pixels, r m = 10 pixels.When ω 0 2 = 0 and ω 8 2 = 1 the mixture model reduces to single mosaic model and single multi-Gamma model, respectively (not shown).Fig. 3 displays the comparison results.

Fig. 6
Fig.6shows simulated images, Fig.6(a) -(g) are corresponding 2 = 0.875 to ω 1 2 = 0.125, respectively.Fig.7shows the results of the estimation of the first-and second-order variograms from simulated images.For ω 7 2 = 0.875 to ω 1 2 = 0.125, the theoretical variograms fits well the experimental variograms with accuracy over 90%.Because of the randomness of Gamma model, when h exceeds the R the values of variograms still have some fluctuations.When the weight of Gamma model is over 80%, the fluctuations become obvious.Table1lists the results of estimated parameter ω 2 of mixture model Z(x).The values of ω 2 are the real weights of the two models, and the values of ω est 2 are the estimated weights.From the results, we can find that from Z 1 (x) to Z 3 (x) the estimation

Fig. 12
Fig.12Result of the estimation of the first-order and second-order variograms from samples (Fig.10)

Table 2 .
Parameter estimation results of polygons in Fig.9