Machine Learning-Aided Sea Ice Monitoring Using Feature Sequences Extracted from Spaceborne GNSS-Reflectometry Data

Two effective machine learning-aided sea ice monitoring methods are investigated using 42 months of spaceborne Global Navigation Satellite System-Reflectometry (GNSS-R) data collected by the TechDemoSat-1 (TDS-1). The two-dimensional delay waveforms with different Doppler spread characteristics are applied to extract six features, which are combined to monitor sea ice using the decision tree (DT) and random forest (RF) algorithms. Firstly, the feature sequences are used as input variables and sea ice concentration (SIC) data from the Advanced Microwave Space Radiometer-2 (AMSR-2) are applied as targeted output to train the sea ice monitoring model. Hereafter, the performance of the proposed method is evaluated through comparing with the sea ice edge (SIE) data from the Special Sensor Microwave Imager Sounder (SSMIS) data. The DTand RF-based methods achieve an overall accuracy of 97.51% and 98.03%, respectively, in the Arctic region and 95.46% and 95.96%, respectively, in the Antarctic region. The DTand RF-based methods achieve similar accuracies, while the Kappa coefficient of RF-based approach is slightly larger than that of the DT-based approach, which indicates that the RF-based method outperforms the DT-based method. The results show the potential of monitoring sea ice using machine learning-aided GNSS-R approaches.


Introduction
Sea ice monitoring shows significant importance because it has notable impacts on the Earth's radiation balance, which affects the global climate significantly. Therefore, having a good knowledge of sea ice extent and distribution is critical for the study of climate change [1].
Sea ice has been monitored with various approaches, such as field observations [2], numerical models [3] and remote sensing [4], the latter of which has been considered as the most efficient which depict the characteristics of DDMs, as input elements may enhance monitoring performance and data processing efficiency. The applications of ML may be categorized into three aspects [35]: classification, developing empirical model and improving computation efficiency. One of the most important parts of sea ice monitoring is to distinguish sea ice from water, which can be regarded as a classification problem. As an ML method, the decision tree (DT) method has been widely applied to sea ice monitoring [41,42]. Another powerful ML algorithm employed for classification is random forest (RF), which creates a variety of individual decision trees that operate as an ensemble [43]. Although the DT and RF algorithms have been applied to monitor sea ice using satellite remote sensing data, such as MODIS and CryoSat-2, there is a lack of information about how DT and RF algorithms can be utilized for monitoring sea ice using spaceborne GNSS-R data. The task of this study is to explore the potential application of spaceborne GNSS-R to distinguish sea ice from water using the DT and RF classifiers. Section 2 firstly gives the description of datasets used in this study and the extraction of features. Then, the sea ice monitoring approaches based on DT and RF algorithms and data processing flow are presented in Section 2. The sea ice monitoring results are presented and discussed in Sections 3 and 4, respectively. Finally, the conclusions are addressed in Section 5.

TDS-1 Mission and Datasets
Spaceborne GNSS-R data from TDS-1 include three different data processing levels, e.g., Level 0 (L0), Level 1 (L1) and Level 2 (L2) [44]. L0 mainly contains the raw data, which are not available to the public except for some sample data. L1 includes L1a and L1b, which are the data converted from the L1a onboard processed DDMs and converted to NetCDF format. The L1b release includes the DDMs and metadata used in this study. Level 2 refers to the wind speed and mean square slope products. DDMs are generated by the Space GNSS Receiver Remote Sensing Instrument (SGR-ReSI) through cross-correlating scattered signals with code replicas generated locally for different time delays and Doppler shifts. When the reflection surface is smooth, most of the scattered power comes from the specular point, and very little from the glistening zone around the specular point [45]. Compared with the sea ice surface, the one of sea water yields a non-coherent reflection with more scattering in the delay and Doppler domains. This distinct characteristic in the spreading from sea ice and water provides the opportunity to monitor sea ice.
The TDS-1 satellite was launched in July 2014 and started its data collection from September 2014. As one of eight payloads onboard on TDS-1, the SGR-ReSI took measurements two days in an eight-day cycle until 2018. The SGR-ReSI was operated in full time mode (7/7 days) during its extension from February to December 2018. The TDS-1 data provide an intense coverage over most of Arctic and Antarctic regions as the satellite runs on a quasi-Sun synchronous orbit with an altitude of~635 km and an inclination of 98.4 • . The TDS-1 data are accessible on the Measurement of Earth Reflected Radio-navigation Signals by Satellite (MERRBys, www.merrbys.co.uk). The available DDMs from TDS-1 consists of 20 Doppler shift bins with an interval of 500 Hz and 128 delay bins with a resolution of 250 ns, which is the length of 0.25 C/A code chips. Figure 1 presents two different DDMs collected over sea water and ice, respectively. It is obvious that the spreading of DDM from sea ice is much less than that of sea water. The reflection of sea ice is more coherent than that of sea water, which results in more scattering in the delay and Doppler domains due to the presence of waves on the open water surface.

Extraction of Features
The scattering components of DDM come from the glistening zone with different delay and Doppler shifts with respect to the specular point. The method proposed in [28] uses the twodimensional delay waveforms generated from DDMs as basic observables for easier data processing. As introduced in [28], the cross section of 20 different Doppler shifts produces 20 delay waveforms, whose summation refers to the integrated delay waveform (IDW) [25] of the DDM over the Doppler domain. The relationship between the power of scattered signals and time delay is illustrated by DDM, which is described by the model proposed in [45] based on the bistatic radar equation: where represents the coherent integration time, τ represents the time delay, represents the function of power antenna footprint, represents the distance from the scattering point to GNSS transmitters, represents the distance from the receiver to the scattering point, Λ is a triangular function as a function of time delay, S is a sinc function in the frequency domain for GPS C/A codes, represents the normalized bistatic radar cross section, represents the Doppler shift frequency and ρ represents the vector from the specular reflection point to the scattering point.
In the TDS-1 mission, the coherent integration time is 1 ms and the Doppler bandwidth can be described by Δ 1 2 ⁄ . If the maximum and minimum Doppler shift of the scattered signal is defined as and , respectively, the width of the glistening zone can be described by . If the Doppler bandwidth is larger than the width of the glistening zone (i.e., ∆ ), the Doppler effects is negligible. The sinc function S is equal to 1 and the cross section with zero Doppler shift (Doppler = 0) is a particularly Central Delay Waveform (CDW) from the DDM. The waveform can be defined as: Another observable termed as differential delay waveform (DDW) was used to describe the degree of difference between CDW and IDW. The DDW between normalized CDW (NCDW) and normalized IDW (NIDW) can be defined as: The IDW is useful to describe the power spreading characteristics due to surface roughness. In order to extract features from DDMs, several data pre-processing schemes presented in the previous study [23,28] should be applied to subtract the noise floor to obtain normalized DDM (NDDM) with

Extraction of Features
The scattering components of DDM come from the glistening zone with different delay and Doppler shifts with respect to the specular point. The method proposed in [28] uses the two-dimensional delay waveforms generated from DDMs as basic observables for easier data processing. As introduced in [28], the cross section of 20 different Doppler shifts produces 20 delay waveforms, whose summation refers to the integrated delay waveform (IDW) [25] of the DDM over the Doppler domain. The relationship between the power of scattered signals and time delay is illustrated by DDM, which is described by the model proposed in [45] based on the bistatic radar equation: where T i represents the coherent integration time, τ represents the time delay, D 2 represents the function of power antenna footprint, R T represents the distance from the scattering point to GNSS transmitters, R R represents the distance from the receiver to the scattering point, Λ is a triangular function as a function of time delay, S is a sinc function in the frequency domain for GPS C/A codes, σ 0 represents the normalized bistatic radar cross section, f D represents the Doppler shift frequency and ρ represents the vector from the specular reflection point to the scattering point.
In the TDS-1 mission, the coherent integration time is 1 ms and the Doppler bandwidth can be described by ∆ f 0 = 1/2T i . If the maximum and minimum Doppler shift of the scattered signal is defined as f max and f min , respectively, the width of the glistening zone can be described by f max − f min . If the Doppler bandwidth is larger than the width of the glistening zone (i.e., ∆ f 0 > f max − f min ), the Doppler effects is negligible. The sinc function S is equal to 1 and the cross section with zero Doppler shift (Doppler = 0) is a particularly Central Delay Waveform (CDW) from the DDM. The waveform can be defined as: Another observable termed as differential delay waveform (DDW) was used to describe the degree of difference between CDW and IDW. The DDW between normalized CDW (NCDW) and normalized IDW (NIDW) can be defined as: The IDW is useful to describe the power spreading characteristics due to surface roughness. In order to extract features from DDMs, several data pre-processing schemes presented in the previous study [23,28] should be applied to subtract the noise floor to obtain normalized DDM (NDDM) with NIDW, NCDW and DDW. Contrary to the previous studies, data with a peak signal-to-noise ratio (SNR) above −3 dB are adopted to increase the amount of data. The more relaxed data filtering strategy is also useful to inspect the applicability and generality of the proposed methods.
There are no effective signals over the several starting and ending delay bins. Therefore, only a part of delay bins from chips −3 to 8.75 (48 delay bins) around the specular point are adopted to extract features. The ground tracks and parts of samples (DDMs and corresponding delay waveforms) of TDS-1 data collected over Baffin Bay on 15 January 2016 are presented in Figure 2. The open water, ice and land are filled with light blue, white and light yellow, respectively. The ground tracks of sea ice and water are depicted by magenta and blue, respectively. The DDMs of sea ice and water are presented with the area marked with cyan and green rectangle respectively. Figure 2a presents the continuous DDMs over the water-ice transition area marked with red rectangle. The corresponding delay waveforms (i.e., NCDW, NIDW and DDW) are shown in Figure 2b. As shown in Figure 2b, the shape of delay waveforms changes from water to ice surface. The largest change in delay waveforms is between DDM 487 and 488.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 22 NIDW, NCDW and DDW. Contrary to the previous studies, data with a peak signal-to-noise ratio (SNR) above -3 dB are adopted to increase the amount of data. The more relaxed data filtering strategy is also useful to inspect the applicability and generality of the proposed methods.
There are no effective signals over the several starting and ending delay bins. Therefore, only a part of delay bins from chips -3 to 8.75 (48 delay bins) around the specular point are adopted to extract features. The ground tracks and parts of samples (DDMs and corresponding delay waveforms) of TDS-1 data collected over Baffin Bay on 15 January 2016 are presented in Figure 2. The open water, ice and land are filled with light blue, white and light yellow, respectively. The ground tracks of sea ice and water are depicted by magenta and blue, respectively. The DDMs of sea ice and water are presented with the area marked with cyan and green rectangle respectively. Figure 2 (a) presents the continuous DDMs over the water-ice transition area marked with red rectangle. The corresponding delay waveforms (i.e., NCDW, NIDW and DDW) are shown in Figure 2 (b). As shown in Figure 2 (b), the shape of delay waveforms changes from water to ice surface. The largest change in delay waveforms is between DDM 487 and 488. A few feature parameters are derived from the delay waveforms to monitor sea ice. Figure 3 presents the NCDW, NIDW and DDW, which can be divided into a left edge (LE) and right edge (RE) according to the point with a delay value of zero. The spaceborne GNSS-R DDMs are generated through cross-correlating scattered signals with code replicas generated locally for different time delays and Doppler shifts [44]. The maximum power is tracked in the Doppler domain to identify the delay value of zero. The earth surface (e.g., sea surface height fluctuations, ice height above the ellipsoid) may affect the geometry and lead to incorrect estimation. This study mainly focuses on the relative change, and not on altimetry applications. Therefore, the impacts of those factors have not been taken into consideration. magenta and blue plots represent the ground tracks of sea ice and water, respectively. The typical DDMs of ice marked with cyan rectangle and water marked with green rectangle are presented. The continuous DDMs from index 481 to 492 for water-ice transition area marked with red rectangle are shown. (b) The continuous delay waveforms of DDM 481 to 492. NCDW, NIDW and DDW are depicted by a blue dotted line, green line and magenta dashed line, respectively.
A few feature parameters are derived from the delay waveforms to monitor sea ice. Figure 3 presents the NCDW, NIDW and DDW, which can be divided into a left edge (LE) and right edge (RE) according to the point with a delay value of zero. The spaceborne GNSS-R DDMs are generated through cross-correlating scattered signals with code replicas generated locally for different time delays and Doppler shifts [44]. The maximum power is tracked in the Doppler domain to identify the delay value of zero. The earth surface (e.g., sea surface height fluctuations, ice height above the ellipsoid) may affect the geometry and lead to incorrect estimation. This study mainly focuses on the relative change, and not on altimetry applications. Therefore, the impacts of those factors have not been taken into consideration. The LE is related to the area above the reflection surface, which results in its insensitivity to the characteristics of reflection surface. Only the RE-related observables are applied to monitor sea ice in this study. Six characteristic parameters termed as RE slope of CDW (RESC), RE slope of IDW (RESI), RE slope of DDW (RESD), RE waveform summation of CDW (REWC), RE waveform summation of IDW (REWI) and RE waveform summation of DDW (REWD) are extracted as features for monitoring sea ice. These features can be computed according to the equations summarized in Table 1. Table 1. The mathematical description of six selected features (i.e., RESC, RESI, RESD, REWC, REWI, REWD). RESC is the right edge slope of CDW. RESI is the right edge slope of IDW. RESD is the right edge slope of DDW. REWC is the right edge waveform summation of CDW. REWI is the right edge waveform summation of IDW. REWD is the right edge waveform summation of DDW.

Features
Mathematical description RESC ∑ ̅ ∑ ̅ The LE is related to the area above the reflection surface, which results in its insensitivity to the characteristics of reflection surface. Only the RE-related observables are applied to monitor sea ice in this study. Six characteristic parameters termed as RE slope of CDW (RESC), RE slope of IDW (RESI), RE slope of DDW (RESD), RE waveform summation of CDW (REWC), RE waveform summation of IDW (REWI) and RE waveform summation of DDW (REWD) are extracted as features for monitoring sea ice. These features can be computed according to the equations summarized in Table 1. Table 1. The mathematical description of six selected features (i.e., RESC, RESI, RESD, REWC, REWI, REWD). RESC is the right edge slope of CDW. RESI is the right edge slope of IDW. RESD is the right edge slope of DDW. REWC is the right edge waveform summation of CDW. REWI is the right edge waveform summation of IDW. REWD is the right edge waveform summation of DDW.

Features
Mathematical Description In the equations in Table 1, n (n ≥ 2) is the number of delay bins for curve fitting and 1 delay bins is equal to 0.25 chips; τ i is the time delay value of the ith point; C R i , I R i and D R i are the waveform values of right edge of CDW, IDW and DDW, respectively; C R , I R and D R are the mean waveform values of points applied for fitting of CDW, IDW and DDW, respectively; τ is the mean of time delay of points applied for fitting. n is set as 5 for RESC, RESI and RESD and 7 for REWC, REWI and REWD.

Validation Data
Two sea ice datasets are used to evaluate the performance of the proposed sea ice monitoring approach. The sea ice edge (SIE) data provided by the Ocean and Sea Ice Satellite Application Facility (OSISAF) are used as the reference data [46,47]. The OSISAF SIE data is generated with a grid resolution of 10 km using a Bayesian approach based on the combination of ASCAT (Advanced Scatterometer) and SSMIS (Special Sensor Microwave Imager Sounder) data with different channels (e.g., 19, 37 and 91 GHz). It is worth noting that the OSISAF data has quality flags which indicate the quality of the sea ice products. The data are divided into five levels according to the confidence levels. The confidence level of 0 means unprocessed, 1 means erroneous, 2 means unreliable, 3 means acceptable, 4 means good and 5 means excellent. The data with a minimum confidence data level of 3 are adopted in this study [46].
The sea ice concentration (SIC) data generated through the Arctic radiation and the turbulence interaction study Sea Ice (ASI) algorithm using AMSR-2 (Advanced Microwave Space Radiometer-2) data are also used as the reference data [4]. This SIC map was obtained from the online sea ice data platform www.meereisportal.de [48]. The reference SIC data are used to generate daily maps in the polar stereographic coordinates with a grid resolution of 6.25 km. The TDS-1 DDMs can be matched with the SIC maps using the location of specular point and date of data collection, which are contained in the data Level 1b. The DDM with a SIC value above 15% is regarded as sea ice, otherwise as sea water [25].

Machine Learning-Aided Sea Ice Monitoring Methods
One of the most important tasks of monitoring sea ice is to distinguish sea ice from water. Therefore, the problem of this study can be regarded as a typical binary classification that can be done by using an ML method on a big dataset. The process flow of monitoring sea ice using ML is presented in Figure 4.
The ML-based sea ice monitoring method includes three steps: (1) feature extraction from the TDS-1 data; (2) the learning process with the training dataset using ML algorithms; (3) the automatic discrimination between data collection over sea ice and over water. A total of seven input variables are used, which are the reference SIC maps and the sequence of six features (i.e., RESC, RESI, RESD, REWC, REWI and REWD) extracted from the TDS-1 data. When the reflection is coherent, the footprint of the TDS-1 DDM is about 6 km by 0.4 km along the track and across track, respectively [25,26,30], which is comparable with the reference SIC maps with a grid resolution of 6.25 km. The footprint is much larger for incoherent reflections. The specular point of each DDM is used to match the reference data. In general, the ML is based on two different data sets (training-set and test-set). The training data are pre-labeled using the reference SIC map; thus, the relationship between input parameters and output results can be modeled using suitable ML algorithms. Then, the output results of test data can be obtained using the pre-built model. According to the process of building models, machine learning could be mainly categorized into supervised learning, unsupervised learning and semi-supervised learning. The characteristics of supervised learning is that training data have priori information (results). In this study, the task is to distinguish sea ice from water and the output results of training datasets can be obtained through the reference SIC data. Therefore, two types of supervised learning-DT and RF-are adopted to monitor sea ice.  Figure 4. Flow diagram of the sea ice monitoring using machine learning (ML). In the first stage (marked by rectangle with black dashed line), the TDS-1 data are processed to extract effective features. In the second stage (marked by a rectangle with a blue line), a classifier is developed using the training data, selected feature sequences (i.e., RESC, RESI, RESD, REWC, REWI and REWD) and ML algorithms, e.g., decision tree (DT) and random forest (RF). In the third stage (marked by a rectangle with a magenta dotted line), the classifier is applied to the test data to generate the sea ice monitoring results and evaluate the performance through comparing with the OSISAF and ASI sea ice data.
The ML-based sea ice monitoring method includes three steps: (1) feature extraction from the TDS-1 data; (2) the learning process with the training dataset using ML algorithms; (3) the automatic discrimination between data collection over sea ice and over water. A total of seven input variables are used, which are the reference SIC maps and the sequence of six features (i.e., RESC, RESI, RESD, REWC, REWI and REWD) extracted from the TDS-1 data. When the reflection is coherent, the footprint of the TDS-1 DDM is about 6 km by 0.4 km along the track and across track, respectively [25,26,30], which is comparable with the reference SIC maps with a grid resolution of 6.25 km. The footprint is much larger for incoherent reflections. The specular point of each DDM is used to match the reference data. In general, the ML is based on two different data sets (training-set and test-set). The training data are pre-labeled using the reference SIC map; thus, the relationship between input parameters and output results can be modeled using suitable ML algorithms. Then, the output results of test data can be obtained using the pre-built model. According to the process of building models, machine learning could be mainly categorized into supervised learning, unsupervised learning and semi-supervised learning. The characteristics of supervised learning is that training data have priori information (results). In this study, the task is to distinguish sea ice from water and the output results of training datasets can be obtained through the reference SIC data. Therefore, two types of supervised learning-DT and RF-are adopted to monitor sea ice.

Decision Tree Algorithm
Decision tree (DT) is one of the simplest and most useful algorithms for classification [49,50]. It has been used to various remote sensing applications [51][52][53]. The structure of a decision tree is In the second stage (marked by a rectangle with a blue line), a classifier is developed using the training data, selected feature sequences (i.e., RESC, RESI, RESD, REWC, REWI and REWD) and ML algorithms, e.g., decision tree (DT) and random forest (RF). In the third stage (marked by a rectangle with a magenta dotted line), the classifier is applied to the test data to generate the sea ice monitoring results and evaluate the performance through comparing with the OSISAF and ASI sea ice data.

Decision Tree Algorithm
Decision tree (DT) is one of the simplest and most useful algorithms for classification [49,50]. It has been used to various remote sensing applications [51][52][53]. The structure of a decision tree is constructed upside down with three parts: internal node, branches and leaf. The first internal node is called the root, where classification starts. The internal node stands for a condition that is expressed by the feature parameters. Based on the node, the decision tree splits into branches according to a discriminant function. The tree ends at the leaf, which represents a final classification decision. Distinguishing between sea ice and water can be regarded as a binary classification problem. Thus, the algorithm C4.5 [54] is used, which recursively splits training data into subdivisions using a set of attributes described by input variables. C4.5 builds decision trees from a set of training data using the concept of information entropy. The training data are a set S = s 1 , s 2 , . . . s i of already classified samples. Each sample s i consists of a p-dimensional vector x 1,i , x 2,i , . . . , x p,i , where the x i values represent attribute values or features of the sample, as well as the class in which s i falls. C4.5 uses the information gain ratio to construct a decision tree. The information gain ratio is defined as: where m is the number of categories; v is the number of selected features; D is the number of samples; D j is the jth sample. In this paper, only two categories, i.e., sea ice and sea water, are included, so m = 2.
As six features are selected, v is equal to 6. The information gain ratio can be simplified as: C4.5 has several advantages. It can mitigate overfitting through single pass pruning process, handle both discrete and continuous data and address the problem of incomplete data, which is common in practical applications.

Random Forest algorithm
Another ensemble learning method for classification is random forest (RF), which constructs a collection of DT at training time. RF combines a boosting sampling strategy and Classification and Regression Tree (CART) to overcome the drawback of a single CART, such as overfitting problems. CART uses a Gini index [55] to measure the impurity of training datasets, while C4.5 utilizes the concept of entropy. The Gini index is described by: where s is the number of categories and p l is the proportion of samples belonging to class l. Since sea ice monitoring is a binary classification problem. Thus, the Gini index can be simplified as: where p can be regarded as the probability that samples belong to sea ice. The advantages of CART include that the rules can be interpreted easily and that it provides automatic processing of parameters selection, data missing, outliers, variable interaction and nonlinear relationships. However, one of the biggest shortcomings of a single CART is overfitting. The strategy of bagging can effectively solve the problem through constructing a large number of independent trees and reduce errors that may be caused by some unstable classifiers [56]. Due to its advantages, RF shows great potential in many remote sensing applications [57].

Results
The TDS-1 data collected over the Arctic and Antarctic regions with the latitude above 55 • N and 55 • S from January 2015 to December 2018 are analyzed in this study. Twenty percent of the data is randomly selected to train the ML-based models to distinguish sea ice from water. The remaining 80% of data is used as the test dataset to validate the sea ice monitoring methods developed using ML algorithms.
As aforementioned, the GNSS-R receiver on the TDS-1 satellite was not always in operation. Thus, the TDS-1 data are not accessible every day. Figure 5 presents the situation of data availability from January 2015 to December 2018. The data unavailability from August to October 2017 is probably due to the scheduled shutdown of TDS-1 mission, which was originally set to the end of July 2017. In fact, the TDS-1 mission was extended from February to December 2018. During its extension, the SGR-ReSI was operated every day, rather than the two of eight-day cycle as in the first three years. The coverage and sampling were increased by a factor of four. The data missing for a few days in 2018 may result from statutory holidays, such as Christmas.
from January 2015 to December 2018. The data unavailability from August to October 2017 is probably due to the scheduled shutdown of TDS-1 mission, which was originally set to the end of July 2017. In fact, the TDS-1 mission was extended from February to December 2018. During its extension, the SGR-ReSI was operated every day, rather than the two of eight-day cycle as in the first three years. The coverage and sampling were increased by a factor of four. The data missing for a few days in 2018 may result from statutory holidays, such as Christmas.

Characteristics of GNSS-R Features
The distribution characteristics of six feature parameters (RESC, RESI, RESD, REWC, REWI and REWD) for sea ice and water are shown in Figure 6. The vertical height of the boxes represents the interquartile range of the samples, while the parallel line depicted in red inside the boxes is the median value of the samples for each feature. The green dotted line represents the threshold obtained by the method proposed in [28] for distinguishing sea ice from water. It is clear that sea ice shows a distinct difference between sea water for all the features considered. This is because the reflection over the sea ice surface is usually more coherent. The sea water surface is often rougher than that of sea ice and easily affected by ocean winds, which results in wider scattering. The median values of each parameter for sea ice and water are significantly different. However, the distribution of features of sea ice and water is more or less overlapped. As shown in Figure 6 (a), the RESC values of sea ice range from 0.15 to 1 and those of sea water range from 0.02 to 0.99; the threshold is 0.745. If RESC < 0.745, it is regarded as sea water; if RESC > 0.745, it is regarded as sea ice. However, some points of RESC below 0.745 appear in sea ice, and some points above 0.745 appear in sea water; these points are overlap. This indicates that simple thresholding of each feature may result in some false discrimination between sea ice and water.

Characteristics of GNSS-R Features
The distribution characteristics of six feature parameters (RESC, RESI, RESD, REWC, REWI and REWD) for sea ice and water are shown in Figure 6. The vertical height of the boxes represents the interquartile range of the samples, while the parallel line depicted in red inside the boxes is the median value of the samples for each feature. The green dotted line represents the threshold obtained by the method proposed in [28] for distinguishing sea ice from water. It is clear that sea ice shows a distinct difference between sea water for all the features considered. This is because the reflection over the sea ice surface is usually more coherent. The sea water surface is often rougher than that of sea ice and easily affected by ocean winds, which results in wider scattering. The median values of each parameter for sea ice and water are significantly different. However, the distribution of features of sea ice and water is more or less overlapped. As shown in Figure 6a, the RESC values of sea ice range from 0.15 to 1 and those of sea water range from 0.02 to 0.99; the threshold is 0.745. If RESC < 0.745, it is regarded as sea water; if RESC > 0.745, it is regarded as sea ice. However, some points of RESC below 0.745 appear in sea ice, and some points above 0.745 appear in sea water; these points are overlap. This indicates that simple thresholding of each feature may result in some false discrimination between sea ice and water.
This study uses the combination of six features derived from the delay waveforms of different Doppler spread characteristics to describe the characteristics of reflecting surface. The six features of samples are composited into sequences, which are applied as input variables to train the sea ice monitoring model. The six features are combined into sequences in order. RESC, RESI, RESD, REWC, REWI and REWD values are presented from bottom to top in the y-axis. The feature sequences of samples in the Arctic and Antarctic regions are presented in Figure 7.  As shown in Figures 7 (a) and (b), the feature sequences of 40,000 samples for sea ice (upper plot) and water (lower plot) show distinct differences, which provide the opportunities of monitoring sea ice. Moreover, the feature sequences can describe the characteristics more accurately than individual features. As shown in Figure 7a,b, the feature sequences of 40,000 samples for sea ice (upper plot) and water (lower plot) show distinct differences, which provide the opportunities of monitoring sea ice. Moreover, the feature sequences can describe the characteristics more accurately than individual features.

Sea Ice Monitoring Performance
The sea ice monitoring models based on DT and RF algorithms are quantitatively assessed using confusion matrices [58] through a comparison with the OSISAF SIE data using the test data. In the field of machine learning and specifically the problem of statistical classification, a confusion matrix [56], also known as an error matrix, is a specific table layout that allows visualization of the performance of a supervised learning algorithm. The confusion matrix is a table with two rows and two columns that reports the number of false positives, false negatives, true positives and true negatives. The error matrices, overall accuracy and kappa coefficient [59] of the agreement are used as indicators to evaluate the performance of the DT and RF models. The performance of the DT and RF models for the Arctic and Antarctic regions are presented in Tables 2 and 3, respectively.
The overall accuracy of DT model is 97.51% and 95.46% for the Arctic and Antarctic, respectively, while the RF model produced an overall accuracy of 98.03% and 95.96% for the Arctic and Antarctic, respectively. The producer and user accuracies of sea water are higher than those of sea ice for both models. This may be because the sea ice with a low SIC is more easily misidentified as sea water. When the surface area with both ice and water is driven by wind field, the surface will become rougher, and the sea ice surface is recognized as sea water. Although the DT and RF models obtain similar overall accuracies, the Kappa coefficient of agreement of RF model is slightly higher than that of DT, which indicates that the performance of the RF algorithm is better than that of DT. Although the overall accuracy obtained in this study is slightly lower than that in the previous study, the dataset used here is much larger and the data filter requirement is lower. This indicates the method developed and applied here is of better applicability and generality. When using data only from the initial mission, as we did in our previous study [25], the overall accuracy of this method is 0.22% better than the REWD method we applied there.
The previous study [37] indicated that the support vector machine (SVM) outperforms the neural network (NN) and convolutional neural network (CNN) methods for detecting sea ice using spaceborne GNSS-R data. SVMs are capable of operating classification tasks by finding a hyperplane that can best distinguish (with the maximum margin) between different types. NNs are extremely flexible in the types of data they can support. NNs do a decent job at learning the important features from basically any data structure, without having to manually derive features. CNNs are much less flexible models compared to a fully connected network, and are biased toward performing well on image. In order to evaluate the performance of proposed methods, the SVM is adopted for comparison. The sea ice monitoring results obtained by SVM based methods are shown in Table 4. The proposed RF-based sea ice monitoring approach shows better accuracy than the SVM-based method, while the SVM-based sea ice monitoring scheme outperforms the DT-based one. The feature sequences applied in this study are extracted from delay waveforms (NCDW, NIDW and DDW) with different doppler shifts.

Discussion
For further analysis, the time series of overall accuracy of sea ice monitoring is computed using all the available data from January 2015 to December 2018 ( Figure 8). The overall accuracy of the Arctic region is significantly lower in September 2016 since the sea ice melts in this season, while the changing trend of the Antarctic region is reverse as the seasonal alternation between the Arctic and Antarctic is opposite.

Discussion
For further analysis, the time series of overall accuracy of sea ice monitoring is computed using all the available data from January 2015 to December 2018 ( Figure 8). The overall accuracy of the Arctic region is significantly lower in September 2016 since the sea ice melts in this season, while the changing trend of the Antarctic region is reverse as the seasonal alternation between the Arctic and Antarctic is opposite. To analyze the impact of each variable, the relative importance of variables for sea ice monitoring is shown in Figure 9. REWD is used at all nodes in the DT algorithm, which results in a relatively high contribution to sea ice monitoring. REWI is useful as it can be used to distinguish sea ice from water with very low error. The REWD is the most important parameter, followed by REWI, RESI, RESD, RESC and REWC in the DT algorithm. Like the DT algorithm, REWC is of the least significance for monitoring sea ice in the RF algorithm.  To analyze the impact of each variable, the relative importance of variables for sea ice monitoring is shown in Figure 9. REWD is used at all nodes in the DT algorithm, which results in a relatively high contribution to sea ice monitoring. REWI is useful as it can be used to distinguish sea ice from water with very low error. The REWD is the most important parameter, followed by REWI, RESI, RESD, RESC and REWC in the DT algorithm. Like the DT algorithm, REWC is of the least significance for monitoring sea ice in the RF algorithm. To analyze the impact of each variable, the relative importance of variables for sea ice monitoring is shown in Figure 9. REWD is used at all nodes in the DT algorithm, which results in a relatively high contribution to sea ice monitoring. REWI is useful as it can be used to distinguish sea ice from water with very low error. The REWD is the most important parameter, followed by REWI, RESI, RESD, RESC and REWC in the DT algorithm. Like the DT algorithm, REWC is of the least significance for monitoring sea ice in the RF algorithm.  The RF-based GNSS-R sea ice monitoring results in March and September 2018 are mapped with the OSISAF SIE data in Figure 10. The white and dark gray edge (partly marked by a rounded rectangle with a red dotted line) represent the minimum and maximum ice extent for March and September in 2018, respectively. In March, the sea ice extent of Arctic region reaches the minimum and maximum on 6 and 14 March, respectively, while the minimum and maximum sea ice extent of Antarctic region appear on 1 and 31 March respectively. As shown in Figure 5d, the data are not available every day in September 2018. From 1 to 17 September, the maximum and minimum sea ice extent of Arctic region occur on 2 and 17 September, respectively, while the sea ice extent of Antarctic region reaches the minimum and maximum on 13 and 17 September, respectively. The scatter points are ground tracks of TDS-1 data with the peak SNR above −3 dB, which results in some gaps in the GNSS-R ground-tracks. In the figures, the presence of sea ice monitored using GNSS-R is illustrated by magenta points, whereas the presence of GNSS-R sea water is depicted by the blue points. As shown in Figure 10b,c, the detected sea ice and water overlaps in some areas. This is because the GNSS-R data span over one month and the ice extent changes rapidly during the melting season in the Arctic and Antarctic regions respectively.  The examples of monitoring sea ice around Greenland using four different methods are presented in Figure 11. The sea ice monitoring results are compared with the ASI SIC data. Two simple thresholding methods based on REWD (i.e., REWD > 0.38 for sea water and REWD < 0.38 for sea ice) and REWI (i.e., REWD > 0.62 for sea water and REWD < 0.62 for sea ice) used in [28] are adopted to monitor sea ice (Figure 11a,b). The simple thresholding methods result in some false monitoring of sea ice; sea ice is identified as sea water or sea water is regarded as sea ice. Although REWD and REWI are considered as useful parameters for distinguishing sea ice from water, simple thresholding based on just one parameter was shown to be insufficient for effectively monitoring sea ice. The results of DTand RF-based approaches are presented in Figure 11c,d, respectively. The false sea ice monitoring of DT-and RF-based methods mainly appear around the sea ice edge areas with a relatively low SIC. The area with a low SIC may be affected by ocean winds, which results in a rougher surface. Then, the sea ice is wrongly identified as sea water. The effects of ocean winds on low SIC have not been analyzed in this study.

Conclusions
In this study, two machine learning-aided GNSS-R methods have been proposed to monitor sea ice using 42 months of TDS-1 data. The sea ice monitoring results are validated with the SIE data from OSISAF. The results showed that the proposed approach successfully distinguishes sea ice from water. The proposed RF-and DT-based sea ice monitoring approaches achieve an overall accuracy of 98.03% and 97.51%, respectively, in the Arctic regions, and 95.96% and 95.46%, respectively, in the Antarctic regions. Another ML-based method (i.e., SVM) used in the previous study [40] is also applied for comparison in this study. The SVM-based method achieves an overall accuracy of 97.62% and 95.61%, respectively, in the Arctic and Antarctic regions with the dataset used in this study.
A total of six features were combined to monitor sea ice, including RESC, RESI, RESD, REWC, REWI and REWD. Although these features have been applied to sense sea ice individually in the Figure 11. Examples of sea ice monitoring results validated against ASI SIC (sea ice concentration) maps from AMSR-2 data on the southwest side of Greenland on 14 March 2018 using four different methods: (a) the REWD thresholding approach, (b) the REWI thresholding approach, (c) the DT-based method in this study and (d) the RF-based method in this study. The land and sea water are represented as light brown and white, respectively. The sea ice concentration (SIC) is demonstrated by the color bar. The green and blue points represent the detected sea ice and sea water, respectively, while the red points represent the false detection.

Conclusions
In this study, two machine learning-aided GNSS-R methods have been proposed to monitor sea ice using 42 months of TDS-1 data. The sea ice monitoring results are validated with the SIE data from OSISAF. The results showed that the proposed approach successfully distinguishes sea ice from water. The proposed RF-and DT-based sea ice monitoring approaches achieve an overall accuracy of 98.03% and 97.51%, respectively, in the Arctic regions, and 95.96% and 95.46%, respectively, in the Antarctic regions. Another ML-based method (i.e., SVM) used in the previous study [40] is also applied for comparison in this study. The SVM-based method achieves an overall accuracy of 97.62% and 95.61%, respectively, in the Arctic and Antarctic regions with the dataset used in this study.
A total of six features were combined to monitor sea ice, including RESC, RESI, RESD, REWC, REWI and REWD. Although these features have been applied to sense sea ice individually in the previous study, the combination of these six features is firstly adopted to monitor sea ice. Compared to the single observable method, the feature sequences can represent the characteristics of reflecting surface more accurately. Therefore, the ML-based approaches achieve higher accuracies than the single observable thresholding method. It would be worth noting that the input features to ML-based methods are different from the single observable thresholding method. Moreover, the spaceborne GNSS-R dataset used here spans 42 months of the TDS-1 mission, which is larger than those applied in the previous studies. The results from this study are encouraging for the GNSS-R applications of machine learning algorithms. Further research on the effects of oceans winds in the low SIC regions will benefit monitoring sea ice. In addition, the combination of multiple ML-based methods (e.g., DT, RF and SVM) will be explored in our future work.