Delineation of Intermittent Rivers and Ephemeral Streams Using a Hybrid Method

Abstract

Intermittent rivers and ephemeral streams are crucial for the water cycle and ecosystem services, yet they are often neglected by managers and researchers, especially in headwater areas. This oversight has caused a lack of comprehensive basemaps for these vital river systems. In headwater regions, water bodies are typically sparse and disconnected, with narrow and less distinct channels. Therefore, we propose a novel hybrid method that integrates topographic data and remote sensing imagery to delineate river networks. Our method reestablishes connectivity among sparsely distributed water bodies through topographic pairs, enhances less distinct channel features using the gamma function, and converts topographic and water indices data into a weighted graph to determine optimal channels with the A* algorithm. The topographic and water indices data are derived from the Multi-Error-Removed Improved-Terrain DEM (MERIT DEM) and an average composite of the Modified Normalized Difference Water Index (MNDWI), respectively. In the upper Lancang-Mekong River basin, our method outperformed five publicly available DEM datasets, achieving over 91% positional accuracy within a 30 m buffer. This hybrid method enhances positional accuracy and effectively connects sparse water bodies in headwater areas, offering promising applications for delineating intermittent rivers and ephemeral streams and providing baseline information for these river systems.

1. Introduction

Intermittent rivers and ephemeral streams (IRES) are river systems characterized by discontinuous water flow throughout the year [1,2,3]. These systems [2] are crucial components of global river networks, comprising approximately 89% of the total river length [4]. Far from being exceptions, IRES are prevalent [5], particularly in arid and semi-arid regions [6], and play a vital role in natural ecosystems and human societies within drylands [7,8]. They provide indispensable ecosystem services such as fresh water and food supply [9], pollutant and nutrient cycling [10], catalyzing biogeochemical reactions [11], sediment transport [12], surface and groundwater exchange [13], stormwater management [14], and recreational opportunities. These services are sensitive and fluctuate in response to direct human activities and climate change [15,16,17]. Unfortunately, IRES often suffer from long-term mismanagement and neglect by researchers and managers [18]. The baseline characteristics of these river systems are poorly understood due to a significant lack of attention and the challenges posed by their sparse and intermittent distribution.

To achieve a comprehensive understanding of IRES, researchers have extensively mapped these water bodies, generating numerous basin-scale or national-scale river network datasets [19,20,21,22]. However, detailed maps of these river systems remain unavailable for many developing countries [23], particularly in mountain headwater areas. Even in regions with previous surveys, results are often imprecise [20]. Accurate delineation of IRES is crucial for precise research, with remote sensing emerging as a promising approach for large-scale, high-resolution mapping [24,25], especially in remote and inaccessible areas such as the High Mountain Asia region. Therefore, substantial methods for delineating river networks have been developed, broadly categorized into digital elevation model-based methods, remote sensing image-based methods, and hybrid methods.

In digital elevation model-based methods, DEM data are acquired through various techniques such as stereo photogrammetry, Synthetic Aperture Radar (SAR), or Light Detection and Ranging (LiDAR) from different satellite missions. These data are extensively used in hydrologic modeling and the delineation of IRES [26,27]. Parameters like direction, slope, accumulation, aspect, and watershed sink can be automatically extracted pixel-wise from the DEM using algorithms like D8 [28] or D-infinity [29], contributing to the generation of river networks. The quality of delineated river networks heavily relies on the resolution and precision of the DEM. Commonly used DEM products include the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) at 30 m resolution and the Shuttle Radar Topography Mission (SRTM) at 1 or 3 arc-second resolutions [30]. Higher resolution DEM products, such as the AHN 0.5 m for the Netherlands [31], the 5 m DEM for Australia, the RGE ALTI 1 m for France, and the USGS 3DEP 10 m for the United States, derived from LiDAR technology, have emerged in recent years. However, the coverage of these datasets is limited. Global or national-scale DEM datasets are seldom updated, which means that the derived river networks may not accurately reflect the dynamics of rivers, particularly smaller ones in headwater areas that are more sensitive to climate change. Anthropogenic features are also inadequately considered in DEM algorithms. Moreover, errors in DEM in specific regions can lead to computational inaccuracies in water sinks or flow directions [32].

In remote sensing image-based methods, the delineation of IRES relies heavily on satellite images, including both optical imagery and SAR. SAR is effective for detecting water bodies due to its sensitivity to backscatter from the water surface [33]. Numerous research endeavors have explored these methods across various scales and landscapes. However, the intricate nature of processing SAR signals has constrained its application in large-scale water body extraction [34]. Mapping water extents from remote sensing images involves leveraging spectral features for pixel-wise applications and object features for object-based methods, achieved through either unsupervised or supervised methods [35]. The water index, a highly efficient category of unsupervised methods for optical imagery, finds widespread application in large-scale basins [36]. Examples include the Normalized Difference Water Index (NDWI [37]), Modified Normalized Difference Water Index (MNDWI [38,39]), Multispectral Water Index (MuWI [40]), Automated Water Extraction Index (AWEI [41]), and Land Surface Water Index (LSWI [42]). Following water body maps, water pixels can be identified and connected to form multiple river lines, which are then organized into river networks. However, a significant limitation of satellite imagery-based delineation is the incomplete nature of derived river centerlines for each basin due to the absence of water points. This lack hinders the connectivity required for constructing river networks in headwater areas.

Given the limitations of existing methods, a hybrid approach has emerged that integrates DEM and remote sensing imagery to delineate river networks. This approach fully leverages the topological features derived from DEM and the water body features mapped from remote sensing imagery, balancing the requirements for topological correctness and positional accuracy. For instance, Wang et al. (2021) proposed a method known as “Remote Sensing Stream Burning”, which integrates a high-resolution river network using 10 m Sentinel-2 imagery at a basin scale. The method extracts water bodies from optical images, identifies flow directions pixel-wise based on the MERIT DEM, and generates a basin-scale river network for the Lancang-Mekong River. It also addresses the challenges posed by IRES in the upper reaches of the river. Similar research can be found [20,43] where SAR images are utilized. These studies focus more on river lines that occur in the middle or lower river reaches. However, a significant challenge remains in the disconnection phenomenon in the upper reaches, potentially leading to measurement errors and an underestimation of length, area, and associated biogeochemical process statistics.

Therefore, this study aims to develop a robust method that integrates remote sensing imagery with MERIT DEM to delineate IRES in the upper reaches. These regions typically feature sparse and disconnected water bodies, with less distinct river channel features. The MERIT DEM dataset effectively reduces data errors [21]. The extensive remote sensing imagery from the Landsat mission series holds the potential to capture intermittent watercourses. We first extract key points from the water bodies to indicate critical junctions in the river network. Then, we build the topological relationships among each pair of key points to generate potential river channels. Finally, the A* algorithm [44] is employed to construct a graph with slope and water indices data to find the optimal channel between each pair of key points. We selected a basin in the upstream region of the Lancang-Mekong River basin to validate our proposed method, given its diverse landscapes with IRES. To assess the accuracy of our method, we compared the extraction results with various ground-truth images. Our research advances river network delineation in three key ways compared to previous studies, as follows: (a) it focuses on delineating IRES in headwater areas, where water bodies are usually sparse and disconnected; (b) it utilizes the gamma function to enhance the visibility of river channel features; (c) it incorporates graph theory to integrate water bodies and topological data, using these as weights to identify optimal river channels. This approach to delineating IRES could establish a methodological baseline for understanding the dynamics of headwater areas.

2. Materials

2.1. Study Area

The Lancang-Mekong River, one of the world’s major river basins, spans multiple countries, with Changdu city marking the dividing line of its upper reaches. As shown in Figure 1, our study area is located in the upstream region of this basin, specifically between 93.8°E and 97.2°E longitude and 31.5°N and 33.8°N latitude. It is bordered by the Yangtze River basin and the Nu-Salween River basin and constitutes a portion of the Three-River-Source National Park in China [45]. Nestled in the southeastern part of the High Mountain Asia region—the highest plateau on the planet—the area encompasses numerous small pools, lakes, intermittent rivers, and ephemeral streams. With elevations ranging from approximately 3485 to 5780 m above sea level, the region has a very low population density. The extensive presence of intermittent and ephemeral watercourses, which typically have a river order of less than 3, results in discontinuous water body scenes within remote sensing imagery. To tackle this challenge, we implemented our proposed method. The boundaries of the study area were established using both the perimeters of the Three-River-Source National Park and global basin boundary data sourced from the HydroBASINS dataset [46].

2.2. Source Data

2.2.1. Water Body Data

A global surface water dataset [47] (https://global-surface-water.appspot.com/#, accessed on 2 January 2024) was employed as one of the input sources for the proposed method. This dataset, developed by the European Commission’s Joint Research Centre, is derived from over three million Landsat satellite images spanning the past 38 years (1984 to 2021) at a high resolution of 30 m. The dataset encompasses various map products, including surface water occurrence, surface water recurrence, surface water seasonality, surface water transition, and maximum water extent. The maximum water extent product, which provides information on all locations identified as water bodies throughout the 38-year period, was selected as the potential water body data to support our experiment.

2.2.2. Water Indices Data

In addition to the water body product, water indices data were necessary as they play a crucial role in indicating water differential information. In our research, interannual changes were excluded in the upstream channel of the river network, and the spatial resolution for our studied river network was 30 m. To mitigate random factors such as temporary water ponds, we employed a multi-year average process using Landsat images from 2013 to 2022. These images were processed on the Google Earth Engine platform, utilizing the MNDWI. The MNDWI is a classical index for surface water detection and has been widely applied in research. It enhances the detection of surface water information, making it suitable for our study of water bodies in the headwater area, where water features may be less distinct. The resulting multi-year water indices data were exported with the same spatial extent and projection information as the global surface water dataset.

$$MNDWI=\frac{b_{green}-b_{MIR}}{b_{green}+b_{MIR}}$$

(1)

Here, $b_{green}$ and $b_{MIR}$ are the surface reflectance values of the Landsat green band and middle infrared band, respectively.

2.2.3. Topography Data

A global topography dataset [48] known as the Multi-Error-Removed Improved-Terrain DEM (MERIT DEM) was utilized as the topography prior data with a resolution of 90 m. Notable for its comprehensive error correction, the MERIT DEM has removed various error components from the existing spaceborne SRTM3 v2.1 and AW3D-30m v1 DEM [49]. These components include absolute bias, stripe noise, speckle noise, and tree height bias. By minimizing technical disturbances during the derivation of river networks, this refined topography data significantly enhances the precision of research analyses. Subsequently, a global hydrography dataset, MERIT Hydro [22], was developed primarily based on the MERIT DEM. This dataset includes crucial features such as flow direction, flow accumulation, hydrologically adjusted elevations, and a river channel width product, which are essential for comprehensive hydrological studies.

2.2.4. Vector Data

The vector data used in our research were primarily sourced from the HydroSHEDS database [46], a comprehensive resource offering global digital layers essential for various research domains, including hydrology, environmental science, conservation, socioeconomics, and health. These layers encompass critical hydrological features such as catchment boundaries, river networks, and lakes, available in diverse resolutions and scales. For our study, we leveraged the HydroBASINS product from HydroSHEDS, which provides a hierarchically nested sub-basin breakdown on a global scale with a resolution of 15 arc-second. Additionally, the boundary of the Three-River-Source National Park was employed to define and generate the spatial extent of our study area. To support our validation, we also incorporated publicly available river networks covering the study area with spatial resolutions finer than 90 m. Given the limited availability of 30 m resolution datasets for our study area, we chose to include one dataset with a 30 m resolution [19] and two with a 90 m resolution [27,50] for river network analysis.

2.2.5. Comparator Data

The comparator data used in our study were sourced from various publicly accessible datasets, all featuring a spatial resolution of 1 arc second (approximately 30 m). These included the SRTM DEM, ASTER DEM, NASA (National Aeronautics and Space Administration) DEM [51], GLO-30 (Copernicus DEM–Global 30 m) DEM [52], and AW3D30 (ALOS World 3D–30 m) DEM [53]. Details regarding the providers and temporal extent of these datasets are provided in

Table 1. Provider and temporal extent for five selected DEM products.

Name	Provider	Temporal Extent
SRTM	NASA & CGIAR	11 February 2000–22 February 2000
ASTER	NASA & METI	1 March 2000–30 November 2013
NASA	NASA & USGS & JPL-Caltech	11 February 2000–22 February 2000
GLO-30	ESA	1 December 2010–31 January 2015
AW3D30	JAXA	24 January 2006–12 May 2011

. It is important to note that the NASA DEM is a reprocessed product derived from multiple data sources, while the others were generated directly from radar or optical missions.

3. Methods

3.1. Overview

As depicted in Figure 2, our proposed method for delineating IRES consists of three main components: deriving key points, building topology, and delineating the river network. This is followed by a result validation step. The process steps are executed sequentially. To begin, key points are identified from the maximum water extent product. These key points, along with the MERIT flow direction data, are utilized as inputs to generate topology pairs. Subsequently, the A* algorithm is employed to plan the optimal river channel between each pair of key points, utilizing the enhanced MNDWI data and the topological data derived from the MERIT dataset. Finally, the detailed delineation of IRES within the study area, particularly in the headwater region, is established. Further details are provided in the following sections.

3.2. Preprocessing

As detailed in Section 2, the maximum water extent product is a composite image compiled from multiple years of data, encompassing all previously detected water bodies within the study area. Since our focus is not on ponds and lakes, and recognizing their potential to affect the precision of our method, we manually excluded them from our study area. To minimize the impact of clouds on the generation of water indices data, a cloud mask was applied to each remote sensing image. Furthermore, to enhance the robustness of our method, we implemented an outlier removal algorithm, which reduces the influence of outlier water bodies typically found in isolation.

The raster products, including the maximum water extent products, the water indices data, and the MERIT DEM, were preprocessed to clip, resample, and reproject them. This ensured they matched the tile size, projection, and metadata of the maximum water extent product, which served as the reference image. Earlier research has often concentrated on river networks wider than three pixels of Landsat imagery, approximately 90 m, due to the significant uncertainties associated with narrower waterways. In this study, we assume that the maximum water extent product provides an accurate representation of the maximum observable water extent at a 30 m resolution for IRES.

3.3. Deriving Key Points

The key points are crucial elements for our proposed method. They represent actual water bodies and highlight critical locations within the potential IRES, such as river confluence points or relay points. They are essential for the accurate delineation of river networks. In our study, the key points were generated from the maximum water extent product and the mean water indices data by employing the following three image-processing algorithms sequentially: the alpha shape [54], skeleton, and local maximum algorithms. As depicted in Figure 3a, a subset of the maximum water extent product shows water bodies sparsely and disconnectedly distributed in the river channels, with water body density gradually increasing as the river flows downstream. However, a common phenomenon is the lack of connectivity among water bodies, and the connectivity from upstream to downstream remains unclear—a typical characteristic of IRES in the headwater area. In this context, these key points could exert a more significant influence on delineating non-perennial rivers than on large perennial rivers.

As shown in Figure 3a, the water bodies, depicted in blue, are dispersed in a fragmented and disconnected manner. To resolve this, we utilized the alpha shape algorithm, which is mathematically well-defined and serves as a generalization of the convex hull and a subgraph of the Delaunay triangulation [55]. This algorithm was employed to envelop the water bodies into multiple alpha shape objects that closely resemble the original shape of the water bodies, as depicted in Figure 3b. Subsequently, these alpha shape objects underwent processing using a morphology labeling algorithm to identify connected components of the alpha shape bodies. The alpha shape objects were then separately regarded as the water body mask and the water body extent, which had been previously detected or classified as water bodies. For the water body mask, the mean water indices derived from the Landsat image series were masked, as shown in Figure 3e. The local maximum algorithm was then applied to generate local peaks in water indices. For the water body extent, a skeleton algorithm was employed to generate skeleton lines, as shown in Figure 3c. The intersections and endpoints derived from these skeleton lines were subsequently relocated to the nearest local peaks of water indices, which were identified as the final key points of the current alpha shape object, as illustrated in Figure 3f.

3.4. Building Topology

Building topology is the second step in delineating IRES, with the goal of establishing a matrix of relationships among all key points derived from Section 3.3. The primary standards for this step are referenced from hydrologic analysis theory, which has been integrated into various software applications and widely applied in numerous research studies over the past few decades. The core source data for topology construction are the flow direction product of MERIT DEM data, as described in Section 2.2.3, and the derived key points from Section 3.3. These are separately illustrated in Figure 4a,b. The MERIT dataset has been refined to address multiple error components and is used solely to indicate the flow direction between each pair of key points. Therefore, the relatively coarse spatial resolution of the MERIT dataset does not influence the accuracy of delineating IRES.

All key points were initially labeled based on the alpha shape objects they belonged to, as identified by the connected component algorithm. First, we computed a list of key points in descending order of their elevation, derived from the MERIT DEM data, to establish the processing order. Subsequently, we systematically searched for connected channel paths for each key point in the ordered list and used the D8 algorithm to find the nearest downstream key point. The acquired downstream key points were organized into a tuple structure with their upstream key points and added to the established relationship matrix if a required downstream key point was found. If a downstream key point could not be found and the search reached the basin boundary, the intersection point was designated as the outlet of that basin. After all key points in the ordered list had been processed, the relative relationships among these key points were completely explored, and the topology matrix was established. The search results can be seen in Figure 4c, where the red-dashed lines indicate key point pairs or topological paths with a topological relationship, including two key points from the same alpha shape object. The purple-dashed lines depict key point pairs involving two key points from different alpha shape objects. The topology pairs within the matrix were derived from the MERIT Direction data and represent only the upstream and downstream relationships. They are used to generate potential channel path buffers, which can reduce computational load, for accurate extraction in the subsequent section, as shown by the solid red and purple lines in Figure 4c.

3.5. Delineating the River Network

This section aims to use the derived auxiliary data from the preceding steps for the accurate delineation of IRES. The features of the available data are primarily determined by their pixel values, with the MERIT DEM and water indices data serving as two important data sources. Hydrological analysis is intended to derive the channel path of the river network with the steepest descent in elevation from the DEM, while the remote sensing image-based method aims to derive the central lines. To fully leverage the strengths of both methods, we have defined a function that integrates features from both the MERIT Slope product and water indices data, depicting the IRES. The formulated equation is shown below:

$$y=f\left(Slope,\Gamma \left(Indices\right)\right)$$

(2)

Here, $Slope$ represents the derived MERIT Slope product, $Indices$ represents the water indices data, and the $\Gamma$ is the gamma function utilized to enhance the water indices data.

In this section, the source data include the Slope data, the water indices data, and the topology matrix data. First, the water indices were processed using the gamma function to enhance their features, given that the water indices data in the headwater regions of the IRES tend to be weak and pose detection challenges compared to midstream and downstream regions. As illustrated in Figure 5b,c, both the derived Slope data and the enhanced water indices data were normalized to the same range, and the desired combined values were calculated using Equation (2). This confirms that our proposed method can utilize both the topology and the water indices information together. Finally, the A* algorithm was employed to delineate the IRES. Originally designed to find the shortest paths between nodes in a weighted graph [56], the topology matrix could be split into several pairs of key points that meet the requirements of the A* algorithm, representing the nodes. The graph’s extent was the buffer area of the channel paths derived from Figure 4c for each pair of key points from upstream to downstream. The constructed values were considered as the weight of each adjacent pixel in eight directions. Implementing these steps in parallel would yield the desired river lines between two key points, extending to cover all key points across the entire study area.

3.6. Positional Accuracy Metric

To quantitatively assess the performance of our proposed method, it is essential to outline the adopted positional accuracy evaluation metrics. In this research, we have chosen completeness as the metric for evaluating positional accuracy, a method that has been utilized in several studies [20,21]. It represents the percentage of the river network derived from our method and various comparator river networks captured by the reference dataset (the ground truth data). Given that there are three types of ground truth basemap layers, multiple comparisons are generated. The equation for assessing completeness is provided below.

$$completeness=\frac{{River}_{network}}{{River}_{ground\_truth}}\times 100\%$$

(3)

Here, ${River}_{network}$ represents the river network obtained by our proposed method or other comparative methods, whereas ${River}_{ground\_truth}$ refers to the buffer zone of the reference river network, defined by a specific buffer size.

4. Experiments and Results

This research focuses on delineating IRES in headwater regions using our proposed method. The delineation results and positional accuracy evaluations are presented in the subsequent sections. Section 4.1 presents the delineation results of our proposed method, along with corresponding visual comparisons to ground truth images. Section 4.2 presents the positional accuracy assessment among various river networks derived from different available DEM datasets.

4.1. Results of Intermittent Rivers and Ephemeral Streams

As described in Section 2.1, the upper Lancang-Mekong River basin was selected as the study area. Initially, our proposed method was employed to extract IRES. Subsequently, three publicly available dynamic basemap layers—produced independently by Esri, Google, and Tianditu corporations for use in ArcGIS 10.6 software—were introduced to serve as the ground truth data source. These dynamic basemap layers encompass a range of remote sensing image tile resolutions globally and have been widely applied.

As depicted in Figure 6, our delineated IRES are presented in full view for visual comparison. We observe that the upper part of the river basin has a higher density of river networks. Statistically, there is a total of 11,693 km of the selected river network, and more than 88.4% are less than 90 m wide. As shown in Figure 7, the three sources of basemaps exhibit overall consistency in terms of landscapes, although some color variations exist. Regardless of the data source used for the ground truth basemap, the river lines consistently lie within the river channel, showing no noticeable deviation. While the derived IRES may deviate from the reference ground truth basemap in the broader riverbed, they consistently remain within the potential riverbed. Additionally, the sparsely located maximum water extent points are all connected by our derived results in both example regions. This example illustrates that our proposed method can accurately extract IRES in headwater areas, especially where water bodies are sparse and disconnected.

4.2. Results of Positional Accuracy Assessment

In addition to visual comparison with the ground truth basemap data, the delineated IRES were also compared with other publicly accessible datasets, as described in Section 2. These datasets include the SRTM DEM, the ASTER DEM, the AW3D30 DEM, the GLO-30 DEM, and the NASA DEM, each with a resolution of 30 m, equivalent to that of our adopted global surface water dataset. The reference data for comparison were obtained from river lines manually delineated from the dynamic basemap layers of Esri, Google, and Tianditu, respectively. Delineation of river networks was performed for each of the five DEM datasets using hydrological analysis, which included steps such as sink filling, flow direction determination, flow accumulation calculation, network extraction, and raster-to-vector river line conversion. All of these processes were implemented in ArcGIS 10.6 software, employing the same parameters and procedures.

For convenience, all derived IRES were converted into vector format. Since comparing the derived topology of river lines is challenging, we employed the buffer of the ground truth data to calculate the related metrics, as mentioned in Section 3.6. Buffer sizes were set at 15 m, 30 m, 45 m, and 60 m. As described by Equation (3), we calculated the intersection between the ground truth buffer zone and the delineated river network within the same basin. We then determined the percentage of the length of IRES that fell within the buffer zone.

A quantitative assessment was conducted to evaluate our proposed method for delineating IRES, with the results presented in

Table 2. Positional accuracy assessment of intermittent rivers and ephemeral streams at different buffer sizes.

Source	Method	15 Meters	30 Meters	45 Meters	60 Meters
Esri	Ours	62.2%	89.5%	95.3%	97.4%
	SRTM	37.6%	66.3%	80.6%	88.9%
	ASTER	27. 5%	52.9%	68.1%	80.4%
	NASA	38.5%	66.7%	80.5%	88.1%
	GLO-30	60.4%	86.8%	93.2%	95.2%
	AW3D30	46. 7%	74.7%	88.1%	91.7%
Google	Ours	69.7%	92.6%	96.5%	98.5%
	SRTM	36.1%	64.4%	79.6%	89.1%
	ASTER	27.3%	51.1%	68.9%	81.3%
	NASA	36.4%	64.3%	79.2%	88.2%
	GLO-30	56.5%	82.5%	92.4%	96.8%
	AW3D30	43.8%	72.9%	86.3%	92.2%
Tianditu	Ours	66.9%	93.3%	97.4%	98.4%
	SRTM	35.6%	67.9%	80.4%	88.7%
	ASTER	27.5%	50.6%	68.3%	80.1%
	NASA	38.5%	66.2%	80.2%	89.3%
	GLO-30	57.4%	85.9%	94.1%	97.1%
	AW3D30	46.9%	75.1%	88.4%	92.9%

. It is evident that our proposed river network delineation method outperforms all other comparison datasets. The comparison river networks were all derived from DEM datasets, given the scarcity of publicly available river network datasets at a 30 m resolution. The GLO-30 DEM dataset performs the best among the comparison datasets, while the ASTER DEM performs the worst across all three types of ground truth and four buffer sizes. Accuracy metrics for all methods increase with the enlargement of the buffer size across all three ground truth sources. The increase is most significant from a buffer size of 15 m to 30 m, and it gradually slows as the buffer size increases further, with the lowest increase observed from a buffer size of 45 m to 60 m. Additionally, the accuracy metrics among different data sources exhibit similar patterns, representing the three types of data sources.

To further validate the performance of various DEM datasets in comparative experiments, we present visual results in Figure 8. As shown in Figure 8, the IRES derived using the hydrological analysis method based on defined thresholds extracted more river lines than our proposed method. Relative to the ground truth line in light blue, only our proposed method and the GLO-30 DEM dataset correctly delineate the river lines, with our method showing better performance than the GLO-30 DEM dataset. The other four datasets derive different river lines and cannot accurately delineate the specified river lines as depicted in the ground truth tile. Additionally, there are some potential river lines in all five comparative subfigures, and the channels differ from one another.

For river sections with significant elevation differences, each of the comparative DEM datasets can relatively accurately delineate IRES, as the actual river lines usually align with the maximum elevation difference. The elevation difference is much larger than the errors, despite potential variations in accuracy. However, for river segments with minor elevation differences, different datasets exhibit noticeable disparities in performance, as elevation errors may dominate the river flow direction and lead to incorrect delineation of IRES. From this perspective, the importance of coupling both remote sensing imagery and DEM datasets to extract IRES becomes evident.

5. Discussion

5.1. Positional Accuracy Assessment and Comparison Scheme

The positional accuracy assessment scheme employed in our experiment has been previously used in studies [20] and [21] to benchmark positional accuracy among comparators. In these works, several freely accessible river networks from previous research that encompass the study area were treated as the ground truth for reflecting positional accuracy. Examples include the OpenStreetMap (OSM) [57] dataset and the Global River Widths from Landsat (GRWL) [19], which include river lines with a 30 m resolution but consider only those wider than 90 m. However, there is a scarcity of publicly accessible river network datasets that meet the required resolution or integrity standards, as our study focuses mainly on IRES. Both OSM and GRWL were not chosen due to their limited contributions and primary coverage focus on only the main stem of the river network within our study area. Other river network datasets derived globally or continentally typically rely on raster images with lower spatial resolutions. Examples include the Basin90m [27] dataset and the HDMA [50] dataset, both of which have a resolution of 90 m. Given that experiments with varying spatial resolutions might result in comparisons that lack significant relevance, existing river networks are utilized solely for visual comparison. Finally, we derived river networks from five types of DEM datasets using hydrologic analysis with a 30 m spatial resolution for the positional accuracy assessment experiment.

In addition, we compiled three types of sources—composite globally sourced remote sensing image layers from Tianditu, Google, and Esri—to serve as the ground truth image. To organize the comparative experiment, the ground-truth river lines were manually mapped according to the three aforementioned dynamic basemaps. We established four buffer groups, each with a size of 15 m, 30 m, 45 m, and 60 m, to evaluate the corresponding experimental results. This is similar to the work in [20], which set three buffer groups of 10 m, 20 m, and 30 m for the Sentinel-1 imagery. As illustrated in Table 2, the positional accuracy assessment results of our method excel across all four buffer size groups and all three dynamic basemaps, validating the merit of our approach. In addition to the results from our method, the results derived from the GLO-30 DEM dataset, produced by the Copernicus program of ESA (European Space Agency), outperform the other comparative DEM datasets. Furthermore, the AW3D30 DEM dataset, generated by JAXA, ranks third. This outcome is predominantly influenced by the generation algorithm and the acquisition time span of each DEM dataset. Both the GLO-30 and AW3D30 DEM datasets were acquired more recently than the SRTM, ASTER, and NASA DEM datasets. They may represent elevations from different time periods; the longer the time span, the greater the elevation differences [49]. Consequently, DEMs captured at different times may introduce varying degrees of error. The impact on delineation methods primarily depends on the time span between DEM acquisition time and ground-truth image capture. Similar circumstances have been documented in previous research regarding the performance of publicly accessible DEM datasets in High Mountain Asia [49]. Additionally, the visual comparison results depicted in Figure 9 reveal that all three previous river networks failed to concentrate on the headwater areas. In contrast, our derived river network effectively delineates these areas, corresponding to the density of the river network.

5.2. Linking Topography with Water Body Position

Our experiment is predicated on two fundamental assumptions. First, we assume that the river lines within our study area on the Tibetan Plateau have remained consistent from 2000 to 2020. This period coincides with the timeframe of our water indices data and the creation of various DEM datasets. This assumption stems from the fact that our study area is characterized by a sparse human population and minimal human-induced modifications to the natural landscape [49]. Particularly in the headwater region within the Three-River-Source National Park, a designated protected area where the influence of human activity is minimal, we anticipate that the river lines have been relatively stable throughout the period under study. Second, we consider the water bodies, which are sparsely distributed across the maximum water extent product, to be authentic water pixels [47]. These pixels denote IRES. Both of these assumptions are deemed equally applicable across all comparators.

According to these assumptions, our goal is to construct detailed river networks in the headwater area, which features plateaus with numerous IRES. The water bodies in this area are typically sparse and disconnected from each other to a large extent, with river channel features that are less distinct than those found in the middle and lower reaches of rivers. These topographical and surface water body features are significantly different from those in previous related study areas. Moreover, these river lines are relatively narrow and susceptible to positional errors. Consequently, conventional river network delineation methods, such as hydrologic analysis [26] or water body extraction [34], prove insufficient for meeting our specific delineation requirements. To address these challenges, we have proposed a novel approach that employs the gamma function to enhance water features within the water indices. We also build topography pairs among the sparsely distributed water bodies using DEM data and have adopted the A* algorithm to identify the optimal river channels.

This hybrid method enables us to integrate topographic information from DEM data with the positional data of water bodies derived from remote sensing imagery. It offers a comprehensive solution for delineating river networks in challenging environments. Compared with the previous hybrid method in research [21], our method is particularly focused on IRES that typically appear in the headwater area, where the harsh natural environment makes manual surveys difficult. It is capable of extracting river lines between sparse and disconnected water bodies, as depicted in high-resolution remote sensing images. Our method provides headwater information that has often been omitted from large-scale river network datasets for an extended period of time. The refined river networks generated can serve as a valuable supplement to existing datasets [19,27,50], establishing a baseline for authentic scenarios and aiding in the creation of river networks over broader areas.

5.3. Limitations and Future Work

Despite the promising performance of our proposed hybrid method within the study area, several limitations require further discussion. First, the accuracy of the generated river channel paths depends on the correctness of key points derived from the global maximum water extent product and the topographic information provided by DEM data. Although numerous sources of DEM data are publicly accessible, a few are as refined as the MERIT DEM, which has been validated by previous research [21]. We utilized the MERIT DEM in our experiment to supply topographic information, opting for its quality over higher-resolution DEM datasets [22]. Additionally, the resolution of the river channel paths largely depends on the resolution of the MNDWI image, with the MERIT DEM exerting minimal influence on this aspect. However, discrepancies compared to real-world conditions persist, potentially leading to erroneous results in some areas. We have chosen to forgo some spatial resolution to enhance vertical resolution. Furthermore, when available, a high-resolution DEM, along with a detailed surface water product, can provide more precise topography data and surface water features. This enhanced information would allow our method to delineate IRES with greater accuracy. Second, computational efficiency is constrained by two main factors: the construction of topography pairs and the delineation of river lines using the A* algorithm. The process of topography pair construction is similar to classical hydrologic analysis, utilizing flow direction data [32]. While hydrologic analysis generates flow accumulation to identify potential river channels, our algorithm creates topography pair information as preliminary data for subsequent path creation. The A* algorithm, though more efficient than Dijkstra’s algorithm, is still computationally intensive.

Future work could focus on improving computational efficiency by conducting distributed processing on more powerful supercomputers, targeting larger study areas rather than individual case studies. Furthermore, our method, specifically tailored for delineating IRES, could be valuable for generating the river network product across broader global regions, particularly in areas with low population density and difficult terrain for manual mapping. This would enhance our understanding of the temporal and spatial dynamics of river channels, particularly in headwater areas.

6. Conclusions

To accurately delineate IRES in headwater areas, we proposed a hybrid method that couples MERIT DEM datasets with surface water imagery. This method focuses on river networks typically found in headwater areas, characterized by sparse and disconnected water bodies, an abundance of narrow river channels, and indistinct channel features. Following the required preprocessing, this hybrid method consists of three main steps: (1) Deriving key points: The maximum water extent product and water indices data are organized into key points, representing the conjunction points in river channels, using the alpha shape, local maximum, and skeleton algorithms. (2) Building topology: Key points are organized into a topology matrix by applying the D8 algorithm to the MERIT Direction data to record the relative topological relationships between key point pairs. (3) Delineating river network: The slope and the enhanced water indices, processed by the gamma function, are integrated and regarded as the weighted graph. This graph is then processed by the A* algorithm to obtain the final optimal river channels in our study area. After developing the method, the derived river networks were compared with previous river networks and potential river networks from five types of DEM datasets. Completeness, as measured by positional accuracy, is our river network’s key performance metric. It surpasses others, achieving over 91% positional accuracy within a 30 m buffer. Consequently, this approach has the potential to generate IRES at larger scales, such as basin or national levels, and provide baselines for understanding surface water bodies in headwater areas.