Method for Extracting Small Water Systems in Mountainous Areas Based on Remote Sensing Images

Abstract

The present invention relates to the technical field of river information interpretation, and discloses a method for extracting small mountain water systems based on remote sensing images, aiming to solve the problem of insufficient accuracy in the existing extraction of small mountain water systems. The solution mainly includes: based on the time series of synthetic aperture radar images of satellites, using the random forest algorithm to preliminarily extract small mountain water systems; solving the correlation between the average slope and catchment area of the watershed based on the digital elevation model, and interpreting the catchment area at the first turning point thereof and taking it as the critical catchment area; using the critical catchment area to correct the preliminarily extracted results to obtain the final small mountain water systems of the watershed. The present invention improves the extraction accuracy of small mountain water systems, provides a scientific basis for simulating the river flow process of small mountain water systems, has the advantages of high calculation accuracy, wide application range and more objective and reasonable results, and is applicable to the extraction of small mountain water systems in mountainous watersheds.

Description

Technical Field

[0001] The invention relates to the technical field of river information interpretation, and in particular to a method for extracting small water systems in mountainous areas based on remote sensing images. Background Art

[0002] Small river systems are an important part of the river network. They help connect the lithosphere, hydrosphere and atmosphere, and play a key role in maintaining biodiversity, biogeochemical processes and the integrity of river system functions. However, due to the huge workload of field annotation, the distribution of small river systems at the basin scale is still unknown. Therefore, it is more feasible to extract small river systems at the basin scale using methods based on remote sensing images and digital elevation models.

[0003] Remote sensing images have excellent temporal and spatial resolution and are a reliable method for extracting river information. Among them, optical remote sensing images have been widely used in the monitoring of surface water systems. Among them, McFeeters et al. (1996) proposed a normalized difference water index based on the contrast between the visible light band and the near-infrared band to quickly extract water information from the image; Xu Hanqiu (2006) used the short-wave infrared band to replace the near-infrared band and proposed a modified normalized difference water index to better suppress vegetation information. Although the water index threshold method is an intuitive and concise water system extraction method, due to the significant differences in the spectral characteristics of water bodies in different regions, the water index threshold method only considers the spectral characteristics of water bodies, and the water body extraction accuracy is often poor. At the same time, optical remote sensing images are easily covered by clouds, and the surface cannot be effectively monitored in areas with dense clouds. River system extraction methods based on digital elevation models have also been widely used. Such methods often consider pixels that meet a certain threshold for characterizing river channel generation as rivers. For example, O'Callaghan and Mark (1984) used catchment area thresholds to determine the drainage system of a watershed, and Peckham (1995) used the Strahler series to determine the drainage system of a watershed. However, the river drainage system extraction algorithm based on the open source digital elevation model is limited by the coarse temporal and spatial resolution, making it difficult to extract small mountainous drainage systems, especially intermittent rivers.

[0004] Therefore, all-weather synthetic aperture radar images are widely used for river extraction in cloudy areas. For example, Klemenjak et al. (2012) combined support vector machine and maximum likelihood classification to extract rivers from synthetic aperture radar images; Zhu et al. (2014) made full use of the grayscale and morphological characteristics of water bodies in synthetic aperture radar images to extract rivers, and their results were significantly better than the traditional grayscale threshold method. Synthetic aperture radar images can penetrate clouds and continuously monitor targets, which makes it possible to extract temporally dynamic water bodies, but the inherent speckle noise problem of synthetic aperture radar images weakens its extraction accuracy for small water systems in mountainous areas. Summary of the invention

[0005] The invention aims to solve the problem of insufficient precision in the existing extraction of small water systems in mountainous areas, and proposes a method for extracting small water systems in mountainous areas based on remote sensing images.

[0006] The technical solution adopted by the present invention to solve the above technical problems is:

[0007] The method for extracting small water systems in mountainous areas based on remote sensing images includes the following steps:

[0008] Step 1, obtaining remote sensing images of a mountainous watershed, and performing multi-temporal filtering on the remote sensing images, wherein the remote sensing images include backscattering time series of satellite ascending and descending orbit images in units of years;

[0009] Step 2, obtaining four time series from the remote sensing image, and respectively calculating the first mathematical statistical features of the four time series, the four time series are the backscattering time series of the ascending orbit image from July to October each year, the backscattering time series of the ascending orbit image from December to March each year, the backscattering time series of the descending orbit image from July to October each year, and the backscattering time series of the descending orbit image from December to March each year;

[0010] Step 3: Draw a small water system based on the Google Earth remote sensing image, and obtain the second mathematical statistical features of the four time series corresponding to the Google Earth remote sensing image, and use the drawing results and their corresponding second mathematical statistical features as training samples to train the random forest model;

[0011] Step 4: input the first mathematical and statistical features into the trained random forest model for classification and identification, and obtain preliminary extraction results of small water systems in mountainous areas;

[0012] Step 5: Use the digital elevation model data to establish a correlation curve between the watershed area and the average slope value in the mountainous area, interpret the watershed area at the first turning point of the curve and use it as the critical watershed area for maintaining the river channel;

[0013] Step 6: Remove the parts of the preliminary extraction results of the small water system that do not meet the critical catchment area to obtain the final extraction results of the small water system in the mountainous area.

[0014] Furthermore, the step 1 specifically includes:

[0015] Step 11: Download radar satellite data products from the Geospatial Data Cloud website to obtain remote sensing images;

[0016] Step 12: Perform multi-temporal filtering on the remote sensing images corresponding to each year. The calculation formula is as follows:

[0017]

[0018] Where (x, y) are the pixel coordinates in the remote sensing image, and I i (x, y) is the pixel value at the (x, y) position in the i-th remote sensing image, N is the total number of remote sensing images, and J k (x, y) is the pixel value at the (x, y) position in the k-th remote sensing image after multi-temporal filtering, is the spatial average of all pixel values in the i-th remote sensing image, is the spatial average of all pixel values in the k-th remote sensing image, is the spatial average of all pixel values in the k-th remote sensing image after multi-temporal filtering. When performing multi-temporal filtering, it satisfies

[0019] Furthermore, in step 2, the first mathematical statistical feature and the second mathematical statistical feature include: median, range, standard deviation, coefficient of variation, kurtosis, and skewness.

[0020] Furthermore, the calculation formula for the first mathematical statistical feature is as follows:

[0021]

[0022]

[0023]

[0024]

[0025]

[0026]

[0027] Where Md t is the median of the t-th pixel, R t is the range of the t-th pixel, SD t is the standard deviation of the t-th pixel, CV t is the coefficient of variation of the t-th pixel, Ku t is the kurtosis of the t-th pixel, Sk t is the skewness value of the t-th pixel, is the value of the t-th pixel in the i-th remote sensing image, i = 1, 2,..., n, where n is the total number of remote sensing images, and sort is to sort the values of the t-th pixel in each remote sensing image from smallest to largest.

[0028] Furthermore, step 3 specifically includes:

[0029] Step 31: Determine the positions of small water systems and other ground objects based on remote sensing images of Google Earth, and convert the positions of small water systems and other ground objects into vector sample points via ArcGIS software. The classification labels corresponding to the vector sample points are 0 and 1 respectively, where 0 indicates the presence of small water systems and 1 indicates the absence of small water systems.

[0030] Step 32: Build a random forest model and K-fold cross-validation based on the machine learning library Scikit-learn on Python software.

[0031] Step 33: Convert the second mathematical statistical features and corresponding classification labels of the four time series corresponding to the vector sample points into a matrix in the form of [1, 25], and use the matrix as training samples to input into the random forest model for training. In the matrix, columns 1 to 24 are the second mathematical statistical features, and column 25 is the classification label.

[0032] Further, the specific steps of step 4 include:

[0033] Step 41: Convert the first mathematical statistical features corresponding to the mountainous watershed into a matrix in the form of [m, 24] and input it into the random forest model. The output format is a matrix in the form of [m, 1], where m is the total number of pixels of the complete remote sensing image, and the first column is the classification label predicted by the random forest model.

[0034] Step 42: Restore the output matrix to a matrix in the form of [i, j], where i and j are the number of rows and columns of the complete remote sensing image respectively.

[0035] Further, the specific steps of step 5 include:

[0036] Step 51: Use Python software to calculate the average slope value corresponding to each catchment area of the mountainous watershed

[0037] Step 52: Plot the relevant relationship curve between the average slope value corresponding to the mountainous watershed and the catchment area A using double logarithmic coordinates, visually interpret the catchment area corresponding to the first turning point in the relevant relationship curve, and take it as the critical catchment area. and the catchment area A are as follows:

[0038] Further, the average slope value corresponding to the mountainous watershed and the catchment area A data are as follows:

[0039] where j = 1, 2,..., M;

[0040] In the formula, is the j-th group of average slope value and catchment area, and M is the number of non-repeated catchment areas;

[0041] The calculation formula is as follows:

[0042]

[0043] In the formula, g is the number of slope value for the catchment area A j , S i is the i-th slope value for the catchment area A j .

[0044] Furthermore, step 6 specifically includes:

[0045] Step 61: Calculate the catchment area distribution of the mountainous watershed based on the digital elevation model, and generate the potential river channel distribution of the mountainous watershed based on the critical catchment area, where the positions satisfying the critical catchment area are marked as 0, and the remaining positions are marked as 1;

[0046] Step 62: Superimpose the preliminary small water system extraction result output by the random forest algorithm with the potential river channel distribution, and use the position distribution finally marked as 0 as the final extraction result of the small water system in the mountainous area.

[0047] The beneficial effects of the present invention are as follows: Based on the statistical characteristics of the satellite remote sensing image time series, the present invention captures small water systems that often exist in the form of mixed pixels, effectively improving the extraction accuracy of small water systems in mountainous areas. And based on the digital elevation model, the preliminary extraction results of small water systems that do not meet the minimum standard for maintaining river channels are eliminated, effectively reducing the interference of the speckle noise of the synthetic aperture radar image itself on the extraction results, improving the extraction accuracy of small water systems in mountainous areas, and providing a scientific basis for simulating the river flow process of small water systems in mountainous areas. In addition, all the data used in the present invention are official open-source data, with reliable data quality and low acquisition difficulty. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 is a schematic flow chart of the method for extracting small water systems in mountainous areas based on remote sensing images according to an embodiment of the present invention;

[0049] Figure 2 is a schematic flow chart of the method for obtaining the preliminary extraction result of small water systems in mountainous areas according to an embodiment of the present invention;

[0050] Figure 3 is a schematic diagram of the preliminary extraction result of small water systems in mountainous areas according to an embodiment of the present invention;

[0051] Figure 4 is a schematic diagram of the correlation curve between the average slope value and the catchment area according to an embodiment of the present invention;

[0052] Figure 5 is a schematic diagram of the final extraction result of small water systems in mountainous areas according to an embodiment of the present invention;

[0053] Figure 6 Schematic diagram for comparing the final extraction result of the small mountain water systems in the embodiments of the present invention with the high-precision satellite images of Google Earth Specific implementation manners

[0054] Taking the Duodigouzi watershed in the Lhasa River Basin of the Tibet Autonomous Region as the research area, the implementation manners of the present invention will be described in detail below. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention

[0055] The satellite used in this embodiment is the domestic high-resolution No. 3 SAR. The 1-meter resolution radar satellite data products can be downloaded from the Geospatial Data Cloud website, and the 30m resolution AW3D30 digital elevation model data can be obtained through Google Earth Engine

[0056] Please refer to Figure 1 and Figure 2 , the method for extracting small mountain water systems based on remote sensing images provided in this embodiment includes the following steps

[0057] Step 1: Obtain the remote sensing images of the mountain watershed and perform multi-temporal filtering on the remote sensing images. The remote sensing images include the backscattering time series of the ascending and descending orbit images of the high-resolution No. 3 satellite on an annual basis

[0058] In this embodiment, Step 1 specifically includes the following steps

[0059] Step 11: Download the radar satellite data products from the Geospatial Data Cloud website to obtain remote sensing images

[0060] In this embodiment, according to the boundary of the Duodigou watershed, a time series set of SAR image products of the corresponding area of the high-resolution No. 3 satellite in the ascending and descending orbits within 2020 is obtained from the Geospatial Data Cloud website, with a total of 61 images

[0061] Step 12: Perform multi-temporal filtering on the remote sensing images corresponding to each year

[0062] In this embodiment, Python software is used to perform multi-temporal filtering on the high-resolution No. 3 SAR image products within 2020. The specific calculation formula is as follows

[0063]

[0064] where (x, y) are the pixel coordinates in the remote sensing image, I i (x, y) is the pixel value at the position (x, y) in the i-th remote sensing image, N is the total number of remote sensing images, J k (x, y) is the pixel value at the position (x, y) in the k-th remote sensing image after multi-temporal filtering is the spatial average of all pixel values in the i-th remote sensing image, is the spatial average of all pixel values in the k-th remote sensing image, is the spatial average of all pixel values in the k-th remote sensing image after multi-temporal filtering, and when performing multi-temporal filtering, it satisfies

[0065] In this embodiment, based on the statistical characteristics of the time series of GF-3 SAR images with a resolution of 1m, small water systems that often exist in the form of mixed pixels are captured, effectively improving the accuracy of extracting small water systems in mountainous areas. And multi-temporal filtering is used to reduce the speckle noise in the original remote sensing image, further improving the accuracy of extracting small water systems in mountainous areas.

[0066] Step 2: Obtain four time series from the remote sensing images, and calculate the first mathematical statistical characteristics of the four time series respectively. The four time series are the backscattering time series of ascending orbit images from July to October each year, the backscattering time series of ascending orbit images from December to March each year, the backscattering time series of descending orbit images from July to October each year, and the backscattering time series of descending orbit images from December to March each year;

[0067] Specifically, by obtaining the backscattering time series of GF-3 satellite ascending and descending orbit SAR images from July to October and from December to March each year, and calculating the first mathematical statistical characteristics of these four time series respectively and using them as the independent variables of the random forest model.

[0068] In this embodiment, step 2 specifically includes the following steps:

[0069] Step 21: Divide the time series of GF-3 SAR images in 2020 into four segments according to the orbit and time period, namely ascending orbit (July to October), ascending orbit (December to March), descending orbit (July to October), and descending orbit (December to March) to capture small water systems;

[0070] Step 22: Calculate the first mathematical statistical characteristics of the image time series corresponding to each pixel in each time period.

[0071] In this embodiment, the first mathematical statistical characteristics include median, range, standard deviation, coefficient of variation, kurtosis, and skewness, and a total of 24 layers (4×6 = 24) of mathematical feature sets are generated. The calculation formulas of the first mathematical statistical characteristics are as follows:

[0072]

[0073]

[0074]

[0075]

[0076]

[0077]

[0078] In the formula, Md t is the median of the t-th pixel, R t is the range of the t-th pixel, SD t is the standard deviation of the t-th pixel, CV t is the coefficient of variation of the t-th pixel, Ku t is the kurtosis of the t-th pixel, Sk t is the skewness value of the t-th pixel, is the value of the t-th pixel in the i-th remote sensing image, i = 1, 2, …, n, where n is the total number of remote sensing images, and sort is to sort the values of the t-th pixel in each remote sensing image from smallest to largest.

[0079] Step 3: Draw small water systems based on Google Earth remote sensing images, and obtain the second mathematical statistical features of four time series corresponding to the Google Earth remote sensing images. Combine the drawing results and their corresponding second mathematical statistical features as training samples to train a random forest model;

[0080] In this embodiment, Step 3 specifically includes the following steps:

[0081] Step 31: Determine the positions of small water systems and other ground objects based on the remote sensing images of Google Earth, and convert the positions of small water systems and other ground objects into vector sample points via ArcGIS software. The classification labels corresponding to the vector sample points are 0 and 1 respectively, where 0 indicates the presence of small water systems and 1 indicates the absence of small water systems;

[0082] Specifically, determine the positions of small water systems and other ground objects (bare land, vegetation, artificial buildings, etc.) based on the high-precision satellite images of Google Earth. Then, input the marked KMZ format file into ArcGIS software to convert it into a vector sample point file, and set their corresponding labels to 1 and 0 respectively. A total of 6784 samples of small water systems and 6076 samples of other ground objects are generated;

[0083] Step 32: Build a random forest model and K-fold cross-validation based on the machine learning library Scikit-learn on Python software, and set the number of decision trees to 1000, the number of features to the square root of the total number of input variables, and the number of folds to 10;

[0084] Step 33: Convert the second mathematical statistical features of the four time series corresponding to the vector sample points and the corresponding classification labels into a matrix in the form of [1, 25], and use the matrix as a training sample to input into the random forest model for training. In the matrix, columns 1 to 24 are the second mathematical statistical features, and column 25 is the classification label.

[0085] Specifically, in this embodiment, Python software is used to read the geographical location information corresponding to the vector sample points and match it with the geographical location information of the pixels in the remote sensing image. The corresponding 24 groups of second mathematical statistical features and classification labels are converted into a matrix in the form of [1, 25] and input into the random forest model for model training. Among them, the first 24 columns are the mathematical feature sets, and the last 1 column is the classification label.

[0086] Step 4: Input the first mathematical statistical features without classification labels into the trained random forest model for classification and recognition to obtain a preliminary extraction result of the small mountainous water systems.

[0087] In this embodiment, Step 4 specifically includes the following steps:

[0088] Step 41: Convert the first mathematical statistical features corresponding to the mountainous watershed into a matrix in the form of [m, 24] and input it into the random forest model. The output format is a matrix in the form of [m, 1], where m is the total number of pixels in the complete remote sensing image, and the first column is the classification label predicted by the random forest model.

[0089] Specifically, based on the number of channels of the image set, the reshape function of Python software is used to convert the mathematical statistical feature set of the Duodi Gully watershed into a matrix in the form of [1085×1059, 24] and input it into the random forest model. The output format is a matrix in the form of [1085×1059, 1], where 1085×1059 = 1149015 is the total number of pixels in the complete remote sensing image, and the first column is the classification label predicted by the random forest model.

[0090] Step 42: Restore the output matrix to a matrix in the form of [i, j], where i and j are the number of rows and columns of the complete remote sensing image, respectively.

[0091] Please refer to Figure 3 , in this embodiment, according to the shape of the initial watershed image, the reshape function of Python software is used to restore the output matrix to a matrix in the form of [1085, 1059], where 1085 and 1059 are the number of rows and columns of the Duodi Gully image, respectively.

[0092] Step 5: Use the digital elevation model data to establish a correlation curve between the catchment area and the average slope value within the mountainous watershed, interpret the catchment area at the first turning point of the curve, and use it as the critical catchment area for maintaining the river channel.

[0093] In this embodiment, step 5 specifically includes the following steps:

[0094] Step 51: Use Python software to calculate the average slope value corresponding to each catchment area in the mountainous watershed

[0095] Specifically, ArcGIS software can be used to import AW3D30 digital elevation model data. After filling depressions, the D8 flow direction algorithm is used to calculate the water flow direction distribution in the Duodi Gully watershed and calculate the catchment area distribution of the watershed accordingly. Then, the slope tool of ArcGIS software is used to calculate the slope value of the raster cells.

[0096] Step 52: Plot the relevant relationship curve between the average slope value corresponding to the mountainous watershed and the catchment area A using double logarithmic coordinates, visually interpret the catchment area corresponding to the first turning point in the relevant relationship curve, and use it as the critical catchment area. Specifically, first, Python software can be used to calculate the average slope value corresponding to each catchment area data in the watershed and sort them to obtain the corresponding average slope value

[0097] and the catchment area A data:

[0098] where is the average slope value and catchment area of the jth group, j = 1, 2,..., M, and M is the number of non-repeated catchment areas;

[0099] Please refer to Figure 4 , the number of catchment area data selected in this embodiment for research is 1000, and the relevant relationship curve between the average slope value and the catchment area within 10 5 m 2 is plotted, which can effectively reduce the amount of calculation and ensure capturing the appearance of the first turning position of the relevant relationship curve.

[0100] Then, the relevant relationship curve between the average slope value and the catchment area A is plotted using double logarithmic coordinates, and the catchment area corresponding to its first turning point is visually interpreted and used as the critical catchment area of the watershed. Please refer to Figure 3 , in this embodiment, it is taken as 4500m 2 .

[0101] Step 6: Remove the part that does not meet the critical catchment area in the preliminary extraction result of the small water systems to obtain the final extraction result of the small water systems in the area.

[0102] In this embodiment, step 6 specifically includes the following steps:

[0103] Step 61: Calculate the catchment area distribution of the mountainous watershed based on the digital elevation model, and generate the potential river channel distribution of the mountainous watershed based on the critical catchment area, where the positions that meet the critical catchment area are marked as 0, and the remaining positions are marked as 1;

[0104] Specifically, the function of conditional raster in ArcGIS software can be used to analyze the catchment area distribution of the watershed, mark the positions that meet the critical catchment area as 0, which is the potential river channel distribution, and mark the remaining positions as 1.

[0105] Step 62: Overlay the preliminary small water system extraction result output by the random forest algorithm with the potential river channel distribution, and use the position distribution finally marked as 0 as the final small water system distribution of the mountainous watershed.

[0106] Please refer to Figure 5 , specifically, the preliminary small water system extraction result output by the random forest model can be overlaid with the potential river channel distribution through ArcGIS software. If the classification label in the preliminary small water system extraction result is 0 and the corresponding mark in the potential river channel distribution is 1, then after the overlay of the preliminary small water system extraction result and the potential river channel distribution, the final mark at the corresponding position is 1, thereby realizing the elimination of the part of the preliminary extraction result of the small water system that does not meet the critical catchment area. Only when the classification label in the preliminary small water system extraction result is 0 and the corresponding mark in the potential river channel distribution is 0, the final mark at the corresponding position is 0. Then, the final small water system distribution of the mountainous watershed can be obtained through the position distribution finally marked as 0.

[0107] Please refer to Figure 6 , by comparing the finally obtained small water system distribution with the high-precision Google Earth image of the corresponding watershed, it can be found that this embodiment effectively captures the distribution of small water systems.

[0108] In summary, the method for extracting small mountain water systems based on remote sensing images provided in this embodiment captures small water systems that often exist in the form of mixed pixels based on the statistical characteristics of the time series of GF-3 SAR satellite remote sensing images with a resolution of 1m, effectively improving the accuracy of small mountain water system extraction. And based on the digital elevation model, it eliminates the preliminary extraction results of small water systems that do not meet the minimum standard for maintaining river channels, effectively reducing the interference of the speckle noise of the synthetic aperture radar image itself on the extraction results, improving the extraction accuracy of small mountain water systems, and providing a scientific basis for simulating the river flow process of small mountain water systems. In addition, all the data applied in the present invention are official open-source data, with reliable data quality and low acquisition difficulty. The GF-3 synthetic aperture radar images and digital elevation model data can cover the whole world.

[0109] As described above, it is only the specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by those skilled in the art within the technical scope disclosed by the present invention should be covered by the protection scope of the present invention.

Claims

1. A method for extracting small water systems in mountainous areas based on remote sensing images, characterized in that, It includes the following steps: Step 1: Obtain remote sensing images of the mountainous watershed and perform multi-temporal filtering on the remote sensing images. The remote sensing images include the backscattering time series of ascending orbit images and descending orbit images of the satellite on an annual basis; Step 2: Obtain four time series from the remote sensing images and calculate the first mathematical statistical features of the four time series respectively. The four time series are the backscattering time series of ascending orbit images from July to October each year, the backscattering time series of ascending orbit images from December to March each year, the backscattering time series of descending orbit images from July to October each year, and the backscattering time series of descending orbit images from December to March each year; Step 3: Draw small water systems based on Google Earth remote sensing images and obtain the second mathematical statistical features of four time series corresponding to the Google Earth remote sensing images. Combine the drawing results and their corresponding second mathematical statistical features as training samples to train a random forest model; Step 4: Input the first mathematical statistical features into the trained random forest model for classification and recognition to obtain a preliminary extraction result of the small water systems in the mountainous area; Step 5: Use digital elevation model data to establish a correlation curve between the catchment area and the average slope value within the mountainous watershed, interpret the catchment area at the first turning point of the curve and use it as the critical catchment area for maintaining the river channel; Step 6: Remove the part that does not meet the critical catchment area from the preliminary extraction result of the small water systems to obtain the final extraction result of the small water systems in the mountainous area.

2. The method for extracting small mountain water systems based on remote sensing images according to claim 1, characterized in that, The specific content of Step 1 includes: Step 11: Download radar satellite data products from the Geospatial Data Cloud website to obtain remote sensing images; Step 12: Perform multi-temporal filtering on the remote sensing images corresponding to each year respectively. The calculation formula is as follows; where (x, y) are the pixel coordinates in the remote sensing image, and I i (x, y) is the pixel value at the (x, y) position in the i-th remote sensing image, N is the total number of remote sensing images, and J k (x, y) is the pixel value at the (x, y) position in the k-th remote sensing image after multi-temporal filtering, is the spatial average of all pixel values in the i-th remote sensing image, is the spatial average of all pixel values in the k-th remote sensing image, is the spatial average of all pixel values in the k-th remote sensing image after multi-temporal filtering. When performing multi-temporal filtering, 3. The method for extracting small mountain water systems based on remote sensing images according to claim 1, wherein, In Step 2, the first mathematical statistical features and the second mathematical statistical features include: median, range, standard deviation, coefficient of variation, kurtosis, and skewness.

4. The method for extracting small mountain water systems based on remote sensing images according to claim 3, characterized in that, The calculation formula of the first mathematical statistical features is as follows: Where Md t is the median of the t-th pixel, R t is the range of the t-th pixel, SD t is the standard deviation of the t-th pixel, CV t is the coefficient of variation of the t-th pixel, Ku t is the kurtosis of the t-th pixel, Sk t is the skewness value of the t-th pixel, is the value of the t-th pixel in the i-th remote sensing image, i = 1, 2, …, n, where n is the total number of remote sensing images, and sort is to sort the values of the t-th pixel in each remote sensing image in ascending order.

5. The method for extracting small mountain water systems based on remote sensing images according to claim 1, characterized in that, The specific content of Step 3 includes: Step 31: Determine the positions of small water systems and other ground objects based on Google Earth remote sensing images, and convert the positions of small water systems and other ground objects into vector sample points through ArcGIS software. The classification labels corresponding to the vector sample points are 0 and 1 respectively, where 0 indicates the presence of small water systems and 1 indicates the absence of small water systems; Step 32: Build a random forest model and K-fold cross-validation based on the machine learning library Scikit-learn on Python software; Step 33: Convert the second mathematical statistical features and the corresponding classification labels of the four time series corresponding to the vector sample points into a matrix in the form of [1,25], and input the matrix as a training sample into the random forest model for training. In the matrix, the first to 24 columns are the second mathematical statistical features, and the 25th column is the classification label.

6. The method for extracting small mountain water systems based on remote sensing images according to claim 1, wherein The specific content of Step 4 includes: Step 41: Convert the first mathematical statistical feature corresponding to the mountainous watershed into a matrix in the form of [m, 24] and input it into the random forest model. The output format is a matrix in the form of [m, 1], where m is the total number of pixels of the complete remote sensing image, and the first column is the classification label predicted by the random forest model; Step 42: Restore the output matrix to a matrix in the form of [i, j], where i and j are the number of rows and columns of the complete remote sensing image respectively.

7. The method for extracting small mountain water systems based on remote sensing images according to claim 1, characterized in that The specific steps of step 5 include: Step 51: Use Python software to calculate the average slope value corresponding to each catchment area in the mountainous watershed Step 52: Plot the correlation curve between the average slope value corresponding to the mountainous watershed and the catchment area A using double logarithmic coordinates, visually interpret the catchment area corresponding to the first turning point in the correlation curve, and take it as the critical catchment area. ​ 8. The method for extracting small mountain water systems based on remote sensing images according to claim 7, characterized in that, The average slope value corresponding to the mountainous watershed and the catchment area A data are as follows: where j = 1, 2, …, M; Wherein, is the average slope value and catchment area of the jth group, and M is the number of non-repeating catchment areas; The calculation formula is as follows: where, g is the number of slope values for the catchment area A j , S i is the i-th slope value for the catchment area A j .

9. The method for extracting small mountain water systems based on remote sensing images according to claim 1, characterized in that The specific steps of step 6 include: Step 61: Calculate the catchment area distribution of the mountainous watershed based on the digital elevation model, and generate the potential river channel distribution of the mountainous watershed based on the critical catchment area, where the positions satisfying the critical catchment area are marked as 0 and the remaining positions are marked as 1; Step 62: Superimpose the preliminary small water system extraction result output by the random forest algorithm on the potential river channel distribution, and use the final position distribution marked as 0 as the final extraction result of the mountainous small water system.

Images