Refined extraction method for small and micro land surface water body based on multi-scale remote sensing

Abstract

The invention discloses a small and micro land surface water body refining extraction method based on multi-scale remote sensing, which comprises the following steps: constructing a high-medium-low resolution image pyramid by using multi-source satellite remote sensing data, realizing self-adaptive wave band combination optimization based on small and micro water body spectral characteristics, and carrying out multi-scale noise suppression and radiation normalization processing; presetting a small and micro water body geometric shape quantification index, constructing a small and micro water body screening mechanism constrained by geometric features, and realizing coarse positioning of a small and micro water body candidate region; a scale weight adaptive distribution algorithm is adopted, contribution weights of data of different scales are dynamically adjusted according to the geometrical characteristics of the candidate areas, a combined discrimination model of spectrum-space-geometrical characteristics of the small and micro water body is established, and accurate extraction of the small and micro water body is achieved; and constructing a boundary optimization method based on a functional energy model, and realizing small and micro water body boundary optimization by combining small-area noise removal, connected region combination and result consistency check of time series data constraint.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of intelligent analysis of remote sensing images, and particularly relates to a multi-scale remote sensing small and micro land surface water body fine extraction method. BACKGROUND

[0002] Small and micro land surface water bodies generally refer to water bodies with an area of 0.01-1 hectares, including small ponds, farmland ditches, small reservoirs and urban landscape water bodies, etc. Such water bodies are an important part of surface water resources, and play a key role in maintaining regional ecological balance, regulating local climate and supporting agricultural irrigation. With the increasing demand for ecological environment governance and fine water resource scheduling management, there is an urgent need for high-precision, large-scale and dynamic monitoring of small and micro water bodies. Traditional ground investigation methods have problems such as high cost and low efficiency, and cannot meet the requirements of large-scale and high-frequency monitoring. Remote sensing technology has become the main technical means for water body monitoring due to its wide coverage, low cost and strong timeliness.

[0003] At present, the water body extraction methods based on remote sensing images mainly include spectral index method, supervised classification method and object-oriented classification method. Although the spectral index method (such as NDWI, MNDWI) is simple to calculate and has high processing efficiency, it is easily affected by mixed pixels and surrounding ground objects, and has limited recognition ability for small and irregularly shaped small and micro water bodies. The supervised classification method relies on a large number of labeled samples and has relatively high classification accuracy, but it still has insufficient ability to distinguish similar spectral features such as shadows and wet soil. The object-oriented classification method can improve the classification effect to some extent through image object segmentation, but its segmentation parameter setting is complex and has low automation, which limits its application in large-scale monitoring.

[0004] In recent years, deep learning technology has made significant progress in semantic segmentation of remote sensing images, and convolutional neural network models such as U-Net and DeepLab have shown good performance in water body extraction tasks. However, the existing methods still have the following problems: first, single scale processing cannot balance the detail preservation and noise suppression of small and micro water bodies; second, there is a lack of special design for the geometric features of small and micro water bodies; third, the boundary extraction accuracy is not high, which cannot meet the requirements of fine application; fourth, the temporal consistency constraint is insufficient, which affects the reliability of long-term monitoring. Therefore, it is urgent to develop a small and micro water body fine extraction method that can effectively process multi-scale information, fully utilize geometric feature constraints and achieve high-precision boundary extraction. SUMMARY

[0005] To solve the technical problems of low recognition accuracy and inaccurate boundary extraction of small and micro water bodies in the prior art, the application provides a small and micro terrestrial surface water body fine extraction method based on multi-scale remote sensing, a small and micro water body fine extraction strategy of three stages of coarse positioning-precise boundary-post optimization is constructed, small and micro water bodies with an area of 0.01-1 hectares can be effectively identified, the boundary positioning accuracy reaches a sub-pixel level, and important support is provided for accurate water resource supervision.

[0006] To achieve the above object, the application provides the following scheme.

[0007] A small and micro terrestrial surface water body fine extraction method based on multi-scale remote sensing, the method comprises the following steps:

[0008] Multi-source satellite remote sensing data is used to construct a high-medium-low resolution image pyramid, adaptive band combination optimization based on small and micro water body spectral characteristics is realized, and multi-scale noise suppression and radiation normalization processing are performed.

[0009] Based on the preprocessed multi-source remote sensing image data, a small and micro water body geometric shape quantization index is preset, a small and micro water body screening mechanism based on geometric feature quantization index is constructed, and coarse positioning of the small and micro water body candidate area is realized.

[0010] Based on the preprocessed multi-source remote sensing image data, a scale weight adaptive distribution algorithm based on water body geometric feature similarity is used, the contribution weight of different scale data is dynamically adjusted according to the geometric feature of the candidate area, a joint discrimination model of small and micro water body spectral-spatial-geometric feature is established, and accurate extraction of small and micro water body is realized.

[0011] A boundary optimization method based on a functional energy model is constructed, and the results of small area noise removal, connected region merging and temporal data constraint consistency test are combined, and small and micro water body boundary optimization is realized.

[0012] Preferably, the method for constructing a high-medium-low resolution image pyramid using multi-source satellite remote sensing data, realizing adaptive band combination optimization based on small and micro water body spectral characteristics, and performing multi-scale noise suppression and radiation normalization processing comprises the following steps:

[0013] The radiation normalization method is used to unify the image element value to the dimensionless reflectivity space, and the processing method is as follows:

[0014]

[0015] Wherein, R i is the normalized reflectivity value of the i-th band; DN i is the original image element value; DN dark and DN white are the element values of dark and bright targets respectively.

[0016] The image is unified to three levels: high resolution layer ≤2m, medium resolution layer 5-30m, and low resolution layer ≥30m using a multi-scale pyramid structure. The spatial resampling uses a bicubic interpolation method, and the formula is:

[0017]

[0018] where f(x, y) represents the pixel value of the image coordinate (x, y) after resampling; DN(x+i, y+j) represents the pixel digital value of (x+i, y+j) in the original image; where i and j are relative coordinate offsets, taking values of -1, 0, 1, 2, corresponding to the 16 neighboring pixels; w(i, j) is the bicubic interpolation weight function, whose value is determined by the pixel offset distance, which can smooth the calculation results and preserve local details; the coordinate (x, y) is the grid coordinate of the target image after resampling, with the upper left corner of the image as the origin, x representing the column direction and y representing the row direction; the offset (i, j) defines the relative position of the neighboring pixels participating in the interpolation in the original image, where i is the horizontal offset and j is the vertical offset;

[0019] For images of different scales, wavelet decomposition and reconstruction methods are used, and the formula is:

[0020] I'(x, y) = ∑ m ∑ n W m,n (x, y)·ψ m,n (x, y);

[0021] where I'(x, y) is the filtered image, ψ m,n (x, y) is the wavelet basis function of scale m and direction n, and W m,n (x, y) is the corresponding wavelet coefficient;

[0022] Combining spectral features and time sequence features for cloud and cloud shadow detection: the Fmask method based on spectral threshold is used to realize automatic identification of cloud and cloud shadow regions, and the cloud mask matrix C(x, y) is obtained; combined with the spectral consistency constraint of time sequence image T t (x, y), the data interpolation of the pixels in the cloud-shaded area is carried out, and the interpolation formula is:

[0023] I''(x, y, t) = α(x, y, t)·I'(x, y, t-1) + (1-α(x, y, t))·I'(x, y, t+1);

[0024] Where I″(x,y,t) represents the image value of the cloud-occluded area after interpolation at time t; I′(x,y,t-1) and I′(x,y,t+1) represent the denoised image values ​​at adjacent times; α(x,y,t) is a weighting factor whose value is dynamically determined based on the spectral similarity of pixels in the time series and the time interval, and can adaptively adjust the contribution of the images before and after to the interpolation result.

[0025] Preferably, the method for coarsely locating candidate regions of small water bodies by constructing a small water body screening mechanism based on the geometric shape quantification index of small water bodies includes:

[0026]

[0027] Among them, F(R) i ) is the joint discriminant function, C i For the comprehensive density index, B i R is the boundary complexity exponent. i Let F(R) be the shape regularity index, and let α, β, γ be the weighting parameters, satisfying α + β + γ = 1; when F(R) i )≥threshold T f At that time, region R i It has been identified as a candidate area for small water bodies.

[0028] Preferably, a scale-weighted adaptive allocation algorithm based on water body geometric similarity is used to dynamically adjust the contribution weights of data at different scales according to the geometric features of candidate regions, and a joint discrimination model of spectral-spatial-geometric features of small water bodies is established to achieve accurate extraction of small water bodies. This method includes:

[0029] Based on the preprocessed multi-source remote sensing image data, denoted as High-Resolution Image I H (x,y), medium resolution image I M (x,y) and low-resolution image I L (x,y) is used to form different scale levels through pyramid decomposition:

[0030] I (s) (x,y)=DownSample(I(x,y),s),s∈{H,M,L};

[0031] Where DownSample(·) represents the downsampling operator, and s is the scale level;

[0032] Standardization is performed using the min-max normalization method:

[0033]

[0034] An improved small water body index W is constructed based on an adaptive band combination optimization method using entropy weighting.s (x,y), this index is achieved by a weighted combination of the Normalized Differential Water Index (NDWI) and the Improved Normalized Differential Water Index (MNNDWI), as shown below:

[0035]

[0036] Among them, W s (x,y) represents the improved micro-water body index of pixel (x,y) at scale s, used to measure the intensity of water body features of that pixel. G (s) NIR (s) With SWIR (s) The green band, near-infrared band, and shortwave infrared band represent the bands at scale s, respectively. Parameters α and β are the band combination weights used to adjust the contribution of each band to water body identification. In the weight determination process, the entropy weight method is introduced, with weights w1 and w2 corresponding to the relative importance of NDWI and MNNDWI, to achieve adaptive optimization of the band combination.

[0037] (α * ,β * ) = (w1, w2);

[0038] Where, α * β represents the band combination coefficient of the optimized NDWI, used to enhance the significance of water body spectral characteristics. Its value is calculated by the entropy weight method, reflecting the relative contribution of NDWI to distinguishing water bodies from non-water bodies; * The band combination coefficients of the optimized MNNDWI are used to supplement the identification ability of small water body edges and complex backgrounds. The values ​​are calculated by the entropy weight method and reflect the relative contribution of MNNDWI to the water body extraction box degree.

[0039] After obtaining water body indices at various scales, a geometric feature quantification index is first constructed by combining the geometric features of candidate water body regions. Simultaneously, a water body probability model P(x,y) is introduced to quantify the probability that each pixel is a water body. Furthermore, a density index and boundary complexity are introduced to define a geometric consistency index Q at scale s. s Based on the geometric consistency index, an adaptive scale weight allocation function is defined, ultimately yielding the multi-scale adaptive water index W. f (x,y):

[0040] W f (x,y)=∑ s∈{H,M,L} ω s ·W s (x,y);

[0041] Among them, W f (x,y) is the multi-scale adaptive water index, ω ss is the adaptive weighting function, where s is the scale level.

[0042] Preferably, boundary optimization methods based on functional energy models include:

[0043]

[0044] Wherein, E(C) is the functional energy of the active contour curve C(s); α is the boundary smoothing control parameter, used to suppress boundary noise and ensure overall contour smoothness; β is the curvature adjustment parameter, used to maintain boundary continuity and prevent excessive bending; λ is the external energy weight, which controls the guiding effect of image information on contour evolution; P(C(s)) is the image gradient potential function, constructed based on the grayscale or spectral feature gradient of water body pixels, which can attract the contour to converge towards the real water body boundary; C′(s) is the first derivative vector, representing the tangent vector of the boundary curve at parameter s, used to describe the local direction change of the curve; and C″(s) is the second derivative vector, representing the rate of change of curvature of the boundary curve at parameter s, used to control the degree of curvature and continuity of the curve.

[0045] Preferred methods for small-area noise removal and connected region merging include:

[0046] Based on boundary refinement, a connected component labeling algorithm is used to analyze small water bodies. Let the set of connected components of small water bodies in the extracted results be {R}. i}, i=1,2,…,N, its area is A(R i ); Perform a threshold operation on each connected region:

[0047] R i ∈Ω if and only if A(R) i )≥A min ;

[0048] Where Ω represents the set of effective connected regions of small water bodies that are retained, that is, all connected regions that satisfy the area threshold, and A min A minimum effective water area threshold is set; areas smaller than this threshold are considered noise and removed. Simultaneously, for spatially adjacent small water bodies with boundaries narrower than the threshold δ, a region merging operation is performed.

[0049]

[0050] Preferably, the methods for verifying the consistency of time-series data constraints include:

[0051] Let R be the extraction result of a certain region at time t. t The result at time t+1 is R. t+1 A temporal similarity constraint function is constructed by integrating spectral similarity, geometric similarity, and spatial neighborhood consistency.

[0052]

[0053] Among them, G(R) t ,R t+1 N(R) represents the consistency index between boundary complexity and shape regularity. t ,R t+1 ) represents the spatial continuity index of the neighboring area, and α, β, and γ are adaptive weighting coefficients determined according to the water body type, area, shape, and background features.

[0054] When S(R) t ,R t+1 When the value is less than the threshold θ, it indicates that there is temporal inconsistency in this region, which requires dynamic optimization and correction.

[0055] Preferably, when S(R) t ,R t+1 When the threshold θ is less than the threshold θ, it indicates that there is temporal inconsistency in this region. Methods for dynamic optimization and correction include:

[0056] Historical results {R1,R2,…,R} are used t-1} and the current result R t The weighted multiple decisions have been updated:

[0057] R t ′=weighted-mode{w1R1,w2R2,...,w t-1 R t-1 ,w t R t};

[0058] Where weighted-mode(·) represents the weighted majority voting function, and the weight w t The allocation is dynamically determined based on the reliability of historical time series, cloud shadow occlusion, and detection confidence.

[0059] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0060] (1) By using multi-scale information fusion and geometric feature constraints, the accuracy of identifying small water bodies has been significantly improved, with a 15-25% improvement in accuracy compared to traditional methods.

[0061] (2) The active contour model is used for boundary refinement, which realizes the precise boundary extraction at the sub-pixel level and controls the boundary positioning error within 1-2 pixels.

[0062] (3) The adaptive weight allocation algorithm can dynamically adjust the processing strategy according to the characteristics of different regions, and has good adaptability to complex terrain and diverse small water body types.

[0063] (4) It can effectively identify ultra-small water bodies with an area of ​​0.01-1 hectare, and the boundary positioning accuracy reaches the sub-pixel level. It is suitable for high-precision intelligent monitoring of small water bodies in fields such as urban waterlogging monitoring, precision agricultural irrigation, dynamic assessment of ecological environment, and refined management of water resources. Attached Figure Description

[0064] To more clearly illustrate the technical solution of the present invention, the drawings used in the embodiments are briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0065] Figure 1 This is a schematic diagram of a method for refining and extracting small-scale land surface water bodies based on multi-scale remote sensing, according to an embodiment of the present invention. Detailed Implementation

[0066] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0067] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0068] Example 1

[0069] like Figure 1 As shown, this invention provides a method for refined extraction of small-scale land surface water bodies based on multi-scale remote sensing, the method comprising:

[0070] Using multi-source satellite remote sensing data, a high-medium-low resolution image pyramid is constructed to achieve adaptive band combination optimization based on the spectral characteristics of small water bodies, and multi-scale noise suppression and radiometric normalization processing are performed.

[0071] Based on preprocessed multi-source remote sensing image data, a geometric shape quantification index for small water bodies, including density index, boundary complexity, and shape regularity, is preset. Based on the geometric shape quantification index of small water bodies, a small water body screening mechanism constrained by geometric features is constructed to achieve coarse localization of candidate areas for small water bodies.

[0072] Based on preprocessed multi-source remote sensing image data, a scale weight adaptive allocation algorithm based on the similarity of water body geometric features is adopted. The contribution weight of data at different scales is dynamically adjusted according to the geometric features of candidate regions. A joint discrimination model of spectral-spatial-geometric features of small water bodies is established to achieve accurate extraction of small water bodies.

[0073] A boundary optimization method based on a functional energy model is constructed, and the results of small-area noise removal, connected region merging, and time-series data constraint consistency verification are combined to achieve boundary optimization of small water bodies.

[0074] In this embodiment, multi-source remote sensing data acquisition:

[0075] 1. The data required for this invention mainly includes the following categories:

[0076] (1) High spatial resolution data

[0077] Acquire sub-meter to meter-level high spatial resolution satellite images from Gaofen-2 (GF-2), Gaofen-6 (GF-6), WorldView series, and GeoEye series covering the study area, with a panchromatic band spatial resolution better than 1 meter.

[0078] (2) Hyperspectral data

[0079] Acquire hyperspectral satellite images of the study area from ZY1-02D, HJ-2A / 2B, Gaofen-5 (GF-5), and EO-1 Hyperion.

[0080] (3) High temporal resolution data

[0081] Acquire multi-temporal resolution satellite images from Sentinel-2A / 2B, Landsat 8 / 9, and domestically produced Gaofen-1 and Gaofen-6 wide-swath satellites within the research period (e.g., 3 months).

[0082] In this embodiment, multi-source remote sensing data preprocessing:

[0083] Before performing refined extraction of small water bodies, this invention first preprocesses multi-source, multi-scale remote sensing data to ensure consistency in spectral, spatial, and radiometric characteristics across images of different resolutions, thus providing reliable input for subsequent geometric feature constraints and boundary refinement. The specific steps are as follows:

[0084] (1) Image radiometric normalization processing

[0085] Due to differences in imaging time, imaging angle, and atmospheric conditions among different satellite sensors, there are significant radiometric deviations between the original images. To achieve comparability of multi-source images, a radiometric normalization method is used to unify the image pixel values ​​into a dimensionless reflectance space. The processing method is as follows:

[0086]

[0087] Among them, R i DN is the normalized reflectivity value for the i-th band; i The original image pixel values; DN dark and DN white These are the pixel values ​​for dark and bright targets, respectively. This normalization process effectively eliminates systematic errors between different data sources.

[0088] (2) Spatial resolution uniformity and pyramid construction

[0089] For high, medium, and low resolution remote sensing data from different sources, this invention proposes a multi-scale image pyramid structure to unify the original remote sensing images into three layers: a high-resolution layer (≤2m), a medium-resolution layer (5–30m), and a low-resolution layer (≥30m), so as to achieve collaborative extraction of information on small water bodies at multiple scales. To achieve spatial resolution consistency, this invention resamples images of different resolutions using a bicubic interpolation method. The resampling formula is as follows:

[0090]

[0091] Where f(x,y) represents the pixel value at image coordinates (x,y) after resampling; DN(x+i,y+j) represents the pixel value at (x+i,y+j) in the original image, where i and j are relative coordinate offsets, ranging from -1, 0, 1, to 2, corresponding to the 16 surrounding neighboring pixels; w(i,j) is the bicubic interpolation weight function, whose value is determined by the pixel offset distance, which can smooth the calculation results and preserve local details; coordinates (x,y) are the grid coordinates after resampling of the target image, usually with the upper left corner of the image as the origin, x representing the column direction and y representing the row direction; offset (i,j) defines the relative position of the neighboring pixels participating in interpolation in the original image, where i is the horizontal offset and j is the vertical offset.

[0092] This method can preserve spatial detail information in multi-scale pyramid levels and avoid excessive amplification of high-frequency noise during high-resolution image resampling, providing a high-quality image foundation for the subsequent fine extraction of small water bodies.

[0093] (3) Multi-scale noise suppression and cloud removal

[0094] To reduce noise and cloud interference during the extraction of small water bodies, this invention employs wavelet transform filtering to suppress high-frequency noise in multi-scale image processing. It also introduces a cloud removal and interpolation mechanism combining spectral and temporal features, thereby achieving effective restoration of missing information and stable detection of water body targets. In the noise suppression stage, this invention uses wavelet decomposition and reconstruction methods for images at different scales, with the following formula:

[0095] I'(x,y)=∑ m ∑ n W m,n (x,y)·ψ m,n (x,y);

[0096] Where I'(x,y) is the filtered image, ψ m,n (x,y) are wavelet basis functions with scale m and direction n, W m,n (x,y) represents the corresponding wavelet coefficients. By setting an adaptive threshold in the high-frequency subband to suppress high-frequency noise and retaining the main spectral features of the water body in the low-frequency subband, the details of small water bodies are preserved while effectively reducing background interference.

[0097] (4) In this invention, to further ensure the continuity and accuracy of monitoring small water bodies in multiple time phases, a cloud and cloud shadow detection and data interpolation method combining spectral and temporal features is proposed. Specifically, it includes the following steps: First, utilizing the spectral features of multi-source remote sensing images, an improved Fmask method based on spectral thresholds is used to automatically identify cloud regions and cloud shadow regions, obtaining the cloud mask matrix C, whose expression is:

[0098]

[0099] Where x and y are the row and column coordinates of the image, C(x, y) = 1 indicates that the pixel at that location is covered by a cloud or cloud shadow, and C(x, y) = 0 indicates that the pixel at that location is a valid observation. Compared with the prior art, this invention not only performs cloud detection based on the spectral threshold of a single image, but also introduces a dynamic threshold adjustment mechanism for small water bodies based on Fmask, to enhance the recognition ability of small shadows and semi-transparent clouds, and achieve accurate masking coverage of small water bodies. Subsequently, for pixel areas obscured by clouds, the time series image T is combined with... i Temporal interpolation of pixel values ​​is performed based on the spectral consistency constraint of (x, y). The interpolation formula is:

[0100] I″(x,y,t)=α(x,y,t)·I′(x,y,t-1)+(1-α(x,y,t))·I′(x,y,t+1);

[0101] Wherein, I″(x,y,t) represents the interpolated image value of the cloud-occluded area at time t; I′(x,y,t-1) and I′(x,y,t+1) represent the denoised image values ​​at adjacent times; α(x,t,y) is a weighting factor, the value of which is dynamically determined based on the spectral similarity of pixels in the time series and the time interval, and can adaptively adjust the contribution of the preceding and following images to the interpolation result. This method can effectively recover the water body information loss caused by cloud interference, ensuring the continuity and integrity of small water body areas in multi-temporal monitoring.

[0102] By using multi-scale wavelet noise suppression and cloud removal and interpolation processing that combines spectral and temporal features, this invention significantly reduces noise interference and cloud occlusion while ensuring the preservation of boundary details of small water bodies and the temporal stability of the extraction results, providing high-quality input data for subsequent geometric feature constraints and boundary refinement.

[0103] The innovation of this invention lies in coupling the weighting factors α, β, and γ with the local spectral statistical characteristics and contrast index of the image to achieve dynamic adaptive optimization, without relying on manual experience to set parameters. Compared with the traditional fixed-parameter MNDWI enhancement method, it can significantly improve the recognition rate and boundary accuracy of small water bodies, and has significant advantages in the detection of small water bodies in complex backgrounds.

[0104] In this embodiment, geometric feature-driven coarse localization of candidate regions for small water bodies is performed:

[0105] To effectively distinguish between small water bodies and pseudo-water bodies (such as shadows, moist soil, and road puddles), a method for locating candidate regions of small water bodies based on geometric feature constraints is proposed. By calculating geometric indices such as the density, boundary complexity, and shape regularity of the target region, rapid screening and optimization of candidate regions for small water bodies are achieved. The specific steps are as follows:

[0106] (1) Extraction of initial water body objects

[0107] To further highlight the spectral differences between water bodies and background features, this invention introduces an adaptive water index enhancement algorithm based on preprocessed multi-source remote sensing data. Using the Normalized Difference Water Index (NDWI) as a foundation and combining it with a weight optimization strategy, it achieves accurate extraction of preliminary water body objects. Specifically, the improved adaptive water index NDWI* expression is defined as follows:

[0108]

[0109] Where G represents the green band reflectance, NIR represents the near-infrared band reflectance, and α, β, and γ are adaptive weighting factors that are dynamically adjusted according to the spectral distribution of the image area.

[0110] In this invention, the initialization of α, β, and γ is globally optimized based on image statistical characteristics. To fully consider the local spectral differences between water bodies and non-water background features in the image, this invention innovatively proposes a method for calculating the global mean and standard deviation based on weighted spatial neighborhoods, rather than traditional simple global statistics. Specifically:

[0111]

[0112] Where SWIR represents the shortwave infrared reflectance, N is the total number of image pixels, and w i The adaptive weight for each pixel is defined as a local spectral similarity weight function of neighboring pixels. Specifically, it is calculated by coupling the spectral gradient with spatial location to enhance the response of the water area and suppress background interference. This neighborhood-weighted statistical method can fully reflect the local spectral characteristics, which is significantly different from the approach in the prior art that only relies on the simple global mean and standard deviation.

[0113] Based on the above global statistics, the initial weights α, β, and γ are set as follows:

[0114]

[0115] Where ∈ is a small constant to prevent division by zero. By dynamically adjusting α, β, and γ, combined with the NDWI* index, the spectral contrast between water bodies and background features can be adaptively enhanced in different image regions, making the initial extraction of water objects more accurate and the boundaries clearer.

[0116] The image is further divided into several local windows (e.g., 32×32 or 64×64 pixels), and for each window, the local water body index contrast C is calculated. loc =μ water -μ background , where μ water and μ background These are the average reflectances of pixels representing water and non-water bodies, respectively. Then, the weights are optimized based on contrast.

[0117]

[0118] Where, σ local denoted as the joint standard deviation of the reflectance of the green band and NIR band within the local window, and k is an adjustment factor used to enhance the stability of the spectral response.

[0119] The weights are then iteratively updated for each local window to maximize the spectral difference between the water body and the background. The iterative formula is as follows:

[0120]

[0121] in, Let C be the contrast gradient, λ be the learning rate, and iterate until C is reached. loc Convergence or reaching the maximum number of iterations.

[0122] This method achieves adaptive weight adjustment, highlights the spectral characteristics of small water bodies, and suppresses pseudo-water body information such as vegetation, building shadows, and soil. Using the Otsu method, it establishes a water body / non-water body segmentation threshold, determines the initial water body region M, and further segments the water body region R into different connected water body regions based on the pixel connectivity principle. i .

[0123] (2) Extraction of high-confidence candidate regions for small water bodies based on geometric features

[0124] For each connected water body region R i Multidimensional geometric feature indices are introduced, including comprehensive compactness index, shape dispersion index, boundary fractal dimension index, radial variance index, etc., to more comprehensively characterize the morphological features of small water bodies, thereby improving the accuracy of small water body extraction.

[0125] Comprehensive tightness index C i for

[0126]

[0127] This index reflects the compactness of the region's shape. When the region is an ideal circle, C... i =1; the narrower or more irregular the shape, the higher C. i The smaller the value.

[0128] Shape dispersion index D i This is used to measure the degree of deviation of the region boundary from the ideal ellipse fit. Let the region R be... i The length of the major axis is L i The length of the minor axis is W i Then the area of ​​its fitted ellipse is Defined as:

[0129]

[0130] When the region contour is close to the fitted ellipse, D i A smaller value indicates a larger value, while a larger value indicates a greater number of complex disturbances at the boundary.

[0131] Boundary fractal dimension exponent F i This is used to characterize boundary complexity. Fractal analysis of the region boundary is performed using box counting. Let N(∈) be the minimum number of boxes required to cover the boundary at scale ∈, defined as:

[0132]

[0133] Theoretically, the boundary dimension of a regular geometric shape is close to 1, while the boundary dimension of a region F is complex and highly tortuous. i Approaching 2.

[0134] Introducing the radial variance exponent V i This is used to reflect the symmetry and regularity of a water body region relative to its centroid. Specifically, let the centroid coordinates of the water body region to be analyzed be (x... c y c The set of boundary points of the region is Where n represents the total number of sampling points on the boundary of the water body region, that is, the number of discrete points used to characterize the boundary. For each boundary point (x j ,y j ), calculate its radial distance r from the centroid. j The formula is:

[0135]

[0136] Based on this, the radial variance exponent V i The definition is as follows.

[0137]

[0138] in, This represents the mean radial distance of all boundary points, i.e.

[0139] This index can quantitatively reflect the geometric regularity of a water body: when the shape of the water body approaches a circle, the radial distance from each boundary point to the centroid changes little, therefore V i The value is relatively small; when the water body area has an irregular shape or obvious protrusions, the radial distance distribution of each boundary point varies significantly, resulting in V i The value increases significantly. Using this index, the regularity of the shape of a region can be quantitatively evaluated during the screening of candidate regions for small water bodies and the constraint of geometric features, thereby further optimizing the extraction binning of water body boundaries.

[0140] By analyzing the above multidimensional geometric feature index {C i D i ,F i V i The present invention uses a combined discrimination method to set multiple threshold constraint rules in the initial screening stage: only when the candidate region simultaneously satisfies C i >T c D i <T d ,F i <T f V i <T v At that time, T c Td T f T v The threshold values ​​are determined using a probability histogram.

[0141] Through the above constraint mechanism, the areas that are ultimately identified as high-confidence small water bodies are effectively removed, thereby effectively eliminating areas with unreasonable shapes or those caused by false water bodies such as shadows or stagnant water.

[0142] (2) Location of candidate areas for small water bodies under the joint constraints of multiple indicators

[0143] After identifying high-confidence candidate areas for small water bodies, this invention designs a multi-index joint constraint mechanism, including a comprehensive density index C. i Boundary complexity index B i and the shape regularity index R i It is used to quantify the spatial morphological characteristics of small water bodies, thereby enabling accurate identification of candidate areas for small water bodies.

[0144] 1. Comprehensive density index C i

[0145] The Comprehensive Density Index measures the compactness of the morphology within a region, effectively distinguishing between bodies of water and pseudo-water bodies such as narrow shadows or roads. It is defined as follows:

[0146]

[0147] Among them, A i Representing region R i area, P i Representing region R i The boundary length. For circular or approximately regular bodies of water, C i The value is close to 1; for elongated or irregularly shaped regions, C i The value decreased significantly.

[0148] 2. Boundary complexity index B i

[0149] The boundary complexity index is used to characterize the tortuosity of water body boundaries and can effectively identify regions with high geometric irregularities. It is defined as follows:

[0150]

[0151] The denominator is related to the area A. i The same circumference of a circle. If the shape of a body of water is nearly circular, then B i ≈1; If the water body boundary is complex or elongated, then B i The value increased significantly.

[0152] 3. Shape regularity index Ri

[0153] The shape regularity index is used to describe the regularity of the morphology of water bodies, and is particularly suitable for distinguishing between regular small water bodies and naturally formed irregular shapes. It is defined as follows:

[0154]

[0155] Among them, L i Representing region R i spindle length, W i R represents the length of the secondary axis. When the water body is nearly circular or square, R... i The value is close to 1; when the region is a long and narrow strip, R i The value is close to 0.

[0156] Based on the above three indices, a joint discriminant function is constructed:

[0157]

[0158] Where α, β, γ are weight parameters, satisfying α + β + γ = 1. When F(R i )≥T f At that time, region R i It has been identified as a candidate area for small water bodies.

[0159] This multi-index joint constraint model (joint discriminant function) enables precise segmentation of candidate regions for small water bodies, significantly improving the accuracy and reliability of small water body extraction.

[0160] In this embodiment, based on lightly processed multi-channel wave image data, a scale-weighted adaptive fusion algorithm based on water body geometric similarity is proposed to address the problem of small water body extraction. The specific steps are as follows:

[0161] 1. Targeting a specific region C in the image i (Obtained through selection based on geometric constraints), calculate its geometric similarity index on images at different resolutions. Set the image scale to L. σ Where σ∈{high resolution, medium resolution, low resolution}, geometric similarity can be defined as a measure of the similarity between the candidate region and the shape of the real water body boundary. For example, normalized Hausdorff distance or shape features can be used:

[0162]

[0163] in, For the candidate region boundary, For boundary estimation of water bodies at scale σ, d H (·) represents the Hausdorff distance. The higher the geometric similarity, the closer the value is to 1.

[0164] 2. Based on this similarity index, an adaptive weighting formula is defined for data at different scales:

[0165]

[0166] Through this weight ω σ The influence of images at different scales on the extraction of small water bodies is dynamically adjusted to maintain high-resolution accuracy while also considering the smoothness and edge suppression effect of low-resolution images. Next, a three-stage self-evolutionary refinement extraction strategy is adopted for candidate regions, implementing boundary extraction and refinement in stages: 1. Coarse-grained stage: Based on weighted fused images at various scales, preliminary water body discrimination is performed on candidate regions, generating coarse-grained boundaries. Stable positioning mainly depends on spectral thresholds and geometric constraints. The formula is expressed as:

[0167]

[0168] Among them, f pixel (x,L σ ) represents the water spectral determination function for pixel x at scale σ, and τ0 is the coarse localization threshold. 2. Refining Stage: Based on the coarse localization, the boundary is refined using feature constraints and geometric constraints of the candidate region. Through optimization of local gradient, edge energy, and geometric consistency, a refined boundary is generated.

[0169]

[0170] Among them, E spectral E is the spectral energy term. geometric E represents the geometric shape error term. spatial λ1, λ2, and λ3 are the energy terms constrained by spatial continuity, and λ1, λ2, and λ3 are the weighting coefficients.

[0171] 3. Subsequent Stages: Based on the refined results, and combining time-series data consistency constraints and small circle noise removal, the boundaries are refined at the sub-pixel level to finally obtain the boundary.

[0172]

[0173] Among them, E final For the combined spectral-spatial-geometric energy, E temporal For timing consistency constraints, E noiseThe small circle noise term is represented by α and β, which are adjustment coefficients. In this invention, the establishment process of the spectral-spatial-geometric feature joint decision model is as follows: Spectral features: Extracting spectral values ​​such as multi-channel band reflectance and single water body indices (ENDW1, MNDWI, etc.) to distinguish between water and non-water body pixels. Spatial features: Including neighborhood relationships, connectivity, and shape similarity, used to improve boundary continuity and consistency between regions. Geometric features: Using quantitative indicators such as shape moments, compactness, length, and brightness to constrain candidate regions and prevent misidentification of non-water bodies. The joint decision model is expressed in the form of a geometrically constrained loss function as follows:

[0174] F(x) = λ4f pixel (x)+λ5f spatial (x)+λ6f geometric (x);

[0175] Among them, λ4, λ5, and λ6 are weighting coefficients.

[0176] Pixel classification and boundary delineation are performed by minimizing threshold energy. The advantages of the model are that it not only uses a single feature for judgment, but also uses multiple features for joint decision-making, which significantly improves the extraction accuracy and boundary location accuracy of small water bodies. At the same time, it overcomes the problem of inaccurate identification of large and small water bodies by single-scale and single-feature extraction methods in existing technologies, and achieves sub-pixel level boundary refinement.

[0177] In this embodiment, the adaptive multi-scale fusion algorithm accurately extracts small water bodies:

[0178] To address the mismatch in identification accuracy of small water bodies caused by differences in spatial resolution among multi-scale remote sensing data, this study employs multi-scale image modeling, adaptive optimization, and a scale weighting mechanism to enhance the fine delineation of small water body boundaries. The specific implementation process includes:

[0179] (1) Multiscale image construction and normalization processing

[0180] First, multi-scale construction and normalization are performed on the preprocessed multi-source remote sensing image data. Let the high-resolution image be I. H (x,y), the medium-resolution image is I M (x,y), the low-resolution image is I L (x,y), where (x,y) represents the spatial pixel coordinates of the image; I H The spatial resolution is typically less than 10m, I M The spatial resolution is in the range of 10-30m, I L The spatial resolution is greater than 30m.

[0181] To achieve multi-scale analysis, the pyramid decomposition method is used to downsample the images at different resolutions, constructing image representations at different scale levels. Specifically, for any scale s∈{H,M,L}, the multi-scale representation of the image is denoted as I. (s) (x,y) is obtained by using the downsampling operator DownSample(·):

[0182] I (s) (x,y)=DownSample(I s (x,y),s);

[0183] Here, DownSample(·) represents the spatial downsampling operation, used to generate a lower-resolution image layer, and s represents the image scale level, corresponding to high resolution H, medium resolution M, and low resolution l, respectively. Due to radiometric differences between multi-source images, a minimum-maximum normalization method is used to standardize each scale image to eliminate spectral biases between different images. The normalized image is denoted as . The calculation formula is as follows:

[0184]

[0185] Wherein, min(I (s) ) and max(I (s) ) respectively represent image I (s) The minimum and maximum values ​​across all pixels. Through this normalization process, grayscale values ​​from images of different scales and sources are mapped to the range of [0,1], eliminating the influence of radiation differences on subsequent identification of small water bodies.

[0186] (2) Enhancement of spectral features of small water bodies based on adaptive band combination optimization

[0187] To address the sensitivity of small water bodies to spectral characteristics, this invention proposes an adaptive band combination optimization method based on entropy weighting to construct an improved small water body index W. s (x,y), this index is achieved through a weighted combination of the Normalized Differential Water Index (NDWI) and the Modified Normalized Differential Water Index (MNNDWI). Specifically, it is represented as follows:

[0188]

[0189] The parameters are defined as follows: W s (x,y) represents the improved micro-water body index for pixel (x,y) at scale s, used to measure the intensity of water body features at that pixel. G (s) (x,y) represents the green band reflectance value of pixel (x,y) at scale s; NIR (s)(x,y) represents the near-infrared reflectance value of pixel (x,y) at scale s, used to calculate the Normalized Differential Water Index (NDWI); SWIR (s) (x,y) represents the shortwave infrared reflectance value of pixel (x,y) at scale s, used to calculate the improved normalized differential water index (MNNDWI); α and β represent the band combination weight coefficients, which are adaptively determined by the entropy weight method, used to adjust the contribution of different bands in the identification of small water bodies to satisfy α+β=1; s represents the image scale level, representing different levels of the high, medium and low resolution image pyramid, used for multi-scale feature extraction.

[0190] In the weight determination process, this invention introduces the entropy weight method to fully utilize the dispersion and contribution capacity of each band's information. Weights w1 and w2 correspond to the relative importance of NDWI and MNNDWI, thereby achieving adaptive optimization of band combinations.

[0191] (α * ,β * ) = (w1, w2);

[0192] Where, α * β represents the band combination coefficient of the optimized NDWI, used to enhance the significance of water body spectral characteristics. Its value is calculated by the entropy weight method, reflecting the relative contribution of NDWI to distinguishing water bodies from non-water bodies; * The band combination coefficients of the optimized MNNDWI are used to supplement the identification ability of small water body edges and complex backgrounds. The values ​​are calculated by the entropy weight method and reflect the relative contribution of MNNDWI to the water body extraction box degree.

[0193] The final improved small water body index W s (x,y) can maximize the spectral distinction between water bodies and non-water bodies, improve the accuracy of identifying small water bodies, and is applicable to remote sensing image processing at different scales and under complex backgrounds.

[0194] (3) Construction of multi-scale adaptive water index based on geometric constraints

[0195] Obtaining water body indices W at various scales H (x,y), W M (x,y), W L After (x,y), this invention first combines the geometric features of the candidate water body region to construct a geometric feature quantification index (including density index, boundary complexity, etc.), providing a constraint basis for subsequent multi-scale fusion. Simultaneously, a water body probability model P(x,y) is introduced to quantify the probability that each pixel is a water body, thereby enhancing the ability to distinguish the edges and complex shapes of small water bodies. Its complete expression is:

[0196]

[0197] Among them, W s (x,y) represents the water index value at scales s∈{H,M,L}, α is a coefficient adjusting the steepness of the curve, and β is the water body discrimination threshold. This probabilistic model quantifies the probability of each pixel being identified as a water body into a probability value between 0 and 1, effectively improving the recognition accuracy for complex boundaries and small water bodies.

[0198] Introducing a density index Where A is the area of ​​the water body, P is the length of the water body boundary, and the closer C is to 1, the more regular the shape, and the boundary complexity. This is used to reflect the tortuosity and morphological complexity of water body boundaries. Based on the above geometric characteristics, a geometric consistency index Q is defined at scale s. s :

[0199] Q s =f(C s B s ,P s );

[0200] Where s∈{H,M,L} represents the high, medium, and low resolution image scales; Cs, Bs, and Ps represent the density index, boundary complexity, and boundary length of the water body region at scale s, respectively; the function f(Cs,Bs,Ps) can be specifically expanded as:

[0201]

[0202] Where α, β, and γ are empirical adjustment coefficients used to balance the contributions of shape regularity, boundary complexity, and area influence to the geometric consistency index. Based on the geometric consistency index, an adaptive scale weight allocation function is defined:

[0203]

[0204] Where s′ represents the scale index corresponding to s, s′∈{H,M,L}; -Qs′ represents the geometric consistency index at scale s′; λ is a semi-constant parameter used to adjust the sensitivity of the weight distribution. The final result is the multi-scale adaptive water index W. f The expression for (x,y) is:

[0205] W f (x,y)=∑ s∈{H,M,L} ω s ·W s (x,y);

[0206] The integrated results automatically enhance the contribution of high-resolution imagery to water body boundary regions to ensure boundary positioning accuracy. By combining geometric feature constraints and water body probability models, the enhancement of water body spectral features and information fusion under multi-scale, multi-source imagery are achieved. The Otsu method is used to determine the water body and non-water body segmentation threshold (segmentation threshold = the optimal threshold T*T^*T* automatically determined on the multi-scale adaptive water index map using the Otsu method, used for segmenting water body and non-water body pixels). Under the constraint of the candidate region range of small water bodies, small water bodies are automatically extracted and generated into a vector object set N.

[0207] In this embodiment, the boundary optimization of small water bodies is performed using both morphological and temporal consistency methods:

[0208] To address the issues of inaccurate boundary treatment, incomplete morphology, and interference from mixed pixels in small water bodies, a boundary optimization method based on a functional energy model is constructed. This method, combined with morphological processing, connected component analysis, and temporal consistency verification, achieves high-precision optimization of small water body boundaries. The implementation process includes the following steps:

[0209] (1) Initial boundary generation and active contour model construction

[0210] For a set of N small water body vector objects, an active contour model is used to refine the boundaries of each small water body object to achieve sub-pixel-level boundary localization. The specific steps are as follows:

[0211] First, for small water body vector objects, an initial boundary curve C(s) = (x(s), y(s)) is generated, where s∈[0,1] represents the parameterized position of the curve; x(s) represents the horizontal coordinate function of the curve in the two-dimensional image coordinate system, and its value range is the continuous real number interval of the image column coordinates, which can achieve sub-pixel level position description between pixel grids; y(s) represents the vertical coordinate function of the curve in the two-dimensional image coordinate system, and its value range is the continuous real number interval of the image row coordinates, which can also achieve sub-pixel level precise positioning; the initial boundary can be fitted by the water body contour in the coarse positioning stage or the region boundary obtained by filtering based on geometric features, ensuring that the initial contour covers the main water body area.

[0212] Secondly, the functional energy of the activity profile curve C(s) is defined as follows:

[0213] E(C)=∫0 1 (α|C'(s)| 2 +β|C””s)| 2 )ds+λ∫0 1 P(C(s))ds;

[0214] Wherein: α is the boundary smoothing control parameter, used to suppress boundary noise and ensure overall contour smoothness; β is the curvature adjustment parameter, used to maintain boundary continuity and prevent excessive bending; λ is the external energy weight, controlling the guiding effect of image information on contour evolution; P(C(s)) is the image gradient potential function, constructed based on the grayscale or spectral feature gradient of water body pixels, which can attract the contour to converge towards the real water body boundary; based on this, the present invention defines C′(s).C″(s) as follows:

[0215] First derivative vector:

[0216]

[0217] This represents the tangent vector of the boundary curve at the parameter s, used to describe the local directional changes of the curve.

[0218] Second derivative vector:

[0219]

[0220] This represents the rate of change of curvature of the boundary curve at parameter s, and is used to control the curvature and continuity of the curve.

[0221] In this invention, the refined extraction of the boundary of small water bodies is achieved by iteratively minimizing the functional energy E(C). Specifically, the active contour curve evolves gradually under the drive of the energy functional, causing the initial boundary to gradually approach the boundary of the real water body. The basic principle is as follows: First, the functional energy...

[0222] E(C)=∫0 1 (α|C′(s)| 2 +β|C″(s)| 2 )ds+λ∫0 1 Pp(C(s))ds;

[0223] By performing variational differentiation, the gradient direction of energy descent is obtained. By deriving the evolution equation of curve C(s) with respect to the time parameter τ, we can obtain:

[0224]

[0225] Expanding it into a specific form:

[0226]

[0227] Where -αC″(s) is a first-order smoothing constraint term used to suppress high-frequency noise at the boundary; βC″′(s) is a second-order curvature adjustment term used to maintain curve continuity and avoid excessive bending. As an external potential energy term, it guides the curve to gradually converge toward the boundary where the image gradient is strongest.

[0228] In practical calculations, this invention employs a discrete iteration strategy, representing the coordinates of the boundary curve C(s) at the k-th iteration as (x... k (s),y k (s)), then the update formula for the (k+1)th iteration is:

[0229] x k+1 (s)=x k (s)+Δτ·F x (x k (s),y k (s));

[0230] y k+1 (s)=y k (s)+Δτ·F y (x k (s),y k (s));

[0231] Where Δτ is the iteration step size; F x With F y These are the components of the gradient descent direction in the horizontal and vertical directions, respectively, determined by the variational results of the energy functional mentioned above.

[0232] (2) Small-area noise removal and connected region merging

[0233] Based on boundary refinement, a connected component labeling algorithm is used to analyze small water bodies. Let the set of connected components of small water bodies in the extracted results be {R}. i}, i=1,2,…,N, its area is A(R i Apply a threshold operation to each connected region:

[0234] R i ∈Ω if and only if A(R) i )≥A min ;

[0235] Among them, A min The minimum effective water body area threshold (e.g., 3 pixel units) is used; regions smaller than this threshold are considered noise and removed. Ω represents the set of retained effective small water body connected regions, i.e., all connected regions that satisfy the area threshold.

[0236] Ω={R i |A( R i)≥A min ,i=1,2,…,N;

[0237] For connected regions that do not meet the area threshold This is considered noise and removed. Furthermore, for adjacent or geometrically close connected regions, a distance threshold D can be used...merge Merging operations are performed to avoid boundary fragmentation and improve the continuity and integrity of small water body extraction:

[0238] R i R j ∈Ω, if d(R) i R j )≤D merge Then R i ∪R j →R ij ;

[0239] Wherein d(R) i R j ) represents region R i With R j The shortest boundary distance, R ij This represents the new connected region after merging. Additionally, for small water bodies that are spatially adjacent and have boundaries narrower than the threshold δ, a region merging operation is performed:

[0240]

[0241] (3) Optimization of temporal consistency of small water bodies

[0242] In long-term remote sensing monitoring, to ensure the stability, continuity, and boundary consistency of the extraction results for small water bodies, this invention proposes a dynamic optimization mechanism based on temporal consistency constraints. Let R be the extraction result for a certain region at time t. t The extraction result at time t+1 is R. t+1 To evaluate the consistency of water extraction results at different time points, a temporal similarity constraint function is introduced, which integrates spectral similarity, geometric similarity, and spatial neighborhood continuity.

[0243]

[0244] in:

[0245] 1. Spectral similarity index

[0246]

[0247] It represents the proportion of spatial overlap of water bodies at adjacent times, and is used to measure the spectral consistency of a region.

[0248] 2. Geometric consistency index G(R) t ,R t+1 ) is defined as:

[0249]

[0250] Among them, C(R) tF(R) represents the complexity of the water body boundary (such as the ratio of perimeter to area). t ) represents shape regularity (such as shape compactness or roundness). This index measures the degree to which the boundary complexity and shape of a water body remain consistent at different times. The value ranges from [0,1]. The larger the value, the more consistent the geometry.

[0251] 3. Spatial neighborhood continuity index N(R) t ,R t+1 ) is defined as:

[0252]

[0253] Among them, Neighborhood (R t+1 ) represents R t+1 The set of neighboring pixels of the region, 1(·) is an indicator function, which counts the consistency of adjacent pixels in the water body region in space, reflecting the continuity of the water body boundary and spatial extension.

[0254] 4. The adaptive weighting coefficients α, β, γ are dynamically determined based on the water body type, area, morphology, and background features to balance the importance of spectral, geometric, and spatial neighborhood under different temporal conditions.

[0255] When S(R) t ,R t+1 When the threshold θ is less than the threshold, it indicates that there is significant temporal inconsistency in the region, requiring dynamic optimization and correction. This invention utilizes historical results {R1, R2, ..., R...} t-1} and the current result R t The weighted multiple-decision mechanism has been updated:

[0256] R t ′=weighted-mode{w1R1,w2R2,...,w t-1 R t-1 ,w t R t};

[0257] Where weighted-mode(·) represents the weighted majority voting function, and the weight w t It can be dynamically allocated based on the reliability of historical time series, cloud shadow occlusion, and detection confidence.

[0258] This optimization mechanism can effectively correct false water body extraction caused by noise, cloud shadows, seasonal interference, or sudden events, while ensuring the boundary continuity and area stability of small water bodies in long-term series monitoring, thus achieving high precision and high reliability in monitoring dynamic changes of water bodies.

[0259] This method achieves dynamic adaptive optimization of water body extraction results, significantly improving the reliability and accuracy of long-term monitoring of small water bodies. Through the above optimization process, the final output of small water body extraction results is obtained, where water bodies are represented by 1 and non-water bodies by 0, and stored in raster format.

[0260] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made to the technical solutions of the present invention by those skilled in the art without departing from the spirit of the present invention should fall within the protection scope defined by the claims of the present invention.

Claims

1. A method for refined extraction of small-scale land surface water bodies based on multi-scale remote sensing, characterized in that, The method includes: Using multi-source satellite remote sensing data, a high-medium-low resolution image pyramid is constructed to achieve adaptive band combination optimization based on the spectral characteristics of small water bodies, and multi-scale noise suppression and radiometric normalization processing are performed. Based on preprocessed multi-source remote sensing image data, a geometric shape quantification index for small water bodies is preset. Based on the geometric shape feature quantification index for small water bodies, a small water body screening mechanism constrained by geometric features is constructed to achieve coarse localization of candidate areas for small water bodies. Based on preprocessed multi-source remote sensing image data, a scale weight adaptive allocation algorithm based on the similarity of water body geometric features is adopted. The contribution weight of data at different scales is dynamically adjusted according to the geometric features of candidate regions. A joint discrimination model of spectral-spatial-geometric features of small water bodies is established to achieve accurate extraction of small water bodies. A boundary optimization method based on a functional energy model is constructed, and the results of small-area noise removal, connected region merging, and time-series data constraint consistency verification are combined to achieve boundary optimization of small water bodies.

2. The method according to claim 1, characterized in that, The method for constructing a high-, medium-, and low-resolution image pyramid using multi-source satellite remote sensing data, achieving adaptive band combination optimization based on the spectral characteristics of small water bodies, and performing multi-scale noise suppression and radiometric normalization processing includes: The image pixel values ​​are unified to a dimensionless reflectance space using a radiometric normalization method. The processing method is as follows: Among them, R i DN is the normalized reflectivity value for the i-th band; i The original image pixel values; DN dark and DN white These are the pixel values ​​for dark and bright targets, respectively; A multi-scale pyramid structure is adopted to unify the imagery into three levels: high-resolution layer ≤2m, medium-resolution layer 5–30m, and low-resolution layer ≥30m. Spatial resampling uses bicubic interpolation, with the following formula: Where f(x,y) represents the pixel value at image coordinates (x,y) after resampling; DN(x+i,y+j) represents the pixel value at (x+i,y+j) in the original image; where i and j are relative coordinate offsets, ranging from -1, 0, 1, 2, corresponding to the 16 surrounding neighboring pixels; w(i,j) is the bicubic interpolation weight function, whose value is determined by the pixel offset distance, which can smooth the calculation results and preserve local details; coordinates (x,y) are the grid coordinates after resampling of the target image, with the upper left corner of the image as the origin, x representing the column direction and y representing the row direction; offset (i,j) defines the relative position of the neighboring pixels participating in interpolation in the original image, where i is the horizontal offset and j is the vertical offset; Wavelet decomposition and reconstruction methods are applied to images of different scales, with the following formula: I'(x,y)=∑ m ∑ n W m,n (x,y)·ψ m,n (x,y); Where I'(x,y) is the filtered image, ψ m,n (x,y) are wavelet basis functions with scale m and direction n, W m,n (x,y) represents the corresponding wavelet coefficients; Cloud and cloud shadow detection combining spectral and temporal features: The Fmask method based on spectral thresholds is used to automatically identify cloud and cloud shadow regions, obtaining the cloud mask matrix C(x,y); this is combined with time-series imagery T... t The spectral consistency constraint of (x,y) is used to interpolate data for pixels in cloud-occluded areas. The interpolation formula is as follows: I″(x,y,t)=α(x,y,t)·I′(x,y,t-1)+(1-α(x,y,t))·I′(x,y,t+1); Where I″(x,y,t) represents the image value of the cloud-occluded area after interpolation at time t; I′(x,y,t-1) and I′(t,y,t+1) represent the denoised image values ​​at adjacent times; α(x,y,t) is a weighting factor whose value is dynamically determined based on the spectral similarity of pixels in the time series and the time interval, and can adaptively adjust the contribution of the images before and after to the interpolation result.

3. The method according to claim 1, characterized in that, Based on the geometric shape quantification index of small water bodies, a screening mechanism for small water bodies constrained by geometric features is constructed. Methods for coarsely locating candidate regions of small water bodies include: Among them, F(R) i ) is the joint discriminant function, C i For the comprehensive density index, B i R is the boundary complexity exponent. i Let F(R) be the shape regularity index, and let α, β, γ be the weighting parameters, satisfying α + β + γ = 1; when F(R) i )≥threshold T f At that time, region R i It has been identified as a candidate area for small water bodies.

4. The method according to claim 1, characterized in that, A scale-weighted adaptive allocation algorithm based on water body geometric similarity is adopted. This algorithm dynamically adjusts the contribution weights of data at different scales according to the geometric features of candidate regions, establishing a joint discrimination model of spectral, spatial, and geometric features for small water bodies. The method for accurately extracting small water bodies includes: Based on the preprocessed multi-source remote sensing image data, denoted as High-Resolution Image I H (x,y), medium resolution image I M (x,y) and low-resolution image I L (x,y) is used to form different scale levels through pyramid decomposition: I (s) (x,y)=DownSample(I(x,y),s),s∈{H,M,L}; Where DownSample(·) represents the downsampling operator, and s is the scale level; Standardization is performed using the min-max normalization method: An improved small water body index W is constructed based on an adaptive band combination optimization method using entropy weighting. s (x,y), this index is achieved by a weighted combination of the Normalized Differential Water Index (NDWI) and the Improved Normalized Differential Water Index (MNNDWI), as shown below: Among them, W s (x,y) represents the improved micro-water body index of pixel (x,y) at scale s, used to measure the intensity of water body features of that pixel. G (s) NIR (s) With SWIR (s) The green band, near-infrared band, and shortwave infrared band represent the bands at scale s, respectively. Parameters α and β are the band combination weights used to adjust the contribution of each band to water body identification. In the weight determination process, the entropy weight method is introduced, with weights w1 and w2 corresponding to the relative importance of NDWI and MNNDWI, to achieve adaptive optimization of the band combination. (a * ,b * )=(w1,w2); Where, α * β represents the band combination coefficient of the optimized NDWI, used to enhance the significance of water body spectral characteristics. Its value is calculated by the entropy weight method, reflecting the relative contribution of NDWI to distinguishing water bodies from non-water bodies; * The band combination coefficients of the optimized MNNDWI are used to supplement the identification ability of small water body edges and complex backgrounds. The values ​​are calculated by the entropy weight method and reflect the relative contribution of MNNDWI to the water body extraction box degree. After obtaining water body indices at various scales, a geometric feature quantification index is first constructed by combining the geometric features of candidate water body regions. Simultaneously, a water body probability model P(x,y) is introduced to quantify the probability that each pixel is a water body. Furthermore, a density index and boundary complexity are introduced to define a geometric consistency index Q at scale s. s Based on the geometric consistency index, an adaptive scale weight allocation function is defined, ultimately yielding the multi-scale adaptive water index W. f (x,y): W f (x,y)=∑ s∈{H,M,L} ω s ·W s (x,y); Among them, W f (x,y) is the multi-scale adaptive water index, ω s s is the adaptive weighting function, where s is the scale level.

5. The method according to claim 1, characterized in that, Methods for constructing boundary optimization based on functional energy models include: Wherein, E(C) is the functional energy of the active contour curve C(s); α is the boundary smoothing control parameter, used to suppress boundary noise and ensure overall contour smoothness; β is the curvature adjustment parameter, used to maintain boundary continuity and prevent excessive bending; λ is the external energy weight, which controls the guiding effect of image information on contour evolution; P(C(s)) is the image gradient potential function, constructed based on the grayscale or spectral feature gradient of water body pixels, which can attract the contour to converge towards the real water body boundary; C′(s) is the first derivative vector, representing the tangent vector of the boundary curve at parameter s, used to describe the local direction change of the curve; and C″(s) is the second derivative vector, representing the rate of change of curvature of the boundary curve at parameter s, used to control the degree of curvature and continuity of the curve.

6. The method according to claim 1, characterized in that, Methods for small-area noise removal and connected region merging include: Based on boundary refinement, a connected component labeling algorithm is used to analyze small water bodies. Let the set of connected components of small water bodies in the extracted results be {R}. i }, i=1,2,…,N, its area is A(R i ); Perform a threshold operation on each connected region: R i ∈Ω if and only if A(R) i )≥A min ; Where Ω represents the set of effective connected regions of small water bodies that are retained, that is, all connected regions that satisfy the area threshold, and A min A minimum effective water area threshold is set; areas smaller than this threshold are considered noise and removed. Simultaneously, for spatially adjacent small water bodies with boundaries narrower than the threshold δ, a region merging operation is performed.

7. The method according to claim 1, characterized in that, Methods for verifying the consistency of time-series data constraints include: Let R be the extraction result of a certain region at time t. t The result at time t+1 is R. t+1 A temporal similarity constraint function is constructed by integrating spectral similarity, geometric similarity, and spatial neighborhood consistency. Among them, G(R) t ,R t+1 N(R) represents the consistency index between boundary complexity and shape regularity. t ,R t+1 ) represents the spatial continuity index of the neighboring area, and α, β, and γ are adaptive weighting coefficients determined according to the water body type, area, shape, and background features. When S(R) t ,R t+1 When the value is less than the threshold θ, it indicates that there is temporal inconsistency in this region, which requires dynamic optimization and correction.

8. The method according to claim 7, characterized in that, When S(R) t ,R t+1 When the threshold θ is less than the threshold θ, it indicates that there is temporal inconsistency in this region. Methods for dynamic optimization and correction include: Historical results {R1,R2,…,R} are used t-1 } and the current result R t The weighted multiple decisions have been updated: R t ′=weighted-mode{w1R1,w2R2,...,w t-1 R t-1 ,w t R t }; Where weighted-mode(·) represents the weighted majority voting function, and the weight w t The allocation is dynamically determined based on the reliability of historical time series, cloud shadow occlusion, and detection confidence.

Images