Multi-source data-driven flood storage and detention area ecological toughness detection method and system

Abstract

The invention belongs to the technical field of image processing, and particularly relates to a multi-source data-driven flood storage and detention area ecological toughness detection method and system, and the method comprises the steps: 1, carrying out the geometric correction and radiation correction of an optical satellite remote sensing image, and forming a time sequence basic data set; 2, a first axis and a second axis of the ecological state matrix correspond to space grids of the flood storage and detention area respectively, and a third axis corresponds to time; 3, mapping the ecological state matrix to a fractal space to form a time coupling matrix; 4, for each time slice of the time coupling matrix, setting a neighborhood with a fixed size around each pixel, and performing time window segmentation and fractal complexity estimation on the multi-scale morphological gradient map; and step 5, calculating the toughness value of the connected region according to the area and the toughness values of all positions in the connected region. According to the method, the accuracy and practicability of ecological toughness detection are remarkably enhanced.

Description

Technical Field

[0001] The present invention belongs to the field of image processing technology, and in particular relates to a multi-source data-driven method and system for detecting the ecological resilience of flood storage and detention areas. Background Art

[0002] Against the backdrop of the increasingly pressing dual demands for flood regulation and ecological protection, flood storage and detention areas, as crucial reservoirs during flood season, not only fulfill flood control functions by reducing flood peaks and mitigating downstream disasters, but are also increasingly being entrusted with the ecological function of restoring and maintaining regional ecosystem stability. Particularly in the context of frequent extreme weather events and increased hydrological fluctuations in river basins, comprehensively evaluating the ecological resilience of flood storage and detention areas before and after flood disturbances has become a crucial interdisciplinary research topic, encompassing water resources, ecology, remote sensing, and geography. Ecological resilience is generally defined as the ability of an ecosystem to recover its structure and function after a sudden disturbance. Its assessment relies on a detailed characterization of the disturbance response, the system's recovery process, and spatial variability. In this context, remote sensing technology, with its wide-area coverage, high timeliness, and non-contact observational capabilities, has become a crucial data source for assessing ecosystem changes before and after flood events, particularly the evolution of water distribution and vegetation status.

[0003] Existing research on analyzing ecological changes in flood storage and detention areas primarily focuses on the following technical approaches: First, single- or dual-temporal change detection methods based on optical remote sensing data. For example, these methods utilize vegetation and water-related indices such as the Normalized Difference Vegetation Index (NDVI) and the Normalized Water Index (NDWI) to extract spatial distribution maps before and after floods. These methods then analyze the areas of change through methods such as image differencing and change thresholding. While computationally simple, these methods lack the ability to respond to continuous temporal changes and cannot effectively capture the dynamic response and internal structural changes of the system. Second, methods employ time series remote sensing analysis. For example, using multi-temporal remote sensing data from MODIS and Sentinel-2, continuous NDVI is used to analyze ecosystem disturbances and recovery processes. However, these methods often rely on pixel-level averaging or trendline fitting, lacking the ability to capture multi-scale structural changes and making it difficult to capture the characteristics of localized, drastic responses to flood disturbances. Third, distributed models or physical process models are used to simulate flood processes and overlay land use or vegetation type data to analyze ecological impacts. However, these models rely on a large number of hydrological parameters and field data, resulting in poor regional adaptability and difficulty in widespread application. Fourth, local studies use manual field surveys and low-altitude drone imagery for detailed analysis. While this approach offers high accuracy, it is costly, time-consuming, and lacks the ability to dynamically monitor large areas. Summary of the Invention

[0004] In view of this, the main purpose of the present invention is to provide a multi-source data-driven method and system for detecting the ecological resilience of flood storage areas. By collecting optical remote sensing images before, during, and after a flood, extracting water masks and vegetation indices to construct a three-dimensional ecological state matrix, and mapping it to fractal space to enhance the ability to express spatiotemporal structures, the method then calculates the temporal coupling matrix and multi-scale morphological gradient map, and estimates the fractal complexity within a sliding time window, ultimately completing pixel-level resilience value labeling and connected region resilience evaluation. This method can comprehensively characterize the ecological response and recovery process under flood disturbances, improve the ability to identify and quantify multi-scale and temporal changes in ecosystems, and significantly enhance the accuracy and practicality of ecological resilience detection.

[0005] The technical solution adopted in the present invention is as follows: A multi-source data-driven method for detecting ecological resilience of flood storage and detention areas, the method comprising: Step 1: Collect optical satellite remote sensing images at different times before, during, and after a flood in the flood storage area; perform geometric and radiometric corrections on the optical satellite remote sensing images to form a basic time series dataset; Step 2: Extract two types of indicators from the time series basic dataset: water mask and vegetation index. Stack the obtained water mask and vegetation index in chronological order to form a three-dimensional ecological state matrix. The first and second axes of the ecological state matrix correspond to the spatial grid of the flood storage and detention area, and the third axis corresponds to time. Step 3: Map the ecological state matrix to a fractal space. Use a hierarchical iterative fractal function to perform self-similar transformation and scale compression on each layer of data. Within the fractal space, couple the fractal patterns at adjacent time points and calculate the similarity of adjacent fractal patterns using a sliding window to form a temporal coupling matrix. Step 4: For each time slice of the temporal coupling matrix, a fixed-size neighborhood is set around each pixel, and the difference between the maximum and minimum values ​​of the indicator in the neighborhood is calculated to obtain the local difference; the local differences of different scales are merged to generate a multi-scale morphological gradient map; the multi-scale morphological gradient map is segmented into time windows and the fractal complexity is estimated; Step 5: By comparing the fractal complexity of the same position over time, mark the toughness value of each position; then perform connected domain analysis on all positions to obtain multiple connected regions; for each connected region, calculate the toughness value of the connected region based on its area and the toughness values ​​of all the positions in it.

[0006] Furthermore, the optical satellite remote sensing images are multispectral images with visible light and near-infrared bands, and are acquired with the same revisit period at different time periods before, during and after the flood; the spatial resolution of the optical satellite remote sensing images is no more than ten meters.

[0007] Furthermore, in step 2, in the optical satellite remote sensing image of each time slice, a fixed threshold range is set based on the fact that water bodies have lower visible light reflectance and higher near-infrared absorption characteristics compared to surrounding objects, and the positions where the pixel reflectance falls within the threshold range are determined to be water areas; and morphological opening and closing operations are performed on the water area determination results to eliminate isolated noise and fill small-scale holes, thereby obtaining a binary water area mask; for each pixel, its near-infrared band reflectance and red band reflectance are obtained, the difference between the two is calculated, and then the sum of the two is calculated, and the ratio of the difference to the sum is used as the vegetation index of the pixel.

[0008] Furthermore, in step 2, the water mask adopts the 8-adjacent connectivity rule, retaining only water patches with an area of ​​not less than 3 pixels, and performing morphological operations of dilation and erosion on isolated patches to reduce pseudo-water noise; all optical satellite remote sensing images are reprojected and unified to the UTM coordinate system based on the WGS-84 ellipsoid, with a fixed pixel size of 10m×10m; the first and second axes of the ecological status matrix correspond to the row and column numbers of the UTM grid, respectively, and the third axis is arranged in order from early to late according to the image acquisition time.

[0009] Furthermore, in step 3, the process of mapping the ecological state matrix to the fractal space includes: setting the iteration depth , the value range is 4 to 6; for each pixel vector in the ecological state matrix Execute the fractal mapping function as follows: ;in, , for Scale compression matrix, is a two-dimensional translation vector; The corresponding pixel grayscale is written into the same position in the fractal space to complete one layer of iteration; the above steps are repeated for all time slices in sequence until a complete multi-layer fractal pattern stack is generated; the scale compression matrix of each layer is Take diagonal form ,in ; Translation vector According to the row and column number of the pixel Calculated as ;in, is the row number, is the column number; is the transpose operation.

[0010] Furthermore, in the fractal space, for adjacent time points and When coupling the fractal pattern of The sliding window of pixels has a window step size of 1 pixel. is an integer subscript index; for each window position , is the coordinate value of the first axis of the window, The coordinate value of the second axis of the window; calculate the two-frame window sub-block and The structural similarity index : ,in and are the average grayscale of the two window sub-blocks, are the grayscale variances of the two window sub-blocks, is the covariance of the two window sub-blocks, , ; For all window positions The structural similarity index Perform pixel-level aggregation to obtain the coupling coefficients of adjacent time points ,in is the total number of valid pixels covered by the window; Fill in the two-dimensional matrix in time order , we get the time coupling matrix, whose row index and column index correspond to the time points ;Time coupling matrix Stored in a symmetrical form, with values ​​of 1 assigned to the diagonal positions, the matrix elements satisfy ,in ; The time limit; is an integer subscript index.

[0011] Furthermore, in step 4, for each time slice of the time coupling matrix, a square neighborhood with a side length of 3 pixels is established for each pixel in the matrix, and the difference between the maximum and minimum index values ​​in the neighborhood is calculated to obtain the local difference map of the first scale; in step 4, in addition to the first scale, square neighborhoods with side lengths of 5 pixels and 7 pixels are further used to calculate the local difference in the same way as the first scale, and the local difference maps of the three scales are weighted linearly superimposed according to weights 0.5, 0.3 and 0.2 to obtain the fused multi-scale morphological gradient map.

[0012] Furthermore, in step 4, for each time slice of the time coupling matrix, a square neighborhood with a side length of 3 pixels is established for each pixel in the matrix, and the difference between the maximum and minimum values ​​of the index value in the neighborhood is calculated to obtain the local difference map of the first scale; in addition to the first scale, square neighborhoods with side lengths of 5 pixels and 7 pixels are further used to calculate the local difference in the same way as the first scale, and the local difference maps of the three scales are weighted linearly superimposed with weights of 0.5, 0.3 and 0.2 to obtain a fused multi-scale morphological gradient map; pixel-by-pixel maximum selection is performed on the weighted multi-scale morphological gradient map, that is, the maximum value of the local differences of the three scales is selected at the same pixel position to generate the final multi-scale morphological gradient map.

[0013] Furthermore, the final multi-scale morphological gradient image sequence is segmented along the time dimension using a sliding time window of length 5, with a window step of 1 time slice, so that there is an overlap of 4 time slices between adjacent time windows to ensure the temporal continuity of the fractal complexity estimation; for the multi-scale morphological gradient image in each time window, the side length of the observation grid is gradually reduced to 32 pixels, 16 pixels and 8 pixels, and the number of non-zero grids occupied by the gradient image at different levels of observation grid size is counted. The observation grid at each level and the corresponding number of non-zero grids are plotted in a logarithmic coordinate system. The slope of the line is obtained by least squares straight line fitting, and the absolute value of the slope is taken as the fractal complexity value of the time window.

[0014] A multi-source data-driven flood storage and detention area ecological resilience detection system, the system includes: a data acquisition unit, used to collect optical satellite remote sensing images at different time periods before, during and after a flood in the flood storage and detention area; geometric correction and radiation correction are performed on the optical satellite remote sensing images to form a time series basic data set; an ecological state matrix construction unit, used to extract two types of indicators from the time series basic data set, namely: water mask and vegetation index; the water mask and vegetation index obtained are stacked in chronological order to form a three-dimensional ecological state matrix; the first and second axes of the ecological state matrix correspond to the spatial grid of the flood storage and detention area, respectively, and the third axis corresponds to time; a time coupling matrix construction unit, used to map the ecological state matrix to a fractal space, and use a hierarchical iterative fractal function to perform self-similar transformation and scale compression on each layer of data; in In the fractal space, the fractal patterns of adjacent time points are coupled, and the similarity of adjacent fractal patterns is calculated through a sliding window to form a time coupling matrix; the fractal complexity calculation unit is used to set a neighborhood of a fixed size around each pixel for each time slice of the time coupling matrix, calculate the difference between the maximum and minimum values ​​of the indicator in the neighborhood, and obtain the local difference; the local differences of different scales are merged to generate a multi-scale morphological gradient map; the multi-scale morphological gradient map is segmented into time windows and the fractal complexity is estimated; the toughness assessment unit is used to mark the toughness value of each position by comparing the fractal complexity of the same position over time; then a connected domain analysis is performed on all positions to obtain multiple connected areas; for each connected area, the toughness value of the connected area is calculated based on its area and the toughness values ​​of all positions therein.

[0015] By adopting the above technical solutions, the present invention produces the following beneficial effects: it can achieve standardized processing and index extraction of multi-temporal remote sensing images before, during, and after a flood, forming a three-dimensional ecological state matrix with a time-sequential structure. This matrix is ​​then mapped to a fractal space through a hierarchical iterative fractal mapping function, effectively enhancing the self-similar expression of spatial texture and structural features. A structural similarity evaluation mechanism is introduced into the fractal space to construct a temporal coupling matrix that captures the degree of structural change between adjacent time slices, thereby revealing the temporal continuity of the disturbance and recovery process. In terms of spatial structure analysis, a multi-scale morphological gradient map is generated through multi-scale neighborhood difference calculation, and sliding window segmentation and fractal complexity estimation are performed in the time dimension to achieve a scale-fused expression from local boundary response to overall complexity evolution. Finally, by comparing the changes in fractal complexity on the time axis, the pixel-level ecological resilience value is labeled, and combined with connected domain analysis, the resilience statistics of the connected area are output, achieving a complete transition from pixel-level dynamic monitoring to regional-level comprehensive evaluation. This invention breaks through the existing technology's reliance on single-scale changes, static differences, and single-indicator responses, significantly improving the interpretability and analytical accuracy of the ecological response process of flood storage and detention areas under flood disturbances, and providing a more scientific, stable, and high-resolution support means for ecological restoration assessment, post-disaster compensation accounting, and ecological management of flood storage and detention areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 A schematic diagram of the system structure of a multi-source data-driven flood storage and detention area ecological resilience detection system provided in an embodiment of the present invention; Figure 2 A schematic diagram of the scale change of the observation grid for fractal complexity estimation provided by an embodiment of the present invention; Figure 3 A schematic diagram of the connected domain analysis principle provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0017] All features disclosed in this specification, or all steps in the disclosed methods or processes, except mutually exclusive features and / or steps, can be combined in any manner.

[0018] Unless otherwise stated, any feature disclosed in this specification (including any appended claims and abstract) may be replaced by other equivalent or similar features. In other words, unless otherwise stated, each feature is only an example of a series of equivalent or similar features.

[0019] refer to Figure 1 A multi-source data-driven method for detecting ecological resilience of flood storage and detention areas, the method comprising: Step 1: Collect optical satellite remote sensing images at different times before, during, and after a flood in the flood storage area; perform geometric and radiometric corrections on the optical satellite remote sensing images to form a basic time series dataset; Step 1 aims to construct a time series foundational dataset for cross-temporal comparison of the multi-source data-driven flood storage and detention area ecological resilience assessment method. This approach involves collecting optical satellite remote sensing images from the flood storage and detention area at different time periods before, during, and after a flood. These images are then geometrically and radiometrically corrected to ensure consistent spatial positioning and comparable radiometric measurements for the same features. This ensures temporal and spatial consistency in the subsequent extraction of water mask and vegetation indices. Specifically, the target area is determined based on the management boundary of the flood storage and detention area. Observation periods covering the pre-flood, mid-flood, and post-flood phases are selected based on flood records. Optical satellite remote sensing images in the visible and near-infrared bands are prioritized, and efforts are made to ensure stable revisit periods, similar imaging geometry, and manageable seasonal variations to minimize interference from non-target factors on the time series features. To improve data availability, candidate images are first screened for quality and subjected to cloud shadow masking. Images with high cloud cover, significant banding noise, or sensor anomalies are eliminated, retaining only those that demonstrate complete spatial coverage, balanced temporal distribution, and continuity across key flood stages. Subsequently, geometric correction is carried out. By using orbital attitude information, sensor imaging geometry, terrain undulation information and known ground control points, the images are geo-referenced and distortion corrected, so that optical satellite remote sensing images of different time periods can be aligned at the pixel level within the flood storage area, ensuring that the row and column positions of the same landform in different time slices remain consistent; when necessary, the local mismatched areas are refined and the residual evaluation is performed until the spatial error meets the pixel-level consistency requirements.

[0020] Next, radiation correction is performed to restore the sensor-recorded observations to physically meaningful surface reflectance. Factors such as atmospheric scattering and absorption, differences in solar altitude and azimuth angles, and sensor response heterogeneity are corrected to make the grayscale values ​​between different time slices comparable. To further stabilize the time series characteristics, the images are uniformly calibrated and normalized to ensure that the visible and near-infrared bands maintain consistent radiation scales at different time periods. After completing geometric and radiation corrections, all optical satellite remote sensing images are cropped based on the boundaries of the flood storage area, retaining only the target area data. The time axis is then sorted, and the slices before, during, and after the flood are organized from early to late according to the acquisition time. For time periods with gaps, qualified images from nearby dates are given priority to supplement to maintain the continuity of the time series. On the basis of the above processing, all optical satellite remote sensing images that have undergone geometric correction and radiation correction are organized and recorded according to unified naming and metadata specifications, clarifying the acquisition time, solar geometry, observation geometry and quality identification corresponding to each time slice. Ultimately, a time series basic dataset with a clear structure, reliable source and cross-period comparability is formed. This lays the foundation for data consistency and traceability for the subsequent extraction of water masks and vegetation indices from the time series basic dataset and stacking them in chronological order to form a three-dimensional ecological status matrix.

[0021] Step 2: Extract two types of indicators from the time series basic dataset: water mask and vegetation index. Stack the obtained water mask and vegetation index in chronological order to form a three-dimensional ecological state matrix. The first and second axes of the ecological state matrix correspond to the spatial grid of the flood storage and detention area, and the third axis corresponds to time. The principle of step 2 is to simultaneously extract water masks and vegetation indices from the time series basic dataset and stack the obtained water masks and vegetation indices in chronological order to form a three-dimensional ecological state matrix. The first and second axes correspond to the spatial grid of the flood storage and detention area, and the third axis corresponds to time. This supports subsequent fractal spatial mapping, temporal coupling, and multi-scale morphological gradient analysis. In the specific implementation process, the optical satellite remote sensing image quality is first reviewed for each time slice in the time series basic dataset to ensure that geometric correction and radiometric correction are consistent. The image is then unified into the UTM coordinate system based on the WGS-84 ellipsoid. The pixel size is fixed at 10m×10m, and the target range is obtained by overlapping and clipping with the flood storage and detention area boundary. Then, within each time slice, a fixed threshold range is set to judge the pixels based on the lower visible light reflectance and higher near-infrared absorption characteristics of water bodies compared to surrounding objects, and a preliminary binarization result is obtained. To suppress isolated noise and small-scale voids, morphological opening and closing operations are performed in sequence, and the 8-adjacent connectivity rule is used to retain only water patches with an area of ​​not less than 3 pixels. At the same time, the boundaries of the retained patches are refined using the morphological operations of dilation and erosion, and a water mask with spatial coherence, continuous boundaries, and stable topological relationships is output.

[0022] The vegetation index is calculated based on the difference between the reflectance of the near-infrared band and the reflectance of the red band. For each pixel, the reflectance values ​​of the two bands are first obtained, and the difference and sum of the two are calculated. The ratio of the difference to the sum is then used as the vegetation index for that pixel. Consistent scaling and invalid value identification are performed across the entire panorama to ensure a one-to-one correspondence with the water mask in both time and space. After completing single-phase processing, the water mask and vegetation index for the same time slice are aligned pixel-wise using row and column numbers, and a data consistency check is performed to eliminate inconsistent pixels caused by residual cloud shadows, sensor striping, or local mismatches. Data locations confirmed to be invalid are uniformly coded to maintain consistent statistical caliber in the temporal dimension in subsequent analysis. All time slices were then organized from morning to night according to image acquisition time. The water mask and vegetation index for each time slice were paired and stacked on the grid to construct a three-dimensional ecological state matrix with row number as the first axis, column number as the second axis, and time as the third axis. During the matrix generation process, the acquisition time and quality indicator of each time slice were recorded to ensure that the matrix could be traced back to the original slice in the time series basic dataset. This transformed the time series basic dataset into a three-dimensional ecological state matrix that combines information on both water body and vegetation changes. This provides structured input and a stable indicator source for mapping the ecological state matrix into a fractal space, coupling fractal patterns at adjacent time points within the fractal space, and estimating multi-scale morphological gradients and segmented fractal complexity on the temporally coupled matrix.

[0023] Step 3: Map the ecological state matrix to a fractal space. Use a hierarchical iterative fractal function to perform self-similar transformation and scale compression on each layer of data. Within the fractal space, couple the fractal patterns at adjacent time points and calculate the similarity of adjacent fractal patterns using a sliding window to form a temporal coupling matrix. The principle of step 3 is to map the three-dimensional ecological state matrix to a fractal space slice by slice without changing the spatial topology of the ecological state matrix. Using a hierarchical iterative fractal function, a self-similar transformation and scale compression are performed on each layer of data to extract the cross-scale structural features formed by the water mask and vegetation index driven by flood processes in the flood storage area. The fractal patterns of adjacent time points are coupled within the fractal space, and the similarity of adjacent fractal patterns is calculated using a sliding window to form a temporal coupling matrix. In the specific implementation process, each time slice of the three-dimensional ecological state matrix is ​​read from early to late, treated as a layer of input, and the mapping result of the fractal space is initialized. The iteration depth range is set to 4 to 6. In each iteration, the positions of all pixels in the time slice are self-similarly transformed and scale compressed according to the hierarchical iterative fractal function. This allows the texture and boundaries represented by the water mask and vegetation index in the same time slice to reappear at smaller scales, thereby gradually generating a fractal pattern stack with obvious self-similar characteristics.

[0024] To ensure spatial consistency, the mapping process maintains a one-to-one correspondence between the first and second axes of the ecological state matrix and the corresponding positions in the fractal space. This ensures that the mapped fractal patterns are strictly aligned with the original pixels in terms of row and column numbers. After each layer is completed, the fractal pattern for that layer is written to the same position in the fractal space until all iterations for that time slice are complete. This process is then repeated for all time slices, forming a complete multi-layer fractal pattern stack. After the fractal space mapping is completed, two layers of fractal patterns at adjacent time points are selected as a pair of inputs in the fractal space. A sliding window of 11 × 11 pixels is used to traverse the first and second axes with a step of 1 pixel. At each window position, the similarity of the two fractal patterns is calculated to measure the structural consistency and difference between adjacent time points within that local range. After the traversal, the similarities obtained for all window positions are summarized at the pixel level to obtain the coupling coefficient for the pair of adjacent time points. These coupling coefficients are then added to the corresponding positions in a two-dimensional matrix in chronological order, gradually constructing the temporal coupling matrix. The temporal coupling matrix is ​​organized symmetrically, with diagonal values ​​assigned 1. This allows it to reflect the coupling strength of fractal patterns between adjacent time points and, through its overall structure, to depict the stage-by-stage differences and continuity before, during, and after a flood. By mapping the ecological state matrix into a fractal space, utilizing a hierarchical iterative fractal function to enhance self-similarity, and calculating the similarity of adjacent fractal patterns within the fractal space using a sliding window, the temporal coupling matrix can stably and objectively represent the spatial structural changes of the flood storage area at different time slices. This provides a cross-scale consistent input foundation for subsequent neighborhood calculations at each time slice, the generation of multi-scale morphological gradient maps, and the temporal window segmentation and fractal complexity estimation.

[0025] Step 4: For each time slice of the temporal coupling matrix, a fixed-size neighborhood is set around each pixel, and the difference between the maximum and minimum values ​​of the indicator in the neighborhood is calculated to obtain the local difference; the local differences of different scales are merged to generate a multi-scale morphological gradient map; the multi-scale morphological gradient map is segmented into time windows and the fractal complexity is estimated; The principle of step 4 is to simultaneously characterize the structural mutation and evolutionary continuity of the flood storage area on both spatial and temporal lines based on the time coupling matrix. First, a neighborhood of fixed size is established with the pixel as the center in each time slice. The local difference is obtained by calculating the difference between the maximum and minimum values ​​of the indicators in the neighborhood. The local differences of different scales are then merged to generate a multi-scale morphological gradient map. The multi-scale morphological gradient map is then segmented into time windows and the fractal complexity is estimated. In this way, the spatiotemporal differences of the water mask and vegetation index changes triggered by the flood process after coupling in the fractal space are stably transcribed into comparable intensity and complexity metrics. In the specific implementation process, in each time slice of the temporal coupling matrix, a square neighborhood with a side length of 3 pixels is established for each pixel in the matrix. The maximum and minimum values ​​of the indicators in the neighborhood are extracted and the difference between the two is calculated to obtain the local difference map of the first scale. To improve the perception of spatial heterogeneity at different scales, square neighborhoods with side lengths of 5 pixels and 7 pixels are further used for the same time slice. The local difference maps of the second scale and the third scale are obtained in the same way as the first scale. The local difference maps of the three scales are weighted linearly superimposed with weights of 0.5, 0.3, and 0.2 to output the fused multi-scale morphological gradient map. The result of the pixel-by-pixel maximum selection is retained as another fusion method, that is, the maximum value of the local differences of the three scales is directly selected at the same pixel position to form the final multi-scale morphological gradient map. The two methods maintain consistent data organization and quality identification in the subsequent process, so as to facilitate method comparison and robustness testing in different study areas and different flood processes.

[0026] After the generation of the multi-scale morphological gradient map of a single temporal phase is completed, the multi-scale morphological gradient map sequence is segmented along the time dimension using a sliding time window of length 5, with a window step of 1 time slice, so that adjacent time windows overlap by 4 time slices to ensure temporal continuity and complete coverage of sudden changes. Within each time window, observation grids with side lengths of 32 pixels, 16 pixels, and 8 pixels are set in sequence by gradually reducing the observation grid. The number of non-zero grid cells occupied by the multi-scale morphological gradient map is counted for each level of observation grid. The observation grid size at each level and the corresponding number of non-zero grid cells are formed into a scattered point set in the logarithmic coordinate system. The value representing the fractal complexity of the time window is obtained by least squares linear fitting. The absolute value of this value is used as the fractal complexity estimation result of the time window and recorded at the index position corresponding to the time window. Through the above process from local differences to multi-scale morphological gradient maps, and then to time window segmentation and fractal complexity estimation, it is possible to stably extract strength indicators that are sensitive to the stage-by-stage structural differences before, during and after the flood on each time slice of the time coupling matrix, and at the same time, extract the cross-scale complexity change trajectory within the sliding time window, providing complete and consistent input for subsequent marking of resilience values ​​through time comparison of fractal complexity at the same location, and calculating the resilience values ​​of connected areas in connected domain analysis.

[0027] Step 5: By comparing the fractal complexity of the same position over time, mark the toughness value of each position; then perform connected domain analysis on all positions to obtain multiple connected regions; for each connected region, calculate the toughness value of the connected region based on its area and the toughness values ​​of all the positions in it.

[0028] The principle of step 5 is to use the fractal complexity of the multi-scale morphological gradient map sequence as a stable characterization in the time dimension, and to perform a temporal comparison of the changes in fractal complexity at the same location in consecutive time windows before, during, and after the flood, to obtain a resilience value that can reflect the disturbance response and recovery process. In space, adjacent locations are aggregated into connected areas through connected domain analysis, and a comprehensive assessment of the connected areas is conducted by combining the area and the resilience values ​​of all locations therein to form a regional-level resilience expression for management and decision-making. In the specific implementation process, first, an ordered index is established for each location on the time series of fractal complexity, and paired comparisons are performed from early to late according to the time window. The change amplitude during the flood relative to before the flood and the recovery amplitude after the flood relative to during the flood are identified respectively, and the directionality and continuity of the change are used as the basis for marking the resilience value of the location; in order to reduce the impact of occasional noise, time series smoothing and outlier detection are first performed on the time series of fractal complexity, and the breakpoints caused by residual cloud shadows or invalid data are filled and quality marked to ensure that the time comparison of the same location is based on the available slices; then, according to the key time period of the flood process, the disturbance peak and recovery inflection point are extracted on the time axis of the same location, and the decline amplitude and recovery speed in several windows after the disturbance are used to jointly characterize the resilience characteristics, and they are standardized and mapped to pixel-level resilience values ​​so that different locations can be compared at the same scale.

[0029] After completing the pixel-level resilience value labeling, a connectivity analysis is performed on all locations in the entire flood storage area, and adjacent locations with valid resilience values ​​are aggregated according to the 8-adjacent connectivity rule to obtain multiple connected areas. During the aggregation process, minimum area constraints and morphological corrections are performed on isolated small blocks with small areas, fragmented morphologies or caused by noise to ensure the spatial continuity and boundary rationality of the connected areas. For each connected area, the area of ​​the connected area is first calculated, and then the resilience values ​​of all locations within the area are statistically summarized. The area and pixel-level resilience values ​​are used to jointly determine the resilience value of the connected area, so that connected areas with larger areas, higher internal resilience values ​​and more uniform spatial distribution can obtain higher stability expressions in the comprehensive evaluation. At the same time, the time index range and quality identification of the connected area are retained to facilitate horizontal comparison between different flood processes. The final output includes two types of results: pixel-level resilience value raster and connected area resilience value list. The former directly reflects the time comparison results of the same location, and the latter reflects the regional aggregation characteristics based on connected domain analysis. Both maintain consistency in row numbers, column numbers and time indexes with the time coupling matrix, multi-scale morphological gradient map and three-dimensional ecological status matrix, and can be traced back to the time series basic data set and optical satellite remote sensing images, realizing a closed-loop evaluation of the ecological resilience detection method of flood storage and detention areas driven by multi-source data from before the flood, during the flood to after the flood.

[0030] Furthermore, at the data source level, multispectral images with visible and near-infrared bands are selected as optical satellite remote sensing images, and are acquired with the same revisit period in three stages. On the one hand, this keeps the sampling interval on the time axis constant, reducing seasonal and sunlight angle differences caused by observation time differences. On the other hand, because optical satellite remote sensing images cover both visible and near-infrared bands, they can simultaneously depict water body and vegetation conditions at the same time, forming a dual-index system sensitive to processes such as flooding, water body expansion, and vegetation recovery after water recedes. In terms of spatial quality, products with a spatial resolution of no more than ten meters are preferred. Finer pixel scales are used to improve the recognition of small-scale water surfaces, narrow ditches, mudflat edges, and shoreline transition zones. This also improves the accuracy of vegetation index expression in the case of mixed ground object pixels. To ensure spatial consistency for cross-period comparisons, all optical satellite remote sensing images were reprojected after geometric and radiometric corrections were completed and unified into the UTM coordinate system based on the WGS-84 ellipsoid. The pixel size was fixed at 10 m × 10 m. By unifying the coordinate system and pixel scale, the slices at different imaging dates, orbits, and attitudes achieved a one-to-one correspondence between row and column numbers within the flood storage area, providing a stable spatial framework for the subsequent construction of a three-dimensional ecological status matrix.

[0031] During the index extraction phase, for each time slice, a fixed threshold range is set to identify pixels based on the low visible light reflectance and high near-infrared absorption characteristics of water bodies compared to surrounding terrain, initially identifying candidate water areas. Considering the potential for isolated noise and holes caused by cloud shadows, sensor strips, and terrain shadows, a morphological opening operation is first performed on the binarized result to eliminate discrete small spots, followed by a morphological closing operation to fill small holes, resulting in a spatially continuous binary water mask with smoother boundaries. To further suppress pseudo-water noise and ensure that water patches have regional significance, the 8-neighbor connectivity rule is used to analyze the connectivity of the water mask, retaining only water patches with an area of ​​at least 3 pixels. For any remaining isolated patches, morphological operations of dilation and erosion are performed based on scene characteristics to optimize patch morphology and topology without destroying the true boundaries. Synchronous with the water mask, the near-infrared band reflectance and the red band reflectance are obtained at each pixel location. The difference between the two is first calculated, and then the sum is calculated. Finally, the ratio of the difference to the sum is used as the vegetation index for the pixel. A consistent calculation process and invalid value identification strategy are maintained across all temporal phases to ensure the comparability of the vegetation index in the temporal dimension and its splicing in the spatial dimension. After completing the extraction of the water mask and vegetation index for a single temporal phase, the two are aligned at the pixel level at the raster level and quality checked. Individual inconsistent pixels caused by residual cloud shadows or geometric micro-mismatches are eliminated. The row and column number indexes determined during the reprojection phase of the optical satellite remote sensing image are inherited, so that each pixel has a traceable spatial identification at different time slices. The images are then arranged in order from morning to night according to their acquisition time, and the water masks and vegetation indices of all time slices are stacked in chronological order to form a three-dimensional ecological status matrix, in which the first and second axes correspond to the row and column numbers of the UTM grid, respectively, and the third axis corresponds to time, thereby embedding the spatiotemporal change information on water expansion, retreat, vegetation damage, and recovery into a structured data container.

[0032] In the subsequent process, this container is directly mapped into a fractal space as input. A hierarchical iterative fractal function is used to perform self-similar transformation and scale compression at each layer to enhance the texture boundaries and block patterns jointly shaped by the water mask and vegetation index. The similarity of the fractal patterns at adjacent time points in the fractal space is calculated using a sliding window to form a temporal coupling matrix. For each time slice of the temporal coupling matrix, the difference between the maximum and minimum values ​​of the index within a fixed-size neighborhood is calculated to obtain the local difference, which is then merged at multiple scales to generate a multiscale morphological gradient map. The multiscale morphological gradient map is then segmented along the time dimension using a sliding time window, and the number of non-zero grid cells is counted based on the step-by-step reduction of the observation grid to complete the fractal complexity estimation. Through the above bottom-up and interdependent process, on the one hand, the uniformity of optical satellite remote sensing images in revisit period, spatial resolution, coordinate reference and pixel scale is guaranteed; on the other hand, the consistency of water mask and vegetation index in judgment criteria, morphological processing, connectivity constraints and pixel-level alignment is ensured. Finally, with the three-dimensional ecological state matrix as a bridge, the original observations of multispectral images are stably converted into a structural representation that is comparable across time periods. This provides a reliable data foundation and a clear index system for subsequent temporal comparison of fractal complexity at the same location to mark resilience values, conduct connected domain analysis on all locations, and calculate the resilience values ​​of connected areas. The entire implementation emphasizes data consistency management and process traceability. The results of any time slice can be traced back to the original optical satellite remote sensing image and the corresponding reprojection, cropping and quality identification records, thereby meeting the stability and reusability requirements of the multi-source data-driven flood storage and detention area ecological resilience detection method in engineering deployment and cross-regional migration.

[0033] Furthermore, in the specific implementation process, each time slice is first read from the three-dimensional ecological state matrix in chronological order, the write cache of the corresponding layer of the fractal space is initialized in each processing, and the iteration depth N is set to a value range of 4 to 6 to ensure that the multi-level characterization of the local and overall structure can be achieved within a limited amount of calculation; then the layered iterative fractal function is applied to each pixel vector in the time slice one by one, and a 2×2 scale compression matrix is ​​used to isotropically scale the pixel position during the iteration process, and the scaling coefficient decreases as a power of 2, so that the high-level level converges quickly to the local structure; at the same time, a two-dimensional translation vector is applied to the position, and the translation amount is fractionally shifted by a power of 2 according to the row and column numbers of the pixel, thereby achieving a smooth expansion of the original spatial pattern while scaling, ensuring that the self-similar transformation neither destroys the traceability of the row and column numbers in the fractal space nor produces topological dislocation across the grid. At the end of each layer iteration, the fractal pattern obtained in the current layer is written to the same position in the fractal space according to the original row and column numbers, ensuring strict alignment of the first and second axes, and the correspondence between the layer and the original time slice is recorded; all layers are processed for the same time slice according to the set iteration depth to form a multi-layer fractal pattern for the time slice; the above operation is repeated for all time slices in sequence to finally generate a complete multi-layer fractal pattern stack.

[0034] To avoid boundary effects and invalid value propagation, a unified quality label is used for invalid pixels in the three-dimensional ecological state matrix during the mapping process and masking is performed during iterative writing. Mirroring or constant extension strategies are enabled for pixels close to the boundary to maintain the continuity of the fractal space on the first and second axes. To reduce the impact of numerical instability on subsequent similarity calculations, a unified grayscale standardization can be performed after each layer is completed to ensure that the grayscale dimensions of different time slices and different iterative layers are consistent. After completing the fractal space mapping, the process of coupling the fractal patterns of adjacent time points in the fractal space begins. Specifically, two adjacent fractal patterns on the time axis are taken as a pair of inputs. A sliding window of size 11×11 pixels and a step size of 1 pixel are used to traverse the entire image in the first and second axis directions. At each window position, the structural similarity index is calculated for the two-frame window sub-block as the object. The structural similarity index is jointly characterized by the average grayscale, grayscale variance and covariance of the two window sub-blocks to simultaneously consider the matching degree at the brightness, consistency and structural levels. Invalid pixels are strictly ignored during the calculation, and the results of the position with insufficient valid pixels in the window are eliminated according to the quality mark to avoid the offset of sparse data on the statistics. After the traversal is completed, the structural similarity indices obtained from all window positions are summarized at the pixel level to calculate the coupling coefficient of the pair of adjacent time points. The value range is limited to between 0 and 1 to intuitively represent the continuous degree from complete dissimilarity to complete identity.

[0035] The coupling coefficients of adjacent time points are entered into the corresponding entries of a two-dimensional matrix from early to late in time, gradually constructing a temporal coupling matrix. This matrix is ​​stored symmetrically, with diagonal positions strictly set to 1 to represent identical time points within the matrix. The row and column indices of the temporal coupling matrix correspond to time points t1 through tT, respectively. The matrix elements satisfy interchangeable symmetry and range constraints, facilitating direct access to time slices in subsequent multiscale morphological gradient calculations. To ensure comparability across different flood processes and study areas, a fixed window size and step size strategy of 11 × 11 and 1 are used when aggregating the coupling coefficients. Furthermore, quality metrics such as the total number of valid pixels M and the proportion of invalid pixels are retained for each pair of time slices, allowing for downweighting or elimination of low-quality entries in subsequent evaluations. At the implementation level, the coupling calculation can be performed using block-wise parallelization and vectorization to reduce computational overhead. For large-scale flood storage and detention area data, data can be sliced ​​by row and column number and seamlessly joined in the final stage, ensuring that each entry in the temporal coupling matrix can be traced back to a specific time pair and valid pixel count. Through the above process, the ecological state matrix is ​​mapped to a fractal space and the fractal pattern coupling of adjacent time points is completed. While maintaining spatial alignment and temporal order, the temporal coupling matrix condenses the cross-scale differences in texture, boundary and block pattern before, during and after the flood, providing a strictly consistent input for establishing a fixed-size neighborhood on each time slice, calculating local differences and synthesizing multi-scale morphological gradient maps; at the same time, the symmetrical structure and standardized value range of the temporal coupling matrix enable it to directly participate in the segmentation processing and fractal complexity estimation of the sliding time window, thereby supporting the subsequent steps of marking the resilience value by temporally comparing the fractal complexity at the same location, performing connected domain analysis and calculating the resilience value of the connected area, thereby realizing the feasibility and reusability of the multi-source data-driven ecological resilience detection method for flood storage and detention areas in the entire link.

[0036] Furthermore, within each time slice, a fixed-size neighborhood is established for each pixel in the matrix and the difference between the maximum and minimum values ​​of the index within this neighborhood is calculated. This local difference, which is sensitive to the abrupt boundary and local amplitude, is obtained. This process is then repeated at multiple scales, with weighted linear superposition or pixel-by-pixel maximum selection, to obtain a multiscale morphological gradient map that takes into account both details and the overall picture. This is then segmented along the temporal dimension using a sliding time window of length 5, with a window step of 1 time slice, ensuring an overlap of 4 time slices between adjacent time windows. Finally, within each time window, the observation grid is gradually reduced and the number of non-zero grid cells occupied by the gradient map is counted. The absolute value of the slope of the least-squares linear fit is used as the fractal complexity value. This transforms flood-driven structural changes into comparable complexity trajectories and provides robust input for subsequent resilience value labeling and connected domain analysis. Specifically, in each time slice of the temporal coupling matrix, a square neighborhood with a side length of 3 pixels is established for each pixel in the matrix.

[0037] All index values ​​in the neighborhood are traversed and the difference between the maximum and minimum values ​​is calculated to generate a local difference map of the first scale. This scale prioritizes capturing fine-grained boundaries and narrow texture changes. In order to improve the response to low-frequency structures corresponding to large-scale water body expansion, retreat, and slow recovery of vegetation index, square neighborhoods with side lengths of 5 pixels and 7 pixels are further established on the same time slice. The local difference is calculated in the same way as the first scale to obtain local difference maps of the second and third scales. After completing the local difference calculation of the three scales, in order to balance boundary sensitivity and regional stability in a unified expression, the local differences of the three scales are firstly calculated according to weights 0.5, 0.3, and 0.2. The weighted linear superposition of the local difference maps is performed to obtain a fused multi-scale morphological gradient map, and the fusion result is checked for invalid value inheritance and boundary completeness to ensure that no pseudo gradient is introduced in the processing near the edge of the image and at the position where there is a missing measurement; at the same time, the results of the pixel-by-pixel maximum selection are calculated in parallel on the same data, that is, the maximum value of the three-scale local difference is directly selected at the same pixel position to generate another final multi-scale morphological gradient map, which is used to emphasize the strongest response in the transition zone of complex terrain or the stage with strong sudden disturbance; the two fusion methods are completely consistent in data organization, time index and quality identification. The subsequent process can select one of them according to the research objectives or compare them in parallel to test the robustness.

[0038] After the generation of the multi-scale morphological gradient map of a single phase is completed, it is arranged into a sequence on the time axis in the order of acquisition, and is segmented along the time dimension using a sliding time window of length 5, with a window step of 1 time slice, so that there are 4 time slices overlapping between adjacent time windows. This setting improves the coverage probability of rapid fluctuations and staged platform periods while ensuring time continuity; for the multi-scale morphological gradient map in each time window, in order to avoid fine texture disturbances driven by noise directly at the pixel level, the strategy of gradually reducing the observation grid is adopted to extract cross-scale occupancy features. Specifically, the observation grids with side lengths of 32 pixels, 16 pixels and 8 pixels are set in turn, and the multi-scale morphological gradient map in the time window is segmented at each level of observation. Coverage statistics are performed on the observation grid, the number of non-zero grids occupied by the gradient map is recorded, and quality indicators such as the effective coverage ratio and the invalid pixel mask ratio are saved for each level to facilitate the subsequent removal of abnormal windows. Subsequently, the observation grid size of each level and the corresponding number of non-zero grids are formed into a scattered point set in the logarithmic coordinate system, and the slope of the line is obtained by least squares straight line fitting. The absolute value of the slope is taken as the fractal complexity value of the time window and written into the time index consistent with the start and end slices of the time window. To ensure the consistency of the method under different flood processes, the side length of the observation grid is set to be fixed as well as the window length and window step size, and windows with insufficient non-zero grid counts or too many invalid pixels are downgraded or marked as unavailable.

[0039] The entire implementation takes the time coupling matrix as input, first generating a multi-scale morphological gradient map at each time slice, and then robustly refining the fractal complexity with a sliding time window, thereby expressing the instantaneous intensity of the local difference between the maximum and minimum values ​​in the spatial neighborhood, and condensing it into a complexity metric that can be compared across scales and stages; in the flood storage area scenario, the square neighborhood with a side length of 3 pixels focuses on responding to narrow waterways, shoreline micro-migration and local vegetation damage, while the square neighborhoods with side lengths of 5 and 7 pixels are more sensitive to large-scale flooding of beaches, retreat of lakeside shallows and community-level recovery; weighted linear superposition with weights 0.5, 0.3 and 0.2 strengthens the priority focus on detailed boundaries while retaining the contribution of medium and large-scale morphology, and the pixel-by-pixel maximum selection provides a stronger disturbance indication during extreme events or sudden increases in image contrast. Through the above consistent, traceable and parameter-fixed process, the local differences in space of the temporal coupling matrix are systematically integrated into a multi-scale morphological gradient map, and a fractal complexity sequence is stably output in time through a sliding time window with a length of 5 and a step size of 1. This provides a continuous, comparable and noise-resistant spatiotemporal feature basis for the subsequent temporal comparison of the fractal complexity at the same location to mark the resilience value, and to combine the area and pixel-level resilience values ​​in the connected domain analysis to calculate the resilience value of the connected area. This ensures that the multi-source data-driven ecological resilience detection method for flood storage and detention areas has consistent measurement caliber and engineering feasibility at all levels from single pixels to connected areas.

[0040] Furthermore, in step 5, for the fractal complexity sequence formed along the time sequence of the same pixel position, the average fractal complexity of the first three time windows of the sequence is first taken as the baseline complexity value; then, in the subsequent time windows, the fractal complexity of the pixel is compared one by one. When the fractal complexity recovers to a threshold value not higher than the baseline complexity value plus 10 percent, the window number is recorded as the recovery time; if the recovery condition is still not met within 8 consecutive time windows, the pixel is marked as a low recovery position. In step 5, a discrete resilience level is assigned to each pixel based on the length of the recovery time. Pixels with a recovery time of no more than 5 time windows are assigned a resilience level of 5, pixels with a recovery time of 6 to 8 time windows are assigned a resilience level of 4, pixels with a recovery time of 9 to 12 time windows are assigned a resilience level of 3, pixels with a recovery time of 13 to 16 time windows are assigned a resilience level of 2, and pixels that have not recovered after more than 16 time windows are assigned a resilience level of 1. The eight-connected region marking algorithm was applied to all pixels with marked resilience levels to identify connected domains. Connected domains with an area of ​​less than 9 pixels were directly removed and excluded from subsequent resilience calculations. For the retained connected domains, the number of pixels in their area was recorded. For each retained connected domain, the arithmetic mean of the resilience levels of the pixels within the domain was calculated. This mean was then multiplied by the area ratio of the total number of pixels within the domain divided by the total number of grid cells to obtain a preliminary regional resilience value. If the regional resilience value was higher than 4 and the area exceeded 5 percent of the total grid area, the regional resilience value was adjusted upward by 0.5 to highlight the importance of large areas of high resilience. If the proportion of pixels marked by the water mask within any connected domain exceeded 30 percent of the regional area, the resilience value of the connected domain was reduced by 1. If the proportion was between 15 percent and 30 percent, the resilience value was reduced by 0.5. If the proportion was less than 15 percent, the resilience value was not adjusted. The resilience values ​​of all connected domains are arranged in descending order and labeled as "high resilience zone", "medium-high resilience zone", "medium resilience zone", "medium-low resilience zone" and "low resilience zone" according to the set rules; at the same time, a regional resilience distribution map consistent with the grid is output, and each pixel is filled with the final resilience value of the connected domain to which it belongs, which will be used for subsequent ecological restoration decisions in flood storage and detention areas.

[0041] To facilitate reproducibility, a rectangular experimental area within the flood storage area was selected, and the UTM coordinate system based on the WGS-84 ellipsoid was uniformly adopted. The pixel size was fixed to 10m×10m, and the row and column numbers were used as spatial indexes to construct a three-dimensional ecological state matrix and subsequent mapping and statistics. The grid size of the experimental area was set to (Line No. , column number ). The time points for collecting optical satellite remote sensing images are Before the Flood , in the flood , after the flood The images are multispectral images covering visible and near-infrared bands, with a spatial resolution of 10 meters and the same revisit period. Geometric and radiometric corrections were performed for all time phases, followed by reprojection to the UTM coordinate system and cropping to the test area.

[0042] For each time slice First, a fixed threshold is set according to the characteristics of water bodies with low reflection in visible light and high absorption in near-infrared to obtain a preliminary water body judgment. Then, morphological opening and closing operations are performed, and the 8-adjacent connectivity rule is used to retain only water patches with an area of ​​not less than 3 pixels, and a binary water mask is output. At the same time, the vegetation index is calculated (the ratio of the difference between the red band and the near-infrared band reflectance to their sum): .in, is the reflectivity in the near-infrared band, is the red band reflectivity, are the pixel row and column numbers, is the time index. Example: In the pixel Up, if ,but .Will and Stacked in chronological order, a three-dimensional ecological status matrix is ​​formed, with the first axis being the row number, the second axis being the column number, and the third axis being time.

[0043] Setting iteration depth (Specific values ​​within the range of 4 to 6). For each time slice, the two-dimensional field is mapped to a fractal space by layered iterative fractal function. For any pixel position vector (consisting of row and column numbers), define the iteration: ,in, for Scale compression matrix, For the Isotropic compressibility of the layer, is a two-dimensional translation vector, are pixel row and column numbers, Represents the transpose operation. The grayscale of the pixel at the corresponding position is written to the same position in the fractal space, and the writing of this layer is completed. Repeat for all time slices to obtain a multi-layer fractal pattern stack. The size is The sliding window is traversed with a step size of 1 to calculate the structural similarity index (SSIM form): ,in, is the coordinate of the upper left corner of the window on the first and second axes, are the average grayscale of the two time sub-blocks, is the variance, is the covariance, is the stability constant. The coupling coefficients of adjacent time points are obtained by pixel-level averaging: ,in is the total number of valid windows. Fill in the two-dimensional matrix in chronological order Form a time coupling matrix that satisfies , ,and Example: In a small area, take A window, if , we can get a single window by substituting . For the whole map The average of the valid windows is obtained .

[0044] For each time slice of the time coupling matrix, a square neighborhood with a side length of 3 pixels is established for each pixel in the matrix, and the difference between the maximum and minimum values ​​of the index in the neighborhood is calculated to obtain the local difference map of the first scale. The same method is used to construct the neighborhood with side lengths of 5 and 7, and we get and The fused multi-scale morphological gradient map is obtained by weighted linear superposition with weights of 0.5, 0.3, and 0.2: . Parallel retain pixel-wise maximum selection: Example: In the pixel At , if the local difference of the three scales is ,but .

[0045] Will (or ) Segment the image along the time dimension using a sliding time window of length 5 and step size 1. For each time window, perform fractal complexity estimation (box counting idea): use observation grids with side lengths of 32 pixels, 16 pixels, and 8 pixels to cover the entire image, and count the number of non-zero grids occupied by the gradient map, which is recorded as ,in For the Level observation grid side length, is the corresponding non-zero grid number. Perform the least squares straight line fitting and take the absolute value of the slope as the fractal complexity of the time window Example: to Window, if statistics are obtained , the slope of the least squares regression is ,but . Write by window center time index .

[0046] The fractal complexity of the same position is compared over time to define the pixel-level toughness value. The mean fractal complexity in the pre-flood window is , the maximum fractal complexity that occurs in a flood is , the mean fractal complexity of the post-flood stable period is The toughness value is defined as ,in, Indicates the degree of recovery relative to the pre-flood baseline, with 1 indicating full recovery to the baseline level and 0 indicating an unrecovered state equivalent to the peak disturbance. ,like , (appears in nearby), ,but .

[0047] Get the full image element level toughness value grid Then, we conduct the connected domain analysis. Using the 8-adjacent connected rule, we will Valid pixels are aggregated into several connected areas For each connected region, calculate the area ,in is the area of ​​a single pixel (10m×10m). The area-weighted mean of the toughness values ​​in the region is calculated as the toughness value of the connected region: ,in For the region Number of pixels. Example: If an area Number of pixels , the average pixel-level toughness value is , then the area , regional toughness value Similarly, we can get the List, as the resilience assessment result of the management unit.

[0048] 1) In the fractal space mapping stage, the iteration depth The choice of affects the balance between details and overall structure. It can balance the computational complexity and expressiveness under the data scale of this example. 2) When calculating the structural similarity index, the constant Used to avoid instability caused by the denominator being close to 0; when the ratio of valid pixels in the window is lower than the set threshold (for example ) can be excluded 3) In the fusion of multi-scale morphological gradient maps, the weight Emphasizes fine-scale mutation identification; if you are more concerned about the extreme disturbance stage, you can choose to use Enhance sensitivity. 4) When estimating fractal complexity, observe the grid edge length The time window length is 5 and the step size is 1 to ensure the comparability across periods. When the ratio of invalid pixels in a window exceeds the threshold, the window is closed. Marked as invalid and ignored in the toughness value calculation. 5) In the definition of pixel-level toughness value, and The mean of the window sets before and after the flood can be taken respectively; if and If the denominator is too close to 0, a minimum difference constraint should be applied or a percentile baseline should be used to maintain numerical stability. 6) The output of the connected domain analysis includes both pixel-level toughness value grids and pixel-level toughness value grids. , also includes a list of connected areas Both can be traced back to the row number, column number and time index of the time coupling matrix, multi-scale morphological gradient map and three-dimensional ecological state matrix to meet the needs of engineering auditing.

[0049] Through the above examples, a comprehensive, multi-source data-driven method for assessing ecological resilience in flood storage and detention areas has been implemented: starting with the standardized processing of optical satellite remote sensing imagery, index extraction, and construction of a three-dimensional ecological state matrix. This method characterizes temporal structural relationships through fractal spatial mapping and a temporal coupling matrix. Multi-scale morphological gradients and fractal complexity are used to characterize the intensity and complexity of spatiotemporal variation. Finally, a quantitative assessment of management units is performed using pixel-level resilience values ​​and connected region resilience values. The above parameters and formulas can be slightly adjusted based on regional characteristics and data quality, but the indexing system, calculation caliber, and quality control strategy remain consistent, allowing for migration and reuse across different flood storage and detention areas.

[0050] refer to Figure 2The figure illustrates the technical implementation principle for estimating observation grid scale variations based on fractal complexity in the present invention. This technical principle is the key step in further fractal complexity estimation after the multi-scale morphological gradient image is segmented into time windows in step 4. Specifically, for the multi-scale morphological gradient image within each time window, the present invention employs a technique of progressively reducing the observation grid to obtain fractal feature parameters. This process involves three different observation grid scale levels: a first-level observation grid with a side length of 32 pixels, a second-level observation grid with a side length of 16 pixels, and a third-level observation grid with a side length of 8 pixels. In the first-level 32-pixel observation grid, the entire multi-scale morphological gradient image is divided into 32-pixel square grids, and the number of non-zero grid cells occupied by the gradient image is counted. The left portion of the figure shows the spatial distribution of the grids at this level, where the black-filled areas indicate the non-zero grid locations occupied by the multi-scale morphological gradient image. By traversing each grid individually, the number of non-zero grid cells at this scale is determined to be 8. In the second-level 16-pixel observation grid, the same technical scheme is employed, but the grid side length is reduced to 16 pixels. Due to the reduction in grid size, more detailed spatial feature changes can be captured. Accordingly, the number of non-zero grid cells occupied by the gradient map increases to 18. The center portion of the figure shows the detailed structure of the grid at this level, observing the significant increase in grid density and the spatial details of the non-zero grid distribution. In the third-level eight-pixel observation grid, the observation grid side length is further reduced to eight pixels, at which point the grid resolution reaches the highest level. Due to the further reduction in grid size, more local morphological change features can be identified, and the number of non-zero grid cells occupied by the gradient map further increases to 32. The right portion of the figure shows the distribution characteristics of the grid at this level, with the grid density reaching the maximum value. After obtaining the three different levels of observation grid sizes and their corresponding non-zero grid numbers, the present invention plots these data points in a logarithmic coordinate system to form a scatter plot. The logarithmic coordinate system at the bottom of the figure illustrates the implementation process of this technology: the horizontal axis represents the logarithmic value of the observation grid side length, and the vertical axis represents the logarithmic value of the number of non-zero grid cells. The three data points correspond to the observation results of (8, 32), (16, 18), and (32, 8), respectively.

[0051] Figure 3 The figure is a principle diagram of the connected domain analysis technology, which shows in detail the core technology implementation process of the connected domain analysis in the present invention. Figure 3 The distribution of resilience values, shown on the left, illustrates the spatial organization of the pixel grid within the flood storage and detention area. This grid uses a unified UTM coordinate system, with a fixed pixel size of 10m x 10m. The black dots in the figure represent locations marked with resilience values ​​after temporal comparison. These locations are determined by comparing the fractal complexity of the same location over time. The spatial distribution of these resilience markers is non-uniform, reflecting the spatial heterogeneity of ecological resilience within the flood storage and detention area. Figure 3 The middle section details the technical implementation of the 8-neighbor connectivity rule. This rule establishes connectivity relationships between the center pixel and its eight surrounding pixels, including horizontal, vertical, and diagonal neighbors. The 8-neighbor connectivity rule ensures the integrity of the connectivity analysis and accurately identifies spatially adjacent toughness locations. The center pixel is filled in black, while the eight adjacent pixels are filled in white and labeled with their adjacency numbers, clearly demonstrating the spatial configuration of adjacency relationships. Figure 3 The right side shows the specific appearance of the connected region results. Connected region A contains three adjacent toughness value locations, distributed linearly; connected region B contains seven toughness value locations, forming a larger connected region; connected region C contains two toughness value locations, forming the smallest connected unit. Each connected region has independent spatial boundaries and internal structural characteristics.

[0052] Although specific embodiments of the present invention have been described above, those skilled in the art will appreciate that these specific embodiments are merely illustrative, and that those skilled in the art may omit, substitute, and modify the details of the methods and systems described above without departing from the principles and spirit of the present invention. For example, combining the above method steps to perform substantially the same functions in substantially the same manner to achieve substantially the same results falls within the scope of the present invention. Accordingly, the scope of the present invention is limited solely by the appended claims.

Claims

1. A multi-source data-driven method for detecting ecological resilience of flood storage and detention areas, characterized by: The method comprises: Step 1: Collect optical satellite remote sensing images at different times before, during, and after a flood in the flood storage area; perform geometric and radiometric corrections on the optical satellite remote sensing images to form a basic time series dataset; Step 2: Extract two types of indicators from the time series basic dataset: water mask and vegetation index. Stack the obtained water mask and vegetation index in chronological order to form a three-dimensional ecological state matrix. The first and second axes of the ecological state matrix correspond to the spatial grid of the flood storage and detention area, and the third axis corresponds to time. Step 3: Map the ecological state matrix to a fractal space. Use a hierarchical iterative fractal function to perform self-similar transformation and scale compression on each layer of data. Within the fractal space, couple the fractal patterns at adjacent time points and calculate the similarity of adjacent fractal patterns using a sliding window to form a temporal coupling matrix. Step 4: For each time slice of the temporal coupling matrix, a fixed-size neighborhood is set around each pixel, and the difference between the maximum and minimum values ​​of the indicator in the neighborhood is calculated to obtain the local difference; the local differences of different scales are merged to generate a multi-scale morphological gradient map; the multi-scale morphological gradient map is segmented into time windows and the fractal complexity is estimated; Step 5: By comparing the fractal complexity of the same position over time, mark the toughness value of each position; then perform connected domain analysis on all positions to obtain multiple connected regions; for each connected region, calculate the toughness value of the connected region based on its area and the toughness values ​​of all the positions in it.

2. The multi-source data-driven flood storage and detention area ecological resilience detection method according to claim 1, characterized in that: Optical satellite remote sensing images are multispectral images with visible light and near-infrared bands, and are acquired with the same revisit cycle at different times before, during and after the flood; the spatial resolution of optical satellite remote sensing images is no more than ten meters.

3. The multi-source data-driven flood storage and detention area ecological resilience detection method according to claim 2, characterized in that: In step 2, in the optical satellite remote sensing image of each time slice, a fixed threshold range is set based on the lower visible light reflectance and higher near-infrared absorption characteristics of water bodies compared to surrounding objects, and locations where the pixel reflectance falls within the threshold range are determined to be water areas; Morphological opening and closing operations are performed on the water area determination results to eliminate isolated noise and fill small holes, thereby obtaining a binary water area mask; for each pixel, its near-infrared band reflectance and red band reflectance are obtained, the difference between the two is calculated, and then the sum of the two is calculated, and the ratio of the difference to the sum is used as the vegetation index of the pixel.

4. The multi-source data-driven flood storage and detention area ecological resilience detection method according to claim 3, characterized in that: In step 2, the water mask uses the 8-adjacent connectivity rule to retain only water patches with an area of ​​at least 3 pixels. Morphological operations such as dilation and erosion are performed on isolated patches to reduce pseudo-water noise. All optical satellite remote sensing images are reprojected and unified to the UTM coordinate system based on the WGS-84 ellipsoid, with a fixed pixel size of 10m × 10m. The first and second axes of the ecological status matrix correspond to the row and column numbers of the UTM grid, respectively, and the third axis is arranged in order from early to late according to the image acquisition time.

5. The multi-source data-driven flood storage and detention area ecological resilience detection method according to claim 4, characterized in that: In step 3, the process of mapping the ecological state matrix to the fractal space includes: setting the iteration depth , the value range is 4 to 6; for each pixel vector in the ecological state matrix Execute the fractal mapping function as follows: ;in, , for Scale compression matrix, is a two-dimensional translation vector; The corresponding pixel grayscale is written into the same position in the fractal space to complete one layer of iteration; the above steps are repeated for all time slices in sequence until a complete multi-layer fractal pattern stack is generated; the scale compression matrix of each layer is Take diagonal form ,in ; Translation vector According to the row and column number of the pixel Calculated as ;in, is the row number, is the column number; is the transpose operation.

6. The multi-source data-driven flood storage and detention area ecological resilience detection method according to claim 5, characterized in that: In the fractal space, for adjacent time points and When coupling the fractal pattern of The sliding window of pixels has a window step size of 1 pixel. is an integer subscript index; for each window position , is the coordinate value of the first axis of the window, The coordinate value of the second axis of the window; calculate the two-frame window sub-block and The structural similarity index : ,in and are the average grayscale of the two window sub-blocks, are the grayscale variances of the two window sub-blocks, is the covariance of the two window sub-blocks, , ; For all window positions The structural similarity index Perform pixel-level aggregation to obtain the coupling coefficients of adjacent time points ,in is the total number of valid pixels covered by the window; Fill in the two-dimensional matrix in time order , we get the time coupling matrix, whose row index and column index correspond to the time points ;Time coupling matrix Stored in a symmetrical form, with values ​​of 1 assigned to the diagonal positions, the matrix elements satisfy ,in ; The time limit; is an integer subscript index.

7. The multi-source data-driven flood storage and detention area ecological resilience detection method according to claim 6, characterized in that: In step 4, for each time slice of the time coupling matrix, a square neighborhood with a side length of 3 pixels is established for each pixel in the matrix, and the difference between the maximum and minimum index values ​​in the neighborhood is calculated to obtain the local difference map of the first scale; in step 4, in addition to the first scale, square neighborhoods with side lengths of 5 pixels and 7 pixels are further used to calculate the local difference in the same way as the first scale, and the local difference maps of the three scales are weighted linearly superimposed with weights of 0.5, 0.3 and 0.2 to obtain the fused multi-scale morphological gradient map.

8. The multi-source data-driven flood storage and detention area ecological resilience detection method according to claim 7, characterized in that: In step 4, for each time slice of the time coupling matrix, a square neighborhood with a side length of 3 pixels is established for each pixel in the matrix, and the difference between the maximum and minimum index values ​​in the neighborhood is calculated to obtain the local difference map of the first scale; in addition to the first scale, square neighborhoods with side lengths of 5 pixels and 7 pixels are further used to calculate the local difference in the same way as the first scale, and the local difference maps of the three scales are weighted linearly superimposed with weights of 0.5, 0.3 and 0.2 to obtain the fused multi-scale morphological gradient map; the weighted multi-scale morphological gradient map is subjected to pixel-by-pixel maximum selection, that is, the maximum value of the local differences of the three scales is selected at the same pixel position to generate the final multi-scale morphological gradient map.

9. The multi-source data-driven flood storage and detention area ecological resilience detection method according to claim 8, characterized in that: The final multi-scale morphological gradient image sequence is segmented along the time dimension using a sliding time window of length 5, with a window step of 1 time slice, so that there is an overlap of 4 time slices between adjacent time windows to ensure the temporal continuity of the fractal complexity estimation; for the multi-scale morphological gradient image in each time window, the side length of the observation grid is gradually reduced to 32 pixels, 16 pixels, and 8 pixels, and the number of non-zero grids occupied by the gradient image at different levels of observation grid size is counted. The observation grid at each level and the corresponding number of non-zero grids are plotted in a logarithmic coordinate system. The slope of the line is obtained by least squares straight line fitting, and the absolute value of the slope is taken as the fractal complexity value of the time window.

10. A multi-source data-driven flood storage and detention area ecological resilience detection system for implementing the method according to any one of claims 1 to 9, characterized in that: The system comprises: a data acquisition unit for collecting optical satellite remote sensing images at different time periods before, during and after a flood in a flood storage area; geometric correction and radiation correction are performed on the optical satellite remote sensing images to form a time series basic data set; an ecological state matrix construction unit for extracting two types of indicators from the time series basic data set, namely, a water area mask and a vegetation index; the water area mask and the vegetation index obtained are stacked in time sequence to form a three-dimensional ecological state matrix; the first axis and the second axis of the ecological state matrix correspond to the spatial grid of the flood storage area, respectively, and the third axis corresponds to time; a time coupling matrix construction unit for mapping the ecological state matrix to a fractal space, and performing self-similar transformation and scale compression on each layer of data using a hierarchical iterative fractal function; and in the fractal space, the adjacent time layers are mapped to the ecological state matrix. The fractal patterns of the points between the pixels are coupled, and the similarity of the adjacent fractal patterns is calculated through a sliding window to form a time coupling matrix; the fractal complexity calculation unit is used to set a neighborhood of a fixed size around each pixel for each time slice of the time coupling matrix, calculate the difference between the maximum and minimum values ​​of the indicator in the neighborhood, and obtain the local difference; the local differences of different scales are merged to generate a multi-scale morphological gradient map; the multi-scale morphological gradient map is segmented by time window and the fractal complexity is estimated; the toughness evaluation unit is used to mark the toughness value of each position by performing time comparison on the fractal complexity of the same position; then the connected domain analysis is performed on all positions to obtain multiple connected areas; for each connected area, the toughness value of the connected area is calculated based on its area and the toughness values ​​of all the positions therein.

Images