A multi-source data driven detection method and system for ecological resilience of a flood storage and detention area

Abstract

The present application belongs to the technical field of image processing, and particularly relates to a multi-source data driven detection method and system for ecological resilience of a flood storage and detention area, which comprises the following steps: step 1: performing geometric correction and radiation correction on an optical satellite remote sensing image to form a time series basic data set; step 2: a first axis and a second axis of an ecological state matrix correspond to spatial grids of the flood storage and detention area respectively, and a third axis corresponds to time; step 3: mapping the ecological state matrix to a fractal space to form a time coupling matrix; step 4: for each time slice of the time coupling matrix, a fixed size neighborhood is set around each pixel, and a multi-scale morphological gradient graph is subjected to time window segmentation and fractal complexity estimation; and step 5: according to the area and the resilience values of all positions therein, the resilience value of the connected region is calculated. The present application significantly enhances the accuracy and practicability of ecological resilience detection.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of image processing, and particularly relates to a multi-source data driven detection method and system for ecological resilience of a flood storage and detention area. BACKGROUND

[0002] Under the background of increasing dual needs of flood storage and ecological protection, the flood storage and detention area, as an important flood storage unit during the flood season, not only bears the functions of flood peak reduction and downstream disaster reduction, but also is increasingly endowed with the ecological function of restoring and maintaining the stability of the regional ecological system. Especially under the background of frequent extreme weather and enhanced fluctuation of hydrological processes in the basin, how to comprehensively evaluate the ecological resilience of the flood storage and detention area before and after being disturbed by floods has become an important direction of interdisciplinary research of water conservancy, ecology, remote sensing, geography and other disciplines. Ecological resilience is usually defined as the ability of an ecological system to restore its structure and function after being disturbed by a sudden disturbance, and its evaluation relies on the detailed description of the disturbance response characteristics, system recovery process and spatial differences. Under this background, remote sensing technology, as a kind of observation means with wide range of coverage, high timeliness and non-contact, has become an important data source for evaluating the changes of ecological system before and after flood events, especially the evolution of water body distribution and vegetation state.

[0003] The existing research on ecological change analysis of the flood storage and detention area mainly focuses on the following technical paths: first, the single-time or double-time change detection method based on optical remote sensing data, such as using vegetation and water related indexes such as normalized difference vegetation index (NDVI) and normalized difference water index (NDWI) to extract spatial distribution maps before and after the flood, and analyzing the change area through image difference, change threshold determination and other methods. This kind of method is simple to calculate, but it is insufficient to reflect the continuous change in time series, and it cannot effectively express the dynamic response and internal structure change of the system. Second, the method of introducing time series remote sensing analysis, such as using multi-temporal remote sensing data provided by MODIS, Sentinel-2 and other data, to analyze the disturbance and recovery process of the ecological system through continuous NDVI. But this kind of method is based on pixel level average or trend line fitting, and lacks the ability to express multi-scale structure change, making it difficult to capture the local severe response characteristics under the disturbance of flood. Third, the distributed model or physical process model is used to simulate the flood process, and the land use or vegetation type data is analyzed to analyze the ecological impact, but the model relies on a large number of hydrological parameters and field measured data, and has poor regional adaptability, making it difficult to be widely applied. Fourth, local research uses artificial field investigation and unmanned aerial low-altitude image to analyze details, which has high accuracy, but high cost and long cycle, and does not have the ability of large-scale dynamic monitoring. SUMMARY

[0004] In view of this, the main purpose of the present application is to provide a multi-source data driven detection method and system for the ecological resilience of a flood storage area, which extracts a water area mask and a vegetation index from optical remote sensing images before, during and after a flood to construct a three-dimensional ecological state matrix, maps it to a fractal space to enhance the spatiotemporal structure expression capability, and then calculates a time coupling matrix and a multi-scale morphological gradient map, estimates the fractal complexity in a sliding time window, and finally completes the pixel-level resilience value labeling and connected region resilience evaluation. This method can comprehensively depict the ecological response and recovery process under flood disturbance, improve the identification and quantification ability of multi-scale and time-series changes of the ecological system, and significantly enhance the accuracy and practicality of ecological resilience detection.

[0005] The technical solutions adopted by the present application are as follows:

[0006] A multi-source data driven detection method for the ecological resilience of a flood storage area, the method comprising:

[0007] Step 1: Collect optical satellite remote sensing images at different time periods before, during and after a flood in the flood storage area; perform geometric correction and radiation correction on the optical satellite remote sensing images to form a time series of basic data sets;

[0008] Step 2: Extract two types of indicators from the time series of basic data sets, namely water area mask and vegetation index; stack the obtained water area mask and vegetation index in time sequence 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 area, and the third axis corresponds to time;

[0009] 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; in the fractal space, couple the fractal patterns of adjacent time points, calculate the similarity of adjacent fractal patterns through a sliding window, and form a time coupling matrix;

[0010] Step 4: For each time slice of the time coupling matrix, set a fixed size neighborhood around each pixel, calculate the difference between the maximum and minimum values of the indicators in the neighborhood to obtain a local difference value; combine the local difference values of different scales to generate a multi-scale morphological gradient map; perform time window segmentation and fractal complexity estimation on the multi-scale morphological gradient map;

[0011] Step 5: Mark the resilience value of each location by comparing the fractal complexity of the same location over time; then perform connected component analysis on all locations to obtain multiple connected regions; for each connected region, calculate the resilience value of the connected region according to its area and the resilience values of all locations in the connected region.

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

[0013] Furthermore, in step 2, in the optical satellite remote sensing image of each time slice, based on the fact that water bodies have lower visible light reflectance and higher near-infrared absorption characteristics compared to surrounding land features, a fixed threshold range is set, and the locations where the pixel reflectance falls within this threshold range are identified as water bodies; morphological opening and closing operations are performed on the water body identification results to eliminate isolated noise and fill small-scale holes, thereby obtaining a binarized water body 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 that pixel.

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

[0015] Furthermore, in step 3, the process of mapping the ecological state matrix to fractal space includes: setting the iteration depth. The value ranges from 4 to 6; for each pixel vector in the ecological state matrix Execute the fractal mapping function as follows: ;in, , for Scale compression matrix It is a two-dimensional translation vector; The corresponding pixel grayscale values ​​are written to the same position in the fractal space, completing one layer of iteration; the above steps are repeated sequentially for all time slices until a complete multi-layer fractal pattern stack is generated; the scale compression matrix of each layer... Adopting a diagonal form ,in Translation vector Based on the row and column number of the cell Calculated as ;in, For line numbers, For column number; This is for the transpose operation.

[0016] Furthermore, within the fractal space, for adjacent time points... and When coupling fractal patterns, a size of [size missing] is used. A sliding window with a step size of 1 pixel. Integer subscript index; for each window position , The coordinate value of the first axis of the window. The second axis coordinate value of the window; calculate the window sub-blocks for two frames. and Structural similarity index : ,in and These represent the average gray levels of the two window sub-blocks. These represent the grayscale variances of the two window sub-blocks, Let covariance be the variance of the two window sub-blocks. , For all window positions Structural similarity index Perform pixel-level aggregation to obtain the coupling coefficients at adjacent time points. ,in The total number of valid cells covered by the window; all Fill the two-dimensional matrix with the data arranged in chronological order. This yields the time coupling matrix, where the row and column indices correspond to specific time points. Time coupling matrix Stored in symmetric form, with a value of 1 assigned to each diagonal position, its matrix elements satisfy... ,in ; Maximum time limit; Integer subscript index.

[0017] 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 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, and the local difference is calculated in the same way as the first scale. The local difference maps of the three scales are then linearly superimposed with weights of 0.5, 0.3 and 0.2 to obtain the fused multi-scale morphological gradient map.

[0018] Further, in step 4, for each time slice of the time-coupled matrix, a 3-pixel square neighborhood 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 a 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, and the local differences are calculated in the same way as the first scale. The local difference maps of the three scales are linearly superimposed with weights of 0.5, 0.3, and 0.2 to obtain a fused multi-scale morphological gradient map. The weighted multi-scale morphological gradient map is subjected to pixel-by-pixel maximum value selection, that is, the maximum value of the three-scale local differences is selected at the same pixel position to generate the final multi-scale morphological gradient map.

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

[0020] The application discloses a multi-source data driven detection system for ecological resilience of a flood storage and detention area, and relates to the technical field of ecological resilience detection. The system comprises a data acquisition unit configured to acquire optical satellite remote sensing images at different time periods before, during and after flood in the flood storage and detention area; the optical satellite remote sensing images are subjected to geometric correction and radiation correction to form a time series basic data set; an ecological state matrix construction unit is configured to extract two types of indexes 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 in time sequence are stacked to form a three-dimensional ecological state matrix; a first axis and a second axis of the ecological state matrix correspond to spatial grids of the flood storage and detention area, and a third axis corresponds to time; a time coupling matrix construction unit is configured to map the ecological state matrix to a fractal space, and utilize a hierarchical iterative fractal function to perform self-similarity transformation and scale compression on each layer of data; in the fractal space, fractal patterns at adjacent time points are coupled, similarity of adjacent fractal patterns is calculated through a sliding window, and a time coupling matrix is formed; a fractal complexity calculation unit is configured to, for each time slice of the time coupling matrix, set a neighborhood of a fixed size around each pixel, calculate a difference between a maximum value and a minimum value of indexes in the neighborhood, and obtain a local difference value; local difference values at different scales are combined to generate a multi-scale morphological gradient graph; the multi-scale morphological gradient graph is subjected to time window segmentation and fractal complexity estimation; a resilience evaluation unit is configured to mark resilience values of each position by comparing fractal complexities at the same position in time; all positions are subjected to connected domain analysis to obtain a plurality of connected regions; for each connected region, an area of the connected region and resilience values of all positions in the connected region are utilized to calculate a resilience value of the connected region.

[0021] By adopting the technical scheme, the following beneficial effects are achieved: the standardized processing and index extraction of multi-time-phase remote sensing images before, during and after floods can be realized, a three-dimensional ecological state matrix with a time sequence structure is formed, and the matrix is mapped to a fractal space through a hierarchical iterative fractal mapping function, so that the self-similarity expression of spatial texture and structure characteristics is effectively strengthened. In the fractal space, a structure similarity evaluation mechanism is introduced, a time coupling matrix is constructed, and the structure change degree between adjacent time slices is captured, so that the time sequence continuity of the disturbance and recovery process is revealed. In terms of spatial structure analysis, a multi-scale morphological gradient graph is generated through multi-scale neighborhood difference calculation, and sliding window segmentation and fractal complexity estimation are performed in the time dimension, so that scale fusion expression from local boundary response to overall complexity evolution is realized. Finally, by comparing the changes of the fractal complexity on the time axis, the pixel-level ecological resilience value marking is completed, and combined with the connected domain analysis, the resilience statistical results of the connected region are output, realizing the complete conversion from pixel-level dynamic monitoring to regional-level comprehensive evaluation. The present application breaks through the dependence on single-scale change, static difference and single-index response in the prior art, significantly improves the explainability and analysis accuracy of the ecological response process of the flood detention basin under flood disturbance, and provides a more scientific, stable and high-resolution support means for ecological restoration evaluation, post-disaster compensation calculation and ecological management of the flood detention basin. BRIEF DESCRIPTION OF DRAWINGS

[0022] Figure 1 A system structure schematic diagram of a multi-source data driven ecological resilience detection system of a flood detention basin is provided for the embodiments of the present application.

[0023] Figure 2 A fractal complexity estimation observation grid scale change schematic diagram is provided for the embodiments of the present application.

[0024] Figure 3 A connected domain analysis principle schematic diagram is provided for the embodiments of the present application. DETAILED DESCRIPTION

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

[0026] Any feature disclosed in this specification, unless stated otherwise, can be replaced by alternative features or equivalents having the same or a similar effect. That is, unless stated otherwise, each feature is one example only of a number of alternative or comparable features.

[0027] REFERENCE Figure 1 A multi-source data driven ecological resilience detection method of a flood detention basin, the method comprising:

[0028] Step 1: Collect optical satellite remote sensing images at different time periods before, during and after the flood in the 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;

[0029] Step 1 aims to build a time series basic data set that can be compared across periods for the multi-source data-driven ecological resilience detection method of the flood storage area. The principle is to collect optical satellite remote sensing images at different time periods before, during and after the flood in the flood storage area, and to perform geometric correction and radiation correction on the images, so that the same ground object has consistency in spatial positioning and comparability in radiation measurement, thereby ensuring the time consistency and spatial consistency when extracting the water area mask and vegetation index. In the specific implementation process, first, determine the target range according to the management boundary of the flood storage area, select the observation period covering the three stages of before, during and after the flood in combination with the flood process record, preferentially select optical satellite remote sensing images with visible and near-infrared bands, and try to ensure stable revisit period, similar imaging geometry and controllable seasonal differences to reduce the interference of non-target factors on the time series characteristics. To improve data availability, first, perform quality screening and cloud shadow mask processing on the candidate images, remove images with high cloud cover, obvious strip noise or sensor abnormalities, and only keep the image set that meets the requirements of complete spatial coverage, balanced time distribution and continuous flood key stage. Then, perform geometric correction, use orbital attitude information, sensor imaging geometry, terrain relief information and known ground control points to perform geographic referencing and distortion correction on the images, so that the optical satellite remote sensing images at different time periods are registered at the pixel level within the flood storage area, ensuring that the same ground object maintains consistent row and column positions in different time slices; when necessary, refine the registration and residual error evaluation of the local mismatch area until the spatial error meets the requirements of pixel-level consistency.

[0030] Then, the radiation correction is performed to restore the observation recorded by the sensor to the surface reflectivity with physical meaning, correct the atmospheric scattering and absorption, the difference of solar elevation and azimuth angle, and the non-uniformity of sensor response, and make the gray values between different time slices comparable. In order to further stabilize the time sequence characteristics, the image is also subjected to unified calibration and normalization processing to ensure that the visible and near-infrared bands remain consistent in radiation scale at different times. After completing the geometric correction and radiation correction, the boundaries of the flood storage and detention area are used to crop all optical satellite remote sensing images, only the target area data is retained, and the time axis is sorted. The slices before, during and after the flood are organized from early to late according to the acquisition time. For the time period with gaps, the qualified images of the adjacent date are preferentially supplemented to maintain the continuity of the time sequence. On the basis of the above processing, all optical satellite remote sensing images subjected to geometric correction and radiation correction are sorted and recorded according to the unified naming and metadata specification. The acquisition time, solar geometry, observation geometry and quality identification corresponding to each time slice are clear, and finally a clear structure, reliable source and cross-period comparable time sequence basic data set is formed. It lays the foundation for data consistency and traceability for subsequent extraction of water area mask and vegetation index from the time sequence basic data set and stacking of three-dimensional ecological state matrix in time sequence.

[0031] Step 2: Extract two types of indicators from the time sequence basic data set, namely: water area mask and vegetation index; stack the water area mask and vegetation index obtained 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 and detention area respectively, and the third axis corresponds to time;

[0032] The principle of step 2 is to extract the water area mask and the vegetation index in the time series basic data set at the same time, and stack the obtained water area mask and the vegetation index in time sequence to form a three-dimensional ecological state matrix, so that the first axis and the second axis correspond to the spatial grid of the flood storage and detention basin respectively, and the third axis corresponds to time, to support subsequent fractal space mapping, time coupling and multi-scale morphological gradient analysis. In the specific implementation process, first, the optical satellite remote sensing image quality of each time slice in the time series basic data set is reviewed to ensure that the geometric correction and the radiation correction are consistent and effective, and are unified to the UTM coordinate system with the WGS-84 ellipsoid as the reference, the pixel size is fixed at 10m x 10m, and the target range is obtained by overlapping and cutting with the boundary of the flood storage and detention basin. Then, in each time slice, according to the characteristics that the water body has lower visible light reflectivity and higher near-infrared absorption than the surrounding ground objects, a fixed threshold range is set to determine the pixels, and a preliminary binary result is obtained; in order to suppress isolated noise and small range holes, morphological open operation and closed operation are performed in turn, and only the water body patches with an area not less than 3 are retained by using the 8-neighbor connected rule, and at the same time, the retained patches are subjected to boundary refinement by morphological operations of expansion and corrosion, and a water area mask that is spatially coherent, boundary continuous and topologically stable is output.

[0033] The calculation of the vegetation index is based on the difference characteristics of the near-infrared band reflectivity and the red band reflectivity. First, the reflectivity values of each pixel in two bands are obtained, the difference value and the sum of the two are calculated, and then the ratio of the difference value to the sum is taken as the vegetation index of the pixel. Scale specification and invalid value identification are performed in the panoramic range to ensure one-to-one correspondence between the water area mask and the vegetation index in time and space. After completing the single time phase processing, the water area mask and the vegetation index of the same time slice are aligned at the pixel level by row number and column number and subjected to data consistency check to exclude inconsistent pixels caused by cloud shadow residues, sensor strips or local mismatches; the invalid data positions are uniformly coded to maintain the consistency of the time dimension statistical caliber in subsequent analysis. Subsequently, all time slices are organized from early to late according to the image acquisition time, the water area mask and the vegetation index of each time slice are 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 identifier of each time slice are recorded to ensure that the matrix can be traced back to the original slice in the time series basic data set. At this point, the time series basic data set is transformed into a three-dimensional ecological state matrix with both water body change and vegetation change information, which provides structured input and stable index source for mapping the ecological state matrix to a fractal space, coupling the fractal patterns of adjacent time points in the fractal space, and carrying out multi-scale morphological gradient and segmented fractal complexity estimation on the time coupling matrix.

[0034] Step 3: mapping the ecological state matrix to a fractal space, using a hierarchical iterative fractal function to perform self-similar transformation and scale compression on each layer of data; in the fractal space, coupling the fractal patterns of adjacent time points, calculating the similarity of adjacent fractal patterns through a sliding window, and forming a time coupling matrix;

[0035] The principle of step 3 is to map the three-dimensional ecological state matrix to a fractal space by time slicing without changing the topological relationship of the ecological state matrix space, to extract the cross-scale structural features formed by the water area mask and the vegetation index under the driving of the flood process of the flood storage and detention area by means of the hierarchical iterative fractal function to perform self-similar transformation and scale compression on each layer of data, and to couple the fractal patterns of adjacent time points in the fractal space, calculate the similarity of adjacent fractal patterns through a sliding window, and form a time coupling matrix. In the specific implementation process, first, read each time slice of the three-dimensional ecological state matrix from early to late according to time, take the time slice as an input layer, initialize the mapping result of the fractal space, set the iteration depth range to 4 to 6, and in each iteration, perform self-similar transformation and scale compression on the position of all pixels in the time slice according to the hierarchical iterative fractal function, so that the texture and boundary represented by the water area mask and the vegetation index in the same time slice repeatedly appear at a smaller scale, thereby gradually generating a fractal pattern stack with obvious self-similarity characteristics.

[0036] To ensure spatial consistency, the first axis and the second axis of the ecological state matrix are kept in one-to-one correspondence with the corresponding positions in the fractal space during the mapping process, so that the mapped fractal pattern can be strictly aligned with the original pixels in row number and column number. After completing each layer, the fractal pattern of the layer is written to the same position in the fractal space until the end of the entire iteration of the time slice. Then the above process is repeated for all time slices in sequence to form a complete multi-layer fractal pattern stack. After completing the mapping of the fractal space, two layers of fractal patterns of adjacent time points are selected as a pair of inputs in the fractal space. A sliding window with a size of 11x11 pixels is used to traverse the first axis and the second axis direction with a step size of 1 pixel. The similarity of the two layers of fractal patterns is calculated at each window position to measure the structural consistency and difference of adjacent time points in the local range. After the traversal is completed, the similarities obtained at all window positions are aggregated at the pixel level to obtain the coupling coefficient of the pair of adjacent time points. The coupling coefficient is filled into the corresponding position of a two-dimensional matrix in time order to gradually build a time coupling matrix. The time coupling matrix is organized in a symmetric form, and the diagonal line position is assigned a value of 1, so that it can not only reflect the coupling strength of the fractal patterns between adjacent time points, but also depict the phase difference and continuity before, during and after the flood through the overall structure of the matrix. Through this process of mapping the ecological state matrix to a fractal space, strengthening the self-similarity feature by using hierarchical iterative fractal function, and calculating the similarity of adjacent fractal patterns in a sliding window in the fractal space, the time coupling matrix can stably and objectively express the spatial structural changes of the flood storage and detention basin in different time slices, providing an input basis with cross-scale consistency for subsequent local difference calculation, multi-scale morphological gradient map generation, and time window segmentation and fractal complexity estimation in each time slice.

[0037] Step 4: For each time slice of the time coupling matrix, a fixed size neighborhood is set around each pixel, and the difference between the maximum and minimum values of the indicators in the neighborhood is calculated to obtain the local difference. Different scales of local differences are combined 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.

[0038] The principle of step 4 is based on the time coupling matrix to depict the structural mutation and evolution continuity of the flood storage and detention area in space and time. First, a fixed size neighborhood is established around each pixel in each time slice, and the local difference is obtained by calculating the difference between the maximum and minimum values of the indicators in the neighborhood. Then, the local difference of different scales is combined to generate a multi-scale morphological gradient map. Finally, the multi-scale morphological gradient map is segmented by time window and the fractal complexity is estimated. Thus, the spatial and temporal differences of the water area mask and the vegetation index change coupled in the fractal space triggered by the flood process are stably converted into comparable intensity and complexity measures. In the specific implementation process, in 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. The maximum and minimum values of the indicators in the neighborhood are extracted and the difference between the two values is calculated to obtain the local difference map of the first scale. To improve the ability to perceive the spatial heterogeneity of different scales, square neighborhoods with side lengths of 5 and 7 pixels are further used in the same time slice. The local difference maps of the second and third scales are obtained in the same way as the first scale. The local difference maps of the three scales are weighted and linearly superimposed according to weights of 0.5, 0.3 and 0.2 to output the fused multi-scale morphological gradient map. The maximum value of each pixel is selected to keep the result as another fusion method, i.e. the maximum value of the three-scale local difference 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 compare the methods and test the robustness under different research areas and different flood processes.

[0039] After the generation of the single time phase multiscale morphological gradient map, the multiscale morphological gradient map sequence is segmented along the time dimension using a sliding time window with a length of 5, and the window step is 1 time slice, so that there are 4 time slices overlapping between adjacent time windows to ensure time continuity and complete coverage of sudden changes. In each time window, the observation grid with a side length of 32 pixels, 16 pixels and 8 pixels is set in turn in the way of gradually reducing the observation grid, the number of non-zero grids occupied by the multiscale morphological gradient map is counted for each level of observation grid, and the size of each level of observation grid and the corresponding number of non-zero grids form a scatter set in the logarithmic coordinate system, and the least square linear fitting is used to obtain a value representing the fractal complexity of the time window, and the absolute value of the value is taken as the fractal complexity estimation result of the time window, which is recorded in the index position corresponding to the time window. Through the above process from local difference to multiscale morphological gradient map, and then to time window segmentation and fractal complexity estimation, the intensity index sensitive to the structural differences before, during and after the flood can be stably extracted on each time slice of the time coupling matrix, and the complexity change trajectory across scales can be extracted in the sliding time window, which provides complete and consistent input for subsequent marking of the resilience value by time comparison of the fractal complexity at the same position, and calculation of the resilience value of the connected region in the connected domain analysis.

[0040] Step 5: Mark the resilience value of each position by time comparison of the fractal complexity at the same position; then perform connected domain analysis on all positions to obtain a plurality of connected regions; for each connected region, calculate the resilience value of the connected region according to the area and the resilience value of all positions in the connected region.

[0041] The principle of step 5 is to use the fractal complexity of the multi-scale morphological gradient image sequence as the stable characterization in the time dimension, and to compare the fractal complexity changes in the same position in the continuous time window before, during and after the flood, to obtain the resilience value that can reflect the disturbance response and recovery process, and to aggregate adjacent positions into connected regions in space through connected component analysis, and to comprehensively evaluate the connected regions by combining the area and the resilience value of all positions in the connected regions, 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 position in the time sequence of fractal complexity, and the time window is paired and compared from early to late, and the change amplitude of the flood relative to the pre-flood and the recovery amplitude of the post-flood relative to 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 position; in order to reduce the influence of incidental noise, first, time series smoothing and outlier detection are performed on the time sequence of fractal complexity, and the breakpoints caused by cloud shadows and invalid data are filled and quality marked, to ensure that the time comparison of the same position is based on available slices; then, according to the key period of the flood process, the disturbance peak and the recovery inflection point are extracted on the time axis of the same position, and the resilience characteristics are described by the fall amplitude and recovery speed in a number of windows after the disturbance, and the resilience value at the pixel level is standardized and mapped, so that different positions can be compared at the same scale.

[0042] After completing the pixel-level resilience value marking, connected component analysis is performed on all positions in the entire flood storage and detention area, and adjacent positions with valid resilience values are aggregated according to the 8-neighbor connectivity rule to obtain multiple connected regions; during the aggregation process, small areas, fragmented shapes or isolated small blocks caused by noise are subjected to minimum area constraint and morphological correction to ensure the spatial continuity and boundary reasonableness of the connected regions; for each connected region, first, the area of the connected region is calculated, and then the resilience values of all positions in the region are statistically summarized, and the area and the pixel-level resilience value are used to determine the resilience value of the connected region, so that the connected regions with larger area, higher internal resilience value and more uniform spatial distribution obtain higher stability expression in comprehensive evaluation; at the same time, the time index range and quality mark of the connected region are retained for horizontal comparison between different flood processes. Finally, two types of results including the pixel-level resilience value raster and the connected region resilience value list are output, the former directly reflects the time comparison results of the same position, and the latter reflects the regional aggregation characteristics based on connected component analysis; both of them are consistent with the row number, column number and time index, and can be traced back to the time sequence basic data set and the optical satellite remote sensing image, to realize the closed-loop evaluation of the multi-source data-driven ecological resilience detection method of the flood storage and detention area from before, during and after the flood.

[0043] Further, multispectral images with visible and near-infrared bands are selected as optical satellite remote sensing images at the data source level, and are acquired at the same revisit period in three stages. On the one hand, the sampling interval on the time axis remains constant, reducing the seasonal and sunlight angle differences caused by observation time differences. On the other hand, since the optical satellite remote sensing images cover both visible and near-infrared bands, the water body and vegetation conditions can be simultaneously described at the same time, forming a dual-index system sensitive to the processes of flood inundation, water body expansion, and vegetation recovery after water recession. In terms of spatial quality, products with a spatial resolution of not more than ten meters are preferentially selected. Through finer pixel scales, the identification ability of small-scale water surfaces, narrow ditches, tidal flat edges, and shoreline transition zones is improved, and the expression accuracy of vegetation index in the case of mixed ground object pixels is also improved. To ensure the spatial consistency of cross-period comparison, all optical satellite remote sensing images are reprojected after geometric correction and radiometric correction, and are unified to the UTM coordinate system with WGS-84 ellipsoid as the reference, with a fixed pixel size of 10m x 10m. Through the unification of the coordinate system and the pixel scale, the slices under different imaging dates, different orbits, and different attitude conditions are one-to-one corresponding in terms of row number and column number within the flood storage and detention area, providing a stable spatial framework for the subsequent three-dimensional ecological state matrix construction.

[0044] In the index extraction link, for each time slice, according to the lower visible light reflectivity and higher near-infrared absorption characteristics of water body compared with surrounding ground objects, a fixed threshold range is set to determine the pixels, and the water body candidate area is obtained initially; considering that cloud shadow, sensor strip and terrain shadow may cause isolated noise and holes, morphological opening operation is first performed on the binary result to eliminate discrete small spots, and then morphological closing operation is performed to fill small holes, so that the binary water area mask which is continuous in space and has smoother boundary is obtained. In order to further suppress false water body noise and ensure that the water body patch has regional significance, the 8-neighbor connectedness rule is used to analyze the connectedness of the water area mask, and only the water body patch with an area of not less than 3 is retained; for the isolated patches still existing, morphological operations of inflation and corrosion are performed on them combined with scene characteristics, so as to optimize the patch morphology and topological structure without damaging the true boundary. Synchronously with the water area mask, the near-infrared band reflectivity and red band reflectivity are obtained at each pixel position, the difference value of the two is calculated first, and then the sum of the two is calculated, finally the ratio of the difference value to the sum is taken as the vegetation index of the pixel, and consistent calculation process and invalid value identification strategy are maintained in the whole time phase, so as to ensure the comparability of the vegetation index in time dimension and the splicing of the vegetation index in space dimension. After the water area mask and the vegetation index of a single time phase are extracted, they are aligned at the pixel level on the grid level and the quality is reviewed, individual inconsistent pixels caused by cloud shadow residual or geometric micro-mismatch are removed, and the row number and column number indexes determined in the re-projection stage of the optical satellite remote sensing image are inherited, so that each pixel has traceable spatial identification on different time slices. Then, according to the image acquisition time from early to late, the water area mask and the vegetation index of all time slices are stacked in time sequence to form a three-dimensional ecological state matrix, wherein the first axis and the second axis correspond to the row number and the column number of the UTM grid respectively, and the third axis corresponds to the time, so that the spatio-temporal change information about water body expansion, shrinkage and vegetation damage and recovery is embedded into a structured data container.

[0045] The container is directly mapped as input in subsequent processes to a fractal space, using hierarchical iterative fractal functions to perform self-similar transformation and scale compression at each layer to enhance the texture boundary and block pattern shaped by the water mask and vegetation index, and to calculate the similarity of fractal patterns at adjacent time points in the fractal space through a sliding window to form a time coupling matrix; each time slice of the time coupling matrix calculates the difference between the maximum and minimum values of the index in a fixed size neighborhood to obtain the local difference, and merges at multiple scales to generate a multi-scale morphological gradient map; the multi-scale morphological gradient map is segmented along the time dimension using a sliding time window, and the number of non-zero grids is counted based on the observation grid, thereby completing the fractal complexity estimation. Through the above bottom-up and layer-by-layer dependent process, on the one hand, it ensures the unity of optical satellite remote sensing images in revisit period, spatial resolution, coordinate reference and pixel scale, and on the other hand, it ensures the consistency of water mask and vegetation index in decision criteria, morphological processing, connectivity constraints and pixel-level alignment. Finally, the three-dimensional ecological state matrix serves as a bridge to stably convert the original observation of multispectral images into cross-period comparable structural representations, and to provide reliable data basis and clear index system for subsequent time comparison of fractal complexity at the same location to mark the resilience value, and connectivity domain analysis and calculation of the resilience value of all locations. The whole implementation emphasizes the consistency management of data and the traceability of the process, and the results of any time slice can be traced back to the original optical satellite remote sensing image and the corresponding re-projection, cropping and quality identification record, thereby meeting the stability and reusability requirements of the multi-source data-driven flood detention basin ecological resilience detection method in engineering deployment and cross-regional migration.

[0046] Further, in the specific implementation process, firstly, each time slice is read in time sequence from the three-dimensional ecological state matrix, the write cache of the corresponding layer of the fractal space is initialized in each processing, and the value range of the iteration depth N is set to 4 to 6, so as to ensure multi-level description of local and global structure in limited calculation amount; subsequently, the layered iterative fractal function is applied to each pixel vector in the time slice one by one, a 2*2 scale compression matrix is used to scale the pixel position in each direction in the iteration process, the scaling coefficient is decreased by the power of 2, and the high level is quickly converged to the local structure; at the same time, a two-dimensional translation vector is applied to the position, the translation amount is shifted by the power of 2 as the denominator of the row number and the column number of the pixel respectively, so that the original space pattern is smoothly expanded while being scaled, and it is ensured that the self-similarity transformation does not destroy the traceability of the row number and the column number in the fractal space, and does not produce topological dislocation across the grid. When each layer iteration is completed, the fractal pattern obtained at the current layer is written into the same position of the fractal space according to the original row number and column number, the strict alignment of the first axis and the second axis is ensured, and the correspondence between the layer and the original time slice is recorded; the processing of all levels is completed for the same time slice according to the set iteration depth, and the multi-level fractal pattern of the time slice is formed; the above operation is repeated in sequence for all time slices, and finally a complete multi-level fractal pattern stack is generated.

[0047] In order to avoid boundary effect and invalid value propagation, a uniform quality identifier is used for invalid pixels in the three-dimensional ecological state matrix in the mapping process, and a mask processing is performed during iterative writing, and a mirror or constant continuation strategy is enabled for pixels near the boundary, so as to maintain the continuity of the fractal space on the first axis and the second axis; in order to reduce the influence of numerical instability on subsequent similarity calculation, the gray scale can be uniformly scaled after each layer is completed, so that the gray scale dimension of different time slices and different iteration layers remains consistent. After the fractal space mapping is completed, the process of coupling the fractal patterns of adjacent time points in the fractal space is entered, the specific method is to take two adjacent fractal patterns on the time axis as a pair of input, use a sliding window with a size of 11*11 pixels and take 1 pixel as a step, traverse the whole image in the first axis and the second axis direction, calculate the structural similarity index at each window position taking two frames of window sub-blocks as objects, the structural similarity index is described by the average gray scale, gray scale variance and covariance of the two window sub-blocks, so as to consider the matching degree of brightness, consistency and structure at the same time; during calculation, invalid pixels are strictly ignored, and in the case that the number of valid pixels in the window is insufficient, the result of this position is excluded according to the quality identifier, so as to avoid the deviation of statistical quantity caused by sparse data; after the traversal is completed, the structural similarity indexes obtained at all window positions are pixel-level summarized, the coupling coefficient of the adjacent time points is calculated, and the numerical range is limited to 0 to 1, so as to intuitively represent the continuous degree from complete dissimilarity to complete identity.

[0048] The coupling coefficients of each adjacent time point are filled in the corresponding entries of the two-dimensional matrix from early to late, the time coupling matrix is gradually constructed, and is stored in a symmetric form, and the diagonal position is strictly set to 1 to reflect the complete consistency of the same time point with itself; the row index and column index of the time coupling matrix correspond to the time points of t1 to tT respectively, the matrix elements meet the constraints of interchange symmetry and value range, which is convenient for reading directly according to the time slice in the subsequent multi-scale morphological gradient calculation. In order to ensure the comparability of the time coupling matrix between different flood processes and different study areas, a fixed window size and step strategy is adopted when summarizing the coupling coefficients, and the settings of 11*11 and 1 are kept unchanged; at the same time, the quality indicators such as the total number of valid pixels M and the proportion of invalid pixels of each pair of time slices are reserved, so as to down-regulate the weight or eliminate the low-quality entries in the subsequent evaluation. On the implementation level, the coupling calculation can use block parallel and vectorization operation to reduce the calculation overhead, and for large-scale detention basin data, it can be processed by row number and column number and spliced seamlessly in the final stage, ensuring that each entry in the time 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. The time coupling matrix maintains spatial alignment and time sequence while condensing the cross-scale differences in texture, boundary and block pattern before, during and after the flood, providing strict consistency input for establishing a fixed size neighborhood on each time slice, calculating the local difference and synthesizing the multi-scale morphological gradient map; at the same time, the symmetric structure and standardized value range of the time 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 comparing the fractal complexity of the same position in time, performing connected component analysis and calculating the connected region resilience value, realizing the realizability and reusability of the multi-source data driven detention basin ecological resilience detection method on the whole link.

[0049] Further, by establishing a fixed-size neighborhood for each pixel in the matrix within each time slice and calculating the difference between the maximum and minimum values of the indicators in the neighborhood, a local difference value sensitive to abrupt boundaries and local amplitudes is obtained. Then, the process is repeated in a multi-scale manner and weighted linear superposition or pixel-by-pixel maximum selection is performed to obtain a multi-scale morphological gradient map that balances details and the whole. Subsequently, a sliding time window with a length of 5 is used along the time dimension for segmentation, and the window step is 1 time slice, so that there are 4 time slices of overlap between adjacent time windows. Finally, in each time window, the observation grid is gradually reduced, and the number of non-zero grids occupied by the gradient map is counted. The absolute value of the slope of the least squares linear fitting is used as the fractal complexity value, so that the structural changes driven by floods are converted into comparable complexity trajectories and provide a robust input for subsequent resilience value labeling and connected component analysis. In the specific implementation process, first, on each time slice of the time-coupled matrix, a square neighborhood with a side length of 3 pixels is established for each pixel in the matrix.

[0050] All indicator values in the neighborhood are traversed, and the difference between the maximum and minimum values is calculated to generate a local difference value map at the first scale, which preferentially captures fine-grained boundaries and narrow texture changes. To improve the response capability to low-frequency structures corresponding to large-scale water body expansion, contraction, and slow vegetation index recovery, square neighborhoods with side lengths of 5 and 7 pixels are further established on the same time slice. Local difference values are calculated in the same way as the first scale to obtain local difference value maps at the second and third scales. After completing the local difference value calculation at the three scales, to balance boundary sensitivity and regional stability in a unified expression, the local difference value maps at the three scales are first weighted and linearly superimposed according to weights of 0.5, 0.3, and 0.2 to obtain a fused multi-scale morphological gradient map. The fusion result is checked for invalid values and boundary completeness to ensure that the processing near the image edge and in the missing location does not introduce false gradients. At the same time, the results of pixel-by-pixel maximum selection are calculated in parallel on the same data, i.e., the maximum value of the local difference values at the three scales 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 complex feature transition zones or stages with strong sudden disturbances. The two fusion methods are completely consistent in data organization, time indexing, and quality identification. The subsequent process can select one or compare them in parallel to test robustness according to the research objectives.

[0051] After the generation of multi-scale morphological gradient maps of single time phase, they are composed into a sequence in the time axis according to the acquisition order, and segmented along the time dimension using a sliding time window with a length of 5. The window step is 1 time slice, so there are 4 overlapping time slices between adjacent time windows. This setting ensures time continuity while improving the coverage probability of rapid fluctuations and stage platforms. For each time window, to avoid the disturbance of fine texture driven by noise on the pixel level, a strategy of gradually reducing the observation grid is used to extract cross-scale occupying features. Specifically, observation grids with side lengths of 32 pixels, 16 pixels, and 8 pixels are set in turn, and the multi-scale morphological gradient maps in the time window are covered on each level of 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 proportion of invalid pixels are saved for each level to remove abnormal windows later. Then the observation grid size and the corresponding non-zero grid number of each level are formed into a scatter set in the log coordinate system, and the slope of the least squares straight line is obtained. 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 fixed and unchanged, and the window length and the window step are also fixed and unchanged. Windows with insufficient non-zero grid counts or excessive invalid pixels are down-weighted or marked as unusable.

[0052] The entire implementation takes the time-coupled matrix as input, first generates multi-scale morphological gradient maps at each time slice, and then robustly extracts fractal complexity with a sliding time window, so as to express the instantaneous intensity of the difference between the maximum and minimum values in the local difference of the spatial neighborhood, and condense it into a complexity measure that can be compared across scales and stages. In the context of flood detention and retention areas, a 3-pixel square neighborhood is particularly sensitive to narrow waterways, shoreline micro-migration, and local damage to vegetation, while 5-pixel and 7-pixel square neighborhoods are more sensitive to sheet flooding, lake shore shallow water recession, and community-level recovery. Weighted linear superposition of weights 0.5, 0.3 and 0.2 enhances the priority attention to detail boundaries, while retaining the contribution of medium and large scale morphology, and the pixel-by-pixel maximum selection provides stronger disturbance indication in extreme events or image contrast sudden rise stages. Through the above consistent, traceable and parameter-fixed process, the local difference of the time-coupled matrix in space is systematically integrated into multi-scale morphological gradient maps, and in time, the fractal complexity sequence is stably output through a sliding time window with a length of 5 and a step of 1, providing a continuous, comparable and noise-resistant spatio-temporal feature basis for subsequent resilience value marking by comparing the fractal complexity at the same location in time, and calculating the resilience value of the connected region in the connected component analysis combined with the area and pixel-level resilience value, ensuring that the multi-source data-driven ecological resilience detection method of flood detention and retention areas has consistent measurement standards and engineering implementability at all levels from single pixel to connected region.

[0053] Further, in step 5, for the fractal complexity sequence formed in time sequence at the same pixel position, firstly, the average value of the fractal complexity of the first 3 time windows in the sequence is taken as the reference 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 reference complexity value plus 10 percent, the window number is recorded as the recovery time; if the recovery condition is still not met in the next 8 time windows, the pixel is marked as a low recovery position. In step 5, according to the length of the recovery time, each pixel is given a discrete resilience level, the pixel with a recovery time not more than 5 time windows is given a resilience level 5, the pixel with a recovery time of 6 to 8 time windows is given a resilience level 4, the pixel with a recovery time of 9 to 12 time windows is given a resilience level 3, the pixel with a recovery time of 13 to 16 time windows is given a resilience level 2, and the pixel with a recovery time more than 16 time windows is given a resilience level 1. For all the pixels marked with resilience levels, an eight-connected region marking algorithm is performed to identify the connected domains; the connected domains with an area less than 9 pixels are directly removed and do not participate in subsequent resilience calculation; for the remaining connected domains, the number of area pixels is recorded. For each remaining connected domain, firstly, the arithmetic mean of the resilience levels of the pixels in the domain is calculated, and then the average value is multiplied by the area ratio of the total number of pixels in the domain to the total number of grids to obtain a preliminary regional resilience value; if the regional resilience value is higher than 4 and the area exceeds 5 percent of the total grid area, the regional resilience value is increased by 0.5 to highlight the importance of large-area high-resilience regions. If any connected domain contains a number of pixels marked by the water area mask that accounts for more than 30 percent of the area of the region, the resilience value of the connected domain is reduced by 1; if the proportion is between 15 percent and 30 percent, the resilience value is reduced by 0.5; if the proportion is less than 15 percent, the resilience value is not adjusted. The resilience values of all connected domains are arranged in descending order, and labels "high resilience area", "medium-high resilience area", "medium resilience area", "medium-low resilience area" and "low resilience area" are assigned in turn according to the set rules; at the same time, a regional-level resilience distribution map consistent with the grid is output, each pixel is filled with the final resilience value of the connected domain to which it belongs, which is used for subsequent ecological restoration decision-making of flood storage and detention areas.

[0054] For ease of reproduction, a rectangular test area in a flood storage and detention area is selected, and the UTM coordinate system with WGS-84 ellipsoid as the reference is uniformly adopted, and the pixel size is fixed at 10m x 10m. The row number and column number are used as spatial indexes to construct a three-dimensional ecological state matrix and subsequent mapping and statistics. The grid size of the test area is set as (row number , column number ). A total of time points of optical satellite remote sensing images are collected: before flood , during flood , and after flood The images are multispectral images with visible and near-infrared bands, and the spatial resolution is 10 m, and the revisit period is the same. After completing the geometric correction and radiation correction of all time phases, reprojecting to the UTM coordinate system and cutting to the test area range.

[0055] For each time slice , a fixed threshold is set according to the characteristics of low reflection of water in visible light and high absorption of near-infrared, to obtain the preliminary water body determination, and then morphological open operation and close operation are performed, and 8-neighbor connected rule is adopted to only retain water body patches with area not less than 3 pixels, and binary water mask is output . At the same time, the vegetation index (the ratio of the difference and the sum of the reflectivity of the red band and the near-infrared band) is calculated: . Wherein, is the reflectivity of the near-infrared band, is the reflectivity of the red band, is the row number and column number of the pixel, is the time index. Example: on the pixel , if , then . Stack and in time order to form a three-dimensional ecological state matrix, whose first axis is the row number, the second axis is the column number, and the third axis is the time.

[0056] Set the iteration depth (the specific value in the range of 4 to 6). For each time slice, map the two-dimensional field to a fractal space according to the hierarchical iterative fractal function. For any pixel position vector (composed of row number and column number), define iteration: , wherein, is the scale compression matrix, is the isotropic compression coefficient of the layer, is the two-dimensional translation vector, are the pixel row and column numbers, indicates the transpose operation. Write the pixel gray value of the corresponding position of after iteration to the same position in the fractal space, complete the writing of this layer, and repeat for all time slices to obtain a stack of multi-layer fractal patterns. In the fractal space, traverse adjacent time points using a sliding window with a size of and a step of 1 to calculate the structural similarity index (SSIM form): , wherein, is the coordinate of the upper left corner of the window in the first axis and the second axis, are the average gray values of the two time sub-blocks, is the variance, Covariance, Covariance, Coupling coefficient at adjacent time points is obtained by averaging at pixel level: where is the total number of effective windows. Fill in the two-dimensional matrix in time order to form the time coupling matrix, satisfying , , and . Example: in a certain small area, take a window at , if , then the single-window coupling coefficient is . Average over the total number of effective windows to obtain .

[0057] For each time slice of the time coupling matrix, establish a 3-pixel square neighborhood for each pixel in the matrix, and calculate the difference between the maximum and minimum values in the neighborhood to obtain the first-scale local difference map . Construct neighborhoods with edge lengths of 5 and 7 pixels in the same way to obtain and . Weighted linear superposition with weights of 0.5, 0.3, and 0.2 to obtain the fused multi-scale morphological gradient map: . Parallelly retain the pixel-by-pixel maximum value selection: . Example: at pixel , if the three-scale local difference is , then .

[0058] Segment (or ) along the time dimension using a sliding time window with a length of 5 and a step of 1. For each time window, perform fractal complexity estimation (box-counting idea): use observation grids with edge lengths of 32 pixels, 16 pixels, and 8 pixels to cover the entire image in turn, and count the number of non-zero grids occupied by the gradient map, denoted as , where is the edge length of the -level observation grid, and is the corresponding number of non-zero grids. Perform least squares linear fitting on , and take the absolute value of the slope as the fractal complexity of the time window. Example: for a window containing to , if the statistics show that , the least squares regression gives a slope of , then . Write according to the center time index of the window.

[0059] The fractal complexity of the same location is compared over time to define the pixel-level resilience value. Let a pixel The mean fractal complexity of the pre-flood window is The maximum fractal complexity occurring during the flood is The mean fractal complexity of the post-flood stable period is The resilience value is defined as where represents the degree of recovery relative to the pre-flood baseline, with 1 indicating full recovery to baseline levels and 0 indicating an unrecovered state equivalent to the peak disturbance. Example: for pixel , if , (occurring near ), then .

[0060] The full image pixel-level resilience value raster is obtained After which, connected component analysis is performed. Using the 8-adjacency connectedness rule, the valid pixels are aggregated into a number of connected regions . For each connected region, the area is calculated, where is the area of a single pixel (10m x 10m). The area-weighted mean of the resilience values within the region is calculated as the resilience value of the connected region: where is the number of pixels in the region . Example: if a region contains pixels, with an average pixel-level resilience value of , then the area of the region is and the resilience value of the region is . In the same way, the list of all regions is obtained, which is the resilience assessment result for the management unit.

[0061] 1) In the fractal space mapping stage, the selection of the iteration depth affects the trade-off between detail and overall structure, which can balance the computational load and expressive power under the data size of this example. 2) In the calculation of the structural similarity index, the constant is used to avoid instability caused by the denominator approaching 0; when the proportion of valid pixels in the window is lower than the set threshold (for example, ), the can not be counted. 3) In the multi-scale morphological gradient map fusion, the weight emphasizes the identification of fine-scale mutations; if more attention is paid to the extreme disturbance stage, the use of Enhance sensitivity. 4) When estimating fractal complexity, observe the side length of the observation grid. The time window length (5) and step size (1) are kept fixed to ensure comparability across periods; when the proportion of invalid pixels in a window exceeds the threshold, the window is closed. Marked as invalid and ignored in toughness value calculation. 5) In the definition of cell-level toughness value, and The mean of the windows before and after the flood can be taken separately; if and If the values ​​are too close, the denominator will approach 0. A minimum difference constraint should be applied, or a percentile baseline should be used to maintain numerical stability. 6) Connected component analysis output includes both pixel-level resilience values ​​and raster values. It also includes a list of connected regions. Both can be traced back to the row numbers, column numbers, and time indexes of the time coupling matrix, multi-scale morphological gradient map, and three-dimensional ecological state matrix, meeting the needs of engineering auditing.

[0062] Through the above examples, a complete process for detecting the ecological resilience of flood storage and detention areas driven by multi-source data has been realized. Starting with the standardized processing of optical satellite remote sensing imagery, index extraction, and the construction of a three-dimensional ecological state matrix, the temporal structural relationships are characterized through fractal spatial mapping and temporal coupling matrices. Multi-scale morphological gradients and fractal complexity are used to characterize the intensity and complexity of spatiotemporal changes. Finally, quantitative assessments for management units are completed using pixel-level resilience values ​​and connected region resilience values. The above parameters and formulas can be slightly adjusted according to regional characteristics and data quality, but the indexing system, calculation methods, and quality control strategies remain consistent, allowing for reuse across different flood storage and detention areas.

[0063] refer to Figure 2, the technical principle of the present application for estimating the observation grid scale change of fractal complexity is shown in the figure. The technical principle is the key link for further implementing the fractal complexity estimation after the time window segmentation of the multi-scale morphological gradient map in step 4. Specifically, for the multi-scale morphological gradient map in each time window, the present application adopts the technical scheme of gradually reducing the observation grid to obtain the fractal feature parameters. The process includes three different observation grid scale levels: the first level is a thirty-two-pixel edge length observation grid, the second level is a sixteen-pixel edge length observation grid, and the third level is an eight-pixel edge length observation grid. In the first level thirty-two-pixel observation grid, the entire multi-scale morphological gradient map is divided according to the thirty-two-pixel edge length square grid, and the number of non-zero grids occupied by the gradient map is counted. The left part of the figure shows the spatial distribution of the grid at this level, where the black filled area represents the non-zero grid position occupied by the multi-scale morphological gradient map. Through the grid-by-grid traversal statistics, the number of non-zero grids at this scale is obtained as 8. In the second level sixteen-pixel observation grid, the same technical scheme is adopted, but the grid edge length is reduced to sixteen pixels. Due to the reduction of grid size, more detailed spatial feature changes can be captured, and accordingly the number of non-zero grids occupied by the gradient map increases to 18. The middle part of the figure shows the detailed structure of the grid at this level, and it can be observed that the grid density is significantly improved and the spatial details of the non-zero grid distribution are observed. In the third level eight-pixel observation grid, the observation grid edge length is further reduced to eight pixels, and the grid resolution reaches the highest level. Due to the further reduction of grid size, more local morphological change features can be identified, and therefore the number of non-zero grids occupied by the gradient map further increases to 32. The right part of the figure shows the distribution characteristics of the grid at this level, and the grid density reaches the maximum value. After obtaining the three different levels of observation grid size and the corresponding number of non-zero grids, the present application draws these data points in the logarithmic coordinate system to form a scatter plot. The logarithmic coordinate system at the bottom of the figure shows the technical implementation process: the horizontal axis represents the logarithmic value of the observation grid edge length, and the vertical axis represents the logarithmic value of the number of non-zero grids. The three data points correspond to the observation results of (8, 32), (16, 18), and (32, 8) respectively.

[0064] Figure 3 The figure is a technical principle diagram for connected component analysis, which shows the core technical implementation process of the connected component analysis in the present application. Figure 3 The left side shows that the resilience value position distribution part shows the spatial organization structure of the pixel grid in the flood storage and detention area. The grid uses a unified UTM coordinate system, and the pixel size is fixed at 10m x 10m. The black dots in the figure represent the marked resilience value positions after time comparison, which are determined by time comparison of the fractal complexity at the same position. The resilience value marking points show non-uniform distribution in space, reflecting the spatial heterogeneity of ecological resilience in the flood storage and detention area.Figure 3 The middle part details the technical implementation of the 8-adjacent connectivity rule. This rule is based on the center pixel and establishes connectivity with its 8 adjacent pixels, including horizontally, vertically and diagonally adjacent pixels. The 8-adjacent connectivity rule ensures the integrity of the connectivity analysis and can accurately identify spatially adjacent resilience value positions. The center pixel is filled with black, and the 8 adjacent pixels are filled with white and labeled with adjacent numbers, clearly showing the spatial configuration of the adjacent relationship. Figure 3 The right side shows the specific morphology of the connected region result. Connected region A contains 3 adjacent resilience value positions, showing linear distribution; connected region B contains 7 resilience value positions, forming a larger connected region; and connected region C contains 2 resilience value positions, constituting the smallest connected unit. Each connected region has independent spatial boundaries and internal structural characteristics.

[0065] Although the specific embodiments of the present application are described above, those skilled in the art should understand that these specific embodiments are only illustrative, and those skilled in the art can make various omissions, substitutions and changes to the details of the above method and system without departing from the principles and essence of the present application. For example, combining the above method steps, so as to perform substantially the same function in substantially the same way to achieve substantially the same result, is within the scope of the present application. Therefore, the scope of the present application is only limited by the appended claims.

Claims

1. A multi-source data-driven detection method for ecological resilience of a flood detention basin, characterized in that, The method comprises: Step 1: collecting optical satellite remote sensing images at different time periods before, during and after flood in the flood storage area; performing geometric correction and radiation correction on the optical satellite remote sensing images to form a time series basic data set; Step 2: extracting two types of indexes from the time series basic data set, which are water area mask and vegetation index; stacking the obtained water area mask and vegetation index 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 of the ecological state matrix corresponds to time; Step 3: mapping the ecological state matrix to a fractal space, using hierarchical iterative fractal function to perform self-similar transformation and scale compression on each layer of data; in the fractal space, coupling the fractal patterns of adjacent time points, calculating the similarity of adjacent fractal patterns through sliding window to form a time coupling matrix; Step 4: for each time slice of the time coupling matrix, setting a fixed size neighborhood around each pixel, calculating the difference between the maximum and minimum values of the indexes in the neighborhood to obtain the local difference value; combining local difference values of different scales to generate a multi-scale morphological gradient map; performing time window segmentation and fractal complexity estimation on the multi-scale morphological gradient map; Step 5: by comparing the fractal complexity of the same position in time, marking the resilience value of each position; then performing connected component analysis on all positions to obtain multiple connected regions; for each connected region, calculating the resilience value of the connected region according to its area and the resilience values of all positions in the connected region.

2. The multi-source data-driven detection method of the ecological resilience of a flood detention basin according to claim 1, wherein, The optical satellite remote sensing image is a multispectral image with visible light and near-infrared bands, and is acquired at the same revisit period before, during and after flood at different time periods; the spatial resolution of the optical satellite remote sensing image is not greater than ten meters.

3. The multi-source data-driven detection method of the ecological resilience of a flood detention basin according to claim 2, wherein, In step 2, in the optical satellite remote sensing image of each time slice, according to the lower visible light reflectivity and higher near-infrared absorption characteristics of water body compared with surrounding ground objects, a fixed threshold range is set, and the positions of the pixel reflectivity falling within the threshold range are determined as water area; And performing morphological opening operation and closing operation on the water area determination result to eliminate isolated noise and fill small range holes, thereby obtaining a binary water area mask; for each pixel, obtaining its near-infrared band reflectivity and red band reflectivity, calculating the difference value, and then calculating the sum, taking the ratio of the difference value to the sum as the vegetation index of the pixel.

4. The multi-source data-driven detection method of the ecological resilience of a flood detention basin according to claim 3, wherein, In step 2, the water area mask adopts 8-neighbor connected rule, only retains water body patches with an area not less than 3, and performs morphological operations of dilation and erosion on isolated patches to reduce false water body noise; all optical satellite remote sensing images are reprojected and unified to UTM coordinate system with WGS-84 ellipsoid as the reference, and the pixel size is fixed at 10m x 10m; the first axis and the second axis of the ecological state matrix correspond to the row number and column number of the UTM grid respectively, and the third axis is arranged in the order from early to late according to the image acquisition time.

5. The multi-source data-driven detection method of Claim 4, wherein, Step 3, the process of mapping the ecological state matrix to fractal space, includes: setting the iteration depth. The value ranges from 4 to 6; for each pixel vector in the ecological state matrix Execute the fractal mapping function as follows: ;in, , for Scale compression matrix It is a two-dimensional translation vector; The corresponding pixel grayscale values ​​are written to the same position in the fractal space, completing one layer of iteration; the above steps are repeated sequentially for all time slices until a complete multi-layer fractal pattern stack is generated; the scale compression matrix of each layer... Adopting a diagonal form ,in Translation vector Based on the row and column number of the cell Calculated as ;in, For line numbers, For column number; This is for the transpose operation.

6. The multi-source data-driven detection method of the ecological resilience of a flood detention basin according to claim 5, wherein, In fractal space, for adjacent time points and When coupling fractal patterns, a size of [size missing] is used. A sliding window with a step size of 1 pixel. Integer subscript index; for each window position , The coordinate value of the first axis of the window. The second axis coordinate value of the window; calculate the window sub-blocks for two frames. and Structural similarity index : ,in and These represent the average gray levels of the two window sub-blocks. These represent the grayscale variances of the two window sub-blocks, Let covariance be the variance of the two window sub-blocks. , For all window positions Structural similarity index Perform pixel-level aggregation to obtain the coupling coefficients at adjacent time points. ,in The total number of valid cells covered by the window; all Fill the two-dimensional matrix with the data arranged in chronological order. This yields the time coupling matrix, where the row and column indices correspond to specific time points. Time coupling matrix Stored in symmetric form, with a value of 1 assigned to each diagonal position, its matrix elements satisfy... ,in ; Maximum time limit; is an integer subscript index.

7. The multi-source data-driven detection method of Claim 6, wherein, In step 4, for each time slice of the time-coupled matrix, a 3-pixel square neighborhood is established for each pixel in the matrix, and the difference between the maximum and minimum values of the index values in the neighborhood is calculated to obtain a local difference value 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, and the local difference values are calculated in the same way as the first scale, and the local difference value maps of the three scales are linearly superimposed according to weights of 0.5, 0.3 and 0.2 to obtain a fused multi-scale morphological gradient map.

8. The multi-source data-driven detection method of Claim 7, wherein, In step 4, for each time slice of the time-coupled matrix, a 3-pixel square neighborhood is established for each pixel in the matrix, and the difference between the maximum and minimum values of the index values in the neighborhood is calculated to obtain a local difference value 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, and the local difference values are calculated in the same way as the first scale, and the local difference value maps of the three scales are linearly superimposed according to weights of 0.5, 0.3 and 0.2 to obtain a fused multi-scale morphological gradient map.

9. The multi-source data-driven detection method of Claim 8, wherein, The final multi-scale morphological gradient map sequence is segmented along the time dimension using a sliding time window with a length of 5, and the window step is 1 time slice, so that there are 4 time slice overlaps between adjacent time windows to ensure the time continuity of the fractal complexity estimation; for the multi-scale morphological gradient map in each time window, the side length of the observation grid is reduced by 32 pixels, 16 pixels and 8 pixels, respectively, the number of non-zero grids occupied by the gradient map under different levels of observation grid size is counted, and the observation grid of each level and the corresponding number of non-zero grids are plotted in a scatter plot in a logarithmic coordinate system, and the slope of the least squares straight line is obtained, and the absolute value of the slope is taken as the fractal complexity value of the time window.

10. A multi-source data-driven detection system for detecting the ecological resilience of a flood detention basin for implementing the method of any one of claims 1 to 9, characterized in that, The system comprises: a data acquisition unit for acquiring optical satellite remote sensing images at different time periods before, during and after flood in the flood storage area; performing geometric correction and radiation correction 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 indexes from the time series basic data set, namely: water area mask and vegetation index; stacking the obtained water area mask and vegetation index 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 of the ecological state matrix corresponds to time; a time coupling matrix construction unit for mapping the ecological state matrix to a fractal space, using a hierarchical iterative fractal function to perform self-similar transformation and scale compression on each layer of data; in the fractal space, coupling the fractal patterns of adjacent time points, calculating the similarity of adjacent fractal patterns through a sliding window to form a time coupling matrix; a fractal complexity calculation unit for each time slice of the time coupling matrix, setting a fixed size neighborhood around each pixel, calculating the difference between the maximum and minimum values of the indexes in the neighborhood to obtain a local difference value; combining local difference values of different scales to generate a multi-scale morphological gradient map; performing time window segmentation and fractal complexity estimation on the multi-scale morphological gradient map; a resilience evaluation unit for marking the resilience value of each location by comparing the fractal complexity of the same location over time; then performing connected component analysis on all locations to obtain multiple connected regions; for each connected region, calculate the resilience value of the connected region according to its area and the resilience values of all locations in the connected region.

Images