A water surface boundary extraction method based on historical frequency and timing correction

Abstract

The application discloses a water surface boundary extraction method based on historical frequency and timing correction, comprising: preprocessing optical remote sensing images and radar remote sensing images; generating a preliminary water body binary mask based on MNDWI and NDVI; constructing a historical water body frequency graph; obtaining T * based on MNDWI and Otsu for the current month; * carrying out water body segmentation based on T * ; taking the union of the region meeting the radar VH polarization band threshold condition and the water body mask to obtain a candidate base map; constructing a screening rule; combining the adjacent month water surface extraction result and the historical water body frequency to correct the candidate base map in timing consistency. The application generates a candidate base map based on the optical remote sensing images of the current month, combines radar data and terrain data, and carries out timing correction through the historical water body frequency and the adjacent month water surface extraction result obtained by the same single month water surface extraction process, so that the accuracy and robustness of the lake and reservoir water surface boundary extraction are improved.

Description

Technical Field

[0001] This invention relates to the fields of hydrological technology and remote sensing image processing, and in particular to a method for extracting water surface boundaries based on historical frequency and time series correction. Background Technology

[0002] Extraction and long-term continuous monitoring of lake and reservoir water surfaces are fundamental tasks in the field of surveying and mapping geographic information technology. The accuracy and efficiency of extraction directly affect the quality of basic water resource surveys, flood emergency response, and wetland ecological protection. Currently, remote sensing technology, due to its advantages in acquiring a series of spatial, temporal, and thematic Earth observation information, provides an effective way to observe surface water dynamics. Commonly used methods include single-band thresholding, multi-band exponential methods, traditional machine learning methods, and deep learning methods.

[0003] Currently, most methods for extracting water surface data in lakes and reservoirs are limited to short-term or coarse-scale monitoring, lacking long-term and fine-scale data support. This makes it impossible to fully grasp the changing trends of water surface at different time and spatial scales, especially the subtle changes under the influence of rapid climate change and human activities. While popular water surface extraction methods such as machine learning have high accuracy, they rely on a large number of high-precision training samples, making the collection cost of such samples too high for long-term monitoring. Thresholding methods are relatively low-cost and suitable for long-term water surface monitoring, but they suffer from the drawback of difficulty in applying universal thresholds to extract water surface boundaries due to spatial heterogeneity. Currently, thresholding methods often rely on optical remote sensing imagery, lacking effective integration with radar and topographic data. Furthermore, their processing is often limited to single-temporal images, neglecting historical water stability and temporal consistency information between adjacent months in the study area. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of the prior art and provide a method for extracting water surface boundaries by historical frequency and time-series correction. Based on the optical remote sensing image of the current month, a candidate base map is generated by combining radar data and topographic data. Time-series correction is performed by using historical water body frequency and water surface extraction results of adjacent months obtained according to the same single-month water surface extraction process, thereby improving the accuracy and robustness of lake and reservoir water surface boundary extraction.

[0005] To achieve the above objectives, the technical solution adopted by this invention is: a water surface boundary extraction method based on historical frequency and time-series correction, comprising: acquiring optical remote sensing images and radar remote sensing images of the target area; preprocessing the optical remote sensing images and radar remote sensing images to obtain monthly composite optical remote sensing reflectance images and monthly composite radar remote sensing images; calculating the improved normalized water index (MNDWI) and normalized vegetation index (NDVI) based on the monthly composite optical remote sensing reflectance images; generating a preliminary water body binary mask corresponding to the time phase based on preset threshold rules of the MNDWI and NDVI; constructing a historical water body frequency map based on the preliminary water body binary mask; and obtaining the optimal threshold T for the current month based on the MNDWI and the global Otsu adaptive thresholding method. * Based on the optimal threshold T of the current month * Water body segmentation is performed to obtain the initial water body mask for the current month. The regions in the current month's composite radar remote sensing image that meet the radar VH polarization band threshold conditions are then joined with the current month's water body mask to obtain the candidate base map for the current month. Threshold filtering rules are constructed to remove noise from the candidate base map for the current month. The candidate base map for the current month is then temporally consistent with the water surface extraction results from the previous month, the current month, the next month, and historical water body frequencies obtained using the same single-month water surface extraction process. The temporally consistent candidate base map for the current month is then morphologically repaired, and isolated small patches are removed based on connected component size or patch morphology characteristics to obtain the lake and reservoir water surface boundary extraction results.

[0006] Preferably, preprocessing the optical remote sensing image and radar remote sensing image to obtain the monthly composite optical remote sensing reflectance image and the monthly composite radar remote sensing image includes: using the median composite method to perform monthly-scale composite of the optical remote sensing image band reflectance and the radar remote sensing image VH polarization band; and performing RefinedLee filtering on the monthly composite radar remote sensing image VH polarization band to facilitate the suppression of speckle noise and the preservation of boundary information.

[0007] Preferably, the improved Normalized Difference Water Index (MNDWI) is expressed by the following formula:

[0008] MNDWI=(Green-SWIR) / (Green+SWIR)

[0009] In the formula: Green represents the reflectivity of the green light band; SWIR represents the reflectivity of the short-wave infrared band;

[0010] The Normalized Difference Vegetation Index (NDVI) is expressed by the following formula:

[0011] NDVI = (NIR - Red) / (NIR + Red)

[0012] In the formula: NIR represents the reflectivity of the near-infrared band; Red represents the reflectivity of the red band.

[0013] Preferably, the preset threshold rules for MNDWI and NDVI include: MNDWI > 0 and NDVI < 0;

[0014] Pixel values ​​of the preliminary water body binary mask in the optical remote sensing image It can be expressed as follows:

[0015] .

[0016] Preferably, the historical water body frequency H(x) is expressed by the following formula:

[0017]

[0018] In the formula: N represents the number of historical images used. This represents the initial binary mask pixel value of pixel x at time t, where water pixels are set to 1 and non-water pixels are set to 0.

[0019] Preferably, the optimal threshold obtained by the Otsu adaptive thresholding method is expressed by the following formula:

[0020]

[0021] In the formula: Here, T represents the inter-class variance, and T represents the candidate threshold.

[0022] It can be expressed as follows:

[0023]

[0024] In the formula: and These represent the proportions of the two types of pixels under threshold T; and These represent the mean values ​​of the two types of pixels, respectively. This represents the population mean.

[0025] Preferably, the initial water mask pixel value for the current month is expressed by the following formula:

[0026]

[0027] In the formula: x represents a pixel; MNDWI(x) represents the improved normalized water index of pixel x; T * Indicates the optimal threshold;

[0028] The region in the current month's composite radar remote sensing image that meets the radar VH polarization band threshold condition is represented by the following formula:

[0029]

[0030] In the formula: This is the backscattering intensity threshold.

[0031] Preferably, the threshold screening rules include slope control, historical frequency map control, radar control, MNDWI index control, and single-band reflectivity control;

[0032] The slope control is as follows:

[0033]

[0034] In the formula: Take 15 degrees;

[0035] The historical frequency control is as follows:

[0036]

[0037] In the formula: T h The value is 0.02;

[0038] The MNDWI index is controlled as follows:

[0039]

[0040] In the formula: Take -0.1;

[0041] The radar control is as follows:

[0042]

[0043] In the formula: Take -15dB;

[0044] The shortwave infrared band reflectivity control is as follows:

[0045]

[0046] In the formula: Take 0.12;

[0047] The near-infrared band reflectivity control is as follows:

[0048]

[0049] In the formula: Take 0.18.

[0050] Preferably, the timing consistency correction includes:

[0051] If a pixel is represented as a water body in the current month, and as a non-water body in the previous and next months, and the historical water body frequency of the pixel is less than the first threshold T1, the pixel is corrected to be a non-water body.

[0052] If a pixel is non-water body in the current month, water body in the previous month and the following month, and the historical water body frequency of the pixel is greater than the second threshold T2, then the pixel is restored to water body.

[0053] If a pixel is not a water body in the current month, is only a water body in the previous month or only in the next month, and the historical water body frequency of the pixel is greater than the third threshold T3, then the pixel is restored to a water body.

[0054] Among them, the first threshold T1 is less than the second threshold T2, and the second threshold T2 is less than the third threshold T3.

[0055] Preferably, morphological repair of the modified candidate base map for the current month includes: performing neighborhood dilation, neighborhood mode filtering, and neighborhood erosion on the candidate base map to fill small holes inside the lake or reservoir and smooth the boundaries; further, based on the connected component size threshold, isolated small patches are removed, and target areas that conform to the surface water characteristics of the lake or reservoir are retained.

[0056] Compared with the prior art, the present invention has the following advantages:

[0057] 1. Based on the optical remote sensing image of the current month, this invention generates candidate base maps by combining radar data and topographic data, and performs time-series correction by using historical water body frequency and water surface extraction results of adjacent months obtained according to the same single-month water surface extraction process, thereby improving the accuracy and robustness of lake and reservoir water surface boundary extraction.

[0058] 2. This invention uses the filtered VH polarization band as auxiliary constraint information, which reduces the uncertainty and complexity of optical and radar fusion while retaining the radar's anti-cloud advantage.

[0059] 3. The water surface boundary extraction method of the present invention can be implemented on a remote sensing cloud computing platform without sample labeling and model training. The calculation method is convenient and suitable for the automated processing of multi-monthly and long-term series of lake and reservoir boundaries, thus improving the efficiency of water surface boundary extraction.

[0060] The present invention will now be described in further detail with reference to the accompanying drawings and embodiments. Attached Figure Description

[0061] Figure 1 This is a flowchart of the present invention;

[0062] Figure 2 This is a schematic diagram illustrating the extraction results of the extraction method of the present invention for lakes in different months;

[0063] Figure 3 This is a schematic diagram illustrating the extraction results of the extraction method of the present invention for different months in a reservoir;

[0064] Figure 4 This is a schematic diagram of the extraction results of the reservoir under cloudy conditions. Detailed Implementation

[0065] like Figure 1 As shown, this invention discloses a method for extracting water surface boundaries based on historical frequency and time-series correction, comprising:

[0066] Acquire optical and radar remote sensing images of the target area;

[0067] The water surface boundary extraction in this application mainly focuses on the extraction of water surface boundaries of lakes and reservoirs, obtaining long-term optical remote sensing images and radar remote sensing images of the target area for subsequent analysis.

[0068] The optical remote sensing images and radar remote sensing images are preprocessed to obtain a monthly composite optical remote sensing reflectance image and a monthly composite radar remote sensing image.

[0069] Preprocessing the optical and radar remote sensing images yields a monthly composite optical remote sensing reflectance image and a monthly composite radar remote sensing image, including:

[0070] The median composite method was used to synthesize the band reflectance of optical remote sensing images and the VH polarization band of radar remote sensing images on a monthly scale, resulting in monthly composite optical remote sensing reflectance images and monthly composite radar remote sensing images.

[0071] Further refined Lee filtering was applied to the VH polarization band of the lunar synthetic radar remote sensing image to suppress speckle noise and preserve boundary information.

[0072] Based on the aforementioned monthly composite optical remote sensing reflectance image, the improved normalized water index MNDWI and normalized vegetation index NDVI were calculated.

[0073] The improved normalized water index MNDWI is expressed by the following formula:

[0074] MNDWI=(Green-SWIR) / (Green+SWIR) (1)

[0075] In the formula: Green represents the reflectivity of the green light band; SWIR represents the reflectivity of the short-wave infrared band; compared with NDVI, MNDWI has better suppression of non-aquatic interference such as bare land and buildings, and is suitable for the initial extraction of water bodies such as lakes and reservoirs;

[0076] The Normalized Difference Vegetation Index (NDVI) is expressed by the following formula:

[0077] NDVI=(NIR-Red) / (NIR+Red) (2)

[0078] In the formula: NIR represents the reflectivity of the near-infrared band; Red represents the reflectivity of the red band.

[0079] NDVI is mainly used to exclude vegetation disturbance, and vegetated areas often have higher NDVI.

[0080] A preliminary water body binary mask for the corresponding time phase is generated based on the preset threshold rules of MNDWI and NDVI.

[0081] Based on long-term optical remote sensing images, the MNDWI and NDVI indices for each month were calculated. Following the preset rules for MNDWI and NDVI: MNDWI > 0 and NDVI < 0, the pixel values ​​of each pixel in the optical remote sensing image were determined. The following formula represents the generation of a preliminary binary water mask for the corresponding time phase based on each pixel value:

[0082] (3)

[0083] In one possible embodiment, all optical remote sensing images of the target area for a single month are acquired. The median composite method is used to perform monthly-scale composite optical remote sensing images for that single month, resulting in a monthly composite optical remote sensing reflectance image for that month. The improved normalized water index (MNDWI) and normalized vegetation index (NDVI) are calculated for each pixel on the monthly composite optical remote sensing reflectance image. Based on the MNDWI and NDVI of each pixel, combined with a preset threshold: MNDWI > 0 and NDVI < 0, when MNDWI is greater than 0 and NDVI is less than 0, the pixels are considered to be water bodies; otherwise, the pixels are considered non-water bodies. When the pixels are considered to be water bodies... Set to 1 when the pixels in this part are not water bodies. Set it to 0, and then generate the preliminary binary mask for the water body for that month.

[0084] Construct a historical water body frequency map based on the preliminary water body binary mask;

[0085] The historical water body frequency H(x) is expressed by the following formula:

[0086] (4)

[0087] In the formula: N represents the number of historical images used. This represents the initial binary mask pixel value of pixel x at time t, where water pixels are set to 1 and non-water pixels are set to 0.

[0088] In one possible embodiment, all optical remote sensing images of the target area for a single year are acquired. The median composite method is used to perform monthly-scale composite optical remote sensing images for each month of that single year, resulting in monthly composite optical remote sensing reflectance images for each month. The improved normalized water index (MNDWI) and normalized vegetation index (NDVI) for each pixel are calculated sequentially on each monthly composite optical remote sensing reflectance image, ultimately yielding 12 sets of MNDWI and NDVI. Based on the MNDWI and NDVI of each pixel, and combined with preset thresholds: MNDWI > 0 and NDVI < 0, the water body binary mask pixel values ​​for each pixel for the 12 months are obtained. The total number of times water appears in each pixel over 12 months is counted. The total number of times water appears in each pixel is divided by the total number of images to obtain the water frequency H(x) of each pixel for that year. The water frequency H(x) ranges from [0,1]. The water frequency H(x) of each pixel for that year constitutes the water frequency map of the target area for that year.

[0089] In another possible embodiment, all optical remote sensing images of the target area for two consecutive years are acquired. The median composite method is used to perform monthly-scale composite optical remote sensing images for each month of the two years, resulting in monthly composite optical remote sensing reflectance images for each month. The improved normalized water index (MNDWI) and normalized vegetation index (NDVI) are calculated for each pixel on each monthly composite optical remote sensing reflectance image, ultimately yielding 24 sets of MNDWI and NDVI. Based on the MNDWI and NDVI of each pixel, and combined with a preset threshold: MNDWI > 0 and NDVI < 0, the 24 monthly binary mask pixel values ​​for water bodies are obtained for each pixel. The total number of times water appears in each pixel over 24 months is counted. The total number of times water appears in each pixel is divided by the total number of images to obtain the water frequency H(x) of each pixel over two years. The water frequency H(x) ranges from [0,1]. The water frequency H(x) of each pixel over two years constitutes the water frequency map of the target area over two years.

[0090] Furthermore, all optical remote sensing images of the target area for three consecutive years were acquired. The median composite method was used to synthesize the optical remote sensing images for each month of the three years at a monthly scale, resulting in monthly composite optical remote sensing reflectance images for each month. The Improved Normalized Difference Water Index (MNDWI) and Normalized Difference Vegetation Index (NDVI) were calculated for each pixel on each monthly composite optical remote sensing reflectance image, ultimately yielding 36 sets of MNDWI and NDVI. Based on the MNDWI and NDVI of each pixel, and combined with preset thresholds: MNDWI > 0 and NDVI < 0, the binary mask pixel values ​​for water bodies for each pixel in the 36 months were obtained. The total number of times water appears in each pixel over 36 months is counted. The total number of times water appears in each pixel is divided by the total number of images to obtain the three-year water frequency H(x) of each pixel. The value range of water frequency H(x) is [0,1]. The three-year water frequency H(x) of each pixel constitutes the three-year water frequency map of the target area.

[0091] The historical water body frequency H(x) ranges from [0,1]. When H(x) is closer to 1, it indicates that the pixel frequently appears as a water body in each year and has high temporal stability. When H(x) is closer to 0, it indicates that the pixel frequently appears as a non-water body in each year, usually corresponding to non-water body areas or temporary water accumulation areas. The historical water body frequency map formed by the historical water body frequency H(x) can be used as temporal prior information for subsequent water body identification to improve the reliability of the discrimination of stable water body areas.

[0092] The optimal threshold T for the current month is obtained based on the MNDWI and the global Otsu adaptive thresholding method. * ;

[0093] Each month, MNDWI is used as the primary feature for initial water body extraction. The MNDWI histogram of the entire target area is statistically analyzed, and the optimal threshold T is obtained using the Otsu method. * The Otsu method determines the threshold by maximizing the inter-class variance. The optimal threshold obtained by the Otsu adaptive thresholding method is expressed by the following formula:

[0094] (5)

[0095] In the formula: Here, T represents the inter-class variance, and T represents the candidate threshold.

[0096] It can be expressed as follows:

[0097] (6)

[0098] In the formula: and These represent the proportions of the two types of pixels under threshold T; and These represent the mean values ​​of the two types of pixels, respectively. This represents the population mean.

[0099] The two types of pixels are water pre-selected pixels with MNDWI values ​​above the threshold and non-water pre-selected pixels with MNDWI values ​​below the threshold.

[0100] Based on the optimal threshold T of the current month * Perform water body segmentation to obtain the initial water body mask for the current month;

[0101] The initial water mask pixel value for the current month is expressed by the following formula:

[0102] (7)

[0103] In the formula: x represents a pixel; MNDWI(x) represents the improved normalized water index of pixel x; T * This represents the optimal threshold.

[0104] In one possible embodiment, all optical remote sensing images of the target area for the current month are acquired. The median composite method is used to perform monthly-scale composite optical remote sensing images for the current month, resulting in a monthly composite optical remote sensing reflectance image. The MNDWI index value for each pixel is calculated on this monthly composite optical remote sensing reflectance image. The MNDWI histogram for each pixel is then calculated. Finally, the optimal threshold T for the current month is obtained using the Otsu method. * According to the optimal threshold T * Combining the MNDWI(x) exponent value of each pixel, if the MNDWI(x) exponent value of a certain pixel is greater than the optimal threshold T * Then the initial water mask pixel value for the current month. If the value is 1, it represents a water body; otherwise, it represents the initial water body mask pixel value for the current month. Setting the value to 0 indicates that the object is not water. After processing each pixel, the initial water mask for the current month is obtained.

[0105] The candidate base map for the current month is obtained by taking the union of the regions in the current month's composite radar remote sensing image that meet the radar VH polarization band threshold condition with the initial water mask for the current month.

[0106] The region in the current month's composite radar remote sensing image that meets the radar VH polarization band threshold condition is represented by the following formula:

[0107] (8)

[0108] In the formula: This is the backscattering intensity threshold.

[0109] To avoid deleting real water bodies due to cloud cover or local anomalies in the current month, regions meeting the VH polarization band threshold conditions are used as supplementary candidate masks. The union of the supplementary candidate masks and the initial water body mask is taken to obtain the candidate base map for the current month. The backscattering intensity threshold T is used. vh1 Use -20 dB.

[0110] In one possible embodiment, all radar remote sensing images of the target area for the current month are acquired. The VH polarization bands of the radar remote sensing images are then synthesized at a monthly scale using the median composite method to obtain a monthly composite radar remote sensing image. Next, the VH polarization bands of the monthly composite radar image are subjected to Refined Lee filtering to suppress speckle noise and preserve boundary information. The backscattering intensity threshold is set to... , here Using a value of -20dB, the VH polarization band in the lunar synthetic radar image is filtered. If the backscattering intensity of a pixel in the VH polarization band of the lunar synthetic radar image is less than... If the pixel is not a water body, it is marked as a supplementary candidate water body pixel and marked as 1; otherwise, it is marked as a non-water body candidate pixel and marked as 0. Then, a supplementary candidate mask for the current month is generated. The supplementary candidate mask for the current month is merged with the initial water body mask for the current month to obtain the candidate base map for the current month. This avoids the problem of insufficient optical remote sensing image data due to excessive cloud cover or rainy weather in the current month, which could lead to abnormal water surface boundary extraction.

[0111] Construct threshold filtering rules to remove noise from the candidate base map for the current month;

[0112] By constructing threshold screening rules, it is easy to eliminate noise that is misidentified as water. The threshold screening rules include slope control, historical frequency map control, radar control, MNDWI index control, and single-band reflectivity control.

[0113] The slope control is as follows:

[0114] (9)

[0115] In the formula: Take 15 degrees, that is ;

[0116] The historical frequency control is as follows:

[0117] (10)

[0118] In the formula: T h The value is 0.02, that is... ;

[0119] The MNDWI index is controlled as follows:

[0120] (11)

[0121] In the formula: Take -0.1, that is ;

[0122] The radar control is as follows:

[0123] (12)

[0124] In the formula: Take -15dB, that is ;

[0125] The shortwave infrared band reflectivity control is as follows:

[0126] (13)

[0127] In the formula: Take 0.12, that is ;

[0128] The near-infrared band reflectivity control is as follows:

[0129] (14)

[0130] In the formula: Take 0.18, that is .

[0131] In one possible embodiment, some noisy data in the candidate base map for the current month is removed. The filtering conditions include: controlling the slope data to within 15 degrees, that is, retaining data with a slope less than 15 degrees and removing the rest of the slope data; if the historical water body frequency of the candidate pixel meets the preset conditions, it is used as the basis for retention or replenishment; marking pixels with MNDWI less than -0.1 as non-water bodies; marking pixels with radar backscattering intensity greater than -15dB as non-water bodies; marking pixels with shortwave infrared band reflectivity greater than 0.12 as non-water bodies; marking pixels with near-infrared band reflectivity greater than 0.18 as non-water bodies. Through the above filtering, noise that is misidentified as water bodies in the candidate base map for the current month can be eliminated, improving the accuracy of the candidate base map.

[0132] By combining the water surface extraction results of the previous month, the current month, and the next month obtained according to the same single-month water surface extraction process, as well as the historical water body frequency, the candidate base map for the current month is corrected for temporal consistency.

[0133] The timing consistency correction includes:

[0134] If a pixel is represented as a water body in the current month, and as a non-water body in the previous and next months, and the historical water body frequency of the pixel is less than the first threshold T1, the pixel is corrected to be a non-water body.

[0135] Expressed as a formula:

[0136] (15)

[0137] In one possible implementation, assuming the current month is March, the pixel value of a certain pixel in the candidate base image for March is... A value of 1 indicates that the pixel was detected as a body of water in March; the pixel value of this pixel in the candidate base map in February of the same year... The value was 0, representing the pixel value in the candidate base map of the same year in April. Similarly, if the value is 0, it means that the pixel value in the candidate base map of both the month before and the month after March is marked as non-water body; and in the historical water body frequency map, the frequency of the pixel value being marked as 1 is less than the first threshold T1, which can be 0.08. Based on the above judgments, this pixel is abnormal noise, so the pixel value in the candidate base map of March is... The marker is changed to 0, meaning that the pixel is not water in the candidate base map for March.

[0138] If a pixel is non-water body in the current month, water body in the previous month and the following month, and the historical water body frequency of the pixel is greater than the second threshold T2, then the pixel is restored to water body.

[0139] Expressed as a formula:

[0140] (16)

[0141] In one possible implementation, assuming the current month is May, the pixel value of a certain pixel in the candidate base image for May is... A value of 0 indicates that the pixel was detected as a non-water body in May; the pixel value of this pixel in the candidate base map in April of the same year... The value is 1, representing the pixel value in the candidate base map in June of the same year. The value is also 1, meaning that the pixel's water mask pixel value in the candidate base map of both the month before and the month after May is marked as water; and in the historical water frequency map, the frequency of the pixel value being marked as 1 is greater than the second threshold T2, which can be 0.18. Based on the above judgment, this pixel is a missed pixel, so the pixel value in the candidate base map of May is... The marker is changed to 1, meaning that the pixel represents water in the candidate base map for May.

[0142] If a pixel is not a water body in the current month, is only a water body in the previous month or only in the next month, and the historical water body frequency of the pixel is greater than the third threshold T3, then the pixel is restored to a water body.

[0143] Expressed as a formula:

[0144] (17)

[0145] In one possible implementation, assuming the current month is October, the pixel value of a certain pixel in the candidate base image for October... A value of 0 indicates that the pixel was detected as a non-water body in October; the pixel value of this pixel in the candidate base map in September of the same year... The value is 1, representing the pixel value in the candidate base map of the same year, in November. A value of 0 indicates that the pixel was labeled as a water body in the month preceding October and as a non-water body in the month following October; or it could be the pixel value in the candidate base map of the same year in September. A value of 0, representing the pixel value in the candidate base map of the same year, in November. A value of 1 indicates that the pixel was marked as a non-water body in the month preceding October and as a water body in the month following October. In other words, the pixel was only marked as a water body in the month preceding or following the current month. In the historical water body frequency map, the frequency of this pixel being marked as 1 is greater than the third threshold T3, where T3 can be 0.3. Based on the above, this pixel is considered a missed detection pixel, and its pixel value in the candidate base map for October is... The marker is changed to 1, meaning that the pixel in the candidate base image for October is a body of water.

[0146] Among them, the first threshold T1 is less than the second threshold T2, and the second threshold T2 is less than the third threshold T3.

[0147] Morphological repair is performed on the candidate base map of the current month after temporal consistency correction, and isolated small patches are removed based on the scale of connected regions or patch morphology features to obtain the lake and reservoir water surface boundary extraction results.

[0148] Morphological repair of the revised candidate base map for the current month includes: performing neighborhood dilation, neighborhood mode filtering, and neighborhood erosion on the candidate base map to fill small holes inside lakes or reservoirs and smooth the boundaries; further, isolated small patches are removed based on the connected component size threshold, while retaining the target area that conforms to the surface water body characteristics of lakes and reservoirs.

[0149] Small holes in the candidate base map are filled in, isolated small patches are removed, and the final water surface boundary extraction results of lakes and reservoirs are obtained. The extraction results are then visually verified and their accuracy evaluated using high-resolution imagery.

[0150] The above method was used to extract the water surface of a lake in Qindu District, Xianyang City, Shaanxi Province. This water body is stable but has a complex surrounding environment. Sentinel-2 multi-band imagery (10m) and Sentinel-1 radar imagery (10m) from 2024 were used to extract the water surface of the study area. Figure 2 As shown, Figure 2 (a) is a true-color image of the lake's water. Figure 2 (b) The results of water body boundary extraction of the lake in March obtained using the method of this application. Figure 2(c) The results of water body boundary extraction of the lake in July obtained using the method of this application. Figure 2 (d) is the result of water body boundary extraction of the lake in November obtained using the method of this application.

[0151] The above method was used to extract the water surface of a mountainous reservoir in Qian County, Xianyang City, Shaanxi Province. The reservoir experiences drastic water changes and frequent cloud and fog interference. Sentinel-2 multi-band imagery (10m) and Sentinel-1 radar imagery (10m) from 2024 were used to extract the water surface of the reservoir study area. Figure 3 As shown, Figure 3 (a) is a true-color image of the reservoir's water body. Figure 3 (b) The results of water body boundary extraction for the reservoir in March obtained using the method described in this application. Figure 3 (c) The results of water body boundary extraction for the reservoir in July obtained using the method described in this application. Figure 3 (d) shows the water body boundary extraction results of the reservoir in November obtained using the method of this application. Figure 4 The results of water surface boundary extraction from multi-cloud imagery. Figure 4 (a) is a true-color image of the reservoir's water body under cloudy conditions. Figure 4 (b) is the result of extracting the reservoir water surface boundary under cloudy conditions using the method of this application.

[0152] The above description is merely a preferred embodiment of the present invention and does not constitute any limitation on the present invention. Any simple modifications, alterations, or equivalent structural transformations made to the above embodiments based on the technical essence of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A method for extracting water surface boundaries based on historical frequency and time-series correction, characterized in that, include: Acquire optical and radar remote sensing images of the target area; The optical remote sensing images and radar remote sensing images are preprocessed to obtain a monthly composite optical remote sensing reflectance image and a monthly composite radar remote sensing image. Based on the aforementioned monthly composite optical remote sensing reflectance image, the improved normalized water index MNDWI and normalized vegetation index NDVI were calculated. A preliminary water body binary mask for the corresponding time phase is generated based on the preset threshold rules of MNDWI and NDVI. Construct a historical water body frequency map based on the preliminary water body binary mask; The optimal threshold T for the current month is obtained based on the MNDWI and the global Otsu adaptive thresholding method. * ; Based on the optimal threshold T of the current month * Perform water body segmentation to obtain the initial water body mask for the current month; The candidate base map for the current month is obtained by taking the union of the regions in the current month's composite radar remote sensing image that meet the radar VH polarization band threshold condition with the water body mask for the current month. Construct threshold filtering rules to remove noise from the candidate base map for the current month; By combining the water surface extraction results of the previous month, the current month, and the next month obtained according to the same single-month water surface extraction process, as well as the historical water body frequency, the candidate base map for the current month is corrected for temporal consistency. Morphological repair is performed on the candidate base map of the current month after temporal consistency correction, and isolated small patches are removed based on the scale of connected regions or patch morphology features to obtain the lake and reservoir water surface boundary extraction results.

2. The water surface boundary extraction method based on historical frequency and time series correction according to claim 1, characterized in that, Preprocessing the optical and radar remote sensing images yields a monthly composite optical remote sensing reflectance image and a monthly composite radar remote sensing image, including: The median composite method was used to synthesize the band reflectance of optical remote sensing images and the VH polarization band of radar remote sensing images on a monthly scale. Refined Lee filtering is applied to the VH polarization band of the lunar synthetic radar remote sensing image to suppress speckle noise and preserve boundary information.

3. The water surface boundary extraction method based on historical frequency and time series correction according to claim 1, characterized in that, The improved normalized water index MNDWI is expressed by the following formula: MNDWI=(Green-SWIR) / (Green+SWIR); In the formula: Green represents the reflectivity of the green light band; SWIR represents the reflectivity of the short-wave infrared band; The Normalized Difference Vegetation Index (NDVI) is expressed by the following formula: NDVI = (NIR - Red) / (NIR + Red); In the formula: NIR represents the reflectivity of the near-infrared band; Red represents the reflectivity of the red band.

4. The water surface boundary extraction method based on historical frequency and time series correction according to claim 1, characterized in that, The preset threshold rules for MNDWI and NDVI include: MNDWI > 0 and NDVI < 0; Pixel values ​​of the preliminary water body binary mask in the optical remote sensing image It can be expressed as follows: 。 5. A method for extracting water surface boundaries based on historical frequency and time series correction as described in claim 1, characterized in that, The historical water body frequency H(x) is expressed by the following formula: ； In the formula: N represents the number of historical images used. This represents the initial binary mask pixel value of pixel x at time t, where water pixels are set to 1 and non-water pixels are set to 0.

6. A method for extracting water surface boundaries based on historical frequency and time series correction as described in claim 1, characterized in that, The optimal threshold obtained by the Otsu adaptive thresholding method is expressed by the following formula: ； In the formula: Here, T represents the inter-class variance, and T represents the candidate threshold. It can be expressed as follows: ； In the formula: and These represent the proportions of the two types of pixels under threshold T; and These represent the mean values ​​of the two types of pixels, respectively. This represents the population mean.

7. A method for extracting water surface boundaries based on historical frequency and time series correction as described in claim 1, characterized in that, The initial water mask pixel value for the current month is expressed by the following formula: ； In the formula: x represents a pixel; MNDWI(x) represents the improved normalized water index of pixel x; T * Indicates the optimal threshold; The region in the current month's composite radar remote sensing image that meets the radar VH polarization band threshold condition is represented by the following formula: ； In the formula: This is the backscattering intensity threshold.

8. A method for extracting water surface boundaries based on historical frequency and time series correction according to claim 1, characterized in that, The threshold screening rules include slope control, historical frequency map control, radar control, MNDWI index control, and single-band reflectivity control. The slope control is as follows: ； In the formula: Take 15 degrees; The historical frequency control is as follows: ； In the formula: T h The value is 0.02; The MNDWI index is controlled as follows: ； In the formula: Take -0.1; The radar control is as follows: ； In the formula: Take -15dB; The shortwave infrared band reflectivity control is as follows: ； In the formula: Take 0.12; The near-infrared band reflectivity control is as follows: ； In the formula: Take 0.

18.

9. A method for extracting water surface boundaries based on historical frequency and time series correction as described in claim 1, characterized in that, The timing consistency correction includes: If a pixel is represented as a water body in the current month, and as a non-water body in the previous and next months, and the historical water body frequency of the pixel is less than the first threshold T1, the pixel is corrected to be a non-water body. If a pixel is non-water body in the current month, water body in the previous month and the following month, and the historical water body frequency of the pixel is greater than the second threshold T2, then the pixel is restored to water body. If a pixel is not a water body in the current month, is only a water body in the previous month or only in the next month, and the historical water body frequency of the pixel is greater than the third threshold T3, then the pixel is restored to a water body. Among them, the first threshold T1 is less than the second threshold T2, and the second threshold T2 is less than the third threshold T3.

10. A method for extracting water surface boundaries based on historical frequency and time series correction as described in claim 1, characterized in that, Morphological repair of the revised candidate base map for the current month includes: performing neighborhood dilation, neighborhood mode filtering, and neighborhood erosion on the candidate base map to fill small holes inside lakes or reservoirs and smooth the boundaries; further, isolated small patches are removed based on the connected component size threshold, while retaining the target area that conforms to the surface water body characteristics of lakes and reservoirs.

Images