Method and system for reconstructing disturbance history of natural forest based on satellite image analysis

Abstract

The application provides a natural forest disturbance history reconstruction method and system based on satellite image analysis, and relates to the technical field of satellite remote sensing forestry. The method comprises the following steps: according to multi-temporal optical satellite images of a natural broad-leaved mixed forest area, the gap boundary of the forest canopy is extracted to obtain a gap profile; according to the gap profile, the extreme points of the boundary curvature are calculated to construct an initial polygon, the center line of the inscribed circle is extracted to generate a skeleton line, the branch length and the included angle of the skeleton line are calculated to obtain a main skeleton line; according to the main skeleton line, V-shaped geometric edges with the tip pointing to the upper side of the slope are identified to obtain a disturbance candidate area; and according to the disturbance candidate area, the pixel-by-pixel drop amplitude of the red-green band reflection difference along the normal direction of the V-shaped edge is calculated to obtain the mixed features of the slope surface shadow and bare soil. The application improves the spatial coverage, recognition accuracy and time continuity.

Description

Technical Field

[0001] This invention relates to the field of satellite remote sensing forestry technology, and in particular to a method and system for reconstructing the history of disturbances in natural forests based on satellite image analysis. Background Technology

[0002] In natural broadleaf mixed forest areas of southern China, landslides and windbreaks induced by typhoons and heavy rains are among the main types of natural disturbances. After disturbances occur, the canopy is damaged, forming strip-shaped or patchy voids. Fallen trees, debris, and exposed soil are interspersed, creating unique micro-topographic features. The spatiotemporal distribution information of historical disturbance events is of great value for forest succession research, carbon cycle assessment, and ecological restoration planning. When investigating the disturbance history of typical natural broadleaf mixed forest areas in southern China, traditional field sampling methods are sometimes used: researchers traverse fallen tree accumulation areas and landslide gullies along pre-set transects, use GPS to record disturbance boundaries, and infer the year of disturbance by cross-dating tree-ring samples. Such surveys typically cover large areas, are time-consuming, and require significant manpower.

[0003] Traditional field survey methods are limited by terrain accessibility, making it difficult to fully cover steep gullies and dense forest areas with fallen logs, thus compromising the integrity of spatial coverage. Tree-ring dating methods rely heavily on living tree samples; when tracing earlier disturbance events, the continuity of temporal reconstruction may be affected by the limited number of surviving samples. Furthermore, survey results can vary to some extent due to differences in the experience of the interpreters. Therefore, methods for automatically extracting canopy gaps, identifying micro-topographic features caused by disturbances such as landslides, and reconstructing disturbance history based on optical satellite imagery still have room for further research and improvement. Summary of the Invention

[0004] This invention provides a method and system for reconstructing the history of disturbance in natural forests based on satellite image analysis. This invention improves spatial coverage, identification accuracy, and temporal continuity.

[0005] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows:

[0006] Firstly, a method for reconstructing the history of disturbance in natural forests based on satellite image analysis, the method comprising:

[0007] Based on multi-temporal optical satellite images of natural broad-leaved mixed forest areas, the canopy gap boundaries are extracted to obtain the gap contours; based on the gap contours, the extreme points of boundary curvature are calculated to construct initial polygons, the lines connecting the centers of the inscribed circles are extracted to generate skeleton lines, and the branch lengths and included angles of the skeleton lines are calculated to obtain the main skeleton lines;

[0008] Based on the main skeleton line, identify the V-shaped geometric edge with its tip pointing upwards on the slope to obtain the interference candidate area; based on the interference candidate area, calculate the pixel-by-pixel drop amplitude of the red and green band reflection difference along the normal direction of the V-shaped edge to obtain the mixed features of slope shadow and bare soil.

[0009] Based on the mixed characteristics of slope shadow and bare soil, the striped texture with alternating high and low short-wave infrared reflectance is extracted to obtain the micro-topographic features of fallen trees and fallen debris; based on the main framework line and the micro-topographic features of fallen trees and fallen debris, the morphological boundary of the ditch is determined.

[0010] Based on the morphological boundaries of the gullies, the candidate areas of disturbance are divided into subsets of disturbance pixels to obtain the disturbance type classification results. Based on the disturbance type classification results, the periodic changes in shortwave infrared reflectance are accumulated in the subsets of disturbance pixels to obtain the cumulative amount of canopy disturbance recovery.

[0011] Based on the cumulative amount of canopy disturbance recovery, the critical time phase for vegetation spectral recovery and stability is determined, and the disturbance occurrence time window is obtained. Based on the disturbance occurrence time window, the disturbance type classification results are mapped to generate a natural forest disturbance history reconstruction sequence, and the reconstruction results are obtained.

[0012] The second aspect is a natural forest disturbance history reconstruction system based on satellite image analysis, including:

[0013] The extraction module is used to extract the canopy gap boundary from multi-temporal optical satellite images of natural broad-leaved mixed forest areas to obtain the gap outline; based on the gap outline, the boundary curvature extreme points are calculated to construct an initial polygon, the inscribed circle center line is extracted to generate the skeleton line, and the branch length and included angle of the skeleton line are calculated to obtain the main skeleton line;

[0014] The identification module is used to identify the V-shaped geometric edge with its tip pointing upwards on the slope based on the main skeleton line, and obtain the interference candidate area; based on the interference candidate area, the pixel-by-pixel drop amplitude of the red and green band reflection difference along the normal direction of the V-shaped edge is calculated to obtain the mixed features of slope shadow and bare soil.

[0015] The module is used to extract the striped texture with alternating high and low shortwave infrared reflectance based on the mixed characteristics of slope shadow and bare soil, and obtain the micro-topographic features of fallen trees and fallen debris; the morphological boundary of the ditch is determined based on the main framework line and the micro-topographic features of fallen trees and fallen debris.

[0016] The segmentation module is used to divide the candidate area of ​​disturbance into a subset of disturbance pixels according to the morphological boundary of the ditch, and obtain the disturbance type classification result; based on the disturbance type classification result, the periodic change of shortwave infrared reflectance is accumulated in the subset of disturbance pixels to obtain the cumulative amount of canopy disturbance recovery.

[0017] The mapping module is used to determine the critical time phase for vegetation spectral recovery stability based on the cumulative amount of canopy disturbance recovery, and obtain the disturbance occurrence time window; based on the disturbance occurrence time window, the disturbance type classification results are mapped to generate a natural forest disturbance history reconstruction sequence, and the reconstruction results are obtained.

[0018] The above-described solution of the present invention has at least the following beneficial effects:

[0019] By automatically extracting canopy gaps and constructing main framework lines from multi-temporal optical satellite imagery, large-scale interference area identification is achieved, overcoming the limitations of traditional surveys due to terrain accessibility and improving spatial coverage integrity and survey efficiency. By calculating the drop amplitude of red and green band reflectance along the normal direction of the V-shaped edge, the mixed features of slope shadow and bare soil are extracted. Combined with the striped texture of high and low shortwave infrared reflectance, the alternating micro-topography of fallen trees and fallen debris is identified, thus accurately depicting the landslide and wind-blown interference characteristics induced by typhoons and rainstorms, improving the objectivity and accuracy of identification.

[0020] By accumulating the changes in shortwave infrared reflectance over time, the cumulative amount of canopy disturbance recovery is obtained, and the critical time phase for vegetation spectral recovery and stability is determined accordingly. This allows for the inversion of the disturbance occurrence time, effectively overcoming the dependence of tree-ring dating on surviving samples and ensuring the continuity of time reconstruction. By mapping the disturbance occurrence time window with the disturbance type classification results, a natural forest disturbance history reconstruction sequence is automatically generated, reducing the subjective differences in human interpretation. Attached Figure Description

[0021] Figure 1 This is a flowchart illustrating the method for reconstructing the history of natural forest disturbance based on satellite image analysis, provided by an embodiment of the present invention.

[0022] Figure 2 This is a schematic diagram of a natural forest disturbance history reconstruction system based on satellite image analysis provided by an embodiment of the present invention. Detailed Implementation

[0023] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0024] The application scenario focuses on a natural broad-leaved mixed forest area in southern China. Located in a subtropical monsoon climate zone, this area is frequently hit by typhoons during the summer and autumn seasons. The strong winds and torrential rains brought by typhoons trigger landslides and windbreaks within the forest. After the canopy is damaged, strip-shaped or patchy voids are formed, with fallen trees and debris accumulating inside these voids, exposing bare soil and creating a unique micro-topographic combination. Multispectral satellite imagery with a spatial resolution of 0.5 meters is used as the data source, encompassing five spectral bands: red, green, blue, near-infrared, and shortwave infrared. This is supplemented by a digital elevation model (DEM) of the same region, with a spatial resolution consistent with the satellite imagery, providing information on topographic slope and aspect. The imaging time span of the satellite imagery covers the complete time series of the past 10 years within the study area, with the imaging time interval between images determined by the satellite revisit cycle.

[0025] like Figure 1 As shown, embodiments of the present invention propose a method for reconstructing the history of disturbance in natural forests based on satellite image analysis. The method includes the following steps:

[0026] Step 1: Based on multi-temporal optical satellite images of natural broad-leaved mixed forest areas, extract the canopy gap boundary to obtain the gap outline; based on the gap outline, calculate the boundary curvature extreme points to construct an initial polygon, extract the inscribed circle center line to generate the skeleton line, calculate the branch length and included angle of the skeleton line to obtain the main skeleton line;

[0027] Step 2: Based on the main skeleton line, identify the V-shaped geometric edge with the tip pointing upwards on the slope to obtain the interference candidate area; based on the interference candidate area, calculate the pixel-by-pixel drop amplitude of the red and green band reflection difference along the normal direction of the V-shaped edge to obtain the mixed features of slope shadow and bare soil.

[0028] Step 3: Based on the mixed characteristics of slope shadow and bare soil, extract the striped texture with alternating high and low shortwave infrared reflectance to obtain the micro-topographic features of fallen trees and fallen debris; determine the morphological boundary of the ditch based on the main framework line and the micro-topographic features of fallen trees and fallen debris.

[0029] Step 4: Based on the morphological boundary of the trench, divide the candidate area of ​​disturbance into a subset of disturbance pixels to obtain the disturbance type classification result; based on the disturbance type classification result, accumulate the periodic change of shortwave infrared reflectance within the subset of disturbance pixels to obtain the cumulative amount of canopy disturbance recovery.

[0030] Step 5: Based on the cumulative amount of canopy disturbance recovery, determine the critical time phase for vegetation spectral recovery and stability, and obtain the disturbance occurrence time window; based on the disturbance occurrence time window, map the disturbance type classification results, generate the natural forest disturbance history reconstruction sequence, and obtain the reconstruction results.

[0031] Based on multi-temporal high-resolution optical satellite imagery of natural broadleaf mixed forest areas, canopy gap boundaries were extracted to obtain gap contours. Specifically, this involved selecting images from the multi-temporal satellite image set that met the following two preset conditions: the first condition was that the image coverage of the study area reached 100%; the second condition was that the cloud cover ratio of the image within the study area was lower than a preset cloud cover threshold. The preset cloud cover threshold was set to 10%, based on the fact that when the cloud cover ratio of a single image exceeded 10%, the undetected area of ​​canopy gaps caused by cloud cover exceeded an acceptable range. The correlation between the missed detection area and the cloud cover ratio was determined by statistical results of the canopy gap extraction integrity rate of satellite images of the forest area under different cloud cover conditions over the past 5 years. 10% was considered an acceptable inflection point on the statistical curve. For each satellite image obtained after screening, radiometric calibration and atmospheric correction operations were performed sequentially to convert the image's digital count values ​​into surface reflectance. The specific process of radiometric calibration was as follows: radiometric calibration parameters were read from the satellite image metadata. These parameters included a gain coefficient and an offset coefficient. The physical meaning of the gain coefficient was the increment of spectral radiance at the sensor's entrance pupil corresponding to each unit increase in the digital count value, with the unit being the digital count value. The physical meaning of the offset coefficient was the equivalent spectral radiance output value generated by the sensor's dark current under no-light conditions when the digital count value was zero, with the same unit as spectral radiance. For each pixel in the image, the value corresponding to the offset coefficient is subtracted from the digital count value of that pixel, and the difference is multiplied by the gain coefficient. The resulting product is the spectral radiance value of that pixel at the entrance pupil of the sensor. The radiometric calibration operation is performed individually for each pixel and each band, with each band using its own gain coefficient and offset coefficient.

[0032] The specific process of atmospheric correction is as follows: Atmospheric correction is performed using the FLAASH correction method. The FLAASH correction method is based on the MODTRAN radiative transfer model and calculates surface reflectance pixel-by-pixel by solving the atmospheric radiative transfer equation. The input parameters for the FLAASH correction method are configured as follows: the latitude and longitude coordinates of the image center are used to determine the solar zenith angle and azimuth angle; the imaging date and time are used to calculate the solar position parameters at the imaging time; the atmospheric model type is set to tropical type based on the imaging season and the latitude of the image center, suitable for low-latitude regions between 30 degrees north and south latitude; the aerosol model type is set to rural type based on the surface cover characteristics of the study area, suitable for non-urban areas with predominantly vegetation cover, and the default visible light band value for aerosol optical thickness is approximately between 0.1 and 0.3; the initial visibility is set based on meteorological records at the imaging time, with a default value of 40 kilometers used when no meteorological records are available. These input parameters are then input into the FLAASH correction processing module, which runs two stages of forward simulation and inverse inversion of the MODTRAN radiative transfer model pixel-by-pixel. In the forward simulation phase, the atmospheric transmittance, upward atmospheric radiation, and downward atmospheric radiation components for each band are calculated based on the input atmospheric parameters. In the inverse retrieval phase, the upward atmospheric radiation component is subtracted from the spectral radiance value at the sensor's entrance pupil. The difference is then divided by the product of the atmospheric transmittance and the downward atmospheric radiation component after surface reflection. The quotient is the surface reflectance value of that pixel in that band. Atmospheric correction is performed pixel-by-pixel and band-by-band. After atmospheric correction is completed for all bands and all pixels, the surface reflectance image for each band is obtained.

[0033] After radiometric calibration and atmospheric correction, the image was differentiated into canopy and gap areas. Within the canopy area, the spongy mesophyll tissue cell walls of healthy broad-leaved trees undergo multiple scattering of near-infrared radiation at the air gap interface, resulting in high near-infrared reflectivity. In the gap area, the ground surface is mainly covered by bare soil, fallen tree trunks, and litter, leading to lower near-infrared reflectivity compared to healthy canopies. Based on these spectral differences, the Otsu method was used to automatically determine the segmentation threshold for the near-infrared surface reflectivity image. The specific processing steps of the Otsu method for determining the segmentation threshold are as follows: The surface reflectivity values ​​of all pixels in the near-infrared image are used as the processing object. All possible candidate segmentation threshold values ​​are iterated. The search range is from the minimum to the maximum value of the surface reflectance of image pixels, and the search step size is the smallest unit of precision of the reflectance value. For each segmentation threshold candidate value... All pixels are divided into two categories: those with near-infrared surface reflectance greater than 100%. The pixels are classified as high reflectivity. Near-infrared surface reflectance is less than or equal to The pixels are classified as low reflectivity. Statistical analysis Class and Number of pixels for each class , and the average surface reflectance of pixels , Calculate the segmentation threshold Corresponding inter-category differences The calculation formula is:

[0034] ;

[0035] In the formula: express The corresponding inter-category difference is dimensionless; the larger the value, the greater the degree of separation between the two categories in the surface reflectance space. express Percentage of class pixels , dimensionless, with values ​​ranging from 0 to 1; express Percentage of class pixels Dimensionless, with values ​​ranging from 0 to 1. ; express The average near-infrared surface reflectance of all pixels in the class, dimensionless; express The average near-infrared surface reflectance of all pixels in the class, dimensionless, is selected after traversal. Reaching the maximum value As the final segmentation threshold ,in accordance with Binarize the image: Near-infrared surface reflectance is greater than... The pixel is labeled as the canopy category (value is 1), less than or equal to The cells are labeled as gap categories (with a value of 0).

[0036] The eight-neighbor edge tracking algorithm is used to extract boundaries for gap categories in the binarized image. The specific execution process of the eight-neighbor edge tracking algorithm is as follows: Scan the binarized image and locate the first pixel with a value of 0 and at least one value of 1 within its eight neighbors. This pixel is used as the starting pixel for edge tracking. Starting from the starting pixel, search for the next pixel with a value of 0 and at least one value of 1 within its eight neighbors according to a preset clockwise neighborhood search order (using the right direction as the 0° reference, searching clockwise in the order of right, bottom right, bottom left, bottom left, top left, top right, and top right, for a total of 8 neighborhood positions). After finding a neighboring pixel that meets the conditions, record it as the boundary pixel and use it as the new current pixel. Rotate 45° clockwise in the opposite direction from the previous boundary pixel to the current pixel as the new starting direction for the search, and continue the search. Repeat the above process until the current pixel returns to the starting pixel, completing the tracking of a closed boundary line. Perform the above eight-neighborhood edge tracking operation on each gap category connected region in the binarized image to obtain a set of closed gap contour lines. Each gap contour line is represented by an ordered sequence of pixel coordinates, where each pixel records its row and column coordinates in the image. This set of closed gap contour lines is the gap contour.

[0037] Based on the gap contour, the tangent direction deflection angle of adjacent pixels is calculated along the boundary direction to obtain a set of boundary curvature extrema points. Based on the set of boundary curvature extrema points, adjacent extrema points are connected to construct an initial polygon, resulting in the inscribed polygon of the gap. Based on the inscribed polygon of the gap, the lines connecting the centers of the inscribed circles of each side are extracted to generate an initial skeleton line network. The length of each branch line segment and the intersection angle are calculated to remove short branch line segments with a length less than a preset length threshold and obtuse branch line segments with an angle greater than a preset angle threshold, resulting in the main skeleton line. The specific implementation process is as follows:

[0038] To calculate the extreme points of boundary curvature, for each ordered pixel coordinate sequence of the gap contour line, sampling points are uniformly extracted along the boundary direction at a preset sampling interval of 2 pixels. The preset sampling interval of 2 pixels is determined based on the fact that in an image with a spatial resolution of 0.5 meters, 2 pixels correspond to an actual ground distance of 1 meter. This interval ensures the spatial resolution required for curvature calculation while reducing the interference of random position noise of adjacent pixels on the estimation of tangent direction. For each sampling point in the sequence, 5 consecutive sampling points with a preset window radius are taken in both the forward and backward directions. A straight line is fitted using the least squares method. The angle between the fitted line and the due north direction is taken as the forward tangent direction and the backward tangent direction, respectively. The preset window radius of 5 sampling points is determined based on the fact that at a sampling interval of 2 pixels, 5 sampling points cover an actual ground distance of approximately 5 meters, which can smooth out local small-scale fluctuations in the canopy gap boundary (such as the serrated undulations of a single tree canopy edge, which are usually less than 5 meters) while preserving large-scale directional change characteristics. The absolute value of the difference between the forward and backward tangent direction angles is calculated to obtain the tangent direction deflection angle at the sampling point. A preset curvature threshold of 45° is used, based on the fact that in canopy gaps formed by landslides, the boundary tangent direction deflection angles at the V-shaped tip and gully sidewall transition points range from 45° to 135°, with 45° serving as the boundary value distinguishing between gradual directional changes and abrupt directional changes. When the tangent direction deflection angle is greater than 45°, the sampling point is marked as a boundary curvature extremum. All sampling points are traversed, and arranged according to their original order on the contour line to obtain the set of boundary curvature extremum points.

[0039] Construct the inscribed polygon within the gap. The extreme points in the boundary curvature extreme point set are arranged in chronological order to form an ordered sequence. Connect the i-th and (i+1)-th extreme points (i from 1 to N1) sequentially, completing the connection of N1 edge segments with straight line segments; finally, connect the N-th and first extreme points to form a closed loop, obtaining the inscribed polygon within the gap, generating the initial skeleton line network. For each edge of the inscribed polygon within the gap, draw a perpendicular line through the midpoint of that edge; the line segment of this perpendicular line inside the polygon is the candidate trajectory for the inscribed circle center. Select candidate center points along the candidate trajectory with a preset step size of 0.5 pixels. The preset step size of 0.5 pixels is based on the fact that at a spatial resolution of 0.5 meters, it corresponds to an actual ground distance of 0.25 meters, which is sufficient to ensure the accuracy of the inscribed circle center positioning. For each candidate center point, calculate the shortest perpendicular distance from it to all edges of the inscribed polygon within the gap as the upper limit of the inscribed circle radius; take the candidate center point corresponding to the maximum radius upper limit as the inscribed circle center of that edge. Traverse all edges to obtain a set of inscribed circle centers. Connect adjacent inscribed circle centers that satisfy the connectivity condition (the ratio of the distance between two inscribed circle centers to the sum of their radii is less than the connectivity threshold of 1.5) with straight line segments to form the initial skeleton network. The connectivity threshold of 1.5 is determined based on the following: when the distance between the centers does not exceed 1.5 times the sum of their radii, two inscribed circles overlap or are adjacent in space, and the skeleton branches can extend continuously.

[0040] The main framework lines are obtained through screening. For each branch segment of the framework line, the Euclidean distance between adjacent pixels (the square root of the sum of the squares of the differences in row coordinates and column coordinates of adjacent pixels) is accumulated pixel by pixel from one end to the other to obtain the length of the branch segment. The preset length threshold is 15 pixels (corresponding to an actual ground distance of approximately 7.5 meters). This value is based on the fact that in southern natural broad-leaved mixed forest areas, the short axis span of small-scale gaps formed by the wind-fallen single or several trees is usually within 7.5 meters. Short branches of the framework lines within such small gaps are unrelated to the main extension direction of typhoon-induced large-scale landslide gullies. Short branch segments with a length of less than 15 pixels are discarded. For each intersection node of the skeleton lines, calculate the included angle between each pair of branch segments that intersect there (take the absolute value of the difference between the two direction angles; if it is greater than 180°, subtract 360°). The preset angle threshold is 120°. The basis for this value is: when the included angle exceeds 120°, the two branches extend approximately in the same direction, and the shorter one is a redundant branch at the edge of the trench. Branches with included angles greater than 120° are discarded, and the shorter one is discarded. After the above two eliminations, the remaining branch segments are reconnected at the intersection node to obtain the main skeleton line.

[0041] Based on the main skeleton line, extract the endpoints and intersection nodes of the skeleton line to obtain the skeleton line topology node set. For each pixel on the main skeleton line, check the number of pixels belonging to the main skeleton line within its eight-neighborhood: if there is only one, it is an endpoint; if there are three or more, it is an intersection node. Summarize all endpoints and intersection nodes to obtain the skeleton line topology node set. Based on the skeleton line topology node set, calculate the azimuth angle of each branch segment extending from the intersection node to the endpoint, obtaining the skeleton line azimuth angle set containing the azimuth angles of each branch segment. For each branch segment, the first 5 consecutive pixels adjacent to the intersection node are selected, and a straight line is fitted using the least squares method. The angle between this straight line and true north is taken as the orientation azimuth (range 0° to 360°). The preset number of 5 pixels is based on the fact that at a spatial resolution of 0.5 meters, it corresponds to an actual ground extension distance of 2.5 meters, which is sufficient to stably reflect the initial spatial orientation of the branch after leaving the intersection node. Based on the set of orientation azimuths of the skeleton lines, branch segments with orientation azimuths pointing upwards on the terrain slope and whose azimuth values ​​are within the preset slope aspect range are selected to obtain the upslope pointing branches. The terrain slope aspect is extracted from the digital elevation model. The slope aspect is calculated using the third-order inverse distance squared weighted difference method. Taking the currently processed center pixel e as the reference, its eight neighboring pixels are arranged as follows: The neighboring pixel directly above e is denoted as b, the neighboring pixel directly below e is denoted as h, the neighboring pixel directly to the left is denoted as d, the neighboring pixel directly to the right is denoted as f, the diagonally opposite pixel to the upper left is denoted as a, the diagonally opposite pixel to the upper right is denoted as c, the diagonally opposite pixel to the lower left is denoted as g, and the diagonally opposite pixel to the lower right is denoted as i. The elevation value of each pixel is represented by the corresponding letter, in meters. The pixel spacing (i.e., the spatial resolution of the digital elevation model) is denoted as Δs. In this embodiment, Δs is taken as 0.5 meters, and the elevation change rate in the east-west direction is... The calculation formula is:

[0042] ;

[0043] In the formula: The rate of change of elevation in the east-west direction is dimensionless. , , These are the elevation values ​​of the three pixels in the left column, in meters. , , These are the elevation values ​​of three pixels in the right column, in meters. The spatial resolution of the digital elevation model is set to 0.5 meters. and Weighting is applied to adjacent pixels, with weight 2 reflecting that adjacent pixels contribute more than diagonal pixels; 8 is the weight normalization factor, representing the north-south elevation change rate. The calculation formula is:

[0044] ;

[0045] according to and Calculate slope aspect value :

[0046] ;

[0047] In the formula, This is the topographic slope value, in degrees, with due north as the baseline and increasing clockwise, ranging from 0° to 360°. It is a two-parameter arctangent function, according to and The symbol automatically determines the quadrant; for each skeleton line branch segment, the average slope aspect of all pixels it covers is taken as the representative value of the terrain slope aspect. The azimuth angle of the branch is compared with the representative value of the terrain slope aspect. When the absolute value of the difference between the two is less than the preset slope aspect deviation threshold of 30°, it is determined that the branch points upwards on the terrain slope. The preset slope aspect deviation threshold of 30° is based on the following: statistics of 156 landslide events cataloged in the study area show that the 95th percentile of the slope aspect deviation angle distribution of the main extension direction of the landslide gully is 28.6°, and 30° is a conservative rounding value.

[0048] When a branch points upwards on the terrain slope and its azimuth angle falls within the preset slope aspect range of 0° to 180°, the branch is marked as an upslope pointing branch. The preset slope aspect range of 0° to 180° is based on statistical data from the cataloging of typhoon-induced landslides in the study area over the past 20 years. 67.3% of the landslides were located on southeast to south-facing windward slopes. Based on the upslope pointing branches, adjacent branch segments on both sides of the upslope pointing branch with symmetrically diverging azimuth angles are extracted as flank constraint branches. The symmetrical divergence determination method is as follows: taking the azimuth angle of the upslope pointing branch as the axis of symmetry, the absolute value of the difference between the deviations of the azimuth angles of the left and right branches from the axis of symmetry is calculated. When this absolute value is less than the preset symmetry tolerance of 20°, it is determined to be a symmetrical divergence. The preset symmetry tolerance of 20° is based on the following: based on the measurement and statistics of the geometric characteristics of the confirmed landslide gully cross section in the study area, the asymmetric deviation of the V-shaped side wall development angle usually does not exceed 20°. A V-shaped geometric edge is constructed with the intersection node as the apex and the side lines (side constraint branches) converging towards the apex. The area enclosed by this V-shaped geometric edge is the interference candidate area.

[0049] Based on the interference candidate region, the pixels on both sides of the V-shaped geometric edge are extracted to obtain the edge pixel set. For the left and right sides of the V-shaped geometric edge, the row and column coordinates of all pixels are extracted pixel by pixel along the tracking direction from the tip vertex, and summarized into the edge pixel set. Based on the edge pixel set, the vertical direction of the local tangent at the pixels on both sides is calculated to obtain the edge normal direction pointing into the candidate region. Specifically, for each edge pixel, three adjacent edge pixels are taken in the forward and backward directions along the edge direction (a total of 7). A straight line is fitted using the least squares method, and the direction of the straight line rotated 90° clockwise is taken as the edge normal direction. The preset normal window radius of 3 pixels is based on the fact that at a spatial resolution of 0.5 meters, it corresponds to an actual ground distance of 1.5 meters. This window width can smooth the discrete sampling noise of the edge pixel coordinates. Based on the edge normal direction, the algorithm moves pixel by pixel along the edge normal direction into the interference candidate area at a preset step size (1 pixel) until it reaches the opposite side edge or reaches a preset maximum step distance of 50 pixels. The preset maximum step distance of 50 pixels is based on the fact that the maximum width of the head of a landslide gully induced by a typhoon in southern natural broad-leaved mixed forest areas is usually within 25 meters, corresponding to 50 pixels. The row and column coordinates of each step pixel are recorded sequentially along each normal direction to obtain the normal step pixel sequence.

[0050] Based on the normal step pixel sequence, the pixel-by-pixel change in the difference between the red and green band reflectance of adjacent pixels is calculated to obtain the red-green reflectance difference. The set of pixels where the red-green reflectance difference jumps from negative to positive and the absolute value of the change is greater than a preset jump threshold is extracted to obtain the mixed feature of slope shadow overlay and bare soil. Specifically, the calculation process is as follows: For two adjacent pixels in the normal step pixel sequence, reflectance values ​​are extracted from the red and green band surface reflectance images, respectively. For each pixel, the difference between the red band surface reflectance and the green band surface reflectance is calculated and denoted as the red-green reflectance difference. The difference between the red-green reflectance difference of the next pixel and the red-green reflectance difference of the previous pixel is calculated. This difference represents the change in red-green reflectance difference at the step position (pixel-by-pixel drop). The ground feature indication of the red-green reflectance difference is as follows: In the area of ​​the slope shadow formed by the V-shaped edge sidewall, the reflectance of the green band is higher than that of the red band, and the red-green reflectance difference is negative; In the area of ​​exposed bare soil, iron oxides in the soil have strong reflectance in the red band, and the reflectance of the red band is higher than that of the green band, and the red-green reflectance difference is positive.

[0051] As the surface transitions from shadow to bare soil along the normal direction from the edge inwards, the difference in red and green reflectance changes from negative to positive. When the absolute value of the change in red and green reflectance at a certain step location exceeds a preset threshold of 0.05, the pixel at that location is marked as a transition pixel between shadow and bare soil. The preset threshold of 0.05 (dimensionless reflectance difference) is determined as follows: Ten known landslide bare soil patches and adjacent shadow transition zones were selected within the study area. 200 sets of spectral samples were simultaneously collected at 0.5-meter intervals using an ASDFieldSpec spectroradiometer. The 5th percentile of the amplitude distribution of the change in red and green reflectance was 0.047, and the 10th percentile was 0.052. 0.05 corresponds to the boundary distinguishing between land cover type transformation and random noise fluctuations at approximately a 95% confidence level. Pixels marked as transition pixels between shadow and bare soil were extracted from all normal step pixel sequences, and the resulting feature of the slope shadow overlay with bare soil was obtained.

[0052] Based on the slope shadow overlay and bare soil mixed features obtained in step 2, the shortwave infrared reflectance within the interference candidate area is extracted to obtain the spatial distribution of shortwave infrared reflectance. The minimum row coordinate, maximum row coordinate, minimum column coordinate, and maximum column coordinate of all pixels in the slope shadow overlay and bare soil mixed features are used to construct the bounding rectangle of the extracted area. The shortwave infrared reflectance values ​​(dimensionless, ranging from 0 to 1) of all pixels within this rectangle are read from the radiometrically calibrated and atmospherically corrected shortwave infrared surface reflectance image, and organized into a two-dimensional grid data according to the pixel row and column coordinates to obtain the spatial distribution of shortwave infrared reflectance.

[0053] Based on the spatial distribution of shortwave infrared reflectance, the pixel-by-pixel difference in shortwave infrared reflectance is calculated along the extension direction of the main skeleton line to obtain a pixel-by-pixel fluctuation sequence of shortwave infrared reflectance. The specific process is as follows: Using the main skeleton line obtained in step 1 as a spatial reference line, each pixel is extended 5 pixels to each side perpendicular to the main skeleton line to form a sampling strip. The preset width of 5 pixels is based on the fact that the projected diameter of a single mature broad-leaved tree crown is approximately 8 to 12 meters (equivalent to 16 to 24 pixels). Extending 5 pixels (equivalent to 2.5 meters) to each side of the main skeleton line is sufficient to cover the staggered distribution of fallen logs and debris at the bottom of the trench. Within the sampling strip, along the extension direction, the average shortwave infrared reflectance of all pixels at each cross-section is taken as the representative value of shortwave infrared reflectance at that cross-section. The representative value of shortwave infrared reflectance at each subsequent cross-section is subtracted from the representative value of the previous cross-section to obtain the pixel-by-pixel fluctuation sequence of shortwave infrared reflectance.

[0054] Based on the pixel-by-pixel fluctuation sequence of shortwave infrared reflectance, continuous pixel segments in which the shortwave infrared reflectance exhibits alternating positive and negative changes and the alternation period is within a preset period range are extracted to obtain the micro-topographic features of alternating fallen trees and dead debris. Specifically, the process involves scanning element by element from the beginning of the wave sequence, recording the position where the symbol of each element in the sequence is reversed, and recording each reversal as an alternation event. The pixel span between two adjacent alternation events is recorded as the alternation period (physically meaning the spatial combination length of a section of fallen logs and a section of debris accumulation along the direction of the gully). The preset period interval is 3 to 15 pixels (corresponding to actual ground distances of 1.5 meters to 7.5 meters). This value is based on the statistical analysis of the axial length of 168 groups of fallen log and debris spatial combinations measured in the field at 12 confirmed landslide accumulation areas within the study area. The average value is 4.2 meters (equivalent to 8.4 pixels), the standard deviation is 1.8 meters (equivalent to 3.6 pixels), the 5th percentile is 1.8 meters (equivalent to 3.6 pixels), and the 95th percentile is 7.1 meters (equivalent to 14.2 pixels). The preset period interval covers 90% of the distribution range. When a continuous pixel segment in a wave sequence satisfies both of the following conditions, the continuous pixel segment is marked as a micro-topographic feature of alternating fallen trees and dead debris: (1) at least 3 positive and negative alternation events occur within the pixel segment (corresponding to at least 2 pairs of complete alternation sequences, which can reliably distinguish between regular textures and random noise); (2) the alternation period of each alternation event falls within 3 to 15 pixels. The spatial range of the continuous pixel segment that meets the conditions is mapped back to the image space to obtain the micro-topographic feature of alternating fallen trees and dead debris.

[0055] Based on the main framework line and the alternating micro-topographic features of fallen trees and fallen debris, the morphological boundaries of the trenches are determined, and the center lines of the alternating micro-topographic features are extracted to obtain the micro-topographic centerline. For the binary spatial range corresponding to each spatial unit of the alternating micro-topographic features of fallen trees and fallen debris, a morphological thinning algorithm is used to extract the centerline. The thinning process is as follows: a set of preset structural element templates are repeatedly applied, and a hit-and-miss transformation is performed on each pixel with a value of 1 (when the value corresponding to the hit position is 1 and the value corresponding to the miss position is 0, the pixel is marked as peelable). At the end of one round, all marked boundary pixels are set from 1 to 0. This process is repeated until no pixels can be further peeled. The remaining pixels with a value of 1 constitute the micro-topographic centerline with a width of 1 pixel. The beginning and end endpoints of the centerlines of adjacent spatial units are connected by straight line segments when the Euclidean distance is within 10 pixels of the preset connection distance. The preset connection distance of 10 pixels is based on the fact that the axial distance between two adjacent alternating areas of fallen trees and fallen debris is usually within 5 meters.

[0056] Calculate the deviation angle between the micro-topography centerline and the main framework line. Extract the local tangent direction at each position on the micro-topography centerline, and simultaneously extract the corresponding local tangent direction on the main framework line. Calculate the absolute value of the angle between the two, which is the deviation angle. Extract continuous line segments with a deviation angle less than a preset deviation threshold of 25° to obtain the ditch morphology boundary. The preset deviation threshold of 25° is determined as follows: when the deviation between the micro-topography centerline and the main framework line exceeds 25°, the micro-texture arrangement direction of fallen logs and debris deviates from the main ditch direction, indicating that these fallen logs were formed by non-landslide factors. For each continuous line segment with a deviation angle less than 25°, step outwards vertically from both ends of the line segment. When the shortwave infrared reflectance at several consecutive step positions no longer exhibits an alternating high and low pattern, it is determined that the area has exceeded the outer edge of the micro-topography feature space. Connect the outer edge points extending to both sides from the start and end points of the line segment to form a closed polygon. The strip-shaped space enclosed by this polygon is the ditch morphology boundary.

[0057] Based on the trench's morphological boundary, pixels inside the boundary are extracted to obtain a candidate pixel set. A ray-mapping method is used for judgment: a ray is emitted from each pixel in a due east direction, and the number of intersections between this ray and each side of the polygon is counted. An odd number of intersections indicates the pixel is inside the polygon, and an even number indicates it is outside. Internal pixels are recorded in the candidate pixel set. Based on the candidate pixel set, the vertical distance of each pixel from the centerline of the trench's morphological boundary is calculated to obtain the pixel spatial distribution distance. The centerline is extracted by taking the average coordinates of pixels at the same extension position on both edge lines to obtain the centerline pixel. Connecting all centerline pixels forms the centerline. For each candidate pixel, the minimum Euclidean distance to all pixels on the centerline is calculated to obtain the pixel spatial distribution distance. Based on the pixel spatial distribution distance, pixels with a vertical distance less than a preset distance limit of 10 pixels are selected to obtain the core interference pixel subset. The preset distance limit of 10 pixels is based on the following criteria: the area within 5 meters on both sides of the center line of the ditch is the core disturbance area of ​​the landslide, and pixels beyond this distance belong to the canopy area above the sidewall of the ditch that has not been directly disturbed by the landslide.

[0058] Based on the core interference pixel subset and combined with the alternating high and low reflectance striped texture of shortwave infrared imagery (i.e., the micro-topographic feature of fallen trees superimposed with debris), the core interference pixel subset is marked as the interference pixel subset, resulting in the interference type classification. The classification method is as follows: check whether each core pixel is located within the spatial range of the alternating micro-topographic feature of fallen trees superimposed with debris. Those within the range are marked as landslide interference pixels, while those within the core interference pixel subset but not within the range of the alternating micro-topographic feature are marked as wind-fall interference pixels. The classification is based on the following: the typical surface feature of landslide interference is that fallen trees are arranged in a regular alternating pattern along the direction of the gully, appearing as a striped texture of alternating high and low reflectance on shortwave infrared images; wind-fall interference is mainly characterized by concentrated lodging, and does not show an alternating striped texture on shortwave infrared images. The core interference pixel subset that has been classified and marked is used as the interference pixel subset. Each pixel carries three attributes: row coordinates, column coordinates, and interference category label, resulting in the interference type classification.

[0059] Based on the interference type classification results, the shortwave infrared reflectance of the interference pixel subset at each temporal phase is extracted to obtain a multi-temporal shortwave infrared reflectance sequence. All available multi-temporal satellite images of the study area over the past 10 years (with revisit cycles of 5 to 16 days) are acquired, and preprocessing is performed on each temporal image according to radiometric calibration and atmospheric correction. From the shortwave infrared surface reflectance image of each temporal image, the shortwave infrared reflectance values ​​of all pixels in the interference pixel subset are extracted, and their average value is calculated as the representative value of the shortwave infrared reflectance of that interference pixel subset at that temporal phase. The temporal phases are arranged chronologically according to imaging time, and the corresponding representative values ​​are arranged in the same chronological order to obtain a multi-temporal shortwave infrared reflectance sequence. Based on the multi-temporal shortwave infrared reflectance sequence, the difference in shortwave infrared reflectance between two adjacent temporal phases is calculated to obtain the change in shortwave infrared reflectance in a single period. The calculation method is: subtract the representative value of the previous temporal phase from the representative value of the later temporal phase, iterating through the data. At each time phase, we obtain Based on the single-period shortwave infrared reflectance variation values, the cumulative canopy disturbance recovery is obtained by sequentially accumulating the shortwave infrared reflectance variation values ​​along the time axis. The accumulation method is as follows: the first time phase on the time axis is set as the baseline time phase, and its cumulative canopy disturbance recovery value is zero. For the second time phase... Each phase ( The cumulative recovery of its canopy disturbance is equal to that of the first to the second period. Total period The algebraic sum of the changes in shortwave infrared reflectance in a single period. The unit of the cumulative recovery of canopy disturbance is the dimensionless cumulative reflectance, which physically means the total net cumulative change in shortwave infrared reflectance in the disturbed area from the reference time phase to the current time phase.

[0060] Based on the cumulative recovery of canopy disturbance, the difference in cumulative recovery between adjacent time phases is calculated to obtain the periodic variation value of the cumulative amount. The calculation method is as follows: subtract the cumulative value of the previous time phase from the cumulative value of the later time phase. Based on the periodic variation value of the cumulative amount, time phases with periodic variation values ​​less than a preset variation limit are selected to obtain the time phases where vegetation recovery has stagnated. The preset variation limit is set at 0.3 times the standard deviation of the cumulative periodic variation value series. The basis for this value is that after vegetation spectral recovery enters a stable stage, the periodic fluctuations mainly come from seasonal phenological changes and random sensor measurement noise, and the amplitude of the periodic variation value usually does not exceed 30% of the overall standard deviation of the entire recovery period. Based on the time series analysis of shortwave infrared reflectance of multiple confirmed restored landslide patches in the study area, the average ratio of the standard deviation of the change value during the stable restoration phase to the standard deviation of the complete restoration phase is 0.27, leaving a margin of 0.3 times the standard deviation. According to the vegetation restoration stagnation phase, the preceding phase is extracted to obtain the critical phase for vegetation spectral stabilization. Combining this critical phase with the start date of the typhoon season in natural broadleaf mixed forests (June 1st each year), a time interval containing the critical phase and the start date of the typhoon season is generated, thus obtaining the disturbance occurrence time window. The basis for choosing June 1st as the start date of the typhoon season is: according to statistics from the China Meteorological Administration's Tropical Cyclone Yearbook, from 1949 to 2023, typhoons forming in June and subsequent months accounted for 87.3% of the total number of typhoons affecting the southern coastal areas of China (18°N to 25°N, 108°E to 122°E).

[0061] Based on the time window of disturbance occurrence, all subsets of disturbance pixels within the time interval are extracted to obtain a temporal disturbance pixel set. For each subset of disturbance pixels in the disturbance type classification results, its corresponding disturbance occurrence time window is obtained and arranged in chronological order. For multiple subsets of disturbance pixels with overlapping time windows, the endpoint pixels of the gully direction morphological boundary closest to the slope top are arranged along the slope surface from the slope top to the slope bottom as a secondary sorting criterion. Based on the temporal disturbance pixel set, the temporal disturbance pixel set is sorted according to the chronological order of the disturbance occurrence time window to obtain a temporally arranged disturbance pixel set. Based on the temporal arranged disturbance pixel set and combined with the disturbance type classification results, a corresponding disturbance type label (landslide disturbance or wind-blown disturbance) and occurrence year (the lower boundary of the disturbance occurrence time window, i.e., the year of the start of the typhoon season) are assigned to each temporal arranged disturbance pixel set. The disturbance type label, disturbance occurrence year, gully direction morphological boundary polygon vertex coordinate sequence, and canopy disturbance recovery cumulative amount sequence are arranged in chronological order of the disturbance occurrence year to generate a natural forest disturbance history reconstruction sequence and obtain the reconstruction result.

[0062] Satellite imagery data from the past 10 years in a natural broadleaf mixed forest area in southern China was processed to reconstruct the historical sequence of natural forest disturbances in the region. The reconstruction results show that this method can automatically identify typhoon-induced landslide disturbance areas and wind-blown disturbance areas from multi-temporal optical satellite imagery, and invert the year of occurrence of each disturbance event through the time series changes of shortwave infrared reflectance. Compared with traditional field plot survey methods, the spatial coverage of this method extends from the transect reachable area to the entire satellite imagery coverage area, the continuity of temporal reconstruction is not limited by the number of surviving standing tree samples, and the determination of disturbance type does not depend on the subjective experience differences of the interpreters.

[0063] like Figure 2 As shown, embodiments of the present invention also provide a natural forest disturbance history reconstruction system based on satellite image analysis, comprising:

[0064] The extraction module is used to extract the canopy gap boundary from multi-temporal optical satellite images of natural broad-leaved mixed forest areas to obtain the gap outline; based on the gap outline, the boundary curvature extreme points are calculated to construct an initial polygon, the inscribed circle center line is extracted to generate the skeleton line, and the branch length and included angle of the skeleton line are calculated to obtain the main skeleton line;

[0065] The identification module is used to identify the V-shaped geometric edge with its tip pointing upwards on the slope based on the main skeleton line, and obtain the interference candidate area; based on the interference candidate area, the pixel-by-pixel drop amplitude of the red and green band reflection difference along the normal direction of the V-shaped edge is calculated to obtain the mixed features of slope shadow and bare soil.

[0066] The module is used to extract the striped texture with alternating high and low shortwave infrared reflectance based on the mixed characteristics of slope shadow and bare soil, and obtain the micro-topographic features of fallen trees and fallen debris; the morphological boundary of the ditch is determined based on the main framework line and the micro-topographic features of fallen trees and fallen debris.

[0067] The segmentation module is used to divide the candidate area of ​​disturbance into a subset of disturbance pixels according to the morphological boundary of the ditch, and obtain the disturbance type classification result; based on the disturbance type classification result, the periodic change of shortwave infrared reflectance is accumulated in the subset of disturbance pixels to obtain the cumulative amount of canopy disturbance recovery.

[0068] The mapping module is used to determine the critical time phase for vegetation spectral recovery stability based on the cumulative amount of canopy disturbance recovery, and obtain the disturbance occurrence time window; based on the disturbance occurrence time window, the disturbance type classification results are mapped to generate a natural forest disturbance history reconstruction sequence, and the reconstruction results are obtained.

[0069] It should be noted that this system is a system corresponding to the above method. All implementation methods in the above method embodiments are applicable to this embodiment and can achieve the same technical effect.

[0070] Embodiments of the present invention also provide a computing device, including: a processor and a memory storing a computer program, wherein the computer program, when executed by the processor, performs the method described above. All implementations in the above method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0071] Embodiments of the present invention also provide a computer-readable storage medium storing instructions that, when executed on a computer, cause the computer to perform the method described above. All implementations in the above method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0072] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for reconstructing the history of disturbance in natural forests based on satellite image analysis, characterized in that, The method includes: Based on multi-temporal optical satellite images of natural broad-leaved mixed forest areas, the canopy gap boundaries are extracted to obtain the gap contours; based on the gap contours, the extreme points of boundary curvature are calculated to construct initial polygons, the lines connecting the centers of the inscribed circles are extracted to generate skeleton lines, and the branch lengths and included angles of the skeleton lines are calculated to obtain the main skeleton lines; Based on the main skeleton line, identify the V-shaped geometric edge with its tip pointing upwards on the slope to obtain the interference candidate area; based on the interference candidate area, calculate the pixel-by-pixel drop amplitude of the red and green band reflection difference along the normal direction of the V-shaped edge to obtain the mixed features of slope shadow and bare soil. Based on the mixed characteristics of slope shadow and bare soil, the striped texture with alternating high and low short-wave infrared reflectance is extracted to obtain the micro-topographic features of fallen trees and fallen debris; based on the main framework line and the micro-topographic features of fallen trees and fallen debris, the morphological boundary of the ditch is determined. Based on the morphological boundaries of the gullies, the candidate areas of disturbance are divided into subsets of disturbance pixels to obtain the disturbance type classification results. Based on the disturbance type classification results, the periodic changes in shortwave infrared reflectance are accumulated in the subsets of disturbance pixels to obtain the cumulative amount of canopy disturbance recovery. Based on the cumulative amount of canopy disturbance recovery, the critical time phase for vegetation spectral recovery and stability is determined, and the disturbance occurrence time window is obtained. Based on the disturbance occurrence time window, the disturbance type classification results are mapped to generate a natural forest disturbance history reconstruction sequence, and the reconstruction results are obtained.

2. The method for reconstructing the history of disturbance in natural forests based on satellite image analysis according to claim 1, characterized in that, Based on multi-temporal optical satellite imagery of natural broadleaf mixed forest areas, canopy gap boundaries were extracted to obtain gap contours. Based on these gap contours, extreme points of boundary curvature were calculated to construct initial polygons. The centers of the inscribed circles were extracted to generate skeleton lines. The branch lengths and included angles of these skeleton lines were calculated to obtain the main skeleton lines, including: Based on multi-temporal high-resolution optical satellite images of natural broad-leaved mixed forest areas, the canopy gap boundaries were extracted to obtain the gap contours; Based on the gap profile, calculate the tangent direction deflection angle of adjacent pixels along the boundary direction to obtain the set of boundary curvature extreme points; Based on the set of extreme points of boundary curvature, connect adjacent extreme points to construct an initial polygon, and obtain the inscribed polygon of the gap. Based on the inscribed polygons in the gaps, the initial skeleton line network is generated by connecting the centers of the inscribed circles on each side. The length of each branch line segment and the intersection angle are calculated to remove short branch line segments with a length less than a preset length threshold and obtuse branch line segments with an angle greater than a preset angle threshold, thus obtaining the main skeleton line.

3. The method for reconstructing the history of disturbance in natural forests based on satellite image analysis according to claim 2, characterized in that, Based on the main skeleton line, identify the V-shaped geometric edge with its tip pointing upwards on the slope to obtain interference candidate regions, including: Based on the main skeleton line, extract the nodes where the endpoints of the skeleton line intersect with the skeleton line to obtain the set of topological nodes of the skeleton line; Based on the set of topological nodes of the skeleton line, calculate the azimuth angle of each branch line segment extending from the intersection node to the endpoint, and obtain the set of skeleton line azimuth angles containing the azimuth angles of each branch line segment. Based on the set of azimuth angles of the skeleton lines, filter the branch lines whose azimuth angles point upwards on the terrain slope and whose azimuth angle values ​​are within the preset slope direction range to obtain the upslope pointing branches. Based on the uphill pointing branch, the adjacent branch segments on both sides of the uphill pointing branch with symmetrical divergence of azimuth angle are extracted as flank constraint branches. A V-shaped geometric edge is constructed with the intersection node as the apex and the side lines converge toward the apex, thus obtaining the interference candidate region.

4. The method for reconstructing the history of disturbance in natural forests based on satellite image analysis according to claim 3, characterized in that, Based on the interference candidate region, the pixel-by-pixel drop in the red-green band reflection difference along the V-shaped edge normal direction is calculated to obtain the mixed characteristics of slope shadow and bare soil, including: Based on the interference candidate region, extract the edge pixel points of both sides of the V-shaped geometric edge to obtain the edge pixel point set; Based on the set of edge pixels, calculate the vertical direction of the local tangent at the pixels on both sides of the edge to obtain the edge normal direction pointing into the candidate region. Based on the edge normal direction, the pixel sequence of normal stepping is obtained by stepping pixel by pixel along the edge normal direction into the interference candidate area at a preset step size. Based on the normal step pixel sequence, the pixel-by-pixel change in the difference between the red band reflectance and the green band reflectance of adjacent pixels is calculated to obtain the red-green reflectance difference. The set of pixels where the red-green reflectance difference jumps from negative to positive and the absolute value of the change is greater than the preset jump threshold is extracted to obtain the slope shadow superimposed bare soil mixed feature.

5. The method for reconstructing the history of disturbance in natural forests based on satellite image analysis according to claim 4, characterized in that, Based on the mixed characteristics of slope shadow and bare soil, the striped texture with alternating high and low short-wave infrared reflectance was extracted to obtain the micro-topographic features of fallen trees and deadwood. Based on the main framework line and the alternating micro-topographic features of fallen trees and debris, the morphological boundaries of the ditch are determined, including: Based on the characteristics of the mixed bare soil and the superimposed shadow on the slope, the short-wave infrared reflectance within the candidate interference area is extracted to obtain the spatial distribution of short-wave infrared reflectance. Based on the spatial distribution of short-wave infrared reflectance, the pixel-by-pixel difference of short-wave infrared reflectance is calculated along the extension direction of the main skeleton line to obtain the pixel-by-pixel fluctuation sequence of short-wave infrared reflectance. Based on the pixel-by-pixel fluctuation sequence of shortwave infrared reflectance, continuous pixel segments in which the shortwave infrared reflectance shows alternating positive and negative changes and the alternation period is within a preset period range are extracted to obtain the micro-topographic features of alternating fallen trees and dead debris. Based on the alternating micro-topographic features of fallen trees and dead debris, the centerline of the alternating micro-topographic features is extracted to obtain the micro-topographic centerline. The deviation angle between the micro-topographic centerline and the main framework line is calculated, and continuous line segments with a deviation angle less than a preset deviation threshold are extracted to obtain the morphological boundary of the trench.

6. The method for reconstructing the history of disturbance in natural forests based on satellite image analysis according to claim 5, characterized in that, Based on the morphological boundaries of the trenches, the interference candidate region is divided into a subset of interference pixels, resulting in an interference type classification, including: Based on the groove orientation and shape boundary, extract the pixels inside the boundary to obtain a candidate pixel set; based on the candidate pixel set, calculate the vertical distance of each pixel from the center line of the groove orientation and shape boundary to obtain the pixel spatial distribution distance. Based on the spatial distribution distance of pixels, pixels with a vertical distance less than a preset distance limit are selected to obtain a subset of core interference pixels. Based on the core interference pixel subset combined with the striped texture of alternating high and low shortwave infrared reflectance, the core interference pixel subset is marked as the interference pixel subset, and the interference type classification result is obtained.

7. The method for reconstructing the history of disturbance in natural forests based on satellite image analysis according to claim 6, characterized in that, Based on the classification results of interference types, the periodic changes in shortwave infrared reflectance are accumulated within the subset of interfering pixels to obtain the cumulative amount of canopy disturbance recovery, including: Based on the classification results of interference types, the shortwave infrared reflectance of the interference pixel subsets under each time phase is extracted to obtain a multi-time phase shortwave infrared reflectance sequence. Based on the multi-temporal shortwave infrared reflectance sequence, the difference in shortwave infrared reflectance between two adjacent temporal phases is calculated to obtain the change value of shortwave infrared reflectance in a single period. Based on the changes in shortwave infrared reflectance in a single period, the changes in shortwave infrared reflectance in each period are accumulated sequentially along the time axis to obtain the cumulative amount of canopy disturbance recovery.

8. The method for reconstructing the history of disturbance in natural forests based on satellite image analysis according to claim 7, characterized in that, Based on the cumulative amount of canopy disturbance recovery, the critical time phase for vegetation spectral stabilization is determined, and the time window for disturbance occurrence is obtained, including: Based on the cumulative recovery of canopy disturbance, the difference in the cumulative recovery of canopy disturbance between adjacent time phases is calculated to obtain the periodic variation value of the cumulative amount; Based on the cumulative change value over a period of time, the time phases in which the cumulative change value over a period of time is less than the preset change limit are selected to obtain the time phases in which vegetation restoration stagnates. Based on the stagnation phase of vegetation recovery, the preceding phase of vegetation recovery stagnation phase is extracted to obtain the critical phase for stable vegetation spectral recovery. By combining the critical time for vegetation spectrum stabilization with the start time of the typhoon-prone season in natural broadleaf mixed forests, a time interval containing the critical phase and the start time of the typhoon-prone season is generated, thus obtaining the time window for disturbance occurrence.

9. The method for reconstructing the history of disturbance in natural forests based on satellite image analysis according to claim 8, characterized in that, Based on the time window of disturbance occurrence, the disturbance type classification results are mapped to generate a historical reconstruction sequence of natural forest disturbance, resulting in the reconstruction results, including: Based on the time window of interference occurrence, extract all subsets of interference pixels within the time interval to obtain the temporal interference pixel set; Based on the temporal interference pixel set, the temporal interference pixel set is sorted according to the order of the interference occurrence time windows to obtain the temporal arrangement interference pixel set. Based on the temporal arrangement of the interference pixel set and the interference type classification results, each temporal arrangement of the interference pixel set is assigned a corresponding interference type label and the year of occurrence, thus obtaining the historical reconstruction sequence of natural forest interference.

10. A natural forest disturbance history reconstruction system based on satellite image analysis, wherein the system implements the method as described in any one of claims 1 to 9, characterized in that, include: The extraction module is used to extract the canopy gap boundary and obtain the gap outline based on multi-temporal optical satellite images of natural broad-leaved mixed forest areas; Based on the gap profile, calculate the extreme points of boundary curvature to construct an initial polygon, extract the line connecting the centers of the inscribed circles to generate the skeleton line, calculate the branch length and included angle of the skeleton line to obtain the main skeleton line; The identification module is used to identify the V-shaped geometric edge with its tip pointing upwards on the slope based on the main skeleton line, and obtain the interference candidate area; based on the interference candidate area, the pixel-by-pixel drop amplitude of the red and green band reflection difference along the normal direction of the V-shaped edge is calculated to obtain the mixed features of slope shadow and bare soil. The module is used to extract the striped texture with alternating high and low shortwave infrared reflectance based on the mixed characteristics of slope shadow and bare soil, and obtain the micro-topographic features of fallen trees and fallen debris; the morphological boundary of the ditch is determined based on the main framework line and the micro-topographic features of fallen trees and fallen debris. The segmentation module is used to divide the candidate area of ​​disturbance into a subset of disturbance pixels according to the morphological boundary of the ditch, and obtain the disturbance type classification result; based on the disturbance type classification result, the periodic change of shortwave infrared reflectance is accumulated in the subset of disturbance pixels to obtain the cumulative amount of canopy disturbance recovery. The mapping module is used to determine the critical time phase for vegetation spectral recovery stability based on the cumulative amount of canopy disturbance recovery, and obtain the disturbance occurrence time window; based on the disturbance occurrence time window, the disturbance type classification results are mapped to generate a natural forest disturbance history reconstruction sequence, and the reconstruction results are obtained.

Images