A method for extracting water area from remote sensing images combining terrain features with observation data

Abstract

The method for extracting water areas from remote sensing images, combining topographic features with observational data, includes the following steps: Step 1: Collect and preprocess images; Step 2: Generate DEM data to obtain the first water extraction result; Step 3: Calculate the vegetation coverage (FVC) of the study area to obtain the second water extraction result; Step 4: Use the maximum inter-class variance method, using the threshold obtained when the error is minimized as the final threshold to obtain the third water extraction result; Step 5: Eliminate the impact of ships in the water body to obtain the fourth water extraction result; Step 6: Correct the fourth water extraction result to obtain the final water extraction result; Step 7: Divide the river into sections according to the characteristics of the river section, and establish a functional relationship between water level and water surface area for each section. By combining multiple sources of data, including topographic features, spectral characteristics, water level observations, and cross-section observations, and performing multiple and detailed processing, a long-term water area monitoring method is constructed, thereby improving the accuracy of water area monitoring.

Description

Technical Field

[0001] The present invention relates to the technical field of remote sensing water surface area extraction, and in particular to a remote sensing image water body area extraction method combining terrain features and observation data. Background Art

[0002] The middle and upper reaches of the Yangtze River are characterized by complex terrain, high water vapor content in the surrounding air, and frequent cloud and fog obstructions. This results in a very limited supply of available optical imagery, making it difficult to construct long-term water area monitoring data. Radar imagery, which penetrates clouds and fog and is unaffected by weather, effectively compensates for this shortcoming. However, in radar imagery, the backscatter coefficient of mountain shadows is similar to that of water bodies, which can lead to mountain shadows being misidentified as water. When using threshold methods to extract water bodies, the threshold value depends solely on statistical characteristics and is subject to uncertainty. Furthermore, the frequent passage of ships on the surface of navigable rivers and reservoirs can cause water surfaces obscured by ships to be classified as non-water bodies. The affected area is often larger than the actual projected area of ​​the ships, severely interfering with the extraction results. Due to the limited revisit period of satellite remote sensing, it is difficult to obtain high-precision daily-scale water surface area changes. Existing studies have mostly estimated this by constructing a functional relationship between water level and water surface area. However, the cross-sectional morphology of water surfaces in actual rivers, reservoirs, etc. varies greatly and is affected by various factors such as season, water storage and release, and flow patterns. It is difficult for traditional single functional relationships to accurately describe them. By determining the functional relationship and fitting the monthly scale parameter set through terrain segmentation analysis, multi-source data and long-series remote sensing water surface extraction, the analysis accuracy of daily scale water surface changes and the accuracy of water area monitoring will be effectively improved. Summary of the Invention

[0003] The technical problem to be solved by the present invention is to provide a method for extracting water body areas from remote sensing images that combines terrain features with observation data. In areas with complex terrain and variable river shapes, during the water body extraction process based mainly on radar SAR images, easily accessible multi-source data such as terrain, water level, and cross-sectional observations are used to eliminate the influence of mountain shadows and ships, thereby improving the spatial accuracy of water area extraction over large areas such as rivers and reservoirs, and achieving daily-scale monitoring.

[0004] In order to solve the above technical problems, the technical solution adopted by the present invention is:

[0005] The method for extracting water area from remote sensing images by combining terrain features with observation data includes the following steps:

[0006] Step 1: Collect multiple synthetic aperture radar (SAR) images of the study area and perform filtering, decibelization, and cropping preprocessing to serve as the basic base map for water surface extraction;

[0007] Step 2: Generate the latest digital elevation model (DEM) data using synthetic aperture radar interferometry (InSAR) technology, and calculate the nearest neighbor river channel relative height model (HAND) through hydrological confluence analysis. Combine HAND with the water body range determined by the maximum water depth observation value of the remote sensing image on the corresponding date in Step 1 to obtain the first water body extraction result.

[0008] Step 3: Use spectral data to calculate the vegetation coverage (FVC) of the study area. Take the FVC mean value of pixels in typical vegetation-covered areas around the river water body as the threshold. Based on the first water body extraction results, remove the high vegetation coverage areas with FVC greater than the threshold to obtain the second water body extraction results.

[0009] Step 4: Use the two-threshold method to extract water bodies. Based on the second extraction result data, use the maximum inter-class variance method, i.e., the OTSU algorithm, to determine the initial threshold for water body / non-water body classification. Use the water level observation data and river section data to determine the cross-section water body boundary position points. Based on the initial threshold, calibrate the water body boundary. The threshold obtained when the error is minimized is used as the final threshold to obtain the third water body extraction result.

[0010] Step 5: Using the Global Surface Water dataset (GSW), which is the permanent water body data, the pixels in Step 4 that are located within the permanent water body and are incorrectly classified as non-water bodies due to the influence of ships and floating objects are defined as water bodies, thus obtaining the fourth water body extraction result.

[0011] Step 6. Calculate the water body index MuWI using multispectral data and set the water body threshold so that the water body boundary includes the main and tributary boundaries of the fourth water body extraction result, excluding the fragments that are far away from the main and tributary but mistakenly classified as water bodies, thereby correcting the final water body extraction result.

[0012] Step 7. Divide the river channel into zones according to the cross-sectional characteristics of the river channel. Establish a functional relationship between water level and water surface area in each zone. Consider the influence of seasonal changes and reservoir storage and release, and fit a monthly parameter set to infer the water area on a daily scale using daily water level observations.

[0013] In Step 1 above, the radar image is taken from the Sentinel-1 satellite GRD data as an example. There are four polarization modes: VV, VH, HH, and HV. The VH polarization mode image data is used, and the GammaMap algorithm and a 7×7 filter window are used to filter the VH data:

[0014]

[0015] Where: represents the backscatter coefficient of VH image data, It is the unit of decibel.

[0016] In Step 2 above, two phases of Sentinel-1 satellite SLC data and the Shuttle Radar Topography Mission (SRTM) 30-meter resolution DEM data downloaded from the U.S. Geological Survey website were used as reference DEMs. The ENVI 5.6 software Sarscape Insardem Workflow module was used to generate updated DEMs for areas where erosion and deposition of reservoirs and rivers have caused significant topographic changes. The DEMs were then used to calculate flow direction, discharge, and river networks, and finally generate the HAND.

[0017] In the above Step 3, the vegetation coverage FVC average of the pixels in the typical vegetation coverage area around the river water body in the study area is used as the threshold, and the high vegetation coverage areas with FVC greater than the threshold are eliminated to complete the second water body extraction; the specific process is as follows:

[0018] The calculation formula of vegetation index NDVI is:

[0019]

[0020] Where NIR stands for near-infrared band, and RED stands for red band, which correspond to band 8 and band 4 in Sentinel-2 images respectively.

[0021]

[0022] Where NDVI soil and NDVI veg Represents the pixel value of zero vegetation coverage and the value of complete vegetation coverage, respectively. The 5% and 95% of NDVI data are used instead of NDVI. soil and NDVI veg , when NDVI is less than NDVI soil When FVC=0, when NDVI is greater than NDVI veg When FVC=1;

[0023] The average FVC of the typical vegetation-covered area around the river water body is 0.8. The area with FVC>0.8 is defined as the high vegetation coverage area and removed from the first water body extraction result image to obtain the second water body extraction result.

[0024] The above Step 4 is based on the VV image of the second water body extraction result. The maximum inter-class variance algorithm OTSU is first used to obtain the initial threshold t for water body / non-water body segmentation; the boundary position points of the cross-section water body determined by the water level observation data and the river section data are used to calibrate the boundary of the water body based on the initial threshold. The threshold obtained when the error is minimum is used as the final threshold T0 to obtain the third water body extraction result; assuming that the grayscale range of the image is [0, T], the corresponding grayscale level The number of pixels is , the total number of pixels N is:

[0025]

[0026] The pixels in the image are divided into two parts a and b according to the gray level t, then the probability of occurrence of pixels a and b is for:

[0027]

[0028]

[0029] The average gray value of A and B and for:

[0030]

[0031]

[0032] The gray value u of the entire image is:

[0033]

[0034] Between-class variance for:

[0035]

[0036] In the range [0, T], the threshold t is increased in steps of 1, and finally When it reaches the maximum value, t is the optimal threshold.

[0037] In the third water body extraction result obtained by the threshold method in Step 5 above, the part blocked by the ship is classified as non-water body. The non-water body is judged. If the non-water body location is within the range of permanent water body, it is judged as water body. The connectedPixelCount() function is used on Google Earth Engine (GEE) to perform 8-neighborhood connectivity analysis to identify broken patches in the image. The focal_mode() function is then used to perform modal filtering, that is, the value of each pixel is replaced by the value that appears most frequently in its neighborhood.

[0038] The above Step 5 calculates the water index MuWI of multiple phases of L2A images of Sentinel-2 satellite after atmospheric correction, and uses the ee.ImageCollection.max() function in Google Earth Engine GEE to synthesize the maximum value of MuWI;

[0039] The calculation formula of water mass index MuWI is:

[0040]

[0041] Where: ρ is the reflectivity, Green is the green light band, Blue is the blue light band, Nir is the near-infrared band, Mir1 is the mid-infrared band 1, and Mir2 is the mid-infrared band 2, corresponding to Band 3, Band 2, Band 8, Band 11, and Band 12 of the Sentinel-2 satellite, respectively.

[0042] The above Step 6 is divided according to the cross-sectional shape of the river. The cross-sectional shape of the river can be divided into four types: irregular, triangular, trapezoidal, and parabolic. Different functional relationships are used to construct the relationship between the ground observed water level and the water surface area of ​​the river in the study area.

[0043] The water surface area in Step 6 above For water bodies such as rivers and reservoirs, the length L along the flow direction is relatively stable and is a constant. The water surface area depends on the water surface width y. The relationship between the water surface width and the water level of the river is related to the shape of the river section and can be classified and fitted. Therefore, the water surface area uses a quadratic function in the parabolic river section. For fitting, the triangular and trapezoidal river sections can be classified into one category, using a linear function Fitting: The stepped river channel is fitted using a piecewise function. Taking the three-step form as an example, it can be expressed as:

[0044]

[0045] Where S is the water surface area, h is the water level, and the formula is divided into sections according to the cross-sectional shapes at different heights; when the water level is When the cross section of the river is parabolic, the formula is When the water level is to The cross section of the river is approximately trapezoidal, and the formula is When the water level is higher than When , the river cross section is parabolic, and the formula is .

[0046] This invention provides a method for extracting water areas from remote sensing images that combines topographic features with observational data. GammaMap filtering is applied to Sentinel-1 VV polarimetric radar data, followed by decibelization and cropping for preprocessing. This serves as the base map for water surface extraction. Accurate water distribution is obtained through five-pass analysis of multi-source data. The multispectral water index (MuWI) is calculated using the time series of Sentinel-2 imagery. Monthly maximum filtering is performed to eliminate areas with MuWI values ​​less than 0.2 from the water extraction results, ultimately yielding the fifth refined water extraction result. According to the cross-sectional characteristics of the river channel, the entire river channel under study is divided into zones, and functional relationships between water level and water surface area are established respectively. Taking into account the influence of seasonal changes and reservoir water storage and release, a monthly parameter set is formed, and the daily water level observation value can be used to infer the water area on a daily scale. The present invention solves the influence of mountain shadows and ships on water extraction from radar images in areas with complex terrain. By combining terrain characteristics, spectral characteristics, water level observation data, cross-sectional observations and other multi-source data for multiple fine processing, a long-term water area monitoring method is constructed, thereby improving the accuracy of water area monitoring. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] The present invention will be further described below with reference to the accompanying drawings and examples:

[0048] Figure 1 This is a flowchart of the calculation process of the present invention Figure 1 ;

[0049] Figure 2 This is a flowchart of the calculation process of the present invention Figure 2 ;

[0050] Figure 3 Schematic diagram of the parabolic river channel cross section in the river channel end face of the present invention;

[0051] Figure 4 Schematic diagram of a triangular river channel cross section in the river channel end face of the present invention;

[0052] Figure 5 Schematic diagram of the trapezoidal river channel cross section in the river channel end face of the present invention;

[0053] Figure 6 Schematic diagram of the stepped river channel cross section at the end face of the river channel in the present invention;

[0054] Figure 7 This is the VH polarization image of Sentinel-1 after pre-processing in the embodiment;

[0055] Figure 8 The left image is a VV polarization image in the embodiment, and the right image is a VH polarization image;

[0056] Figure 9 This is the first water extraction result in the embodiment;

[0057] Figure 10 This is the second water extraction result in the embodiment;

[0058] Figure 11 This is the water body extraction result using OTSU threshold in the example;

[0059] Figure 12 This is the water body extraction result after adjusting the OTSU threshold value based on the river section observation points in the embodiment;

[0060] Figure 13 This is the third water extraction result in the embodiment;

[0061] Figure 14 This is the fourth water extraction result in the embodiment;

[0062] Figure 15 This is the fifth water extraction result in the embodiment;

[0063] Figure 16 It is the change of water surface area on a daily basis in the embodiment. DETAILED DESCRIPTION

[0064] The technical solution of the present invention is described in detail below with reference to the accompanying drawings and embodiments.

[0065] The method for extracting water area from remote sensing images by combining terrain features with observation data includes the following steps:

[0066] Step 1: Collect multiple synthetic aperture radar (SAR) images of the study area and perform filtering, decibelization, and cropping preprocessing to serve as the basic base map for water surface extraction;

[0067] Step 2: Generate the latest digital elevation model (DEM) data using synthetic aperture radar interferometry (InSAR) technology, and calculate the nearest neighbor river channel relative height model (HAND) through hydrological confluence analysis. Combine HAND with the water body range determined by the maximum water depth observation value of the remote sensing image on the corresponding date in Step 1 to obtain the first water body extraction result.

[0068] Step 3: Use spectral data to calculate the vegetation coverage (FVC) of the study area. Take the FVC mean value of pixels in typical vegetation-covered areas around the river water body as the threshold. Based on the first water body extraction results, remove the high vegetation coverage areas with FVC greater than the threshold to obtain the second water body extraction results.

[0069] Step 4: Use the two-threshold method to extract water bodies. Based on the second extraction result data, use the maximum inter-class variance method, i.e., the OTSU algorithm, to determine the initial threshold for water body / non-water body classification. Use the water level observation data and river section data to determine the cross-section water body boundary position points. Based on the initial threshold, calibrate the water body boundary. The threshold obtained when the error is minimized is used as the final threshold to obtain the third water body extraction result.

[0070] Step 5: Using the Global Surface Water dataset (i.e., the permanent water body data in GSW), define the pixels in Step 4 that are located within the permanent water body and are incorrectly classified as non-water bodies due to the influence of ships and floating objects as water bodies, thus obtaining the fourth water body extraction result.

[0071] Step 6. Calculate the water body index MuWI using multispectral data and set the water body threshold so that the water body boundary includes the main and tributary boundaries of the fourth water body extraction result, excluding the fragments that are far away from the main and tributary but mistakenly classified as water bodies, thereby correcting the final water body extraction result.

[0072] Step 7. Divide the river channel into zones according to the cross-sectional characteristics of the river channel. Establish a functional relationship between water level and water surface area in each zone. Consider the influence of seasonal changes and reservoir storage and release, and fit a monthly parameter set to infer the water area on a daily scale using daily water level observations.

[0073] In Step 1 above, radar images are taken from the Sentinel-1 satellite GRD data as an example. There are four polarization modes: VV, VH, HH, and HV. The VH polarization mode image data is used, and the GammaMap algorithm and a 7×7 filter window are used to filter the VH data:

[0074]

[0075] Where: represents the backscatter coefficient of VH image data, It is the unit of decibel.

[0076] In Step 2 above, two phases of SLC data from the Sentinel-1 satellite and 30-meter-resolution DEM data from the Shuttle Radar Topography Mission (SRTM) downloaded from the U.S. Geological Survey website were used as reference DEMs. The ENVI 5.6 software (Survey Insardem Workflow) module was used to generate updated DEMs for areas where erosion and deposition of reservoirs and rivers have caused significant topographic changes. The DEMs were then used to calculate flow direction, discharge, and river networks, and finally generate the HAND.

[0077] The first water extraction is accomplished by setting a Hand and Area (HAND) threshold based on actual ground water level observations from the corresponding date in the remote sensing image. Specifically, the Hand and Area threshold is the maximum observed water depth in the remote sensing image for that month. This is used to identify water and non-water areas. A mask analysis is then performed on the VV data basemap to remove non-water areas and most mountain shadows. For example, if the SAR image was taken on January 8, 2020, the average water depth for January was 140.94m, with a maximum depth of 143.7m. In this case, a Hand and Area threshold of 143.7m is set, and the corresponding water area is used as a mask to extract the VV data basemap.

[0078] In the above Step 3, the high vegetation coverage area is eliminated by setting a threshold through the vegetation coverage FVC of the study area, completing the second water body extraction; the specific process is as follows:

[0079] The value of vegetation coverage FVC in the study area is between 0 and 1. The closer it is to 1, the higher the vegetation coverage. The average FVC of the typical area with vegetation coverage around the river water body is 0.8. The area with FVC>0.8 is defined as a high vegetation coverage area and is removed from the first water body extraction result image to obtain the second water body extraction result.

[0080] For example, if the SAR image was taken on January 8, 2020, the NDVI of the Sentinel-2 image for the study area in January is calculated, and then the ee.ImageCollection.max() function in GEE is used to synthesize the monthly maximum value of the NDVI. The vegetation index FVC is calculated using the synthesized monthly maximum value NDVI, and finally the areas with FVC greater than 0.8 in the SAR image on January 8 are eliminated.

[0081] The calculation formula of vegetation index NDVI is:

[0082]

[0083] Where NIR stands for near-infrared band, and RED stands for red band, which correspond to band 8 and band 4 in Sentinel-2 images respectively.

[0084]

[0085] Where NDVI soil and NDVI veg Represents the pixel value of zero vegetation coverage and the value of complete vegetation coverage, respectively. The 5% and 95% of NDVI data are used instead of NDVI. soil and NDVI veg , when NDVI is less than NDVI soil When FVC=0, when NDVI is greater than NDVI veg When , take FVC=1.

[0086] In Step 4 above, based on the VV image of the second water body extraction result, the maximum inter-class variance algorithm (OTSU) is first used to obtain the initial threshold t for water body / non-water body segmentation. The water body boundary position points of the cross-section determined by the water level observation data and the river section data are used to calibrate the water body boundary based on the initial threshold. The threshold obtained when the error is minimized is used as the final threshold T0 to obtain the third water body extraction result.

[0087] The maximum inter-class variance algorithm OTSU algorithm is also known as the Otsu algorithm. Its principle is to divide the entire image into target pixels and background pixels; assuming that the grayscale range of the image is [0, T], the corresponding grayscale level The number of pixels is , the total number of pixels N is:

[0088]

[0089] The pixels in the image are divided into two parts a and b according to the gray level t, then the probability of occurrence of pixels a and b is for:

[0090]

[0091]

[0092] The average gray value of A and B and for:

[0093]

[0094]

[0095] The gray value u of the entire image is:

[0096]

[0097] Between-class variance for:

[0098]

[0099] In the range [0, T], the threshold t is increased in steps of 1, and finally When it reaches the maximum value, t is the optimal threshold.

[0100] The above Step 5 uses the Permanent water layer in the JRC Global Surface Water YearlyHistory 1.4 dataset from 2013 to 2021, and identifies raster permanent water bodies with a history of more than 8 years as permanent water bodies to eliminate the impact of ships in water body extraction;

[0101] Specifically, in the third water body extraction result obtained by the threshold method, the part blocked by the ship is classified as non-water body. If the non-water body is located near a permanent water body, it is determined to be a water body. The connectedPixelCount() function is used on Google Earth Engine (GEE) to perform 8-neighborhood connectivity analysis to identify broken patches in the image. The focal_mode() function is then used for modal filtering, that is, the value of each pixel is replaced by the value that appears most frequently in its neighborhood.

[0102] The above Step 5 calculates the water index MuWI of multiple Sentinel-2 satellite images, and uses the ee.ImageCollection.max() function in Google Earth Engine GEE to synthesize the maximum value of MuWI;

[0103] The MuWI is a water index developed based on Sentinel-2 satellite data, which can generate 10m high-precision water body maps. According to the time of the SAR image, the MuWI maximum value filtering result of the corresponding month is selected to eliminate the grids with MuWI less than 0.2 in the fourth extraction result. For example: the imaging time of the radar SAR image is January 8, 2020. After the fourth water body extraction, the MuWI of the Sentinel-2 image in January for many years is calculated and maximum value filtering is performed. The area with a value less than 0.2 in the result is non-water body. The corresponding pixels are removed from the fourth water body extraction result to complete the final extraction of the water body.

[0104] The calculation formula of water mass index MuWI is:

[0105]

[0106] Where: ρ is the reflectivity, Green is the green light band, Blue is the blue light band, Nir is the near-infrared band, Mir1 is the mid-infrared band 1, and Mir2 is the mid-infrared band 2, corresponding to Band 3, Band 2, Band 8, Band 11, and Band 12 of the Sentinel-2 satellite, respectively.

[0107] The above Step 6 divides the river into sections according to the cross-sectional morphology. The river sections can be divided into four types: irregular, triangular, trapezoidal, and parabolic. Different functional relationships are used to construct the relationship between the ground-observed water level and the water surface area of ​​the river in the study area. Considering the influence of seasonal changes and reservoir storage and release, the parameters in different months may be different. To reduce the error, a monthly-scale parameter set is formed by monthly fitting. The daily water level observation value can be used to infer the daily water area, which is more accurate.

[0108] The water surface area in Step 6 above For water bodies such as rivers and reservoirs, the length L along the flow direction is relatively stable and is a constant. The water surface area mainly depends on the water surface width y. The relationship between the water surface width and the water level of the river is related to the shape of the river section and can be classified and fitted. Therefore, the water surface area uses a quadratic function in the parabolic river section. For fitting, the triangular and trapezoidal river sections can be classified into one category, using a linear function Fitting: The stepped river channel has a complex shape, so a piecewise function is used for fitting. Taking the three-step form as an example, it can be expressed as:

[0109]

[0110] Where S is the water surface area, h is the water level, and the formula is divided into sections according to the cross-sectional shapes at different heights; when the water level is When the cross section of the river is parabolic, the formula is ;according to Figure 6 The cross section shown, when the water level is to The cross section of the river is approximately trapezoidal, and the formula is When the water level is higher than When , the river cross section is parabolic, and the formula is ; If it is divided into more levels, the same can be applied.

[0111] The representative river sections in the middle and upper reaches of the Yangtze River were selected as cases. According to the cross-sectional characteristics of the upstream and downstream river channels, they can be divided into four types: parabolic, trapezoidal, triangular, and stepped. The relationship curve between water surface area and water level was fitted. The water surface area was calculated from December 27, 2019 to January 2, 2021. During this period, a total of 27 SAR images were collected by the Sentinel-1 satellite. The parabolic river sections are as follows: Figure 3 As shown, the fitting function is: , ; Triangular river section such as Figure 4 As shown, the fitting function is: , 64; ladder-shaped river sections such as Figure 5 As shown, the fitting function is: , ; Step-shaped river channel section such as Figure 6 As shown, the water level is higher than 160m. At this time, the river channel presents a parabolic shape, and the fitting function is: ,

[0112] Example:

[0113] This experiment used the middle and upper reaches of the Yangtze River as a case study. Using Sentinel-1 and Sentinel-2 remote sensing imagery from 2019 to 2022, data from the Global Water Survey dataset (GSW), water levels observed at surface hydrological stations and water level stations, and river cross-section data, water body extraction was performed on the GEE platform. The specific steps of the experiment were based on a Sentinel-1 radar image taken on September 4, 2020.

[0114] Step 1: Collect images and pre-process them. In the backscatter coefficient image of VV polarization mode, the ship appears as a lot of star-shaped bright spots, which even affect the water body around the ship. Therefore, the image of VH polarization mode is selected as the base map for water body extraction to avoid the influence of bright spots. Figure 7 and 8 As shown in .

[0115] Step 2: Generate a HAND model. Combine the remote sensing image time and ground observation water level data in step 1 to determine the HAND threshold and generate the first water body extraction result. Taking the Sentinel-1 image on September 4, 2020 as an example, the maximum water depth in September was 132.51m, and the maximum elevation difference between the upstream and downstream of the river was about 3m. The HAND threshold was set to 135.51m; the first water body extraction result is as follows Figure 9 shown.

[0116] Step 3: Calculate the vegetation coverage FVC of the study area to obtain the second water extraction result; calculate the NDVI of all Sentinel-2 remote sensing images in September over the past four years and perform maximum synthesis to calculate the FVC for September, thus avoiding the impact of fog on optical images and causing data loss. The average FVC of typical areas with vegetation coverage around river water bodies is 0.8. Areas with FVC>0.8 are defined as high vegetation coverage areas. The high vegetation coverage areas in the first water extraction results are eliminated, and the second water extraction results are as follows: Figure 10 shown.

[0117] Step 4: Use the maximum inter-class variance method and the water boundary points in the river cross-section diagram to determine the optimal segmentation threshold: First, use the maximum inter-class variance algorithm to determine the initial segmentation threshold, and apply this threshold to water body segmentation, such as Figure 11As shown; secondly, check the boundary point position of the water body boundary and the water body cross-section monitoring data in the figure, and fine-tune the threshold until the distance between the water body boundary and the cross-section observation data boundary is minimized. This step uses ground observation data to improve the accuracy of remote sensing image analysis using the maximum inter-class variance method. The results are as follows Figure 12 and 13 shown.

[0118] Step 5: Remove the influence of ships in the water body and obtain the fourth water body extraction result; the backscatter coefficient of the ship is much higher than that of the water body and close to the non-water body pixels, so the threshold method cannot be used to eliminate its interference and it will be mistakenly classified as a non-water body; add GSW permanent water body data to eliminate the influence of ships; if the permanent water body includes non-water body pixels, it is defined as a water body, thus obtaining the fourth water body extraction result, the result is as follows Figure 14 shown.

[0119] Step 6: Correct the fourth water extraction result to obtain the final water extraction result: The fourth water extraction result contains many misclassified patches outside the permanent water body, which need to be corrected; calculate the MuWI water index of all September images of Sentinel-2 in the past four years and perform maximum synthesis to avoid the problem of limited available data and low accuracy caused by the influence of clouds on optical images. Finally, eliminate the areas where MuWI is less than 0.2 to complete the correction of the fourth water extraction result. The results are as follows Figure 15 shown.

[0120] Step 7. Divide the river into different zones according to the cross-sectional characteristics, and establish a functional relationship between water level and water surface area for each zone: the middle and upper reaches of the Yangtze River from the Three Gorges Dam to Badong are divided into four types: parabolic, trapezoidal, triangular, and stepped; based on the water area extracted from each pair of Sentinel-1 images in 2020, a functional relationship is constructed. The parabolic river section is as follows: Figure 3 As shown, the fitting function is: , ; Triangular river section such as Figure 4 As shown, the fitting function is: , 64; ladder-shaped river sections such as Figure 5 As shown, the fitting function is: , ; Step-shaped river channel section such as Figure 6 As shown, the water level is higher than 160m. At this time, the river channel presents a parabolic shape, and the fitting function is: , Finally, the daily scale water area is calculated based on the water level, and the daily scale water surface area of ​​the entire segment is obtained by accumulating each segment. The results are as follows Figure 16 shown.

Claims

1. A method for extracting water area from remote sensing images by combining terrain features with observation data, characterized in that: The following steps are involved: Step 1: Collect multiple synthetic aperture radar (SAR) images of the study area and perform filtering, decibelization, and cropping preprocessing to serve as the basic base map for water surface extraction. Step 2: Generate the latest digital elevation model (DEM) data using synthetic aperture radar interferometry (InSAR) technology, and calculate the nearest neighbor river channel relative height model (HAND) through hydrological confluence analysis. Combine HAND with the water body range determined by the maximum water depth observation value of the remote sensing image on the corresponding date in Step 1 to obtain the first water body extraction result. Step 3: Use spectral data to calculate the vegetation coverage (FVC) of the study area. Take the FVC mean value of pixels in typical vegetation-covered areas around the river water body as the threshold. Based on the first water body extraction results, remove the high vegetation coverage areas with FVC greater than the threshold to obtain the second water body extraction results. Step 4: Use the two-threshold method to extract water bodies. Based on the second extraction result data, use the maximum inter-class variance method, i.e., the OTSU algorithm, to determine the initial threshold for water body / non-water body classification. Use the water level observation data and river section data to determine the cross-section water body boundary position points. Based on the initial threshold, calibrate the water body boundary. The threshold obtained when the error is minimized is used as the final threshold to obtain the third water body extraction result. Step 5: Using the Global Surface Water dataset (GSW), which is the permanent water body data, the pixels in Step 4 that are located within the permanent water body and are incorrectly classified as non-water bodies due to the influence of ships and floating objects are defined as water bodies, thus obtaining the fourth water body extraction result. Step 6. Calculate the water body index MuWI using multispectral data and set the water body threshold so that the water body boundary includes the main and tributary boundaries of the fourth water body extraction result, excluding the fragments that are far away from the main and tributary but mistakenly classified as water bodies, thereby correcting the final water body extraction result. Step 7. Divide the river channel into zones according to the cross-sectional characteristics of the river channel. Establish a functional relationship between water level and water surface area in each zone. Consider the influence of seasonal changes and reservoir storage and release, and fit a monthly parameter set to infer the water area on a daily scale using daily water level observations.

2. The method for extracting water area from remote sensing images by combining terrain features with observation data according to claim 1, characterized in that: The radar image in Step 1 is taken as an example of Sentinel-1 satellite GRD data, which has four polarization modes: VV, VH, HH, and HV. The image data of VH polarization mode is used, and the GammaMap algorithm and 7×7 filter window are used to filter the VH data: Where: represents the backscatter coefficient of VH image data, is the backscatter coefficient of the VH image data expressed in decibels.

3. The method for extracting water area from remote sensing images by combining terrain features with observation data according to claim 1, characterized in that: In Step 2, two phases of Sentinel-1 satellite SLC data and the Shuttle Radar Topography Mission (SRTM) 30-meter resolution DEM data downloaded from the U.S. Geological Survey website were used as reference DEMs. The ENVI 5.6 software Sarscape Insardem work flow module was used to generate updated DEMs for areas where reservoirs and river erosion and deposition have caused significant changes in topography. The DEMs were then used to calculate flow direction, discharge, and river networks, and finally generate the HAND.

4. The method for extracting water area from remote sensing images by combining terrain features with observation data according to claim 1, characterized in that: In Step 3, the vegetation coverage FVC average of the pixels in the typical vegetation coverage area around the river water body in the study area is used as a threshold, and the high vegetation coverage areas with FVC greater than the threshold are eliminated to complete the second water body extraction; the specific process is as follows: The calculation formula of vegetation index NDVI is: Where NIR stands for near-infrared band, and RED stands for red band, which correspond to band 8 and band 4 in Sentinel-2 images respectively. Where NDVI soil and NDVI veg Represents the pixel value of zero vegetation coverage and the value of complete vegetation coverage, respectively. The 5% and 95% of NDVI data are used instead of NDVI. soil and NDVI veg , when NDVI is less than NDVI soil When FVC=0, when NDVI is greater than NDVI veg When FVC=1; The average FVC of the typical area with vegetation coverage around the river water body is 0.

8. The area with FVC > 0.8 is defined as the high vegetation coverage area and removed from the first water body extraction result image to obtain the second water body extraction result.

5. The method for extracting water area from remote sensing images by combining terrain features with observation data according to claim 1, characterized in that: The Step 4 is based on the VV image of the second water body extraction result. The maximum inter-class variance algorithm OTSU is first used to obtain the initial threshold t for water body / non-water body segmentation; the boundary position points of the cross-section water body determined by the water level observation data and the river section data are used to calibrate the boundary of the water body on the basis of the initial threshold. The threshold obtained when the error is minimum is used as the final threshold T0 to obtain the third water body extraction result; assuming that the grayscale range of the image is [0, T], the corresponding grayscale level The number of pixels is , the total number of pixels N is: The pixels in the image are divided into two parts a and b according to the gray level t, then the probability of occurrence of pixels a and b is for: The average gray value of A and B and for: The gray value u of the entire image is: Between-class variance for: In the range [0, T], the threshold t is increased in steps of 1, and finally When it reaches the maximum value, t is the optimal threshold.

6. The method for extracting water area from remote sensing images by combining terrain features with observation data according to claim 1, characterized in that: In the third water body extraction result obtained by the threshold method in Step 5, the part blocked by the ship is classified as a non-water body. The non-water body is judged. If the non-water body is located within the range of a permanent water body, it is judged as a water body. The connectedPixelCount() function is used on Google Earth Engine (GEE) to perform 8-neighborhood connectivity analysis to identify broken patches in the image. The focal_mode() function is then used to perform modal filtering, that is, the value of each pixel is replaced by the value that appears most frequently in its neighborhood.

7. The method for extracting water area from remote sensing images by combining terrain features with observation data according to claim 1, characterized in that: The Step 5 calculates the water index MuWI of multiple phases of L2A images of Sentinel-2 satellite after atmospheric correction, and uses the ee.ImageCollection.max() function in Google Earth Engine GEE to perform MuWI maximum value synthesis; The calculation formula of water mass index MuWI is: Where: ρ is the reflectivity, Green is the green light band, Blue is the blue light band, Nir is the near-infrared band, Mir1 is the mid-infrared band 1, and Mir2 is the mid-infrared band 2, corresponding to Band 3, Band 2, Band 8, Band 11, and Band 12 of the Sentinel-2 satellite, respectively.

8. The method for extracting water area from remote sensing images by combining terrain features with observation data according to claim 1, characterized in that: The Step 6 is divided into four types according to the cross-sectional shape of the river. The cross-sectional shape of the river can be divided into four types: irregular, triangular, trapezoidal, and parabolic. Different functional relationships are used to construct the relationship between the ground observation water level and the water surface area of ​​the river in the study area.

9. The method for extracting water area from remote sensing images by combining terrain features with observation data according to claim 1, characterized in that: The water surface area in Step 6 , for the water body of river channel and reservoir, the length L along the flow direction is relatively stable and constant, and the water surface area depends on the water surface width y; the relationship between the water surface width and the water level of the river channel is related to the cross-sectional shape of the river channel and can be classified and fitted. Therefore, the water surface area uses a quadratic function in the parabolic river channel cross-section For fitting, the triangular and trapezoidal river sections can be classified into one category, using a linear function Fitting: The stepped river channel is fitted using a piecewise function, which is divided into three steps and can be expressed as: Where S is the water surface area, h is the water level, and the formula is divided into sections according to the cross-sectional shapes at different heights; when the water level is When the cross section of the river is parabolic, the formula is When the water level is to The cross section of the river is approximately trapezoidal, and the formula is When the water level is higher than When , the river cross section is parabolic, and the formula is .

Images