Incident Angle Dependence of Sentinel-1 Texture Features for Sea Ice Classification

Robust and reliable classification of sea ice types in synthetic aperture radar (SAR) images is needed for various operational and environmental applications. Previous studies have investigated the class-dependent decrease in SAR backscatter intensity with incident angle (IA); others have shown the potential of textural information to improve automated image classification. In this work, we investigate the inclusion of Sentinel-1 (S1) texture features into a Bayesian classifier that accounts for linear per-class variation of its features with IA. We use the S1 extra-wide swath (EW) product in ground-range detected format at medium resolution (GRDM), and we compute seven grey level co-occurrence matrix (GLCM) texture features from the HH and the HV backscatter intensity in the linear and logarithmic domain. While GLCM texture features obtained in the linear domain vary significantly with IA, the features computed from the logarithmic intensity do not depend on IA or reveal only a weak, approximately linear dependency. They can therefore be directly included in the IA-sensitive classifier that assumes a linear variation. The different number of looks in the first sub-swath (EW1) of the product causes a distinct offset in texture at the sub-swath boundary between EW1 and the second sub-swath (EW2). This offset must be considered when using texture in classification; we demonstrate a manual correction for the example of GLCM contrast. Based on the Jeffries–Matusita distance between class histograms, we perform a separability analysis for 57 different GLCM parameter settings. We select a suitable combination of features for the ice classes in our data set and classify several test images using a combination of intensity and texture features. We compare the results to a classifier using only intensity. Particular improvements are achieved for the generalized separation of ice and water, as well as the classification of young ice and multi-year ice.


Introduction
Synthetic aperture radar (SAR) is a primary tool for monitoring of sea ice conditions in the polar regions [1][2][3]. A radar system is an active device that both transmits and receives electromagnetic radiation in the microwave region, and is thus independent of sunlight and cloud conditions. The resulting continuous imaging capability of the SAR is important for operational ice services worldwide [3,4]. The analysis and interpretation of the SAR images and the production of ice charts is at present carried out manually and therefore subject to the expertise of the individual ice analyst [5,6]. Furthermore, while timeliness of ice charts is a critical requirement, the manual image analysis is a time-consuming process [7]. In combination with an increasing volume of available SAR imagery, this underlines the need for automated or computer-assisted classification of sea ice. The backscatter signature of sea ice in radar images, however, depends on a variety of different factors, including sea ice, environmental, and radar parameters [8,9]. Despite multiple efforts and various approaches, robust and automated classification of ice types therefore remains a challenging task [3].
The important radar parameters that influence the signal include radar frequency, polarization, and incident angle (IA) [9]. While frequency and polarization are fixed for a given sensor and operation mode, IA varies across the image. The backscatter intensity (also referred to as image brightness or tone) from a homogeneous surface varies with IA and decreases across a SAR image from near-range (low IA) to far-range (high IA). In the logarithmic domain, i.e., with the intensity given in dB, the decrease is approximately linear with a constant slope per surface class [10][11][12][13]. Despite many studies showing that the rate of decrease differs between different surface types, most classification approaches apply a global IA correction during pre-processing, using a constant slope for the entire SAR image [14][15][16][17][18]. Although such approaches can achieve good results, they neglect the known physical differences in backscatter behavior with IA for different surface types. Lohse et al. [19] recently introduced a method to directly incorporate the class-dependent effect of IA into a supervised algorithm. The method is based on a Bayesian classifier (i.e., it maximizes probabilities) with multi-variate Gaussian probability density functions (PDFs), where the constant mean value is replaced with a linearly varying mean value. The linear slopes are class-dependent and thus directly included in the classifier. Since it is based on underlying Gaussian PDFs and accounts for a per-class IA effect, the classification method is referred to as the GIA (Gaussian incident angle) classifier. The approach achieves improved classification results compared to intensity-based methods with a global IA correction. However, ambiguities in backscatter intensity remain for individual classes at some IA ranges [19]. In particular, changes in the sea surface roughness (the sea surface state), caused for example by varying wind conditions, ocean currents, or natural and oil slicks, can complicate the reliable classification of open water (OW).
Previous studies have shown that in many cases textural information can help to resolve ambiguities in sea ice classification, both for the binary problem of ice-water classification and for the multi-class separation of different ice types [18,[20][21][22][23][24][25][26][27]. Texture generally refers to the local spatial variation of tone or brightness within an image at a given scale [28]. While many different ways of extracting texture features exist [24], the most common texture features used for sea ice classification are based on the grey level co-occurrence matrix (GLCM) [28]. A straightforward way to directly utilize information from the GLCM in a pixel-based classifier is to extract scalar features from the matrix. Such GLCM-based texture features have been used in a variety of studies and algorithms and improved overall classification accuracy (CA) of OW areas vs. sea ice [14,18,[21][22][23][24][25]27,29]. Texture extraction, and in particular computation of the GLCM, requires several input parameters such as window size, quantization levels, displacement distance, and displacement direction (see Section 3 for details). Many of these parameters have been investigated in previous studies. The optimal choice, however, differs between studies (Table 1) and depends on class definitions and data properties. To our best knowledge, a systematic investigation of the dependence of common GLCM texture features on IA for different classes has not been performed prior to this study. All approaches presented in the literature that use the GLCM for the analysis of sea ice imagery either apply a global IA correction or no correction.
In this study, we therefore investigate the per-class IA dependence of different texture features. We do so with the intention of incorporating texture directly into the GIA classifier, which accounts for per-class IA effects. The main prerequisite for this incorporation is that there must be a clearly defined relationship between texture parameters and IA. Ideally, for the linear GIA, this relationship should be constant or linear, and the distribution of the individual features around the linear function can be approximated as Gaussian. These ideal conditions allow for the direct use of the texture features in the existing algorithm, while a more complicated IA relationship or a clearly non-Gaussian distribution would require changing the underlying model of the GIA classifier. After we confirm that common texture features fulfill these conditions, we select a useful set of features and demonstrate the benefits of including them in the classification process.
This paper is organized as follows. Section 2 gives an overview of the data set, including the standard pre-processing steps that were applied. The outline of the investigations in this study is presented in Section 3, followed by a detailed description of the computation and selection of the texture features as well as the tested parameter settings. Section 4 presents the results of these investigations. We discuss our findings in Section 5, pointing out implications and limitations for the usage of texture in sea ice type classification from SAR data, and in particular in the GIA classifier. In Section 6, we summarize the main findings and conclusions. Table 1. Overview of GLCM computation parameters from selected studies that investigate GLCM texture features for sea ice classification. The parameters are window size w, co-occurrence distance d, displacement angle α, and number of grey level quantization intervals k. w and d are given in number of pixels, α is given in degrees. Bold type indicates the preferred choice selected by these studies, where applicable. The dB column indicates whether texture was computed from the linear or from the logarithmic intensity. The question mark indicates that this information is not explicitly given in the study. All SAR imagery in this study is Sentinel-1 (S1) data acquired in extra-wide swath (EW) operation mode [30]. S1 operates at C-band (5.4 GHz) providing either single-or dualpolarization data. As part of the European Copernicus Earth observation program, all S1 data are freely available (e.g., through the Copernicus Open Access Hub). The data used in this study are acquired at dual-polarization (HH and HV) and downloaded in groundrange detected format at medium resolution (GRDM). The EW GRDM product comes at a pixel spacing of 40 × 40 m with an actual spatial resolution of approximately 93 × 87 m; its values are multi-looked intensities with 18 looks in the first sub-swath EW1 and 12 looks in the remaining sub-swaths EW2 to EW5. As standard pre-processing, we apply the thermal noise correction implemented in ESA's Sentinel Application Platform (SNAP) and calibrate the data to obtain the normalized radar cross-section σ 0 . All processing is performed in the ground-range detected image geometry.

Training and Validation Data
We use the training and validation data set introduced by Lohse et al. [19] in 2020. This data set is based on the visual inspection and expert analysis of overlapping SAR (S1) and optical remote sensing data acquired during winter conditions between 2015 and 2019. The main identified classes in the data set are open water (OW), leads with OW or newly formed ice (NFI), brash or pancake ice, young ice (YI), level first-year ice (LFYI), deformed first-year ice (DFYI), and multi-year ice (MYI). A detailed description of the data set, image locations and acquisition times, and the selection of classes and training polygons is given in Lohse et al. [19]. For parts of this study, we have added new images and test polygons to the existing data set. Generally, the training data for individual classes used in this study are collected from a large number of scenes. In some cases we show training data from a single image for a particular class. The image IDs (product unique IDs) are given in the respective subsections of this article.

Domain-Dependent Texture Extraction and Calculation of Separability
After pre-processing, we extract various texture features from both the HH and the HV channel of the S1 images. Initially, we compute all texture features from both the backscatter intensity in the linear domain (normalized radar cross-section σ 0 ) and the backscatter intensity in the logarithmic domain. In the logarithmic domain, the intensities are given in decibels (dB): HH dB = 10 · log 10 (σ 0 HH ) HV dB = 10 · log 10 (σ 0 HV ) The relationship of intensity with IA is approximately exponential in the linear domain, and thus in turn approximately linear in the logarithmic domain ( Figure 1). These differences in IA dependence, in combination with the change of variance with IA in the linear domain, are expected to translate into differences in IA dependence of the texture features extracted from the respective intensity domains. The initial extraction of texture from both linear and logarithmic intensity allows us to find the preferred domain to compute the texture features. For the preferred domain, we test a variety of parameter settings (see Sections 3.2.1 and 3.2.2 for details), and investigate the variation of the extracted features with IA for these different settings. We adjust the borders of the training regions according to the size of the texture windows, such that little or no mixing of ice classes occurs within the texture windows. We then use the Jeffries-Matusita (JM) distance [31] to evaluate all features and parameter settings in terms of class separability for various two-class cases.
The JM distance is an established separability measure between class distributions that returns values between zero (no separability) and two (perfect separability). It is given by where D B is the Bhattacharrya distance according to Features with a JM value above one are commonly considered useful for classification [31].
On the basis of this evaluation, we finally select a suitable feature set and demonstrate the benefits of incorporating textural information into the GIA classifier. In particular, we present examples for the classification of sea ice against OW, as well as YI against MYI, and compare against results obtained from a classifier based on intensities only.

Calculation of Texture Features
Since many previous studies suggest that the GLCM provides a useful method to generate texture features that can improve sea ice classification results [24,26], we focus on this method (Section 3.2.1). However, calculation of the GLCM is computationally expensive and time-consuming, which can impede its use in operational applications. As an example for a simpler and more easily calculated texture feature, we therefore also extract the variance (Section 3.2.2) and assess its usability compared to the GLCM features.

GLCM Texture Features
The GLCM provides a second-order statistic for extraction of texture features [25,28]. It calculates the probability of a pixel with grey level value i occurring at a certain distance and angle from another pixel with grey level j within a given window. The key parameters that must be set are the window size w for which to calculate the GLCM, the co-occurrence distance d, the displacement angle α, and the number of grey level quantization intervals k. Algebraically, the GLCM can be expressed as where S w,d,α,k is an element of the GLCM for a given window size, direction, co-occurrence distance, and grey level quantization; P w,d,α,k is the frequency of occurrence of grey levels i and j; and k is the number of quantized grey levels. The size of the GLCM depends on the number of grey levels. In order to neglect effects from ice floe rotation and changes in the angle between the radar-look direction and the physical structures on the ice, the GLCM is often calculated for different directions (0 • , 45 • , 90 • , 135 • ) and then averaged before feature extraction [14,25]: The resulting averaged GLCM S still includes the effects of directional structures, but the effects are diluted by the averaging. The specific orientation of the structures is not reflected any more in the averaged GLCM. The remaining parameters are usually chosen manually or optimized for a particular study and then kept fixed. Individual scalar texture features can be calculated from the GLCM according to Equations (6)- (12): Correlation : Energy : Ene = √ ASM (10) Different tested parameter settings and feature choices from selected studies that use GLCM texture features for the interpretation of sea ice SAR imagery are summarized in Table 1. If provided, the preferred/optimal choice of each study is indicated with bold type. It is evident that there is no consensus in the literature on an optimal set of GLCM features and parameters. The choices that lead to the best classification results for the individual studies differ from window sizes between 5 and 64 pixels, co-occurrence distances between 1 and 8, and grey level quantization levels between 8 and 64. The optimal features and parameter set depend on the class definitions, the image pre-processing steps (in particular multi-looking and re-sampling), and the data properties (in particular frequency, spatial resolution, polarization). In this study, we therefore test a variety of GLCM parameter settings (Table 2), covering a reasonable range of settings that is based on the literature values in Table 1. Our goal is to assess the effect of different settings on potential IA dependence of the features, and to find a suitable set of features and parameters for the specific data set that we use.
To be consistent and to ensure identical quantization for all images, we choose a uniform quantization with equally spaced grey level intervals. We clip the minimum and maximum grey levels at −35 and +5 dB for HH and −40 and 0 dB for HV, respectively. To avoid directional effects, we average the GLCMs obtained for four directions (0 • , 45 • , 90 • , 135 • ) before computation of individual scalar features. Table 2. Summary of GLCM computation parameters tested in this study (w: window size, d: GLCM co-occurrence distance, k: GLCM grey levels). The window is applied on the original 40 × 40 m pixel spacing of the S1 EW GRDM product. Quantization is performed with equally spaced intervals between −35 and +5 dB for HH and −40 and 0 dB for HV, respectively. The GLCM is calculated for four different directions (0 • , 45 • , 90 • , 135 • ) and then averaged before feature extraction. Computation of the GLCM is time-consuming and depends on multiple different parameters. Hence, it is interesting to also test other texture features that can be calculated faster and more easily, and investigate if they can be used instead of GLCM features. In this study, we investigate the variance (Var) as an example for such simpler texture features. The only required input parameter is the window size. We calculate variance from the logarithmic intensity for the same window sizes as the GLCM features (Table 2). Unlike the GLCM, variance does not depend on a distance parameter inside the defined texture window. Additionally, it is not sensitive to the spatial orientation of physical structures on the ice. However, many of the physical structures that we are interested in, for example leads or pressure ridges, have some specific spatial orientation. Even though we are not interested in this specific orientation itself, looking in different directions can be necessary to detect the physical structures' presence. Hence, the larger computational effort of the GLCM can be beneficial to detect certain structures that the variance may miss, depending on the applied window size and the physical size of the structures on the ground. A comparison between the different approaches to quantify texture in terms of class separability is therefore useful.

Results
In this section, we present the results of our study. As we have tested a large number of different texture parameter settings and individual features, it is not feasible to show a full overview of all tested combinations. We therefore present and discuss representative examples for the various experiments performed in this study.

Texture and IA
We begin by comparing the influence of IA on GLCM texture features computed from intensity in dB against GLCM texture features computed from linear intensity. An example including both intensity domains and three selected GLCM features is shown in Figure 2. When computed from the logarithmic intensity, the GLCM features showed no significant per-class variation with IA ( Figure 2, upper panel); when computed from the linear intensity, there was an evident variation ( Figure 2, lower panel). For some features (for example entropy and homogeneity) the relationship appeared to be approximately linear over part of the shown IA range. Overall, however, the IA dependence of the GLCM texture features was significantly more complicated when computed from linear intensity compared to logarithmic intensity. Furthermore, the OW class in the example of Figure 2 showed some internal variability in intensity, caused by changes in the sea surface state across the image. The GLCM texture from the logarithmic intensity appeared to be less sensitive to such internal class variation than the GLCM texture from the linear intensity. Density distribution of HH intensity and three selected GLCM texture features with IA for OW training data selected over the full range of a single S1 image (Image ID: 77 BA). The upper panel shows HH intensity in the logarithmic domain (in dB) and the GLCM features extracted from it; the lower panel shows intensity in the linear domain with its respective GLCM features. Note that the OW class displays some internal class variation in intensity due to varying sea surface state across the image. GLCM parameter settings: w = 11, d = 4, k = 16.
All tested texture features revealed a significant offset at the boundary between the first and the second sub-swath of the image, which is located at an IA of approximately 28.5 • . This offset was observed for all tested parameter settings, and it occured independently of the intensity metric. Since the offset will affect the performance of the texture features in any classifier, it requires further investigation. It is reasonable to assume that the offset in texture values at the sub-swath boundary is at least partly caused by the different number of looks in the respective sub-swaths of the S1 EW GRDM product. The image contrast, for example, is directly linked with the variance of the image brightness within a given window. A change in the number of looks causes a corresponding change in the magnitude of variance. Therefore, we can presumably correct the contrast values in sub-swath EW1 manually by multiplying with the square root of the ratio of number of looks between EW1 and EW2. Figure 3 shows the distribution of HH contrast with IA before and after the correction. The sub-swath boundaries are visible in the distribution and additionally marked in the figure. The proposed correction successfully removed the offset. Note that the correction may be less straightforward for other texture features, depending on the formula for their computation.

EW1
EW2 EW3 EW4 Figure 3. Density distribution of GLCM contrast (computed from HH intensity in dB) with IA for OW training data selected over the full range of a single S1 image (Image ID: F2FE). The left side shows contrast computed directly from the GRDM product; the right side shows contrast manually corrected by multiplying with the square root of the ratio of number of looks in the sub-swaths. The sub-swath boundaries (indicated with red arrows) are clearly visible in the texture profile. In particular the boundary between EW1 and EW2 is characterized by a distinct offset that requires a correction. GLCM parameter settings: w = 21, d = 4, k = 16.
On the basis of findings from the comparison of texture from linear intensity against texture from logarithmic intensity, all of the following calculations were performed with intensity given in dB. Further tests with different GLCM parameter settings (57 different settings in total, Table 2) confirmed that the selected GLCM features were not (or only weakly) dependent on the IA. Figure 4 shows examples for two selected features and four representative parameter settings for the MYI training data. The changes in the different parameters (w, d, and k) clearly affected the numerical values of the texture features. The variance of the distributions around the mean value decreased with increasing window size w (Figure 4 from left to right). Therefore, the linear trend in the feature distribution with IA was more easily visible for larger window sizes (>∼21 pixels). The dependency of texture with IA was linear (and almost constant) for all the tested parameter settings; hence, the parameter settings did not affect the general inclusion of GLCM texture into the concept of the GIA classifier.
All texture features shown so far have been extracted from the HH channel of the data. Figure 5 shows an example of three selected texture features extracted from both HH and HV channel for LFYI training data. For a weak signal close to or below the nominal noise floor of the sensor, the distribution of texture with IA was strongly affected by the noise profile. This problem is of particular importance in the HV channel, as the signal at HV polarization was generally weaker than the signal at HH polarization. These noise floor artifacts in the texture features will complicate the inclusion of the HV texture features in the GIA classifier. Additionally, texture profiles for both channels displayed artifacts at the sub-swath boundaries between EW2 and EW5. The stronger these artifacts are, the more they will affect the use of texture across sub-swath boundaries in the classification. Note that we have only applied the calibration and nominal noise-floor corrections provided by ESA in the SNAP software. Improved noise-floor corrections and more accurate data calibration between the sub-swaths may remedy this problem in the future.  Our tests for the variance as an example of a more simply calculated texture feature give similar results (not shown). When computed from the HH channel in dB, variance was approximately constant over the full range of the image, except for a significant offset between sub-swaths EW1 and EW2. When extracted from the HV channel, which often has a signal strength close to the nominal noise floor, the noise profile was clearly visible in the IA relationship of the variance. The assumptions needed for the GIA classifier (linear IA dependence, approximately Gaussian distribution) are then violated.

Texture and Different Sea Surface States
One of the main challenges in sea ice classification is to automatically separate sea ice and OW. As spatial and temporal variations in sea surface state affect the backscatter intensity from OW, a purely intensity-based classifier will often struggle with a generalized separation of sea ice and water, unless additional information is included. We now confirm that texture can help to overcome the issue of different backscatter intensity for varying sea surface states. Figure 6 shows training data of OW areas collected from multiple S1 images. The left column shows the density distributions of HH intensity in dB with IA. Clearly, the intensity levels differed between the images; assuming that the intensity level can be explained by Bragg scattering, its variations can be explained by differences in OW surface roughness. The numerical values of the texture features were consistent over the selected images and wind states ( Figure 6, column 2, 3, 4). In agreement with previous results (Figure 4), the numerical texture values and the width of the distributions differed between different GLCM parameter settings (not shown here). Furthermore, the distributions of the texture features remained constant over the IA range, except for the offset between EW1 and EW2. Figure 6. Density distribution of HH intensity in dB and three selected GLCM texture features with IA for OW training data selected from three different images (image IDs: C113, 6346, D4FB). The lowest panel shows the data from all three images combined. The intensity values clearly differ between the images, and the texture values are consistent. The offset in texture between sub-swaths EW1 and EW2 is caused by the different number of looks in the respective sub-swaths. The texture artifacts at the remaining sub-swath boundaries are caused by errors in the calibration and noise correction between the sub-swaths. GLCM parameter settings: w = 21, d = 4, k = 16.

Separability of Different Classes
The results presented so far show that the tested GLCM texture features can be directly incorporated into the concept of the GIA classifier, given that the signal is strong enough to avoid noise floor artifacts. The GIA classifier assumes a linear relationship of its features with IA, and an approximately Gaussian feature distribution; thus, the incorporation of GLCM texture features calculated from intensity in dB is straightforward. For the S1 EW GRDM product, the first sub-swath should be ignored or corrected according to the number of looks. However, for the inclusion of the texture features to be useful in terms of improving classification results, we need to investigate the separability of different classes for individual features and varying parameter settings. Note that at this point we are interested in the individual separability of different pairs of classes, as shown on the left side of Table 3. The results of this kind of analysis allow a better understanding of the specific gain achieved by single texture parameters, and they are useful, for example, for designing decision trees such as described in [32].  We performed this analysis by computing the JM distance between the different feature distributions for all settings. Since the HV texture is strongly affected by the noise, we only evaluate texture features extracted from the HH channel at this point. In total, we have analyzed 57 different parameter settings (Table 2) for eight separate features (seven GLCM features (Equations (6)- (12)) and variance). Table 3 presents the JM distance for eight selected parameter settings. Corresponding class distributions for four of these settings are shown in Figure 7. For the given data set and the parameter settings tested in this study, we found that window size was the parameter with the largest influence on class separability. If a feature offers any separability at all between two tested classes, the separability improves with increasing window size (Table 3 and Figure 7, rows 1, 2, and 5). In agreement with previous studies, we found that many of the tested GLCM features allowed for partial separation of OW from thicker ice types (that is LFYI, DFYI, and MYI). In particular, HH ASM, HH contrast, HH dissimilarity, HH energy, HH entropy, and HH homogeneity all displayed distinct class distributions that revealed at least some separability of LFYI/DFYI/MYI and OW, with JM distance values close to one (for w = 21) and significantly above one (for w = 51). Figure 8 presents more histogram examples of HH dissimilarity and HH energy for the classes OW and MYI. Again, significant improvement in separability with increasing window size was clearly visible. Separation between the aforementioned thicker ice types (LFYI, DFYI, MYI) is not possible based on the investigated texture features only. This can be seen by the overlapping histograms for LFYI and MYI in Figure 7 (row 4) and was confirmed by the low JM distances between the distributions for all parameter settings and features (Table 3). Partial separation of LFYI and MYI may be possible using the GLCM correlation from the HV channel (Figure 9, right-hand side); however, the HV signal is often close to the nominal noise floor and the channel must be treated cautiously. None of the tested features separated well between YI and OW. The distributions of these two classes overlapped significantly for all features and all parameter settings (Figure 7, row 3), resulting in low JM distances. However, HH texture showed the potential to improve the classification of MYI against YI in refrozen leads. For several features, JM distances between these two classes were close to and exceeded one at large window sizes (w = 21, 51). Histogram examples of HH dissimilarity and HH energy for YI and MYI are shown on the left-hand side of Figure 9. Histograms of HH dissimilarity, HH energy, and HV correlation for YI against MYI (left) and LFYI against MYI (right). YI and MYI can be partly separated using HH texture. LFYI and MYI are inseparable in HH texture, but can be partly separated in HV correlation, although significant overlap of the distributions remains.

Classification Result Examples
Finally, we demonstrate the potential of including useful texture features in the GIA classifier to improve classification results. On the basis of class separability indicated by the JM distances, we chose different feature combinations and trained the GIA algorithm to classify various test images. Three examples are shown in Figures 10 and 11. The features used for the classification of the presented images were HH intensity, HV intensity, HH contrast, HH dissimilarity, and HH energy. Because of the offset in numerical texture values between sub-swaths EW1 and EW2, we masked out the results for the first sub-swath. Figure 10 shows all used features (a, b, c, and d) together with the classification result (e) and a map of sea ice concentration (SIC) (f), which can be directly obtained from the classification result. OW and sea ice were well separated in the image, and the ice edge can be successfully detected in both the ice type map as well as the SIC. Some areas within the pack ice were classified as OW; without complementary information, it is difficult to assess whether these areas are in fact OW or YI. Classification of an image containing large OW areas (such as the example in Figure 10) is challenging based on intensity only, as one would need to know the sea surface state to select the correct intensity level and training data for the OW class. Without that additional information, it is only the inclusion of textural information that makes the classification of the image with a generalized classifier feasible. We therefore do not present a result based on intensity only for this case.  Figure 11 shows false-color images (a and b) and classification results (c, d, e, and f) for two different images that were almost entirely covered by sea ice. The classification results in the middle column (Figure 11c,d) were obtained with a GIA classifier based on intensities only. The same generalized classifier was used for both images. While the classifier captured the YI areas in the lower image correctly (ID B7A9), there was significant mis-classification of YI areas as MYI in the upper image (ID 89A6). This mis-classification occurred despite the fact the YI areas appeared visually similar in both images. A minor change in the properties of the YI areas and thus in the backscatter intensity from the surface can result in the confusion of YI and MYI. The texture signatures of the two classes can help to solve this problem. When the classification was performed including the textural information (Figure 11e,f), the YI areas were classified correctly in both images. The examples demonstrate the superior capability of the GIA classifier using both intensity and texture to separate YI and MYI in situations where the purely intensity-based classifier fails.

Discussion
We have investigated the class-dependent variation of different texture features with IA, in order to include them into the GIA classifier. The GIA classifier accounts for perclass variation of its features with IA, and it requires approximately linear relationships and Gaussian distributions. Our results show that it is possible to directly incorporate GLCM texture into the GIA classifier. When computed from intensity in the logarithmic domain (that is in dB), the tested GLCM features do not depend on IA or reveal only a weak, approximately linear dependency. When computed from intensity in the linear domain, the tested GLCM features show considerable variation with IA. For some features the variation is approximately linear over part of the IA range, while for other features it is more complicated. Generally, we therefore recommend to compute texture from the logarithmic intensity; the slope of texture with IA is then constant (approximately zero), and the features can be used in the GIA classifier that assumes a linear relationship. Furthermore, for any other application in a different classification algorithm, as long as texture features are computed from logarithmic intensity, no global correction of the IA effect during pre-processing is needed. We have tested a variety of different GLCM parameter settings and find that the results regarding the IA dependence of the features are independent of the parameter choice. They hold for all tested window sizes, grey levels, and distances. For classes with a weak intensity signal, the S1 noise profile is clearly visible in the distribution of texture features as a function of IA. While this noise pattern will cause problems for any classification algorithm, it is particularly challenging for the GIA classifier, as the basic assumptions of linear IA dependence and Gaussian distributions are violated. The HV channel is more problematic than the HH channel in this regard because the signal at HV polarization is often weak and close to the nominal noise floor. HV texture should thus be used carefully, and it will benefit from an improved thermal noise correction of the S1 data.
Furthermore, we find that there is a distinct offset in all texture features at the subswath boundary between EW1 and EW2 of the S1 EW GRDM product. This offset is caused by the different number of looks applied to the individual sub-swaths; the GRDM product is delivered with 18 looks in EW1 and 12 looks in EW2 to EW5. We demonstrate how to correct this offset for the example of GLCM contrast, by multiplying with the square root of the ratio of number of looks in the different EW GRDM sub-swaths. Corrections can also be applied for other GLCM features, but it will depend on the formula for the texture computation. The exact correction factors for all individual features require further investigation and are not part of this study. However, the offset must be considered whenever using texture features extracted from S1 wide-swath products in GRDM format across the full range of the image. Generally, it must be kept in mind that texture is dependent on the speckle contribution in the intensity, and thus on the number of looks in the data.
When integrated into the GIA classifier, GLCM texture features help to resolve some of the inherent ambiguities found with intensity-only classifiers. For example, the general separation of OW and thicker sea ice types, such as LFYI, DFYI, and MYI, is significantly improved by the inclusion of texture. While the backscatter intensity from OW is dependent on the sea surface state and may require the training of several OW classes for different sea surface conditions, the texture signature of OW is nearly independent of the sea state. It suffices to train one OW class for a smooth water surface (where the signal will be close to the nominal noise floor), and one OW class for all rough surface conditions ( Figure 6). Another challenge of a classifier based on intensity only is the separation of YI and MYI. Especially for YI with frost flowers or a snow crust, which increases the small-scale roughness and causes strong backscatter from the YI surface, mis-classification of YI as MYI can occur [14,33]. We demonstrate in this study that the inclusion of texture features can significantly improve the separation of YI and MYI.
While some ambiguities of the sea ice type classification can be improved or solved by adding texture features, other classes are not separable in the texture feature space alone. Their per-class texture distributions overlap significantly for all tested features and parameter settings. Hence, backscatter intensity remains an important feature for the classification of these ice types. In particular, this is true for the separation between the thicker ice types (LFYI, DFYI, and MYI), and for the separation of YI and OW. FYI and MYI can be distinguished quite well based on their intensity, as the less saline MYI will cause more volume scattering [34], which results in a stronger signal at both HH and HV polarization. We therefore recommend to always include intensity as a feature in ice type classification. YI and OW, on the other hand, can be more challenging. When the sea surface state, and thus the intensity level of OW, is known, backscatter intensity can be used to overcome this ambiguity. However, since the sea surface conditions are not usually known a priori, the separation of YI and OW remains difficult at this point. We will explore possible solutions such as the inclusion of SIC from passive microwave observations [35] or the results of SAR wind retrieval algorithms [36] in our future work.
The evaluation of features and parameter settings in this study is based on JM distances between different class distributions. Generally, we find that larger window sizes improve the separability between the classes in the used data set. It must be kept in mind, however, that large texture windows result in smoothing and thus in a coarser effective resolution of the results. Hence, there is a trade-off between spatial resolution of the results and separability of sea ice and OW, as well as YI and MYI.
The main focus of our analysis has been on the use of GLCM texture features, as the GLCM method has been widely used for sea ice image analysis [18,[20][21][22][23][24][25][26][27]. However, the calculation of the GLCM requires substantial computational resources and is time consuming. For operational ice charting, where timeliness of the results is a main requirement, this can be problematic. We have therefore also tested variance as an example of a more simply calculated texture feature. We find that variance from intensity in dB is roughly independent of IA and can potentially be used as a faster and alternative to GLCM texture. Further investigation of other computationally cheap texture features, such as for example the coefficient of variation or wavelets, is needed in future work and may help to speed up the process of automated sea ice classification, making application of the algorithms feasible in operational ice charting.

Conclusions
In this study, we have investigated the IA dependence of seven commonly used GLCM texture features extracted from the S1 EW GRDM product, and we assessed their potential to be included in IA-sensitive sea ice classification (the GIA classifier). When calculated from intensity in dB, the GLCM features are found to be almost independent of IA and can thus be directly included in the GIA classifier, with the estimated slope being approximately zero. Particular attention must be paid to classes with a weak signal, which will lead to noise artifacts in the texture parameters, and to the first sub-swath of the EW GRDM product, as the different number of looks in the sub-swaths results in an offset of texture at the sub-swath boundary. We have shown how this offset can be corrected for the example of GLCM contrast.
We have tested a large number of GLCM parameters and evaluated the resulting features in terms of class separability. Using per-class histograms of feature distributions in combination with the JM distance, we have selected meaningful combinations of texture features and backscatter intensities to train a classifier and demonstrate the improvement in ice-water classification as well as the separation of YI and MYI on various examples compared to classification based on intensity only. Our analysis shows that larger texture windows (up to 51 × 51 pixels) generally result in better class separability, albeit at the cost of reduced spatial resolution of the image.  Acknowledgments: The presented work contains primary and altered Sentinel data products (@Copernicus data). We would like to thank the anonymous reviewers for their constructive comments and feedback.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: