A method for active and passive collaborative remote sensing extraction of surface water resources across climate zones

Abstract

The application relates to a kind of passive and active cooperative remote sensing extraction methods of cross-climate zone surface water resources, to satellite historical image as the basis preprocessing obtains the optical synthetic satellite image and radar synthetic satellite image of target research area, and trains to obtain optical image water body identification model, and contact the best water body classification threshold under radar synthetic satellite image, water body distribution identification is realized respectively through model identification and classification threshold, and finally the identification results of two kinds of ways are fused, obtain the water body distribution corresponding to each historical acquisition time under target research area respectively;Design scheme will optical and radar data depth fusion, realizes large-scale, long time sequence fast automation water body mapping, significantly improves the precision and reliability of extraction result, not only can make up the limitation of single data source, also can provide important support for the dynamic monitoring and scientific management of watershed water resources, especially suitable for the water resource monitoring of cross-climate zone complex watershed.

Description

TECHNICAL FIELD

[0001] The present application relates to a kind of passive and active cooperative remote sensing extraction method of surface water resources across climate zone, belong to surface water body detection technical field. BACKGROUND

[0002] Surface water resources are an important part of global ecosystems, with significant environmental regulation functions and ecological benefits, closely related to human health, living environment, agricultural production and economic activities. In the Nile River Basin in Africa, an important region spanning tropical, subtropical and arid climate zones, surface water resources are not only a key factor in maintaining ecological balance and social development, but also an important foundation for coordinated development among countries in the basin. With the intensification of global climate change and human activities, the distribution and dynamic changes of water resources in this basin have increasingly attracted attention.

[0003] Efficient detection of the spatiotemporal distribution characteristics of surface water bodies and dynamic monitoring of their change rules are of great significance in complex regions like the Nile River Basin across climate zones. On the one hand, such monitoring can provide data support for the scientific management of water resources and help stabilize and restore regional ecosystems; on the other hand, it can provide scientific basis for climate change impact assessment and transnational water resources management strategies.

[0004] Remote sensing technology has been widely used in the extraction and monitoring of surface water resources due to its real-time, periodic, large-scale coverage and multi-source data fusion advantages. Optical remote sensing technology, such as Sentinel-2 imagery, excels in low, medium and high resolution surface water body extraction with its rich spectral information, but it is easily affected by cloud shadows, complex terrain and snow and other interference environments. Radar remote sensing technology, such as Sentinel-1 imagery, can penetrate clouds and has all-weather imaging capability, making it suitable for surface water extraction in complex weather conditions. However, radar imagery is heavily affected by speckle noise, and its image grayscale is not as continuous as optical imagery. In addition, building and mountain shadows can also cause some misjudgments in radar data scattering characteristics.

[0005] In recent years, the method of active and passive remote sensing data cooperation has gradually attracted attention. Some researchers try to combine the complementary advantages of optical and radar remote sensing data to improve the accuracy and reliability of surface water extraction. For example, Li et al. (2022) used the random forest (RF) algorithm to achieve high-precision classification of surface water in Sentinel-2 images, with producer's accuracy, user's accuracy and F1 value reaching 91.208%, 88.706% and 89.940%, respectively. Bioresita et al. (2019) combined Sentinel-1 and Sentinel-2 images and used decision-level data fusion technology to study the extraction method of permanent and temporary water bodies in the central region of Ireland, achieving good results. Cao et al. (2019) used the local histogram characteristics of SAR images to suppress the influence of speckle noise on water extraction results through a double-peak method. Saghafi et al. (2021) used multiple features of Landsat-8, Sentinel-1 and Sentinel-2 images to train multiple classifiers and perform decision-level fusion to improve the recognition accuracy of surface water.

[0006] Although the above research has made significant progress in the application of remote sensing technology in surface water extraction, there are still the following limitations: (1) Many methods rely on a large number of manually drawn samples for model training, while the technology for automatically obtaining high-quality samples still needs to be further optimized; (2) The collaborative application of optical and radar data is still not deep enough, and most studies only use a single data source or do not fully exploit the fusion potential of multi-source data; (3) Water extraction research is mostly focused on a single time period, lacking long-term sequence data support, making it difficult to accurately reflect the spatiotemporal dynamic changes of water bodies in cross-climate zones, which is not conducive to in-depth analysis of the comprehensive impact of climate and environmental change. SUMMARY

[0007] The technical problem to be solved by the present application is to provide a cross-climate zone surface water resource active and passive collaborative remote sensing extraction method, which deeply fuses optical and radar data, combines automated sample generation technology and multi-temporal data processing, and accurately performs temporal monitoring and analysis of surface water resources.

[0008] The present application adopts the following technical solution to solve the above technical problems: The present application designs a cross-climate zone surface water resource active and passive collaborative remote sensing extraction method for extracting the distribution of water bodies in the target study area, comprising the following steps:

[0009] Step A. Obtain optical satellite historical images and radar satellite historical images corresponding to each preset historical acquisition time in the target study area within the preset historical period, and obtain optical composite satellite images and radar composite satellite images corresponding to the target study area through preprocessing, and then enter step B;

[0010] Step B. Based on the optical synthetic satellite images corresponding to the target study area, and the optical synthetic satellite images respectively corresponding to each preset first water body index layer and each preset vegetation index layer, the classification labels of each pixel position in the satellite images about the preset water body classification and each non-water body classification are obtained, and the target classification model is trained to obtain an optical image water body recognition model; then step C is entered;

[0011] Step C. According to the classification labels of each pixel position about the preset categories, for the radar synthetic satellite images corresponding to the target study area, and the preset second water body index values respectively corresponding to each pixel position in the radar synthetic satellite images, the best water body classification threshold is obtained, and then step D is entered;

[0012] Step D. For each local image of each pixel position in each optical satellite historical image, the optical image water body recognition model is applied for processing to obtain the classification labels of each pixel position in each optical satellite historical image, and then the water body distribution in each optical satellite historical image is obtained;

[0013] At the same time, for each pixel position in each radar satellite historical image, the best water body classification threshold is compared and judged to determine whether each pixel position in each radar satellite historical image is a water body area, and then the water body distribution in each radar satellite historical image is obtained;

[0014] Then step E is entered;

[0015] Step E. For each historical acquisition time, each intersection region between the water body distribution in the optical satellite historical image and the water body distribution in the radar satellite historical image at the historical acquisition time is obtained, and the water body distribution of the target study area at the historical acquisition time is formed by each intersection region, and then the water body distribution of the target study area corresponding to each historical acquisition time is obtained.

[0016] As a preferred technical solution of the present application: in step A, for the obtained optical satellite historical images and radar satellite historical images corresponding to the target study area respectively corresponding to each preset historical acquisition time in the preset historical period, steps A1-1 to A1-3 and steps A2-1 to A2-2 are simultaneously performed:

[0017] Step A1-1. For each optical satellite historical image, the optical satellite historical images in which the area ratio of the cloud area to the whole image area is less than a preset ratio threshold are obtained as each optical primary satellite image, and then step A1-2 is entered;

[0018] Step A1-2. For each optical primary satellite image, remove the cloud layer local area covering at least the preset lower limit number of pixels in the optical primary satellite image, update the optical primary satellite image, and then enter step A1-3;

[0019] Step A1-3. For each pixel position in the satellite image, obtain the median value of the pixel value corresponding to each optical primary satellite image at the pixel position, to form the pixel median value corresponding to the pixel position; and then combine the pixel median values corresponding to each pixel position to form the optical synthetic satellite image.

[0020] Step A2-1. For each radar satellite historical image, first perform boundary noise removal processing to update, and then use the Refine Lee filtering algorithm to remove the scattering noise in the radar satellite historical image to form a radar primary satellite image, and then enter step A2-2.

[0021] Step A2-2. For each pixel position in the satellite image, obtain the median value of the pixel value corresponding to each radar primary satellite image at the pixel position, to form the pixel median value corresponding to the pixel position; and then combine the pixel median values corresponding to each pixel position to form the radar synthetic satellite image.

[0022] As a preferred technical solution of the present application: the step B includes the following steps B1 to B3.

[0023] Step B1. Obtain the preset first water body index values corresponding to each pixel position in the optical synthetic satellite image of the target research area, including the normalized water body index value NDWI, the normalized difference water body index value MNDWI, the land water body extraction index value LSWI, and the AWEIsh value and AWEInsh value in the automatic water body extraction index, to form the first water body index layer corresponding to the optical synthetic satellite image, including the normalized water body index layer, the normalized difference water body index layer, the land water body extraction index layer, and the automatic water body extraction index layer.

[0024] At the same time, obtain the preset vegetation index values corresponding to each pixel position in the optical synthetic satellite image of the target research area, including the normalized vegetation index value NDVI and the vegetation enhancement index value EVI, to form the vegetation index layer corresponding to the optical synthetic satellite image, including the normalized vegetation index layer and the vegetation enhancement index layer.

[0025] Then enter step B2.

[0026] Step B2. According to the optical synthetic satellite image, and the preset each first water body index layer and each vegetation index layer corresponding to the optical synthetic satellite image, combined with the dynamic world land classification data set, the classification labels of each pixel position in the satellite image about the preset water body classification and each non-water body classification are determined, and a preset number of pixel positions under each classification label are randomly selected to form each sample pixel position. The local image of the optical synthetic satellite image corresponding to the sample pixel position is used to form a sample in combination with the classification label of the sample pixel position, and then a sample data set is formed by each sample, and then step B3 is entered;

[0027] Step B3. Based on the sample data set, the local image in the sample is used as the input, and the classification label in the sample is used as the output. The target classification model is trained to obtain an optical image water body recognition model.

[0028] As a preferred technical solution of the present application: in step B2, for each sample pixel position obtained by random selection, the maximum inundation range layer in the JRC global land surface water dynamic data set is applied to update the mask processing of the local image of the optical synthetic satellite image corresponding to each sample pixel position corresponding to the water body classification label, and then the local image of the optical synthetic satellite image corresponding to the sample pixel position is used to form a sample in combination with the classification label of the sample pixel position.

[0029] As a preferred technical solution of the present application: the non-water body classification includes vegetation classification label, bare land classification label, cloud classification label and snow classification label.

[0030] As a preferred technical solution of the present application: step C includes the following steps C1 to C3 to obtain the best water body classification threshold value.

[0031] Step C1. Obtain a preset second water body index value of each pixel position in the radar synthetic satellite image corresponding to the target research area corresponding to the dual-polarization radar water body index value SDWI, and then enter step C2.

[0032] Step C2. According to the classification label of each pixel position about the preset each category, the dual-polarization radar water body index value SDWI of each pixel position in the radar synthetic satellite image is statistically determined to determine the water body initial classification threshold T1 of the dual-polarization radar water body index value SDWI of each pixel position in the radar synthetic satellite image for distinguishing the water body classification from each non-water body classification.

[0033] Meanwhile, statistics are made on the vertical polarization backscattering coefficient value VV involved in the calculation of the dual-polarized radar water body index value SDWI of each pixel position in the radar synthetic satellite image, to determine the water body initial classification threshold T2 of the vertical polarization backscattering coefficient value VV of the pixel position in the radar synthetic satellite image for distinguishing the water body classification from each non-water body classification;

[0034] Then, step C3 is entered;

[0035] Step C3. Based on the water body initial classification threshold T1 and the water body initial classification threshold T2 respectively, the optimal SDWI water body classification threshold of the dual-polarized radar water body index value SDWI of the pixel position in the radar synthetic satellite image for distinguishing the water body classification from each non-water body classification and the optimal VV water body classification threshold of the vertical polarization backscattering coefficient value VV of the pixel position in the radar synthetic satellite image for distinguishing the water body classification from each non-water body classification are determined according to the water body edge detection of the radar synthetic satellite image in combination with the Otsu algorithm.

[0036] As a preferred technical solution of the present application: in step C3, based on the water body initial classification threshold T1, steps C3a-1 to C3a-3 are performed as follows:

[0037] Step C3a-1. According to the water body initial classification threshold T1, the Canny edge detection operator is applied to determine the edge information of each water body in the radar synthetic satellite image, and each water body with an edge perimeter greater than a preset perimeter threshold is obtained as each water body to be analyzed, and then step C3a-2 is entered;

[0038] Step C3a-2. For each water body to be analyzed, the edge of the water body to be analyzed is extended in the water direction and the water body outside direction by a preset distance respectively to form an edge buffer zone area, and then the edge buffer zone area corresponding to each water body to be analyzed is obtained, and then step C3a-3 is entered;

[0039] Step C3a-3. The dual-polarized radar water body index value SDWI of each pixel position in the edge buffer zone area corresponding to each water body to be analyzed is statistically analyzed, and the Otsu algorithm is applied to determine the optimal SDWI water body classification threshold of the dual-polarized radar water body index value SDWI of the pixel position in the radar synthetic satellite image for distinguishing the water body classification from each non-water body classification by maximizing the inter-class variance;

[0040] Meanwhile, based on the water body initial classification threshold T2, steps C3b-1 to C3b-3 are performed as follows:

[0041] Step C3b-1. According to the initial classification threshold T2 of the water body, the Canny edge detection operator is applied to determine the edge information of each water body in the radar synthetic satellite image, and each water body with an edge perimeter greater than a preset perimeter threshold is obtained as each water body to be analyzed, and then step C3b-2 is entered.

[0042] Step C3b-2. For each water body to be analyzed, the edge buffer region is formed by extending the edge of the water body to be analyzed in the water direction and the water outside direction by a preset distance, and then the edge buffer region corresponding to each water body to be analyzed is obtained, and then step C3b-3 is entered.

[0043] Step C3b-3. The vertical polarization backscatter coefficient value VV of each pixel position in the edge buffer region corresponding to each water body to be analyzed is statistically analyzed, and the Otsu algorithm is applied to maximize the inter-class variance to determine the optimal VV water body classification threshold of the vertical polarization backscatter coefficient value VV of the pixel position in the radar synthetic satellite image for distinguishing water body classification from each non-water body classification.

[0044] As a preferred technical solution of the present application: in step D, for each pixel position in each radar satellite historical image, it is determined whether the pixel position is a water body region or a non-water body region by judging whether the dual-polarization radar water body index value SDWI of the pixel position is greater than the optimal SDWI water body classification threshold and the vertical polarization backscatter coefficient value VV of the pixel position is less than the optimal VV water body classification threshold.

[0045] As a preferred technical solution of the present application: it further includes the following step F, which is entered after step E is executed;

[0046] Step F. For each union region between the water body distribution in the optical satellite historical image and the water body distribution in the radar satellite historical image at each historical acquisition time, the following operations are performed to update the union region and further update the water body distribution of the target research region corresponding to each historical acquisition time;

[0047] First, according to the region with a slope greater than a preset slope threshold as a mountain shadow area, the union region is masked and processed to remove the mountain shadow area in the union region, and the union region is updated. Then, according to the JRC global land surface water data set global human settlement area layer, the union region is subjected to secondary mask processing, and the building shadow area and the misextracted building area in the union region are removed, and the union region is updated.

[0048] Compared with the prior art, the cross-climate zone surface water resource active and passive cooperative remote sensing extraction method has the following technical effects:

[0049] (1) The present application designs a cross-climate zone surface water resource active and passive collaborative remote sensing extraction method. First, based on satellite historical images, the optical synthetic satellite images and radar synthetic satellite images of the target research area are obtained by preprocessing, then the optical image water body recognition model is trained, and the best water body classification threshold is obtained by contacting the radar synthetic satellite images, the water body distribution is recognized by the model recognition and classification threshold, and finally the recognition results of the two ways are fused to obtain the water body distribution of the target research area corresponding to each historical acquisition time. The design scheme deeply fuses optical and radar data, realizes large-scale, long-time sequence rapid automation water body mapping, significantly improves the accuracy and reliability of the extraction result, not only makes up for the limitations of single data source, but also provides important support for dynamic monitoring and scientific management of watershed water resources, especially suitable for water resource monitoring of complex cross-climate zone basins, and provides important data support and technical support for water resource management, ecological protection and hydrological dynamic change analysis of the African Nile River basin in application. BRIEF DESCRIPTION OF DRAWINGS

[0050] Figure 1 is a schematic diagram of a cross-climate zone surface water resource active and passive collaborative remote sensing extraction method designed by the present application.

[0051] Figure 2 is the RGB composite image of the median composite image Sentinel 2 data of the Blue Nile River basin in 2021 obtained in the embodiment of the present application.

[0052] Figure 3 is a water body extraction flowchart based on Sentinel 1 remote sensing image in the embodiment of the present application.

[0053] Figure 4 is the monthly water body extraction result of the Blue Nile River basin based on optical and radar data collaboration realized by the method. DETAILED DESCRIPTION

[0054] The specific embodiments of the present application will be further described in detail below in combination with the drawings of the specification.

[0055] The present application designs a cross-climate zone surface water resource active and passive collaborative remote sensing extraction method, which is used to extract the water body area distribution in the target research area. In actual application, as shown in Figure 1 , the specific design and execution steps are as follows.

[0056] Step A. Obtain optical satellite historical images and radar satellite historical images corresponding to each preset historical acquisition time in the target research area within the preset historical period, and obtain optical synthetic satellite images and radar synthetic satellite images corresponding to the target research area by preprocessing, then enter step B.

[0057] The above step A in practical application, for the target research area respectively corresponding to each preset historical acquisition time in the preset historical period of optical satellite historical image, step A1-1 to step A1-3 are executed, and the optical synthetic satellite image corresponding to the target research area is obtained.

[0058] Step A1-1. For each optical satellite historical image, the cloud area area ratio of each optical satellite historical image is obtained, which is less than 20% of the preset ratio threshold, as each optical primary satellite image, and then step A1-2 is entered.

[0059] Step A1-2. For each optical primary satellite image, the cloud local area covering at least the preset lower limit pixel number is removed, the optical primary satellite image is updated, and then step A1-3 is entered.

[0060] Step A1-3. For each pixel position in the satellite image, the median value of the pixel value corresponding to each optical primary satellite image is obtained, and the pixel median value corresponding to the pixel position is constructed; then the pixel median value corresponding to each pixel position is combined into an optical synthetic satellite image.

[0061] At the same time, for the target research area respectively corresponding to each preset historical acquisition time in the preset historical period of radar satellite historical image, step A2-1 to step A2-2 are executed, and the radar synthetic satellite image corresponding to the target research area is obtained.

[0062] Step A2-1. For each radar satellite historical image, boundary noise removal processing is first executed to update, and then Refine Lee filtering algorithm is used to remove scattering noise in the radar satellite historical image to construct a radar primary satellite image, and then step A2-2 is entered.

[0063] Step A2-2. For each pixel position in the satellite image, the median value of the pixel value corresponding to each radar primary satellite image is obtained, and the pixel median value corresponding to the pixel position is constructed; then the pixel median value corresponding to each pixel position is combined into a radar synthetic satellite image.

[0064] Step B. Based on the optical synthetic satellite image corresponding to the target research area, and the optical synthetic satellite image corresponding to each preset first water body index layer and each preset vegetation index layer, the classification label of each pixel position in the satellite image about the preset water body classification and each non-water body classification is obtained, and the target classification model is trained to obtain the optical image water body recognition model; then step C is entered.

[0065] In practical applications, the above step B is specifically designed to perform the following steps B1 to B3.

[0066] Step B1. The following calculation formula is used:

[0067]

[0068]

[0069] AWEIsh = BLUE + 2.5 x GREEN - 1.5 x (NIR + SWIR1) - 0.25 x SWIR2

[0070] AWEInsh = 4 x (GREEN - SWIR1) - (0.25 x NIR + 2.75 x SWIR2)

[0071] The target research area corresponds to the optical synthetic satellite image, and each pixel position corresponds to a preset first water body index value, including the normalized water body index value NDWI, the normalized difference water body index value MNDWI, the land water body extraction index value LSWI, and the AWEIsh value and AWEInsh value in the automatic water body extraction index, which constitutes the first water body index layer corresponding to the optical synthetic satellite image, including the normalized water body index layer, the normalized difference water body index layer, the land water body extraction index layer, and the automatic water body extraction index layer. Wherein, GREEN represents the pixel value of the green light band corresponding to the pixel position in the optical synthetic satellite image, NIR represents the pixel value of the near-infrared band corresponding to the pixel position in the optical synthetic satellite image, BLUE represents the pixel value of the blue light band corresponding to the pixel position in the optical synthetic satellite image, SWIR1 and SWIR2 respectively represent the pixel value of the first short-wave infrared band corresponding to the pixel position in the optical synthetic satellite image and the pixel value of the second short-wave infrared band.

[0072] In practical applications, the normalized water body index value NDWI can enhance the water body in the satellite image, and is mainly used for detecting and monitoring the small changes in water content. The disadvantage of this index is that it is sensitive to building structures and cannot effectively suppress building noise and mountain shadows in mountainous areas, which may lead to overestimation of water bodies; the normalized difference water body index value MNDWI can effectively suppress background noise such as buildings and bare land, and has higher precision in extracting water bodies in densely built urban areas; the land water body extraction index value LSWI is used to evaluate the distribution and changes of surface water bodies, which has important significance for water resources management and environmental monitoring; the automatic water body extraction index AWEI can consistently improve the extraction accuracy of water bodies in the presence of various environmental noise, while providing stable threshold values, which have two forms, AWEIsh value and AWEInsh value.

[0073] Meanwhile, the following calculation formula is used:

[0074]

[0075] Obtaining the preset each vegetation index value corresponding to each pixel position in the optical synthetic satellite image corresponding to the target research area, including the normalized vegetation index value NDVI, the vegetation enhancement index value EVI, to constitute each vegetation index layer corresponding to the optical synthetic satellite image, including the normalized vegetation index layer, the vegetation enhancement index layer; wherein, RED represents the pixel value of the red light band corresponding to the pixel position in the optical synthetic satellite image.

[0076] Then enter step B2.

[0077] In practical application, the normalized vegetation index value NDVI can show the information of plant growth, ecological system vitality and productivity; the vegetation enhancement index value EVI can improve the sensitivity of vegetation, reduce the influence of soil background and atmosphere, and has higher sensitivity and superiority in monitoring vegetation change, and is widely used in grassland degradation monitoring, grassland resource quantitative analysis and other researches.

[0078] Step B2. According to the optical synthetic satellite image, and the preset each first water body index layer and each vegetation index layer corresponding to the optical synthetic satellite image, combined with the dynamic world land cover data set (Dynamic World Land Cover), determine the classification label of each pixel position in the satellite image under the preset water body classification and each non-water body classification, and randomly select a preset number of pixel positions under each classification label to constitute each sample pixel position; further, for each sample pixel position, first apply the maximum submerged range layer in the JRC global land surface water dynamic data set (JRC GSW), respectively for the local image of the optical synthetic satellite image corresponding to each sample pixel position corresponding to the water body classification label, to update the mask processing, and then with the local image of the optical synthetic satellite image corresponding to the sample pixel position, contact the classification label of the sample pixel position, to constitute a sample, and then by each sample to constitute a sample data set, and then enter step B3. In practical application, the design of each non-water body classification includes vegetation classification label, bare land classification label, cloud classification label and snow classification label.

[0079] Step B3. Based on the sample data set, taking the local image in the sample as input and the classification label in the sample as output, the target classification model such as random forest classification model is trained, and the optical image water body recognition model is obtained.

[0080] Step C. According to the classification label of each pixel position in the satellite image under the preset each category, for the radar synthetic satellite image corresponding to the target research area, and the preset second water body index value corresponding to each pixel position in the radar synthetic satellite image, the best water body classification threshold is obtained, and then enter step D.

[0081] The above step C is specifically designed to perform the following steps C1 to C3 in practical application to obtain the optimal water body classification threshold.

[0082] Step C1. The following calculation formula is used:

[0083] SDWI = ln(10 x VV x VH) - 8

[0084] A preset second water body index value corresponding to the dual-polarized radar water body index value SDWI of each pixel position in the radar synthetic satellite image corresponding to the target research area is obtained, and then step C2 is entered, wherein VV and VH represent the vertical polarization backscatter coefficient value and the horizontal polarization backscatter coefficient value corresponding to the pixel position in the radar synthetic satellite image, respectively.

[0085] Step C2. According to the classification labels of each pixel position in the satellite image with respect to the preset categories, the dual-polarized radar water body index value SDWI of each pixel position in the radar synthetic satellite image is counted to determine the water body initial classification threshold T1 of the dual-polarized radar water body index value SDWI of the pixel position in the radar synthetic satellite image for distinguishing the water body classification from each non-water body classification.

[0086] At the same time, the vertical polarization backscatter coefficient value VV involved in the calculation of the dual-polarized radar water body index value SDWI of each pixel position in the radar synthetic satellite image is counted to determine the water body initial classification threshold T2 of the vertical polarization backscatter coefficient value VV of the pixel position in the radar synthetic satellite image for distinguishing the water body classification from each non-water body classification.

[0087] Then step C3 is entered.

[0088] The water body initial classification threshold T1 and the water body initial classification threshold T2 are obtained to ensure high precision of water body sample classification and adapt to the differences in feature characteristics of regional objects across climate zones.

[0089] Step C3. Based on the water body initial classification threshold T1 and the water body initial classification threshold T2, respectively, the water body edge detection in the radar synthetic satellite image under the Canny edge detection operator is determined in combination with the Otsu algorithm (OTSU) to determine the optimal SDWI water body classification threshold of the dual-polarized radar water body index value SDWI of the pixel position in the radar synthetic satellite image for distinguishing the water body classification from each non-water body classification, and the optimal VV water body classification threshold of the vertical polarization backscatter coefficient value VV of the pixel position in the radar synthetic satellite image for distinguishing the water body classification from each non-water body classification.

[0090] In practical applications, to obtain the optimal SDWI water classification threshold and the optimal VV water classification threshold, the following steps C3a-1 to C3a-3 are performed based on the water initial classification threshold T1.

[0091] Step C3a-1. According to the water initial classification threshold T1, the Canny edge detection operator is applied to determine the edge information of each water body in the radar synthetic satellite image, and each water body with an edge perimeter greater than a preset perimeter threshold is obtained as each water body to be analyzed, and then step C3a-2 is entered.

[0092] In step C3a-1, the Canny operator obtains the transition boundary between water and non-water through gradient detection, retains key edge features, and eliminates small-area noise or invalid boundaries by designing a greater-than-pre-set perimeter threshold to generate optimized water edges.

[0093] Step C3a-2. For each water body to be analyzed, the edge buffer region is formed by extending the edge of the water body to be analyzed in the water direction and the water body outside direction by a preset distance, and then the edge buffer region corresponding to each water body to be analyzed is obtained, and then step C3a-3 is entered.

[0094] Step C3a-3. The dual-polarized radar water index value SDWI of each pixel position in the edge buffer region corresponding to each water body to be analyzed is counted, and the Otsu algorithm is applied to maximize the inter-class variance to determine the optimal SDWI water classification threshold of the dual-polarized radar water index value SDWI of the pixel position in the radar synthetic satellite image for distinguishing water classification from each non-water classification.

[0095] At the same time, the following steps C3b-1 to C3b-3 are performed based on the water initial classification threshold T2 to obtain the optimal VV water classification threshold.

[0096] Step C3b-1. According to the water initial classification threshold T2, the Canny edge detection operator is applied to determine the edge information of each water body in the radar synthetic satellite image, and each water body with an edge perimeter greater than a preset perimeter threshold is obtained as each water body to be analyzed, and then step C3b-2 is entered.

[0097] Step C3b-2. For each water body to be analyzed, the edge buffer region is formed by extending the edge of the water body to be analyzed in the water direction and the water body outside direction by a preset distance, and then the edge buffer region corresponding to each water body to be analyzed is obtained, and then step C3b-3 is entered.

[0098] Step C3b-3. Statistics are performed on the vertical polarization backscattering coefficient values VV of each pixel position in the corresponding edge buffer zone region of each water body to be analyzed, and the Otsu algorithm is applied to determine the optimal VV water body classification threshold value of the pixel position in the radar synthetic satellite image for distinguishing the water body classification from each non-water body classification by maximizing the inter-class variance.

[0099] Step D. The optical satellite historical image water body recognition model is applied to process each pixel position in each local image of each optical satellite historical image to obtain the classification label of each pixel position in each optical satellite historical image, and then the water body distribution in each optical satellite historical image is obtained.

[0100] Meanwhile, the optimal water body classification threshold value is compared and judged for each pixel position in each radar satellite historical image. If the dual-polarization radar water body index value SDWI of the pixel position is greater than the optimal SDWI water body classification threshold value, and the vertical polarization backscattering coefficient value VV of the pixel position is less than the optimal VV water body classification threshold value, it is determined that the pixel position is a water body region, otherwise it is determined that the pixel position is a non-water body region. Whether each pixel position in each radar satellite historical image is a water body region is determined, and then the water body distribution in each radar satellite historical image is obtained; then step E is entered.

[0101] Step E. For each historical acquisition time, each union region between the water body distribution in the optical satellite historical image and the water body distribution in the radar satellite historical image at the historical acquisition time is obtained, and the water body distribution of the target research region at the historical acquisition time is formed by each union region. The water body distribution of the target research region corresponding to each historical acquisition time is obtained, and then step F is entered.

[0102] Step F. For each union region between the water body distribution in the optical satellite historical image and the water body distribution in the radar satellite historical image at each historical acquisition time, the following operations are performed to update the union region and then update the water body distribution of the target research region corresponding to each historical acquisition time.

[0103] Firstly, the terrain slope is calculated by using the digital elevation model (DEM) data, and the region with a slope greater than a preset slope threshold value is marked as a mountain shadow area. The union region is masked and processed to remove the mountain shadow area in the union region, reduce the noise of the misextracted water body, and update the union region. Then, according to the global human settlement area layer (GHSL) in the JRC global land surface water data set (JRC GSW), the union region is subjected to secondary mask processing to remove the building shadow area and the misextracted building area in the union region, and the union region is updated.

[0104] The application is applied to the actual situation, and is specifically based on a Google Earth Engine (GEE) cloud platform, utilizes annual long-time series image sets of Sentinel-1 GRD (radar satellite images) and Sentinel-2 MSL (optical satellite images), and is applied to surface water extraction research in the African Blue Nile Basin. The advantages of optical and radar remote sensing images are fully combined, and through active and passive collaborative data processing technology, accurate extraction and dynamic monitoring of surface water in a cross-climate zone are realized.

[0105] In the selection of experimental data, Sentinel-1 (radar satellite images) and Sentinel-2 (optical satellite images) data of the African Blue Nile in 2021 are obtained, as shown in Figure 2 In the processing of Sentinel-2 image data, based on the cloud probability layer (S2_CLOUD_PROBABILITY) provided in the image set, the cloud-covered area of the annual image of the Blue Nile Basin is de-clouded. Through the ee.ImageCollection.median() function in GEE, the monthly images are median-synthesized to generate a monthly median image set, that is, an optical synthetic satellite image and a radar synthetic satellite image. At the same time, the areas with missing data in the monthly median image and the cloud-covered areas are marked as Nodata areas, so as to supplement the water body information of these areas by Sentinel-1 images in the subsequent steps.

[0106] In the processing of Sentinel-1 (radar satellite image) image, first, based on the GEE platform, the formula ee.Image(10).pow(img.divide(10)) is used to convert the radar backscattering coefficient into dimensionless data form, and then RefineLee filtering processing is performed to remove the speckle noise in the radar image. Subsequently, based on the GEE platform, the formula ee.Image(10).multiply(img.log10()) is used to convert the dimensionless data into backscattering coefficient again, to ensure the physical consistency of the image data. Subsequently, the ee.ImageCollection.median() function in GEE is used to generate monthly median synthetic images, to provide high-quality radar image data for subsequent analysis.

[0107] To enhance the ability to distinguish water bodies from other features, five first water indices (NDWI, MNDWI, LSWI, AWEIsh, AWEInsh) and two vegetation indices (NDVI, EVI) were selected. In terms of sample generation, the Dynamic World Land Cover dataset and the JRC Global Surface Water dataset were used as the data source for automatic sample generation. First, the Dynamic World Land Cover dataset was synthesized monthly using the ee.Reducer.mode() function in GEE. Then, the ee.Image.remap() function was used to reclassify the land cover classification samples in the Dynamic World Land Cover dataset into water, vegetation, land, cloud, and snow. For the reclassified dataset, the ee.Image.stratifiedSample() function was used to randomly select 12000 samples from each category as the initial training samples.

[0108] To further improve the quality of water samples, the maximum inundation range layer of the JRC Global Surface Water dataset (JRC GSW) was used to mask the water samples in the initial training samples, and finally 10560 high-quality water sample points were retained. After combining the water samples with other land cover classification samples, the backscattering coefficient of VV polarization and the eigenvalue of SDWI were calculated. According to the distribution of eigenvalues of different land cover categories, the initial classification threshold of water body in VV polarization was set to -17.2, and the initial classification threshold of SDWI water index was set to 0.19, which was used as the initial classification standard for the subsequent classification model.

[0109] Random forest classification method is widely used in remote sensing classification tasks due to its strong anti-overfitting ability, excellent ability to handle large-scale data, and high robustness. In this study, the automatically generated sample data was input into the random forest model to classify the Sentinel-2 optical image and obtain the monthly water distribution results in the study area.

[0110] To supplement the water information in the Nodata area of Sentinel-2 image, Sentinel-1 data was used for water extraction, such as Figure 3The Edge-OTSU algorithm was used to process the radar image, which automatically determines the best threshold based on the initial threshold to improve the accuracy of water body classification. The VV polarization data of Sentinel-1 radar image can effectively reflect the scattering characteristics of water body due to its strong penetration; the SDWI index combines the VV and VH polarization characteristics, further enhancing the discrimination ability of water body and other ground objects.

[0111] The initial classification value of VV and SDWI obtained by statistics was used as the threshold, and the ee.Algorithms.CannyEdgeDetector() function in GEE was used to extract the water body edge in Sentinel-1 image, which preliminarily represented the water body boundary in the image. Then, the edge strength threshold was set to 0.05, and the water body edge information greater than the threshold was retained, and the noise and weak boundary information was removed, to generate the optimized Sentinel-1 water body boundary. The optimized boundary was processed by buffer zone, and the buffer zone width was set to 500 meters, and the histogram distribution in the buffer zone was counted. Based on the OTSU algorithm, the water body classification was carried out on the VV polarization and SDWI water body index respectively, and the intersection of the extraction results of the two was taken, to finally obtain the water body extraction result of Sentinel-1 image.

[0112] In the water body extraction result fusion stage, the water body extraction result of Sentinel-1 image was merged with the water body extraction result of Sentinel-2 image to generate seamless surface water distribution data of the Blue Nile Basin. To further improve the accuracy of the results, the slope threshold was set to 5 degrees, and the DEM data of the Shuttle Radar Topography Mission (SRTM) was used to calculate the slope by using the ee.Terrain.slope() function in GEE, and the area with slope greater than 5 degrees was masked as mountain shadow area, so as to remove the possible mis-extracted water body results. In addition, the water body extraction result was masked again combined with the global human settlement layer (GHSL), to remove the mis-extracted area caused by building shadow.

[0113] Finally, the high-precision monthly water body extraction result of Sentinel-1 and Sentinel-2 data cooperation was obtained, as shown in Figure 4 The results achieved seamless and dynamic water body distribution monitoring in complex environments across climate zones, providing solid data support and technical support for water resource management and dynamic change analysis in the African Blue Nile Basin.

[0114] The application designs a cross-climate zone surface water resource active-passive cooperative remote sensing extraction method. First, based on satellite historical images, optical synthetic satellite images and radar synthetic satellite images of a target research area are obtained by preprocessing, then a water body recognition model of the optical images is trained, and the best water body classification threshold is obtained by connecting the radar synthetic satellite images, the water body distribution is recognized from the model recognition and the classification threshold, and finally the recognition results of the two ways are fused to obtain the water body distribution of the target research area corresponding to each historical acquisition time; the optical and radar data are deeply fused, large-scale and long-time sequence rapid automatic water body mapping is realized, the accuracy and reliability of the extraction result are significantly improved, the limitations of a single data source can be made up, important support can be provided for dynamic monitoring and scientific management of the water resources of a basin, and the method is especially suitable for water resource monitoring of a complex cross-climate zone basin, and important data support and technical support are provided for water resource management, ecological protection and hydrological dynamic change analysis of the African Nile basin in application.

[0115] The embodiments of the application are described in detail above with reference to the drawings, but the application is not limited to the above-described embodiments, and various changes can be made within the knowledge of those skilled in the art without departing from the purpose of the application.

Claims

1. A method for cross-climate zone surface water resource active and passive collaborative remote sensing extraction, for extracting water body area distribution in a target study area, characterized in that, Comprising the following steps: Step A. Obtain optical satellite historical images and radar satellite historical images corresponding to each preset historical acquisition time within a preset historical period for a target study area, respectively, and obtain optical composite satellite images and radar composite satellite images corresponding to the target study area through preprocessing, and then enter step B; Step B. Based on the optical composite satellite images corresponding to the target study area, and each first water body index layer and each preset vegetation index layer corresponding to the optical composite satellite images, obtain the classification labels of each pixel position in the satellite images with respect to the preset water body classification and each non-water body classification, and train the target classification model to obtain an optical image water body recognition model; then enter step C; Step C. According to the classification labels of each pixel position in the satellite images with respect to each preset category, obtain the best water body classification threshold value for the radar composite satellite images corresponding to the target study area and the preset second water body index value corresponding to each pixel position in the radar composite satellite images, and then enter step D; The above step C comprises the following steps C1 to C3 to obtain the best water body classification threshold value; Step C1. Obtain the preset second water body index value corresponding to each pixel position in the radar synthetic satellite image corresponding to the target study area, and then enter Step C2. Step C2. Obtain the water body index value corresponding to each pixel position in the radar synthetic satellite image corresponding to the target study area, and then enter Step C3. Step C2. According to the classification label of each pixel position in the satellite image about the preset each category, the dual-polarization radar water body index value of each pixel position in the radar synthetic satellite image is calculated Statistics are performed to determine the dual-polarization radar water body index value of the pixel position in the radar synthetic satellite image Water body initial classification threshold for distinguishing water body classification from each non-water body classification ; Meanwhile, the value of the vertical polarization backscattering coefficient corresponding to the position of each pixel in the radar synthetic satellite image is determined The value of the vertical polarization backscattering coefficient involved in the calculation Statistics are performed to determine the value of the vertical polarization backscattering coefficient corresponding to the position of each pixel in the radar synthetic satellite image A water body initial classification threshold value for distinguishing the water body classification from each non-water body classification ; Then enter step C3; Step C3. Based on the water body initial classification threshold , the water body initial classification threshold , according to the water body edge detection in the radar synthetic satellite image, combined with the Otsu algorithm, the dual-polarization radar water body index value corresponding to the pixel position in the radar synthetic satellite image is determined The best water body classification threshold for distinguishing water body classification from each non-water body classification , and the vertical polarization backscatter coefficient value corresponding to the pixel position in the radar synthetic satellite image The best water body classification threshold for distinguishing water body classification from each non-water body classification ; Step D. Apply the optical image water body recognition model to process each local image of each pixel position in each optical satellite historical image to obtain the classification labels of each pixel position in each optical satellite historical image, and further obtain the water body distribution in each optical satellite historical image; At the same time, compare and determine whether each pixel position in each radar satellite historical image is a water body area by contacting the best water body classification threshold value, and further obtain the water body distribution in each radar satellite historical image; Then enter step E; In the above step D, for each pixel position in each radar satellite historical image, it is determined whether the dual-polarization radar water index value of the pixel position is Greater than optimal Water body classification threshold, and the vertical polarization backscatter coefficient value of the pixel position Less than optimal The water body classification threshold is: if it is, the pixel position is determined to be a water body area; otherwise, the pixel position is determined to be a non-water body area. This determines whether each pixel position in each radar satellite historical image is a water body area. Step E. Obtain each union region between the water body distribution in the optical satellite historical image and the water body distribution in the radar satellite historical image at each historical acquisition time, and form the water body distribution of the target study area at the historical acquisition time from each union region, and further obtain the water body distribution of the target study area at each historical acquisition time.

2. The method according to claim 1, wherein the method is characterized in that: In step A, steps A1-1 to A1-3 and steps A2-1 to A2-2 are simultaneously performed for the obtained optical satellite historical images and radar satellite historical images corresponding to each preset historical acquisition time within a preset historical period for a target study area: Step A1-1. Screen each optical satellite historical image to obtain each optical primary satellite image with a cloud area ratio less than a preset ratio threshold in the entire image, and then enter step A1-2; Step A1-2. Remove the cloud local area covering at least a preset lower limit number of pixels in each optical primary satellite image, update the optical primary satellite image, and then enter step A1-3; Step A2-1. Screen each radar satellite historical image to obtain each radar primary satellite image with a cloud area ratio less than a preset ratio threshold in the entire image, and then enter step A2-2; Step A2-2. Remove the cloud local area covering at least a preset lower limit number of pixels in each radar primary satellite image, update the radar primary satellite image, and then enter step A2-3; Step A1-3. For each pixel position in the satellite image, obtain the median value of the pixel values of the pixel position corresponding to each optical primary satellite image, to form the pixel median value corresponding to the pixel position; then combine the pixel median values corresponding to each pixel position to form the optical synthetic satellite image. Step A2-1. For each radar satellite historical image, first perform boundary noise removal processing update, then use the Refine Lee filtering algorithm to remove the scattering noise in the radar satellite historical image to form the radar primary satellite image, and then proceed to step A2-2. Step A2-2. For each pixel position in the satellite image, obtain the median value of the pixel values of the pixel position corresponding to each radar primary satellite image, to form the pixel median value corresponding to the pixel position; then combine the pixel median values corresponding to each pixel position to form the radar synthetic satellite image.

3. The method according to claim 1, wherein the method is characterized in that: The step B includes the following steps B1 to B3; Step B1. Obtain the preset first water body index values corresponding to each pixel position in the optical synthetic satellite image corresponding to the target study area, including the normalized water body index value , the normalized difference water body index value , the land water extraction index value , and the value and the value in the automatic water body extraction index, to form each first water body index layer corresponding to the optical synthetic satellite image, including the normalized water body index layer, the normalized difference water body index layer, the land water extraction index layer, and the automatic water body extraction index layer. Meanwhile, preset each vegetation index value corresponding to each pixel position in the optical synthetic satellite image of the target research area is obtained, including normalized vegetation index value , vegetation enhancement index value , and each vegetation index layer corresponding to the optical synthetic satellite image is constituted, including normalized vegetation index layer and vegetation enhancement index layer. Then proceed to step B2; Step B2. According to the optical synthetic satellite image, and the preset each first water body index layer and each vegetation index layer corresponding to the optical synthetic satellite image, combined with the dynamic world land classification data set, determine the classification labels of each pixel position in the satellite image under the preset water body classification and each non-water body classification, and randomly select a preset number of pixel positions under each classification label to form each sample pixel position. Form a sample by combining the local image of the sample pixel position in the optical synthetic satellite image with the classification label of the sample pixel position, and then form a sample data set by combining each sample, and then proceed to step B3; Step B3. Based on the sample data set, take the local image in the sample as the input and the classification label in the sample as the output, train the target classification model to obtain the optical image water body recognition model.

4. The method according to claim 3, wherein the method is characterized in that: In step B2, for each sample pixel position randomly selected, first apply the maximum inundation range layer in the JRC global land surface water dynamic data set to perform mask processing update on the local image of the optical synthetic satellite image corresponding to each sample pixel position corresponding to the water body classification label, and then form a sample by combining the local image of the sample pixel position in the optical synthetic satellite image with the classification label of the sample pixel position.

5. The method according to claim 3, wherein the method is characterized in that: The non-water body classification includes vegetation classification label, bare land classification label, cloud classification label, and snow classification label.

6. The method of claim 1, wherein the method is characterized by: the step C3a-1 to step C3a-3 are performed as follows: , the step C3a-1 to step C3a-3 are performed as follows: Step C3a-1. According to the initial classification threshold of water body , the Canny edge detection operator is applied to determine the edge information of each water body in the radar synthetic satellite image, and each water body with an edge perimeter greater than a preset perimeter threshold is obtained as each water body to be analyzed, and then step C3a-2 is entered. Step C3a-2. For each water body to be analyzed, extend the edge of the water body by a preset distance in the water direction and the water body outside direction to form an edge buffer region, and then obtain the edge buffer region corresponding to each water body to be analyzed, and then proceed to step C3a-3; Step C3a-3. Bipolarization radar water index value of each pixel position in the edge buffer zone region corresponding to each water body to be analyzed Statistics are performed, and the Otsu algorithm is applied to determine the pixel position corresponding bipolarization radar water index value in the radar synthetic satellite image by maximizing the inter-class variance The best water body classification threshold value for distinguishing water body classification from each non-water body classification Water body classification threshold value At the same time, based on the initial classification threshold of the water body , the following steps C3b-1 to C3b-3 are performed; Step C3b-1. According to the initial classification threshold of water body , the Canny edge detection operator is applied to determine the edge information of each water body in the radar synthetic satellite image, and each water body with an edge perimeter greater than a preset perimeter threshold is obtained as each water body to be analyzed, and then step C3b-2 is entered. Step C3b-2. For each water body to be analyzed, extend the edge of the water body by a preset distance in the water direction and the water body outside direction to form an edge buffer region, and then obtain the edge buffer region corresponding to each water body to be analyzed, and then proceed to step C3b-3; Step C3b-3. The vertical polarization backscattering coefficient values of each pixel position in the edge buffer zone region corresponding to each water body to be analyzed are determined Statistics are performed, and the Otsu algorithm is applied to determine the vertical polarization backscattering coefficient values of the pixel positions in the radar synthetic satellite image by maximizing the inter-class variance The optimal threshold value for distinguishing the water body classification from each non-water body classification Water body classification threshold value.

7. The method of claim 1, wherein the method is characterized by: Further comprising the following step F, after step E is executed, proceed to step F; Step F. After step E is executed, proceed to step F. Step F. For each of the union regions between the water body distribution in the historical optical satellite image and the water body distribution in the historical radar satellite image at each historical acquisition time, the following operations are performed to update the union region and further update the water body distribution of the target study area at each historical acquisition time: First, the union region is masked and the mountain shadow region in the union region is removed to update the union region according to the region with a slope greater than a preset slope threshold. Then, the union region is subjected to secondary masking to remove the building shadow region and the misextracted building region in the union region according to the global human settlement region layer in the JRC global land surface water data set, and the union region is updated.

Images