A method for automatic mapping of seasonal rivers based on digital terrain constraints and satellite remote sensing

Abstract

The present invention discloses a method for automatically drawing seasonal rivers based on digital terrain constraints and satellite remote sensing. Based on the watershed network derived from the digital terrain model DEM and remote sensing images, the automatic extraction and drawing of seasonal rivers in a specified period of time is realized. This method uses the watershed network derived from the digital terrain model DEM as a spatial constraint condition, uses the optical remote sensing image water index method combined with adaptive threshold segmentation to realize the automatic water body extraction of each scene image within the study period, and overcomes the problem of monitoring difficulties caused by the dynamic and changeable seasonal rivers by calculating the frequency of water body occurrence within the period. Subsequently, the center line of the river mask is generated by object-based morphological operations, and the length, width, water area and other indicators of the seasonal river under various water body states are estimated. This method takes into account the continuous change state of the seasonal river within the time period, and can quickly and automatically draw the spatial distribution of the river, providing method support for understanding the spatiotemporal evolution characteristics of the river hydrological situation in areas with missing data.

Description

Technical Field

[0001] The invention belongs to the field of hydrological remote sensing, and in particular relates to a seasonal river automatic drawing method based on digital terrain constraints and satellite remote sensing. Background Art

[0002] Rivers originate from the complex interactions between various systems in the Earth's surface sphere. They play a key role in the Earth's water cycle by providing basic ecosystem services, influencing climate change and providing freshwater resources. In recent years, human activities and climate change have changed the global hydrological cycle, causing most rivers in the world to become fragile and unstable. These unstable rivers control material transport channels, affecting regional drainage patterns, hydrological cycle processes and ecological service functions.

[0003] Understanding the spatial pattern of rivers is an important basis for assessing river connectivity, revealing the hydrological cycle process and ecological and environmental impact. However, the spatial distribution of global rivers has not been fully recorded at this stage, especially seasonal rivers, whose dynamic and changeable characteristics make them extremely difficult to be finely monitored. Field measurement is the most accurate method to map river distribution and obtain river hydrological parameters, but large-scale and wide-area surveys cannot break through the limitations of economic and practical factors. The development of remote sensing earth observation technology has not only improved the spatial refinement of river observation, but also increased the frequency of river observation. Remote sensing mapping of river distribution is mainly achieved through two ways. One is remote sensing digital products derived from terrain, such as digital elevation models (DEM). DEM data provides relevant river information including potential catchment network and elevation, but the highest spatial resolution of DEM data publicly available in the world is only 30 meters. Its accuracy and spatial resolution limit the ability to extract narrow seasonal rivers, and cannot reflect the real water body flooding and time dynamic information. The other is to map the flooding range of rivers through remote sensing images, but it is difficult to automatically distinguish rivers from other types of water bodies in two-dimensional images. In addition, even if a multi-satellite collaborative observation method is adopted, the final result is the instantaneous state of the river at multiple times. It is still unable to break through the technical bottleneck of continuous and refined expression of dynamic and changeable seasonal rivers, and it is difficult to promote and use it on a large regional scale. Summary of the invention

[0004] The purpose of the present invention is to address the problems existing in the above-mentioned prior art and to provide a method for automatic drawing of seasonal rivers based on digital terrain constraints and satellite remote sensing. The method uses a watershed network derived from digital elevation model (DEM) data and satellite remote sensing images to extract the flooded range, reflect the continuous state of river flooding, and realize the refined drawing of dynamic and changeable seasonal rivers.

[0005] A seasonal river automatic drawing method based on digital terrain constraints and satellite remote sensing, comprising:

[0006] Obtain watershed network data and satellite remote sensing image data derived from DEM of the study area;

[0007] Generate buffer zones of preset distances on both sides of the catchment network line data to form potential flooding areas of seasonal rivers;

[0008] Extract the water area in each satellite remote sensing image, spatially aggregate the water areas extracted from all images, and calculate the frequency of water inundation in potential flooded areas within a period of time;

[0009] Obtain the water area range corresponding to different flooding conditions, including maximum, average and instantaneous flooding conditions;

[0010] Convert raster data of water areas with different flooding conditions into N Connect the spatial objects and obtain the theoretical river range through the intersection of the water body object and the line data;

[0011] Generate the centerline of the river channel through morphological operation on the theoretical river channel range;

[0012] The river centerline is used as a reference to the river waterside line, and the width of the river is estimated section by section through morphological distance transformation operations, completing the automatic drawing of the seasonal river.

[0013] As a preferred implementation, it also includes performing spatial intersection judgment on the buffer zone after the buffer zone is generated, and aggregating spatially overlapping components to form a potential flooding area of ​​the seasonal river.

[0014] As a preferred implementation, it also includes calculating the normalized water index after invalid pixels are removed from the satellite remote sensing image. Further, the invalid pixels include cloud, snow, and mountain shadow pixels; the cloud and snow pixels are obtained based on the cloud probability data of the remote sensing image at the same time, and the mountain shadow pixels are obtained based on the DEM data of the remote sensing image at the same time.

[0015] As a preferred embodiment, the water area is extracted in the following manner:

[0016] The normalized water index NDWI of each image is calculated, and the frequency distribution characteristics of each NDWI grayscale image are statistically analyzed. The segmentation threshold with the largest inter-class variance is searched using an equally spaced step size. The NDWI grayscale image is segmented based on the segmentation threshold to extract the water area.

[0017] As a preferred implementation, the frequency segmentation thresholds of water bodies in different flooding states are determined according to the sampling survey results, and the water area ranges corresponding to the different flooding states are obtained.

[0018] As a preferred implementation, the theoretical river range is obtained in the following manner:

[0019] The intersection between the water body object and the potential flooding area is judged. The object that meets the intersection condition is a high-confidence river channel object, otherwise it is an initial low-confidence river channel object;

[0020] Traverse the initial low-confidence river channel object and make an adjacent judgment with the high-confidence river channel object. If they are adjacent, the initial low-confidence river channel object is determined as a medium-confidence river channel object, otherwise it is a low-confidence river channel object;

[0021] The medium and high confidence river objects are subjected to conditional judgment, and the judgment condition is set as "the total number of pixels is greater than a " or "Maximum flooding frequency is greater than b %" or "The object's aspect ratio is greater than c ", if the conditions are met, it will be retained as a river channel object, otherwise it will be deleted; a , b , c is the preset threshold.

[0022] As a preferred implementation, it also includes performing a morphological closing operation on the retained river channel object to obtain a theoretical river channel range.

[0023] As a preferred implementation, the erosion operation and the morphological opening operation are cyclically performed on the theoretical river channel range until the river channel grid image is completely eroded to obtain the river channel center pixel.

[0024] As a preferred implementation, the estimating the width of the river section by section by using the morphological distance transformation operation includes:

[0025] The morphological distance transformation function is used to calculate the shortest Euclidean distance between all pixels in the river grid and the water edge, and the distance value corresponding to the center pixel of the river is extracted;

[0026] Traverse the center pixel of the river channel and count the number of river channel pixels in the 3*3 neighborhood of each pixel. If the number is ≤3, the river channel at that location is determined to be a single pixel width. Based on the remote sensing image resolution, the single pixel width attribute is assigned a value L Meters; if the number is greater than 3, then twice the Euclidean distance from the pixel to the water edge is taken as the river width

[0027] As a preferred implementation, the raster data of the water area in different flooding states is converted into an eight-connected spatial object, and the theoretical river range is obtained through the intersection of the water body object and the potential flooded area. Usually, four-connected or eight-connected spatial objects are used for water processing. In this application, the eight-connected processing has higher judgment accuracy than the four-connected one. For example, seasonal rivers with narrow widths are mostly represented by 1 to 2 pixels on a 10-meter resolution image, and there are many cases of diagonal connections between pixels. If four-connected is used, the two diagonally connected pixels will be broken up into independent objects, while the two diagonally connected pixels can still be connected as a whole when eight-connected is selected to avoid accidental deletion during judgment.

[0028] As a preferred implementation, the satellite remote sensing image data is Sentinel-2 satellite image data.

[0029] The present invention has the following beneficial effects:

[0030] (1) The seasonal river automatic mapping method proposed in the present invention greatly simplifies the workload of field measurements. It can be completed based on publicly available remote sensing digital products and images, and is particularly suitable for areas with insufficient data.

[0031] (2) The algorithm of the present invention is simple to implement and can quickly map the spatial distribution of rivers. It can be applied to the monitoring of dynamic changes of rivers on a large spatial scale. It has good applicability to rivers with different flow patterns and different geomorphic environments. It can be extended to large regional and even global scale research, providing methodological support for revealing the spatiotemporal evolution characteristics of regional river hydrological conditions under the background of global change. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 It is the sample area map and the DEM watershed constraint range of the embodiment of the present invention.

[0033] Figure 2 It is a flow chart of the method of the present invention.

[0034] Figure 3 It is a schematic diagram of the process of river channel extraction, centerline extraction, and river width estimation of the present invention.

[0035] Figure 4 This is the river centerline extraction result of the sample area in the embodiment of the present invention.

[0036] Figure 5 This is the verification result of the river width estimation in the embodiment of the present invention. DETAILED DESCRIPTION

[0037] The specific implementation of the present invention is further described in detail below in conjunction with the accompanying drawings and examples. The following implementation cases are used to illustrate the present invention, but are not used to limit the scope of the present invention.

[0038] The embodiment uses an inland area with an arid and semi-arid climate as a sample area to illustrate the technical solution of the present invention. The average altitude of the inland area is more than 4,000 meters. The plateau climate is cold and relatively dry. Most of the rivers in the inland area are seasonal rivers. Precipitation and glacier water are their main sources of supply. Most of them flow into lakes or dry up in deserts and salt marshes. Considering the characteristics of rivers in the inland area, as well as the challenges of difficulty in actual measurement and insufficient monitoring data in the inland area, it is representative as an embodiment of the present invention.

[0039] The main process of this embodiment is as follows Figure 2 As shown, the following steps are included:

[0040] Step 1: Preprocess the DEM-derived watershed network (Global River Topology (GRIT) product of GLO-30 DEM), and clip the watershed network data according to the sample area boundary of the embodiment. Generate a one-kilometer buffer on both sides of the watershed network line data, perform spatial intersection judgment on the generated buffers, and aggregate the spatially overlapping components to generate the potential flooding area of ​​the seasonal river. Spatial aggregation can be achieved using the dissolve operation of ArcGIS. Figure 1 It is a partial illustration of the sample area map and the DEM watershed constraint range.

[0041] Step 2: Obtain valid observation pixels of Sentinel-2 satellite remote sensing images. Obtain Sentinel-2 images, Sentinel-2 cloud probability data, and GLO-30 DEM data from July to October 2019, 2020, and 2021 from the Google Earth Engine (GEE) platform (https: / / earthengine.google.com / ). Pixels with Sentinel-2 cloud probability data raster values ​​greater than or equal to 20% and pixels with Sentinel-2 SCL band raster values ​​equal to 3 / 8 / 9 / 11 are determined as invalid observation pixels interfered by clouds, cloud shadows, and snow. Input GLO-30 DEM data, solar azimuth angle, solar zenith angle and other parameters of each image into the terrain analysis tool to predict invalid observation pixels blocked by mountain shadows. Eliminate invalid pixels blocked by clouds, cloud shadows, snow, and mountain shadows in Sentinel-2 satellite remote sensing images.

[0042] Step 3: Synthesize the flooding frequency. The normalized water index NDWI is calculated using the green band and near-infrared band within the effective observation range of each Sentinel-2 satellite image. The calculation formula is shown in (1). The frequency distribution characteristics of the water index grayscale image generated for each scene are statistically analyzed. The segmentation threshold with the largest inter-class variance is searched using an equally spaced step size. The water index NDWI grayscale image is segmented based on this threshold to extract the water area in the image. The water area extracted from all images is spatially aggregated to calculate the water flooding frequency of the potential flooded area within the time period (the proportion of the number of times an area is flooded by detected water bodies within a certain period of time to the total number of observations). The calculation formula is shown in (2).

[0043] NDVI= (Green-NIR) / (Green+NIR) (1)

[0044] Where Green represents the reflectivity of the green light band, and NIR represents the reflectivity of the near-infrared band.

[0045] Water Frequency=N water / N vo (2)

[0046] Where N water Indicates the total number of times the water body is detected to be flooded during the period, N vo Represents the total number of valid observations in the time period.

[0047] Step 4. Extract various water areas. Randomly select river sections within the sample area, analyze the temporal variation trend of the river width of the sampled river sections during the observation period, and deduce the optimal water body frequency segmentation threshold corresponding to the river width under different flooding conditions. Specifically: randomly select several river sections with a length of 100 meters, obtain the time series of the width of each river section (manually measured from July to October 2019-2021), and count the maximum and average values ​​of the time series; then adjust different water body frequency thresholds so that the maximum and average values ​​of the segmented water width correspond to the results obtained from the time series.

[0048] Based on the results of the river section time series analysis, in this embodiment, a 10% threshold is used to represent the maximum water range (maximum flooding state), and a 50% threshold is used to represent the average water range (average flooding state).

[0049] Step 5: River channel extraction. Convert the raster data of the water area expressing various flooding states (maximum, average, instantaneous) into eight-connected spatial objects (such as Figure 3In a), the water body object and the water collection line data are intersected and judged. The object that meets the intersection condition is temporarily determined as a high-confidence river channel object, otherwise it is temporarily determined as a low-confidence river channel object. Traverse the low-confidence river channel object to determine whether it is adjacent to the high-confidence river channel object. If it is adjacent, it is a medium-confidence river channel object, otherwise it is a low-confidence river channel object. The medium and high-confidence river channel objects are conditionally judged, and the judgment conditions are set as "the total number of pixels is greater than 3" or "the maximum flooding frequency is greater than 50%" or "the object aspect ratio is greater than 3". If the conditions are met, it is retained as a river channel object, otherwise it is deleted, and the river mask is obtained, such as Figure 3 b. And calculate the water area of ​​seasonal rivers as the attribute index of seasonal rivers.

[0050] Perform morphological closing operations (dilation operation first, then erosion operation) on the retained river channel objects to merge them into one object, and fill the center area of ​​the object less than or equal to 0.1 km 2 The theoretical river channel range is obtained by the holes, such as Figure 3 Middle c.

[0051] Step 6: Extract the center line of the river. Perform morphological corrosion operation on the raster image of the theoretical river range to obtain a refined river image. Perform morphological opening operation on the corroded raster image (corrosion operation first, then dilation operation), and mark the pixels that need to be deleted in the opening operation as the center pixels of the river. Repeat the corrosion operation and morphological opening operation until the raster image of the river is completely corroded. Vectorize the pixels marked as the center of the river and convert them into spatial point objects. Match the nearest points one by one based on the Euclidean distance less than 15 meters, and connect the paired points to obtain the center vector line of the river (such as Figure 3 (d) The seasonal river length is estimated as an attribute indicator of the seasonal river. Figure 4 This is a partial illustration of the results of extracting the center line of the sample area.

[0052] Step seven, river channel width estimation. First, the morphological distance transformation function is used to calculate the shortest Euclidean distance from all pixels in the river channel grid to the water's edge. The river channel center pixel layer is spatially superimposed to extract the distance value corresponding to the river channel center pixel. Then, the center pixels of the river channel are traversed one by one, and the number of river channel pixels in the 3*3 neighborhood of each pixel is counted. If the number is less than or equal to 3, it is determined that the river channel at that location is only a single pixel width. In this embodiment, a 10-meter resolution Sentinel-2 satellite remote sensing image is used, so the single pixel width attribute is assigned a value of 10 meters; if the number is greater than 3, then twice the Euclidean distance from the pixel to the water's edge is taken as the river channel width (such as Figure 3 Repeat the vectorization of the center pixel of the river channel and the connection operation of the paired points in step 6, and take the average of the widths of the paired points as the width value of the river section (such as Figure 3 f), as the attribute index of seasonal rivers. The verification results of river width estimation are shown in Figure 5 。

Claims

1. A method for automatic drawing of seasonal rivers based on digital terrain constraints and satellite remote sensing, characterized in that: include: Obtain watershed network data and satellite remote sensing image data derived from DEM of the study area; Generate buffer zones of preset distances on both sides of the catchment network line data to form potential flooding areas of seasonal rivers; Extract the water area in each satellite remote sensing image, spatially aggregate the water areas extracted from all images, and calculate the frequency of water inundation in potential flooded areas within a period of time; Obtain the water area range corresponding to different flooding conditions, including maximum, average and instantaneous flooding conditions; Convert raster data of water areas with different flooding conditions into N Connect the spatial objects and obtain the theoretical river range through the intersection of the water body object and the line data; The river centerline is generated by morphological operation on the theoretical river range; the theoretical river range is obtained as follows: The intersection between the water body object and the potential flooding area is judged. The object that meets the intersection condition is a high-confidence river channel object, otherwise it is an initial low-confidence river channel object; Traverse the initial low-confidence river channel object and make an adjacent judgment with the high-confidence river channel object. If they are adjacent, the initial low-confidence river channel object is determined as a medium-confidence river channel object, otherwise it is a low-confidence river channel object; The medium and high confidence river objects are conditionally judged, and the judgment condition is set as "the total number of pixels is greater than a "or" Maximum flooding frequency is greater than b %" or "The object's aspect ratio is greater than c ", if the conditions are met, it will be retained as a river channel object, otherwise it will be deleted; a , b , c is the preset threshold; By referring to the river waterline from the river centerline, the river width is estimated section by section through morphological distance transformation operation, and the seasonal river is automatically drawn, including: The morphological distance transformation function is used to calculate the shortest Euclidean distance between all pixels in the river grid and the water edge, and the distance value corresponding to the center pixel of the river is extracted; Traverse the center pixel of the river channel and count the number of river channel pixels in the 3*3 neighborhood of each pixel. If the number is ≤3, the river channel at that location is determined to be a single pixel width. Based on the remote sensing image resolution, the single pixel width attribute is assigned a value L If the number is greater than 3, then twice the Euclidean distance between the pixel and the water edge is taken as the river width.

2. The method according to claim 1, characterized in that It also includes, after the buffer zone is generated, making a spatial intersection judgment on the buffer zone, and aggregating spatially overlapping components to form a potential flooding area of ​​the seasonal river.

3. The method according to claim 1, characterized in that It also includes calculating the normalized water index after invalid pixels are removed from the satellite remote sensing image; the invalid pixels include cloud, snow, and mountain shadow pixels; the cloud and snow pixels are obtained based on the cloud probability data of the remote sensing image in the same period, and the mountain shadow pixels are obtained based on the DEM data of the remote sensing image in the same period.

4. The method according to claim 1, characterized in that: The extraction method of the water area is: The normalized water index NDWI of each image is calculated, and the frequency distribution characteristics of each NDWI grayscale image are statistically analyzed. The segmentation threshold with the largest inter-class variance is searched using an equally spaced step size. The NDWI grayscale image is segmented based on the segmentation threshold to extract the water area.

5. The method according to claim 1, characterized in that According to the sampling survey results, the frequency segmentation threshold of water bodies in different flooding states is determined to obtain the water area range corresponding to different flooding states.

6. The method according to claim 1, characterized in that It also includes performing morphological closing operations on the retained river channel objects to obtain the theoretical river channel range.

7. The method according to claim 1, characterized in that The erosion operation and morphological opening operation are cyclically performed on the theoretical river channel range until the river channel grid image is completely eroded to obtain the center pixel of the river channel.

8. The method according to claim 1, characterized in that The raster data of the water area range in different flooding states are converted into eight-connected spatial objects, and the theoretical river channel range is obtained through the intersection of the water body object and the potential flooding area.

Images