An Adaptive Unmixing Method Based on Iterative Multi-Objective Optimization for Surface Water Fraction Mapping (IMOSWFM) Using Zhuhai-1 Hyperspectral Images

Abstract

Surface water fraction mapping is an essential preprocessing step for the subpixel mapping (SPM) of surface water, providing valuable prior knowledge about surface water distribution at the subpixel level. In recent years, spectral mixture analysis (SMA) has been extensively applied to estimate surface water fractions in multispectral images by decomposing each mixed pixel into endmembers and their corresponding fractions using linear or nonlinear spectral mixture models. However, challenges emerge when introducing existing surface water fraction mapping methods to hyperspectral images (HSIs) due to insufficient exploration of spectral information. Additionally, inaccurate extraction of endmembers can result in unsatisfactory water fraction estimations. To address these issues, this paper proposes an adaptive unmixing method based on iterative multi-objective optimization for surface water fraction mapping (IMOSWFM) using Zhuhai-1 HSIs. In IMOSWFM, a modified normalized difference water fraction index (MNDWFI) was developed to fully exploit the spectral information. Furthermore, an iterative unmixing framework was adopted to dynamically extract high-quality endmembers and estimate their corresponding water fractions. Experimental results on the Zhuhai-1 HSIs from three test sites around Nanyi Lake indicate that water fraction maps obtained by IMOSWFM are closest to the reference maps compared with the other three SMA-based surface water fraction estimation methods, with the highest overall accuracy (OA) of 91.74%, 93.12%, and 89.73% in terms of pure water extraction and the lowest root-mean-square errors (RMSE) of 0.2506, 0.2403, and 0.2265 in terms of water fraction estimation. This research provides a reference for adapting existing surface water fraction mapping methods to HSIs.

1. Introduction

Surface water, a vital component of land water storage, is crucial for maintaining diverse and complex social and ecological systems [1]. Climate change and increasing anthropogenic activities are causing deteriorating water quality, declining water levels, and shrinking surface water areas, which has widespread impacts on human production and daily activities [2,3,4,5]. Changes in surface water distribution affect both the terrestrial and global water cycles, leading to severe floods or droughts [6,7,8]. Therefore, precisely mapping and monitoring the spatial distribution of surface water is essential for water resource management, agricultural irrigation, ecological conservation, disaster alerting, as well as the sustainable development of society and the economy [9,10,11,12].

In recent years, the rapid advancement of remote sensing (RS) technology has made it an efficient tool for comprehensively monitoring surface water bodies, thanks to its macrographic, periodic, and cost-effective nature [13]. RS data with various temporal, spatial, and spectral resolutions have been widely used for surface water extraction. Data from the Moderate-Resolution Imaging Spectroradiometer (MODIS), Landsat series, and Sentinel-2 are particularly popular on account of their relatively high spatial or temporal resolutions and free access [14,15,16,17]. For example, Xiao et al. [16] proposed a remote sensing framework for detailed monitoring of surface water dynamics using Landsat series in which a gap-filling method was developed to restore water bodies obscured by clouds and an innovative index was designed to quantify the intermittency of surface water bodies. Li et al. [17] proposed a high spatiotemporal surface water mapping framework combining Sentinel-1 and Sentinel-2 satellite imagery with ancillary products to generate 10 m surface water maps at a temporal resolution of 15 days. However, these multispectral images (MSIs) contain only a limited number of discrete spectral bands. The insufficient spectral information restricts their capability to distinguish surface water bodies from other background features with similar spectral characteristics [18]. In the past few decades, the hyperspectral image (HSI), which contains hundreds of continuous spectral channels from visible light to infrared, has witnessed enticing prospects in extracting water body information [19]. Extensive research has been conducted on applying HSIs for surface water mapping [20,21,22,23]. Zhao et al. [20] proposed a PCA-NDWI model for surface water extraction in urban areas using uncrewed aerial vehicle (UAV) hyperspectral data. In this model, the traditional normalized difference water index (NDWI) was calculated with the principal components obtained by principal component analysis (PCA) of the green and near-infrared bands to exploit the abundant spectral information of HSIs. Xie et al. [21] developed a hyperspectral difference water index (HDWI) for water extraction by integrating the reflectance of the red and near-infrared bands and calculating the reflectance difference between these two spectral regions. Considering the difficulties of selecting optimal threshold parameters for the traditional spectral water index, Chen et al. [22] deeply analyzed the characteristics of surface water spectra collected from HSIs and constructed inequality constraints as well as physical magnitude constraints to identify water from urban scenes. Qin et al. [23] proposed an improved U-Net to extract small water bodies from Zhuhai-1 HSIs, where the structure of the traditional U-Net was deepened to learn the nuances of different spectral bands in consideration of the large number of bands in HSIs. Actually, the primary advantage of HSIs lies in their richness of spectral information. However, directly applying traditional surface water extraction methods to HSIs often leads to inadequate leverage of the abundant spectral information. Effectively exploiting the spectral information contained in HSIs is a central focus in hyperspectral surface water mapping.

Representative methods for extracting surface water bodies from MSIs or HSIs mainly include the following: (1) Single-band thresholding. This kind of method extracts water bodies by comparing the reflectance of a single band (usually the infrared band) to a predefined threshold [24,25]. (2) Multi-band water index. This index is constructed by performing numerical operations between different bands to enhance the spectral characteristics of water bodies, and a corresponding threshold needs to be selected to extract water pixels [26,27,28]. (3) Machine learning (ML) and deep learning (DL). Traditional ML methods, such as support vector machine (SVM), random forest (RF), and gradient boosted machines (GBM), have exhibited effective performance in water extraction [2,29,30]. Additionally, newly developed DL-based methods, especially convolutional neural networks (CNNs), have become a research focus and demonstrated more robust results than traditional ML algorithms due to their capability to capture high-level abstractions of the original data [31,32,33,34]. (4) Object-oriented classification. Unlike pixel-level extraction methods, object-oriented classification considers spatially adjacent pixels with similar spectral features as a single object, serving as the basic processing unit [35]. Traditional methods such as band thresholding or ML are then applied to detect surface water objects based on their spectral and geometric features [36,37].

Most surface water extraction methods perform binary classification, dividing pixels into water or non-water categories. However, due to limited spatial resolution and ground feature heterogeneity, many mixed pixels exist in remote sensing images [38]. Mixed pixels create uncertainty in binary classification-based methods, particularly at water–land transitions. To achieve more accurate water extraction, subpixel mapping (SPM) has been introduced [39,40]. SPM determines the spatial distribution of each land cover class within a mixed pixel by breaking the mixed pixel into subpixels and categorizing each subpixel based on coarse fraction maps [41,42,43]. However, the accuracy of SPM is significantly influenced by the quality of the estimated fraction maps [44], which serve as critical inputs for SPM of land cover types, including water bodies [45]. Thus, an accurate water fraction map is vital for subpixel-level surface water mapping.

Existing methods for water fraction estimation mainly involve soft classification-based approaches [46], regression-based approaches [16,47], and spectral mixture analysis (SMA) [48,49,50,51]. SMA is a physically based technique for solving the mixed pixel problem [52,53], which considers mixed pixels to be linear or nonlinear combinations of endmember signatures weighted by their correspondent abundance fractions [54]. It is widely employed in water fraction estimation, primarily due to its clear physical interpretation and its advantage of not necessitating a substantial number of training samples [49]. For instance, Xie et al. [48] proposed an automated subpixel surface water mapping (ASWM) method for Landsat 8 OLI imagery, where mixed water–land pixels were determined by a double threshold segmentation method based on the histogram of the NDWI. Water fractions of each mixed pixel were then estimated via a linear spectral mixture model (LSMM), with endmembers iteratively selected from a local window. Xiong et al. [49] introduced a subpixel surface water extraction (SSWE) method for Landsat 8 OLI data, using a new water index, the all bands water index (ABWI), to extract pure water pixels. Mixed water–land pixels were obtained by the morphological dilation of pure pixels, and water fractions were calculated by local multiple endmember spectral mixture analysis. Jiang et al. [50] improved surface water fraction mapping using Sentinel-2 images by combining morphological dilation with morphological erosion to determine mixed pixels, and solved spectral variability by locally generating multiple water and land endmembers for unmixing. Lv et al. [51] proposed a spatial–spectral subpixel surface water extraction method (SSEM), combining spatial and spectral attributes to refine and decompose mixed water–land pixels.

Despite extensive research on water fraction estimation, challenges remain, particularly when the application scenarios transfer from MSIs to HSIs. Conventional water indexes used for segmenting pure water, such as the normalized difference water index (NDWI) [26] and the modified normalized difference water index (MNDWI) [27], utilize only a few bands of RS images, missing vast quantities of spectral information in HSIs. Strategies such as band averaging or summing multiple bands to construct water indexes [20,49] fail to capture the intrinsic spectral characteristics of water bodies, which are manifested as nuanced reflectance fluctuations across different bands in HSIs. Furthermore, careful consideration should be given to the accuracy of endmember extraction. Widely used local endmember extraction methods [48,49,50,51] may incorrectly identify endmembers, particularly in two cases: (1) Pure water pixels are inaccurately classified. If pure water pixels around a mixed pixel are mistakenly classified to the land class, they may be extracted as land endmembers, and vice versa. (2) Few pure pixels surround a mixed pixel. In this case, neighboring mixed pixels may be extracted as land endmembers. In the above two cases, incorrectly classified pure pixels or neighboring mixed pixels may be extracted as endmembers, leading to inaccurate estimation of water fractions. Therefore, local endmember extraction is not effective in all scenarios, and optimal endmembers should be adaptively extracted for mixed pixels under complex local conditions.

To address these challenges, this research proposes an adaptive unmixing method based on iterative multi-objective optimization for surface water fraction mapping (IMOSWFM) using Zhuhai-1 HSI. The main contributions of this work are summarized from two aspects: (1) A modified normalized difference water fraction index (MNDWFI) tailored for hyperspectral surface water extraction is designed, which elaborately exploits the spectral information of each band in HSIs through an unmixing process, accurately reflecting water distribution details and facilitating fine extraction of water bodies. (2) An iterative water fraction estimation framework based on reconstruction error optimization is proposed, where a multi-objective optimization algorithm is utilized to adaptively find the global optimal endmember combination for mixed pixels from different regions, significantly improving water fraction estimation accuracy.

2. Materials and Methods

2.1. Study Area

Nanyi Lake, located in Xuancheng City, Anhui Province (Figure 1), was selected as the primary study area. It extends between longitudes 118°50′E and 119°3′E and latitudes 31°1′N and 31°10′N. As the largest natural freshwater lake in Southern Anhui, Nanyi Lake plays an important role in regulating the water level of the Shuiyang River, which is connected to the Yangtze River, and maintaining the stability of the regional ecosystem [55,56].

Three distinct areas around Nanyi Lake, each with diverse surface water types and ground features, were meticulously selected for the experiments. The first area lies in the southwest, where two slender rivers converge. Small water bodies are distributed throughout the eastern land region, and several narrow ditches exist in the west. Surface water bodies are characterized by their clarity in this area. The second area, located in the northeast, contains a wide river surrounded by lots of block-shaped artificial ponds. The turbidity of the river and ponds varies greatly, resulting in noticeable spectral variabilities among the pure water pixels. The third area, situated in the southeast, consists of large tracts of paddy fields divided into small blocks by narrow ridges, both horizontally and vertically. These narrow ridges have similar spectra to the pure water pixels around them, posing huge challenges for the precise extraction of pure water pixels and estimation of the water fractions. In addition, we carried out experiments on three other areas located in different parts of China, Area 4, Area 5, and Area 6, to verify the applicability of IMOSWFM to diverse regions, which also contain different types of surface water. The locations of all the study areas are illustrated in Figure 1. The divergent characteristics of the surface water environments in these six study areas render them excellent candidates for evaluating the proposed surface water fraction mapping method.

2.2. Zhuhai-1 OHS Data

The RS data used in this research were acquired by the Zhuhai-1 Orbita hyperspectral satellites (OHS), as displayed in Figure 1. The Zhuhai-1 hyperspectral satellites are the first commercial hyperspectral satellites launched and integrated into a network in China. Each pixel of the Zhuhai-1 OHS image is 10 m × 10 m in size and there are 32 bands ranging from 0.443 μm to 0.940 μm, with a spectral resolution of 0.25 μm [57]. The even distribution of spectral bands from visible light to near-infrared as well as the high spatial resolution render the Zhuhai-1 OHS image an ideal data source for distinguishing surface water bodies from other similar ground features, which require high discrimination capabilities in the spectral domain [23]. The details of the Zhuhai-1 OHS images are summarized in

Table 1. Summary of the Zhuhai-1 OHS and reference images used in the present study.

Study Area	Product ID	Size of OHS Image	Acquisition Date of OHS Image	Acquisition Date of Reference Image
Area 1	HEM2_20230227235680_0014_L1B_CMOS3		27 February 2023	27 January 2023
Area 2		300 × 300 pixels		9 January 2023
Area 3				9 January 2023
Area 4	HEM2_20230219224117_0011_L1B_CMOS3		19 February 2023	/
Area 5	HGM2_20230215234218_0004_L1B_CMOS2	500 × 500 pixels	15 February 2023	/
Area 6	HGM2_20230203235828_0008_L1B_CMOS2		3 February 2023	/

.

2.3. Reference Data

The reference water fraction maps used to evaluate the performance of IMOSWFM and other compared algorithms were derived from fine spatial resolution (<1 m) RGB images downloaded from the GEOVIS Earth DataCloud “https://datacloud.geovisearth.com/ (accessed on 11 May 2024)”, acquired on 27 January 2023 (Area 1) and 9 January 2023 (Area 2 and 3). The reference surface water boundaries were manually digitized from the GEOVIS Earth (Beijing, China) images through visual interpretation. The Zhuhai-1 OHS images were co-registered with the GEOVIS Earth images using a second order polynomial. For each area, we selected at least 15 ground control points (GCPs), resulting in root-mean-square errors (RMSEs) of 0.18, 0.17, and 0.21 pixels, respectively. The reference water fraction maps were finally generated by calculating the proportion of surface water area in each OHS pixel, with the same spatial resolution as the OHS image.

2.4. Framework of IMOSWFM

The preprocessing of the Zhuhai-1 OHS image, which consists of layer stacking, radiometric calibration, and 6S atmospheric correction [58], is automatically carried out by the “OHS data processing tool” provided by Orbita. At the beginning of IMOSWFM, the improved multi-objective discrete particle swarm optimization algorithm (IMODPSO) [59] is utilized to extract endmembers from the original image. Three classes of endmembers are extracted: water, vegetation, and impervious surface. Then the fraction of water and other endmembers are estimated by the fully constrained least squares (FCLS) method.

Based on the estimated fraction maps, the MNDWFI is calculated to distinguish pure water, mixed water–land, and pure land pixels through a double threshold segmentation method. In the main loop, the water fractions of the mixed pixels are adaptively estimated through an iterative unmixing process. Firstly, endmembers are specifically extracted by IMODPSO. Secondly, the mixed pixels are unmixed by FCLS, and the reconstruction error of each mixed pixel is calculated. The estimated water fractions are assigned to the corresponding mixed pixels with a reconstruction error less than a predefined threshold r. The remaining mixed pixels, for which the water fractions are not assigned, will be unmixed in the next iteration.

IMOSWFM terminates when the number of mixed pixels obtaining water fractions is less than 1000 in two successive iterations or when the number of remaining mixed pixels to be unmixed is less than 5% of the total number of mixed pixels. Once the termination condition is met, endmembers are finally extracted by IMODPSO for all remaining mixed pixels, and their water fractions are estimated.

The framework of IMOSWFM is presented in Figure 2. In the following sections, we will describe the details of IMOSWFM.

2.5. Related Work

2.5.1. IMODPSO

IMODPSO simultaneously optimizes two objective functions: the volume inverse and the reconstruction error, which is the RMSE between the original image and its reconstructed image [59]. The volume inverse objective function is based on convex geometry theory, which assumes that endmembers are located at the vertices of the data simplex, forming a larger volume than that of mixed pixels [60]. Therefore, a larger volume indicates a higher likelihood that the pixels are endmembers. Thus, in IMODPSO, the volume inverse is employed as one objective function. In addition, the RMSE objective function is utilized, following the approach of many evolutionary algorithm-based endmember extraction methods [61,62]. A smaller RMSE indicates a better result of endmember extraction. IMODPSO aims to simultaneously minimize both objectives, which can be expressed as follows:

$$\begin{array}{c}\min :F\left(\mathbf{x}\right)=\left(f_1\right(\mathbf{x}),f_2(\mathbf{x}\left)\right)\\ f_1(\mathbf{x})={\displaystyle \frac{1}{\mathrm{volume}\left(\mathbf{A}\right(\mathbf{x}\left)\right)}}={\displaystyle \frac{(P-1)!}{\left|\det \left[\begin{array}{cccc}1& 1& \cdots & 1\\ {\mathbf{a}}_1& {\mathbf{a}}_2& \cdots & {\mathbf{a}}_P\end{array}\right]\right|}}\\ f_2(\mathbf{x})=\mathrm{RMSE}(\mathbf{Y},\hat{\mathbf{Y}})={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}\sqrt{{\displaystyle \frac{1}{L}}{‖{\mathbf{y}}_i-{\hat{\mathbf{y}}}_i‖}^2}\end{array}$$

(1)

where x denotes an individual of the swarm, representing a combination of endmembers, P denotes the number of endmembers, L and N denote the number of bands and the number of pixels in the image, respectively. A(x) is the endmember matrix and a_i (i = 1, 2, …, P) is a vector with P-1 dimensions obtained from the original endmember spectral vector by minimum noise fraction transformation [63]. Y = (y₁, y₂, …, y_N) and $\hat{\mathbf{Y}}=({\hat{\mathbf{y}}}_1,{\hat{\mathbf{y}}}_2,...,{\hat{\mathbf{y}}}_N)$ represent the original image matrix and the reconstructed image matrix, respectively. The reconstructed image $\hat{\mathbf{Y}}$ is obtained by the linear spectral mixture model (LSMM) using the extracted endmember spectra and the estimated abundance maps produced by the unconstrained least squares (UCLS).

The reason for optimizing both the volume inverse and the RMSE instead of only one of them lies in the inherent contradiction between the two objective functions. Existing research [64] has demonstrated that applying the volume inverse is in favor of capturing rare endmembers, but the results are easily influenced by outliers and noises in the real image. On the contrary, the RMSE is robust against interference from outliers and noises but possibly neglects some prospective rare endmembers. A tough dilemma in extracting optimal endmembers arises from the phenomenon that one certain endmember combination cannot achieve best values for both objectives concurrently [65]. Specifically, the improvement of one objective function is usually accompanied by the deterioration of the other, rendering it difficult to determine which objective function to prioritize. In multi-objective optimization, if m objective functions are for minimization, solution x is said to dominate y if and only if f_i(x) ≤ f_i(y) for i = 1, …, m and f_j(x) < f_j(y) for at least one objective function j [66]. A solution is called a Pareto optimal or nondominated solution if no other solution dominates it [67]. The goal of multi-objective optimization is identifying a series of nondominated solutions. By applying multi-objective optimization methodologies, a better trade-off can be achieved between the two objective functions, leading to more robust and accurate endmember extraction. This, in turn, greatly benefits the unmixing of mixed water–land pixels in the subsequent process.

The IMODPSO algorithm, proposed by Tong et al. [59], enhances the performance of the original multi-objective discrete particle swarm optimization (MODPSO) [64] algorithm by eliminating the predefined random selection probability as well as adopting an archive strategy. Two user-defined parameters, namely, the number of endmembers and the maximum number of iterations, need to be input to IMODPSO. At the initializing stage, IMODPSO randomly initializes the position of each particle. The global best archive (GBA) is initialized with the nondominated solutions of the initial population. In the main loop of IMODPSO, firstly, the best local guide gbest of each particle is chosen from the GBA by utilizing the Sigma method [68]. Secondly, the velocity and position of each particle are either updated based on self and social experience or randomly updated. Then, the pbest of each particle is updated, and these eliminated pbest are preserved in an archive. Subsequently, the GBA is updated by selecting all the nondominated solutions of the combination set (pbest ∪ A). At last, if the predefined maximum number of iterations is achieved, the algorithm will terminate and output the GBA; otherwise, the process of IMODPOS returns to the stage of choosing the gbest for each particle.

2.5.2. Selection of the Optimal Endmember Combination

Given that the output of IMODPSO is a set of nondominated solutions, and each solution corresponds to a combination of endmembers, it is crucial to reasonably select one to achieve robust water fraction estimation results. Due to fact that no single nondominated solution can simultaneously have better values than all others for both objective functions, determining the optimal solution is challenging. To find the solution that best trades off the two objective functions, we first adopt the min-max normalization to map the raw values of the two objective functions to the range of [0, 1] as listed below:

$${f_1}^{\prime }\left(\mathbf{x}\right)={\displaystyle \frac{f_1\left(\mathbf{x}\right)-f_{1\min }}{f_{1\max }-f_{1\min }}}, {f_2}^{\prime }\left(\mathbf{x}\right)=\frac{f_2\left(\mathbf{x}\right)-f_{2\min }}{f_{2\max }-f_{2\min }}$$

(2)

where f_1min and f_1max denote the minimum and maximum values of the first objective function among the solutions in GBA, respectively. The same definition applies to f_2min and f_2max. Then, the solution with the minimum sum of the two normalized objective functions is selected as the optimal solution as follows:

$${\mathbf{x}}_{optimal}=\underset{{\mathbf{x}}_i}{\mathrm{argmin}}{f_1}^{\prime }\left({\mathbf{x}}_i\right)+{f_2}^{\prime }\left({\mathbf{x}}_i\right)$$

(3)

To describe the strategy of selecting the optimal solution more clearly, the sketch map of the distribution of the selected optimal solution in the objective space is delineated in Figure 3. It can be seen that the selected solution is located at the bottom left of the objective space, which has relatively small values on both the volume inverse and RMSE compared with other unselected trivial solutions.

2.6. Adaptive Unmixing Framework for Surface Water Fraction Mapping

2.6.1. Spectral Characteristics of Surface Water

The output of IMODPSO is a combination of endmember spectra consisting of water, vegetation, and impervious surface. Before estimating the water fraction, it is necessary to accurately distinguish the water endmember from other land classes. To ascertain spectral signatures of typical ground features in Zhuhai-1 OHS images, we manually selected a set of water, vegetation, and impervious surface pixels from the image covering Nanyi Lake. The corresponding spectra are displayed in Figure 4. Notice that the last four bands (bands 29 to 32) were removed due to noise. As shown in Figure 4, the reflectance of water is significantly lower than that of vegetation and impervious surface in the infrared bands, consistent with previous studies [26]. Therefore, the spectrum with the lowest average reflectance in the infrared bands (from bands 20 to 28) is considered the water spectrum.

2.6.2. Extraction of Pure and Mixed Pixels

The abundant spectral information recorded by HSIs facilitates the precise identification of diverse ground features. However, most existing research on subpixel surface water mapping uses multiple-band spectral water indices to distinguish pure and mixed water pixels, relying on only a few bands of the image [45,48]. Cai et al. [69] proposed a normalized difference water fraction index (NDWFI) for the spatiotemporal mapping of surface water based on Landsat images, defined as follows:

$$\mathrm{NDWFI}={\displaystyle \frac{A_{\mathrm{Dark}}-A_{\mathrm{Others}}}{A_{\mathrm{Dark}}+A_{\mathrm{Others}}}}$$

(4)

where A_Dark denotes the fractional abundance of dark endmembers with a reflectance of 0 in all bands, and A_Others denotes the sum of the fractional abundance of other land endmembers. Both A_Dark and A_Others are obtained by the fully constrained least squares (FCLS). Actually, NDWFI leverages the low reflectance of surface water to estimate water fraction by using dark endmembers. However, the water spectra in Figure 4 showcase non-negligible fluctuations in reflectance from visible to near-infrared bands, significantly differing from the spectra of dark endmembers. This unveils the intrinsic reflection characteristics of water bodies, making it challenging for the NDWFI of a pure water pixel to reach 1, since the abundance of dark endmembers estimated by FCLS cannot approach 1, even if no other ground objects are present. Under this condition, the NDWFI values of pixels in one image may distribute in a narrow range, complicating the segmentation of pure and mixed pixels. Actually, in real scenes, surface water spectra are influenced by various factors such as illumination, atmospheric condition, terrain effect, and their intrinsic characteristics [70]. Spectra of surface water pixels from different RS images or different regions within a single image can vary significantly. For these water pixels with evident spectral variability, the estimated abundance of dark endmembers can also vary greatly, leading to fluctuated NDWFI values across diverse water bodies. For instance, the NDWFI of some clean water pixels with spectra similar to those of dark pixels can reach 1, while some turbid water pixels may have a lower NDWFI, such as 0.8. The varying distribution of NDWFI values among different water pixels poses great challenges to finding the appropriate threshold without a priori knowledge. For this reason, the NDWFI may exhibit uncertainty in extracting surface water when applied to RS data under complex imaging conditions.

To address this, we propose a modified normalized difference water fraction index named the MNDWFI to fully exploit the spectral information of HSIs and enhance surface water features. The MNDWFI is defined as follows:

$$\mathrm{MNDWFI}={\displaystyle \frac{A_{\mathrm{Water}}-A_{\mathrm{Others}}}{A_{\mathrm{Water}}+A_{\mathrm{Others}}}}$$

(5)

where A_Water denotes the fractional abundance of water endmembers, adaptively estimated with endmembers extracted from the image. The MNDWFI better reflects the water content since the dark endmember is replaced by the water endmember, which is dynamically extracted and represents the common spectral signatures of water bodies in one certain image. Therefore, the proposed MNDWFI stably distributes between −1 and 1, with larger values indicating a higher water fraction. Based on the MNDWFI, pixels can be divided into pure water, mixed water–land, and pure land pixels using the double threshold method as shown in Figure 5. The land threshold t₁ is automatically determined by the Otsu segmentation algorithm [71], while the water threshold t₂ is manually selected based on the maximum local change rate at the junction between the valley and the right peak of the histogram. The MNDWFI of pure water pixels falls within the range of (t₂, 1], and the MNDWFI of mixed water–land pixels falls within the range of [t₁, t₂].

2.6.3. Iterative Water Fraction Estimation for the Mixed Pixels

IMOSWFM applies an iterative framework to adaptively estimate the water fraction of mixed pixels in different regions. The pivotal motivation behind this iterative framework is to address the problem of inaccurate endmembers extraction arising from local endmember extraction methods. The local endmember extraction methods assume that mixed water–land pixels are usually located in the transitions from water to land. Pure land or water pixels (i.e., endmembers) lie around these mixed water–land pixels. However, large errors exist in the linear SMA (LSMA) results, partly because some endmembers are misclassified during the process of endmember selection [45]. If the segmentation results are inaccurate, land pixels may be mistakenly extracted as water endmembers. Similarly, water pixels may be extracted as land endmembers. Mixed water–land pixels classified as pure pixels also have the chance to be extracted as endmembers. These misclassifications can cause certain errors in the LSMA results. The constituents of land pixels are complicated. A pure land pixel may contain diverse ground features. Thus, these pure land pixels are not entirely pure from a more detailed perspective. Land endmembers extracted from these mixed land pixels can have negative impacts on the abundance estimation results.

Therefore, optimal endmembers may not be successfully extracted from only a local extent. Expanding the searching space can sometimes help to find the potential optimal endmember combination. Meanwhile, endmembers should not be extracted under the strict restrictions of existing classification results, since errors are inevitably generated during the classification process. To solve these issues, an iterative unmixing framework is adopted in IMOSWFM to automatically find the global optimal endmembers from the whole image for each mixed pixel by minimizing the volume inverse and the reconstruction error, leveraging the powerful global search capability of multi-objective optimization algorithms, which helps improve the water fraction estimation results.

As displayed in Figure 6, the process begins with the IMODPSO algorithm extracting an initial set of endmembers, which consists of water, vegetation, and impervious surface. Abundance maps of these endmembers are then estimated by FCLS. Next, the MNDWFI map of the original image is calculated by Equation (5). The pixels of the whole image are classified into pure water pixels, mixed water–land pixels, and pure land pixels using the double threshold method described in Section 2.6.2. The water fractions for pure water pixels and pure land pixels are set to 1 and 0, respectively, and these pixels are subsequently removed from the original image.

Following this, endmembers are extracted specifically for the mixed water–land pixels by IMODPSO. In this step, the objective function RMSE, as defined in Equation (1), is calculated only from the mixed pixels. Abundances for the extracted endmembers of these pixels are estimated by FCLS, and the spectra of reconstructed pixels are generated by the LSMM. Pixels with an RMSE value less than a user-defined threshold r are considered accurately unmixed, and their water fractions are determined.

Once the water fractions of the mixed pixels are determined, these pixels are removed from the current image and do not participate in the next iteration of water fractions. Note that the searching space of IMODPSO remains the whole image, while the calculation of the objective function RMSE adaptively changes based on the remaining mixed pixels. With each iteration, the number of remaining mixed pixels gradually decreases.

When less than 1000 pixels are removed in two successive iterations, or the number of remaining pixels is less than 5% of the total number of the mixed pixels, all remaining pixels will be unmixed in the final iteration regardless of their RMSE values. This approach ensures the completion of the water fraction map.

2.6.4. Restart Mechanism

Due to the stochastic property of particle swarm optimization, IMODPSO may occasionally fail to extract the endmember for ground features in some iterations. For example, it might extract two vegetation spectra and one impervious surface spectrum, or two water spectra and one vegetation spectrum. This can lead to inaccurate water fraction estimation, especially if the water spectrum is not extracted or if multiple water spectra are extracted.

To tackle this problem, a restart mechanism is adopted to avoid producing invalid endmember combinations. At each iteration, the endmember spectrum extracted by IMODPSO with the lowest average reflectance in the infrared bands is initially considered the water spectrum, as depicted in Section 2.6.1. The remaining two spectra are classified as land spectra. Then, the NDWI of the three endmembers are calculated. If the NDWI of the water spectrum is negative, or if one of the land spectra has a positive NDWI, indicating either no successful extraction of water spectrum or the extraction of more than one water spectrum, this endmember combination is deemed invalid, and IMODPSO is executed again.

If IMODPSO fails to obtain a valid endmember combination after three independent runs in a single iteration, the endmembers will be inherited from the last iteration. Specifically, if the final water spectrum has a negative NDWI, it will be replaced by the water spectrum extracted in the last iteration. If one of the land spectra has a positive NDWI, both of the two land spectra will be replaced by the land spectra extracted in the last iteration.

2.7. Performance Metrics

Quantitative evaluations were carried out on the accuracy of both the classification of pure water and the estimation of water fractions. The overall accuracy (OA) and kappa coefficient were calculated to assess the classification results. To calculate OA and kappa coefficient, the estimated and reference water fraction maps were transformed into binary images in which values of pure water pixels remained 1 and values of the pure land and mixed pixels were set to 0. Further, the root-mean-square error (RMSE) and the systematic error (SE) were adopted to quantify the accuracy of estimated water fractions. The two evaluation metrics are listed as follows:

$$\begin{array}{c}\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{N}}\left({\displaystyle \sum _{i=1}^N}{\left(f_i\left(ref\right)-f_i\left(esti\right)\right)}^2\right)}\\ \mathrm{SE}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}\left(f_i\left(ref\right)-f_i\left(esti\right)\right)\end{array}$$

(6)

where f_i(ref) and f_i(esti) denote the reference water fraction and the estimated water fraction of the ith pixel, respectively. N denotes the number of pixels in the image.

3. Results

3.1. Comparison Algorithms and Parameter Settings

Several SMA-based algorithms for surface water fraction estimation, namely the automatic subpixel water mapping method (ASWM) [48], subpixel surface water mapping method based on morphological dilation and erosion (DESWM) [50], and the spatial–spectral extraction method (SSEM) [51], are chosen as methods to compare to validate the effectiveness of IMOSWFM. All of the compared algorithms utilize the traditional multi-band water index, exploiting only a few bands to extract pure water, and most of them extract local endmembers to unmix mixed water–land pixels. For IMOSWFM, the user-defined RMSE threshold r is set to 0.01 while the maximum number of iterations for IMODPSO is configured to 100 for all areas.

3.2. Qualitative Evaluation of Accuracy

The surface water fraction maps obtained by IMOSWFM and other compared methods are displayed in Figure 7. It can be discerned that the water fraction maps obtained by IMOSWFM were closest to the reference maps both in overall representation and detailed accuracy. Firstly, IMOSWFM achieved more accurate identification of pure and mixed water pixels. The three compared methods mistakenly classified a large number of mixed or pure land pixels located at the edge of water bodies as pure water, especially in Area 2 and Area 3. This discrepancy stems from the limited ability of the NDWI to capture the spectral characteristics of surface water by only using the green and infrared band, since ASWM, DESEW, and SSEM extract pure and mixed pixels based on the histograms of the NDWI. In contrast, the MNDWFI fully exploits the abundant spectral information provided by HSIs through an unmixing process. This theoretical advantage enables IMOSWFM to accomplish various surface water extraction tasks that demand high spectral discrimination capabilities. Secondly, it is noteworthy that a considerable amount of small water bodies was not detected by DESWE and SSEM in Area 1 and Area 2 compared with ASWM and IMOSWFM, indicating that the automatic threshold selection methods adopted by DESWE and SSEM are not sufficiently adaptable to different scenes. From the zoomed-in figures of three selected regions in Figure 7, IMOSWFM captured the detailed variations in surface water fractions better than other algorithms. Ridges between adjacent paddy fields were accurately recognized and manifested as mixed pixels with low water fractions by IMOSWFM, particularly in Region 2 and Region 3. The compared methods failed to identify such nuanced characteristics. The water fraction maps estimated by IMOSWFM well matched the reference maps, clearly showing the transitions from water to non-water regions. DESWE excessively blurred the boundaries between water and non-water pixels because DESWE utilized both morphological dilation and erosion to detect the mixed pixels (Region 1 and Region 2). Consequently, pure water or non-water pixels near the boundaries may have been erroneously divided into mixed pixels. SSEM performed the worst among all the methods. ASWM outperformed DESWE and SSEM. However, the results of ASWM were inferior to IMOSWFM due to the incorrect extraction of pure and mixed pixels.

The surface water fraction maps of Area 4, Area 5, and Area 6 located in different parts of China estimated by IMOSWFM and the three comparison methods with the Zhuhai-1 OHS images are illustrated in Figure 8. Note that no reference water fraction map was generated for the three areas. It can be seen from Figure 8 that IMOSWFM outperformed other methods visually both in pure water extraction and water fraction estimation. In Area 4, the three compared methods erroneously extracted large amounts of impervious surface pixels as surface water pixels. Conversely, some water pixels within the paddy fields were not extracted, resulting in numerous cavities inside the water bodies. In comparison, IMOSWFM was able to extract the paddy fields completely and distinguish the boundaries between adjacent water blocks as shown in Region 4. Also, few urban pixels were classified as water pixels by IMOSWFM. In terms of Area 5, the main streams can be accurately extracted by all the methods. However, the three compared methods failed to extract small water bodies in urban regions. Specifically, IMOSWFM successfully extracted the narrow tributary located in the west of Area 5, which was disregarded by other methods. It is noteworthy that some shadows were extracted as water bodies by IMOSWFM, particularly in the low-albedo urban regions situated in the northwest of this area. The water fraction maps estimated by the four methods of Area 6 were similar. However, from the zoomed-in figures of Region 6, we can observe that IMOSWFM best delineated the geometric shapes of the paddy fields with clearly visible boundaries. In general, benefiting from the capability of the MNDWFI to fully exploit the spectral information provided by HSIs, IMOSWFM can well capture the details of water distribution and obtain robust surface water mapping results when applied to different HSIs.

To further elucidate why ASWM, DESWE, and SSEM produced poor water fraction mapping results in Area 4, we compared the NDWI (Figure 9b) and MNDWFI (Figure 9c) maps of Area 4 and analyzed the spectra of four representative pixels (Figure 9d). As displayed in Figure 9b, the NDWI map of Area 4 shows noticeable contrast between the paddy fields and their surrounding cropland. However, there are fewer differences between the NDWI values of the paddy fields and the urban areas, making it a challenge to appropriately distinguish between the two land types through a global NDWI threshold. Compared with the NDWI, the MNDWFI greatly enhances the contrast between surface water and other kinds of ground features (Figure 9c). Surface water bodies appear as bright tones in the MNDWFI map, with their internal boundary details clearly visible. Figure 9d presents the spectra of four typical pixels in which P1 and P4 are impervious and located in the northern and southeastern parts of Area 4, respectively. P2 and P3 are both surface water pixels. Notably, P2 is located in a slender river, which was accurately classified as pure water by the three traditional methods, while P3 is located within the paddy fields at the center of Area 4, which was misclassified as land by the three traditional methods. From Figure 9d, spectra of the two water pixels, P2 and P3, are not consistent with traditional water spectra. Their reflectance of NIR (band 20) is unexpectedly higher than that of the green band (band 10), resulting in negative NDWI values, similar to those observed for P1 and P4. P3 even has a smaller NDWI than P1. As a result, it was not extracted by the three traditional surface water fraction mapping methods. The spectral anomalies of surface water in Area 4 possibly stem from the inappropriate atmospheric correction [19] or the eutrophication of water bodies [22], which can lead to underestimation of reflectance in the visible bands or an increase in reflectance of the NIR bands. Despite the presence of spectral anomalies, it is obvious that water and impervious surfaces each possess distinguishable spectral characteristics in terms of both spectral curve shape and physical reflectance magnitude. Nevertheless, the NDWI, which is calculated from only two bands, cannot fully make use of this useful spectral information. In contrast, the MNDWFI takes each band of HSIs into consideration through spectral unmixing and adaptively extracts water endmembers from HSIs with varying spectral features, which demonstrates the more robust capability of hyperspectral surface water extraction, regardless of spectral variability.

3.3. Quantitative Evaluation of Accuracy

Table 2. OA and kappa coefficient of IMOSWFM and the compared methods in Area 1, Area 2, and Area 3 with the Zhuhai-1 OHS images.

	Overall Accuracy			Kappa Coefficient
	Area 1	Area 2	Area 3	Area 1	Area 2	Area 3
ASWM	89.03%	91.49%	83.16%	0.6163	0.7741	0.6579
DESWE	89.28%	91.22%	82.90%	0.5931	0.7572	0.6438
SSEM	88.55%	90.74%	79.49%	0.5968	0.7450	0.5970
IMOSWFM	91.74%	93.12%	89.73%	0.6026	0.8080	0.7632

The bold entities represent the optimal values.

displays the OA and kappa coefficient of IMOSWFM and the compared methods. IMOSWFM obtained the highest overall accuracy in each area, indicating the effectiveness of the MNDWFI in extracting pure water pixels. In Area 1 and Area 2, IMOSWFM had higher OA values, mainly due to the precise extraction of small water bodies which widely existed in the two areas. The greatest enhancement in overall accuracy by IMOSWFM was observed in Area 3, surpassing the second ranked ASWM by 6.57%. The reason is that the three compared methods failed to extract the narrow ridges distributed horizontally and vertically between the large tracts of paddy fields in Area 3. Most of the pixels located at these ridges were mistakenly classified as pure water by the three compared methods while IMOSWFM divided them into mixed pixels, as shown in Figure 7. In terms of the kappa coefficient, IMOSWFM notably outperformed other methods in Area 2 and Area 3. ASWM had the highest kappa coefficient in Area 1. Nevertheless, the kappa coefficient obtained by IMOSWFM was close to ASWM and was higher than DESWE and SSEM.

To quantitatively evaluate the surface water extraction accuracy of the four methods in Area 4, Area 5, and Area 6, we manually delineate the pure water pixels in the three small zoomed-in regions displayed in Figure 8 through visual interpretation using the original 10 m Zhuhai-1 OHS HSIs due to the lack of ground truth water fraction maps. The corresponding OA and kappa coefficient of IMOSWFM and the compared methods are presented in

Table 3. OA and kappa coefficient of IMOSWFM and the compared methods in Region 4, Region 5, and Region 6 with the Zhuhai-1 OHS images.

	Overall Accuracy			Kappa Coefficient
	Region 4	Region 5	Region 6	Region 4	Region 5	Region 6
ASWM	55.18%	91.06%	86.48%	0.1164	0.7738	0.7217
DESWE	51.64%	90.84%	83.91%	0.0843	0.7677	0.6795
SSEM	57.79%	95.11%	85.44%	0.1642	0.8814	0.6935
IMOSWFM	91.09%	96.61%	89.84%	0.8214	0.9191	0.7919

The bold entities represent the optimal values.

. As shown in Table 3, IMOSWFM obtained the highest OA and kappa in all regions. In Region 4, IMOSWFM increased the OA by over 30% and increased the kappa by more than 0.6 compared with the other three methods. The reason behind the extremely low OA and kappa generated by ASWM, DESWE, and SSEM is rooted in the fact that they can hardly distinguish surface water and impervious pixels, consistent with the results displayed in Figure 8. Although the water maps of the four methods in Region 5 were visually similar, IMOSWFM still produced the highest OA and kappa in Region 5, indicating more accurate spatial positioning of river edges. With regard to Region 6, IMOSWFM also achieved the highest OA of 89.84% and kappa of 0.7919, followed by ASWM with an OA of 86.48% and kappa of 0.7217, which demonstrated the outstanding performance of IMOSWFM in the detection of nuanced variations in water fractions, as shown in Figure 8.

The RMSE and SE of IMOSWFM and the compared methods are shown in

Table 4. RMSE and SE of IMOSWFM and the compared methods in Area 1, Area 2, and Area 3 with the Zhuhai-1 OHS images.

	RMSE			SE
	Area 1	Area 2	Area 3	Area 1	Area 2	Area 3
ASWM	0.2700	0.2433	0.2648	0.0461	0.0314	0.0875
DESWE	0.2662	0.2572	0.2782	0.0108	−0.0049	0.0639
SSEM	0.2946	0.2858	0.3144	−0.0165	−0.0341	0.1184
IMOSWFM	0.2506	0.2403	0.2265	−0.0023	−0.0083	0.0516

The bold entities represent the optimal values.

. IMOSWFM produced the minimum RMSE in all areas, revealing that the water fractions estimated by IMOSWFM were most consistent with the reference water fraction maps. The SE value represents the general deviation trend between estimated and reference water fractions, where a positive SE means overesitmation and vice versa. We can see from the SE results that all of the methods overestimated the water fractions in Area 3 and most of the methods underestimated the water fractions in Area 2. The minimum SE values were obtained by IMOSWFM in Area 1 and Area 3. In Area 2, DESWE had the smallest SE value. However, the SE of IMOSWFM in Area 2 was much smaller than ASWM and SSEM, though comparable to DESWE, demonstrating the robust performance of the proposed iterative multi-objective optimization framework for surface water fraction estimation.

4. Discussion

4.1. Comparison Among NDWI, NDWFI, and MNDWFI

To further verify the effectiveness of MNDWFI in the extraction of pure water pixels, here we visually compare the histograms of the NDWI, NDWFI, and MNDWFI and calculate the corresponding OA and kappa coefficient of classification on the three areas around Nanyi Lake. The histograms of the three water indices are illustrated in Figure 10. In general, all of the water indices exhibit the characteristics of bimodal distribution where the two peaks from left to right represent land and water pixels, respectively. The advantages of the MNDWFI mainly embody two aspects: (1) Wider span of values. The values of the MNDWFI have a wider span (from −1 to 1) than the NDWI (from −0.8 to 0.6) and NDWFI (from −1 to 0.7), rendering a longer transition between the two peaks, thus facilitating the precise division of pure and mixed pixels. The range of the NDWI displayed in the results is not between −1 and 1 because only when one pixel has a reflectance of 0 in the NIR band can its NDWI reach 1. Similarly, the NDWI of one pixel can achieve −1 only when it has a reflectance of 0 in the green band. However, the reflectance of pixels of major land cover types in the green and NIR bands is usually not 0 [28,72]. Even for pure water pixels, they have very low reflectance in the NIR band. Therefore, the range of the NDWI can be narrower in real images; (2) More concentrated distribution of pure water pixels. The MNDWFI values of large amounts of pure water pixels distribute around 1, which provides more evident characteristics of numerical distribution for extracting pure water.

Table 5. Comparison of OA and kappa coefficient of NDWI, NDWFI, and MNDWFI on the three areas around Nanyi Lake.

	Overall Accuracy			Kappa Coefficient
	Area 1	Area 2	Area 3	Area 1	Area 2	Area 3
NDWI	89.35%	90.51%	81.99%	0.5923	0.7502	0.6245
NDWFI	91.53%	92.72%	89.64%	0.6533	0.7988	0.7613
MNDWFI	91.74%	93.12%	89.73%	0.6026	0.8080	0.7632

The bold entities represent the optimal values.

shows the OA and kappa coefficient of the NDWI, NDWFI, and MNDWFI across the three areas around Nanyi Lake. The threshold of pure water and land pixels was determined from the histogram by manually selecting the value with the maximum local change rate at the junction between the valley and the right peak. From Table 3 it can be observed that the MNDWFI consistently demonstrated superior performance in extracting pure water compared to the NDWI and NDWFI. The NDWI produced the minimum overall accuracy and kappa coefficients in the three areas since it only exploited few bands of the hyperspectral images. The NDWFI and MNDWFI improved the accuracy of the NDWI by leveraging the full spectral information available in hyperspectral images. Moreover, the MNDWFI further enhanced the classification accuracy of the NDWFI by using the spectrum of water endmembers instead of dark endmembers, which increased the separability between pure water and other pixels, as visually depicted in Figure 10. The outstanding capability of the MNDWFI in pure water extraction makes it a recommended choice for arduous surface water mapping tasks.

4.2. Ablation Experiments

To further verify the performance of each component of IMOSWFM, we carried out several ablation experiments. Specifically, four different combinations of components were considered in this experiment: (1) IMOSWFM. This was the proposed complete algorithm in this paper. (2) IMOSWFM_NDWFI. In this combination, the NDWFI was used, instead of the MNDWFI, to extract pure water pixels, and the iterative unmixing framework was used to adaptively estimate water fractions. (3) IMOSWFM_{non-iterative}. This combination used the MNDWFI to extract pure water pixels, but removed the iterative unmixing framework. Water fractions of all the mixed pixels were estimated in one iteration. (4) IMOSWFM_{NDWFI_non-iterative}. This combination used the NDWFI to extract pure water pixels, and meanwhile removed the iterative unmixing framework. The RMSE and SE of the four combinations on the three areas are illustrated in Figure 11. In terms of RMSE, Figure 11a shows that IMOSWFM and IMOSWFM_{non-iterative} produced smaller RMSE values than IMOSWFM_NDWFI and IMOSWFM_{NDWFI_non-iterative} in all three areas, respectively, indicating the effectiveness of the MNDWFI in accurately segmenting water bodies. The performance of IMOSWFM deteriorated when the iterative unmixing framework was not adopted, since IMOSWFM_{non-iterative} had lager RMSE values than IMOSWFM and IMOSWFM_{NDWFI_non-iterative} had larger RMSE values than IMOSWFM_NDWFI. Therefore, the effectiveness of our proposed iterative unmixing framework can be verified. Within the realm of SE, Figure 11b shows that the SE values generally became larger when the MNDWFI was not used or the iterative unmixing framework was not applied in Area 1 and Area 2. In Area 3, IMOSWFM had the maximum SE. However, given that SE represents the general deviation trend between estimated and reference water fractions, a small SE does not necessarily indicate a high accuracy of water fraction estimation. From Figure 11a we can see that IMOSWFM had the minimum RMSE in Area 3. Thus, the indispensability of each component of IMOSWFM can be proved.

4.3. Convergence Analysis

To analyze the convergence of IMOSWFM, convergence curves of the two objective functions optimized by IMODPSO in Area 1 of one certain run are illustrated in Figure 12. Five iterations were finally conducted in this run. As the evaluation times increased, the values of both objective functions decreased rapidly in the early stage and gradually converged to a small value at the end of optimization, unveiling the efficient convergence of IMOSWFM, thus ensuring reliable endmember combinations for accurately estimating water fractions. Figure 13 complements this analysis by depicting the number of remaining mixed pixels at different iterations. It can be discerned that the number of remaining mixed pixels reduced constantly as the iteration increased, since the mixed water–land pixels, which were accurately unmixed, were removed from the original image at the end of each iteration. There was a sharp decline after the first iteration because the pure water and non-water pixels were also removed from the original image in the first iteration in addition to the accurately unmixed mixed pixels. After four iterations, the number of remaining mixed pixels had been reduced to a small value which met the termination condition, and the algorithm finally stopped.

4.4. Error Analysis

Although IMOSWFM outperformed other methods in the three areas, it is undeniable that differences still existed between the water fraction maps generated by IMOSWFM and the reference maps. As shown in Figure 7, the main source of errors stems from the extraction of pure water and mixed water–land pixels. Firstly, quantities of land pixels located at the boundaries of adjacent water bodies were mistakenly divided into mixed pixels. Secondly, a lot of small water bodies were not extracted by IMOSWFM. The reason is that only one set of thresholds was utilized to segment the whole image. IMOSWFM adopts a histogram-based double threshold segmentation method on the assumption that pure water and land pixels can be expressed by the two distinct peaks in the histogram of the MNDWFI. However, due to the variation in illumination and atmospheric conditions, the scales of the MNDWFI in different regions change dramatically, resulting in a wide overlapping range between the pure and mixed pixels in the histogram. Specifically, we observed that pure land pixels bordered by a large water area generally have higher MNDWFI values than the pure water pixels located in some isolated small water bodies. If the land threshold t₁ is set high enough to distinguish the former, the latter will not be successfully extracted since the MNDWFI values of these pixels are lower than t₁. On the contrary, if the land threshold t₁ is set low enough to extract the latter, the former will not be properly divided into land pixels since the MNDWFI values of these pixels are higher than t₁. The same principle applies to the water threshold t₂. Therefore, errors occurred inevitably when pure water pixels were extracted globally by the double threshold segmentation method in IMOSWFM. A local thresholding method which applies a varying threshold for different regions of the image may be an effective way to address this problem [73,74,75].

Residual errors also exist in the estimation of water fractions. IMOSWFM adopts an iterative framework to adaptively unmix mixed water–land pixels. In each iteration, only the mixed pixels with RMSE smaller than the user-defined threshold r obtained the water fractions. A smaller RMSE usually means a better unmixing result [61,76,77]. Theoretically, setting a smaller r can lead to more accurate water fraction estimation results for mixed pixels in each iteration, as the RMSE of all these pixels would be smaller than r. However, the improvement of accuracy comes at the expense of efficiency. An r too small can lead to an extremely high computational time because few pixels are able to be successfully unmixed in each iteration. Therefore, it is of vital importance to find an appropriate r which can trade off both the efficiency and accuracy. Another issue arises from the determination of the number of endmembers. IMOSWFM assumes that each mixed pixel consists of three types of endmembers, including water, vegetation, and impervious surface. This assumption may not apply to some pixels containing complex ground features, which gives rise to the residual errors between the estimated and reference water fraction maps. Multitasking optimization can be used in IMOSWFM to simultaneously optimize solutions containing different numbers of endmembers and adaptively search the optimal number for different mixed pixels [78].

5. Conclusions

In this research, we introduced IMOSWFM, an adaptive unmixing method based on iterative multi-objective optimization for surface water fraction mapping from hyperspectral images. In IMOSWFM, a modified normalized difference water fraction index, named the MNDWFI, was proposed to exploit each spectral band of HSIs for the accurate extraction of pure water. Also, an iterative unmixing framework was designed to obtain a more precise water fraction map where optimal endmembers were adaptively extracted for different mixed water–land pixels, leveraging the robust global search ability of multi-objective evolutionary algorithms. Experiments on the three Zhuhai-1 HSIs demonstrated the effectiveness of IMOSWFM in surface water fraction mapping. From the obtained water fraction maps, IMOSWFM was able to preserve the details of water fraction variation and obtain robust water fraction estimation results in diverse image scenes regardless of spectral variability. Specifically, IMOSWFM had the highest OA of 91.74%, 93.12%, and 89.73%, and the lowest root-mean-square errors of 0.2506, 0.2403, and 0.2265 in the three primary study areas, indicating the exceptional capability of IMOSWFM in both pure water extraction and water fraction estimation. In our future work, we will strive to solve the problem of determining multiple thresholds for pure water extraction and automatically estimating the number of endmembers. As discussed in Section 4.4, a single threshold for the segmentation of pure water may result in errors because the scales of the water index vary greatly. A novel pure water extraction method needs to be proposed to find optimal thresholds for pixels from different regions. Moreover, to further improve the accuracy of water fraction estimation, it is also important to dynamically find the optimal number of endmembers in each iteration. An advanced spectral mixture model can also be considered in the unmixing process to deal with the complex interactions between multiple constituent components within the mixed pixels.