Method for long-period monitoring of ecological restoration of substation construction based on time-series remote sensing images

Abstract

The application relates to the technical field of remote sensing monitoring, and discloses a long-period monitoring method for ecological restoration of a transformer substation construction based on time-series remote sensing images. The method comprises the following steps: extracting each pixel annual peak NDVI to obtain a recovery rate time-series matrix; performing radial gradient statistics to obtain a gradient slope sequence; calibrating an underground damage index and identifying a gradient solidification unit to generate a repair priority map; and calculating a spatial variance ratio to output an acceptance conclusion. The application solves the problems that the existing ecological restoration monitoring method based on time-series remote sensing images cannot identify the distribution of the underground damage degree by using the spatial heterogeneity of the recovery rate in the construction disturbance area and cannot determine whether the ecological function is spatially homogeneously restored.

Description

Technical Field

[0001] This application relates to the field of remote sensing monitoring technology, and in particular to a long-term monitoring method for ecological restoration during substation construction based on time-series remote sensing images. Background Technology

[0002] As a core hub node of the power system, the construction of substations inevitably disturbs the surrounding ecological environment. Activities such as vegetation clearing, land leveling, heavy machinery compaction, and underground foundation excavation all damage the surface vegetation cover and underground soil structure to varying degrees. With the continuous improvement of my country's ecological compliance supervision requirements for power infrastructure construction, long-term quantitative monitoring of the ecological restoration status of substation construction disturbance areas after completion has become a necessary part of environmental protection acceptance for projects. In existing methods, ecological restoration monitoring based on time-series remote sensing imagery mainly extracts the normalized vegetation index (NDVI) time-series data of each pixel in the construction disturbance area. The NDVI value of each pixel after completion is used as the criterion for judgment. This is combined with time-series change detection algorithms such as LandTrendr and CCDC to assess the vegetation restoration status of the disturbance area. These methods have been widely used in fields such as mining area reclamation and post-disaster recovery.

[0003] However, existing ecological restoration monitoring methods based on time-series remote sensing imagery have fundamental technical flaws. These methods independently assess whether the NDVI value has reached the recovery threshold within a single pixel dimension, treating spatial differences in recovery rates between pixels within the construction disturbance area as noise that needs to be smoothed out, rather than utilizing effective information. This results in monitoring conclusions that only reflect the apparent characteristic of whether vegetation cover has recovered, failing to distinguish between two fundamentally different states: an overall NDVI recovery but persistent soil compaction damage in the construction core area, and a truly homogeneous spatial restoration of ecological functions. Specifically, after substation construction, construction companies commonly use temporary greening methods such as planting fast-growing grasses or crops in the disturbance area. Such measures can quickly raise the NDVI value to near-background levels, causing existing single-pixel assessment methods to misjudge recovery as achieved. Meanwhile, the structural damage to the underground soil in the construction core area caused by long-term heavy machinery compaction and deep excavation cannot be identified using the NDVI value index.

[0004] Existing methods completely discard the spatial relationship information of the recovery rate of each pixel, focusing only on the absolute value of the NDVI of a single pixel, which leads to a series of progressive technical problems. The intensity of construction pressure at different spatial locations within the construction disturbance zone essentially exhibits a structural gradient that decreases from the core area to the periphery. This gradient leaves an observable imprint on the spatial distribution of surface vegetation recovery rate. However, existing methods cannot extract this spatial gradient information, and therefore cannot map the spatial heterogeneity of the observable surface recovery rate to the distribution of the invisible underground damage level, thus failing to locate hotspots of persistent underground damage. In the absence of hotspot location, existing methods also lack quantitative criteria for determining whether the spatial heterogeneity caused by construction has faded to the natural background level over time, resulting in a lack of tracking evidence for the fading process of underground damage structure in ecological restoration acceptance conclusions. Ultimately, this leaves the ecological compliance supervision of substation construction at the level of superficial assessment of surface green volume, failing to achieve spatial structural identification and long-term quantitative tracking of the degree of construction ecological damage. Summary of the Invention

[0005] This application provides a long-term monitoring method for ecological restoration of substations based on time-series remote sensing images. It solves the problems of existing ecological restoration monitoring methods based on time-series remote sensing images being unable to identify the distribution of underground damage degree by utilizing the spatial heterogeneity of restoration rate in the construction disturbance area and unable to determine whether the ecological function has achieved spatial homogeneous restoration. It realizes the inversion mapping from the surface observable restoration rate gradient to the underground invisible damage degree, as well as the quantitative determination of long-term ecological homogeneous restoration based on the spatial variance ratio sequence.

[0006] This application provides a long-term monitoring method for ecological restoration during substation construction based on time-series remote sensing imagery. The method includes: Step S1: Perform annual peak NDVI extraction on each pixel in the multi-temporal remote sensing image of the substation construction disturbance area to obtain the pixel-by-pixel recovery rate time series matrix. Step S2: Perform radial gradient statistics on the pixel-by-pixel recovery rate time series matrix with the centroid of the construction red line as the origin to obtain the radial gradient curve of the recovery rate and the corresponding gradient slope sequence. Step S3: Using the intercept value of the radial gradient curve of the recovery rate as a proxy index for the degree of underground damage, and combining it with engineering parameters, the underground damage index is obtained by linear regression calibration; the Mann-Kendall trend test is performed on the recovery rate deviation time series of each spatial unit in the gradient slope sequence, and spatial units whose deviation values ​​are consistently significantly lower than the surrounding mean are marked as gradient solidification units. The underground damage index is superimposed with the gradient solidification units to obtain a repair priority map. Step S4: Calculate the spatial variance ratio of the recovery rates of the per-pixel recovery rate time series matrix between the perturbed area and the reference area on an annual basis. When the spatial variance ratio meets the preset threshold for three consecutive years, output the ecological homogeneity restoration acceptance conclusion.

[0007] The technical solution provided in this application uses a pixel-by-pixel recovery rate time-series matrix as the core data carrier, completely preserving the spatial relationship of recovery rates between pixels within the construction disturbance zone within the matrix structure. This overcomes the inherent limitations of existing solutions that independently evaluate NDVI values ​​in a single pixel dimension and discard spatial relationship information between pixels. Based on this, radial gradient statistics are performed on the recovery rate time-series matrix with the centroid of the construction redline as the origin. This transforms the damage pressure structure within the disturbance zone, decreasing from the construction core area to the outer natural buffer zone, into a quantitatively describable radial gradient curve and gradient slope sequence of the recovery rate. This transforms spatial heterogeneity information, which was previously treated as noise and smoothed out in existing solutions, into an effective signal carrying the spatial distribution law of underground damage levels. This achieves a technological leap from observable surface measurements to invisible underground damage quantities—a leap not achieved by existing solutions. The intercept of the radial gradient curve of the recovery rate is used as a proxy indicator of the degree of underground damage in the construction core area. Combined with the engineering parameters in the construction archive, the underground damage index is obtained by linear regression calibration. The macroscopic observation of surface vegetation by satellite remote sensing and the measured data of soil compaction at the construction site are integrated in the same inversion framework, giving the underground damage index a traceable and calibrable physical basis.

[0008] Furthermore, the targeted application of the Mann-Kendall nonparametric trend test in the ecological restoration monitoring scenario of substation construction allows for the statistically reliable identification of a solidified trend in which the recovery rate deviation of each spatial unit is consistently significantly lower than the surrounding mean. The restoration priority map generated by superimposing the underground damage index with the gradient solidification units accurately locates hotspots with severe damage and insufficient natural recovery capacity, providing a spatial basis with geographic coordinate boundaries for targeted restoration decisions after completion. In terms of long-term compliance acceptance judgment, the spatial variance ratio sequence of the recovery rate between the disturbed area and the reference area replaces the single-pixel NDVI threshold judgment in the existing scheme. This elevates the evaluation dimension of ecological restoration from "whether the green volume of a single point meets the standard" to "whether the spatial heterogeneity of the disturbed area has receded to the natural background level." The judgment condition of meeting the preset threshold for three consecutive years further eliminates the interference of single-year climate anomalies on the acceptance conclusion, ensuring that the ecological homogeneous restoration acceptance conclusion is based on the statistical foundation of stable trends over many years. This fundamentally solves the technical defect of the existing scheme, which neglects spatial structure information, leading to acceptance conclusions remaining at the level of apparent green volume assessment. Attached Figure Description

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

[0010] Figure 1 This is a schematic diagram of an embodiment of the long-term monitoring method for ecological restoration of substation construction based on time-series remote sensing images in this application. Figure 2 This is a schematic diagram of the spatial distribution of the Mann-Kendall Z statistic and the identification results of the gradient solidification unit in the embodiments of this application; Figure 3 This is a schematic diagram showing the comparison between the spatial variance ratio sequence of 5 consecutive years after completion and a preset threshold in an embodiment of this application. Detailed Implementation

[0011] This application provides a long-term monitoring method for ecological restoration during substation construction based on time-series remote sensing imagery. The terms first, second, third, fourth, etc. (if present) in the specification, claims, and accompanying drawings are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments described herein can be implemented in a sequence other than that illustrated or described herein. Furthermore, the terms include or have, and any variations thereof, are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or device that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or devices.

[0012] For ease of understanding, the specific process of the embodiments of this application is described below. Please refer to [link / reference]. Figure 1 One embodiment of the long-term monitoring method for ecological restoration of substation construction based on time-series remote sensing imagery in this application includes: Step S1: Perform annual peak NDVI extraction on each pixel in the multi-temporal remote sensing image of the substation construction disturbance area to obtain the pixel-by-pixel recovery rate time series matrix. Specifically, the per-pixel recovery rate time-series matrix is ​​constructed using each raster pixel within the construction disturbance area as the row index and each natural year after completion as the column index. Each element in the matrix stores the interannual NDVI increment of that pixel in the corresponding year, i.e., the difference between the annual peak NDVI of the current year and the annual peak NDVI of the previous year. The annual peak NDVI is taken as the maximum value of the effective observations during the growing season from June to September each year. This time window is set based on the fact that the peak of temperate vegetation in the Northern Hemisphere is concentrated in this interval during the growing season, eliminating the interference of bare soil spectra on vegetation signals during the non-growing season. The completion and acceptance date is taken as zero time, ensuring that the annual increment is measured from the vegetation state after construction, eliminating the confusion of the baseline value caused by construction disturbance.

[0013] Step S2: Perform radial gradient statistics on the per-pixel recovery rate time series matrix with the centroid of the construction red line as the origin to obtain the radial gradient curve of the recovery rate and the corresponding gradient slope sequence. Specifically, the pressure intensity of substation foundation construction on the surface ecology is highest in the core excavation area and decreases towards the outer natural buffer zone, forming a concentric circle damage distribution structure. The pixels in the disturbed area are statistically divided into rings according to radial distance, with a ring size of 10m, which matches the 10m spatial resolution of Sentinel-2 imagery. This ensures that each distance ring contains a sufficient number of pixels to support the stability of the mean statistics within the ring. The mean recovery rate within the same distance ring represents the ecological recovery level at that distance. Linear regression is performed on the mean of each distance ring with radial distance as the independent variable to obtain a gradient slope sequence. A slope greater than zero indicates that the recovery rate increases from the inside to the outside, confirming the existence of the spatial structure of underground damage.

[0014] Step S3: Using the intercept value of the radial gradient curve of the recovery rate as a proxy index for the degree of underground damage, and combining it with engineering parameters, the underground damage index is obtained by linear regression calibration; the Mann-Kendall trend test is performed on the recovery rate deviation time series of each spatial unit in the gradient slope sequence, and spatial units whose deviation values ​​are consistently significantly lower than the surrounding mean are marked as gradient solidification units. The underground damage index is superimposed with the gradient solidification units to obtain the repair priority map. Specifically, the intercept value is taken as the function value of the radial gradient curve at the point where the radial distance is zero. It reflects the deviation of the recovery rate of the construction core area relative to the overall disturbance area. The larger the deviation, the more severe the compaction or backfilling damage to the underground soil in the core area. This value is calibrated by performing multiple linear regression with the frequency of heavy machinery operation and the foundation excavation depth in the construction file. The regression coefficients are fitted with least squares using the measured soil compaction data as the dependent variable to obtain the underground damage index. The determination of the gradient solidification unit adopts the standard of a Mann-Kendall test Z statistic less than -1.96, corresponding to a two-sided significance level of 0.05. The repair priority index thresholds of 0.6 and 0.3 are determined according to the principle of trisection of the underground damage index range from 0 to 1, which spatially distinguishes the three categories of high, medium and low damage levels.

[0015] Step S4: Calculate the spatial variance ratio of the recovery rates of the per-pixel recovery rate time series matrix between the perturbed area and the reference area on an annual basis. When the spatial variance ratio meets the preset threshold for three consecutive years, output the ecological homogeneity restoration acceptance conclusion.

[0016] Specifically, the preset threshold for the spatial variance ratio is the upper limit of the average interannual variation of the variance of the reference area over five consecutive years. This upper limit represents the maximum normal fluctuation range of the spatial difference in vegetation restoration rate under natural conditions. The reference area is selected from at least three areas with the same landform type as the construction disturbance area and an area of ​​not less than 5 hectares each. The average of multiple reference areas is used to eliminate the interference of the micro-topographical differences of climate in a single reference area on the determination of the threshold. The judgment condition of meeting the threshold for three consecutive years excludes the occasional occurrence of indicators caused by single-year climate anomalies. When the spatial variance ratio is not greater than the preset threshold for three consecutive years, the conclusion of ecological homogeneous restoration acceptance is output. Otherwise, the conclusion of delayed acceptance or continuous damage is determined by combining whether the proportion of the first-level restoration priority area in the total area of ​​the disturbance area in the restoration priority map exceeds 30%.

[0017] In one specific embodiment, step S1 includes: Based on the construction start date and completion acceptance date recorded in the substation construction archive, atmospheric correction and topographic radiation correction were performed on multi-temporal satellite remote sensing images with a spatial resolution of no less than 10m for five consecutive natural years before construction and five consecutive natural years after completion to obtain a time-series data stack of surface reflectance. The difference between the near-infrared band reflectance and the red band reflectance of each pixel in the surface reflectance time series data stack is divided by the sum of the two. The normalized vegetation index is calculated for each pixel in the construction disturbance area hourly to obtain the pixel-level NDVI time series set. Peak synthesis is performed on the effective observations from June to September of each calendar year in the pixel-level NDVI time series to obtain the annual peak NDVI of each pixel. Based on the annual peak NDVI of each pixel, the annual increment is calculated by taking the completion and acceptance date as the zero time point and calculating the difference between the annual peak NDVI of two adjacent years. The annual NDVI increment of each pixel for five consecutive years after completion is assembled in the spatial dimension to obtain the per-pixel recovery rate time series matrix.

[0018] Specifically, the surface reflectance time-series data stack is a three-dimensional data structure formed by sorting all pre-processed images from five consecutive calendar years before construction and five consecutive calendar years after completion by timestamp and then cropping them to the extended range of the construction boundary. The three dimensions are pixel column coordinates, pixel row coordinates, and image acquisition time. The data from the five years before construction is used to establish a baseline reference for vegetation growth, while the data from the five years after completion covers the complete long-term process of ecological restoration. The symmetrical setting of the two time windows ensures that the baseline and restoration period data are comparable in terms of climate background. The requirement of a spatial resolution of not less than 10m is determined based on the actual situation that the area of ​​the substation construction disturbance zone is usually between several hectares and tens of hectares. A resolution lower than this will lead to severe pixel mixing at the disturbance boundary, making it impossible to support subsequent pixel-level gradient analysis. Atmospheric correction uses the 6S radiative transfer model to eliminate the systematic bias of atmospheric scattering and absorption on surface reflectance. Topographic radiative correction uses the C-correction method to eliminate the uneven irradiance caused by differences in slope aspect in mountainous areas. Both corrections are necessary preprocessing steps for accurate extraction of pixel-level surface reflectance. The difference between the corrected near-infrared reflectance and the red reflectance, divided by the sum of the two, yields the normalized vegetation index for each pixel. This index has high sensitivity to green vegetation coverage in the low to medium coverage range and is suitable for monitoring scenarios where vegetation gradually recovers from zero coverage after construction disturbance.

[0019] The time window from June to September within each calendar year in the pixel-level NDVI time series is determined based on the pattern that temperate vegetation in the Northern Hemisphere completes the full growth process from leaf expansion to canopy closure during this period. The maximum effective observation value within this window is taken as the annual peak NDVI, which eliminates low-value anomalies introduced by cloud pollution imagery and filters out interference from bare soil and litter spectra during the non-growing season on vegetation signals. The annual peak NDVI of each pixel is calculated by taking the completion and acceptance date as zero and comparing the difference between two adjacent years. The difference represents the vegetation recovery rate of that pixel in that year. Positive values ​​represent net increase in vegetation cover, negative values ​​represent degradation, and zero values ​​represent no change. The interannual NDVI increments of each pixel for five consecutive years after completion are arranged in the row direction by pixel coordinates and in the column direction by year to form a two-dimensional matrix, which yields the pixel-by-pixel recovery rate time series matrix. This matrix completely preserves the recovery rate values ​​of each spatial location within the perturbed area in each year. The spatial relative relationships between the elements in the matrix are the information basis for subsequent radial gradient statistics.

[0020] In one specific embodiment, step S2 includes: For each pixel in the per-pixel recovery rate time series matrix, a local deviation is calculated by comparing the average recovery rate of all pixels in a circular neighborhood with a radius of 3 pixels with the recovery rate of the corresponding pixel itself, and the annual recovery rate deviation value of each pixel is obtained. Based on the annual recovery rate deviation of each pixel, the pixels in the disturbance area are divided into rings according to radial distance with the centroid of the construction red line as the origin and a step size of 10m. The recovery rate deviation of all pixels in the same distance ring is statistically analyzed to obtain the radial gradient curve of the recovery rate. A linear regression was performed on the radial gradient curve of the recovery rate with radial distance as the independent variable and the mean of the recovery rate deviation loop as the dependent variable to obtain the gradient slope sequence; The gradient decay rate is obtained by performing a first-order linear regression on the gradient slope sequence for five consecutive years after completion, with the year as the independent variable.

[0021] Specifically, the annual recovery rate deviation value for each pixel is obtained by taking the arithmetic mean of the recovery rates of all pixels within a circular neighborhood of 3 pixels centered on the pixel in the pixel-by-pixel recovery rate time-series matrix, and then subtracting the pixel's own recovery rate from this mean. The neighborhood of 3 pixels corresponds to 30 meters above ground in a 10m resolution image. This scale matches the typical operational unit scale of substation construction disturbances, capturing the local recovery rate differences between the core construction area and the surrounding buffer zone without mixing in the recovery signals from distant natural vegetation due to an excessively large neighborhood. A negative deviation value indicates that the pixel's recovery rate is lower than the local neighborhood mean, while a positive deviation value indicates that it is higher than the local neighborhood mean. The spatial distribution of the deviation value directly reflects the degree of local heterogeneity of the recovery rate within the disturbance area. Using the centroid of the construction red line as the origin, each pixel in the disturbance area is divided into rings according to its Euclidean distance from the centroid, with a step size of 10m, consistent with the spatial resolution of the Sentinel-2 image. This ensures that each distance ring contains at least dozens of pixels, supporting the robustness of the mean statistics within the ring. The arithmetic mean of the recovery rate deviation values ​​of all pixels within the same distance ring is taken to obtain the mean value within the ring at that distance. The mean values ​​of each distance ring are arranged and connected according to the radial distance to form a radial gradient curve of the recovery rate.

[0022] A univariate linear regression was performed on the radial gradient curve of the recovery rate, with radial distance as the independent variable on the horizontal axis and the mean deviation of the recovery rate of each distance loop as the dependent variable on the vertical axis. The regression slope is the gradient slope. A slope greater than zero indicates that the recovery rate increases with increasing distance from the construction core area, and the recovery rate in the construction core area is consistently lower than that in the surrounding areas, confirming the existence of spatial structural damage in the subsurface. A slope not greater than zero indicates that the spatial distribution of the recovery rate in the disturbed area tends to be homogeneous in that year. The gradient slopes calculated year by year for five consecutive years after completion form a gradient slope sequence. A first-order linear regression was performed on this sequence with the year as the independent variable. The regression slope is the gradient decay rate. A gradient decay rate less than zero indicates that the gradient weakens over time and the ecosystem evolves towards spatial homogeneity. A gradient decay rate not less than zero indicates that the gradient has not receded and the structural spatial heterogeneity caused by construction continues to exist. This indicator transforms the temporal change of the observable recovery rate on the surface into a quantitative measure of the subsurface damage decay process, which is fundamentally different from the existing approach that only focuses on the absolute value change of NDVI of a single pixel and cannot reflect the spatial structural decay law of damage.

[0023] In one specific embodiment, step S3 uses the intercept of the radial gradient curve of the recovery rate as a proxy index for the degree of subsurface damage, and obtains the subsurface damage index through linear regression calibration in combination with engineering parameters, including: Extract the function intercept value of the radial gradient curve of the recovery rate at the point where the radial distance is zero to obtain the underground damage proxy quantity in the construction core area; Multiple linear regression calibration was performed based on the underground damage proxy volume in the construction core area and the frequency of heavy machinery operation and foundation excavation depth recorded in the construction archives. The regression coefficients were fitted with least squares using the measured soil compaction data on site as the dependent variable to obtain the regression coefficient set. Based on the underground damage proxy quantity and regression coefficient set in the construction core area, the underground damage degree inversion calculation is performed on each spatial unit in the disturbance zone to obtain the underground damage index.

[0024] Specifically, the radial gradient curve of the recovery rate is fitted with radial distance as the horizontal axis and the mean deviation of the recovery rate of each distance loop as the vertical axis. A straight line is obtained through linear regression. The intercept of this line at the point where the horizontal axis is zero is the surrogate quantity of underground damage in the construction core area. The physical meaning of the intercept value is the predicted value of the recovery rate deviation at the centroid of the construction boundary. The greater the negative bias of this value, the stronger the inhibition of the recovery rate of the construction core area relative to the overall disturbed area, corresponding to the more severe structural damage to the underground soil caused by heavy machinery compaction or construction waste landfill. The rationale for using the intercept value as a surrogate quantity of underground damage is that underground soil compaction directly inhibits root growth and soil permeability, resulting in the surface vegetation recovery rate being consistently lower than that of the surrounding undamaged areas. There is a monotonic correspondence between the two, and the absolute value of the intercept value quantifies the degree of inhibition. The amount of underground damage in the core construction area is taken as one of the independent variables. Together with the frequency of heavy machinery operation and the foundation excavation depth recorded in the construction archive, they form a set of independent variables for multiple linear regression. The measured soil compaction data collected on site after construction is taken as the dependent variable. The regression coefficients corresponding to each independent variable are obtained by minimizing the sum of squared regression residuals through least squares fitting. The three independent variables correspond to three regression coefficients plus a constant term, which together form a set of regression coefficients.

[0025] The underground damage index is calculated by substituting the underground damage proxy quantity of the construction core area of ​​each spatial unit within the disturbance zone into the linear equation determined by the regression coefficient set. The calculation result is normalized and mapped to a dimensionless interval of 0 to 1, where 0 represents no damage to the underground ecological structure and consistency with the natural background, and 1 represents the most severe underground damage, corresponding to extreme working conditions such as long-term heavy machinery compaction or deep excavation without backfilling. The normalization process uses the minimum and maximum values ​​of the output values ​​of the linear equations of all spatial units within the disturbance zone as boundaries, and performs max-min normalization on the output values ​​of each unit to ensure that the relative damage level ranking of the index within the disturbance zone is consistent with the absolute calculated value. The spatial distribution of the underground damage index directly reflects the spatial heterogeneity of the damage to the underground ecological structure caused by construction disturbance. The high value area of ​​the damage index is concentrated in the core area of ​​foundation excavation where heavy machinery is frequently used, and decreases with the decrease of operation intensity in the outer periphery. This spatial distribution feature is consistent with the structure of the recovery rate from low in the inside to high in the outside as shown by the radial gradient curve. The consistency of the spatial pattern of the two independent sources of indicators verifies the rationality of the technical approach of using the intercept value as the underground damage proxy.

[0026] In one specific embodiment, step S3 involves performing a Mann-Kendall trend test on the recovery rate deviation time series of each spatial unit in the gradient slope sequence, and marking spatial units with deviation values ​​that are consistently significantly lower than the surrounding mean as gradient solidification units, including: A time-series vector was constructed for the recovery rate deviation values ​​of each spatial unit in the gradient slope sequence over five consecutive years after completion. The Mann-Kendall nonparametric trend test was performed on the time-series vector with a significance level of 0.05 as the test threshold to obtain the Z statistic for each spatial unit. Based on the Z-statistic, spatial cells whose Z-statistic consistently satisfies a value less than -1.96 are identified as spatial cells whose recovery rate deviation is consistently significantly lower than the surrounding mean. The identification results are then spatially connected and merged to obtain gradient solidification cells.

[0027] Specifically, the recovery rate deviation values ​​of each spatial unit for five consecutive years after completion are arranged in chronological order to form a time-series vector. This vector has a length of 5, and its elements are the recovery rate deviation values ​​of the spatial unit from year 1 to year 5 after completion. The Mann-Kendall nonparametric trend test compares the pairwise relationships of all elements in the time-series vector, and calculates the difference between the logarithm of a later value being greater than the logarithm of a later value being less than the logarithm of a later value. This difference is then standardized to obtain the Z-statistic. A negative Z-statistic indicates that the time-series vector shows a downward trend, meaning that the recovery rate deviation value of the spatial unit continues to decrease over time. The larger the absolute value of the Z-statistic, the more significant the downward trend. Using a significance level of 0.05 as the test threshold, the corresponding two-sided critical value for a normal distribution is -1.96. A Z-statistic less than -1.96 indicates that the trend of the recovery rate deviation value of the spatial unit being significantly lower than the surrounding mean is statistically valid. This critical value is a universally accepted statistical standard, ensuring that the conclusion is reproducible and verifiable. The Mann-Kendall test does not require time series data to follow a specific distribution and is applicable to short time series with a sample size of 5. This is the basis for choosing this method instead of the nonlinear regression significance test.

[0028] All spatial units with a Z-statistic less than -1.96 are spatially connected and merged. Specifically, a binary mask is constructed for pixels meeting the criteria. An 8-neighborhood connectivity analysis is performed on pixels with a value of 1 in the mask. Spatially adjacent pixels meeting the criteria are grouped into the same patch. Isolated patches with an area smaller than a minimum mapping unit (i.e., one pixel area) are removed, and continuous patches with an area not less than one pixel are retained. The set of all retained patches is the gradient solidification unit. The physical meaning of a gradient solidification unit is a spatial set where the recovery rate within the construction disturbance zone is consistently and systematically lower than the surrounding area, and this suppression trend is statistically significant. This set spatially corresponds to the area with the most severe underground soil compaction or backfill damage and serves as the spatial carrier for subsequent overlay of the underground damage index to generate a repair priority map. The connectivity merging operation integrates the discrete single-pixel judgment results into continuous patches with geographical boundaries, giving the output results spatial accuracy that directly guides on-site repair operations.

[0029] Figure 2This diagram illustrates the spatial distribution of the Mann-Kendall Z-statistic and the identification results of gradient solidification units in this embodiment. The shades of gray in the diagram represent the magnitude of the Z-statistic for the recovery rate deviation over a continuous five-year period after completion of each spatial unit within the construction disturbance zone. Darker colors indicate a larger absolute value of the Z-statistic and a more significant trend of the recovery rate deviation consistently falling below the surrounding mean. The dashed boundary in the diagram represents the critical contour line where the Z-statistic equals -1.96. The area inside the critical line (i.e., the light-colored rectangular area in the center of the diagram) represents the gradient solidification units identified after connecting and merging units whose Z-statistic consistently meets the criterion of being less than -1.96. This area spatially coincides with the heavy machinery operation range in the core construction area, confirming the rationality of the technical approach of determining the area of ​​persistent underground damage based on the trend of the recovery rate deviation.

[0030] In one specific embodiment, step S3 involves superimposing the underground damage index with the gradient solidification unit to obtain a repair priority map, including: A spatial overlay operation is performed on the underground damage index and the gradient solidification unit. The underground damage index values ​​of each spatial unit within the gradient solidification unit are retained, while the underground damage index values ​​of non-gradient solidification units are set to zero, thus obtaining the spatial distribution of the repair priority index. Based on the spatial distribution of the restoration priority index, spatial units with a restoration priority index greater than 0.6 are designated as first-level restoration priority areas, spatial units with a restoration priority index between 0.3 and 0.6 are designated as second-level restoration priority areas, and spatial units with a restoration priority index not greater than 0.3 are designated as natural recovery monitoring areas. Spatial merging and boundary extraction are performed on the three types of areas to obtain a restoration priority map.

[0031] Specifically, spatial overlay operations are performed at the pixel level. For each spatial unit within the disturbance zone, it is determined whether it belongs to a gradient solidification unit. Spatial units belonging to gradient solidification units retain their original subsurface damage index values, while those not belonging to gradient solidification units have their subsurface damage index values ​​forcibly set to zero. The basis for this zeroing operation is that the recovery rate deviation trend of non-gradient solidification units has not statistically reached the standard of sustained significant inhibition. The subsurface damage of these units is still within the range of natural recovery and does not need to be included in the priority sequence for artificial intervention. Setting their index to zero allows repair resources to be concentrated on areas where damage is continuously solidified. After the overlay operation is completed, each spatial unit within the disturbance zone holds a repair priority index value. Gradient solidification units hold the original values ​​from the subsurface damage index, while non-gradient solidification units hold zero values. The continuous spatial distribution of the index values ​​of all units constitutes the spatial distribution of the repair priority index.

[0032] The thresholds for the restoration priority index, 0.6 and 0.3, are determined based on the principle of trisecting the normalized subsurface damage index range of 0 to 1. The 0-1 range is evenly divided into three segments: below 0.3 corresponds to low damage, 0.3 to 0.6 to medium damage, and above 0.6 to high damage. These three segments correspond to three management response intensities: natural recovery monitoring area, secondary restoration priority area, and primary restoration priority area. Spatial merging merges adjacent spatial units of the same category into continuous patches, eliminating single-pixel-level category islands. Boundary extraction generates vector boundary lines for the outer contours of each merged patch, enabling the output restoration priority map to have the spatial accuracy to directly overlay engineering drawings and guide the delineation of on-site restoration work areas. The spatial location and area of ​​the primary restoration priority area in the restoration priority map directly correspond to the work areas requiring manual soil improvement or waste removal. The secondary restoration priority area corresponds to areas requiring assisted vegetation replanting. The natural recovery monitoring area only requires continuous remote sensing tracking without human intervention. The delineation of these three types of areas transforms the temporal remote sensing analysis conclusions into operable spatial management units.

[0033] In one specific embodiment, step S4 includes: Spatial variance calculation is performed on the annual recovery rate values ​​of all pixels in the construction disturbance area in the pixel-by-pixel recovery rate time series matrix to obtain the annual spatial variance sequence of the disturbance area; Spatial variance calculation is performed on the annual recovery rate values ​​of all pixels in the reference area in the pixel-by-pixel recovery rate time series matrix. The arithmetic mean of the annual spatial variances of at least 3 reference areas is taken to obtain the annual spatial variance sequence of the reference area. Based on the annual spatial variance sequence of the disturbance area and the annual spatial variance sequence of the reference area, the ratio of the two is calculated for each of the five consecutive years after completion to obtain the spatial variance ratio sequence. Based on the spatial variance ratio sequence, determine whether there is a spatial variance ratio in the sequence that is not greater than a preset threshold for three consecutive years. If so, output the ecological homogeneity restoration acceptance conclusion; otherwise, output the delayed acceptance conclusion based on the restoration priority map.

[0034] Specifically, the calculation of the annual spatial variance sequence of the disturbed area is performed for each natural year after completion. The recovery rate value of all pixels within the construction boundary is taken for that year. The arithmetic mean of the squares of the differences between the recovery rate values ​​of all pixels in that year and the arithmetic mean of the recovery rate of pixels in the disturbed area is calculated; this is the spatial variance of the disturbed area for that year. The annual spatial variances for five consecutive years after completion are arranged by year to obtain the annual spatial variance sequence of the disturbed area. The annual spatial variance sequence of the reference area is calculated using the same rules for at least three reference areas. Reference areas are selected as areas with the same landform type as the construction disturbed area, unaffected by construction activities before and after the completion and acceptance date, and with an area of ​​at least 5 hectares each. The lower area limit ensures that the number of pixels in the reference area is sufficient to support the robustness of the spatial variance statistics. The requirement of at least three reference areas eliminates statistical bias caused by local climatic and micro-topographical differences in a single reference area. The arithmetic mean of the annual spatial variances of three or more reference areas is used to obtain the annual spatial variance sequence of the reference area. For each of the five consecutive years following completion, the spatial variance of the disturbed area is divided by the spatial variance of the reference area to obtain the spatial variance ratio for that year. The spatial variance ratios for the five years are arranged by year to form a spatial variance ratio sequence. This sequence reflects the dynamic change over time of the amplification factor of the spatial heterogeneity of the recovery rate of the disturbed area relative to the natural background.

[0035] The preset threshold is determined by taking the upper limit of the average interannual variation of the annual spatial variance sequence of the reference area over five consecutive years. This upper limit represents the maximum normal fluctuation range of the spatial difference in vegetation restoration rate under natural conditions without construction disturbance. This is used as the quantitative boundary to determine whether the spatial heterogeneity of the disturbed area has receded to the level of the natural background. This threshold is derived from the measured data of the reference area rather than being a fixed value preset by humans, ensuring that the threshold is compatible with the natural background characteristics of the specific construction area. A yearly scan of the spatial variance ratio sequence is performed to determine whether there are three consecutive years in which the spatial variance ratio is not greater than a preset threshold. The three-year consecutive three-year criterion excludes occasional occurrences caused by single-year climate anomalies (such as abnormally dry or wet years), ensuring that the acceptance conclusion is based on a stable trend over many years. When the three-year consecutive three-year compliance condition is met, an ecological homogeneous restoration acceptance conclusion is output. This conclusion indicates that the spatial heterogeneity of the restoration rate in the construction disturbance area has receded to a level that is not substantially different from the natural background. When the condition is not met, the proportion of the primary restoration priority area to the total area of ​​the disturbance area in the restoration priority map is considered. If it exceeds 30%, a persistent damage delayed acceptance conclusion is output; if it does not exceed 30%, a normal delayed acceptance conclusion is output. The two types of delayed conclusions correspond to differentiated follow-up monitoring and intervention plans.

[0036] Figure 3This is a schematic diagram comparing the spatial variance ratio sequence over five consecutive years after completion with a preset threshold in an embodiment of this application. The broken line in the diagram represents the spatial variance ratio sequence formed by the ratio of the annual spatial variance of the disturbed area to the annual spatial variance of the reference area. The horizontal dashed line represents the preset threshold of 1.20, determined based on the upper limit of the average annual variance variation range of the reference area. The diagonally filled intervals represent the years in which the spatial variance ratio exceeds the preset threshold, and the dotted intervals represent the years in which the spatial variance ratio meets the preset threshold. The starting point for continuous compliance judgment marked in the diagram corresponds to the 4th year. The spatial variance ratios for both the 4th and 5th years are not greater than the preset threshold. Considering that the 3rd year still exceeds the threshold, the continuous compliance judgment condition in this embodiment is only met for two years between the 4th and 5th years. This indicates that monitoring needs to be extended to the 6th year to confirm whether the acceptance judgment condition of three consecutive years of compliance is met.

[0037] In one specific embodiment, in the calculation of the annual spatial variance sequence of the disturbance area, the average of the squares of the differences between the recovery rate values ​​of all pixels within the boundary and the mean recovery rate of the disturbance area is taken, with the construction red line boundary as the range, to obtain the annual spatial variance sequence of the disturbance area.

[0038] In the calculation of the annual spatial variance sequence of the reference area, at least three areas with the same landform type as the construction disturbance area, which were not affected by construction activities during the construction period and have an area of ​​not less than 5 hectares are selected as the reference area. The arithmetic mean of the annual spatial variance of each reference area is calculated to obtain the annual spatial variance sequence of the reference area.

[0039] In determining the preset threshold, the annual variation amplitude of the spatial variance sequence of the reference area is statistically analyzed. The upper limit of the mean of the annual variation amplitude of the reference area itself is used as the upper limit of the acceptable spatial heterogeneity amplification factor of the perturbation area. The upper limit is determined as the preset threshold of the spatial variance ratio.

[0040] Specifically, the spatial range of the annual spatial variance sequence of the disturbance area is based on the construction red line vector boundary registered in the substation construction archives. The construction red line boundary is the legally confirmed outer contour of the construction land in the engineering permit documents. Using this boundary as the calculation range ensures that the spatial variance statistics cover all plots directly affected by construction activities and do not include the recovery rate data of natural vegetation outside the boundary. The mean recovery rate of the disturbance area is taken as the arithmetic mean of the recovery rate values ​​of all pixels within the construction red line boundary in that year. The arithmetic mean of the squares of the differences between the recovery rate values ​​of each pixel and this mean is the spatial variance of the disturbance area in that year. This calculation rule is the standard form of unbiased variance estimation in statistics. The lower limit of the reference area of ​​5 hectares is set based on the fact that in a 10m spatial resolution image, 5 hectares corresponds to approximately 500 pixels. This number can reduce the standard error of the spatial variance estimation to within 5% of the mean, meeting the statistical robustness requirement for comparison with the spatial variance of the disturbed area. The requirement of no less than 3 locations ensures that the reference area sample covers natural vegetation patches in different directions around the construction disturbance area, eliminating the systematic bias of single-direction climate micro-topographic differences on the spatial variance of the reference area.

[0041] The derivation of the preset threshold is based on the annual spatial variance sequence of the reference area. The absolute value of the difference between two consecutive years of annual spatial variance in the reference area is calculated, yielding four interannual variation amplitude values. The arithmetic mean of these four values ​​is taken as the mean interannual variance variation amplitude of the reference area. This mean reflects the normal interannual fluctuation level of the spatial difference in vegetation restoration rate under natural conditions without construction disturbance. The upper limit is set by adding its own standard deviation to the mean interannual variance variation amplitude of the reference area. This upper limit corresponds to the statistical upper bound of the normal fluctuation range of the reference area. The criterion is that the ratio of the spatial variance of the disturbed area to the spatial variance of the reference area should not exceed this upper limit. That is, when the amplification factor of the spatial heterogeneity of the disturbed area does not exceed the upper bound of the natural fluctuation of the reference area, the spatial distribution of the restoration rate of the disturbed area is judged to have receded to the natural background level. This threshold is gradually derived from measured data of the reference area and is directly linked to the natural background characteristics of specific construction areas. Substation construction areas with different geographical locations and different vegetation types will obtain preset thresholds adapted to their natural backgrounds, rather than uniformly using fixed values.

[0042] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A long-term monitoring method for ecological restoration during substation construction based on time-series remote sensing imagery, characterized in that, The method includes: Step S1: Perform annual peak NDVI extraction on each pixel in the multi-temporal remote sensing image of the substation construction disturbance area to obtain the pixel-by-pixel recovery rate time series matrix. Step S2: Perform radial gradient statistics on the pixel-by-pixel recovery rate time series matrix with the centroid of the construction red line as the origin to obtain the radial gradient curve of the recovery rate and the corresponding gradient slope sequence. Step S3: Using the intercept value of the radial gradient curve of the recovery rate as a proxy index for the degree of underground damage, and combining it with engineering parameters, the underground damage index is obtained by linear regression calibration; the Mann-Kendall trend test is performed on the recovery rate deviation time series of each spatial unit in the gradient slope sequence, and spatial units whose deviation values ​​are consistently significantly lower than the surrounding mean are marked as gradient solidification units. The underground damage index is superimposed with the gradient solidification units to obtain a repair priority map. Step S4: Calculate the spatial variance ratio of the recovery rates of the per-pixel recovery rate time series matrix between the perturbed area and the reference area on an annual basis. When the spatial variance ratio meets the preset threshold for three consecutive years, output the ecological homogeneity restoration acceptance conclusion.

2. The long-term monitoring method for ecological restoration of substations based on time-series remote sensing imagery according to claim 1, characterized in that, Step S1 includes: Based on the construction start date and completion acceptance date recorded in the substation construction archive, atmospheric correction and topographic radiation correction were performed on multi-temporal satellite remote sensing images with a spatial resolution of no less than 10m for five consecutive natural years before construction and five consecutive natural years after completion to obtain a time-series data stack of surface reflectance. The difference between the near-infrared band reflectance and the red band reflectance of each pixel in the surface reflectance time series data stack is divided by the sum of the two. The normalized vegetation index is calculated for each pixel in the construction disturbance area hourly to obtain the pixel-level NDVI time series set. Peak synthesis is performed on the effective observations from June to September of each natural year in the pixel-level NDVI time series to obtain the annual peak NDVI of each pixel. Based on the annual peak NDVI of each pixel, the annual increment is calculated by taking the completion and acceptance date as zero time and calculating the difference between the annual peak NDVI of two adjacent years. The annual NDVI increments of each pixel for five consecutive years after completion are assembled in the spatial dimension to obtain the per-pixel recovery rate time series matrix.

3. The long-term monitoring method for ecological restoration of substation construction based on time-series remote sensing imagery according to claim 1, characterized in that, Step S2 includes: For each pixel in the per-pixel recovery rate time series matrix, a local deviation calculation is performed on the difference between the average recovery rate of all pixels in a circular neighborhood with a radius of 3 pixels and the recovery rate of the corresponding pixel itself, to obtain the annual recovery rate deviation value of each pixel. Based on the annual recovery rate deviation value of each pixel, the pixels in the disturbance area are divided into rings according to radial distance with the centroid of the construction red line as the origin and a step size of 10m. The recovery rate deviation value of all pixels in the same distance ring is statistically analyzed to obtain the radial gradient curve of the recovery rate. A linear regression was performed on the radial gradient curve of the recovery rate with radial distance as the independent variable and the mean of the recovery rate deviation loop as the dependent variable to obtain the gradient slope sequence. Based on the gradient slope sequence, a first-order linear regression was performed on the gradient slope over five consecutive years after completion, with the year as the independent variable, to obtain the gradient decay rate.

4. The long-term monitoring method for ecological restoration of substation construction based on time-series remote sensing imagery according to claim 1, characterized in that, In step S3, the intercept value of the radial gradient curve of the recovery rate is used as a proxy index for the degree of subsurface damage. Combined with engineering parameters, a subsurface damage index is obtained through linear regression calibration, including: Extract the function intercept value of the radial gradient curve of the recovery rate at the point where the radial distance is zero to obtain the underground damage proxy quantity in the construction core area; Based on the underground damage proxy amount in the construction core area and the frequency of heavy machinery operation and foundation excavation depth recorded in the construction archive, multiple linear regression calibration was performed. The regression coefficients were fitted with least squares using the measured soil compaction data as the dependent variable to obtain the regression coefficient set. Based on the underground damage proxy quantity of the construction core area and the regression coefficient set, the underground damage degree inversion calculation is performed on each spatial unit in the disturbance area to obtain the underground damage index.

5. The long-term monitoring method for ecological restoration of substations based on time-series remote sensing imagery according to claim 4, characterized in that, In step S3, the Mann-Kendall trend test is performed on the recovery rate deviation time series of each spatial unit in the gradient slope sequence, and spatial units with deviation values ​​that are consistently significantly lower than the surrounding mean are marked as gradient solidification units, including: A time-series vector is constructed from the recovery rate deviation values ​​of each spatial unit in the gradient slope sequence for five consecutive years after completion. The Mann-Kendall nonparametric trend test is performed on the time-series vector with a significance level of 0.05 as the test threshold to obtain the Z statistic of each spatial unit. Based on the Z statistic, spatial cells whose Z statistic consistently satisfies a value less than -1.96 are identified as spatial cells whose recovery rate deviation is consistently significantly lower than the surrounding mean. The identification results are then spatially connected and merged to obtain gradient solidification cells.

6. The long-term monitoring method for ecological restoration of substation construction based on time-series remote sensing imagery according to claim 5, characterized in that, In step S3, the underground damage index is superimposed with the gradient solidification unit to obtain a repair priority map, including: A spatial superposition operation is performed on the underground damage index and the gradient solidification unit. The underground damage index values ​​of each spatial unit in the gradient solidification unit are retained, and the underground damage index values ​​of non-gradient solidification units are set to zero to obtain the spatial distribution of the repair priority index. Based on the spatial distribution of the repair priority index, spatial units with a repair priority index greater than 0.6 are designated as first-level repair priority areas, spatial units with a repair priority index between 0.3 and 0.6 are designated as second-level repair priority areas, and spatial units with a repair priority index not greater than 0.3 are designated as natural recovery monitoring areas. Spatial merging and boundary extraction are performed on the three types of areas to obtain a repair priority map.

7. The long-term monitoring method for ecological restoration of substation construction based on time-series remote sensing imagery according to claim 1, characterized in that, Step S4 includes: Spatial variance calculation is performed on the annual recovery rate values ​​of all pixels in the construction disturbance area in the pixel-by-pixel recovery rate time series matrix to obtain the annual spatial variance sequence of the disturbance area; Based on the annual recovery rate values ​​of all pixels in the reference area in the pixel-by-pixel recovery rate time series matrix, spatial variance calculation is performed, and the arithmetic mean of the annual spatial variances of no less than 3 reference areas is taken to obtain the annual spatial variance sequence of the reference area. Based on the annual spatial variance sequence of the disturbance region and the annual spatial variance sequence of the reference region, the ratio of the two is calculated for each of the five consecutive years after completion to obtain the spatial variance ratio sequence. Based on the spatial variance ratio sequence, determine whether there is a spatial variance ratio in the spatial variance ratio sequence that is not greater than a preset threshold for three consecutive years. If so, output an acceptance conclusion for ecological homogeneous restoration; otherwise, output a delayed acceptance conclusion based on the restoration priority map.

8. The long-term monitoring method for ecological restoration of substation construction based on time-series remote sensing imagery according to claim 7, characterized in that, In the calculation of the annual spatial variance sequence of the disturbance area, the average of the squares of the differences between the recovery rate values ​​of all pixels within the boundary and the average recovery rate of the disturbance area is taken, with the construction red line boundary as the range, to obtain the annual spatial variance sequence of the disturbance area.

9. The long-term monitoring method for ecological restoration of substation construction based on time-series remote sensing imagery according to claim 8, characterized in that, In the calculation of the annual spatial variance sequence of the reference area, at least three areas with the same landform type as the construction disturbance area, which were not affected by construction activities during the construction period and have an area of ​​not less than 5 hectares are selected as the reference area. The arithmetic mean of the annual spatial variance of each reference area is calculated to obtain the annual spatial variance sequence of the reference area.

10. The long-term monitoring method for ecological restoration of substation construction based on time-series remote sensing imagery according to claim 9, characterized in that, In determining the preset threshold, the annual variation amplitude of the spatial variance sequence of the reference area is statistically analyzed, and the upper limit of the mean of the annual variation amplitude of the reference area itself is used as the upper limit of the acceptable spatial heterogeneity amplification factor of the perturbation area. The upper limit value is determined as the preset threshold of the spatial variance ratio.

Images