Automatic pasture land depression reservoir identification method based on digital elevation model

Abstract

The invention provides a pasture land depression reservoir automatic identification method based on a digital elevation model, and the method comprises the steps: collecting the topographic data of a target watershed through employing an unmanned aerial vehicle-mounted LiDAR technology, and generating an original digital elevation model DEM of the target watershed; topographic data noise is eliminated based on the DEM, and a smooth digital elevation model DEMF and potential depression position raster data DEMEX are generated; generating morphological reconstruction elevation data based on DEMF, and obtaining a residual digital elevation model CS3 of the target drainage basin; processing CS3 by using an MATLAB region connectivity analysis algorithm, and identifying a region connection body to generate a preliminary identification result CS4 of the earth embankment protruding out of the earth surface; a final earth embankment recognition result CS7 is obtained based on the CS4; and performing spatial overlay analysis on the CS7 and DEMEX, and outputting three-dimensional coordinates of the depression reservoir. The method can be used for quickly and automatically identifying the depression reservoir in the pasture land in the semi-arid region.

Description

Technical Field

[0001] The present invention relates to the field of soil and water conservation technology, in particular to a method for automatically identifying grassland depression reservoirs based on a digital elevation model. Background Art

[0002] Grassland is an important land resource in my country, mainly distributed in the semi-arid, arid and alpine areas in the west and north of my country. It is an important barrier area for maintaining national ecological security and a key node for the protection and restoration of ecosystems. At the same time, grassland provides basic production and living materials for herders, and is the main battlefield for developing modern animal husbandry and promoting sustainable development of the regional economy. However, in some semi-arid and arid areas, due to insufficient precipitation, some grasslands are facing the dilemma of water shortage and degradation of water conservation functions, causing grassland degradation and a decline in its ecological service functions. In order to make full use of precipitation resources and provide water for poultry and livestock, based on the regional terrain characteristics, herders spontaneously built simple earth embankments in some low-lying areas in the basin, and the depressions surrounded by the earth embankments were transformed into poultry and livestock water pools (such as Figure 2 Compared to specialized water storage and irrigation systems, depression reservoirs are constructed on-site at a lower cost. While storing rainwater, they also intercept eroded sediment from the slopes above the depression. This intercepted sediment not only reduces sediment transport within the watershed but, due to its high nutrient content, further enhances the local habitat for plant growth in the depression and promotes regional ecological function. Therefore, depression reservoirs play a vital role in the comprehensive utilization of grassland water resources and ecological restoration. However, due to the limited capacity of the depressions, their relatively simple overall structure, and the lack of ongoing maintenance and dredging, depression reservoirs and their surrounding earthen embankments are easily damaged by rainfall, runoff, and human and livestock disturbance. The sediment accumulated in the embankments and depressions surrounding the reservoirs can easily become a new source of sediment, exacerbating regional soil erosion and ecological degradation (Nichols et al., 2018; Nichols et al., 2021). These reservoirs are key targets for targeted soil erosion prevention and control and subsequent soil and water conservation design in the new era.

[0003] Accurately identifying the locations of depression-type reservoirs scattered across vast grasslands is a prerequisite for their subsequent management and maintenance, but it is also a highly challenging task. The difficulties lie in four main aspects: First, there was no unified technical standard for the construction of depression-type reservoirs, resulting in a certain degree of randomness in the height and shape of the earthen embankments surrounding the depressions. Second, the long-term operation of existing depression-type reservoirs has further altered the height and shape of the earthen embankments surrounding the depressions due to rainfall, runoff, and human and livestock interference, exacerbating the difficulty of identifying depression-type reservoirs through morphological parameters. Third, the location of depression-type reservoirs is often based on herders' subjective judgment of the basin's topography and the need for livestock and poultry water replenishment at different times, resulting in a high degree of randomness in their locations. Fourth, it is difficult to directly distinguish depression-type reservoirs from naturally occurring depressions within the basin; only depressions with artificially constructed earthen embankments are considered depression-type reservoirs.

[0004] Traditional methods for identifying the location of depression-type reservoirs primarily include field surveys and remote sensing visual interpretation. The former involves field investigations to record the location, shape, and operating conditions of depression-type reservoirs; the latter uses remote sensing data to manually delineate existing depression-type reservoirs and determine their locations. Both methods require significant human and material resources and have long identification cycles, making them difficult to effectively identify and accurately locate scattered depression-type reservoirs across large pastures. Furthermore, even with remote sensing interpretation, it can be difficult to objectively determine whether the identified depressions are naturally occurring or have been the product of human intervention. Therefore, a rapid and automated depression-type reservoir identification technology is urgently needed. Summary of the Invention

[0005] The purpose of the present invention is to provide a fast and automatic method for identifying water reservoirs in grassland depressions in semi-arid areas, so as to support sustainable utilization of grassland water resources and watershed management in semi-arid areas.

[0006] In order to achieve the above purpose, a technical solution is provided: a method for automatically identifying grassland depression reservoirs based on a digital elevation model, comprising the following steps.

[0007] S10. Use UAV-mounted LiDAR technology to collect topographic data of the target watershed and generate the original digital elevation model (DEM) of the target watershed;

[0008] S20, based on the original digital elevation model DEM of the target watershed, the terrain data noise is eliminated by the spatial filtering algorithm of the MATLAB platform to generate a smoothed digital elevation model DEM_F, and the potential depression location raster data DEM_EX of the target watershed is further generated by the depression filling algorithm of the MATLAB platform;

[0009] S30, based on the smoothed digital elevation model DEM_F, constructing a planar disk-shaped structural element and corroding the watershed elevation data to generate corroded elevation data, then using the watershed elevation data and the corroded elevation data to generate morphologically reconstructed elevation data, and further obtaining a residual digital elevation model CS3 of the target watershed;

[0010] S40, using the MATLAB regional connectivity analysis algorithm to process the residual digital elevation model CS3 of the target watershed, identifying regional connected bodies, and screening potential earth embankment areas in combination with a preset connected body length threshold, thereby generating a preliminary identification result CS4 of earth embankments protruding from the surface;

[0011] S50, based on the preliminary recognition result CS4 of the earth embankment protruding from the ground, using the MATLAB regional connectivity analysis algorithm, and setting the image boundary threshold, earth embankment length threshold, and earth embankment eccentricity threshold, to obtain the final earth embankment recognition result CS7;

[0012] S60 , spatially overlaying and analyzing the final identification result CS7 of the earth embankment and the grid data DEM_EX of the potential depression position in the region, and outputting the three-dimensional coordinates of the depression reservoir.

[0013] Furthermore, the step S20 specifically includes:

[0014] S21. Call the filter function in the MATLAB platform to perform a convolution operation on the original digital elevation model DEM through a 3-cell × 3-cell moving window to eliminate the high-frequency noise components in the terrain data and generate the target watershed smooth terrain data DEM_F with continuous and smooth characteristics;

[0015] S22, based on the target watershed smoothed terrain data DEM_F, call the fillsinks function to obtain the target watershed filled depression terrain data DEM_FF;

[0016] S23. Subtract the smoothed terrain data DEM_F of the target watershed from the filled terrain data DEM_FF of the target watershed to obtain the grid data DEM_EX of the potential depression position of the target watershed.

[0017] Furthermore, the step S30 specifically includes:

[0018] S31, using the MATLAB platform, storing the terrain elevation in the target watershed smooth terrain data DEM_F in matrix form, recorded as CS1;

[0019] S32, calling the strel function in the MATLAB platform to construct a planar disc-shaped structural element, denoted as SE;

[0020] S33, calling the imerode function in the MATLAB platform, using the planar disc-shaped structuring element SE to erode the elevation matrix CS1, and obtaining the eroded elevation data CSmaker;

[0021] S34, using the MATLAB platform, based on the eroded elevation data CSmaker and the elevation matrix CS1, calling the imreconstruct function to construct the morphologically reconstructed elevation data CS2;

[0022] S35. Using the MATLAB platform, subtract the morphologically reconstructed elevation data CS2 from the elevation data CS1 to obtain a residual digital elevation model CS3 of the target watershed.

[0023] Furthermore, the step S40 specifically includes:

[0024] S41, calling the bwconncomp function in the MATLAB platform to identify the connected area in the residual digital elevation model CS3 of the target watershed, which is recorded as CSC1;

[0025] S42, calling the regionprops function in the MATLAB platform, counting the major axis of the ellipse of each connected body in the connected region CSC1, numbering the connected bodies as 1, 2, 3, ... N, and recording the major axis of the ellipse corresponding to them as L1, L2, L3, ... L N , data storage, denoted as CSstats;

[0026] S43. Using the MATLAB platform, set the major axis threshold of the earth embankment ellipse, call the ismember function, exclude the connected bodies in the connected area CSC1 whose ellipse major axis is smaller than the set major axis threshold, store the remaining connected bodies, and obtain a preliminary recognition image CS4 of the earth embankment protruding from the surface.

[0027] Furthermore, the step S50 specifically includes:

[0028] S51, using the MATLAB platform, setting an image boundary threshold, excluding the connected bodies located at the periphery of the image in the preliminary identification result CS4 of the earth embankment protruding from the ground, storing the position information of the remaining connected bodies, and obtaining the earth embankment image CS5 of the identified earth embankment protruding from the ground;

[0029] S52. Using the MATLAB platform, based on the earth embankment image CS5 that identifies the protruding surface, call the imfill function to fill the connected body holes in the earth embankment image CS5 to obtain the earth embankment image CS6 that identifies the protruding surface;

[0030] S53, using the MATLAB platform, based on the earth embankment image CS6 where the protruding surface is identified, calling bwconncomp to calculate the connected area CSC2 in the earth embankment image CS6 where the protruding surface is identified;

[0031] S54, calling the regionprops function in the MATLAB platform, counting the ellipse major axis and ellipse eccentricity of each connected body in the connected region CSC2, numbering the connected bodies as 1, 2, 3, ... N, and recording the corresponding ellipse major axes as L1, L2, L3, ... L N , the corresponding ellipse eccentricity is F1, F2, F3, ... F N , data storage, recorded as CSstats2;

[0032] S55. Using the MATLAB platform, set the ellipse major axis threshold and the ellipse eccentricity threshold, call the ismember function, and exclude the connected bodies in the connected area CSC2 whose ellipse major axis is smaller than the set major axis threshold and whose ellipse eccentricity is smaller than the set eccentricity threshold. That is, exclude the straight-line earth embankments protruding from the ground surface, and store the remaining connected bodies to obtain the final identified earth embankment image CS7.

[0033] Furthermore, the step S60 specifically includes:

[0034] S61, calling the bwconncomp function in the MATLAB platform to identify the connected area in the earth embankment image CS7 protruding from the ground surface, which is recorded as CSC3;

[0035] S62. Call the labelmatrix function in the MATLAB platform to number each independent connected body in the connected area CSC3 as 1, 2, 3, ..., N, and the corresponding elevation is recorded as 1, 2, 3, ..., N. The data is stored as CS7_lable;

[0036] S63, replacing the elevation matrix in the original digital elevation model DEM with CS7_lable to obtain a new digital elevation model DEM_new;

[0037] S64, based on the newly created digital elevation model DEM_new, calling the GRIDobj2polygon function to obtain the boundary coordinates S of the earth embankment;

[0038] S65 , spatially overlay analysis is performed on the earth embankment boundary coordinates S and the potential depression location grid data DEM_EX in the region, and the three-dimensional coordinates of the depression reservoir are output.

[0039] The advantage of this method is that, compared to traditional manual methods of directly identifying water-retention depressions, this method utilizes the MATLAB platform to focus on the key characteristic of the elliptical shape of the earthen embankments surrounding water-retention depressions. This method first identifies easily recognizable earthen embankments protruding from the ground, and then determines whether the surrounding area is a water-retention depression by determining whether the embankment is linear or elliptical. By determining the presence of earthen embankments and their elliptical shape, it can directly identify depressions and reservoirs that have been manually interfered with, thus avoiding the possibility that the identified depressions are natural depressions rather than human-interferenced water-retention depressions. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 This is a flow chart of a method for automatically identifying grassland depression reservoirs based on a digital elevation model according to an embodiment of the present invention;

[0041] Figure 2 is a schematic diagram of a water storage depression to be identified within a watershed according to an embodiment of the present invention;

[0042] Figure 3 Schematic diagram of a smoothed digital elevation model DEM_F according to an embodiment of the present invention;

[0043] Figure 4 Schematic diagram of the residual digital elevation model CS3 of the target watershed according to an embodiment of the present invention;

[0044] Figure 5 is a schematic diagram of a preliminary recognition result CS4 according to an embodiment of the present invention;

[0045] Figure 6 1 is a schematic diagram of the final recognition result CS7 of the earth embankment according to an embodiment of the present invention;

[0046] Figure 7 It is a schematic diagram of outputting the coordinate position of the identified depression according to an embodiment of the present invention. DETAILED DESCRIPTION

[0047] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0048] The present invention utilizes the MATLAB platform and, based on terrain data, can identify grassland water-retention depressions in batches, quickly, with high precision, and automatically, and output their coordinate positions. During the identification process, there is no need to save intermediate calculation data, thereby saving hard disk storage space. The entry point of the present invention's technology is that simple earthen embankments are built around the water-retention depressions, and the earthen embankments surrounding the depressions have elliptical features. If earthen embankments are built around the depressions identified in the area, they must be water-retention depressions that have been artificially intervened, rather than natural depressions. The eccentricity of an object is the ratio of the focal length of the ellipse to the major axis of the ellipse after the object is approximated as an ellipse. The eccentricity is generally between 0 and 1. An object with an eccentricity of 0 can be considered a circular body, while an object with an eccentricity of 1 can be considered a straight line segment. Therefore, if the eccentricity of the earthen embankment found in a certain area is less than a certain eccentricity threshold, the earthen embankment has an elliptical feature, and the area surrounded by the earthen embankment is a water-retention depression. In other words, instead of directly looking for water-retention depressions, we first look for more easily identifiable earthen embankments protruding from the ground, and then determine whether the shape of the earthen embankments is linear or elliptical, so as to determine whether the area around the earthen embankments is a water-retention depression.

[0049] like Figure 1 As shown, the embodiment of the present invention takes a certain target small watershed as an example ( Figure 2 ), provides a method for automatically identifying grassland water storage depressions based on a digital elevation model, comprising the following steps:

[0050] S10, using drone-borne LiDAR technology to obtain original digital elevation model DEM data; the data acquisition in step S10 uses drone LiDAR technology, including flight route design, point cloud data collection, trajectory solution, point cloud fusion and coordinate conversion, point cloud filtering and construction of digital elevation model DEM and other steps.

[0051] S20. Based on the original digital elevation model DEM of the target watershed, the terrain data noise is eliminated through the spatial filtering algorithm of the MATLAB platform to generate a smoothed digital elevation model DEM_F. The potential depression location raster data DEM_EX of the target watershed is further generated through the depression filling algorithm of the MATLAB platform.

[0052] Step S20 specifically includes:

[0053] S21. Call the filter function in the MATLAB platform to perform convolution operation on the original digital elevation model DEM through a 3 cell × 3 cell moving window to eliminate the high-frequency noise components in the terrain data and generate the target watershed smooth terrain data DEM_F ( Figure 3 ).

[0054] S22, based on the target watershed smoothed terrain data DEM_F, call the fillsinks function to obtain the target watershed filled depression terrain data DEM_FF;

[0055] S23, subtracting the smoothed terrain data DEM_F of the target watershed from the filled depression terrain data DEM_FF of the target watershed to obtain the potential depression location raster data DEM_EX of the target watershed;

[0056] S31, using the MATLAB platform, storing the terrain elevation in the target watershed smooth terrain data DEM_F in matrix form, recorded as CS1;

[0057] S32. Call the strel function in the MATLAB platform to construct a planar disc-shaped structural element, denoted as SE;

[0058] S33, calling the imerode function in the MATLAB platform, using the planar disc-shaped structuring element SE to erode the elevation matrix CS1, and obtaining the eroded elevation data CSmaker;

[0059] S34, using the MATLAB platform, based on the eroded elevation data CSmaker and the elevation matrix CS1, calling the imreconstruct function to construct the morphologically reconstructed elevation data CS2;

[0060] S35, using the MATLAB platform, subtracting the morphologically reconstructed elevation data CS2 from the elevation data CS1 to obtain the residual digital elevation model CS3 of the target watershed ( Figure 4 ).

[0061] S40. Use the MATLAB regional connectivity analysis algorithm to process the residual digital elevation model CS3 of the target watershed, identify regional connected bodies, and use a preset connected body length threshold to screen potential earth embankment areas, generating a preliminary identification result CS4 of earth embankments protruding from the surface.

[0062] Step S40 specifically includes:

[0063] S41, calling the bwconncomp function in the MATLAB platform to identify the connected area in the residual digital elevation model CS3 of the target watershed, which is recorded as CSC1;

[0064] S42, calling the regionprops function in the MATLAB platform, counting the major axis of the ellipse of each connected body in the connected region CSC1, numbering the connected bodies as 1, 2, 3, ... N, and recording the major axis of the ellipse corresponding to them as L1, L2, L3, ... L N , data storage, denoted as CSstats;

[0065] S43, using the MATLAB platform, set the major axis threshold of the earth embankment ellipse, call the ismember function, exclude the connected bodies in the connected area CSC1 whose ellipse major axis is smaller than the set major axis threshold, store the remaining connected bodies, and obtain the preliminary recognition image CS4 of the earth embankment protruding from the surface ( Figure 5 ).

[0066] S50 , based on the preliminary recognition result CS4 of the earth embankment protruding from the ground, using the MATLAB regional connectivity analysis algorithm, and setting the image boundary threshold, earth embankment length threshold and earth embankment eccentricity threshold, to obtain the final earth embankment recognition result CS7.

[0067] Step S50 specifically includes:

[0068] S51, using the MATLAB platform, setting an image boundary threshold, excluding the connected bodies located at the periphery of the image in the preliminary identification result CS4 of the earth embankment protruding from the ground, storing the position information of the remaining connected bodies, and obtaining the earth embankment image CS5 of the identified earth embankment protruding from the ground;

[0069] S52. Using the MATLAB platform, based on the earth embankment image CS5 that identifies the protruding surface, call the imfill function to fill the connected body holes in the earth embankment image CS5 to obtain the earth embankment image CS6 that identifies the protruding surface;

[0070] S53, using the MATLAB platform, based on the earth embankment image CS6 where the protruding surface is identified, calling bwconncomp to calculate the connected area CSC2 in the earth embankment image CS6 where the protruding surface is identified;

[0071] S54, calling the regionprops function in the MATLAB platform, counting the ellipse major axis and ellipse eccentricity of each connected body in the connected region CSC2, numbering the connected bodies as 1, 2, 3, ... N, and recording the corresponding ellipse major axes as L1, L2, L3, ... L N , the corresponding ellipse eccentricity is F1, F2, F3, ... F N , data storage, recorded as CSstats2;

[0072] S55. Using the MATLAB platform, set the ellipse major axis threshold and the ellipse eccentricity threshold, call the ismember function, and exclude the connected bodies in the connected area CSC2 whose ellipse major axis is smaller than the set major axis threshold and whose ellipse eccentricity is smaller than the set eccentricity threshold. That is, exclude the straight-line earth embankments protruding from the ground surface, and store the remaining connected bodies to obtain the final earth embankment image CS7 ( Figure 6 ).

[0073] S60 , spatially overlaying and analyzing the final identification result CS7 of the earth embankment and the grid data DEM_EX of the potential depression position in the region, and outputting the three-dimensional coordinates of the depression reservoir.

[0074] Step S60 specifically includes:

[0075] S61, calling the bwconncomp function in the MATLAB platform to identify the connected area in the earth embankment image CS7 protruding from the ground surface, which is recorded as CSC3;

[0076] S62. Call the labelmatrix function in the MATLAB platform to number each independent connected body in the connected area CSC3 as 1, 2, 3, ..., N, and the corresponding elevation is recorded as 1, 2, 3, ..., N. The data is stored as CS7_lable;

[0077] S63, replacing the elevation matrix in the original digital elevation model DEM with CS7_lable to obtain a new digital elevation model DEM_new;

[0078] S64, based on the newly created digital elevation model DEM_new, calling the GRIDobj2polygon function to obtain the boundary coordinates S of the earth embankment;

[0079] S65, spatially overlaying and analyzing the earth embankment boundary coordinates S and the potential depression location grid data DEM_EX in the region, and outputting the three-dimensional coordinates of the water storage depression ( Figure 7 ).

[0080] The present invention can identify grassland water-storage depressions in batches, quickly, with high precision, and automatically, and output their coordinate positions. Taking the above case as an example, the present invention can accurately locate three water-storage depressions within a range of 2.4km×4.2km, requiring 15.4 seconds and achieving 100% recognition accuracy. If these three water-storage depressions are identified and their coordinates are output manually based on the ArcGIS platform, the time required is approximately 25 minutes; and during the recognition process, there is no need to save intermediate calculation data, thus saving hard disk storage space. Taking the above case as an example, if manual interpretation is used to find the depression, it is necessary to first locate the position of the identified depression, then construct a vector file for the position of the depression, then edit it, and then store it; however, the present invention can directly output the depression position in the last step, without the need to save any intermediate process data. The recognition results can support the sustainable use of grassland water resources and watershed management in semi-arid areas.

[0081] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in the present invention should be included in the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be based on the scope of protection of the claims.

Claims

1. A method for automatically identifying grassland depression reservoirs based on digital elevation models, characterized in that: The steps include: S10. Use UAV-mounted LiDAR technology to collect the topographic data of the target watershed and generate the original digital elevation model (DEM) of the target watershed; S20, based on the original digital elevation model DEM of the target watershed, the terrain data noise is eliminated by the spatial filtering algorithm of the MATLAB platform to generate a smoothed digital elevation model DEM_F, and the potential depression location raster data DEM_EX of the target watershed is further generated by the depression filling algorithm of the MATLAB platform; S30, based on the smoothed digital elevation model DEM_F, constructing a planar disk-shaped structural element and corroding the watershed elevation data to generate corroded elevation data, then using the watershed elevation data and the corroded elevation data to generate morphologically reconstructed elevation data, and further obtaining a residual digital elevation model CS3 of the target watershed; S40, using the MATLAB regional connectivity analysis algorithm to process the residual digital elevation model CS3 of the target watershed, identifying regional connected bodies, and screening potential earth embankment areas in combination with a preset connected body length threshold, thereby generating a preliminary identification result CS4 of earth embankments protruding from the surface; S50, based on the preliminary recognition result CS4 of the earth embankment protruding from the ground, using the MATLAB regional connectivity analysis algorithm, and setting the image boundary threshold, earth embankment length threshold, and earth embankment eccentricity threshold, to obtain the final earth embankment recognition result CS7; S60. Using the MATLAB platform, spatially overlay and analyze the final identification result CS7 of the earth embankment and the grid data DEM_EX of the potential depression position in the region, and output the three-dimensional coordinates of the water storage depression.

2. The method for automatically identifying grassland depression reservoirs based on a digital elevation model according to claim 1, wherein: The step S20 specifically includes: S21. Call the filter function in the MATLAB platform to perform a convolution operation on the original digital elevation model DEM through a 3-cell × 3-cell moving window to eliminate the high-frequency noise components in the terrain data and generate the target watershed smooth terrain data DEM_F with continuous and smooth characteristics; S22, based on the target watershed smoothed terrain data DEM_F, call the fillsinks function to obtain the target watershed filled depression terrain data DEM_FF; S23. Subtract the smoothed terrain data DEM_F of the target watershed from the filled terrain data DEM_FF of the target watershed to obtain the grid data DEM_EX of the potential depression position of the target watershed.

3. The method for automatically identifying grassland depression reservoirs based on a digital elevation model according to claim 1, wherein: The step S30 specifically includes: S31, using the MATLAB platform, storing the terrain elevation in the target watershed smooth terrain data DEM_F in matrix form, recorded as CS1; S32, calling the strel function in the MATLAB platform to construct a planar disc-shaped structural element, denoted as SE; S33, calling the imerode function in the MATLAB platform, using the planar disc-shaped structuring element SE to erode the elevation matrix CS1, and obtaining the eroded elevation data CSmaker; S34, using the MATLAB platform, based on the eroded elevation data CSmaker and the elevation matrix CS1, calling the imreconstruct function to construct the morphologically reconstructed elevation data CS2; S35. Using the MATLAB platform, subtract the morphologically reconstructed elevation data CS2 from the elevation data CS1 to obtain a residual digital elevation model CS3 of the target watershed.

4. The method for automatically identifying grassland depression reservoirs based on a digital elevation model according to claim 1, wherein: The step S40 specifically includes: S41, calling the bwconncomp function in the MATLAB platform to identify the connected area in the residual digital elevation model CS3 of the target watershed, which is recorded as CSC1; S42, calling the regionprops function in the MATLAB platform, counting the major axis of the ellipse of each connected body in the connected region CSC1, numbering the connected bodies as 1, 2, 3, ... N, and recording the major axis of the ellipse corresponding to them as L1, L2, L3, ... L N , data storage, denoted as CSstats; S43. Using the MATLAB platform, set the major axis threshold of the earth embankment ellipse, call the ismember function, exclude the connected bodies in the connected area CSC1 whose ellipse major axis is smaller than the set major axis threshold, store the remaining connected bodies, and obtain a preliminary recognition image CS4 of the earth embankment protruding from the surface.

5. The method for automatically identifying grassland depression reservoirs based on a digital elevation model according to claim 1, wherein: The step S50 specifically includes: S51, using the MATLAB platform, setting an image boundary threshold, excluding the connected bodies located at the periphery of the image in the preliminary identification result CS4 of the earth embankment protruding from the ground, storing the position information of the remaining connected bodies, and obtaining the earth embankment image CS5 of the identified earth embankment protruding from the ground; S52. Using the MATLAB platform, based on the earth embankment image CS5 that identifies the protruding surface, call the imfill function to fill the connected volume holes in the earth embankment image CS5 to obtain the earth embankment image CS6 that identifies the protruding surface; S53, using the MATLAB platform, based on the earth embankment image CS6 where the protruding surface is identified, calling bwconncomp to calculate the connected area CSC2 in the earth embankment image CS6 where the protruding surface is identified; S54, calling the regionprops function in the MATLAB platform, counting the ellipse major axis and ellipse eccentricity of each connected body in the connected region CSC2, numbering the connected bodies as 1, 2, 3, ... N, and recording the corresponding ellipse major axes as L1, L2, L3, ... L N , the corresponding ellipse eccentricity is F1, F2, F3, ... F N , data storage, recorded as CSstats2; S55. Using the MATLAB platform, set the ellipse major axis threshold and the ellipse eccentricity threshold, call the ismember function, and exclude the connected bodies in the connected area CSC2 whose ellipse major axis is smaller than the set major axis threshold and whose ellipse eccentricity is smaller than the set eccentricity threshold. That is, exclude the straight-line earth embankments protruding from the ground surface, and store the remaining connected bodies to obtain the final identified earth embankment image CS7.

6. The method for automatically identifying grassland depression reservoirs based on a digital elevation model according to claim 1, wherein: The step S60 specifically includes: S61, calling the bwconncomp function in the MATLAB platform to identify the connected area in the earth embankment image CS7 protruding from the ground surface, which is recorded as CSC3; S62. Call the labelmatrix function in the MATLAB platform to number each independent connected body in the connected area CSC3 as 1, 2, 3, ..., N, and the corresponding elevation is recorded as 1, 2, 3, ..., N. The data is stored as CS7_lable; S63, replacing the elevation matrix in the original digital elevation model DEM with CS7_lable to obtain a new digital elevation model DEM_new; S64. Based on the newly created digital elevation model DEM_new, call the GRIDobj2polygon function to obtain the boundary coordinates S of the earth embankment; S65 , spatially overlay analysis is performed on the earth embankment coordinate boundary S and the potential depression location grid data DEM_EX in the region, and the three-dimensional coordinates of the depression reservoir are output.

Images