Isometric mapping algorithm based GNSS-R sea ice detection

1.   Introduction

Sea ice is a critical component of the Earth's climate system and plays a significant role in global ocean and atmospheric circulations ^[1]. Sea ice information provides knowledge and plays a beneficial role in ocean voyages and natural resource exploration. Sea ice and open water (OW) show different characteristics in roughness ^[2], which is one important parameter for sea ice detection. Compared to OW, sea ice has higher albedo and roughness, which play an important role in energy exchanging between sea and air. The reduced extent of Arctic sea ice and thinner ice cover have been indicated in some previous studies ^[3]. The roughness and dielectric constant of ocean surface change with the variation of the ground surface type, such as the appearance of sea ice. The surface of sea ice presents more smoother than that of OW. These characteristics of OW and sea ice are the basis to detect sea ice.

During recent research, global navigation satellite systems reflectometry (GNSS-R) has played a powerful role in using L-band signals scattered from the Earth's surface to sense bio-geophysical features ^[4]. The initial GNSS-R application was ocean altimetry detection ^[5] after the concept of GNSS-R was proposed in 1988 ^[6]. Subsequently, the applications of GNSS-R have been extended to wind speed retrieval ^[7,8,9], snow depth estimation ^[10], soil moisture sensing ^[11,12], sea target detection ^[13] and sea ice detection ^{[14,15,16,17]}.

Since the greatest amplitude of specular scattering is presented from the Earth's surface GNSS reflected signals, the specular scattering geometry can be used in the applications of GNSS-R. As one of the most important GNSS-R observables, the Delay-Doppler Map (DDM) ^[18] is a function of time delays and Doppler shifts. The function is to describe the power scattering from the reflected surface. Through integrating DDM in the Doppler domain, the integration delay waveform (IDW) in Figure 1 can be obtained as another GNSS-R observable. The reflection over open water (OW), whose surface is rough, often has a continuous pixel number jump like the blue line in Figure 1. The sea ice surface is often considered as relatively smooth, and the IDW is shown as the red line in Figure 1, which has only one pixel number jump.

Figure 1.  Classical IDWs of sea ice and OW.

DownLoad: Full-Size Img PowerPoint

In a recent GNSS-R sea ice detection development, Yan et al. ^[16] successfully utilized the number of DDM pixels with signal-to-noise ratio (SNR) above a threshold for sea ice remote sensing. Zhu et al.^[19] recognized the transitions of ice and water with the differential DDM observable. Yan et al. ^[20] classified the sea ice and open water with convolutional neural networks. Similarly, Hu et al. ^[21] used IDW to detect sea ice and OW with the linear discriminant algorithm (LDA) method and analyzed the noise impact in sea ice detection. DDM has been widely applied as reference data to detect sea ice. sea ice studies using IDW as reference data are scarce, the observable of IDW has more potential to sense sea ice. IDW is utilized as research data to detect sea ice in this study, and the isometric mapping (ISOMAP) is proposed to obtain the observable of IDW and use four different classifiers to classify the low-dimensional features into sea ice and OW. These classifiers are based on support vector machine (SVM), gradient boosting decision tree (GBDT), LDA and K-nearest neighbors (KNN). Compared with previous studies, the whole IDW is used as data in this paper, instead of leading edge slope (LES) or trailing edge slope (TES), as observations for sea ice detection.

This paper is structured as follows: Section 2 introduces data sources and methods. Section 3 presents the results of the ISOMAP-based sea ice detection in this paper. The discussion is presented in Section 4, and Section 5 is the conclusions.

2.   Data sources and ISOMAP-based method

2.1.   Data sources

2.1.1.   TechDemoSat-1 IDW data

TechDemoSat-1 (TDS-1) can provide the spaceborne GNSS-R DDM for spaceborne GNSS-R study. In 2014, TDS-1 was launched as a national technology demonstration satellite. In the all regions of the world, the related datasets can be found at the website (http://www.merrbys.co.uk).

In ^[18], the power of scattered signals was described as time delay function, the function was proposed as theoretical integrated delay waveform (IDW):

where $\tau$ and ${f}_{D}$ are the time delay value and the Doppler shift frequency, $\overrightarrow{\rho }$ and ${T}_{i}$ are the scattering area pixels and the coherent integration time, $\lambda$ is the wavelength of the L1 signal, ${P}_{t}$ is the GNSS transmitting power, $\chi$ is Woodward's ambiguity function (WAF) ^[22], ${G}_{t}$ and ${G}_{r}$ are transmitter antenna gain and receiver antenna gain, ${R}_{t}$ is the distance from the transmitter to the surface point (SP), ${R}_{r}$ is the distance from the receiver to SP, ${D}^{2}$ is the function of power antenna footprint, and ${\sigma }_{0}\left(\overrightarrow{\rho }\right)$ is the normalized bistatic radar cross section. The reflection over rough open water (OW) surfaces often has a continuous pixel number jump. The IDW of the considered relatively smooth sea ice surface has only one pixel number jump.

Because of the opportunity in GNSS-R sea ice classification, 625555 IDWs were utilized in this paper. In detail, a consecutive period from February to April 2018 was chosen to select IDWs for detecting sea ice continuously. Due to the presence of massive amounts of sea ice at high latitudes, the experiment was chosen for areas above 70°N latitude. These measures can provide the benefits in sea ice detection study. The DDM metadata provides each DDM's specular point position, which is also utilized to calculate the correlation with the surface data from the National Snow and Ice Data Center (NSIDC). In ^[23], details about TDS-1 DDM metadata can be found for the detection experiment.

2.1.2.   Surface type data from NSIDC

The surface type data are provided from the NSIDC. The surface type data are utilized as true type of IDWs in the ISOMAP-based sea ice detection experiment. NSIDC focuses on Earth's atmospheric and oceanic changes, catastrophic weather warnings, and nautical and aeronautical charts. The NSIDC website (https://nsidc.org) provides the available ground data. Based on latitude and longitude, two-dimensional matrix is the present of the surface-type. The NSIDC spatial resolution is 1 km.

2.2.   ISOMAP-based method

The ISOMAP-based method is based on the isometric mapping and classifier. The low-dimensional observed features of IDW are extracted by using ISOMAP, and then the classifiers such as SVM, GBDT, LDA, KNN are used to classify the low-dimensional features of IDW into sea ice or OW.

2.2.1.   Isometric mapping

Isometric mapping ^[24] is an improved function of multidimensional scaling (MDS) ^[25], overcoming the difficulties of MDS through defining a new metric and substituting that for Euclidean distance. The IDWs can preserve the relative stability of the feature when mapping to low-dimensional features. ISOMAP implements a neighborhood graph that connects the closer feature points. ISOMAP takes the following steps:

In the first step, the neighbors of each sample x_i on the low dimensional manifold M based on some appropriate distance metrics d_x(x_i, x_j) in input space X is determined. A k nearest neighbor algorithm is used by ISOMAP to determine neighbors. These neighborhood relationships are represented in a weighted graph G in which d_G(x_i, x_j) = d_X(x_i, x_j), if x_i and x_j are neighbors, and ${d}_{X}\left({x}_{i}, {x}_{j}\right) = \infty$ otherwise.

Estimate distance d_M(x_i, x_j) between any pair of points on the manifold M is the second step. Due to the unknown embedding manifold, d_M(x_i, x_j) approximates the shortest path between x_i and x_j on G, which is calculates by the Floyd-Warshall algorithm ^[26]:

The shortest paths between any two points are represented in a matrix D where ${D}_{ij} = {d}_{G}({x}_{i}, {x}_{j})$. The last step is to apply classical MDS to obtain the matrix of distance, and the output is the low-dimensional features of IDW.

2.2.2.   Support vector machine

As a classical machine learning algorithm, support vector machine (SVM) ^[27] is often used to determine the type boundaries and the basis of SVM is statistical learning theory. The binary classification problems are solved with SVM. The input data can be presented as

where x_i are the classification features of input data, and y_j is the type label of x_i. The SVM follows the following rule in the linear classification problem:

where ${\omega }^{T}x+b = 0$ is a hyperplane, $\omega$ is the coefficient parameter, and b is the coefficient and bias parameter. The classification algorithm is based on maximum interval, which can be expressed as:

Then, the parameters $\omega$ and b can be obtained to classify the input data.

2.2.3.   Gradient boosting decision tree

With improving the boosting tree algorithm, Gradient boosting was proposed successfully^[28]. It is composed with boosting and gradient descent. In order to get the optimal model, the gradient descent is performed in function space.

The first step is to initialize the weak learner. The weak learner is shown as follows:

where ${y}_{i}\in \left\{\mathrm{0, 1}\right\}$. The next step is to calculate the negative gradients r_mi at each gradient m ($m = 1, \cdot \cdot \cdot, m$).

r_mi obtained in first step is utilized as the new true sample value to compute a new regression tree ${f}_{m}\left(x\right)$ and corresponding leaf node area ${R}_{jm}$, . Then, the best-fit values of each corresponding leaf node area for the leaf region will be calculated as follows:

where $j = 1, \cdot \cdot \cdot, J$, and J is the number of leaf nodes of the regression tree. The strong learner will be updated.

After m times of updating the learner, the final learner can be obtained.

2.2.4.   Linear discriminant algorithm

The dimensionality of data can be reduced by the linear discriminant algorithm (LDA) ^[29], and the types of data can be determined with these low-dimensional features. The projected vector is calculated as follows:

where $x$ is n-dimensional column vectors, A is an $m\cdot n$ vector, $y$ is the m-dimensional projected vector. A can consist of 0 matrix and experiment data, such as $A = [O, IDW{]}^{T}$.The best projection vector $x$ will be obtained with the Fisher projection criterion $J\left(x\right)$. The total scatter of the projected samples can be characterized by the trace of the covariance matrix of the projected features vector ^[30] in the theory of the Fisher projection criterion. $J\left(x\right)$ is calculated as follows:

where $x$ is the unitary column vector. ${S}_{b}$ is the maximum inter-class scattering matrix, ${S}_{w}$ is the minimum intra-class scattering matrix. Due to the nonsingular property of ${S}_{w}$, a generalized eigenvalue problem ^[30] is considered to present the optimization problem. ${S}_{b}$ and ${S}_{w}$ can be calculated as follows:

where ${\mu }_{0}$ and ${\mu }_{1}$ denote sample means of two different type data sets, ${X}_{0}$ and ${X}_{1}$ are IDW sample sets of two different types, and ${x}_{s}$ denotes IDW sample in ${X}_{0}$ or ${X}_{1}$. When unitary vector $x$ maximizes $J\left(x\right)$, the Fisher optimal projection axis can be obtained. Find $x$ such that it equals 0 and the relationship between $J\left(x\right)$ and $x$, $x$ can be obtained as follows:

With maximizing criterion, the optimal projection ${x}_{opt}$ and the input optimal feature vector can be obtained as follows:

Euclidean distance is employed to classify the features into sea ice and OW. The Euclidean distance ^[31] is utilized to characterize similarity.

where ${Y}_{1} = \left[{y}_{{1}_{1}}, \cdots, {y}_{{1}_{d}}\right]$ and ${Y}_{2} = \left[{y}_{{2}_{1}}, \cdots, {y}_{{2}_{d}}\right]$ are the input feature matrices of different types. Each IDW is categorized as either ${T}_{1}$ or ${T}_{2}$.

When $d(Y, {Y}_{l}) = \mathit{min}(d({Y}_{i}, {Y}_{l}))$, the classification can be calculated as follows:

2.2.5.   K-nearest neighbors

K-nearest neighbors (KNN) ^[32] is a basic algorithm of machine learning which is usually used in applications of classification and regression. KNN uses the distance between different feature values to classify the data. The first step for KNN is determine a distance measurement method, such as Euclidean distance. Then, the k nearest samples of $x$ in training sets, $T = \left\{\left({x}_{1}, {y}_{1}\right), \cdot \cdot \cdot, \left({x}_{N}, {y}_{N}\right)\right\}$, are found to construct a new set N_k(x). The n-dimensional vector of a sample is y_i, and ${y}_{i}\in Y = \{{c}_{1}, \cdot \cdot \cdot, {c}_{K}\}$. The type of the sample $x$ is determined with the principle of majority voting like follows:

where $i = 1, \cdot \cdot \cdot, N$, $j = 1, \cdot \cdot \cdot, K$, and $I$ is the indication function

2.2.6.   Sea ice detection based on ISOMAP and classifiers

In this paper, the type of the ground surface is sea ice and OW. The classifiers of SVM, GBDT, LDA and KNN were proposed to classify observations such as DDM features and SNR as sea ice or OW. Four classifiers with different mechanisms are used to classify IDW features acquired by ISOMAP. SVM is a supervised classification discriminative algorithm, GBDT is a supervised greedy discriminative algorithm, LDA is a supervised dimensionality reduction algorithm, and KNN is a supervised lazy learning algorithm. The performance evaluations are compared based on the evaluated quantities in Table 1. The process flow of the experiment is presented in Figure 2.

Table 1.  The definitions of evaluation metrics.

Evaluation Metric	Function
Accuracy (%)	$\frac{{TP + TN}}{{TP + TN + FP + FN}}$
Precision (%)	$\frac{{TP}}{{TP + FP}}$
Recall (%)	$\frac{{TP}}{{TP + FN}}$
F1-value	$\frac{{\text{Re} call \cdot \Pr ecision \cdot (1 + {\beta ^2})}}{{\beta \cdot (\text{Re} call + \Pr ecision)}}, \beta = 1$
Kappa coefficient	$\frac{2·(TP·TN-FN·FP)}{(TP+FP)·(FP+TN)+(TP+FN)·(FN+TN)}$

---

 | Show Table

DownLoad: CSV

Figure 2.  Process flow of ISOMAP-based method.

DownLoad: Full-Size Img PowerPoint

In step 2, low-dimensional features are extracted from IDWs, and 30% of samples are randomly selected as the training set. The remaining 70% of samples are used as the test set after low-dimensional feature extraction.

In step 3, the low-dimensional features are applied to training the classifiers of SVM, GBDT, LDA and KNN, and the remaining samples' low-dimensional features are used as test sets to classify sea ice and OW.

In step 4, the detection performance is analyzed and evaluated. The accuracy, precision, recall, F1 value, kappa coefficient ^[33] and confusion matrix are analyzed in the results.

3.   Results

As illustrated in Section 2.2.6, the sea ice detection is conducted with a two-step method. The first step aims to extract low-dimensional features from IDWs using ISOMAP. After extracting low-dimensional features, the low-dimensional features are employed for sea ice and OW detection using SVM, GBDT, LDA and KNN classifiers. There, 30% of the low-dimensional features of previously selected samples are used as the training set, and the remaining 70% of low-dimensional features are used as the test set to classify sea ice and OW.

The ISOMAP-based confusion matrices of RF, SVM, GBDT and LDA classifiers are presented in Figure 3, which presents the classification results of each class. The evaluation metrics are computed with the equations listed in Table 1 and shown in Table 2.

Figure 3.  Confusion matrices of ISOMAP-based method: (a) SVM, (b) GBDT, (c) LDA, (d) KNN.

DownLoad: Full-Size Img PowerPoint

Table 2.  Evaluation metrics of ISOMAP-based sea ice detection.

Evaluation metric	ISOMAP-based detection
	SVM	GBDT	LDA	KNN
Accuracy (%)	99.44	85.58	91.88	98.82
Precision (%)	97.42	58.43	71.59	94.66
Recall (%)	99.88	99.93	99.60	99.81
F1-value	0.98	0.74	0.78	0.96
Kappa coefficient	0.98	0.65	0.83	0.97

---

 | Show Table

DownLoad: CSV

The evaluation metrics are shown in Table 2. The predicted results are validated with the test dataset, whose true types are selected from the NSIDC surface data. In Figure 3 (a), (b), (c) and (d), 437889 samples are selected for testing, in which there are 89038 OW and 348851 sea ice samples. With SVM, GBDT, LDA and KNN classifiers, the ISOMAP-based sea ice detection obtains accuracy values of 99.44%, 85.58%, 91.88% and 98.82%. The accuracy is comparable to previous studies ^[34,35,36].

The ISOMAP-based OW-sea ice detection has comparable performance through SVM, GBDT, LDA and KNN classifiers. The detection results for February 2018 are shown in Figure 4 to present the overall space distribution.

Figure 4.  Sea ice detection of ISOMAP-based method for February 2018: (a) SVM, (b) GBDT, (c) LDA, (d) KNN.

DownLoad: Full-Size Img PowerPoint

The overall spatial distribution of the February 2018 classification results is shown in Figure 4, which demonstrates the distribution of predicted and reference types.

4.   Discussion

In order to analyze the errors in ISOMAP-based sea ice detection, the low-dimensional features classified with classifiers are presented in Figure 5. The errors appear on the distribution boundaries of low-dimensional data, which is caused by the noise distribution in the waveform of IDW. The greater the degree and number of data fluctuations in IDW, the more difficult the IDW waveform is for ISOMAP-based sea ice detection.

Figure 5.  Low-dimensional features of ISOMAP-based method: (a) SVM, (b) GBDT, (c) LDA, (d) KNN.

DownLoad: Full-Size Img PowerPoint

Another experiment is implemented to analyze the relationship between the misclassifications and the degree and number of data fluctuations in IDW. First, the number of data fluctuations in IDW is calculated and counted. Within three consecutive time delay sequences of IDW, if the difference between the maximum pixel power and the minimum pixel power is greater than 1000 pixels, it will be counted as one data fluctuation.

As shown in Figure 6, the data fluctuation numbers in predicted IDWs are counted to present dramatic fluctuations in the IDW waveform. From Figure 6 (b), (c), the data fluctuation values of wrong predicted results are mainly concentrated in the range of 5 to 30 and the range of 0 to 1. However, the data fluctuation values of wrong predicted results are mainly concentrated only in the range of 5 to 30 for Figure 6(a), (d), and the data fluctuation values of correct predicted results are mainly concentrated in the range of 0 to 5 for Figure 6. This shows that the data with too many dramatic fluctuations can be filtered out to improve the quality of data and the GNSS-R sea ice detection. Therefore, future experiments need to not only screen at SNR but also use data fluctuations to ensure the quality of the DDMs or IDWs.

Figure 6.  Data fluctuation counts of ISOMAP-based method: (a) SVM, (b) GBDT, (c) LDA, (d) KNN.

DownLoad: Full-Size Img PowerPoint

5.   Conclusions

In this paper, the ISOMAP-based method is proposed to use IDWs for GNSS-R sea-ice detection. The experiments show that the ISOMAP-based method can be used to conduct sea ice detection and obtain great classification accuracy. Above 70°N and during February-April 2018, the feasibility of ISOMAP-based GNSS-R sea ice detection with selected IDWs is verified with the accuracy of 99.44%. The errors in ISOMAP-based sea ice detection have been analyzed, and the analysis shows that the data fluctuation numbers can reflect the quality of the data. The data predicted incorrectly are concentrated in the range of 5 to 30, and the data predicted correctly are concentrated in the range of 0 to 5. Therefore, the data fluctuation number can be used as another important parameter to improve the GNSS-R sea ice detection by selecting the higher quality data.

Acknowledgements

This work was sponsored by the National Natural Science Foundation of China (52071199).

Use of AI tools declaration

The authors declare they have not used artificial intelligence (AI) tools in the creation of this article.

Conflict of Interest

The authors declare that they have no conflict of interest.