Superpixel Classification with the Aid of Neighborhood for Water Mapping in SAR Imagery

1. Introduction

Extracting water regions in radar images has been an important field in remote sensing, given the fact that water segmentation has many real-world applications such as ice/snow/oceans monitoring and change detection (CD) in water in dams, lakes, and rivers. It is significant to map water regions promptly and accurately and provide the information about the water regions, especially when a natural disaster like flooding occurs. Nowadays Synthetic Aperture Radar (SAR) data are widely used for providing imagery for water mapping, as the radar technology offers customers weather condition independent, 24/7 available, and high-resolution satellite imagery.

Multi-temporal SAR imagery, where images of a region are made on several diffrent dates [1], can utilize CD as a methodology for flood detection; it uses not only flood imagery but also its pre-/post-flood imagery as reference for flood mapping [2,3]. The behavior of surface water, which shows high temporal variability (the standard deviation of the backscattered intensities) and low minimum backscatter in time series, enables us to use a simple threstholding approach for segmenting water accurately [4]. It is equally true, however, that multi-temporal imagery is not available to be used at any time; we need at least two images (ideally produced under the same conditions except for date) for water mapping. Single-date SAR imagery’s larger uncertainty leaves it behind multi-temporal in accurate extraction, but it is more economical of resources and time; it has therefore an advantage in a situation which needs fast action.

In general, electromagetic (plane) waves have a degree of freedom, e.g. a rotation around the propagation axis. This allows us to have linearly polarization which is eigher vertical (V) or horizontal (H). Considering both transmitted waves and received, we have combinations of V and H such as the vertical-horizontal (Vertical transmit and Horizontal received, VH) mode. Likewise we have VV, HV, and HH. It is thought that the HH mode is more suitable for distinguishing flood regions from non-flood [5,6]. The VV mode is more sensitive than the HH and cross-polarizations (VH, HV) to small-scale roughness of waves on the water surface. The VH component, which makes it easier for us to notice differences between a flood region and a non-flood, is more useful than the VV for detecting flood regions [7,8].

As to algorithms for land/water segmentation (river, lake, flood regions etc) in SAR imagery, thresholding-based approaches have been traditionally used [5]. Its underlying assumption is that the backscattering amplitudes which differ between the Earth’s surface types allows us to classify the points of imagery by setting thresholds in the amplitude histogram. For instance, the surface layer backscatter coefficents of high moistured soil tend to be higher than those of low moistured soil [9]; the amplitude of smooth surface such as calm water regions is low, while rough surface yields a higher amplitude. As a result, water regions are generally darker in SAR imagery [1], and thus a threshold for land regions can be set. At the same time, the thresholds depend on several conditions, that is, wind, incident angle, polarization, speckle noize, and therefore segmentation results are affected by these conditions.

There have been various nonneuralnet machine learning (ML) approaches used for several data source: random forest for Sentinel-1 C-band [10], K-means clustering for Sentinel-1A [11], Fuzzy C-Means clustering for Gaofen-3 and Sentinel-1 [12], Support Vector Machine (SVM) [13], Markov Random Field (MRF) for TerraSAR-X [14], XGBoost for Resourcesat-2 [15]. In addition, several statistical models have been used for Radarsat-2 [16], COSMO-Sky-Med [17], and Sentinel-1 [18].

Also, superpixel-based segmentation approaches have been applied for water mapping. In the superpixel framework, an image which consists of standard pixels is partitioned into groups according to algorithms, and each of the groups has statistically or morphologically meaningful information. In terms of image description, it could be argued, the superpixel representation maintains the same level of detail as the pixel representation, while the latter uses larger computational resources [19]. The framework threfore has potential capability to organize data from imagery faster and focus on more relevent regions.

Superpixel-based approaches have a range of applications: change detection [20], image segmentation [21], automatic target recognition (ATR), extracting water regions and coastlines [22], etc. For segmentation within the framework of superpixel, there has been a variety of superpixel methods especially for optical images; simple linear iterative clustering (SLIC), an edge-based SLIC (ESLIC) [23], pixel intensity and location similarity (PILS) [24], quick shift (QS) [25], and turbo pixels (TP) [26,27]. When it comes to SAR imagery, mixture-based superpixel (MISP) method is robust to SAR imegery’s speckle noise, which facilitates more practical use of the model [28].

Not only non-neuralnet ML models but neuralnet-based models have been applied for water mapping. CNN-based models especially U-Net models have shown good performances in typical evaluation metrics [29] - besides, they have been applied to near-real-time flood detection because of their speed [8].

In relation to encoder-decoder-type architecture, we have seen that the course of evolution of Seq2Seq [30] into Transformer was directed by the alignment model [31] (commonly known as the attention mechanism [32,33]). Compared to CNN and RNN, Transformer has smaller complexity per layer, minimum number of sequential operations, and maxium per lengths [33]. The attention-based models, which make it easier to find distant relations between words in sentences, have shown a brilliant performance in Natural Language Processing and accordingly been widely used in the field. The attention-based models are used in not just NLP but Computer Vision subject. Vision Transformer (ViT) converts an image into a sequence of patches which are used for embeddings, together with positional encoding, so that it can leverage the contextual information about the image [34].

The attention mechanism has been used for improving a neuralnet model’s human interpretability as well. For instance, real driving datasets were analyzed by constructing a visual attention model and computing attention maps to find which parts of the video image were cared by the model [35]. For chest radiograph classification, an attention-based model outperformed GradCAM (gradient-weighted class activation mapping) – it was concluded that, compared to GradCAM, the attention-based model could be more helpful for radiologists in terms of decision making [36]. Also, audio data were classified by an attention-based deep sequence model, where the attention scores were taken advantaged of to visualize which parts of the specrogram were cared [37].

The objective of the present study is to propose novel segmentation schemes whereby contextual information about superpixels is utilized for the input of machine learning to perform land/water segmentation of SAR imagery. This study also assesses our model’s performances in land/water classification, applying our models for Sentinel-1’s flood data [38].

The rest of this paper is organized as follows. Setcion 2 describes the satellite data and its data prosessing for the present study, and presents the method for generating and annotating superpixels. The section also provides schemes which are designed for unstructured/structured input for ML, showing performance indicators for land/water segmentation as well. Section 3 describes the setting of parameters; the section does feature selection for the structured scheme, presenting the results of numerical experiments. The section also makes a qualitative analysis by using attention rollout scores. Section 4 interprets the findings from Section 3. The final section is devoted to our concluding remarks.

2. Materials and Methods

2.1. Sen1Floods11 Dataset

Sen1floods11 is a dataset which processed the data of 11 flood events during 2016-2019, providing global coverage of 4831 chips of 512x512 pixels [38]. The flood imagery data of the dataset were from Sentinel-1 and Sentinel-2 – the former has VV and VH, and the latter has 13 bands, where the ground resolution is sampled to 10 meters on all bands. In the dataset, 4370 chips were automatically labeled; Otsu’s thresholding algorithm for the histogram from the VH band was used for the Sentinel-1 flood mapping reference, and thresholds of Normalized Difference Vegetation Index (NDVI) and Modified Normalized Difference Water Index (MNDWI) were used for the Sentinel-2. The Sen1floods11 has also 446 chips which were hand labeled by trained remote sening analysis, using a GUI of Google Earth Engine.

The hand labeled chips were treated as high quality data, and their test data were used for performance assessment [38]. The breakdown of the hand-labeled dataset is shown in Table 1, where the dataset contains 252 chips for train, 89 for validation, and 90 for test. Also 15 chips for Bolivia data is included in the hand-labeled dataset. Among them we used train, validation, and test chips.

Here it should be noted that we have 3 types of water area: permanent water (PW), temporaly water (TW), and all water (AW). PW was determined by using the JRC (European Commission Joint Research Centre) surface water dataset [39] which provides surface water observations at 50m resolusion on a monthly basis. The permanent water regions were given as the regions which were detected as water at both the start 1984 and the end 2018 of the dataset [38]. TW is an area temporalily covered by water, and AW is given by the union of TW and PW.

2.2. Data Processing

We now provide a brief description about the whole picture of processing the data. We use Sentinel-1 SAR imagery (chips) and hand-labeled data for training and assessing our models. The values of a chip is converted from dB scale into amplitude. The chip is then divided into four tiles, each of which is 256x256 pixels and used for generation of superpixels. Each superpixel is labeled a land/water superpixel acccording to percentages of land/water pixels in the superpixel. In the present unstructrued scheme, the generated superpixel data are further processed so that the data can be used as the input of our neuralnet model, where either our single or neighborhood framework is selected. Under the single framework, only a target superpixel is given in a training/testing process, while the neighbor provides both a target and its neighbor superpixels. The model is trained after the relating hyperparameters are tuned. Then each superpixel is classified by the model, which leads to segmentation of the chip.

When it comes to our structured scheme, features are extracted from each of the generated superpixels and their values (together with land/water class values) are stored as their corresponding columns in a file. Also our single/neighborhood scheme is adopted; the number of columns under the neighbor is larger than those under the single. Then hyperparameter tuning is performed and a model is trained together with feature selection in order to build an optimal model. The trained model is used for superpixel land/water classification for each superpixel, or equivalently, segmentation is performed. The flow diagram of the data processing is illustrated in Figure 1.

2.3. Superpixel Generation

Suppose a SAR image consists of N pixels and it is divided into K superpixels. Then we assume the conditional probability distribution of the amplitude $a_n$ of the n-th pixel given the k-th superpixel is modelled by Nalagami distribution [28]:

$$\begin{array}{c}\hfill p\left(a_n\right|{\mu }_k,{\nu }_k)={\displaystyle \frac{2}{\Gamma \left({\nu }_k\right)}}{\left({\displaystyle \frac{{\nu }_k}{{\mu }_k}}\right)}^{{\nu }_k}{\left(a_n\right)}^{2{\nu }_k-1}{e}^{-{\displaystyle \frac{{\nu }_k{a}_n^2}{{\mu }_k}}}\end{array}$$

(1)

where ${\mu }_k,{\nu }_k$ are the Nakagami parameters for the k-th superpixel, and $\Gamma \left(x\right)$ is the Gamma function. As to the conditional distribution of the position ${\mathbf{q}}_n={(x_n,y_n)}^T$ of the n-th pixel given the k-th superpixel, it is expressed by divariate normal distribution:

$$\begin{array}{c}\hfill p\left({\mathbf{q}}_n\right|{\mathbf{m}}_k,{\Sigma }_k)={\displaystyle \frac{1}{\left(2\pi \right)\sqrt{|{\Sigma }_k|}}}{e}^{-{\displaystyle \frac{1}{2}}{({\mathbf{q}}_n-{\mathbf{m}}_k)}^T{\Sigma }_k^{-1}({\mathbf{q}}_n-{\mathbf{m}}_k)}\end{array}$$

(2)

where the centroid and covariance are ${\mathbf{m}}_k$ and ${\Sigma }_k$, respectively. We introduce a handy parameter ${\theta }_k=\{{\mu }_k,{\nu }_k,{\mathbf{m}}_k,{\Sigma }_k\}$ for superpixels.

Assuming that the amplitude and position of a pixel is statistically independent, we introduce a vector ${\mathbf{f}}_n={(a_n,{\mathbf{q}}_n)}^T$ for the n-th pixel. After that, we introduce a K-dimensional label vector ${\mathbf{z}}_n\in {\mathbf{R}}^K$ defined by

$$\begin{array}{cc}\hfill & {\mathbf{z}}_n\in \{{(1,0,\dots ,0)}^T,{(0,1,0,\dots ,0)}^T,\dots ,{(0,\dots ,0,1)}^T\},\hfill \\ \hfill & (n=1,\dots ,N),\hfill \end{array}$$

(3)

where each vector is normalized as ${\sum }_{k=1}^K{z}_{n,k}=1$. This vector connects a pixel to a superpixel; for example, ${\mathbf{z}}_8=(0,1,0,\dots ,0)$ means that the 8th pixel belongs to the 2nd superpixel. The prior distribution for this vector is represented by a multinominal distribution

$$\begin{array}{c}\hfill p\left({\mathbf{z}}_n\right|{\omega }_{1:K})=\prod _{k=1}^K{\omega }_k^{z_{n,k}}\end{array}$$

(4)

where ${\omega }_1,\dots ,{\omega }_K$ are the parameters for the multinominal distribution. The prior of ${\omega }_{1:K}$ is

$$\begin{array}{c}\hfill p\left({\omega }_{1:K}\right)={\displaystyle \frac{\Gamma \left(\alpha K\right)}{{(\Gamma \left(\alpha \right))}^K}}\prod _{k=1}^K{\omega }_k^{\alpha -1},\end{array}$$

(5)

where $\alpha$ is the concentration parameter.

Then we optimize these parameters ${\left\{{\theta }_k\right\}}_{1,\dots ,K}$, ${\left\{{\omega }_k\right\}}_{1,\dots ,K}$, ${\left\{{\mathbf{z}}_n\right\}}_{1,\dots ,N}$ to obtain the posterior distribution. We deal with it by splitting the divariate distribution-related part and concentration-related: $p({\mathbf{z}}_{1:N},{\theta }_{1:K},{\omega }_{1:K}|{\mathbf{f}}_{1:N})\propto p\left({\mathbf{f}}_{1:N}\right|{\mathbf{z}}_{1:N},{\theta }_{1:K})p\left({\mathbf{z}}_{1:N}\right|{\omega }_{1:K})$. Given that $a_n$ and ${\mathbf{q}}_n$ are statistically independent, the former is futher splitted by the position-related part and the amplitude-related [28]: ${\displaystyle p\left({\mathbf{f}}_{1:N}\right|{\mathbf{z}}_{1:N},{\theta }_{1:K})=\prod _{n=1}^N\prod _{k=1}^K{\left[p\left(a_n\right|{\theta }_k)p\left({\mathbf{q}}_n\right|{\theta }_k)\right]}^{z_{n,k}}}$

2.3.1. Block ICM

As to maximum a posteriori (MAP), we have two assumptions: (1) each pixel in a likelihood function is independent (2) each pixel is Markovian. In this case, we can iteratively update the parameters to see the convergence for them. This method, called Iterated Conditional Modes (ICM) [40], allows us to obtain the local minima of the parameters at low computational cost. In the present study, we use a block ICM approach which is faster than conventional ICM, and parameters are updated in the following order [28]:

$$\begin{array}{c}\hfill {\mathbf{z}}_n^{\left(t\right)}=\underset{{\mathbf{z}}_n}{max}p\left({\mathbf{f}}_n\right|{\mathbf{z}}_n,{\theta }_{1:K}^{(t-1)})p\left({\mathbf{z}}_n\right|{\omega }_{1:K}^{(t-1)})\end{array}$$

(6)

$$\begin{array}{c}\hfill {\theta }_k^{\left(t\right)}=\underset{{\theta }_k}{max}p\left({\mathbf{f}}_{1:N}\right|{\mathbf{z}}_{1:N}^{\left(t\right)},{\theta }_k)\end{array}$$

(7)

$$\begin{array}{c}\hfill {\omega }_k^{\left(t\right)}=\underset{{\omega }_k}{max}p\left({\mathbf{z}}_n^{\left(t\right)}\right|{\omega }_{1:K})p\left({\omega }_{1:K}\right)\end{array}$$

(8)

where t means pseudo-time. That is, we first obtain ${\mathbf{z}}_n^{\left(t\right)}$, and accordingly we analytically determine $\left\{{\theta }_k\right\}$ and $\left\{{\omega }_k\right\}$. Specifically we have

$${\mathbf{m}}_k^{\left(t\right)}={\displaystyle \frac{1}{|S_k|}}\sum _{k\in S_k}{\mathbf{q}}_n$$

for the centroid parameter, and

$${\Sigma }_k^{\left(t\right)}={\displaystyle \frac{1}{|S_k|}}\sum _{n\in S_k}({\mathbf{q}}_n-{\mathbf{m}}_k^{\left(t\right)}){({\mathbf{q}}_n-{\mathbf{m}}_k^{\left(t\right)})}^T$$

for the covariance. As to Nakagami scale parameter, we have ${\mu }_k^{\left(t\right)}={\displaystyle \frac{1}{|S_k|}}{\sum }_{n\in S_k}{a}_n^2$. When it comes to the Nakagami shape parameter ${\nu }_k^{\left(t\right)}$, we have the derivative of the likelihood function:

$$log{\displaystyle \frac{{\nu }_k}{{\mu }_k^{(t-1)}}}-\psi \left({\nu }_k\right)+{\displaystyle \frac{2}{|S_k|}}\sum _{k\in S_k}log{a}_n$$

where $\psi$ is the digamma function. We then search for the root of it to update the shape parameter.

2.3.2. Superpixel Annotation

After generating superpixels from the given data, we label them land/water according to predefined strategies. That is, the superpixel of which land/water pixels compose more than 68 percent is given as a land/water superpixel.

2.4. Superpixel Classification

The generated superpixels are classified into two major categories, (a) water (b) land. To do this, we use supervized machine learning, considering two types of schemes: unstructured and structured data schemes. Under the unstructured scheme, we use shallow variant models (MISP-SDT) of DeiT-Ti [41], where DeiT is the abbreviation of Data-efficient image Transformer; we use boosting algorithms under the structured scheme.

2.4.1. Unstructured Data Scheme

We present our unstructured data model and the architecture. The generated superpixels are put into tiles; each of them has 72x72x3 pixels. These tiles are augmented through random flips and rotations. We divide a tile image $\mathbf{x}\in {\mathbf{R}}^{H\times W\times C}$ ($H=72,W=72,C=3$) into flattened patches ${\mathbf{x}}_{pat}\in {\mathbf{R}}^{N\times \left(P^{2}C\right)}$, where $N(=HW/P^2)$ is the number of patches [34] and $P=6$. The patch sequence and the class label (CLS) emmbedding ${\mathbf{x}}_{CLS}\in {\mathbf{R}}^N$ are put together, and accordingly the length of the sequence is $N+1$. Then they are embedded by $E\in {\mathbf{R}}^{\left(P^{2}C\right)\times D}$, and equipped with positional embedding $E_{pos}\in {\mathbf{R}}^{(N+1)\times D}$, where D is hidden dimension (latent vector size [34]). This yields the input for the transformer block:

$$\begin{array}{c}\hfill {\zeta }_0=[{\mathbf{x}}_{CLS}^T;{\mathbf{x}}_{pat}E]+E_{pos}.\end{array}$$

(9)

In the transformer block, the input enters layer normalization and multi-head self attention (MSA), and the output of the MSA is added to the input to form an intermediate object ${\tilde{\zeta }}_1$. The intermediate goes into layer normalization and subsequently a Multi-layer Perceptron (MLP) which consists of an nn layer with GELU activation and a dropout layer; the output of the MLP is added to the intermediate to obtain the input ${\zeta }_1$ for the next transformer block. Likewise we have

$$\begin{array}{cc}\hfill & {\tilde{\zeta }}_b=MSA\left(LayerNorm\left({\zeta }_{b-1}\right)\right)+{\zeta }_{b-1}(b=1,2,\dots ,B)\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & {\zeta }_b=MLP\left(LayerNorm\left({\tilde{\zeta }}_b\right)\right)+{\tilde{\zeta }}_b\hfill \end{array}$$

(11)

where B is the number of transformer blocks. In the present study, we set $B=8$. Then, the first row (CLS) of ${\zeta }_B$ enters the classification head which contains a layer normalization and a linear; land/water classification is performed.

$$\begin{array}{c}\hfill \mathbf{y}=CH\left({\left({\zeta }_B\right)}_{0,*}\right)\end{array}$$

(12)

where $CH$ is the classification head, and ${\left({\zeta }_B\right)}_{0,*}$ is the first row of ${\zeta }_B$.

Single/Mask/Neighborhood models

In the present study, we leverage the color channels R, G, B of the ViT input imagery ${\mathbf{x}}^{H\times W\times C}$ so that we can incorporate the neighborhood information of a target superpixel into the scheme. To see this, the amplitude data from a raw tif file, its contrast-stretched amplitude, and the information about whether a superpixel is a target one or neghborhood are assigned to $R,G$, and B, respectively. We then call R, G, and B, the SAR channel, the stretched SAR channel, and the role channel, respectively. In putting the role channel to use in the unstructured scheme, we adopt three approches for the input information about the MISP-SDT: (i) a target superpixel only (ii) a target and its neighborhood regions and masking the role channel (iii) a target and its neighborhood regions.

We call (i), (ii), and (iii) the single model, mask model and neighborhood model. In the single model, the role channel is set 0, and in the SAR and streched channels regions except for a target superpixel are set 0. This means that a target’s surrounding region is treated like a padding. In the mask model, the role channel is set 0. In the neighborhood model, the padding region, the neighborhood region, and the target superpixel region in the role channel are set 0, 0.5, 1.0, respectively. An example of the unstructured input for the single/mask/neighborhood models are shown in Figure 2, where the 1st, 2nd, and 3rd rows represent the inputs for the single, mask, and neighborhood schemes, respectively. In the figure, the SAR, stretched, and role channels for these schemes are shown as well.

Considering these things, we show the architecture of the MISP-SDT and an example of its input in Figure 3.

Attention rollout

We also provide a short explanation of MSA and attention score. Suppose $\mathbf{z}\in {\mathbf{R}}^{N\times D}$ is the input sequence, and the query, key, and value $(\mathbf{q},\mathbf{k},\mathbf{v}\in {\mathbf{R}}^{N\times D_h})$ are given by $[\mathbf{q},\mathbf{k},\mathbf{v}]=\mathbf{z}{U}_{qkv}$, where $U_{qkv}\in {\mathbf{R}}^{D\times 3D_h}$. Typically, we set $D_h=D/h$, where h is the number of heads. Then the self-attention is given as $SA\left(\mathbf{z}\right)=W\mathbf{v}$, where $W\in \mathbf{R}$ is an attention weight matrix. The weight matrix is explicitly expressed as

$$\begin{array}{c}\hfill W_{ij}=\mathrm{softmax}{(\mathbf{q}{\mathbf{k}}^T/\sqrt{D_h})}_{ij}={\displaystyle \frac{exp({\sum }_{m=1}^{D_h}{q}_{im}{k}_{mj}/\sqrt{D_h})}{\sum exp({\sum }_{m=1}^{D_h}{q}_{im}{k}_{mj}/\sqrt{D_h})}}\end{array}$$

(13)

Then MSA is given as

$$\begin{array}{c}\hfill MSA\left(\mathbf{z}\right)=[S{A}_1,S{A}_2,\dots ,S{A}_{h-1},S{A}_h]U_{MSA},\end{array}$$

(14)

where $U_{MSA}\in {\mathbf{R}}^{h{D}_h\times D}$. We then move on to attention rollout. Suppose $A^{\left(b\right)}$ is the attention score matrix for the b-th transformer block. Then, its raw attention is defined by $A^{\left(b\right)}=0.5(I+{\langle W_{att,h}\rangle }_h)$, where $W_{att,h}$ is the attention matrix of an h-head attention, and usually it is averaged over the heads. Then, given the number of transformer blocks is B, the attention rollout for the attention matrices is obtained by the following recursive calculation [42]:

$$\begin{array}{c}\hfill {\tilde{A}}^{\left(b\right)}=\left\{\begin{array}{cc}{A}^{\left(b\right)}{\tilde{A}}^{(b-1)}\hfill & (0<b\le B-1)\hfill \\ A^{\left(b\right)}\hfill & (b=0)\hfill \end{array}\right.\end{array}$$

(15)

2.4.2. Structured Data Scheme

We also consider turning unstructured superpixel data into structured data and using them as inputs for machine learning. It is expected that constructed data from neighborhood superpixels have many NaNs since the number of neighbor superpixels differs between target ones. On top of that, noise and outliers in amplitude could often appear in sar imagery. Also, there could be a case where the ground truth is not so accurate. Taking these things into account, we choose XGBoost [43], which handles sparsity to build tree structures and is robust to outliers, as the algorithm for this case.

Now let us provide a brief explanation about the XGBoost framework. Suppose we have a set of N data points $\mathcal{D}=\left\{({\mathbf{x}}_j,y_j)\right|{\mathbf{x}}_j\in {\mathbf{R}}^m,{\mathbf{y}}_j\in \mathbf{R},(j=1,\dots ,N)\}$, and a function $F^{*}$ maps $\mathbf{x}$ to y. Following the typical gradient boosting procedure, we approximate $F^{*}$ by the additive expansion of functions ${\left\{f_k\left(\mathbf{x}\right)\right\}}_{k=1,\dots ,t}$ (aka base learner [44]). In the XG tree boosting method, we have

$$\begin{array}{c}\hfill {\widehat{y}}_i^{\left(t\right)}={\widehat{y}}_i^{(t-1)}+f_t\left({\mathbf{x}}_i\right)=\sum _{k=1}^t{f}_k\left({\mathbf{x}}_i\right),\end{array}$$

(16)

after the t-th iteration that constructs a decision tree [43]. The shrinkage parameter $\eta \in [0,1]$, which each of the base learner functions is multiplied by, is also considered. This parameter beefs up the technical capabilities for training an XGB model to improve its performances, since it disperses the power of each base learner function over future trees [43]. The column subsampling technique is considered as well to run the program faster and address an overfitting issue.

The objective function of the t-th iteration is expressed by the sum of convex loss functions ${\displaystyle \sum _{i=1}^{n}l(y_i,{\widehat{y}}_i^{(t-1)}+f_t\left({\mathbf{x}}_i\right))}$ and a regularization term $\Omega \left(f_t\right)$. Expanding the former with regard to ${\widehat{y}}_i^{(t-1)}$ and retaining by the second-order term, we find

$$\begin{array}{cc}\hfill {\mathcal{L}}^{\left(t\right)}& \simeq \sum _{i=1}^n\left\{l(y_i,{\widehat{y}}_i^{(t-1)})+{\displaystyle \frac{1}{1!}}·{\displaystyle \frac{\partial l}{\partial {\widehat{y}}^{(t-1)}}}{f}_t\left(\mathbf{x}\right)+{\displaystyle \frac{1}{2!}}·{\displaystyle \frac{{\partial }^{2}l}{\partial {\left({\widehat{y}}^{(t-1)}\right)}^2}}{f}_t^2\left(\mathbf{x}\right)\right\}+\Omega \left(f_t\right)\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & =\sum _{i=1}^n\left\{l(y_i,{\widehat{y}}_i^{(t-1)})+{\displaystyle \frac{1}{1!}}{g}_i{f}_t\left(\mathbf{x}\right)+{\displaystyle \frac{1}{2!}}{h}_i{f}_t^2\left(\mathbf{x}\right)\right\}+\Omega \left(f_t\right)\hfill \end{array}$$

(18)

where $g_i=\left({\displaystyle \frac{\partial }{\partial {\widehat{y}}^{(t-1)}}}\right)l(y_i,{\widehat{y}}^{(t-1)}){|}_{{\widehat{y}}^{(t-1)}={\widehat{y}}_i^{(t-1)}}$, and $h_i={\left({\displaystyle \frac{\partial }{\partial {\widehat{y}}^{(t-1)}}}\right)}^{2}l(y_i,{\widehat{y}}^{(t-1)}){|}_{{\widehat{y}}^{(t-1)}={\widehat{y}}_i^{(t-1)}}$. In the XGBoost regime, the regularization term is given by ${\displaystyle \Omega \left(f\right)=\gamma \sum _{j=1}^T+\frac{1}{2}\lambda \sum _{j=1}^T{w}_j^2=\gamma T+\frac{1}{2}\lambda \sum _{j=1}^T{w}_j^2}$, where T is the number of the leaves of a tree. The objective function is therefore written as [43]

$$\begin{array}{c}\hfill {\mathcal{L}}^{\left(t\right)}=\sum _{i=1}^{n}l(y_i,{\widehat{y}}_i^{(t-1)})+\sum _{j=1}^T\left\{\left(\sum _{i\in I_j}{g}_i\right)w_j+{\displaystyle \frac{1}{2}}\left(\sum _{i\in I_j}{h}_i+\lambda \right)w_j^2\right\}+\gamma T,\end{array}$$

(19)

which indicates that the objective function is a convex quadratic function. Obviously it has the minimum at ${\displaystyle w_j={\widehat{w}}_j=-\frac{{\sum }_{i\in I_j}{g}_i}{{\sum }_{i\in I_j}{h}_i+\lambda }}$. Substituting this into equation (19) yields the analytic form of it:

$$\begin{array}{c}\hfill {\mathcal{L}}_{min}^{\left(t\right)}=\sum _{j=1}^T{\displaystyle \frac{-{\left({\sum }_{i\in I_j}{g}_i\right)}^2}{2({\sum }_{i\in I_j}{h}_i+\lambda )}}+\sum _{i=1}^{n}l(y_i,{\widehat{y}}_i^{(t-1)})+\gamma T.\end{array}$$

(20)

In building a tree according to Information Gain(IG), the objective function is utilized as the score function to calculate the IG for determining the splitting point of a node.

When it comes to the structured inputs to the XGB, we extract feature information from a superpixel and represent it as (the row of) a dataframe. For instance, a target superpixel provides us with the 50th percentile, standard deviation of the amplitude, the area of it, and they are named pc50_cent, std_cent, area_cent, respectively. A column is assigned to the values of a feature. Likewise, a neighbor superpixel provides feature columns. Also, we consider the angle between a basis axis and the line segment which is bounded by the x-y median point of a target superpixel and that of a neighbor superpixel. In the process of superpixel generation, a superpixel could become like a set of uncombined islands, and therefore we use the number of the islands as a feature as well. All of the features used in the MISP-XGB models are shown in Table A1 of the appendix.

Feature selection is an important part of building machine learning models, allowing us to train a model faster, have better interpretability, and extend models more easily. There exist several ways to perform feature selection: filter, wrapper, embedded, and hybrid. To some extent, feature selection is embedded into the XGBoost as the tree boosting algorithm automatically picks the sets of features in real terms during the program running – the algorithm uses the feature which maximizes the XGBoost’s information gain for each splitting of a tree. Having said that, it is more desirable to construct an optimal model by conducting feature selection explicitly, because it improves model performance.

Specifically for XGBoost, rank-based feature selection methods were developed [45,46]. In the present study, the feature importance score for each of the features is calculated according to the weight which counts how many times a feature is used for tree splitting. Then the features are arranged in descending order of the importance score. Provided that x is a threshold, the top x percent of the features are then used to train a model to get the validation score. For several different thresholds, the validation scores are calculated to find out an optimal model.

2.5. Performance Metrics

Now that we consider a land/water binary classification problem, we introduce 4 cases in terms of segmentation: True Positive (TP), True Negative (TN), False Positive (FP), and False Negative (FN). TP is the case where water is correctly predicted, and TN is the case where land is correctly predicted. FP is the case where water is wrongly predicted, and FN is the case where land is wrongly predicted. These terms are important especially when there is a significant imbalance between land and water. We threfore use typical segmentation performance indicators: Accuracy, which is defined by ${\displaystyle \frac{TP+TN}{TP+TN+FP+FN}}$, Recall (or Sensitivity) ${\displaystyle \frac{TP}{TP+FN}}$, Presicion ${\displaystyle \frac{TP}{TP+FP}}$, and the F-1 score (the harmonic mean of Recall and Precision). Intersection over Union (IoU) ${\displaystyle \frac{TP}{TP+FP+FN}}$, which is known as Critical Success Index (CSI) and commonly used in this sort of segmentation [29], is considered as well.

In the present study, we define IoU (at the chip level) and its average by [47]

$$\begin{array}{cc}\hfill & Io{U}_s={\displaystyle \frac{T{P}_s}{T{P}_s+F{P}_s+F{N}_s}}\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & mIoU={\displaystyle \frac{1}{S}}\sum _{s=1}^{S}Io{U}_s\hfill \end{array}$$

(22)

where S is the total number of the samples, and $T{P}_s,F{P}_s,F{N}_s$ mean the true positive, the false positive, and the false negative for the s-th sample, respectively. We also consider IoU at the pixel level, which is given by $IoU={\displaystyle \frac{{\sum }_{s=1}^{S}T{P}_s}{{\sum }_{s=1}^S(T{P}_s+F{P}_s+F{N}_s)}}$ [47].

3. Results

3.1. Parameter Tuning

Numerical experiments were conducted using hand-labeled train data, where the VV mode and VH were used. We basically genarated around 1444 superpixels for a chip, setting the concentration parameter $\alpha =1000000$ [28], and iteration parameter 25 for the block ICM.

In carrying out numerical experiments, we tuned relating parameters for the unstructured/structured models. For the unstructured models, the learning rate and weight decay were set $0.0005\times {\displaystyle \frac{BatchSize}{512}}$ and 0.05, respectively. AdamW was chosen as the optimizer [41].

For the structured models, the boosting count and the max depth of a tree were set 1450 and 10, respectively. The learning rate was set 0.0035, and early stopping round was 250. Also, the land-to-water ratio was set 11.319:1 for the weight in order to address the class imbalance issue. The parameter choices and values for the models are presented in Table 2.

It is noted that we excluded NaN region of SAR imagery from performance assessment as our schemes were designed to use SAR imagery, not optical imagery. Also, the numerical experiments were performed within the framework of Tensorflow (2.15.1) under the unstructured scheme, together with NVIDIA Tesla T4 GPU. Under the structured scheme, XGBoost (2.0.3) was used for the numerical experiments.

3.2. Feature Selection

Under the MISP-XGB schemes, we had 19 features and 342 features in the single and neighborhood models, respectively. Setting several different percentage thresholds for the features, we calculated the validation scores to determine an optimal set of the features.

Figure 4 shows how the validation scores converge in relationship to the top x percent of the whole features. The x-axes of the figures represent the top x percent of all features for the single models (upper subfigures) and the top x percent of all features for the neighborhood models (lower); the y-axes show the validation scores which are mIoU (left subfigures) and IoU (right). Each of the subfigures presents the results from the VV mode (red-colored) and the VH (blue).

It can be seen from the upper subfigures that the validation scores start to converge at around the top 50 percent of the single features for the single model in the VV mode. Similarly, convergence starts at around the top 50 percent for the single model in the VH mode. Therefore, the optimal threshold for the feature selection was set 50 percent; namely this amounted to top 9 features to be used for the single model in both VV and VH modes. Likewise, we estimated that the threshold is at around the top 20 percent of the whole features – or equivalently we obtained top 40 features and top 38 features in the VV mode and the VH, respectively, in the neighborhood models.

The features used under the MISP-XGB scheme are presented in Figure A1 in Appendix. In the figure, bar charts in the VV and VH modes are filled with red and blue, respectively – the results from the single (upper subfigure) and neighborhood (lower) models are shown together. From the figure, it can be seen that neighborhood-superpixel features accounted for 60 percent of the top 10 features in the VV mode. They accounted for 90 percent of the top 10 features in the VH modes, which suggests that neighborhood-superpixel features were more frequently used in the VH mode, compared to the VV.

3.3. Experimental Study

The results of the numerical experiments are presented in Table 3. We firstly look at the results from the unstructured scheme. From the table, it can be seen that the neighborhood model generally showed better performances than the single model in the segmentation performance indicators in both polariation modes, except that the recall value of the neighborhood model was worse than that of the single in the VV mode, and the recall value of the neighborhood model was comparable to that of the single model in the VH mode. For both VV and VH, the mask model underperformed the single and neighborhood models except for recall.

Then we go on to the structured scheme. In the VV mode, the neighborhood model generally had better scores than the single model except that the recall value of the neighborhood model was lower than that of the single model. This pattern emerged in the VH mode as well – in this sense the results from both unstructured and structured schemes were in line with each other.

A few of examples of the results from the MISP-SDT scheme are shown in Figure 5, where the SAR imagery (1st subfigure), the segmentation results from the single model (2nd), the neighborhood model (3rd), and the ground truth (4th) for Somalia_94102 (top) and USA_504150 (bottom) are seen. In the figure, land, water, and unknown regions are filled with green, blue, and grey respectively.

When it comes to polarization, the VH mode produced on the whole more encouraging results than the VV; for instance under the unstructured scheme the f1 score and mIoU in the single model were increased from 0.376, 0.292 (VV) to 0.431, 0.338 (VH), respectively. A previous study showed that the precision of the VV component is generally larger than that of the VH – the co-polarized mode got lower scores than the cross-polarized in recall [8]. The finding was in line with our results in terms of precision and recall. These results indicate that the cross-polarized components are more suitable for flood/non-flood segmentation. This is consistent with previous studies.

We then compare our models to others in terms of AW. Table 4 displays a table which draws a comparison between our models and others. In the table, FCNN (fully convolusional neural network) is a CNN-based model [38], where FCNN (HL) and FCNN (S1W) are the models trained on respectively the hand-labeled and Sentinel-1 Weak datasets. Also, AlbuNet-34 (AN-34) is a model architecture developed from the LinkNet which links the input of its encoder block to the output of its conresponding decoder block to retain spatial information [48]. AN-34’s encoder is based on pretrained ResNet-34 [49]. Using summation for the information flow from the encoder to the decoder, AN-34 has fewer parameters than U-Net to gain faster training and inference speed [50]. From the table, it can be seen that the MISP-SDT and MISP-XGB models outperformed these CNN-based models in IoU and mIoU.

3.4. Qualitative Analysis

Also, we consider how the existence of neighborhood superpixels makes a difference in the unstructured scheme, by utlizing the attention map, where the attention rollout scores relating to CLS label embedding are visualized.

Figure 6 shows several examples of sets of stretched SAR data channel (the 1st column), role channel (2nd), attention map for single (3rd), attention map for mask (4th), attention map for neighborhood (5th), and class label map (the ground truth) (6th). For the class label map, land, water, and unknown regions are represented by green color, blue, and black, respectively. From the subfigures, we see that high score regions (red-colored regions) are mainly overlapped with the target region and its surrounding regions in the single scheme – this suggests that the target and its surrounding regions are more closely looked at by the single model. In terms of the distribution of attention map scores, the mask model showed patterns which were similar to those of the single model, although there exist some exceptions such as a more duffused type or almost uniformly distributed type. Meanwhile, in the neighborhood model, high and reltively high score (light-green-colored regions) regions are spread to the neighborhood regions, suggesting that not only the target and its surrounding regions but the neghborhood regions are cared by the model. Having acknowledged this, it is also crucial to recognize that these are examples of the numerous results of attention rollout scores that were calcuated by the single, mask, and neighborhood models. We look at this more carefully in the discussion section.

4. Discussion

We have seen that the neighborhood model generally outperforms the single model in both VV and VH modes for the MISP-SDT/XGB schemes. In the VV mode, the neighborhood model had better scores than the single model except that the recall values (0.459 in the unstructured scheme and 0.533 in the structured) for the single models were higher than those (0.425 in the unstructured and 0.514 in the structured, respectively) for the neighborhood models. These differences in recall between the single model and the neighborhood were mitigated in the VH mode: the recall values were 0.542 and 0.580 in the single models, compared to those (0.543 and 0.570) in the neighborhood models. That is, they are comparable for the unstructured scheme, and the gap was narrowed for the structured scheme in the VH mode.

As to the polarization, it has been found that the VV models underperformed the VH for both MISP-SDT and MISP-XGB, regardless of whether the model is single or neighbohood. These findings were in line with those of previous studies, supporting the notion that cross-polarized SAR imagery is more suitable for flood detection than the co-polarized. It was also observed that under the MISP-XGB, the results from the neighborhood model showed 6 neighbor features were ranked in the top 10 features in the VV mode; 9 neighbor features were ranked in the top 10 in the VH mode. It could be argued that the VH mode made it easier for the algorithm to utilize neighbor features for the classification, compred to the VV.

We have carried out the qualitative analysis as well; we now introduce a concept of Near Neighbor region to evaluate the results from the qualitative analysis more carefully. Roughly, the Near Neighbor region of a target superpixel is defined to be the vicinity of the target superpixel. More specifically, we define the circle area whose center is the centroid of the target region of a superpixel and the radius is equal to 1.1 times the semi major axis of the target region. Then, we define the Near Neighbor region of the superpixel by the intersection of the inside region of the circle area and the outside region of the target superpixel.

Figure 7 illustrates a few examples of part of SAR imagery (left subfigures), its target/neighbor superpixels (center), and corresponding Target/Near Neighbor/Far Neighbor regions (right). In the right subfigures, Target, Near neighbor, and Far Neighbor regions are filled with white, indigo, and gray, respectively.

Then, we calculated the average of the attention rollout scores for each region (i.e. Padding/Far Neighbor/Near Negihbor/Target), using the single, mask, and neighborhood models for the test VH data.

Figure 8 shows boxplots of the averages of attention rollout scores for Pad (1st subfigure), Far Neighbor (2nd), Near Neighbor (3rd), and Target (4th) regions, where attention rollout scores were min-max normalized on a scale of 0 to 1. In the figure, each subfigure shows the results from the single (red), mask (green), and neighborhood (blue) models. From the figure, it is apparent that for the neighborhood model the attention rollout scores in Target, Near Neighbor, and Far Neighbor regions are the highest, second highest, and third highest, respectively. This pattern emerges from the single and mask models.

It is noticeable that the score for Target region was the highest in the neighborhood model among the four regions – likewise the scores were in the mask and single models. As to the comparison between the three models, the neighborhood, single, and mask models had the highest, 2nd highest, 3rd highest scores, respectively for Target region, and this pattern is similar to those for Far Neighbor and Near Neighbor regions. It could be argued that these findings are related to the fact that the neighborhood model is given more information than the single and mask, and that the mask is not given information about the role of a superpixel.

It is thought that the neighborhood model gains an advantage in terms of segmentation because the model is provided more information about SAR and its stretch, and the roles of regions. These things suggest that the role channel has positive effects on land/water segmentation, allowing the neighborhood model to get an advantage. These findings suggest that contextual information about a target superpixel and its neighbor ones has overall positive effects on classification performanes of the algorithm.

5. Conclusions

We have devised novel superpixel schemes which use the mixture-based superpixel framework for single-date SAR imagery, incorporating the unstructured/structured data of both target and neighbor superpixels into the input for the ML algorithms. Designing the unstructured and structured schemes to annotate and segment target superpixels, we have designed MISP-SDT and MISP-XGB schemes, respectively. We have also developed single/mask/neighborhood models and single/neighborhood models for the unstructured and structured schemes, respectively. Accordingly, our schemes have enabled us to take advantage of contexual information about SAR imagery to segment superpixels.

Our schemes were applied to Sentinel-1 SAR data to examine segmentation performances. The results from the numerical experiments demonstrated that the cross-polarized mode produced better results than the co-polarized, which were consistent with previous studies. Under our MISP-SDT/XGB schemes, where we did not reply on Sentinel-2 optical imagery, the neighborhood models outperformed the FCNN model in IoU and mIoU.

As to models, the neighborhood model as a whole gave better performances than the single; this pattern of numerical results emerged regardless of whether the scheme was unstructured or structured.

We have used attention maps and feature importance scores, demonstrating that neighbor regions were looked at or used by the ML algorithms in the neighborhood models. Our findings suggest that under the unstructured/structured schemes contextual information provided by neighbor superpixels helps the ML algorithms to improve land/water classification performance. As part of future work, we mention that the interpretability of ML could be enhanced by further developing this kind of schemes. Moreover, it will be worthwhile to consider how this sort of schemes could be applied to different research fields and topics.