Forestry ecosystem health assessment method based on multi-source indexes

Abstract

The invention relates to the technical field of ecological environment monitoring, in particular to a forestry ecosystem health assessment method based on multi-source indexes. The method comprises the following steps: obtaining forest infrared data and forest point cloud data; performing crown energy flux calculation according to the forest infrared data to obtain crown energy flux data; constructing an earth surface roughness model according to the forest region point cloud data to obtain an earth surface roughness model; performing surface energy reflection processing according to the surface roughness model to obtain surface energy reflection data; performing double-domain energy gradient processing according to the crown energy flux data and the earth surface energy reflection data to obtain double-domain energy gradient data; and according to the double-domain energy gradient data, ecological sub-health region identification is carried out to obtain ecological sub-health region data. According to the method, the crown energy flux and the earth surface reflection features are fused, the vertical energy conduction abnormal region in the ecological system is identified, and the space identification capability and the discrimination precision of the forestry ecological sub-health region are improved.

Description

Technical Field

[0001] This invention relates to the field of ecological environment monitoring technology, and in particular to a method for assessing the health of forestry ecosystems based on multi-source indicators. Background Technology

[0002] With the intensification of global climate change and the continuous increase in the frequency and intensity of human activities, forest ecosystems, as a key support for terrestrial ecological security, are facing increasingly severe challenges to their structural stability, functional integrity, and energy cycle coordination. Multiple factors, including changes in hydrothermal patterns, extreme climate events, land use alteration, and local disturbances (such as logging, pests and diseases, and road construction), are accelerating the heterogeneous evolution and degradation of ecosystems. Traditional forestry ecological health assessment methods are mostly based on single-point-in-time vegetation indices (such as NDVI and EVI) or static land cover information for modeling, often failing to capture micro-scale disturbances and lacking the ability to dynamically express energy flow mechanisms and spatial structural evolution. Summary of the Invention

[0003] To address the aforementioned technical problems, this invention proposes a forestry ecosystem health assessment method based on multi-source indicators, thereby resolving at least one of the aforementioned technical issues.

[0004] This application provides a method for assessing the health of forest ecosystems based on multi-source indicators, including the following steps: Step S1: Acquire infrared data and point cloud data of the forest area; calculate the canopy energy flux based on the infrared data of the forest area to obtain the canopy energy flux data; Step S2: Construct a surface roughness model based on forest area point cloud data to obtain a surface roughness model; perform surface energy albedo processing based on the surface roughness model to obtain surface energy albedo data; Step S3: Perform dual-domain energy gradient processing based on canopy energy flux data and surface energy albedo data to obtain dual-domain energy gradient data; Step S4: Identify ecologically sub-healthy regions based on dual-domain energy gradient data to obtain ecologically sub-healthy region data.

[0005] This invention achieves layered modeling of the energy structure characteristics of the forest canopy and surface layers through multi-source data fusion of infrared thermal imaging and lidar point clouds. The acquisition of canopy energy flux reflects vegetation evapotranspiration and ecological metabolic levels, while the surface roughness model characterizes micro-topographic disturbance and structural fragmentation. Furthermore, surface energy feedback characteristics are formed through coupled albedo processing. The construction of dual-domain energy gradients not only reveals the vertical energy coordination between the canopy and the surface but can also be used to identify potential ecological fault zones or areas of abnormal energy flow. This invention can accurately identify ecologically unhealthy areas without relying on manual annotation, enabling ecological health early warning based on energy characteristics and improving the interpretability and predictive capabilities of forestry remote sensing monitoring.

[0006] Preferably, the calculation of canopy capacity flux is as follows: Radiation brightness temperature correction was performed based on infrared data from the forest area to obtain canopy surface temperature data; Net radiation data of the canopy is obtained by calculating the net radiation based on the canopy surface temperature data. Meteorological monitoring data was acquired, and sensible heat flux was estimated based on canopy net radiation data and meteorological monitoring data to obtain canopy sensible heat flux data; Latent heat flux was estimated based on canopy net radiation data and canopy sensible heat flux data to obtain canopy evapotranspiration heat dissipation flux data; By integrating the canopy sensible heat flux data and the canopy evapotranspiration heat flux data, the tree canopy energy flux data is obtained.

[0007] The canopy energy flux calculation described in this invention obtains the canopy surface temperature by performing radiation brightness temperature correction on infrared remote sensing data. This, combined with meteorological monitoring data, constructs a canopy energy balance model to systematically estimate key parameters such as net radiation, sensible heat flux, and latent heat flux. Compared to traditional remote sensing indicators that rely on a single vegetation index, this thermodynamic approach more directly reflects vegetation evapotranspiration capacity and energy exchange intensity, offering higher physical interpretability and dynamic response capabilities. By integrating sensible and latent heat fluxes into canopy energy flux data, a comprehensive characterization of vegetation physiological state and ecological metabolic activity can be achieved.

[0008] Preferably, the surface roughness model is constructed as follows: An elevation raster model was constructed based on forest area point cloud data to obtain the elevation raster model. Local feature extraction is performed on the elevation raster model to obtain local feature data; Roughness calculation is performed on local feature data to obtain roughness data; The surface scattering characteristic parameters are calculated based on the roughness data to obtain the surface roughness model.

[0009] This invention constructs a high-precision elevation raster model from forest point cloud data and extracts local undulation, aspect variation, and structural fragmentation features to achieve multi-dimensional perception of surface micro-topographic disturbances. Through roughness calculation, it not only captures traditional elevation fluctuation information but also integrates structural disturbance factors such as directional distortion and spatial connectivity, thereby generating roughness data with spatial texture and structural sensitivity features. By coupling roughness with energy reflection mechanisms, surface scattering characteristic parameters are obtained, effectively improving the physical realism and spatial resolution of surface energy response modeling.

[0010] Preferably, local feature extraction specifically involves: Local relief is calculated for the elevation grid model to obtain local relief data; The slope direction variation data is obtained by calculating the slope direction variation based on the local undulation data; Based on the slope cutting direction variation data, a patch structure connectivity map is constructed to obtain patch structure connectivity map data; The patch structure fragmentation index is calculated based on the patch structure connectivity graph data to obtain the patch structure fragmentation index data. Local feature data are obtained by integrating local undulation data, slope direction variation data, patch structure connectivity map data, and patch structure fragmentation index data.

[0011] This invention effectively overcomes the limitations of traditional terrain parameter extraction, which relies solely on undulation, by extracting local features at multiple scales and dimensions from an elevation raster model. By combining four types of information—local undulation, slope aspect distortion, patch structural connectivity, and fragmentation index—a structural disturbance-sensitive feature system is constructed, capable of precisely revealing the micro-disturbance patterns and spatial heterogeneity evolution characteristics of forest surfaces. This invention can identify potential landslide zones, wind erosion path potential zones, and shrub degradation patches, demonstrating good adaptability and discrimination capabilities against ecological disturbances.

[0012] Preferably, the roughness calculation is as follows: Local elevation relief layers are constructed based on local feature data and elevation raster models to obtain elevation relief map data. The aspect distortion layer is extracted from the elevation raster model to obtain the aspect distortion layer data; Based on local feature data, elevation relief map data, and slope aspect distortion layer data, patch structure fragmentation layer is extracted to obtain patch structure data. Roughness data is obtained by performing roughness weighting on the patch structure data.

[0013] This invention achieves roughness modeling with directional, structural, and scale-aware capabilities by utilizing multi-source disturbance indices such as elevation undulation, aspect distortion, and patch structure fragmentation during roughness calculation. It identifies areas of significant topographic fluctuation by identifying local elevation changes, characterizes the discontinuity of disturbance direction by combining aspect statistical dispersion, and reflects the integrity of spatial organization and ecological connectivity by using patch structure fragmentation, effectively breaking through the traditional roughness expression method centered on elevation standard deviation. By weighting roughness and fusing multi-dimensional features, the resulting roughness data not only enhances the sensitivity of surface geometry but also improves the overall quality of the roughness data.

[0014] Preferably, the extraction of the slope aspect distortion layer specifically involves: Slope aspect layer data is obtained by extracting the slope aspect layer from the elevation raster model. Calculate the neighboring slope aspect value from the slope aspect layer data to obtain the neighboring slope aspect value data; Local aspect dispersion data is obtained by calculating the local aspect dispersion based on the neighboring aspect value data. Based on the local dispersion data of slope aspect, the flat areas of the slope aspect layer data are removed to obtain the slope aspect distortion layer data.

[0015] In this invention, during the extraction of the aspect distortion layer, a discrete metric method for local aspect disturbances is constructed by performing neighborhood statistical analysis on the aspect layer generated from the elevation raster model. This method effectively identifies the complexity and discontinuity of the slope direction distribution. A circular statistical method is used to calculate the dispersion of neighborhood aspect values, avoiding the jump errors that occur when traditional methods handle 360-degree circumferential angles. By removing flat edge areas and filtering regions with minimal slope variation or consistent direction, the extracted distortion index is ensured to have actual disturbance directionality.

[0016] Preferably, the surface energy albedo treatment specifically includes: Obtain surface spectral reflectance data; The surface roughness model and surface spectral reflectance data are registered to obtain joint surface data. Roughness reflectance coupling processing is performed on the surface joint data to obtain surface coupled data; The surface coupling data is processed by albedo perturbation mapping to obtain surface perturbation energy data; Surface energy albedo is obtained by performing surface energy perturbation data.

[0017] This invention achieves precise control over surface albedo characteristics by constructing a surface roughness model and spatially registering it with surface spectral reflectance. Compared to traditional albedo estimation methods that rely solely on spectral information, this method uses surface roughness as a structural perturbation factor and performs coupled modeling based on joint data, effectively improving the response to energy scattering inhomogeneities. Through albedo perturbation map processing, areas of anomalous albedo caused by microscale topographic undulations, structural fragmentation, or bare land exposure can be identified. The generated surface energy albedo data not only preserves spectral physical characteristics but also reflects the spatial variability of actual energy reflection behavior.

[0018] Preferably, step S3 specifically includes: Vertical energy channel data are obtained by spatial projection alignment based on canopy energy flux data and surface energy albedo data. The energy difference between the vertical energy channels is calculated in two domains to obtain the energy difference data between the two domains. Gradient fields are constructed based on the energy difference data in the two domains to obtain the energy gradient data in the two domains.

[0019] In this invention, during the dual-domain energy gradient processing, a vertically consistent energy channel data structure is constructed by aligning the spatial projections of canopy energy flux data and surface energy albedo data, ensuring the spatial accuracy of energy conduction analysis between the canopy and the surface. Through dual-domain energy difference calculation, energy imbalances between different ecological levels in forest areas can be effectively revealed, such as the coupled anomalies of weakened evapotranspiration and enhanced surface scattering. The constructed energy gradient field not only reflects the spatial abrupt changes in energy distribution but also has the ability to identify ecological fault zones, areas of abnormal energy exchange, and potential degraded patches. This invention breaks through the limitations of traditional ecological health assessments based solely on single-layer indicators, providing a technical foundation for monitoring the vertical functional coordination of ecosystems through canopy-surface dual-domain synergistic physical field modeling.

[0020] Preferably, step S4 specifically includes: A two-domain energy gradient map is constructed based on the two-domain energy gradient data to obtain the two-domain energy gradient map data. Local anomaly detection is performed on the dual-domain energy gradient map data to obtain local anomaly data; Gradient continuity disruption regions are identified in local anomaly data to obtain disrupted region data; Based on the data of the damaged areas, ecological sub-health cluster analysis was performed to obtain data on the types of ecological sub-health. An ecological sub-health region layer is generated from the ecological sub-health type data to obtain ecological sub-health region data.

[0021] In this invention, the identification of ecologically sub-healthy areas utilizes a dual-domain energy gradient map and local anomaly detection to achieve precise spatial localization of vertical energy disturbance anomalies in forest areas. The local anomaly detection step identifies high-risk points such as abrupt changes in energy gradients, conduction discontinuities, and direction reversals, providing structural support for boundary extraction of ecologically damaged areas. Gradient continuity disruption identification helps identify areas where energy conduction chains are interrupted due to ecological degradation, disturbance, or fragmented terrain. Combined with cluster analysis, sub-healthy areas are categorized, enabling the system not only to identify risk areas but also to classify and analyze their formation mechanisms and characteristics.

[0022] The beneficial effects of this invention are as follows: Step S1 acquires canopy radiation temperature data through infrared thermal imaging, and estimates the sensible and latent heat fluxes of the canopy by combining meteorological parameters, forming an energy metabolism flux index for the canopy, which directly reflects the vegetation's evapotranspiration capacity and physiological activity. Step S2 constructs a high-precision terrain model based on laser point cloud data, and further integrates slope aspect dispersion and structural fragmentation index to extract surface roughness features, which are used to represent the degree of surface disturbance and energy albedo. Step S3 achieves vertical coupling modeling of energy in both the canopy and surface domains through spatial registration, calculates the energy difference and gradient distribution between the two, and accurately reflects the energy coordination and conduction continuity within the ecosystem. Step S4, based on gradient anomaly detection and damaged area identification, combined with spatial clustering methods, achieves the classification and spatial positioning of ecologically sub-healthy areas, providing clear guidance for ecological risk management and resource intervention. Attached Figure Description

[0023] Other features, objects, and advantages of this application will become more apparent from the following detailed description of the non-limiting embodiments, taken with reference to the accompanying drawings: Figure 1 A flowchart illustrating the steps of a forestry ecosystem health assessment method based on multi-source indicators is shown in one embodiment. Figure 2 A flowchart illustrating the steps of a method for calculating canopy energy flux according to an embodiment is shown. Figure 3 A flowchart illustrating the steps of a method for constructing a surface roughness model according to an embodiment is shown. Figure 4 A flowchart illustrating the steps of a dual-domain energy gradient processing method according to one embodiment is shown. Figure 5 A flowchart illustrating the steps of an embodiment of a method for identifying ecologically sub-healthy areas is shown. Detailed Implementation

[0024] The technical method of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without inventive effort are within the scope of protection of this invention.

[0025] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0026] It should be understood that although the terms "first," "second," etc., may be used herein to describe various units, these units should not be limited by these terms. These terms are used merely to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, a first unit may be referred to as a second unit, and similarly, a second unit may be referred to as a first unit. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0027] Please see Figures 1 to 5 This application provides a method for assessing the health of forest ecosystems based on multi-source indicators, including the following steps: Step S1: Acquire infrared data and point cloud data of the forest area; calculate the canopy energy flux based on the infrared data of the forest area to obtain the canopy energy flux data; Specifically, airborne / vehicle-mounted thermal infrared sensors are used to image the target forest area, acquiring apparent brightness and temperature data. A pre-set atmospheric correction model is then used to perform radiometric correction on the original images, adjusting parameters such as atmospheric transmittance, surface emissivity, and albedo to obtain corrected results reflecting the actual surface temperature of the tree canopy. Combined with collected meteorological station monitoring data (including wind speed, air temperature, air humidity, and solar radiation flux), the system sequentially calculates the net radiation, sensible heat flux, and latent heat flux of the canopy layer. For example, when calculating the net radiation of the canopy... ,in Net radiation of the canopy, Canopy albedo The incident solar radiation flux represents the intensity of shortwave solar radiation reaching the top of the canopy. Canopy emissivity represents the ability of the canopy surface to emit long-wave radiation. Stefan-Boltzmann constant represents the energy constant radiated per unit area of ​​a blackbody per unit time. The surface temperature of the canopy. Given the air temperature; calculate the sensible heat flux. ,in For sensible heat flux, air density, For isobaric specific heat capacity, The surface temperature of the canopy. For air temperature, Aerodynamic drag represents the resistance of airflow to heat transfer; latent heat flux is calculated. ,in For latent heat flux, Net radiation of the canopy, Sensible heat flux is defined as follows: net radiation represents the difference between the energy absorbed and emitted by the canopy; sensible heat flux represents the intensity of heat exchange caused by the temperature difference between the canopy and the air; and latent heat flux represents the energy consumed by the canopy per unit area due to transpiration. The system aggregates the sensible and latent heat fluxes of the canopy according to spatial location, generating a canopy energy flux data raster layer, expressed in the form of power density per unit area (watts per square meter).

[0028] Step S2: Construct a surface roughness model based on forest area point cloud data to obtain a surface roughness model; perform surface energy albedo processing based on the surface roughness model to obtain surface energy albedo data; Specifically, the lidar point cloud is filtered and classified to extract ground point sets and generate a digital elevation model (DEM) with a resolution of 1–5 meters. A 3×3 sliding window is used to calculate surface structure features on the DEM, including local undulation (characterizing the strength of elevation undulation within the window), aspect distortion (characterizing the dispersion and rotation of aspect within the window), and patch structure fragmentation (characterizing the complexity and connectivity of surface boundaries). These features are normalized and summarized according to preset weights or directly summarized to form a structure-enhanced surface roughness index / vector.

[0029] In energy albedo processing, multispectral images strictly registered with the digital elevation model are acquired. Surface reflectance is retrieved using red, green, and near-infrared bands and coupled with the aforementioned roughness index through same-pixel registration to construct a joint input layer. An albedo correction model is established based on this joint layer: while maintaining the physical meaning of the original reflectance, a preset roughness influence factor is set so that the corrected albedo increases appropriately with increasing surface roughness, thus more closely reflecting the actual scattering and multiple reflection effects. Spatial perturbation detection (e.g., calculating local coefficients of variation) is performed on the corrected albedo layer to identify and balance small-scale anomalous fluctuations, outputting a surface energy albedo layer for energy balance calculations, expressed using power density per unit area (W / m²).

[0030] Step S3: Perform dual-domain energy gradient processing based on canopy energy flux data and surface energy albedo data to obtain dual-domain energy gradient data; Specifically, spatial registration is performed on the canopy energy flux raster and the surface energy albedo raster to ensure consistency in resolution, coordinate system, and projection method. A vertical energy channel is established for each pixel, and the energy flux difference between the canopy layer and the surface layer at the same spatial location is calculated. Spatial gradient calculations are performed on the energy difference layer using edge detection operators (such as Sobel or Prewitt operators) to extract the amplitude and direction information of energy changes. The resulting dual-domain energy gradient layer can visually reveal the uneven energy distribution between the canopy and the surface.

[0031] Step S4: Identify ecologically sub-healthy regions based on dual-domain energy gradient data to obtain ecologically sub-healthy region data.

[0032] Specifically, local anomaly detection is performed based on a dual-domain energy gradient layer. Standardized analysis of energy gradient values ​​at different spatial locations is conducted to identify high-gradient regions with prominent energy changes. Spatial clustering analysis is performed on these anomalies to identify and construct interconnected sets of anomaly regions. The system identifies gradient continuity breaks, detecting abrupt changes in energy gradients and selecting regions with clear energy abrupt boundary boundaries as candidate areas for ecological damage or degradation. The system uses clustering algorithms (such as density clustering or Gaussian mixture models) to classify these anomaly regions, categorizing them into different sub-health types based on their energy characteristics. These include metabolic decline (manifested as weakened energy exchange capacity), energy breakage (manifested as interrupted vertical conduction), and structural fragmentation (manifested as impaired spatial continuity). For example, by statistically analyzing characteristic parameters such as the average energy gradient value, gradient direction consistency index, spatial connectivity, and energy fluctuation frequency of each region / cluster, the clusters are matched with known ecological degradation patterns. For example, when the overall energy gradient of a region is low and the fluctuation range is gentle, it is classified as "metabolic decline type"; when there is a significant abrupt change in the vertical gradient and the energy conduction between upper and lower layers is interrupted, it is classified as "energy fracture type"; when the spatial structure is fragmented and the energy direction is disordered, it is classified as "structural fragmentation type". An ecological sub-health classification layer is output to form a distribution map of ecological sub-health areas in forest regions.

[0033] Preferably, the calculation of canopy capacity flux is as follows: Step S11: Perform radiation brightness temperature correction based on forest area infrared data to obtain canopy surface temperature data; Specifically, thermal infrared remote sensing image data of the forest area is acquired, and atmospheric correction and surface emissivity correction are performed on the original brightness temperature image. The system can use a radiative transfer model (such as MODTRAN) and combine it with parameters such as measured / historical temperature, humidity and atmospheric transmittance from regional meteorological stations to correct the deviations caused by atmospheric absorption and scattering, thereby retrieving the canopy surface temperature layer.

[0034] Step S12: Calculate the net radiation based on the canopy surface temperature data to obtain the canopy net radiation data.

[0035] Specifically, the system combines parameters such as canopy surface temperature, total incident solar radiation flux, and canopy albedo to calculate the difference between the radiant energy absorbed and released per unit area of ​​the canopy, i.e., the net radiation. ,in Net radiation of the canopy, Canopy albedo The incident solar radiation flux represents the intensity of shortwave solar radiation reaching the top of the canopy. Canopy emissivity represents the ability of the canopy surface to emit long-wave radiation. Stefan-Boltzmann constant represents the energy constant radiated per unit area of ​​a blackbody per unit time. The surface temperature of the canopy. The air temperature is used as the reference point. The calculation process consists of two parts: first, determining the actual absorbed energy based on the canopy's absorption of shortwave solar radiation; and second, estimating the released longwave radiation energy based on the temperature difference between the canopy and the air, and its emission characteristics. The difference between the two is the net radiation value of the canopy.

[0036] Step S13: Obtain meteorological monitoring data, and estimate sensible heat flux based on canopy net radiation data and meteorological monitoring data to obtain canopy sensible heat flux data; Specifically, the system acquires meteorological monitoring data within the study area, including parameters such as wind speed, air temperature, air density, and specific heat capacity. Combined with canopy net radiation data, and based on the stability-corrected Monin-Obukhov similarity theory, the heat exchange between the canopy and the atmosphere is calculated. By considering factors such as airflow turbulence characteristics and aerodynamic drag, the energy transfer efficiency between the canopy and the air is reflected. The system estimates the sensible heat flux transferred per unit area of ​​the canopy by incorporating variables such as wind speed, the difference between canopy surface temperature and air temperature, air specific heat capacity, and aerodynamic drag. ,in For sensible heat flux, air density, For isobaric specific heat capacity, The surface temperature of the canopy. For air temperature, Aerodynamic drag represents the resistance of airflow to heat transfer. This flux indicates the intensity of heat transferred from the canopy layer to the air via thermal convection and is an important indicator for assessing the energy balance and microclimate regulation capacity of forest areas.

[0037] Step S14: Estimate the latent heat flux based on the net radiation data and sensible heat flux data of the canopy to obtain the evapotranspiration heat dissipation flux data of the canopy; Specifically, the system calculates the energy balance at the canopy level. Based on the obtained data on net canopy radiation and sensible heat flux, the portion of net radiation not used for sensible heat transfer is considered as energy consumed by vegetation evapotranspiration, i.e., latent heat flux. ,in For latent heat flux, Net radiation of the canopy, This refers to sensible heat flux. This flux represents the potential heat energy released by the canopy per unit area to the outside world through water evaporation and stomatal transpiration, reflecting the vegetation's comprehensive ability to regulate water and energy. A higher latent heat flux value indicates stronger evapotranspiration and more active physiological metabolism in the forest area; conversely, a lower value may indicate that the ecosystem is under water stress or in a state of functional decline.

[0038] Step S15: Integrate the canopy sensible heat flux data and the canopy evapotranspiration heat dissipation flux data to obtain the canopy energy flux data.

[0039] Specifically, the system registers and aligns the calculated sensible heat flux and evapotranspiration flux of the canopy using a unified spatial grid. The sensible and latent heat fluxes at the same spatial location are then integrated to obtain a canopy energy flux dataset reflecting the total energy exchange across the canopy levels. This data is measured in power density per unit area (watts per square meter).

[0040] Preferably, the surface roughness model is constructed as follows: Step S21: Construct an elevation raster model based on forest area point cloud data to obtain the elevation raster model; Specifically, the acquired raw lidar point cloud data is preprocessed, including the removal of noise points, outliers, and non-ground reflection points. A progressive triangulation filtering algorithm is used to extract ground points from the point cloud, separating a set of points representing the landform. The system performs spatial interpolation calculations on the ground point set, using methods such as inverse distance weighted interpolation or kriging interpolation, to generate a digital elevation model with a resolution of 0.5 to 2 meters. In the output elevation raster model, each pixel corresponds to a ground elevation value at a center point, labeled in meters.

[0041] Step S22: Extract local features from the elevation raster model to obtain local feature data; Specifically, the system uses a fixed-size sliding window (e.g., 5×5 pixels) as the analysis unit to extract various local indicators reflecting micro-topographic changes from the digital elevation model. The standard deviation and range of elevation values ​​within the window are calculated, i.e., the undulation index. Based on the aspect information derived from the elevation model, the system statistically analyzes the dispersion of the aspect angle within the window, using a circular standard deviation to measure the strength of aspect change, thus obtaining the aspect distortion index. For areas with high undulation or slope, spatial connectivity analysis is performed, calculating the ratio of the boundary length to the area of ​​each topographic patch, thus obtaining the patch structure fragmentation index. These three local indicators are normalized and fused according to preset weights, or directly spliced ​​after normalization to form the local feature vector corresponding to each pixel.

[0042] Step S23: Calculate the roughness of the local feature data to obtain the roughness data; Specifically, the system calculates the surface roughness index based on the aforementioned extracted local feature indicators using a pre-defined structure-enhanced surface roughness index calculation model. This model uses three indicators—undulation, aspect distortion, and patch structure fragmentation—as the main input parameters and performs a linear weighted summation calculation according to empirically set weights. The weights can be determined empirically based on the topographic features of different forest areas, or they can be adaptively inverted through training with sample data. The system performs the calculation pixel-by-pixel across the entire region, generating a roughness index layer and normalizing the results to a range of 0 to 1. A higher roughness index value indicates more severe surface undulation, more complex structure, and poorer spatial continuity in the region.

[0043] Step S24: Calculate the surface scattering characteristic parameters based on the roughness data to obtain the surface roughness model.

[0044] Specifically, the system uses the roughness index calculated above as an input factor to construct a model for adjusting the effect of micro-topography on surface energy reflection characteristics. This model corrects the original spectral reflectance through roughness influence parameters, ensuring that the reflectance reflects the multiple scattering and reflection enhancement effects brought about by surface microstructures. ,in This is the corrected surface albedo. The original spectral reflectance represents the surface reflectance value calculated from multispectral or hyperspectral remote sensing data; this is a preset value. The roughness adjustment factor represents an empirical coefficient used to control the degree of energy reflection enhancement caused by surface roughness. It can be set based on measured samples from typical areas or the results of radiation balance experiments. The surface roughness index represents the surface roughness index value at the corresponding spatial pixel location. The corrected reflectance value reflects the true energy scattering characteristics of the surface, and its variation is positively correlated with the roughness index; that is, the rougher the surface, the stronger the energy scattering. The roughness adjustment factor can be set based on measured data or empirical parameters of typical sample areas to control the adjustment range of the albedo. The output surface roughness model is a composite layer, containing the spatial location of each pixel, the roughness index, and the corresponding scattering enhancement parameters.

[0045] Preferably, local feature extraction specifically involves: Local relief is calculated for the elevation grid model to obtain local relief data; Specifically, the system uses a fixed-size sliding window (e.g., 5) on the digital elevation model (DEM). 5 or 7 Using 7 pixels as the analysis unit, elevation data within a window range is extracted. The degree of fluctuation in surface elevation within this neighborhood is calculated, including the standard deviation of elevation (representing the dispersion of elevation value distribution) and the range of elevation (representing the magnitude of the elevation difference between the highest and lowest points). Based on the above, the system generates a local undulation layer, where each pixel records the intensity of elevation change in its neighborhood, in meters (m).

[0046] The slope direction variation data is obtained by calculating the slope direction variation based on the local undulation data; Specifically, the system calculates the aspect layer based on the Digital Elevation Model (DEM), with each cell corresponding to a slope angle value ranging from 0 to 360 degrees. The system uses a fixed-size sliding window (e.g., 5... 5 or 7 The system uses 7 pixels as the analysis unit to statistically analyze the changes in all aspect angles within the window. It employs circular statistics to perform dispersion analysis on aspect angles, representing the degree of concentration and directional consistency of the aspect distribution. ,in The standard deviation is circular, representing the intensity of the change in aspect angle at the current position (i,j) within the sliding window. It is the natural logarithm function. For the vector composition magnitude, The number of pixels in the neighborhood. The cosine values ​​of the slope angles of all cells within the window. The slope angle represents the slope angle value of the pixel located at coordinates (i,j). The data represents the sine of the aspect angles of all pixels within the window. The upper limit is the total number of pixels in the sliding window, and the lower limit is the center pixel within the sliding window. If the aspect angles within the window change little, the dispersion is low, indicating stable slope direction and continuous terrain morphology. If the aspect angles are scattered and the direction changes drastically, the dispersion is high, indicating large slope distortion and complex terrain directionality. The system generates a slope direction variation layer, and each pixel records the intensity of its local aspect change.

[0047] Based on the slope cutting direction variation data, a patch structure connectivity map is constructed to obtain patch structure connectivity map data; Specifically, the system normalizes the slope direction variation layer, limiting its values ​​to between 0 and 1. Based on a preset threshold (e.g., values ​​greater than 0.7 are considered high-direction disturbance areas), the region is divided into high-disturbance and low-disturbance areas, generating a binarized disturbance mask map. The system uses eight-neighbor connected component labeling to perform spatial connectivity analysis on the mask map, aggregating mutually contacting or adjacent high-disturbance pixels into independent surface patches. Each patch is assigned a unique number in the connectivity layer, and its geometric attributes such as area, boundary length, principal axis direction, and shape index (the system calculates the shape index value of each patch by comparing the actual boundary length with the circumference of an equivalent circle) are calculated and recorded.

[0048] The patch structure fragmentation index is calculated based on the patch structure connectivity graph data to obtain the patch structure fragmentation index data. Specifically, based on the previously generated patch connectivity layer, the system calculates the structural features of all identified connected surface patches. It statistically analyzes the boundary length and area of ​​each patch, and calculates the ratio of boundary length to area to obtain the boundary complexity; the more tortuous and irregular the boundary, the higher this ratio. The system calculates the patch number density, i.e., the distribution of the number of patches per unit area. The system calculates the average patch shape index. Combining multiple indicators such as patch number, boundary complexity, and area, the system constructs a patch structure fragmentation index / fragmentation index value. ,in This is the breakage index value. For the number of plaques, The average boundary length per unit area. The total area of ​​the region. This is an area stabilization correction term, introduced to prevent numerical anomalies caused by extremely small areas. A higher fragmentation index value indicates a denser distribution of patches, more complex boundaries, and a more discontinuous spatial structure within the region. The system assigns the fragmentation index value corresponding to the patch to each pixel, generating a patch structure fragmentation index layer.

[0049] Local feature data are obtained by integrating local undulation data, slope direction variation data, patch structure connectivity map data, and patch structure fragmentation index data.

[0050] Specifically, the system performs unified registration processing on various spatial layers to ensure consistency in spatial resolution, coordinate system, and projection method across different data sources. Using the pixel as the smallest unit of analysis, the system performs raster alignment on the undulation layer, aspect disturbance layer, patch connectivity layer, and patch fragmentation index layer, ensuring that each layer corresponds to the same pixel coordinates at the same spatial location. The system constructs a local feature vector for each pixel containing multi-dimensional terrain features, including the intensity of elevation undulation, the intensity of local aspect change, the patch number, and the corresponding fragmentation index value.

[0051] Preferably, the roughness calculation is as follows: Local elevation relief layers are constructed based on local feature data and elevation raster models to obtain elevation relief map data. Specifically, based on the existing elevation raster model (DEM) and combined with previously extracted local feature data, the system performs scale optimization and stability enhancement processing on local undulation to improve its response to micro-undulation disturbances and suppress misjudgments of noise-dominated areas. During processing, the system sets a dynamic sliding window range based on the local slope aspect change intensity and fragmentation index values ​​in the local feature data; that is, a smaller window is used for areas with higher values ​​to enhance local response, and a larger window is used to smooth out minor fluctuations. Elevation data is extracted within each dynamic window, and the standard deviation and range of elevation are calculated. Simultaneously, an undulation volatility index / surface undulation intensity value (such as the ratio of elevation difference to the elevation of the central pixel) is introduced. The system normalizes and performs local smoothing filtering on the preliminary calculation results to construct a stability-enhanced local elevation undulation layer. Each pixel records its surface undulation intensity value after scale optimization and anomaly correction, resulting in elevation undulation map data.

[0052] The aspect distortion layer is extracted from the elevation raster model to obtain the aspect distortion layer data; Specifically, a multi-scale aspect difference layer is constructed on the aspect layer, including a radial difference layer: calculating the aspect difference between the central pixel and its four directional (north, south, east, west) pixels; a diagonal difference layer: calculating the degree of abrupt change in aspect angle with its diagonal neighbors; and a maximum direction jump layer: recording the direction and magnitude of the largest aspect difference between the current pixel and its neighbors. The system sets an angle stability index, constructing a composite index reflecting the trend of local directional abrupt changes by statistically analyzing the mean, maximum, and variance of the angle change magnitude in the above difference maps. This index shows a low value in areas with good aspect continuity and a high value in areas with multiple overlapping directions, concentrated inflection points, or significant breakage trends. The system also performs a removal operation on pixels with slopes below a set threshold (e.g., 5°) to prevent misidentification of aspect disturbances in micro-slope areas. The system assigns this index to the center of each pixel, generating aspect distortion layer data.

[0053] Based on local feature data, elevation relief map data, and slope aspect distortion layer data, patch structure fragmentation layer is extracted to obtain patch structure data. Specifically, the system normalizes the elevation relief layer and the aspect distortion layer to make data from different sources comparable on the same numerical scale. Following preset weights (e.g., 0.5 for each, or dynamically adjusted weights based on sample data training), the system performs linear weighted fusion on the two types of feature data to generate a perturbation potential map. The system then performs threshold segmentation on the perturbation potential map, identifying areas exceeding a set threshold as high-perturbation zones. For these zones, an eight-neighbor connected component analysis method is used to aggregate spatially continuous high-perturbation pixels into a complete connected patch. For each connected patch, the system calculates its structural characteristic indices, including the ratio of boundary length to area (representing boundary complexity), shape factor (representing the patch's roundness and compactness), and the average distance between neighboring patches (representing fragmentation density and spatial uniformity). The system then performs weighted fusion / combination / stitching of these indices back to the elevation model to form a patch fragmentation index layer, obtaining patch structure data.

[0054] Roughness data is obtained by performing roughness weighting on the patch structure data.

[0055] Specifically, the system performs weighted calculations on surface structural features based on three input layers: elevation undulation layer, aspect distortion layer, and patch fragmentation layer. The system performs numerical normalization and spatial registration on each input layer to ensure consistent data scale and spatial correspondence. Using pre-set or sample-trained weighting coefficients (e.g., undulation weight 0.4, aspect distortion weight 0.3, patch fragmentation weight 0.3), the three types of features are weighted and fused to construct a structure-enhanced roughness index model. This model combines surface elevation undulation, slope orientation disturbance, and structural fragmentation features to represent the irregularity and disturbance intensity of regional micro-topography from multiple perspectives. The system calculates the corresponding roughness value for each pixel.

[0056] Preferably, the extraction of the slope aspect distortion layer specifically involves: Slope aspect layer data is obtained by extracting the slope aspect layer from the elevation raster model. Specifically, the system generates an aspect layer based on a digital elevation model (DEM) using an aspect calculation algorithm (such as the Horn algorithm or the Zevenbergen & Thorne algorithm). The aspect calculation analyzes the elevation difference relationship between each pixel and its eight neighboring pixels to determine the maximum descent direction of the slope at that location. The calculation result is output as an aspect angle value, ranging from 0 to 360 degrees, representing the azimuth of the slope orientation. The generated aspect layer maintains consistency with the original elevation model in spatial resolution and coordinate system, with each pixel recording its corresponding aspect angle in degrees.

[0057] Calculate the neighboring slope aspect value from the slope aspect layer data to obtain the neighboring slope aspect value data; Specifically, the system sets a fixed-size sliding window (e.g., 5) on the slope aspect layer. The system uses 5 pixels as the center and performs neighborhood analysis sequentially around each pixel. For each center pixel, the system extracts the aspect angle values ​​of all pixels within the window, forming a set of neighborhood aspect values ​​for that location. To avoid numerical jumps in aspect angles at the 360-degree cycle boundary (e.g., 5 degrees and 355 degrees should be considered close in direction rather than differing by 350 degrees), the system converts all aspect angles into unit vector form to obtain aspect direction information. The neighborhood aspect value data output by the system contains the aspect direction information within the neighborhood of each pixel.

[0058] Local aspect dispersion data is obtained by calculating the local aspect dispersion based on the neighboring aspect value data. Specifically, the system extracts the corresponding set of neighborhood aspect vectors for each pixel to statistically analyze the consistency of aspect distribution within that area. The system averages all aspect vectors within the neighborhood to obtain an average vector representing the overall dominant slope direction. The system calculates the magnitude of this average vector to measure the concentration of aspect directions within the neighborhood: a larger magnitude indicates concentrated aspect and stable terrain direction; a smaller magnitude indicates dispersed aspect distribution and complex directional changes. The system calculates a circular standard deviation based on the magnitude of the average vector. A higher circular standard deviation indicates more drastic aspect changes and a more complex surface structure. The system writes the circular standard deviation value corresponding to each central pixel into an aspect dispersion layer to generate local aspect dispersion data.

[0059] Based on the local dispersion data of slope aspect, the flat areas of the slope aspect layer data are removed to obtain the slope aspect distortion layer data.

[0060] Specifically, the system calculates the slope value of each pixel based on a digital elevation model (DEM) to represent the degree of surface tilt. A lower slope threshold is set (e.g., 5 degrees, or dynamically determined based on the statistical percentile of the regional topographic distribution) to identify flat areas. When the slope value of a pixel is below this threshold, or when the slope standard deviation in its neighborhood is small and the overall variation is not significant, the system determines that the area is a flat area. For such pixels, the system sets the corresponding aspect dispersion value (CSD value) to zero or marks it as invalid in the aspect dispersion layer. The remaining pixels retain the true aspect dispersion results, thus obtaining the aspect distortion layer data. Preferably, the surface energy albedo treatment specifically includes: Obtain surface spectral reflectance data; Specifically, the system acquires multispectral remote sensing data consistent with the spatial extent of the study area. Data sources can include satellite imagery, Landsat 8 / 9 imagery, or multispectral imaging systems mounted on UAVs. For the acquired raw imagery, the system processes it using atmospheric correction methods (such as Sen2Cor or FLAASH modules) to convert the raw radiance values ​​into surface reflectance. In the corrected surface reflectance imagery, the system selects key bands for surface energy analysis and vegetation structure modeling, including the red band, near-infrared band (NIR), and shortwave infrared band (SWIR). The system outputs a spectral reflectance layer expressed as dimensionless reflectance values ​​(range 0 to 1), whose spatial resolution must be consistent with the roughness layer, or unified to a specified standard through resampling.

[0061] The surface roughness model and surface spectral reflectance data are registered to obtain joint surface data. Specifically, the system acquires the constructed surface roughness model layer, where each pixel contains a structure-enhanced roughness index. The system performs spatial registration between the surface spectral reflectance layer and the roughness layer. After registration, the system combines the roughness value and spectral reflectance value at each location, pixel by pixel, to construct a joint input data structure. This joint data includes both surface geometric features and surface physical spectral response information. The system can also set other auxiliary variables, such as Normalized Difference Vegetation Index (NDVI), slope, and surface moisture, according to actual application needs. The generated surface joint data layer is a multi-feature composite raster structure.

[0062] Roughness reflectance coupling processing is performed on the surface joint data to obtain surface coupled data; Specifically, based on the aforementioned constructed surface joint data layer, the system uses a roughness adjustment model to correct the original surface reflectance. The system supports various coupling model forms, including linear coupling models (such as...). ,in The effective albedo after coupling is divided by the corrected original spectral reflectance. The original surface spectral reflectance, This is an empirical adjustment coefficient. (Surface roughness index / roughness value) and exponential coupling model (e.g.) ,in The effective albedo after coupling is divided by the corrected original spectral reflectance. The original surface spectral reflectance, is the base of the natural logarithm. This is an empirical adjustment coefficient. The roughness index (or roughness value) can be flexibly selected based on the actual application scenario. The roughness index serves as a structural perturbation factor, while an empirical adjustment coefficient controls the sensitivity of roughness to reflectance correction. This parameter is set between 0.2 and 0.5 and can be obtained by fitting measured data from typical sample areas. During the coupling process, the system combines the roughness value with the original spectral reflectance value pixel-by-pixel, coupling the surface data. The corrected albedo layer has the same units and spatial resolution as the original layer, with values ​​ranging from 0 to 1.

[0063] The surface coupling data is processed by albedo perturbation mapping to obtain surface perturbation energy data; Specifically, based on the coupled surface albedo layer, the system employs a sliding window analysis method at the pixel level to extract variation characteristics within local areas. Within each window, the system calculates the standard deviation of the local albedo; it applies first-order gradient operations (such as the Sobel operator) to identify the spatial direction of albedo value changes and abrupt change edges; and it calculates second-order variation indices (such as the Laplacian operator) to represent the severity and fluctuation trend of albedo changes within the region. Through the combination of these multiple indices, the system constructs an albedo perturbation layer to represent the spatial heterogeneity of surface energy reflection. Higher values ​​in the layer indicate more drastic reflectivity changes and a more complex surface energy reflection structure in the region. The system can also set quantile thresholds based on the distribution of perturbation values ​​to classify perturbation intensity into weak, moderate, and strong levels. The surface perturbation energy data output by the system serves as an albedo spatial response map.

[0064] Surface energy albedo is obtained by performing surface energy perturbation data.

[0065] Specifically, the system acquires solar incident irradiance data matching the study period and region. Data sources may include ground-based meteorological station observations, ERA5 meteorological reanalysis data, or satellite radiation products such as MODIS. Irradiance data can be total irradiance across all wavelengths, or shortwave radiant flux consistent with the remote sensing reflectance band, measured in watts per square meter. Based on the acquired incident radiant flux, the system combines the coupled surface albedo layer to calculate the surface reflected energy for each pixel. , Energy reflected from the Earth's surface For solar incident irradiance data, This is the surface albedo layer, which represents the corrected original surface reflectance. It estimates the reflected energy under current solar radiation conditions based on the albedo level of a local area. The system uses previously generated perturbation factor or perturbation index layers as weighting terms and sets an albedo perturbation adjustment mechanism to better reflect the impact of surface micro-topographic differences and structural perturbations on energy reflection. The output surface energy albedo layer is expressed in watts per square meter, representing the actual energy reflection distribution characteristics of different regions under micro-scale perturbations.

[0066] Preferably, step S3 specifically includes: Step S31: Spatial projection alignment is performed based on canopy energy flux data and surface energy albedo data to obtain vertical energy channel data; Specifically, the system acquires the required input layers, including a canopy energy flux layer (watts per square meter) and a surface energy albedo layer, used to describe the energy output characteristics of the canopy and surface layers, respectively. The system performs a consistency check on the spatial features of the two layers to ensure alignment in terms of spatial resolution (e.g., 10 meters or 30 meters), projection coordinate system (e.g., WGS84 or UTM projection), map boundaries, and layer size. When there are resolution or projection differences between the two layers, the system prioritizes the layer with the higher resolution or the dominant coordinate system as the reference, using bilinear interpolation or nearest neighbor interpolation to resample the other layer, ensuring a one-to-one correspondence between the two layers at pixel locations. After resampling and spatial registration, the system merges the two layers according to the pixel correspondence, combining the canopy energy flux value and the surface energy albedo value at each spatial location to construct a vertical energy channel unit. In the output vertical energy channel layer, each pixel contains two energy components, representing the energy expression of the canopy and the surface, respectively.

[0067] Step S32: Calculate the dual-domain energy difference for the vertical energy channel data to obtain dual-domain energy difference data; Specifically, based on the registered vertical energy channel layer, the system calculates the energy difference between the canopy layer and the surface layer pixel by pixel. The difference calculation can be performed in two ways: one is an absolute difference model, which directly calculates the numerical difference between the canopy energy flux and the surface energy albedo; the other is a proportional difference model, which normalizes the difference by the ratio of the energy difference to its total amount. To avoid division by zero, the system sets a small constant in the denominator as a compensation factor. Depending on the actual needs, the system can normalize the calculation results, limiting the difference to the range of 0 to 1. The generated dual-domain energy difference layer can be stored in the original energy units (e.g., watts per square meter) or processed into a dimensionless form, depending on the selected model.

[0068] Step S33: Construct the gradient field based on the dual-domain energy difference data to obtain dual-domain energy gradient data.

[0069] Specifically, the system performs first-order spatial derivative operations on the dual-domain energy difference layer obtained from the aforementioned calculations. The system employs edge enhancement operators commonly used in image processing for gradient calculation, such as the Sobel operator, or, depending on accuracy requirements, the Prewitt, Roberts, or Scharr operators. Through derivative calculations in the horizontal and vertical directions, the system generates two output layers: a gradient magnitude layer, representing the intensity of energy difference changes in space; and an optional gradient direction layer, indicating the main direction of energy change in degrees. The output gradient magnitude layer serves as the dual-domain energy gradient map.

[0070] Preferably, step S4 specifically includes: Step S41: Construct a dual-domain energy gradient map based on the dual-domain energy gradient data to obtain dual-domain energy gradient map data; Specifically, the system analyzes the degree of gradient direction shift within the local pixel cluster corresponding to the dual-domain energy gradient data. If the direction change exceeds a preset angle threshold, the system performs a compression and smoothing operation on the gradient magnitude in that region based on a preset direction consistency penalty factor. This effectively suppresses misleading boundary responses caused by high-frequency noise and improves the coherence and physical interpretability of the gradient map in ecological fracture identification. The system performs structure-sensitive multi-scale filtering operations, and presets multiple gradient convolution kernels of different scales (e.g., 3D). 3, 5 5, 7 The system performs weighted response calculations on the dual-domain energy difference layers (7 levels). For each pixel, the system evaluates its structural boundary response intensity at each scale, prioritizing the scale result with the most obvious performance as the output. Specifically, it calculates the mean or variance of the local gradient magnitude for each pixel as the structural boundary response intensity evaluation index; the scale corresponding to the largest index value is selected as the output value for that pixel. After gradient convolution, the system overlays an external digital elevation model or slope data. Background adaptive adjustment is implemented in areas with significant terrain undulations. In areas with abrupt slope changes, the system sets a suppression factor based on the possibility of aspect projection, attenuating the gradient values ​​accordingly to avoid misidentifying energy differences caused by terrain shading or orientation changes as ecological anomalies. The system constructs a background energy gradient surface as a reference benchmark. This benchmark surface is calculated based on the overall map energy distribution trend, such as by constructing a continuous and smooth background model through moving average or low-pass filtering. The system repositions the energy gradient value of the current pixel relative to this background trend to obtain a relative evaluation result of the disturbance amplitude, resulting in dual-domain energy gradient map data.

[0071] Step S42: Perform local anomaly detection on the dual-domain energy gradient map data to obtain local anomaly data; Specifically, based on the constructed dual-domain energy gradient layer, the system performs statistical analysis on the gradient intensity of each pixel and sets a dynamic detection threshold to identify energy anomaly regions. This process can be achieved through a local Z-score method or a quantile threshold strategy: the former identifies anomalous pixels significantly above the average level by statistically analyzing the mean and standard deviation of all pixels in the image; the latter extracts high-value regions located in the upper layers of the gradient value distribution (e.g., the top 10%) as potential anomalies. The system performs spatial clustering tests on the initially identified anomalous pixels. Spatial autocorrelation indices (such as Moran's I or Ripley's K function) are used to analyze the spatial distribution characteristics of the anomalies. The system outputs local anomaly data, with each pixel assigned only a value of 0 or 1, representing normal areas and significantly anomalous areas, respectively.

[0072] Step S43: Identify the gradient continuity failure area of ​​the local abnormal data to obtain the failure area data; Specifically, the system identifies spatially adjacent anomalous pixels in the local anomaly layer, extracts connected patches using the 8-neighborhood connection rule, and constructs a preliminary set of anomalous regions. For each connected region, the system analyzes the spatial consistency of its internal gradient directions and calculates the gradient angle differences between pixels within the patch using the previously generated gradient direction map. If a significant abrupt change in gradient direction is found within a region (e.g., the angle change exceeds a set threshold), the region is determined to have a break in the continuity of energy flux transfer direction and is considered a potential area of ​​ecological structural damage. The system extracts key morphological attributes from these damaged regions, including indicators such as area, boundary length, and shape complexity, and sets a minimum ecological patch threshold (e.g., areas with an area less than 100 square meters are excluded) to remove unrepresentative noise patches from the analysis. The system outputs a mask layer of damaged regions, i.e., damaged region data, where the value of each pixel indicates whether it belongs to the identified gradient break zone. A value of 1 represents structural damage, and a value of 0 represents normal gradient continuity.

[0073] Step S44: Perform ecological sub-health cluster analysis based on the damaged area data to obtain ecological sub-health type data; Specifically, the system extracts various feature variables from the identified damaged areas, including basic indicators such as average energy difference, average gradient intensity, and regional morphological and structural characteristics (e.g., compactness, boundary tortuosity). It can also integrate ecological remote sensing auxiliary information such as the vegetation index (NDVI), bare land index (NDBSI), and soil moisture content to construct a multidimensional feature vector representing the ecological state of the patches. Based on unsupervised clustering algorithms (such as K-means, density clustering, or Gaussian mixture models), the system performs ecological feature clustering analysis on these patches, automatically identifying the differential distribution patterns of sub-healthy ecological areas in terms of energy characteristics, morphological structure, and vegetation response. The number of clusters can be adaptively selected based on model evaluation indicators (such as information criteria or silhouette coefficient) to ensure the stability and discriminative power of the classification results. The output is an ecological sub-healthy type layer, with each damaged area assigned a corresponding ecological type label, such as "weakened canopy metabolism type," "surface albedo imbalance type," or "structural fracture driven type."

[0074] Step S45: Generate an ecological sub-health region layer from the ecological sub-health type data to obtain ecological sub-health region data.

[0075] Specifically, the system maps the ecological sub-health type labels obtained from the aforementioned clustering analysis back to their original geospatial locations, constructing a complete raster layer containing the ecological type of each pixel. Image processing techniques such as morphological closing operations can be applied to smooth the regional boundaries. The system can overlay auxiliary layers such as administrative divisions, forest stand types, or ecological reserve boundaries, performing spatial cropping or zoning weighting to meet the accuracy requirements of specific regional governance or multi-level management. The output ecological sub-health region layer can serve as data support for ecological restoration priority zoning, forest degradation level classification, and ecological health evolution monitoring.

[0076] Therefore, the embodiments should be regarded as exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended application documents rather than the foregoing description. Thus, it is intended that all variations falling within the meaning and scope of the equivalents of the application documents be incorporated into the invention.

[0077] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. A method for assessing the health of forestry ecosystems based on multi-source indicators, characterized in that, Includes the following steps: Step S1: Acquire infrared data and point cloud data of the forest area; calculate the canopy energy flux based on the infrared data of the forest area to obtain the canopy energy flux data; Step S2: Construct a surface roughness model based on forest area point cloud data to obtain a surface roughness model; perform surface energy albedo processing based on the surface roughness model to obtain surface energy albedo data; Step S3: Perform dual-domain energy gradient processing based on canopy energy flux data and surface energy albedo data to obtain dual-domain energy gradient data; Step S4: Identify ecologically sub-healthy regions based on dual-domain energy gradient data to obtain ecologically sub-healthy region data.

2. The method according to claim 1, characterized in that, The calculation of canopy capacity flux is as follows: Radiation brightness temperature correction was performed based on infrared data from the forest area to obtain canopy surface temperature data; Net radiation data of the canopy is obtained by calculating the net radiation based on the canopy surface temperature data. Meteorological monitoring data was acquired, and sensible heat flux was estimated based on canopy net radiation data and meteorological monitoring data to obtain canopy sensible heat flux data; Latent heat flux was estimated based on canopy net radiation data and canopy sensible heat flux data to obtain canopy evapotranspiration heat dissipation flux data; By integrating the canopy sensible heat flux data and the canopy evapotranspiration heat flux data, the tree canopy energy flux data is obtained.

3. The method according to claim 1, characterized in that, The surface roughness model is constructed as follows: An elevation raster model was constructed based on forest area point cloud data to obtain the elevation raster model. Local feature extraction is performed on the elevation raster model to obtain local feature data; Roughness calculation is performed on local feature data to obtain roughness data; The surface scattering characteristic parameters are calculated based on the roughness data to obtain the surface roughness model.

4. The method according to claim 3, characterized in that, Local feature extraction specifically involves: Local relief is calculated for the elevation grid model to obtain local relief data; The slope direction variation data is obtained by calculating the slope direction variation based on the local undulation data; Based on the slope cutting direction variation data, a patch structure connectivity map is constructed to obtain patch structure connectivity map data; The patch structure fragmentation index is calculated based on the patch structure connectivity graph data to obtain the patch structure fragmentation index data. Local feature data are obtained by integrating local undulation data, slope direction variation data, patch structure connectivity map data, and patch structure fragmentation index data.

5. The method according to claim 3, characterized in that, Roughness calculation is as follows: Local elevation relief layers are constructed based on local feature data and elevation raster models to obtain elevation relief map data. The aspect distortion layer is extracted from the elevation raster model to obtain the aspect distortion layer data; Based on local feature data, elevation relief map data, and slope aspect distortion layer data, patch structure fragmentation layer is extracted to obtain patch structure data. Roughness data is obtained by performing roughness weighting on the patch structure data.

6. The method according to claim 5, characterized in that, The extraction of the aspect distortion layer is specifically as follows: Slope aspect layer data is obtained by extracting the slope aspect layer from the elevation raster model. Calculate the neighboring slope aspect value from the slope aspect layer data to obtain the neighboring slope aspect value data; Local aspect dispersion data is obtained by calculating the local aspect dispersion based on the neighboring aspect value data. Based on the local dispersion data of slope aspect, the flat areas of the slope aspect layer data are removed to obtain the slope aspect distortion layer data.

7. The method according to claim 1, characterized in that, The specific steps of surface energy albedo processing are as follows: Obtain surface spectral reflectance data; The surface roughness model and surface spectral reflectance data are registered to obtain joint surface data. Roughness reflectance coupling processing is performed on the surface joint data to obtain surface coupled data; The surface coupling data is processed by albedo perturbation mapping to obtain surface perturbation energy data; Surface energy albedo is obtained by performing surface energy perturbation data.

8. The method according to claim 1, characterized in that, Step S3 is as follows: Vertical energy channel data are obtained by spatial projection alignment based on canopy energy flux data and surface energy albedo data. The energy difference between the vertical energy channels is calculated in two domains to obtain the energy difference data between the two domains. Gradient fields are constructed based on the energy difference data in the two domains to obtain the energy gradient data in the two domains.

9. The method according to claim 1, characterized in that, Step S4 is as follows: A two-domain energy gradient map is constructed based on the two-domain energy gradient data to obtain the two-domain energy gradient map data. Local anomaly detection is performed on the dual-domain energy gradient map data to obtain local anomaly data; Gradient continuity disruption regions are identified in local anomaly data to obtain disrupted region data; Based on the data of the damaged areas, ecological sub-health cluster analysis was performed to obtain data on the types of ecological sub-health. An ecological sub-health region layer is generated from the ecological sub-health type data to obtain ecological sub-health region data.

Images