A forestry resource dynamic monitoring system and method based on remote sensing images

Abstract

This invention relates to the field of forestry resource monitoring technology, and discloses a dynamic monitoring system and method for forestry resources based on remote sensing imagery. The system includes: generating temporal remote sensing reflectance data within the same grid; generating a continuous canopy occupancy map; generating an infill set; determining the intrinsic neck scale; calculating the topological spectrum intensity of the continuous canopy infill; calculating the interlayer decoupling topological index; generating monitoring judgment results; and outputting change patches. This invention constructs an adaptive multi-scale erosion sequence based on the intrinsic neck scale, combines Eulerian features to quantify the infill topological fragmentation process, and integrates the degree of infill fragmentation with the stability state of the outer envelope through the interlayer decoupling topological index. This achieves accurate determination of change units and precise spatial patch location, avoiding the omissions or misjudgments of deep forest understory structure changes by traditional techniques, ensuring that the monitoring results accurately reflect the actual mechanisms of dynamic changes in forest resources.

Description

Technical Field

[0001] This invention relates to the field of forestry resource monitoring technology, and more specifically, to a dynamic monitoring system and method for forestry resources based on remote sensing imagery. Background Technology

[0002] Dynamic monitoring of forestry resources is the core support for ecological protection and management. Mixed uneven-aged forests, as an important forest type, generally adopt a continuous cover management model. Their upper canopy often appears continuous and intact in remote sensing images, without obvious large-scale breaks or loss of green. The dynamic changes of understory resources in this type of forest are complex, including selective felling of target trees, expansion of the natural regeneration layer, and reconstruction of microcavities, which directly affect the forest's ecological function and management benefits, and urgently require the support of accurate dynamic monitoring technology.

[0003] Current mainstream forestry remote sensing monitoring technologies mostly focus on the overall characteristics of the canopy, analyzing indicators such as canopy coverage, average reflectance, vegetation greenness, and changes in large patch boundaries to determine the status of forest resources. These methods rely on significant changes in the external morphology or overall spectral characteristics of the canopy, employing techniques such as simple spectral difference, fixed threshold segmentation, and comparison after conventional classification. The core logic is to capture abrupt changes in the overall area or spectral characteristics of the canopy.

[0004] However, in continuously covered mixed uneven-aged forests, substantial abrupt changes in understory resources often precede changes in the external morphology of the canopy. This manifests as a relatively stable outer canopy envelope, but with internal microcavities undergoing topological fragmentation, contraction, and splitting. Traditional monitoring techniques, failing to address this deep structural change, cannot distinguish between a stable outer envelope with fragmented internal pores and a lack of substantial change. This leads to missed or misjudged dynamic changes in understory resources, making it difficult to meet the practical needs of refined dynamic monitoring of such forests. Summary of the Invention

[0005] This invention provides a forestry resource dynamic monitoring system and method based on remote sensing imagery, which solves the technical problems mentioned in the background.

[0006] This invention provides a forestry resource dynamic monitoring system based on remote sensing imagery, comprising:

[0007] The first module normalizes multispectral images to generate temporal remote sensing reflectance data in the same grid.

[0008] The second module uses the same grid of temporal remote sensing reflectance data to construct the canopy occupancy field and binarizes it to generate a continuous canopy occupancy map.

[0009] The third module performs hole filling on the continuous canopy occupancy map to extract the continuous canopy outer envelope, calculates the area of ​​the outer envelope, and generates an inner hole set by subtracting the continuous canopy occupancy map from the continuous canopy outer envelope.

[0010] The fourth module extracts the skeleton based on the boundary distance of the inner hole set, and determines the intrinsic neck size based on the local hole width of the skeleton;

[0011] The fifth module divides the discrete erosion scale sequence according to the intrinsic neck scale, generates a residual erosion set sequence by eroding the pore set step by step, and calculates the pore topological spectrum intensity of the continuous canopy by combining the Eulerian features of the residual pore set sequence.

[0012] The sixth module constructs an outer envelope stable modulation term based on the temporal variation of the outer envelope area, and fuses the topological spectrum intensity of the continuous canopy pores with the outer envelope stable modulation term to obtain the interlayer decoupling topological index;

[0013] The seventh module calculates the temporal dynamic change based on the interlayer decoupling topology index, and performs threshold segmentation on the temporal dynamic change to generate monitoring and judgment results;

[0014] The eighth module locates the principal scale position with the largest increase in Eulerian features in the discrete erosion scale sequence based on the monitoring and judgment results. It then compares the differences in the residual pore set before and after the principal scale position to determine the main topological fragmentation zone and outputs the change patch by connecting the domains.

[0015] This invention provides a method for dynamic monitoring of forestry resources based on remote sensing imagery, comprising the following steps:

[0016] Step S1: Normalize the multispectral image to generate temporal remote sensing reflectance data for the same grid.

[0017] Step S2: Construct a canopy occupancy field using the same grid temporal remote sensing reflectance data, and binarize it to generate a continuous canopy occupancy map;

[0018] Step S3: Perform hole filling on the continuous canopy occupancy map to extract the continuous canopy outer envelope, calculate the area of ​​the outer envelope, and generate the inner hole set by subtracting the continuous canopy occupancy map from the continuous canopy outer envelope;

[0019] Step S4: Extract the skeleton based on the boundary distance of the inner pore set, and determine the intrinsic neck size based on the local pore width of the skeleton;

[0020] Step S5: Divide the discrete erosion scale sequence according to the intrinsic neck scale, erode the pore set step by step to generate the residual pore set sequence, and calculate the pore topological spectrum intensity of the continuous canopy by combining the Euler features of the residual pore set sequence.

[0021] Step S6: Construct an outer envelope stable modulation term based on the temporal change of the outer envelope area, and fuse the topological spectrum intensity of the continuous canopy pores with the outer envelope stable modulation term to obtain the interlayer decoupling topological index;

[0022] Step S7: Calculate the temporal dynamic change based on the interlayer decoupling topology index, and perform threshold segmentation on the temporal dynamic change to generate monitoring and judgment results;

[0023] Step S8: Based on the monitoring and judgment results, locate the main scale position with the largest increase in Euler feature in the discrete erosion scale sequence, compare the difference in the residual pore set before and after the main scale position to determine the main topological fragmentation zone, and output the change patch by connecting the domains.

[0024] The beneficial effects of this invention are as follows: By shifting the monitoring object from the overall characteristics of the canopy to the set of pores within the continuous canopy, this invention can capture substantial abrupt changes in understory resources under the stable state of the canopy envelope in continuously covered mixed uneven-aged forests. Based on an intrinsic neck scale, an adaptive multi-scale erosion sequence is constructed. The pore topological fragmentation process is quantified using Eulerian features. The degree of pore fragmentation and the stable state of the outer envelope are integrated through interlayer decoupled topological indices, ultimately achieving the determination of change units and precise spatial location of patches. The entire technical process is adapted to the business logic of forestry unit-based monitoring, avoiding the omissions or misjudgments of deep understory structural changes that are common with traditional techniques. This ensures that the monitoring results not only conform to the actual mechanisms of dynamic changes in forest resources but also directly serve the actual needs of forestry management, improving the targeting and adaptability of dynamic monitoring of continuously covered mixed uneven-aged forests. Attached Figure Description

[0025] Figure 1 This is a flowchart of a method for dynamic monitoring of forestry resources based on remote sensing imagery according to the present invention;

[0026] Figure 2 This is a schematic diagram of the multi-temporal dynamic monitoring of the present invention;

[0027] Figure 3 This is a schematic diagram of the comparative experiment of the present invention. Detailed Implementation

[0028] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, features described in some examples may be combined in other examples.

[0029] It should be noted that, unless otherwise defined, the technical or scientific terms used in one or more embodiments of the present invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in one or more embodiments of the present invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" indicate that the element or object preceding the term encompasses the elements or objects listed following the term and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0030] like Figures 1-3 As shown, a forestry resource dynamic monitoring system based on remote sensing imagery includes:

[0031] The first module normalizes multispectral images to generate temporal remote sensing reflectance data in the same grid.

[0032] The second module uses the same grid of temporal remote sensing reflectance data to construct the canopy occupancy field and binarizes it to generate a continuous canopy occupancy map.

[0033] The third module performs hole filling on the continuous canopy occupancy map to extract the continuous canopy outer envelope, calculates the area of ​​the outer envelope, and generates an inner hole set by subtracting the continuous canopy occupancy map from the continuous canopy outer envelope.

[0034] The fourth module extracts the skeleton based on the boundary distance of the inner hole set, and determines the intrinsic neck size based on the local hole width of the skeleton;

[0035] The fifth module divides the discrete erosion scale sequence according to the intrinsic neck scale, generates a residual erosion set sequence by eroding the pore set step by step, and calculates the pore topological spectrum intensity of the continuous canopy by combining the Eulerian features of the residual pore set sequence.

[0036] The sixth module constructs an outer envelope stable modulation term based on the temporal variation of the outer envelope area, and fuses the topological spectrum intensity of the continuous canopy pores with the outer envelope stable modulation term to obtain the interlayer decoupling topological index;

[0037] The seventh module calculates the temporal dynamic change based on the interlayer decoupling topology index, and performs threshold segmentation on the temporal dynamic change to generate monitoring and judgment results;

[0038] The eighth module locates the principal scale position with the largest increase in Eulerian features in the discrete erosion scale sequence based on the monitoring and judgment results. It then compares the differences in the residual pore set before and after the principal scale position to determine the main topological fragmentation zone and outputs the change patch by connecting the domains.

[0039] In one embodiment of the present invention, the method of normalizing multispectral images to generate temporal remote sensing reflectance data within the same grid includes:

[0040] In each phase Below, the multispectral image is positioned at any pixel location within a predetermined grid. Multispectral reflectance data are generated at the location The formula is:

[0041] ;

[0042] in For phase Multispectral reflectance data, Red light reflectance, For green light reflectance, Red-edge reflectivity, Near-infrared reflectance;

[0043] For each time phase With each band At all cell positions in the same grid The median of the multispectral reflectance data is taken to obtain the median operator result. ,in ;

[0044] For each time phase With each band Constructing the absolute deviation and at all cell positions in the same grid The median absolute deviation is taken from the absolute deviation to obtain the result of the median absolute deviation operator. ;

[0045] For each time phase Each band Each pixel position Calculate the normalized reflectance using the following formula. :

[0046] ;

[0047] in For normalized reflectivity, It is a positive constant;

[0048] For each time phase With each cell position Normalized reflectance is assembled into temporal remote sensing reflectance data of the same grid. The formula is: ;

[0049] in This is remote sensing reflectance data from the same grid at different times.

[0050] It should be noted that multispectral imagery refers to forest canopy spectral imaging data acquired by remote sensing equipment equipped with multispectral sensors. It reflects the radiation and reflectance characteristics of the forest canopy in different spectral bands and can be acquired through satellite multispectral remote sensing systems, UAV multispectral imagers, and airborne multispectral sensors. A grid is a unified pixel raster spatial reference set for multispectral images from different time periods, reflecting the geometric registration standard of multi-source remote sensing images. The preferred value is a square raster with a pixel resolution of 1 to 10 meters. This resolution range balances spatial detail in forestry resource monitoring with data processing efficiency, and is suitable for monitoring the microstructure of mixed-age forest canopies. Pixel location is the geospatial coordinate point corresponding to each pixel in the grid, which can be determined using remote sensing image georegistration technology combined with Global Positioning System (GPS) and Geographic Information System (GIS). Multispectral reflectance data is the set of reflectance values ​​for each band at a specific pixel location in the grid, reflecting the comprehensive spectral reflectance characteristics of the land cover at that location. It can be obtained by extracting the reflectance values ​​for each band after multispectral image radiometric calibration and atmospheric correction. Red reflectance is the surface reflectance value of a pixel location in the red band, reflecting the reflectivity of ground features to red light radiation. Green reflectance is the surface reflectance value of a pixel location in the green band, reflecting the reflectivity of ground features to green light radiation. Red edge reflectance is the surface reflectance value of a pixel location in the red edge band, reflecting the physiological activity and structural characteristics of the vegetation canopy. Near-infrared reflectance is the surface reflectance value of a pixel location in the near-infrared band, reflecting the cellular structure and coverage of the vegetation canopy. The median operator result is a statistical value calculated from the median of the single-band reflectance values ​​of all pixel locations in the same grid, reflecting the central distribution characteristics of reflectance in that band. The absolute deviation is the absolute difference between the band reflectance of a single pixel location and the result of the median operator for that band, reflecting the degree of deviation of the reflectance of a single pixel from the central distribution of the band. The median absolute deviation operator result is a statistical value calculated from the median of the absolute deviation of a certain band across all pixel locations within the same grid, reflecting the discrete distribution characteristics of the reflectance values ​​for that band. Normalized reflectance is the value of the band reflectance after standardization processing using the median operator result and the median absolute deviation operator result, reflecting the relative reflectance characteristics of the pixels after eliminating scale and extreme value interference. Temporal remote sensing reflectance data within the same grid is a numerical set of normalized reflectance values ​​for each band in band order at a specific time phase, reflecting the comprehensive multispectral reflectance characteristics of the forest canopy after standardization at that time phase.

[0051] It should be noted that the preset rule for the pre-defined grid is to spatially divide the monitoring area using equally spaced square grids. The technical specifications are a pixel resolution of 1 to 10 meters, a Gauss-Kruger projection, and the 2000 National Geodetic Coordinate System. For example, when monitoring a mixed uneven-aged forest, the pre-defined grid is set to a 5-meter resolution square grid, using a Gauss-Kruger 3-degree zone projection, with the 2000 National Geodetic Coordinate System as the coordinate reference. All multispectral images from all time phases are registered to this grid reference. Positive constants are preferably set to 10 to the power of -5 to avoid cases where the denominator is zero. The range of all pixel positions in the pre-defined monitoring area for forestry resource dynamic monitoring is defined as all valid pixel positions in the same grid within the pre-defined monitoring area, excluding invalid pixels and pixels without data outside the monitoring area. The time intervals for each phase are selected based on the dynamic monitoring needs of forestry resources: 1 to 3 months for short-term monitoring, 6 to 12 months for medium-term monitoring, and 1 to 3 years for long-term monitoring. The selection criteria are clear skies without clouds and similar solar altitude angles, with the difference in solar altitude angle controlled within 10 degrees to avoid interference from weather and solar altitude angle on reflectance data. The sensor types for multispectral imagery can include satellite sensors, UAV sensors, and airborne sensors. The wavelength range is: red band 620 nm to 680 nm, green band 520 nm to 580 nm, red-edge band 730 nm to 770 nm, and near-infrared band 780 nm to 900 nm.

[0052] It should be noted that this invention addresses the issue of geometric and radiometric scale differences in multispectral images across different time phases. By establishing a unified grid reference for image geometric registration, and after extracting multi-band reflectance data, robust statistical methods based on median and median absolute deviation are used for numerical standardization to eliminate interference from extreme values ​​and scale differences. The data is then assembled in band order to form standardized, grid-based temporal remote sensing reflectance data, providing a unified spectral data foundation for subsequent canopy structure analysis in forestry resource dynamic monitoring. This allows for comparison of multispectral images from different time phases at the same geometric and radiometric scale, reducing the interference of extreme values, improving the statistical robustness of reflectance data, eliminating numerical scale differences between different time phases and bands, preserving the band characteristics of multispectral reflectance, ensuring the consistency and comparability of subsequent canopy structure analysis data, and providing reliable standardized input data for subsequent steps such as canopy occupancy field construction.

[0053] In one embodiment of the present invention, a canopy occupancy field is constructed using temporal remote sensing reflectance data from the same grid, and a continuous canopy occupancy map is generated by binarization, including:

[0054] In each phase With each cell position in the same grid Read the remote sensing reflectance data of the same grid in the same time phase. Obtain the normalized reflectance of red light Green light normalized reflectance Red-edge normalized reflectance Near-infrared normalized reflectance ;

[0055] In each phase With each cell position Calculate the canopy occupancy field :

[0056] ;

[0057] in For the canopy to occupy the field, It is a positive constant;

[0058] Canopy occupying the field As input, determine the binarization threshold. :

[0059] ;

[0060] in For inter-class variance, For candidate thresholds, and These are the proportions of two types of pixels, and These are the mean values ​​of the two types of pixels, This is the population mean;

[0061] In each phase With each cell position Using binarization threshold Canopy Occupied Field Processing to generate continuous canopy occupancy maps :

[0062] ;

[0063] in This is a map showing the continuous canopy occupancy.

[0064] It should be noted that the normalized reflectance of red light, normalized reflectance of green light, normalized reflectance of red edge, and normalized reflectance of near-infrared light reflect the relative reflectance characteristics of the corresponding bands after eliminating scale and extreme value interference. The canopy occupancy field is a continuous value calculated by combining normalized reflectance of specific bands, reflecting the probability that a pixel location is occupied by the upper canopy. The maximization of inter-class variance criterion is a statistical criterion for automatic thresholding, reflecting the judgment logic of achieving pixel classification by maximizing inter-class differences. The binarization threshold is a critical value calculated by maximizing the inter-class variance criterion, reflecting the numerical boundary standard between canopy-occupied and non-occupied pixels. Inter-class variance is the variance statistical value of two classes of pixels under different candidate thresholds, reflecting the degree of numerical difference between the two classes of pixels. Candidate thresholds are alternative critical values ​​used to calculate inter-class variance within the range of canopy occupancy field values. The first value is the identifier value assigned when the canopy occupancy field value is greater than or equal to the binarization threshold, reflecting the canopy-occupied pixel identifier setting, which is 1. The second value is the identifier assigned when the canopy occupancy value is less than the binarization threshold, reflecting the setting of non-occupied pixels in the canopy, which is 0. The continuous canopy occupancy map is a raster image composed of binarized identifier values, reflecting the spatial occupancy distribution characteristics of the upper canopy within the monitoring area.

[0065] It should be noted that the specific steps for calculating the inter-class variance using the maximization of inter-class variance criterion are as follows: First, determine the range of candidate threshold values ​​as the minimum to maximum value of the canopy occupancy field; second, traverse all candidate thresholds with a constant step size of 0.01; third, for each candidate threshold, classify pixels into canopy-occupied and non-occupied classes; fourth, calculate the proportion, mean, and overall mean of the two classes; fifth, calculate the inter-class variance based on the proportion and mean; sixth, select the candidate threshold corresponding to the maximum inter-class variance as the binarization threshold. The selection range of candidate thresholds is the minimum to maximum value of the canopy occupancy field within the monitoring area, and the traversal method is a linear traversal with a constant step size. The step size is initially set to 0.01, which ensures the accuracy of the threshold calculation while controlling the computational load and adapting to the processing efficiency of forestry remote sensing data. The canopy cell structure of healthy vegetation exhibits high reflectivity in the near-infrared band. The red edge band is sensitive to changes in vegetation physiological activity and canopy structure. Canopy-occupied areas show high red edge and near-infrared normalized reflectivity values; summing these values ​​further amplifies the numerical signal of canopy-occupied areas. Red light is easily absorbed by vegetation chlorophyll, resulting in low red reflectivity in canopy-occupied areas, while non-canopy components such as surface soil and dead branches have high red reflectivity. Green light is easily affected by the superposition of canopy shadows and surface background, leading to greater fluctuations in green reflectivity values ​​in non-canopy areas. Combining these two factors effectively characterizes the numerical features of non-canopy areas. Furthermore, continuous canopy occupancy fields only reflect the probability of canopy occupancy and cannot be directly used for subsequent topological analysis. By converting these fields to binary 0s and 1s through comparative assignment, the spatial boundaries between canopy occupancy and non-occupancy can be clearly defined, providing structured foundational data for subsequent extraction of continuous canopy envelopes and pore sets.

[0066] It should be noted that this invention utilizes the normalized reflectance of each band after standardization, combined with the spectral characteristics that the red edge and near-infrared light have strong responses to the canopy, while red and green light are easily affected by non-canopy components. By constructing a canopy occupancy field through specific band combinations, the numerical difference between canopy occupancy and non-occupancy is amplified. Then, the binarization threshold is adaptively determined using the maximization of inter-class variance criterion, converting the continuous canopy occupancy field into a binarized continuous canopy occupancy map. This provides structured spatial distribution data for subsequent extraction of the continuous canopy envelope and pore set. It aligns with the spectral response characteristics of forest canopies, strengthens the numerical representation of canopy occupancy, and the adaptive threshold determination avoids subjective bias from manual pre-setting. The binarized canopy occupancy map clearly defines the spatial boundary of canopy occupancy, adapting to the numerical computation requirements of subsequent topological structure analysis and improving the adaptability of canopy occupancy determination to different monitoring areas.

[0067] In one embodiment of the present invention, a continuous canopy occupancy map is subjected to hole filling to extract the continuous canopy outer envelope, and the area of ​​the outer envelope is calculated. An inner hole set is generated by subtracting the continuous canopy occupancy map from the continuous canopy outer envelope, including:

[0068] In each phase With each cell position in the same grid , with continuous canopy occupancy map Perform hole-filling operation on the input. , thus obtaining a continuous canopy outer envelope The formula is:

[0069] ;

[0070] in For external background collection, For background collection, For the set of cells on the same grid boundary, For the same grid region;

[0071] In each phase With each monitoring unit Construct a continuous canopy outer envelope set within the monitoring unit And calculate the outer envelope area. The formula is:

[0072] ;

[0073] in For monitoring unit In time phase The continuous outer envelope area of ​​the canopy, For set area operator;

[0074] In each phase With each cell position in the same grid , consisting of a continuous canopy outer envelope Subtracting continuous canopy occupancy map Generate an inner hole set The formula is:

[0075] ;

[0076] And construct the internal hole assembly within the monitoring unit. ,in This is a binary representation of the set of internal holes.

[0077] It should be noted that the continuous canopy envelope is the raster data obtained after filling holes in the continuous canopy occupancy map, reflecting the spatial distribution characteristics of the shell formed by the continuous canopy as a whole. The background set is the set of all non-occupied pixels in the continuous canopy occupancy map, reflecting the spatial distribution characteristics of pixels not covered by the canopy within the monitoring area. The set of pixels at the same grid boundary is the set of all pixels located at the geographic boundary within the same grid, reflecting the spatial boundary range characteristics of the same grid. Background connectivity paths are the spatial connections between background pixels, reflecting the spatial connectivity characteristics of non-occupied pixels. The external background set is the set of pixels in the background set that have connectivity paths with the set of pixels at the same grid boundary, reflecting the spatial distribution characteristics of the background outside the canopy within the monitoring area. The hole set is the set of pixels remaining in the background set after removing the external background set, reflecting the spatial distribution characteristics of non-occupied pixels surrounded within the continuous canopy. The monitoring unit is the basic spatial unit for dynamic monitoring of forestry resources, reflecting the zoning spatial characteristics of forestry management. The continuous canopy envelope set within a monitoring unit is the set of pixels within the continuous canopy envelope that fall within the geographic range of the monitoring unit, reflecting the spatial characteristics of the canopy shell within a single monitoring unit. The envelope area is the geographic area corresponding to the continuous canopy envelope set within the monitoring unit, reflecting the spatial coverage scale of the continuous canopy shell within a single monitoring unit. The aperture set is the binary raster set obtained by subtracting the continuous canopy occupancy map from the continuous canopy envelope, reflecting the spatial distribution characteristics of the cavities enclosed within the continuous canopy. The aperture set within a monitoring unit is the set of pixels within the aperture set that fall within the geographic range of the monitoring unit, reflecting the internal cavity spatial characteristics of the continuous canopy within a single monitoring unit.

[0078] It should be noted that the hole filling operation uses the FloodFill algorithm. The execution steps are as follows: First, extract the background set from the continuous canopy occupancy map; second, using the boundary cells of the same grid as seed points, perform flood filling on the background set to mark the outer background set; third, determine the unmarked background cells as the hole set; fourth, assign the cells within the hole set to the canopy occupancy, completing the filling and generating the continuous canopy outer envelope. The background connectivity path adopts the 8-neighbor connectivity determination rule, that is, the cells in the eight directions of a certain background cell (up, down, left, right, left-up, left-down, right-up, right-down) are all considered as connected neighbors; the connectivity analysis uses the seed filling method, using the grid boundary background cells as seeds, traversing all 8-neighbor connected background cells to complete the marking of the outer background set. The numerical determination rule for the binarization of the aperture set is as follows: if the result of subtracting each pixel is 1, it is determined to be an aperture pixel; if the result is 0, it is determined to be a non-aperture pixel. The assignment standard is that aperture pixels are assigned a value of 1, and non-aperture pixels are assigned a value of 0, forming a binarized aperture set raster. In addition, a spatial overlay and clipping method of raster and vector is used to spatially overlay the polygonal vector boundary of the monitoring unit with the binarized aperture set raster, extract all raster pixels within the vector boundary range to form the aperture set within the monitoring unit, and remove all pixels outside the boundary.

[0079] It should be noted that this invention addresses the characteristics of internal holes and external background in continuous canopy occupancy maps. It distinguishes between the two through hole-filling operations, filling only the internal holes, and extracting the continuous canopy outer envelope that characterizes the canopy shell. Then, according to forestry management-related monitoring unit division standards, it trims the outer envelope set within each unit and calculates its area. Finally, through pixel-by-pixel difference calculation between the outer envelope and the original occupancy map, it separates the set of internal holes within the continuous canopy, simultaneously trimming the set of internal holes within each unit. This provides unitized structural data for subsequent internal hole topological feature analysis. This accurately distinguishes the external background of the canopy from internal holes, ensuring that the outer boundary of the canopy outer envelope remains unchanged. It adapts to the operational needs of forestry unit-based monitoring, effectively separating the set of internal holes from the canopy shell. The obtained unitized outer envelope area and internal hole set data provide a precise spatial structural basis for subsequent topological fragmentation analysis, while allowing the structural analysis results to directly match the unit division logic of forestry management.

[0080] In one embodiment of the present invention, extracting the skeleton based on the boundary distance of the inner hole set and determining the intrinsic neck size based on the local hole width of the skeleton includes:

[0081] In each phase With each monitoring unit Within the range, the collection of internal holes within the monitoring unit For each cell position, as input. Calculate the boundary distance of the internal hole assembly The formula is: ;

[0082] in For Euclidean distance operators, For the collection of internal holes The boundary pixel set;

[0083] In each phase With each monitoring unit Within the range, the internal aperture collection within the monitoring unit Extract skeleton from input The formula is:

[0084] ;

[0085] in For the skeleton point set of the internal hole assembly, For skeleton points;

[0086] For each skeleton point Calculate the local hole width at the skeleton point. And calculate the median width of the local holes in the skeleton. The intrinsic neck scale is calculated using the following formula. : ;

[0087] in It is an intrinsic neck scale.

[0088] It should be noted that the boundary distance of the inner aperture set is the spatial distance from each pixel within the inner aperture set to the boundary of the inner aperture set. The skeleton is the set of points formed by the centers of all the largest inscribed circles within the inner aperture set, reflecting the core topological structure and main channel spatial distribution characteristics of the inner aperture set. The set of centers of the largest inscribed circles is the set of centers of the largest inscribed circles at each location within the inner aperture set, reflecting the core constituent point characteristics of the inner aperture set skeleton. A skeleton point is a single center point that constitutes the inner aperture set skeleton. The local aperture width of a skeleton point is twice the boundary distance of the inner aperture set corresponding to the skeleton point, reflecting the local spatial width characteristics of the inner aperture at that skeleton point location. The median of the local aperture widths of the skeleton reflects the central distribution characteristics of the local aperture widths of the inner aperture set skeleton. The intrinsic neck scale is the value obtained by statistically analyzing the median of the local aperture widths of the selected smaller skeleton points, reflecting the typical scale characteristics of the narrow neck positions where the inner aperture set is most prone to topological clipping.

[0089] It should be noted that the boundary of the inner hole set is extracted using the Canny algorithm in edge detection. The criterion for determining a boundary cell is that there is a non-inner hole cell in the neighborhood of the inner hole set; that is, an inner hole cell with a value of 0 in its 8-neighborhood is determined as a boundary cell. The method for solving the maximum inscribed circle is to combine the Euclidean distance transformation result. The radius of the maximum inscribed circle at a certain position of the inner hole is the boundary distance at that position, and the center of the circle is the cell coordinates at that position. The rule for determining the center coordinates is to use the geographic coordinates of the cell's center, based on the geographic coordinate system of the same grid, accurate to the meter. The skeleton extraction uses the median transformation algorithm (MAT) based on distance transformation. The parameter settings require connectivity to be 8-neighborhood and the distance threshold to be set to 1, that is, to remove small skeleton branches with a distance value less than 1 and retain the core skeleton structure. The coordinate extraction method for skeleton points involves extracting the center geographic coordinates of skeleton pixels, outputting two-dimensional coordinate values. The selection criteria for valid skeleton points are to remove isolated single skeleton pixels and retain skeleton points with at least one 8-neighbor skeleton point to avoid interference from invalid skeleton points caused by noise. Invalid data removal rules involve removing values ​​with local aperture widths less than 1. These values ​​correspond to tiny gaps in the inner hole, not actual cavity structures; removing them ensures the statistical results closely match the actual topological features of the inner hole. Furthermore, skeleton smoothing employs a moving average method, using three consecutive skeleton points as a window and calculating the average coordinates within the window as the new coordinates of the intermediate point. Noise reduction uses a length threshold method, removing skeleton branches with a pixel length less than 5 and retaining core skeleton branches with a length greater than or equal to 5, improving the integrity and accuracy of the skeleton.

[0090] It should be noted that this invention addresses the structural scale differences of borehole sets within different monitoring units. Based on the borehole sets within a monitoring unit, it calculates the distance from each pixel to the boundary using Euclidean distance. A median transformation algorithm is employed to extract and preserve the skeleton of the core borehole topology. Twice the distance from the skeleton point boundary is defined as the local borehole width. After two rounds of median statistics and filtering, interference from extremely wide borehole regions is eliminated. The intrinsic neck scale, which characterizes the narrow neck feature of the borehole, is determined, providing an adaptive scale benchmark for subsequent multi-scale erosion analysis of the borehole. This accurately characterizes the spatial location and local width features of borehole pixels. The extracted skeleton completely preserves the core topology of the borehole. Multi-step median processing effectively reduces extreme value interference. The intrinsic neck scale accurately points to the key parts of the borehole topological breakup. Furthermore, the analysis process is adapted to the operational logic of forestry unit-based monitoring, providing an adaptive scale basis that fits the borehole's own structure for subsequent multi-scale erosion analysis targeting topological fragmentation mechanisms. This allows subsequent erosion operations to better align with the actual topological characteristics of the borehole.

[0091] In one embodiment of the present invention, a discrete erosion scale sequence is divided based on the intrinsic neck scale, and a residual erosion set sequence is generated by progressively eroding the pore set. The topological spectrum intensity of the continuous canopy pores is calculated by combining the Eulerian features of the residual pore set sequence, including:

[0092] In each phase With each monitoring unit Within the range, using intrinsic neck scale Calculate the number of discrete erosion scales : ;

[0093] in For discrete erosion scales, The floor operator;

[0094] For each scale label satisfy Calculate the erosion radius :

[0095] ;

[0096] Using radius Circular structure element For internal hole assembly Perform morphological erosion calculations to obtain the set of residual internal pores. : ;

[0097] in For morphological erosion operators; combine all sets of residual internal pores in the layers to generate a sequence of residual internal pore sets. ;

[0098] For each layer of residual internal pores Calculate Euler features :

[0099] ;

[0100] in It is an Euler feature. The number of connected components. This represents the number of internal loops.

[0101] In each phase With each monitoring unit Within the range, utilizing the outer envelope area Calculate the topological spectrum intensity of the continuous canopy pores :

[0102] ;

[0103] in The intensity of the topological spectrum of continuous canopy pores. It is a positive constant.

[0104] It should be noted that the discrete erosion scale sequence is an ordered set of erosion scales divided based on the intrinsic neck scale, reflecting the hierarchical characteristics of the erosion scales adapted to the narrow neck structure of the internal pores. The discrete erosion scale number is a positive integer obtained by rounding up the intrinsic neck scale, reflecting the total number of scale levels for the hierarchical erosion of the internal pore set. The scale label is the sequence number used to identify each level in the discrete erosion scale sequence, reflecting the order of erosion and the position of the scale level. The erosion radius is the radius of action of the morphological erosion operation corresponding to each scale label, reflecting the contraction amplitude of the internal pore set at the corresponding scale. The disk structural element is the circular action template used in the morphological erosion operation, reflecting the spatial action form of the internal pore set contracting uniformly from the boundary to the center. The residual internal pore set is the set of cavities remaining after the internal pore set has undergone the erosion operation at a certain scale label, reflecting the topological structure characteristics of the internal pores at that erosion scale. The Eulerian feature is a value obtained by the difference between the number of connected components and the number of internal loops, reflecting the overall topological structure characteristics of the residual internal pore set at a single scale. The number of connected components is the number of independent connected regions in the residual pore set, reflecting the spatial connectivity distribution characteristics of the residual pore set. The number of internal loops is the number of closed loops in the residual pore set, reflecting the spatial pore nesting characteristics of the residual pore set. The continuous canopy pore topological spectrum intensity is a single-valued index obtained by summing and normalizing the Eulerian characteristics, reflecting the overall topological fragmentation characteristics of the pore set under multi-scale erosion.

[0105] It should be noted that the intrinsic neck scale is the typical narrow neck scale where the inner hole is most prone to topological breakage. The erosion scale sequence, based on this, has its largest scale corresponding to the narrow neck scale. Gradual erosion operates from the lower scales to the narrow neck scales, closely mirroring the actual topological fragmentation process of the inner hole, which begins to break at the narrow neck and then gradually disintegrates. Morphological erosion must be performed step-by-step according to integer scale levels. Continuous intrinsic neck scales cannot directly match this operational logic; the floor operation converts them to positive integers, ensuring the number of erosion scales is integer, allowing the step-by-step erosion operation to be executed in an orderly manner. The spatial action of the disk-shaped structural element is isotropic; during erosion, it contracts simultaneously inward from all directions of the inner hole boundary without directional bias, consistent with the characteristic of the inner hole uniformly fragmenting from the boundary in its natural state. If rectangular structural elements are used, directional contraction will occur, deviating from the actual fragmentation process. Using disk-shaped structural elements to perform morphological erosion achieves uniform contraction of the inner hole assembly from the boundary to the center, restoring the natural process of inner hole topological fragmentation. In addition, each erosion operation causes the pore set to shrink and undergo topological changes to a certain extent. Low-scale erosion only removes tiny edges, while high-scale erosion causes narrow necking and pore fragmentation. The residual pore sets at each scale are combined in sequence, which can fully demonstrate the entire topological change process of the pore from intact to fragmented.

[0106] It should be noted that the scale marker ranges from 1 to the discrete erosion scale number as a continuous positive integer. The numbering rule is to number them sequentially according to the erosion radius, with scale marker 1 corresponding to the minimum erosion radius and scale marker equal to the discrete erosion scale number corresponding to the maximum erosion radius. The generation rule for the disk structure element is to generate an isotropic circular template with the erosion radius as the circular radius. The raster pixel representation is to assign a value of 1 to all pixels within the circular area with the template center as the origin, and assign a value of 0 to pixels outside the area. During pixelation, a rounding rule is used to determine whether a pixel is within the circular area. The morphological erosion operation adopts the basic erosion algorithm of binary morphology. The execution logic is to slide the disk structure element pixel by pixel on the inner hole set raster. Only when the structure element completely falls into the inner hole set pixels is the central pixel retained as 1; otherwise, it is assigned a value of 0. The connectivity determination rule adopts the 8-neighborhood connectivity rule. The method for counting the number of connected components is to use the seed filling method. Starting from any unmarked valid inner hole cell, all its connected neighbor cells are traversed and marked as a connected component. This operation is repeated until all valid inner hole cells are marked. The number of marked cells is the number of connected components. The 8-neighborhood rule is used as the criterion for determining connected neighbors.

[0107] It should be noted that this invention addresses the multi-scale evolution characteristics of internal hole topological fragmentation. Using the intrinsic neck scale as a baseline, a discrete erosion scale number is obtained by rounding up. The discrete erosion scale sequence is then divided, and the erosion radius of each scale is calculated. Disk structural elements are used to perform step-by-step morphological erosion on the internal hole set, generating a sequence of residual internal holes reflecting topological evolution. Eulerian features at each scale are calculated using the number of connected components and internal loops. These features are summed and normalized by combining the discrete erosion scale number and the outer envelope area to obtain a single-valued topological spectrum intensity integrating multi-scale topological features. This provides a quantitative indicator for subsequent inter-layer decoupling analysis. This allows the erosion scale to adapt to the narrow neck structure of the internal hole itself, enabling step-by-step erosion to completely reconstruct the internal hole topological evolution process. Eulerian features effectively quantify the single-scale topological structure. The integration of multi-scale features into a single-valued indicator facilitates subsequent quantitative analysis and inter-unit comparison. Normalization eliminates scale differences between monitoring units, improving the comparability of analysis results from different units. The obtained topological spectrum intensity accurately characterizes the overall topological fragmentation characteristics of the internal hole, providing a reliable quantitative basis for the construction of subsequent inter-layer decoupling topological indices.

[0108] In one embodiment of the present invention, an outer envelope stable modulation term is constructed based on the temporal variation of the outer envelope area, and the interlayer decoupling topological index is obtained by fusing the topological spectrum intensity of the continuous canopy pores with the outer envelope stable modulation term, including:

[0109] In each phase With each monitoring unit Within the range, using the outer envelope area of ​​the current time phase outer envelope area compared to the previous time phase Calculate the relative change in the outer envelope area : ;

[0110] in This represents the relative change in the outer envelope area. For absolute value operators, It is a positive constant;

[0111] In each phase With each monitoring unit Within the range, the relative change in the outer envelope area is used. Constructing an outer envelope-stabilized modulation term :

[0112] ;

[0113] in For outer envelope stable modulation term, It is a natural exponential function;

[0114] In each phase With each monitoring unit Within the range, the topological spectrum intensity of continuous canopy pores With outer envelope stable modulation term Inter-layer decoupling topology index is obtained through fusion calculation. :

[0115] ;

[0116] in This is the interlayer decoupling topology index.

[0117] It should be noted that the relative change in the outer envelope area is the ratio of the absolute difference between the outer envelope area of ​​the current time phase and the previous time phase to the sum of the outer envelope area of ​​the previous time phase and a positive constant, reflecting the degree of relative change in the outer envelope area of ​​the continuous canopy over time. The outer envelope stability modulation term is the result of calculating the negative of the relative change in the outer envelope area using the natural exponential function, reflecting the stability of the continuous canopy outer envelope between time phases. The interlayer decoupling topological index is the result of multiplying the intensity of the continuous canopy inner pore topological spectrum with the outer envelope stability modulation term, reflecting the comprehensive quantitative characteristics of topological fragmentation of the continuous canopy inner pore under the premise of outer envelope stability. The time-correspondence matching rule is that adjacent time phases within the same monitoring unit, season, and imaging conditions are considered the current and previous time phases. The imaging conditions require that the difference in solar altitude angle be within 10 degrees and that all images be clear and cloudless. The time interval is defined according to monitoring needs: 1 to 3 months for short-term dynamic monitoring, 6 to 12 months for medium-term monitoring, and 1 to 3 years for long-term monitoring. Since the outer envelope area base values ​​of different monitoring units are different, the absolute change value cannot objectively reflect the degree of change. The relative change range, with the base value as a reference, can make the degree of area change of units at different scales comparable.

[0118] It should be noted that this invention addresses the problem that the intensity of the topological spectrum within continuous canopy pores cannot accurately reflect the stability of the outer envelope. By calculating the temporal relative variation of the outer envelope area, the influence of unit scale differences on the area change characterization is eliminated. Then, a monotonically decreasing outer envelope stability modulation term is constructed using a natural exponential function to accurately characterize the degree of outer envelope stability. Finally, the modulation term is multiplied by the topological spectrum intensity, internalizing the outer envelope stability condition into the index. This generates an interlayer decoupled topological index that simultaneously integrates the degree of topological fragmentation within the pores and the stability of the outer envelope, thus locking in the target monitoring scenario of outer shell stability and internal fragmentation. This makes the degree of outer envelope area variation comparable across monitoring units of different scales, allows for a detailed characterization of the outer envelope stability, enables the multiplicative fusion to achieve the inherent constraints of the stability condition, and integrates dual core features into the single-value index, making subsequent dynamic change analysis more aligned with the actual monitoring scenario requirements and improving the matching degree between the quantitative index and the target monitoring scenario.

[0119] In one embodiment of the present invention, the temporal dynamic change is calculated based on the interlayer decoupling topology index, and the temporal dynamic change is thresholded to generate a monitoring and judgment result, including:

[0120] In each phase With each monitoring unit Within the range, utilize the inter-layer decoupling topology index of the current time phase. Interlayer decoupling topological index compared to the previous time phase Calculate the dynamic changes in the time phase The formula is: ;

[0121] in For monitoring unit In time phase The temporal dynamic changes;

[0122] For each time phase The set of temporal dynamic changes of all monitoring units Calculate the threshold for the input The formula is: ;

[0123] in For inter-class variance, For candidate thresholds, and The proportions of the two types of monitoring units are respectively. and These represent the average temporal dynamic changes of the two types of monitoring units, respectively. This represents the overall mean of the dynamic changes over time. The total number of monitoring units;

[0124] In each phase With each monitoring unit Within the range, using threshold dynamic changes in phase Perform threshold segmentation and generate monitoring and judgment results. The formula is:

[0125] ;

[0126] in For monitoring and judgment results.

[0127] It should be noted that the temporal dynamic change is the difference between the current temporal interlayer decoupling topological index and the previous temporal interlayer decoupling topological index, reflecting the temporal variation amplitude and trend of the degree of continuous canopy pore topological fragmentation within the monitoring unit. The threshold is a critical value calculated using the maximization of inter-class variance criterion on the set of temporal dynamic changes of all monitoring units, reflecting the numerical boundary standard for distinguishing whether a monitoring unit has undergone target-type structural changes. The ratio of the two types of monitoring units is the ratio of the number of monitoring units with target changes to the number of monitoring units without target changes under a certain candidate threshold, reflecting the quantitative distribution characteristics of the two types of monitoring units. The overall mean of the temporal dynamic change of all monitoring units is the average value of the temporal dynamic change of all monitoring units in the current period, reflecting the average characteristics of the pore topological fragmentation changes of the overall monitoring units in the current period. The monitoring judgment result is the binary identification result obtained after threshold segmentation of the temporal dynamic change, reflecting whether the monitoring unit has undergone target-type structural changes of continuous canopy pore topological fragmentation.

[0128] It should be noted that the interlayer decoupling topology index is a comprehensive quantitative value of the internal porosity fragmentation and the stability of the outer envelope. The difference directly reflects the change in the current period relative to the previous period; a positive increment indicates increased fragmentation, a negative increment indicates decreased fragmentation, and a zero increment indicates no significant change. The degree of internal porosity topological change varies across different monitoring periods, and the data distribution characteristics differ. Determining the threshold based on all data for the current period ensures that the boundary standard aligns with the actual situation and avoids insufficient adaptability of fixed thresholds. The selection range for candidate thresholds is from the minimum to the maximum value of the effective phase dynamic change in the current period; the traversal step size is preferentially set to 0.001 to ensure the accuracy of threshold determination and control the computational load; the screening rule is to remove candidate thresholds with zero inter-class variance during the traversal process, retaining only the thresholds corresponding to the effective variance. The first value is preferentially set to 1 to identify monitoring units that have undergone target-type structural changes; the second value is preferentially set to 0 to identify monitoring units that have not undergone target-type structural changes.

[0129] It should be noted that this invention addresses the core need for dynamic monitoring of forestry resources. It calculates the temporal dynamic change value of the interlayer decoupling topological index to directly characterize the degree of temporal change in porosity topological fragmentation. Based on the dynamic changes of all monitoring units in the current period, a threshold is adaptively determined using the maximization of inter-class variance criterion, avoiding subjective bias from manual pre-setting. Finally, the dynamic changes are binarized through threshold segmentation to generate a monitoring judgment result that clearly defines whether a target change has occurred, thus filtering out monitoring units where porosity topological fragmentation has occurred. This allows the dynamic changes to accurately reflect the temporal trend of porosity fragmentation, the adaptive threshold to fit the actual distribution characteristics of the current data, and the binarized monitoring judgment result to make the selection of changing units clear and explicit. The unit-by-unit execution logic adapts to the operational needs of forestry unit-based monitoring, ensuring the objectivity and repeatability of the judgment result, and providing a clear unit range basis for subsequent location of topological fragmentation.

[0130] In one embodiment of the present invention, based on the monitoring and judgment results, the principal scale position with the largest increase in Eulerian feature in the discrete erosion scale sequence is located. The difference in the residual pore set before and after the principal scale position is compared to determine the main topological fragmentation zone, and the resulting change patch is output via connectivity mapping. This includes:

[0131] In each phase With each monitoring unit Within the range, the monitoring judgment results are met. At that time, within the scale marking range Internal calculation of Euler characteristic amplification and determine the principal scale position. :

[0132] ;

[0133] in Principal scale location For the first Euler features of the residual internal pore set in the layer, The discrete erosion scale number is used; the Euler characteristic calculation formula is:

[0134] ;

[0135] in The number of connected components. This represents the number of internal loops.

[0136] Meeting the monitoring and judgment results At that time, using the principal scale position With the continuous canopy outer envelope set within the monitoring unit Calculate the main generation zone of topological fragmentation :

[0137] ;

[0138] in This is the main zone of topological fragmentation. For set difference operators, For set intersection operators, The set of residual internal holes corresponding to the main scale location. This refers to the collection of continuous canopy outer envelopes within the monitoring unit;

[0139] Meeting the monitoring and judgment results At that time, for the main generation zone of topological fragmentation Perform connected componentization operation and output a set of changed features. :

[0140] ;

[0141] in For a set of changing patches, For the connected domain transformation operator.

[0142] It should be noted that the Euler feature increase is the difference in Euler feature values ​​between the sets of residual pores at adjacent scale markers, reflecting the degree of change in pore topology between adjacent erosion scales. The principal scale position is the scale marker corresponding to the maximum value of the Euler feature increase in the discrete erosion scale sequence, reflecting the core scale level where pore topological fragmentation occurs most significantly. The main pore fragmentation zone is the intersection of the set difference of the residual pore sets before and after the principal scale position and the set of continuous canopy envelopes within the monitoring unit, reflecting the core spatial region where pore topological fragmentation occurs. Variation patches are independent connected components obtained by connecting the main pore fragmentation zone through connected component operations, reflecting independent spatial distribution units of pore topological fragmentation. The set of variation patches is the set of all independent variation patches within the main pore fragmentation zone, reflecting the overall spatial distribution characteristics of pore topological fragmentation within the monitoring unit.

[0143] It should be noted that Eulerian features characterize single-scale topology, and the difference between adjacent scales directly reflects the magnitude of topological changes. A larger difference indicates more severe fracturing and breakage of the pore at that scale. The scale with the maximum increase is the core scale of topological fragmentation. Therefore, the increase in Eulerian features can quantify the degree of change in the pore's topology at adjacent erosion scales, with the maximum increase corresponding to the scale where topological fragmentation is most significant. Topological fragmentation of the pore is a gradual process under progressive erosion, with different scales exhibiting different erosion effects. The principal scale location represents the node with the most severe fragmentation in this process, accurately pinpointing the key scale where topological fragmentation occurs, providing core evidence for locating the spatial occurrence area. The set difference operation may include the background region outside the canopy envelope. This region is unrelated to pore topological fragmentation. By performing an intersection operation with the outer envelope set, only the eroded region within the outer envelope is retained, making the spatial range of the principal occurrence zone more precise. The specific calculation method for the increase in Eulerian feature is to subtract the Eulerian feature of the previous scale from the Eulerian feature of the later scale. No additional determination is needed for the value's sign; a positive value indicates increased topological fragmentation, a negative value indicates decreased fragmentation, and zero indicates no significant change. All values ​​directly participate in determining the principal scale location. When multiple scales have Eulerian feature increases of a maximum value, the scale with the smallest value is selected as the principal scale location, as smaller scales correspond to earlier topological fragmentation and are the initial core scale of the fragmentation process. Furthermore, the sorting rule for the changed patch set is to sort by area from largest to smallest, and for patches with the same area, sort by the x-coordinate of their center geographic coordinates from smallest to largest. The data storage format is SHP vector data, while retaining the patch's unique identifier, area, center coordinates, and other attribute information.

[0144] It should be noted that this invention addresses the problem that unit-level monitoring results cannot reflect the location of spatial changes. For monitoring units determined to have undergone target changes, it locates the most significant principal-scale location of topological fragmentation by calculating the Eulerian feature amplification of adjacent scales. Then, it extracts the erosion reduction region using the set difference of residual pore sets before and after the principal-scale location. Combined with the intersection operation of the canopy outer envelope set, it limits the effective spatial range, obtaining the main generation zone of topological fragmentation. Finally, it divides the main generation zone into independent change patches through connected domain operations, achieving the conversion from unit-level determination to spatialized patches. Focusing on monitoring units that have undergone changes improves the targeting and efficiency of the analysis. The precise location of the principal-scale location points to the core scale of topological fragmentation, the set operation accurately defines the effective spatial area of ​​fragmentation, and the change patches obtained through connected domain operations meet the spatial representation needs of forestry monitoring. This achieves a closed loop from unit-level determination to precise spatial location of topological fragmentation changes, making the monitoring results more aligned with the practical application needs of forestry resource management.

[0145] In one embodiment of the present invention, a method for dynamic monitoring of forestry resources based on remote sensing imagery includes the following steps:

[0146] Step S1: Normalize the multispectral image to generate temporal remote sensing reflectance data for the same grid.

[0147] Step S2: Construct a canopy occupancy field using the same grid temporal remote sensing reflectance data, and binarize it to generate a continuous canopy occupancy map;

[0148] Step S3: Perform hole filling on the continuous canopy occupancy map to extract the continuous canopy outer envelope, calculate the area of ​​the outer envelope, and generate the inner hole set by subtracting the continuous canopy occupancy map from the continuous canopy outer envelope;

[0149] Step S4: Extract the skeleton based on the boundary distance of the inner pore set, and determine the intrinsic neck size based on the local pore width of the skeleton;

[0150] Step S5: Divide the discrete erosion scale sequence according to the intrinsic neck scale, erode the pore set step by step to generate the residual pore set sequence, and calculate the pore topological spectrum intensity of the continuous canopy by combining the Euler features of the residual pore set sequence.

[0151] Step S6: Construct an outer envelope stable modulation term based on the temporal change of the outer envelope area, and fuse the topological spectrum intensity of the continuous canopy pores with the outer envelope stable modulation term to obtain the interlayer decoupling topological index;

[0152] Step S7: Calculate the temporal dynamic change based on the interlayer decoupling topology index, and perform threshold segmentation on the temporal dynamic change to generate monitoring and judgment results;

[0153] Step S8: Based on the monitoring and judgment results, locate the main scale position with the largest increase in Euler feature in the discrete erosion scale sequence, compare the difference in the residual pore set before and after the main scale position to determine the main topological fragmentation zone, and output the change patch by connecting the domains.

[0154] It should be noted that, as Figure 2 As shown, the left side presents multispectral images acquired by a forest farm monitoring unit in June 2023 and June 2024. After reading these images from consecutive time phases, the system automatically generated a canopy occupancy map in the upper right corner to distinguish between the canopy coverage area and the understory gaps. Further combining interlayer decoupling index and topological spectral intensity calculations, the system identified areas where substantial abrupt changes in understory resources occurred. The overlay map of detection results in the lower right corner clearly shows the extracted set of change patches, highlighted in red. The red patches are located within the internal regions of units 4, 9, and 12, with corresponding areas of 65 square meters, 92 square meters, and 78 square meters, respectively. These patches were all converted into vector data with geometric boundaries and geographic coordinates. The entire process was executed automatically without manual intervention, effectively demonstrating the automated processing capabilities and spatial change detection effect of this invention in multi-temporal dynamic monitoring.

[0155] It should be noted that, as Figure 3 As shown in the diagram, a comparative experiment demonstrates the performance difference between the interlayer decoupling index proposed in this invention and the traditional normalized vegetation index (NZVI) difference method in detecting understory mutations. The upper part constructs three typical forest change scenarios: an area with cleared understory shrubs, an area with slight canopy growth, and an area with understory weed growth. The traditional method's detection map in the lower left shows that it is overly sensitive to canopy changes, easily misinterpreting natural canopy growth as mutations, and failing to effectively separate the actual replacement information of understory vegetation. The interlayer decoupling index map in the lower right shows that this invention effectively suppresses interference from canopy state fluctuations by separating the topological changes of the canopy and understory spaces. In the area with cleared understory shrubs, the response value of the interlayer decoupling index is significantly higher. In the area with slight canopy growth, the traditional index shows a high false alarm rate, while the interlayer decoupling index has an extremely low response. In the area with understory weed growth, the new index also maintains excellent detection sensitivity. This comparative result confirms that this invention has significant effectiveness and reliability in overcoming canopy shading and interlayer signal confusion.

[0156] It should be noted that the interval and threshold sizes are set for ease of comparison. The size of the threshold depends on the amount of sample data and the base number set by those skilled in the art for each set of sample data, as long as it does not affect the proportional relationship between the parameter and the quantized value. Furthermore, the above formulas are all dimensionless calculations, and the formulas are derived from software simulations using a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0157] The embodiments of this example have been described above. However, this example is not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms based on the guidance of this example, and all of them are within the protection scope of this example.

Claims

1. A forestry resource dynamic monitoring system based on remote sensing imagery, characterized in that, include: The first module normalizes multispectral images to generate temporal remote sensing reflectance data in the same grid. The second module uses temporal remote sensing reflectance data from the same grid to construct a canopy occupancy field and binarizes it to generate a continuous canopy occupancy map. The third module performs hole filling on the continuous canopy occupancy map to extract the continuous canopy outer envelope, calculates the area of ​​the outer envelope, and generates an inner hole set by subtracting the continuous canopy occupancy map from the continuous canopy outer envelope. The fourth module extracts the skeleton based on the boundary distance of the inner hole set, and determines the intrinsic neck size based on the local hole width of the skeleton; The fifth module divides the discrete erosion scale sequence according to the intrinsic neck scale, generates a residual erosion set sequence by eroding the pore set step by step, and calculates the pore topological spectrum intensity of the continuous canopy by combining the Eulerian features of the residual pore set sequence. The sixth module constructs an outer envelope stable modulation term based on the temporal variation of the outer envelope area, and fuses the topological spectrum intensity of the continuous canopy pores with the outer envelope stable modulation term to obtain the interlayer decoupling topological index. The seventh module calculates the temporal dynamic change based on the interlayer decoupling topology index, and performs threshold segmentation on the temporal dynamic change to generate monitoring and judgment results; The eighth module locates the principal scale position with the largest increase in Eulerian features in the discrete erosion scale sequence based on the monitoring and judgment results. It then compares the differences in the residual pore set before and after the principal scale position to determine the main topological fragmentation zone and outputs the change patch by connecting the domains.

2. The forestry resource dynamic monitoring system based on remote sensing imagery according to claim 1, characterized in that, At each time phase, multispectral reflectance data is generated at any pixel location in the predetermined grid of the multispectral image. The multispectral reflectance data includes red light reflectance, green light reflectance, red edge reflectance, and near-infrared reflectance. For each time phase and each band, the median of the multispectral reflectance data is taken at all pixel locations in the same grid to obtain the median operator result. For each time phase and each band, an absolute deviation is constructed, and the median of the absolute deviation is taken at all pixel positions in the same grid to obtain the median absolute deviation operator result. For each time phase, each band, and each pixel location, the results of the median operator and the median absolute deviation operator are used to perform numerical standardization to calculate the normalized reflectance. For each time phase and each pixel location, the normalized reflectance is assembled into remote sensing reflectance data of the same time phase in band order.

3. The forestry resource dynamic monitoring system based on remote sensing imagery according to claim 1, characterized in that, At each time phase and at each pixel location within the same grid, remote sensing reflectance data for the same time phase of the same grid are read to obtain red light normalized reflectance, green light normalized reflectance, red edge normalized reflectance, and near-infrared normalized reflectance. At each time phase and at each pixel location, the canopy occupancy field is obtained by calculating the ratio of the numerator to the denominator by using the sum of the red edge normalized reflectance and the near-infrared normalized reflectance as the numerator and the red light normalized reflectance and the green light normalized reflectance plus a positive constant as the denominator. Using the canopy occupancy field as input, the binarization threshold is determined by maximizing the inter-class variance criterion. The inter-class variance is calculated using the proportion of the two classes of pixels, the mean of the two classes of pixels, and the overall mean under different candidate thresholds. At each time phase and at each pixel location, the canopy occupancy field is compared with a binarization threshold. When the canopy occupancy field value is greater than or equal to the binarization threshold, a first value is assigned, and when the canopy occupancy field value is less than the binarization threshold, a second value is assigned, thus generating a continuous canopy occupancy map.

4. The forestry resource dynamic monitoring system based on remote sensing imagery according to claim 1, characterized in that, At each time phase and at each pixel position in the same grid, a hole-filling operation is performed with the continuous canopy occupancy map as input to obtain the continuous canopy outer envelope. A background set consisting of all non-occupied pixels is defined, the set of boundary pixels in the same grid is determined, and the external background set connected to the set of boundary pixels in the same grid is extracted according to the background connectivity path. The background pixels in the background set excluding the external background set are determined as the hole set. The continuous canopy outer envelope is generated by filling the hole set. Within each time phase and each monitoring unit, a continuous canopy outer envelope set is constructed within the monitoring unit, consisting of pixels belonging to the continuous canopy outer envelope. The area contained within the continuous canopy outer envelope set within the monitoring unit is calculated to obtain the outer envelope area. At each time phase and at each pixel location in the same grid, the continuous canopy occupancy map is subtracted from the continuous canopy outer envelope to obtain a binary set of internal holes, and the set of internal holes within the monitoring unit is extracted according to the range of the monitoring unit.

5. A forestry resource dynamic monitoring system based on remote sensing imagery according to claim 1, characterized in that, Within each time phase and each monitoring unit, taking the set of inner holes within the monitoring unit as input, the Euclidean distance from each pixel position to the boundary of the set of inner holes is calculated, and the boundary distance of the set of inner holes is obtained. Within each time phase and each monitoring unit, an extraction operation is performed with the set of inner holes within the monitoring unit as input. The set of centers of all the largest inscribed circles within the set of inner holes is determined as the skeleton, and the skeleton points contained in the skeleton are obtained. For each skeleton point, the local aperture width of the skeleton point is defined by twice the distance of the inner aperture set boundary. The local aperture widths of all skeleton points are counted and the median is calculated to obtain the median of the local aperture widths of the skeleton. The local aperture widths of skeleton points with values ​​less than or equal to the median of the local aperture widths of the skeleton are further selected to form a set, and the median is taken to obtain the intrinsic neck scale.

6. A forestry resource dynamic monitoring system based on remote sensing imagery according to claim 1, characterized in that, Within each time phase and each monitoring unit, the intrinsic neck scale is used to perform an up-rounding operation to obtain the discrete erosion scale number. The erosion radius is calculated for each scale marker, and morphological erosion operations are performed on the set of internal pores using disk structural elements with erosion radius to obtain the set of residual internal pores in each layer. The set of residual internal pores in all layers is combined to generate a sequence of residual internal pore sets. For each layer of residual internal holes, the Euler feature is calculated by subtracting the number of internal loops from the number of connected components; The Euler features of the residual pore sets in each layer are summed, and the summation result is divided by the product of the discrete erosion scale number, the outer envelope area, and a positive constant to calculate the topological spectrum intensity of the continuous canopy pores.

7. A forestry resource dynamic monitoring system based on remote sensing imagery according to claim 1, characterized in that, Within each time phase and each monitoring unit, the absolute value of the numerical difference between the outer envelope area of ​​the current time phase and the outer envelope area of ​​the previous time phase is calculated, and the absolute value of the numerical difference is divided by the sum of the outer envelope area of ​​the previous time phase and the positive constant to obtain the relative change range of the outer envelope area. Within each time phase and each monitoring unit, the outer envelope stable modulation term is obtained by performing exponential calculation using the negative of the relative change amplitude of the outer envelope area using the natural exponential function. Within each time phase and each monitoring unit, the topological spectrum intensity of the continuous canopy pores is multiplied with the stable modulation term of the outer envelope to obtain the interlayer decoupling topological index.

8. A forestry resource dynamic monitoring system based on remote sensing imagery according to claim 1, characterized in that, Within each time phase and each monitoring unit, the inter-layer decoupling topology index of the current time phase is subtracted from the inter-layer decoupling topology index of the previous time phase to calculate the dynamic change of the time phase. In each time phase, the set of time-phase dynamic changes of all monitoring units is used as input, and the threshold is determined by maximizing the inter-class variance criterion. The inter-class variance is determined by the ratio of the two types of monitoring units under different candidate thresholds, the mean of the two types of monitoring units, and the overall mean of the time-phase dynamic changes of all monitoring units. Within each time phase and each monitoring unit, a threshold is used to perform a segmentation operation on the dynamic change of the time phase. When the dynamic change of the time phase is greater than or equal to the threshold, the judgment result is determined as the first value, and when the dynamic change of the time phase is less than the threshold, the judgment result is determined as the second value, thus obtaining the monitoring judgment result.

9. A forestry resource dynamic monitoring system based on remote sensing imagery according to claim 1, characterized in that, Within each time phase and within each monitoring unit, under the condition that the monitoring judgment result meets the first value, the Euler feature amplification between the sets of residual inner holes in two adjacent layers is calculated within the scale mark range. The Euler feature amplification is calculated using the numerical difference between adjacent scale marks of the Euler feature determined by the number of connected components and the number of internal loops. The main scale position is located by finding the scale mark corresponding to the maximum value of the Euler feature amplification. Under the condition that the monitoring and judgment results meet the first value, the set of residual internal pores of two adjacent layers is located using the principal scale position. The set difference between the set of residual internal pores of the previous layer corresponding to the principal scale position and the set of residual internal pores corresponding to the principal scale position is calculated. The set difference is then intersected with the set of continuous canopy outer envelope within the monitoring unit to determine the main topological fragmentation zone. Under the condition that the monitoring and judgment results meet the first value, the connectivity operation is performed on the main generation zone of topological fragmentation, dividing the binary set into a set of disconnected connected components, and outputting a set of changed patches.

10. A method for dynamic monitoring of forestry resources based on remote sensing imagery, characterized in that, Implementing a forestry resource dynamic monitoring system based on remote sensing imagery as described in any one of claims 1 to 9 includes the following steps: Step S1: Normalize the multispectral image to generate temporal remote sensing reflectance data for the same grid. Step S2: Construct a canopy occupancy field using the same grid temporal remote sensing reflectance data, and binarize it to generate a continuous canopy occupancy map; Step S3: Perform hole filling on the continuous canopy occupancy map to extract the continuous canopy outer envelope, calculate the area of ​​the outer envelope, and generate the inner hole set by subtracting the continuous canopy occupancy map from the continuous canopy outer envelope; Step S4: Extract the skeleton based on the boundary distance of the inner pore set, and determine the intrinsic neck size based on the local pore width of the skeleton; Step S5: Divide the discrete erosion scale sequence according to the intrinsic neck scale, erode the pore set step by step to generate the residual pore set sequence, and calculate the pore topological spectrum intensity of the continuous canopy by combining the Euler features of the residual pore set sequence. Step S6: Construct an outer envelope stable modulation term based on the temporal change of the outer envelope area, and fuse the topological spectrum intensity of the continuous canopy pores with the outer envelope stable modulation term to obtain the interlayer decoupling topological index; Step S7: Calculate the temporal dynamic change based on the interlayer decoupling topology index, and perform threshold segmentation on the temporal dynamic change to generate monitoring and judgment results; Step S8: Based on the monitoring and judgment results, locate the main scale position with the largest increase in Euler feature in the discrete erosion scale sequence, compare the difference in the residual pore set before and after the main scale position to determine the main topological fragmentation zone, and output the change patch by connecting the domains.

Images