Artificial forest degradation area extraction method and system based on high-resolution remote sensing

Abstract

The application discloses a method and system for extracting artificial forest degradation area based on high-resolution remote sensing, and relates to the technical field of remote sensing image processing, which comprises the following steps: acquiring multi-temporal remote sensing images and extracting the initial boundary of artificial forest; constructing a double-sided buffer zone on both sides of the initial boundary, layering sampling along the normal direction to obtain a feature profile, determining the spatial position of the boundary transition zone by gradient mutation characteristics, and extracting the width distribution data; dividing the initial boundary into multiple boundary segments, calculating the difference degree of the transition zone width values of adjacent boundary segments, and marking the abnormal boundary segment when the difference degree exceeds the preset threshold; establishing a detection path in the artificial forest along the normal direction from the abnormal boundary segment as the starting point, calculating the correlation degree between the vegetation characteristics and the transition zone width value; determining the detection endpoint position according to the correlation degree attenuation characteristics, taking the detection path coverage area as the degradation candidate area, and outputting the degradation area through connectivity merging and morphological screening. The application realizes accurate and reliable extraction of artificial forest degradation area.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing technology, specifically to a method and system for extracting degraded areas of plantations based on high-resolution remote sensing. Background Technology

[0002] Planted forests serve as crucial ecological barriers and timber resource bases, and their health directly impacts ecological security and economic benefits. Due to factors such as pest and disease infestation, climate change, and soil degradation, planted forests often experience localized degradation during their growth. Timely and accurate monitoring and identification of degraded areas are of great significance for forestry resource management and ecological protection.

[0003] Traditional monitoring of plantation degradation mainly relies on ground surveys and manual patrols, which suffers from low efficiency, high cost, and limited coverage, making it difficult to achieve large-scale dynamic monitoring. With the development of remote sensing technology, plantation monitoring methods based on remote sensing imagery have gradually become a research hotspot. Existing remote sensing monitoring methods mainly identify degraded areas through vegetation index change detection and texture feature analysis, but these methods often treat the entire forest area as a unified analysis object, ignoring the special characteristics of the plantation boundary areas.

[0004] As the boundary between plantations and their surrounding environment, the boundary zone of a plantation exhibits distinct transitional ecological characteristics. This boundary area is more directly affected by the external environment and often becomes the first and most pronounced area of ​​degradation. Existing monitoring methods lack in-depth analysis of the characteristics of this boundary transition zone, failing to effectively utilize anomalies in the boundary area to guide the identification of degraded areas. Furthermore, existing methods lack correlation analysis with changes in boundary characteristics when determining the spatial extent of degraded areas, resulting in insufficient accuracy and reliability in degraded area extraction. Summary of the Invention

[0005] The purpose of this invention is to provide a method and system for extracting degraded areas of plantations based on high-resolution remote sensing, aiming to solve at least one of the technical problems existing in the prior art.

[0006] The technical solution of this invention is: a method for extracting degraded areas of plantations based on high-resolution remote sensing, comprising the following steps: Acquire multi-temporal remote sensing images of the target area and perform preprocessing to extract the initial boundaries of the plantation forest; A double-sided buffer zone is constructed on both sides of the initial boundary. Vegetation features are sampled in layers along the normal direction of the initial boundary to obtain a feature profile that crosses the initial boundary. The spatial location of the boundary transition zone is determined by the gradient abrupt change features of the feature profile, and the width distribution data of the boundary transition zone is extracted. The initial boundary is divided into multiple boundary segments. For each boundary segment, the corresponding transition band width value is extracted from the width distribution data. The degree of difference between the transition band width values ​​of adjacent boundary segments is calculated. When the degree of difference exceeds the preset difference threshold, the current boundary segment is marked as an abnormal boundary segment. Starting from the abnormal boundary segment, a detection path is established inside the plantation along the normal direction of the abnormal boundary segment. Vegetation features are extracted on the detection path and the correlation between the vegetation features and the width of the transition zone corresponding to the abnormal boundary segment is calculated. The detection endpoint is determined based on the decay characteristics of the correlation. The area covered by the detection path from the abnormal boundary segment to the detection endpoint is taken as the degradation candidate area. The degradation candidate area is merged based on connectivity and morphological screening to output the degradation area of ​​the plantation forest.

[0007] Acquire multi-temporal remote sensing images of the target area and perform preprocessing to extract the initial boundaries of the plantation, including: Remote sensing image data of the target area at different time phases are acquired and radiometrically corrected and geometrically registered. Vegetation index time series data are extracted from the registered remote sensing image data. Seasonal decomposition of vegetation index time series data yields trend and fluctuation components. The spatial decay direction of the trend component and the spatial amplitude distribution of the fluctuation component are extracted. The directional consistency between the spatial decay direction and the spatial amplitude distribution is calculated. Candidate boundary regions are determined based on the spatial abrupt change locations of the directional consistency. Phase information of vegetation index time series data is extracted within the candidate boundary region, the temporal offset of phase information of adjacent pixels is calculated, a spatial distribution field of temporal offset is constructed, and the gradient extremum position of the spatial distribution field is extracted. Connecting the gradient extrema locations forms a boundary tracing path, which is then used as the initial boundary of the artificial forest.

[0008] A double-sided buffer zone is constructed on both sides of the initial boundary. Vegetation features are sampled hierarchically along the normal direction of the initial boundary to obtain a feature profile spanning the initial boundary. The spatial location of the boundary transition zone is determined by the gradient abrupt change features of the feature profile, and the width distribution data of the boundary transition zone is extracted, including: Using the initial boundary as a baseline, an inner buffer zone and an outer buffer zone are constructed by extending a set distance into and out of the artificial forest, forming a double-sided buffer zone. Topographic factors within the double-sided buffer zones are extracted, and the initial boundary is divided into multiple boundary units according to the rate of change of the angle of the tangent direction of the initial boundary. Within each boundary unit, vegetation features are sampled hierarchically at fixed intervals along the normal direction of the initial boundary within the double-sided buffer zone to extract vegetation indices, and then arranged in spatial order to form a feature profile that spans the initial boundary. Calculate the gradient change rate of the feature profile, perform terrain correction on the gradient change rate in combination with terrain factors, identify the extreme points of the corrected gradient change rate as gradient abrupt change features, and extract the gradient peak intensity. Extreme points with values ​​higher than a preset intensity threshold are selected based on the gradient peak intensity, and inner and outer boundary points are determined based on the positions of the extreme points in the inner and outer buffer zones. The width of the transition zone is obtained by calculating the distance between the inner boundary point and the outer boundary point, and the heterogeneity of the transition zone is obtained by calculating the coefficient of variation of the vegetation index within the area between the inner boundary point and the outer boundary point. The transition band width and transition band heterogeneity of each boundary unit are arranged along the initial boundary to form the width distribution data of the boundary transition band.

[0009] The heterogeneity of the transition zone is obtained by calculating the coefficient of variation of vegetation indices within the region between the inner and outer boundary points, including: Extract vegetation indices from all sampling locations within the area between the inner and outer boundary points, calculate the standard deviation and mean of the vegetation indices, and then calculate the basic coefficient of variation by ratio of the standard deviation and mean. Starting from the inner boundary point, vegetation indices are extracted from each sampling location in spatial order along the normal direction to the outer boundary point to form a spatial sequence of vegetation indices. The spatial sequence of vegetation index is segmented and fitted. The positions where the fitting error exceeds the preset error threshold are identified as segment boundaries. The number and distribution interval of segment boundaries are statistically analyzed to construct spatial differentiation feature parameters. The transition zone heterogeneity is obtained by correcting the basic coefficient of variation based on the spatial heterogeneity coefficient.

[0010] The initial boundary is divided into multiple boundary segments. For each boundary segment, the corresponding transition band width value is extracted from the width distribution data. The difference in the transition band width values ​​between adjacent boundary segments is calculated. When the difference exceeds a preset difference threshold, the current boundary segment is marked as an abnormal boundary segment, including: The initial boundary is divided into multiple boundary segments, and the corresponding transition band width value and transition band heterogeneity are extracted from the width distribution data for each boundary segment. Construct a joint feature space for the transition zone width value and the heterogeneity of each boundary segment, calculate the local density of feature points of each boundary segment in the joint feature space, identify boundary segments with local density lower than the global density mean and perform spatial clustering, and extract isolated boundary segments as anomalous response boundary segments. The feature distance between the abnormal response boundary segment and its adjacent boundary segments in the joint feature space is calculated as the degree of difference. A verification path is established inside the plantation along the normal direction of the abnormal response boundary segment whose degree of difference exceeds the preset difference threshold. The vegetation index spatial sequence is extracted on the verification path and trend decomposition is performed to obtain the decay trend component. The slope of the decay trend component is calculated as the decay rate. A two-dimensional discriminant space for decay rate and degree of difference is constructed. Abnormal response boundary segments in the two-dimensional discriminant space are identified where both decay rate and degree of difference deviate from the normal distribution area and are marked as abnormal boundary segments.

[0011] Starting from the anomalous boundary segment, a detection path is established within the plantation along the normal direction of the anomalous boundary segment. Vegetation features are extracted along the detection path, and the correlation between the vegetation features and the width of the transition zone corresponding to the anomalous boundary segment is calculated, including: Starting from the abnormal boundary segment, sampling points are set up inside the artificial forest along the normal direction of the abnormal boundary segment to form a detection path. Multi-temporal vegetation indices are extracted from each sampling point as vegetation characteristics. Using the width of the transition zone corresponding to the abnormal boundary segment as a reference scale, a sliding analysis window with the reference scale as the window length is constructed on the detection path. The standard deviation of vegetation characteristics at each sampling point within the sliding analysis window is calculated, and the change sequence of the standard deviation along the detection path is extracted. Sampling points where the standard deviation changes abruptly are identified from the change sequence as heterogeneity transition points. The temporal stability difference of vegetation characteristics before and after the heterogeneity transition point is calculated, the spatial distance from the heterogeneity transition point to the anomalous boundary segment is calculated, and the distance normalization coefficient is obtained by calculating the ratio of the spatial distance to the width of the transition zone. By multiplying the distance normalization coefficient with the temporal stability difference, the correlation between vegetation characteristics and the transition zone width value corresponding to the anomalous boundary segment is obtained.

[0012] The detection endpoint location is determined based on the decay characteristics of the correlation. The area covered by the detection path from the anomaly boundary segment to the detection endpoint location is designated as a degradation candidate area. The degradation candidate areas are then merged based on connectivity and filtered by morphology. The output of the plantation degradation areas includes: The correlation degree of each sampling point on each detection path is extracted and arranged into a spatial sequence of correlation degree according to the detection path from the anomaly boundary segment to the interior of the plantation; The second-order difference of the correlation spatial sequence is calculated, and the position where the absolute value of the second-order difference is less than the preset difference threshold is identified as the decay stable position. Based on the decay characteristics of the correlation, the sampling point corresponding to the decay stable position is determined as the detection endpoint position. The spatial region covered by the detection path from the abnormal boundary segment to the detection endpoint position is extracted as the degradation candidate region. Calculate the boundary overlap rate and correlation similarity between each degraded candidate region and its adjacent degraded candidate regions. Based on the boundary overlap rate and correlation similarity, merge the degraded candidate regions according to connectivity to form connected degraded regions. The ratio of the spatial variance to the mean of the correlation degree in each connected degraded region is calculated as an internal consistency index. Based on the internal consistency index, the connected degraded regions are morphologically screened, and the plantation degraded regions are output.

[0013] This invention provides a system for extracting degraded areas of plantations based on high-resolution remote sensing, the system comprising: The preprocessing module is used to acquire multi-temporal remote sensing images of the target area and perform preprocessing to extract the initial boundaries of the plantation forest; The transition zone identification module is used to construct a double-sided buffer zone on both sides of the initial boundary, perform layered sampling of vegetation features along the normal direction of the initial boundary to obtain a feature profile that crosses the initial boundary, determine the spatial location of the boundary transition zone through the gradient abrupt change features of the feature profile, and extract the width distribution data of the boundary transition zone. The anomaly marking module is used to divide the initial boundary into multiple boundary segments, extract the corresponding transition band width value from the width distribution data for each boundary segment, calculate the degree of difference in the transition band width value between adjacent boundary segments, and mark the current boundary segment as an abnormal boundary segment when the degree of difference exceeds a preset difference threshold. The correlation calculation module is used to establish a detection path in the artificial forest starting from the abnormal boundary segment and along the normal direction of the abnormal boundary segment. On the detection path, vegetation features are extracted and the correlation between the vegetation features and the width value of the transition zone corresponding to the abnormal boundary segment is calculated. The degradation output module is used to determine the detection endpoint location based on the decay characteristics of the correlation. The area covered by the detection path from the abnormal boundary segment to the detection endpoint location is used as the degradation candidate area. The degradation candidate area is merged based on connectivity and morphological screening to output the degradation area of ​​the plantation forest.

[0014] One technical solution provided in this embodiment of the invention is an electronic device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps in any of the aforementioned methods.

[0015] One technical solution provided in this embodiment of the invention is a computer-readable storage medium storing computer program instructions, which, when executed by a processor, implement the steps in any of the aforementioned methods.

[0016] This invention identifies potential degradation initiation locations by analyzing anomalous changes in the width of boundary transition zones, enabling earlier detection of degradation signs and improving monitoring timeliness. By establishing a correlation between anomalous boundary segments and the degree of degradation within plantations, targeted detection from boundary features to the extent of internal degradation is achieved, avoiding blind searches and improving the targeting and accuracy of degradation area identification. Utilizing correlation decay characteristics to determine the spatial boundary of degradation impacts ensures that degradation area extraction conforms to the spatial diffusion patterns of ecological degradation and has clear quantitative criteria, improving the reliability of the results. Through connectivity merging and morphological filtering, isolated noise points and irregular fragments are effectively removed, resulting in degradation areas with good spatial continuity and morphological integrity, facilitating subsequent degradation causal analysis and remediation planning. The overall method fully utilizes the sensitivity of boundary areas to changes in the external environment, achieving efficient and accurate monitoring of plantation degradation. Attached Figure Description

[0017] Figure 1 A flowchart illustrating the method for extracting degraded areas of plantations based on high-resolution remote sensing, provided in an embodiment of the present invention; Figure 2 This is a schematic diagram of the structure of the plantation degradation area extraction system based on high-resolution remote sensing according to an embodiment of the present invention. Detailed Implementation

[0018] like Figure 1 As shown, Figure 1 A flowchart of a method for extracting degraded areas of plantations based on high-resolution remote sensing provided in an embodiment of the present invention, the method comprising the following steps: Step 101: Acquire multi-temporal remote sensing images of the target area and perform preprocessing to extract the initial boundary of the artificial forest.

[0019] In some embodiments of the present invention, step 101 may specifically include the following sub-steps: Sub-step 1011: Obtain remote sensing image data of the target area at different time phases and perform radiometric correction and geometric registration; extract vegetation index time series data from the registered remote sensing image data. Sub-step 1012 involves seasonally decomposing the vegetation index time series data to obtain trend components and fluctuation components, extracting the spatial decay direction of the trend components and the spatial amplitude distribution of the fluctuation components, calculating the directional consistency between the spatial decay direction and the spatial amplitude distribution, and determining candidate boundary regions based on the spatial abrupt change locations of the directional consistency. Sub-step 1013: Extract phase information of vegetation index time series data within the candidate boundary region, calculate the time series offset of phase information of adjacent pixels, construct the spatial distribution field of the time series offset, and extract the gradient extreme value position of the spatial distribution field. Sub-step 1014: Connect the gradient extreme value positions to form a boundary tracing path, and use the boundary tracing path as the initial boundary of the artificial forest.

[0020] First, acquire multi-temporal remote sensing imagery of the target area. This data can be obtained from various high-resolution remote sensing satellites, including optical and radar remote sensing satellites. The acquired remote sensing images should cover at least one complete vegetation growing season in the target area, with time intervals typically between 15 and 30 days, in order to capture the seasonal changes in vegetation growth.

[0021] Radiometric correction and geometric registration were performed on the acquired multi-temporal remote sensing images. Radiometric correction mainly eliminates the influence of atmospheric scattering and absorption on the image radiometric values, while geometric registration ensures accurate spatial correspondence between images from different temporal phases. Radiometric correction employs an atmospheric radiative transfer model, which considers factors such as solar radiation, atmospheric path radiation, and background radiation, converting the digital values ​​of the remote sensing images into surface reflectance. Geometric registration uses a control point matching method, selecting prominent ground feature points from images of different temporal phases as control points. Affine transformation is then used to achieve spatial correspondence between images, with registration accuracy controlled within 0.5 pixels.

[0022] Vegetation index time-series data were extracted from the registered remote sensing images. The vegetation index was calculated using either the Normalized Difference Vegetation Index (NDVI) or the Enhanced Vegetation Index (EVI). The NDVI calculation formula is: NDVI = (NIR - R) / (NIR + R), where NIR is the near-infrared reflectance and R is the red light reflectance. The EVI further reduces the influence of soil background and atmosphere by incorporating blue light reflectance: EVI = 2.5 × (NIR - R) / (NIR + 6 × R - 7.5 × B + 1), where B is the blue light reflectance.

[0023] Seasonal decomposition was performed on the extracted vegetation index time-series data to separate trend and fluctuation components. The STL decomposition method was used to decompose the vegetation index time-series data into a trend term, a seasonal term, and a residual term. The trend term reflects the long-term trend of vegetation change, the seasonal term reflects the periodic changes in vegetation growth, and the residual term contains random fluctuations. The spatial decay direction of the trend component was extracted using gradient direction analysis to calculate the gradient vector of the trend component in space; the direction of the gradient vector is the spatial decay direction. The spatial amplitude distribution of the fluctuation component was extracted by calculating the standard deviation of the seasonal term; a larger standard deviation indicates more pronounced seasonal fluctuations.

[0024] Vector field correlation analysis was used to calculate the directional consistency between the spatial attenuation direction and the spatial amplitude distribution. For each pixel, the gradient vectors of the trend component and the fluctuation component within a local window centered on that pixel were calculated, and the directional consistency was evaluated by calculating the cosine of the angle between them. Regions with high directional consistency indicate that the vegetation growth trend and seasonal fluctuations are spatially consistent, usually representing the same type of vegetation cover; while locations where directional consistency changes abruptly may be the boundaries of different vegetation types. Candidate boundary regions were determined based on the spatial abrupt changes in directional consistency, and a threshold method was used to identify pixel locations where the directional consistency gradient value exceeded a set threshold.

[0025] Hilbert transform or discrete Fourier transform is used to extract phase information from vegetation index time-series data within candidate boundary regions. Phase information reflects the temporal rhythm of vegetation growth, and different types of vegetation exhibit different phase characteristics. The temporal offset of phase information between adjacent pixels is calculated. For each pair of adjacent pixels, the cross-correlation function of their vegetation index time-series data is calculated, and the time lag when the cross-correlation function reaches its maximum value is the temporal offset.

[0026] Construct a spatial distribution field of temporal offsets, assign the calculated temporal offsets between each pair of adjacent pixels to the corresponding positions to form a spatially continuous distribution field, extract the gradient extremum positions of the spatial distribution field, and use the Canny edge detection algorithm or level set method to identify the pixel positions in the spatial distribution field of temporal offsets where the gradient value reaches a local maximum.

[0027] Connecting the gradient extrema locations forms the boundary tracing path using a shortest path algorithm, such as Dijkstra's algorithm. A connected graph is constructed at the gradient extrema locations, where nodes represent gradient extrema and edge weights are the weighted sum of spatial distances and temporal offset differences between nodes. The path with the minimum weight is found in the connected graph; this path is the boundary tracing path, and it is used as the initial boundary of the artificial forest.

[0028] This invention achieves accurate extraction of plantation boundaries by analyzing the spatiotemporal variation characteristics of vegetation indices in multi-temporal remote sensing images. It overcomes the limitations of traditional methods that rely on single-temporal images and struggle to distinguish vegetation types with similar spectral characteristics, significantly improving the accuracy and reliability of plantation boundary identification. This method does not rely on a large number of labeled samples, reducing labor costs and enhancing its universality and applicability, providing effective technical support for large-scale plantation monitoring and management.

[0029] Step 102: Construct a double-sided buffer zone on both sides of the initial boundary, perform layered sampling of vegetation features along the normal direction of the initial boundary to obtain a feature profile that crosses the initial boundary, determine the spatial location of the boundary transition zone through the gradient abrupt change features of the feature profile, and extract the width distribution data of the boundary transition zone.

[0030] In some embodiments of the present invention, step 102 may specifically include the following sub-steps: Sub-step 1021: Using the initial boundary as a baseline, extend the artificial forest by a set distance to construct an inner buffer zone and an outer buffer zone, forming a double-sided buffer zone. Extract the topographic factors within the double-sided buffer zone and divide the initial boundary into multiple boundary units according to the angular change rate of the tangent direction of the initial boundary. Sub-step 1022: In each boundary unit, vegetation features are sampled in layers along the normal direction of the initial boundary within the double-sided buffer at fixed intervals to extract vegetation indices, and arranged in spatial order to form a feature profile that spans the initial boundary. Sub-step 1023: Calculate the gradient change rate of the feature profile, perform terrain correction on the gradient change rate in combination with terrain factors, identify the extreme points of the corrected gradient change rate as gradient abrupt change features, and extract the gradient peak intensity. Sub-step 1024: Filter out extreme points that are higher than the preset intensity threshold based on the gradient peak intensity, and determine the inner boundary point and the outer boundary point based on the position of the extreme point in the inner buffer and the outer buffer. Sub-step 1025: Calculate the distance between the inner boundary point and the outer boundary point to obtain the width of the transition zone, and calculate the vegetation index variation coefficient within the area between the inner boundary point and the outer boundary point to obtain the heterogeneity of the transition zone. Sub-step 1026: Arrange the transition band width and transition band heterogeneity of each boundary unit along the initial boundary to form the width distribution data of the boundary transition band.

[0031] First, using the initial boundary as a baseline, extend a predetermined distance inwards and outwards from the plantation to construct an inner buffer zone and an outer buffer zone, forming a double-sided buffer zone. The inner buffer zone points inwards from the plantation, and the outer buffer zone points outwards from the plantation. The buffer zone distance is typically set between 30 and 100 meters, adjusted according to the degree of spatial heterogeneity of the study area. The larger the buffer zone, the more complete the transition features can be captured, but it may also introduce more interference information.

[0032] Topographic factors, including slope, aspect, and elevation, were extracted within the bilateral buffer zones. These factors were extracted based on a digital elevation model and calculated using topographic analysis algorithms. Slope represents the degree of surface inclination, aspect represents the direction of surface inclination, and elevation represents the ground elevation. These topographic factors will be used in subsequent analyses to correct for the influence of topography on vegetation characteristics.

[0033] Calculate the tangent angles between adjacent points on the initial boundary, and then calculate the rate of change of the tangent angles along the boundary. Locations with larger rates of change usually correspond to inflection points or regions with greater curvature on the boundary. Divide the boundary into multiple relatively straight boundary cells to facilitate subsequent analysis along the normal direction. The length of the boundary cells is usually controlled between 50 and 200 meters to ensure that the boundary shape within the cell is relatively simple.

[0034] Within each boundary unit, vegetation characteristics are stratified and sampled at fixed intervals along the normal direction of the initial boundary within the bilateral buffer zones. The normal direction is perpendicular to the boundary tangent and points from the inside of the plantation to the outside. The interval for stratified sampling is typically set at 5-10 meters, and vegetation index values ​​are extracted at each sampling point. Sampling begins at the innermost edge of the inner buffer zone, crosses the initial boundary, and extends to the outermost edge of the outer buffer zone, forming a sampling line that crosses the boundary. The vegetation index can be the Normalized Difference Vegetation Index (NDVI) or the Enhanced Vegetation Index (EDI), and the specific calculation method is the same as in step 101.

[0035] The vegetation index values ​​at each point on the sampling line are arranged in spatial order to form a characteristic profile that spans the initial boundary. The characteristic profile reflects the gradual change of vegetation characteristics from the inside to the outside of the plantation. By analyzing the characteristic profile, the location and width of the boundary transition zone can be identified.

[0036] The gradient rate of change of the feature profile is calculated by dividing the difference in vegetation index between adjacent sampling points by the sampling interval distance, which indicates the speed of spatial change of vegetation features.

[0037] Topographic correction is performed on the gradient change rate by incorporating topographic factors to eliminate the influence of topographic changes on the spatial distribution of vegetation characteristics. The topographic correction employs a multiple regression model, using the gradient change rate as the dependent variable and topographic factors such as slope and aspect as independent variables to establish a regression equation. The residuals are used as the corrected gradient change rate. The correction formula is: G cor =G-(a×Slope+b×Aspect+c×Elevation+d), where G cor G represents the corrected gradient change rate, G represents the original gradient change rate, Slope, Aspect, and Elevation represent the slope, aspect, and elevation values ​​at the sampling point, respectively, and a, b, c, and d are regression coefficients.

[0038] The extreme points of the corrected gradient rate of change are identified as gradient abrupt change features, and the gradient peak intensity is extracted. The gradient rate of change typically exhibits fluctuating changes on the feature profile. Extreme points are the locations where the gradient rate of change reaches a local maximum or minimum, corresponding to the regions where the spatial changes in vegetation features are most drastic or most gradual. The gradient peak intensity is the absolute value of the gradient rate of change at the extreme point, reflecting the significance of the changes in vegetation features.

[0039] Extreme points exceeding a preset intensity threshold are selected based on the gradient peak intensity. This threshold is typically set to 1.5-2 times the average gradient rate of change to identify significant gradient abrupt changes. Inner and outer boundary points are determined based on the location of these extreme points within the inner and outer buffer zones, respectively. The inner boundary point is located within the inner buffer zone where a significant gradient abrupt change occurs, and the outer boundary point is located within the outer buffer zone where a significant gradient abrupt change occurs.

[0040] The distance between the inner and outer boundary points is calculated to obtain the width of the transition zone. The width reflects the transition distance from planted forests to non-planted forest areas; a larger width indicates a slower transition, while a smaller width indicates a more distinct boundary. The coefficient of variation of vegetation indices within the area between the inner and outer boundary points is calculated to obtain the heterogeneity of the transition zone. The coefficient of variation, the ratio of the standard deviation to the mean, reflects the degree of variation in vegetation characteristics within the transition zone; a larger coefficient of variation indicates more uneven vegetation characteristics within the transition zone.

[0041] The width and heterogeneity of the transition zones for each boundary unit are arranged along the initial boundary to form the width distribution data of the boundary transition zones. This data will be used in subsequent steps to determine the degree of degradation of the plantation and the spatial distribution of the degraded areas.

[0042] This invention achieves a detailed characterization of the transition zone at the boundary of plantations by constructing a double-sided buffer zone and using feature profile analysis, overcoming the limitations of traditional boundary extraction methods in identifying complex transition zones. The boundary transition zone identification technology based on gradient abrupt change features can not only accurately locate the boundary but also characterize the width variation and heterogeneous distribution of the boundary, providing crucial information for subsequent identification of degraded areas in plantations and improving the accuracy and reliability of degraded area extraction.

[0043] In sub-step 1025, calculating the vegetation index variation coefficient within the region between the inner and outer boundary points to obtain the heterogeneity of the transition zone includes: Extract vegetation indices from all sampling locations within the area between the inner and outer boundary points, calculate the standard deviation and mean of the vegetation indices, and then calculate the basic coefficient of variation by ratio of the standard deviation and mean. Starting from the inner boundary point, vegetation indices are extracted from each sampling location in spatial order along the normal direction to the outer boundary point to form a spatial sequence of vegetation indices. The spatial sequence of vegetation index is segmented and fitted. The positions where the fitting error exceeds the preset error threshold are identified as segment boundaries. The number and distribution interval of segment boundaries are statistically analyzed to construct spatial differentiation feature parameters. The transition zone heterogeneity is obtained by correcting the basic coefficient of variation based on the spatial heterogeneity coefficient.

[0044] Vegetation index values ​​are extracted from all sampling locations within the region between the established inner and outer boundary points. These sampling locations are the sampling points set at fixed intervals along the normal direction. The resulting set of vegetation index values ​​reflects the spatial distribution characteristics of vegetation conditions within the transition zone.

[0045] Statistical analysis was performed on the extracted set of vegetation index values, calculating their standard deviation and mean. The ratio of the calculated standard deviation to the mean of the vegetation index was then used to obtain the basic coefficient of variation. The basic coefficient of variation reflects the overall degree of variation of vegetation indices within the transition zone; a larger coefficient of variation indicates a more uneven vegetation characteristic within the transition zone.

[0046] Starting from the inner boundary point, vegetation indices are extracted spatially along the normal direction to the outer boundary point, forming a spatial sequence of vegetation indices. This spatial sequence preserves the continuous spatial variation information of the vegetation indices, which helps to identify the spatial heterogeneity characteristics within the transition zone. The spatial sequence of vegetation indices can be represented as a one-dimensional array sorted by distance, where each element corresponds to the vegetation index value of a sampling point.

[0047] Piecewise linear fitting or piecewise polynomial fitting methods are used to perform piecewise fitting analysis on the spatial sequence of vegetation indices. The purpose of piecewise fitting is to identify inflection points in the spatial changes of vegetation indices, which often represent significant changes in vegetation type or growth status. The fitting process uses the sliding window method, with the window size typically set to 3-5 sampling points. Linear or quadratic functions are used for fitting within the window.

[0048] For the fitting results within each window, the fitting error is calculated, which is the difference between the actual vegetation index value and the fitted value. The locations where the fitting error exceeds a preset error threshold are identified as segment boundaries. This error threshold is typically set to 1.5-2 times the overall fitting error to ensure that significant change points are identified. These segment boundaries represent abrupt changes in the spatial variation of the vegetation index, reflecting the discontinuity of vegetation characteristics within the transition zone.

[0049] The number and distribution interval of segment boundaries are statistically analyzed to construct spatial differentiation characteristic parameters, which include two parts: segment boundary density and segment boundary distribution uniformity. Segment boundary density represents the number of segment boundaries per unit distance, calculated as the total number of segment boundaries divided by the width of the transition zone; segment boundary distribution uniformity represents the coefficient of variation of the distance between each segment boundary, calculated as the ratio of the standard deviation to the mean of the distance between adjacent segment boundaries.

[0050] The basic coefficient of variation is corrected based on spatial differentiation characteristic parameters to obtain the heterogeneity of the transition zone. The correction formula is: H=CV base ×(1+α×D+β×U), where H represents the heterogeneity of the transition zone, and CV baseThe base coefficient of variation is represented by , D represents the segment boundary density, U represents the evenness of the segment boundary distribution, and α and β are weighting coefficients, typically ranging from 0.2 to 0.5, adjusted according to regional characteristics. When the segment boundary density is high and unevenly distributed, the heterogeneity value of the transition zone is large, indicating complex spatial variations in vegetation characteristics within the transition zone; when the segment boundary density is low and evenly distributed, the heterogeneity value of the transition zone is small, indicating gentle spatial variations in vegetation characteristics within the transition zone.

[0051] The degree of heterogeneity in the transition zone directly reflects the spatial complexity of vegetation characteristics at the boundary of plantations. At the boundary of healthy and stable plantations, the heterogeneity in the transition zone is usually small, and vegetation characteristics show a regular gradual change; while in degraded areas, due to the influence of human disturbance, pests and diseases or other stress factors, the spatial distribution of vegetation characteristics is more complex, showing a higher degree of heterogeneity in the transition zone.

[0052] This invention, by combining the basic coefficient of variation and spatial differentiation characteristic parameters, achieves precise quantification of the spatial heterogeneity of vegetation characteristics in the transition zone of plantation boundaries, effectively capturing the characteristic differences of different types of degraded areas. The transition zone heterogeneity index not only considers the overall statistical characteristics of vegetation indices but also integrates spatial continuity analysis results, accurately reflecting the spatial differentiation characteristics of plantation boundary degradation and laying the foundation for subsequent identification and classification assessment of degraded areas.

[0053] Step 103: Divide the initial boundary into multiple boundary segments, extract the corresponding transition band width value from the width distribution data for each boundary segment, calculate the degree of difference between the transition band width values ​​of adjacent boundary segments, and mark the current boundary segment as an abnormal boundary segment when the degree of difference exceeds the preset difference threshold.

[0054] In some embodiments of the present invention, step 103 may specifically include the following sub-steps: Sub-step 1031: Divide the initial boundary into multiple boundary segments, and extract the corresponding transition band width value and transition band heterogeneity from the width distribution data for each boundary segment; Sub-step 1032: Construct a joint feature space of the transition zone width value and the heterogeneity of the transition zone for each boundary segment, calculate the local density of feature points of each boundary segment in the joint feature space, identify boundary segments with local density lower than the global density mean and perform spatial clustering, and extract isolated boundary segments as abnormal response boundary segments. Sub-step 1033: Calculate the feature distance between the abnormal response boundary segment and its adjacent boundary segments in the joint feature space as the degree of difference. Establish a verification path in the plantation along the normal direction of the abnormal response boundary segment whose degree of difference exceeds the preset difference threshold. Extract the vegetation index spatial sequence on the verification path and perform trend decomposition to obtain the decay trend component. Sub-step 1034: Calculate the slope of the decay trend component as the decay rate, construct a two-dimensional discrimination space for decay rate and degree of difference, and identify abnormal response boundary segments in the two-dimensional discrimination space where both decay rate and degree of difference deviate from the normal distribution area, and mark them as abnormal boundary segments.

[0055] First, the initial boundary is divided into multiple boundary segments. The division method can be based on the geometric features or spatial location of the boundary. The length of each boundary segment is generally controlled between 100 and 300 meters to ensure relatively uniform terrain and vegetation features within each segment. For each segment, the corresponding transition zone width value and transition zone heterogeneity are extracted from the width distribution data obtained in the previous step. These two indicators are important bases for judging the degree of anomaly in the boundary segment.

[0056] For the extracted boundary segment transition zone width and heterogeneity, a joint feature space is constructed for analysis. The joint feature space is a two-dimensional plane, where the horizontal axis represents the transition zone width and the vertical axis represents the transition zone heterogeneity. Each boundary segment corresponds to a feature point in this space, and the coordinates of the feature point are the transition zone width value and the transition zone heterogeneity value of that boundary segment. Boundary segments with abnormal features can be identified through joint feature space analysis.

[0057] Calculate the local density of feature points on each boundary segment in the joint feature space. Local density reflects the degree of clustering of other feature points around a feature point, and can be represented by the number of other feature points within a specific radius around the feature point. The formula for calculating local density is: , where ρ i Let d represent the local density of feature points in the i-th boundary segment, n represent the total number of boundary segments, and d represent the local density of feature points in the i-th boundary segment. ij d represents the Euclidean distance between the i-th and j-th feature points in the feature space. c This represents the cutoff distance, typically taken as 10%-15% of the maximum distance in the feature space. χ is an indicator function, where d... ij <d c The value is 1 if the condition is met, and 0 otherwise.

[0058] Calculate the global density mean, which is the average of the local densities of all feature points on the boundary segments. Identify boundary segments with local densities lower than the global density mean. These boundary segments are relatively isolated in the joint feature space and may represent anomalies. Perform spatial clustering analysis on these boundary segments with low local density using distance-based clustering methods, such as K-means clustering or density clustering algorithms, to group boundary segments with similar features and close spatial locations into one class.

[0059] Isolated boundary segments are extracted from the clustering results; these are those that do not belong to any cluster or belong to a cluster with a small number of nodes. These segments are marked as anomalous response boundary segments. Isolated distribution means that the transition zone characteristics of this boundary segment are significantly different from the surrounding boundary segments, potentially indicating the presence of degradation.

[0060] The feature distance between the abnormal response boundary segment and its adjacent boundary segments in the joint feature space is calculated as the degree of difference. The feature distance is calculated using Euclidean distance. The calculated degree of difference is compared with a preset difference threshold, which is usually set as the 90th percentile of the feature distances between all boundary segments, or adjusted according to the specific application scenario. When the degree of difference exceeds the preset difference threshold, further verification is performed to confirm whether the abnormal response boundary segment does indeed exhibit degradation.

[0061] A verification path is established within the plantation along the normal direction of the boundary segment of the abnormal response where the difference exceeds a preset difference threshold. The verification path extends from the boundary into the plantation, with a length typically of 200-500 meters, to verify whether the boundary anomaly is related to vegetation degradation within the plantation. Vegetation indices are extracted along the verification path at fixed intervals, forming a spatial sequence of vegetation indices, with intervals typically 10-20 meters.

[0062] Trend decomposition was performed on the spatial sequence of vegetation indices along the validation path to extract the decay trend component. Trend decomposition employed either the moving average method or a multinomial fitting method to decompose the vegetation index sequence into trend and fluctuation components. The trend component reflects the overall trend of vegetation index changes along the validation path, while the decay trend component represents the portion of the vegetation index that decreases with increasing distance.

[0063] The slope of the decay trend component is calculated as the decay rate. The decay rate indicates how quickly the vegetation index decreases towards the interior of the plantation along the validation path, and is an important indicator for judging the degree of correlation between boundary anomalies and interior degradation. When the decay rate is negative and has a large absolute value, it indicates that the vegetation condition deteriorates rapidly from the boundary to the interior, and there may be obvious degradation.

[0064] A two-dimensional discriminant space is constructed, representing the decay rate and the degree of difference, with the horizontal axis representing the degree of difference and the vertical axis representing the decay rate. In this discriminant space, feature points in normal boundary segments are typically concentrated in specific regions, forming normal distribution areas. The center and boundary of the normal distribution area can be determined using a multivariate normal distribution model or the contour envelope method. The multivariate normal distribution model is a statistical method that describes the probability distribution of data by calculating the mean vector and covariance matrix of the feature point distribution; it is often used to identify the normal range in the data. The contour envelope method, on the other hand, visually divides normal and abnormal regions by determining the geometric boundaries (such as convex hulls or α-shapes) that contain the vast majority of feature points.

[0065] In a two-dimensional discrimination space, anomalous response boundary segments are identified where both the decay rate and the degree of difference deviate from the normal distribution area. The degree of deviation can be quantified by calculating the Mahalanobis distance from the feature point to the center of the normal distribution area. When the Mahalanobis distance exceeds a preset threshold, the boundary segment is marked as a determined anomalous boundary segment, serving as the boundary portion of the degraded plantation area.

[0066] This invention achieves accurate identification of degraded areas at the boundaries of plantations by combining multidimensional feature space analysis and verification path detection. This method not only considers the anomalous features of the boundary transition zone but also verifies the correlation between anomalies and internal vegetation degradation, effectively avoiding false positives. The combined use of feature space analysis and spatial clustering techniques improves the accuracy of degraded boundary segment identification, while trend decomposition analysis on the verification path enhances the reliability of the discrimination results, providing technical support for the accurate extraction and subsequent monitoring of degraded areas.

[0067] Step 104: Starting from the abnormal boundary segment, establish a detection path inside the artificial forest along the normal direction of the abnormal boundary segment, extract vegetation features on the detection path, and calculate the correlation between the vegetation features and the width value of the transition zone corresponding to the abnormal boundary segment.

[0068] In some embodiments of the present invention, step 104 may specifically include the following sub-steps: Sub-step 1041: Starting from the abnormal boundary segment, set up sampling points in the artificial forest along the normal direction of the abnormal boundary segment to form a detection path, and extract multi-temporal vegetation indices as vegetation features at each sampling point. Sub-step 1042: Using the width of the transition zone corresponding to the abnormal boundary segment as a reference scale, a sliding analysis window with the reference scale as the window length is constructed on the detection path. The standard deviation of vegetation characteristics at each sampling point within the sliding analysis window is calculated. The change sequence of the standard deviation along the detection path is extracted. The sampling points where the standard deviation changes abruptly are identified from the change sequence as heterogeneity transition points. Sub-step 1043: Calculate the temporal stability difference of vegetation characteristics before and after the heterogeneity transition point, calculate the spatial distance from the heterogeneity transition point to the anomalous boundary segment, and calculate the ratio of the spatial distance to the width of the transition zone to obtain the distance normalization coefficient. Sub-step 1044 involves multiplying the distance normalization coefficient with the temporal stability difference to obtain the correlation between vegetation characteristics and the transition zone width value corresponding to the anomalous boundary segment.

[0069] First, a series of sampling points are established within the plantation along its normal direction, starting from the anomalous boundary segment, to form a detection path. The normal direction refers to the direction perpendicular to the boundary segment and pointing inwards from the plantation. The length of the detection path is generally 3-5 times the width of the transition zone corresponding to the anomalous boundary segment to ensure sufficient coverage of the interior area. The sampling points along the detection path are typically spaced 10-20 meters apart to ensure sufficient sampling density to capture spatial variability in vegetation characteristics.

[0070] Multi-temporal vegetation indices are extracted as vegetation features at each sampling point. "Multi-temporal" refers to remote sensing images acquired at different times, typically covering a complete vegetation growth cycle or multiple seasons. For plantation areas, remote sensing images from multiple time phases, such as the early growing season, vigorous growth period, and late growth period, can be selected. At each sampling point, vegetation index values ​​are extracted from the multi-temporal remote sensing images to form a vegetation feature sequence for that sampling point. Vegetation indices can be standardized vegetation indices, enhanced vegetation indices, or other indices suitable for characterizing vegetation status.

[0071] The width of the transition zone corresponding to the anomalous boundary segment is used as a reference scale, reflecting the characteristic spatial range of vegetation change at the boundary. A sliding analysis window, with a length equal to the reference scale, is constructed along the detection path based on this reference scale. The sliding analysis window is a region that can move along the detection path, used for local analysis of the spatial variability of vegetation characteristics. The window starts at the beginning of the detection path and moves forward one sampling point distance at a time until it covers the entire detection path.

[0072] At each sliding analysis window position, the standard deviation of vegetation characteristics at each sampling point within the window is calculated. The variation sequence of the standard deviation along the detection path is extracted; this sequence reflects the changing trend of spatial heterogeneity of vegetation characteristics. Within normal areas, the spatial heterogeneity of vegetation characteristics should be relatively stable; however, at the transition between degraded and healthy areas, heterogeneity often changes significantly. Sampling points where the standard deviation undergoes abrupt changes are identified from the variation sequence and marked as heterogeneity transition points. The abrupt change criterion can be set as the rate of change of the standard deviation exceeding 2 or 3 times the average rate of change of the sequence.

[0073] The temporal stability difference of vegetation characteristics before and after the heterogeneity transition point is calculated. Temporal stability refers to the degree of fluctuation of vegetation indices over time, which can be characterized by the coefficient of variation. The temporal coefficient of variation is calculated separately for the regions before and after the heterogeneity transition point. The region before the transition point refers to the area from the transition point to the starting point of the detection path, and the region after the transition point refers to the area from the transition point to the ending point of the detection path. The difference in temporal stability is calculated as the absolute difference in the temporal coefficients of variation between the two regions.

[0074] The spatial distance from the heterogeneity transition point to the anomalous boundary segment is calculated. This distance reflects the extent to which degradation extends inward from the boundary. The spatial distance can be calculated by measuring the distance between the transition point's location on the probe path and the starting point of the probe path, taking into account actual geographic coordinates for precise measurement. The ratio of the spatial distance to the transition zone width is calculated to obtain the distance normalization coefficient. The distance normalization coefficient represents the proportional relationship between the degradation extension distance and the boundary transition zone width, and can be used to assess the spatial extent of degradation.

[0075] The correlation between vegetation characteristics and the transition zone width corresponding to the anomalous boundary segment is obtained by multiplying the distance normalization coefficient with the temporal stability difference. The larger the correlation value, the stronger the correlation between the anomalous features of the boundary segment and the vegetation features of the interior region, meaning that the boundary anomaly is likely directly related to the degradation of the interior region.

[0076] Association degree grading can be achieved by calculating the distribution characteristics of association degree values ​​and further classifying them into three levels: high association, medium association, and low association (high association ≥ 75th percentile, medium association between 25th and 75th percentiles, and low association ≤ 25th percentile) using the percentile method. High association indicates that boundary anomalies are closely related to internal degradation and may be different manifestations of the same degradation process; medium association indicates that there is some association between the two, but it may be influenced by other factors; low association indicates that boundary anomalies may be caused by other reasons and are not closely related to internal vegetation degradation.

[0077] Correlation analysis can identify boundary anomalies that are indeed related to vegetation degradation, improving the accuracy of degraded area extraction. Highly correlated anomalous boundary segments and their corresponding internal regions can be considered as key degraded areas, providing a basis for subsequent precise delineation of degraded area boundaries and assessment of degradation levels.

[0078] This invention achieves accurate identification and boundary determination of degraded areas in plantations by establishing detection paths and analyzing the correlation between vegetation and boundary features. Multi-temporal vegetation feature analysis enhances the spatiotemporal dimension of degradation identification, while the sliding window technique captures spatial heterogeneity changes. Heterogeneity transition points reflect the location of degradation boundaries, and the correlation index constructed by temporal stability differences and distance normalization effectively determines the degree of correlation between boundary anomalies and internal degradation. This multi-dimensional degradation area identification method overcomes the limitations of traditional single-indicator identification and can accurately extract morphologically complex and varying degrees of plantation degradation areas.

[0079] Step 105: Determine the detection endpoint location based on the decay characteristics of the correlation. Take the area covered by the detection path from the abnormal boundary segment to the detection endpoint location as the degradation candidate area. Perform connectivity merging and morphological screening on the degradation candidate area to output the degradation area of ​​the plantation forest.

[0080] In some embodiments of the present invention, step 105 may specifically include the following sub-steps: Sub-step 1051: Extract the correlation degree of each sampling point on each detection path, and arrange them in a spatial sequence of correlation degree from the anomaly boundary segment to the interior of the artificial forest according to the detection path; Sub-step 1052: Calculate the second difference of the correlation spatial sequence, identify the position where the absolute value of the second difference is less than the preset difference threshold as the decay stable position, determine the sampling point corresponding to the decay stable position as the detection endpoint position according to the decay characteristics of the correlation, and extract the spatial area covered by the detection path from the abnormal boundary segment to the detection endpoint position as the degradation candidate area. Sub-step 1053: Calculate the boundary overlap rate and correlation similarity between each degraded candidate region and its adjacent degraded candidate regions, and merge the degraded candidate regions based on the boundary overlap rate and correlation similarity to form a connected degraded region; Sub-step 1054: Calculate the ratio of the spatial variance to the mean of the correlation degree in each connected degraded region as an internal consistency index. Based on the internal consistency index, perform morphological screening on the connected degraded regions and output the degraded plantation forest regions.

[0081] The correlation degree obtained in the previous step is an important indicator for measuring the correlation between boundary anomalies and internal degradation. It is necessary to extract the correlation degree values ​​of each sampling point along each detection path. These correlation degree values ​​are arranged according to the detection path from the anomaly boundary segment towards the interior of the plantation, forming a spatial sequence of correlation degree. The spatial sequence of correlation degree reflects the changing trend of degradation correlation from the boundary to the interior, which is crucial for determining the boundary of the degradation area.

[0082] The second difference of the spatial sequence of correlation is calculated to detect inflection points in the trend of correlation changes. The second difference represents the change in the rate of change of correlation, and can accurately identify the position where the correlation changes from rapid change to stability. The formula for calculating the second difference is as follows: Δ 2 C y =(C y+2 -C y+1 )-(C y+1 -C y )=C y+2 -2C y+1 +C y， Where, Δ 2 C y C represents the second difference value of the y-th point in the correlation space sequence. y C y+1 and C y+2 These represent the correlation values ​​of three consecutive sampling points in the spatial sequence.

[0083] By calculating second-order differences, locations where the absolute value of the second-order difference is less than a preset difference threshold are identified as stable decay locations. The preset difference threshold is typically set to 20% to 30% of the average absolute value of the second-order differences, or adjusted according to the specific application scenario. A stable decay location indicates a point where the correlation change tends to level off, meaning the degradation effect has stabilized at this location. Based on the decay characteristics of the correlation, the sampling point corresponding to the stable decay location is determined as the detection endpoint. The detection endpoint is usually the boundary point where the degraded region transitions to the healthy region.

[0084] The spatial region covered by the detection path from the anomaly boundary segment to the detection endpoint is extracted and used as a degradation candidate region. This candidate region represents areas initially identified as potentially degraded and requires further verification and optimization. If multiple detection paths are adjacent or intersect, their corresponding degradation candidate regions may partially overlap; in this case, their connectivity needs to be considered.

[0085] Calculate the boundary overlap rate and correlation similarity between each degraded candidate region and its adjacent degraded candidate regions. The boundary overlap rate represents the degree to which two degraded candidate regions share a boundary, and is calculated as the ratio of the length of the shared boundary to the total length of the two region boundaries. The correlation similarity represents the similarity of the correlation distribution within two degraded candidate regions, and can be characterized by calculating the difference in the mean correlation values ​​between the two regions. The correlation similarity calculation formula is: , of which S te M represents the correlation similarity between the t-th and e-th degenerate candidate regions. t and M e These represent the average correlation between the two candidate regions, ranging from 0 to 1. A higher value indicates a higher similarity.

[0086] Degraded candidate regions are merged based on boundary overlap rate and correlation similarity to form connected degraded regions. When the boundary overlap rate exceeds 30% and the correlation similarity is higher than 0.8, two degraded candidate regions can be merged into one connected region. Connectivity merging can eliminate redundant degraded candidate regions, forming more complete and coherent degraded regions. The merging process is iterative until no candidate regions that meet the conditions can be merged.

[0087] The ratio of the spatial variance to the mean of the connectivity within each connected degraded region is calculated as an internal consistency index. Internal consistency reflects the uniformity of the connectivity distribution within a degraded region and is an important parameter for measuring the reliability of the degraded region. Spatial variance calculates the dispersion of connectivity in spatial distribution, while the mean represents the overall connectivity level. A smaller ratio indicates higher internal consistency and greater reliability of the degraded region.

[0088] Morphological screening of connected degraded regions is performed based on internal consistency indices to output the final degraded plantation areas. Screening criteria can include an internal consistency index below a preset threshold, typically between 0.3 and 0.5. Further screening can be conducted based on morphological features such as area and shape complexity to remove excessively small or irregularly shaped areas, which may be due to noise or misjudgment. Morphological screening can employ mathematical morphological operations, such as opening and closing operations, to smooth boundaries, fill small holes, and remove isolated regions.

[0089] The output degraded plantation areas are the result of multi-level verification and optimization, ensuring high reliability. These areas can be saved in vector form, containing boundary coordinates and attribute information such as degradation degree and correlation distribution, facilitating subsequent analysis and management. The extracted degraded areas can be overlaid on the original remote sensing image for visualization, intuitively presenting the spatial distribution characteristics of plantation degradation.

[0090] The spatial distribution characteristics of degraded areas can be used to analyze the causes and trends of degradation. For example, the shape, size, and spatial distribution patterns of degraded areas may be related to specific environmental factors, pest and disease transmission methods, or human disturbance. This information is of great reference value for developing targeted protection and restoration measures.

[0091] This invention achieves accurate extraction and boundary delineation of degraded areas through correlation decay analysis and region merging screening. This method combines spatial sequence analysis, second-order difference detection, connectivity merging, and morphological screening techniques, effectively solving the identification challenges of blurred boundaries and irregular shapes of degraded areas in complex landscape backgrounds using traditional methods. By determining the degradation range through spatial variation characteristics of correlation, the spatial accuracy and reliability of plantation degradation monitoring are significantly improved.

[0092] like Figure 2 As shown, Figure 2 This is a schematic diagram of a system for extracting degraded areas of plantations based on high-resolution remote sensing, provided in an embodiment of the present invention. The system includes: The preprocessing module 201 is used to acquire multi-temporal remote sensing images of the target area and perform preprocessing to extract the initial boundary of the artificial forest; The transition zone identification module 202 is used to construct a double-sided buffer zone on both sides of the initial boundary, perform layered sampling of vegetation features along the normal direction of the initial boundary to obtain a feature profile that crosses the initial boundary, determine the spatial location of the boundary transition zone through the gradient abrupt change features of the feature profile, and extract the width distribution data of the boundary transition zone. The anomaly marking module 203 is used to divide the initial boundary into multiple boundary segments, extract the corresponding transition band width value from the width distribution data for each boundary segment, calculate the degree of difference between the transition band width values ​​of adjacent boundary segments, and mark the current boundary segment as an abnormal boundary segment when the degree of difference exceeds a preset difference threshold. The correlation calculation module 204 is used to establish a detection path in the artificial forest starting from the abnormal boundary segment and along the normal direction of the abnormal boundary segment, extract vegetation features on the detection path and calculate the correlation between the vegetation features and the width value of the transition zone corresponding to the abnormal boundary segment. The degradation output module 205 is used to determine the detection endpoint location based on the decay characteristics of the correlation. The area covered by the detection path from the abnormal boundary segment to the detection endpoint location is taken as the degradation candidate area. The degradation candidate area is merged for connectivity and filtered for morphology, and the degradation area of ​​the plantation forest is output.

[0093] One technical solution provided in this embodiment of the invention is an electronic device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps in any of the aforementioned methods.

[0094] One technical solution provided in this embodiment of the invention is a computer-readable storage medium storing a computer program, wherein the processor executes the computer program to implement the steps in any of the aforementioned methods.

[0095] The specific embodiments described above are preferred embodiments of the present invention and are not intended to limit the specific scope of the present invention. The scope of the present invention includes, but is not limited to, these specific embodiments. All equivalent changes made in accordance with the shape and structure of the present invention are within the protection scope of the present invention.

Claims

1. A method for extracting degraded areas of plantations based on high-resolution remote sensing, characterized in that, Includes the following steps: Acquire multi-temporal remote sensing images of the target area and perform preprocessing to extract the initial boundaries of the plantation forest; A double-sided buffer zone is constructed on both sides of the initial boundary. Vegetation features are sampled in layers along the normal direction of the initial boundary to obtain a feature profile that crosses the initial boundary. The spatial location of the boundary transition zone is determined by the gradient abrupt change features of the feature profile, and the width distribution data of the boundary transition zone is extracted. The initial boundary is divided into multiple boundary segments. For each boundary segment, the corresponding transition band width value is extracted from the width distribution data. The degree of difference between the transition band width values ​​of adjacent boundary segments is calculated. When the degree of difference exceeds the preset difference threshold, the current boundary segment is marked as an abnormal boundary segment. Starting from the abnormal boundary segment, a detection path is established inside the plantation along the normal direction of the abnormal boundary segment. Vegetation features are extracted on the detection path and the correlation between the vegetation features and the width of the transition zone corresponding to the abnormal boundary segment is calculated. The detection endpoint is determined based on the decay characteristics of the correlation. The area covered by the detection path from the abnormal boundary segment to the detection endpoint is taken as the degradation candidate area. The degradation candidate area is merged based on connectivity and morphological screening to output the degradation area of ​​the plantation forest.

2. The method according to claim 1, characterized in that, Acquire multi-temporal remote sensing images of the target area and perform preprocessing to extract the initial boundaries of the plantation, including: Remote sensing image data of the target area at different time phases are acquired and radiometrically corrected and geometrically registered. Vegetation index time series data are extracted from the registered remote sensing image data. Seasonal decomposition of vegetation index time series data yields trend and fluctuation components. The spatial decay direction of the trend component and the spatial amplitude distribution of the fluctuation component are extracted. The directional consistency between the spatial decay direction and the spatial amplitude distribution is calculated. Candidate boundary regions are determined based on the spatial abrupt change locations of the directional consistency. Phase information of vegetation index time series data is extracted within the candidate boundary region, the temporal offset of phase information of adjacent pixels is calculated, a spatial distribution field of temporal offset is constructed, and the gradient extremum position of the spatial distribution field is extracted. Connecting the gradient extrema locations forms a boundary tracing path, which is then used as the initial boundary of the artificial forest.

3. The method according to claim 1, characterized in that, A double-sided buffer zone is constructed on both sides of the initial boundary. Vegetation features are sampled hierarchically along the normal direction of the initial boundary to obtain a feature profile spanning the initial boundary. The spatial location of the boundary transition zone is determined by the gradient abrupt change features of the feature profile, and the width distribution data of the boundary transition zone is extracted, including: Using the initial boundary as a baseline, an inner buffer zone and an outer buffer zone are constructed by extending a set distance into and out of the artificial forest, forming a double-sided buffer zone. Topographic factors within the double-sided buffer zones are extracted, and the initial boundary is divided into multiple boundary units according to the rate of change of the angle of the tangent direction of the initial boundary. Within each boundary unit, vegetation features are sampled hierarchically at fixed intervals along the normal direction of the initial boundary within the double-sided buffer zone to extract vegetation indices, and then arranged in spatial order to form a feature profile that spans the initial boundary. Calculate the gradient change rate of the feature profile, perform terrain correction on the gradient change rate in combination with terrain factors, identify the extreme points of the corrected gradient change rate as gradient abrupt change features, and extract the gradient peak intensity. Extreme points with values ​​higher than a preset intensity threshold are selected based on the gradient peak intensity, and inner and outer boundary points are determined based on the positions of the extreme points in the inner and outer buffer zones. The width of the transition zone is obtained by calculating the distance between the inner boundary point and the outer boundary point, and the heterogeneity of the transition zone is obtained by calculating the coefficient of variation of the vegetation index within the area between the inner boundary point and the outer boundary point. The transition band width and transition band heterogeneity of each boundary unit are arranged along the initial boundary to form the width distribution data of the boundary transition band.

4. The method according to claim 3, characterized in that, The heterogeneity of the transition zone is obtained by calculating the coefficient of variation of vegetation indices within the region between the inner and outer boundary points, including: Extract vegetation indices from all sampling locations within the area between the inner and outer boundary points, calculate the standard deviation and mean of the vegetation indices, and then calculate the basic coefficient of variation by ratio of the standard deviation and mean. Starting from the inner boundary point, vegetation indices are extracted from each sampling location in spatial order along the normal direction to the outer boundary point to form a spatial sequence of vegetation indices. The spatial sequence of vegetation index is segmented and fitted. The positions where the fitting error exceeds the preset error threshold are identified as segment boundaries. The number and distribution interval of segment boundaries are statistically analyzed to construct spatial differentiation feature parameters. The transition zone heterogeneity is obtained by correcting the basic coefficient of variation based on the spatial heterogeneity coefficient.

5. The method according to claim 1, characterized in that, The initial boundary is divided into multiple boundary segments. For each boundary segment, the corresponding transition band width value is extracted from the width distribution data. The difference in the transition band width values ​​between adjacent boundary segments is calculated. When the difference exceeds a preset difference threshold, the current boundary segment is marked as an abnormal boundary segment, including: The initial boundary is divided into multiple boundary segments, and the corresponding transition band width value and transition band heterogeneity are extracted from the width distribution data for each boundary segment. Construct a joint feature space for the transition zone width value and the heterogeneity of each boundary segment, calculate the local density of feature points of each boundary segment in the joint feature space, identify boundary segments with local density lower than the global density mean and perform spatial clustering, and extract isolated boundary segments as anomalous response boundary segments. The feature distance between the abnormal response boundary segment and its adjacent boundary segments in the joint feature space is calculated as the degree of difference. A verification path is established inside the plantation along the normal direction of the abnormal response boundary segment whose degree of difference exceeds the preset difference threshold. The vegetation index spatial sequence is extracted on the verification path and trend decomposition is performed to obtain the decay trend component. The slope of the decay trend component is calculated as the decay rate. A two-dimensional discriminant space for decay rate and degree of difference is constructed. Abnormal response boundary segments in the two-dimensional discriminant space are identified where both decay rate and degree of difference deviate from the normal distribution area and are marked as abnormal boundary segments.

6. The method according to claim 1, characterized in that, Starting from the anomalous boundary segment, a detection path is established within the plantation along the normal direction of the anomalous boundary segment. Vegetation features are extracted along the detection path, and the correlation between the vegetation features and the width of the transition zone corresponding to the anomalous boundary segment is calculated, including: Starting from the abnormal boundary segment, sampling points are set up inside the artificial forest along the normal direction of the abnormal boundary segment to form a detection path. Multi-temporal vegetation indices are extracted from each sampling point as vegetation characteristics. Using the width of the transition zone corresponding to the abnormal boundary segment as a reference scale, a sliding analysis window with the reference scale as the window length is constructed on the detection path. The standard deviation of vegetation characteristics at each sampling point within the sliding analysis window is calculated, and the change sequence of the standard deviation along the detection path is extracted. Sampling points where the standard deviation changes abruptly are identified from the change sequence as heterogeneity transition points. The temporal stability difference of vegetation characteristics before and after the heterogeneity transition point is calculated, the spatial distance from the heterogeneity transition point to the anomalous boundary segment is calculated, and the distance normalization coefficient is obtained by calculating the ratio of the spatial distance to the width of the transition zone. By multiplying the distance normalization coefficient with the temporal stability difference, the correlation between vegetation characteristics and the transition zone width value corresponding to the anomalous boundary segment is obtained.

7. The method according to claim 1, characterized in that, The detection endpoint location is determined based on the decay characteristics of the correlation. The area covered by the detection path from the anomaly boundary segment to the detection endpoint location is designated as a degradation candidate area. The degradation candidate areas are then merged based on connectivity and filtered by morphology. The output of the plantation degradation areas includes: The correlation degree of each sampling point on each detection path is extracted and arranged into a spatial sequence of correlation degree according to the detection path from the anomaly boundary segment to the interior of the plantation; The second-order difference of the correlation spatial sequence is calculated, and the position where the absolute value of the second-order difference is less than the preset difference threshold is identified as the decay stable position. Based on the decay characteristics of the correlation, the sampling point corresponding to the decay stable position is determined as the detection endpoint position. The spatial region covered by the detection path from the abnormal boundary segment to the detection endpoint position is extracted as the degradation candidate region. Calculate the boundary overlap rate and correlation similarity between each degraded candidate region and its adjacent degraded candidate regions. Based on the boundary overlap rate and correlation similarity, merge the degraded candidate regions according to connectivity to form connected degraded regions. The ratio of the spatial variance to the mean of the correlation degree in each connected degraded region is calculated as an internal consistency index. Based on the internal consistency index, the connected degraded regions are morphologically screened, and the plantation degraded regions are output.

8. A system for extracting degraded areas of plantations based on high-resolution remote sensing, used to implement the method described in any one of claims 1-7, characterized in that, The system includes: The preprocessing module is used to acquire multi-temporal remote sensing images of the target area and perform preprocessing to extract the initial boundaries of the plantation forest; The transition zone identification module is used to construct a double-sided buffer zone on both sides of the initial boundary, perform layered sampling of vegetation features along the normal direction of the initial boundary to obtain a feature profile that crosses the initial boundary, determine the spatial location of the boundary transition zone through the gradient abrupt change features of the feature profile, and extract the width distribution data of the boundary transition zone. The anomaly marking module is used to divide the initial boundary into multiple boundary segments, extract the corresponding transition band width value from the width distribution data for each boundary segment, calculate the degree of difference in the transition band width value between adjacent boundary segments, and mark the current boundary segment as an abnormal boundary segment when the degree of difference exceeds a preset difference threshold. The correlation calculation module is used to establish a detection path in the artificial forest starting from the abnormal boundary segment and along the normal direction of the abnormal boundary segment. On the detection path, vegetation features are extracted and the correlation between the vegetation features and the width value of the transition zone corresponding to the abnormal boundary segment is calculated. The degradation output module is used to determine the detection endpoint location based on the decay characteristics of the correlation. The area covered by the detection path from the abnormal boundary segment to the detection endpoint location is used as the degradation candidate area. The degradation candidate area is merged based on connectivity and morphological screening to output the degradation area of ​​the plantation forest.

9. An electronic device, characterized in that, include: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the computer program, implements the steps of the method as described in any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores computer program instructions that, when executed by a processor, implement the steps of the method as described in any one of claims 1 to 7.

Images