Forest resource dynamic supervision method based on multi-temporal laser radar

Abstract

The invention provides a forest resource dynamic supervision method based on a multi-temporal laser radar, and the method comprises the steps: carrying out the data collection of a target forest in different periods, and obtaining multi-temporal point cloud data, which is used for monitoring the change of the target forest in a time sequence; performing hierarchical processing on point cloud of each period in the multi-temporal point cloud data by taking an echo sequence and echo intensity as constraint conditions to obtain canopy layer point cloud data and ground layer point cloud data, calculating a crown change index and a ground change index, setting a crown change threshold value and a ground change threshold value, and calculating a crown change threshold value and a ground change threshold value according to the crown change index and the ground change index; the method comprises the following steps: identifying regions with significant change of crowns and grounds, realizing classification of crown and ground changes, obtaining crown and ground change analysis results, and drawing a forest resource change distribution map. According to the method, quantitative analysis is carried out by constructing a plurality of change indexes, automatic extraction, classification and mapping of forest crown and ground change information are completed in combination with change threshold segmentation, and forest resource changes are dynamically monitored.

Description

Technical Field

[0001] The present invention belongs to the technical field of dynamic supervision of forest resources, and in particular relates to a dynamic supervision method of forest resources based on multi-temporal laser radar. Background Art

[0002] Forest resource monitoring is crucial in addressing global challenges such as climate change, forest degradation and biodiversity loss. With the development of remote sensing technology, especially the application of technologies such as LiDAR, high-resolution satellite images and drones, modern forest resource monitoring can provide more accurate and timely monitoring data. At the same time, the combination of big data and artificial intelligence technology makes the processing and analysis of remote sensing data more efficient, and can achieve real-time change detection and trend prediction. However, the current dynamic forest resource supervision scheme has the following shortcomings:

[0003] (1) Most of them rely on data collection at a single time point and lack multi-temporal comparative analysis, making it difficult to monitor long-term changes in forest resources and unable to effectively capture long-term trends and short-term fluctuations caused by seasonal or climate change;

[0004] (2) More attention is paid to changes in the canopy layer, while less attention is paid to changes in the ground layer. Changes in the ground layer are usually assessed using relatively rough indicators, lacking detailed analysis of changes in the ground layer, and often ignoring the impact of ground undulations and human disturbance on forest ecosystems;

[0005] (3) Most of them use fixed thresholds for change detection, which is difficult to adapt to the changing needs of different monitoring areas and time periods. In addition, they usually monitor the canopy layer and the ground layer independently, lacking a data analysis and dynamic adjustment mechanism that comprehensively considers the changes in both layers, resulting in low accuracy in monitoring changes in forest resources.

[0006] Therefore, we need to develop a dynamic monitoring method for forest resources based on multi-temporal lidar, which can conduct comprehensive data collection, accurately analyze the point cloud data of the ground layer, and combine the changing characteristics of the canopy and ground layers to comprehensively reflect the dynamic changes of forest ecosystems. Summary of the invention

[0007] The purpose of the present invention is to provide a method for dynamic monitoring of forest resources based on multi-temporal lidar to solve the problems mentioned in the above background technology of the existing dynamic monitoring scheme of forest resources, such as single data time point, lack of ground layer data monitoring, and lack of comprehensive data analysis of canopy and ground layer change characteristics.

[0008] To achieve the above object, the present invention provides a method for dynamic monitoring of forest resources based on multi-temporal laser radar, and the method is specifically as follows:

[0009] Collecting data of the target forest at multiple different time periods, including but not limited to performing two or more laser radar scans at different time points, to obtain multi-temporal point cloud data of the target forest;

[0010] Taking the laser radar echo order and echo intensity as constraints, performing layered processing on the point cloud of each period in the multi-temporal point cloud data, dividing the point cloud into a canopy layer point cloud and a ground layer point cloud, and obtaining canopy layer point cloud data and ground layer point cloud data;

[0011] Based on the canopy point cloud data, a canopy change index of each period of the canopy point cloud is calculated, and based on the canopy change index, a canopy change threshold is set, and a canopy change analysis is performed on the canopy point cloud data of each period to identify areas where the canopy has changed significantly, and based on the combination pattern of the canopy change index, the nature and cause of the significant change in the canopy is judged to achieve classification of the canopy change and obtain a canopy change analysis result;

[0012] Based on the ground layer point cloud data, calculating the elevation distribution index of the ground layer point cloud in each period, so as to characterize the fluctuation change and human disturbance status of the target forest ground;

[0013] Based on the elevation distribution index, a ground change threshold is set, and a ground change analysis is performed on the ground layer point cloud data of each period. Areas with significant changes in ground elevation and ground point cloud density are identified and classified, and the spatial distribution pattern and possible causes of ground changes are obtained as the ground change analysis results.

[0014] Based on the above scheme, the crown change indicators include: crown point cloud density, average crown height, crown coverage and maximum crown height; after calculating the crown change index of the crown layer point cloud for each period, the change amount and change rate of the crown change index in different periods are also calculated to evaluate the structural dynamics of the crown within a certain time scale. The change amount reflects the increase or decrease in the absolute value of the crown change index, and the change rate reflects the rate of change of the crown change index relative to the initial value.

[0015] Based on the above scheme, based on the canopy point cloud data, a raster data layer is generated to intuitively display the spatial change distribution of the canopy change index, and the change amount raster and change rate raster of the canopy change index in different time periods are generated to obtain the spatiotemporal change layer of the canopy structure, which comprehensively displays the dynamic changes of the canopy layer;

[0016] After setting the crown change threshold, threshold segmentation of the crown layer point cloud data is performed to divide the significance levels of the change amount grid and the change rate grid of the crown change index. The crown change threshold includes but is not limited to a point cloud density change amount threshold, a point cloud density change rate threshold, a maximum height change amount threshold, a coverage change amount threshold, a coverage change rate threshold, an average height change amount threshold, and an average height change rate threshold.

[0017] Based on the above scheme, the point cloud density change threshold and the point cloud density change rate threshold are set as follows: for the crown point cloud density, two crown change thresholds are defined: point cloud density change threshold And the point cloud density change rate threshold , when the density change of the crown point cloud in area A of the target forest And the crown point cloud density change rate When , it is determined that a significant drop in the tree crown point cloud density has occurred in area A;

[0018] The setting of the coverage change threshold and the coverage change rate threshold is specifically as follows: for the tree crown coverage, the coverage change threshold is defined as and coverage change rate threshold , when the coverage change of area A And the coverage change rate When , there was a significant drop in tree crown coverage in area A.

[0019] Based on the above scheme, if area A in the target forest meets the threshold conditions of sudden drop in crown point cloud density and sudden drop in crown coverage at the same time, that is, and , and , and when these two sudden drops are concentrated and continuous in space, area A is judged to be a suspected deforestation area.

[0020] Based on the above scheme, the elevation distribution indicators include but are not limited to: average ground elevation, minimum ground elevation, and ground point cloud density;

[0021] The ground change threshold includes but is not limited to: a minimum elevation change threshold, an average elevation change threshold, a ground point cloud density change threshold, and a ground point cloud density change rate threshold;

[0022] The setting of the minimum elevation change threshold is specifically as follows: for the lowest ground elevation, the minimum elevation change threshold is set , which is used to determine whether there is a significant negative elevation anomaly in the lowest depression area on the ground. If the change in the lowest ground elevation in area A is , it is determined that significant negative ground elevation changes have occurred in area A.

[0023] Based on the above scheme, the setting of the average elevation change threshold is specifically as follows: for the average ground elevation, the average elevation change threshold is set , used to determine whether the ground as a whole has a significant uplift change. If the average elevation change in area A is > , it is determined that a significant positive ground elevation change has occurred in area A.

[0024] The setting of the ground point cloud density change amount threshold and the ground point cloud density change rate threshold is specifically as follows: for the ground point cloud density, the ground point cloud density change amount threshold is set. And the ground point cloud density change rate threshold , which is used to determine whether the sampling density of the ground point cloud has increased significantly. If the density change of the ground point cloud in area A is > , ground point cloud density change rate > , it is determined that the ground density in area A has increased significantly.

[0025] Based on the above scheme, if area A satisfies the following conditions at the same time: the change in the lowest ground elevation is less than the lowest elevation change threshold, the change in the ground point cloud density is greater than the ground point cloud density change threshold, and the ground point cloud density change rate is greater than the ground point cloud density change rate threshold, that is, a zero value area appears in the lowest ground elevation and the ground point cloud density increases suddenly, then area A is determined to be a newly added road area;

[0026] If area A simultaneously satisfies the following conditions: the average elevation change is greater than the average elevation change threshold, the ground point cloud density change is greater than the ground point cloud density change threshold, and the ground point cloud density change rate is greater than the ground point cloud density change rate threshold, but the lowest ground elevation is not less than the lowest elevation change threshold, then area A is determined to be a newly added stockpile area.

[0027] Based on the above scheme, combined with the crown change analysis results and the ground change analysis results, the resource change area data of the target forest is extracted and mapped, and the dynamic change process of the target forest resources is comprehensively described;

[0028] The resource change area includes but is not limited to the suspected tree felling area, the suspected tall tree felling area, the suspected tree natural death or damage area, the newly added road area and the newly added stockpile area.

[0029] Compared with the prior art, the present invention has the following advantages and effects:

[0030] (1) Repeated data collection at two or more different times is performed through the airborne LiDAR system. By comparing and analyzing data from multiple periods, the changes in the forest study area over time are monitored, and the dynamic changes of forest resources are accurately captured, which can reflect the long-term trends and short-term fluctuations of the forest ecosystem.

[0031] (2) After completing the canopy analysis, the elevation distribution index of the ground layer point cloud data for each period is calculated to characterize the ground undulations and human disturbance conditions, and to reveal the spatial distribution patterns and possible causes of ground changes, providing effective support for forest management and ecological protection decision-making;

[0032] (3) A dynamic change threshold is set to identify and classify significant change areas in real time based on the amplitude and rate of change in the canopy and ground layers. Combining the change characteristics of the canopy and ground layers can more comprehensively reflect the dynamic changes of forest ecosystems. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] In order to more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the following is a brief introduction to the drawings required for the specific embodiments or the description of the prior art. In all the drawings, similar elements or parts are generally identified by similar reference numerals. In the drawings, the elements or parts are not necessarily drawn according to the actual scale.

[0034] Figure 1 It is a flow chart of a method for dynamic monitoring of forest resources based on multi-temporal lidar provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0035] In order to more clearly explain the purpose, technical solutions and advantages of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all of the embodiments. The example implementation methods can be implemented in various forms and should not be understood as being limited to the examples described herein. On the contrary, these implementation methods are provided to make the present invention more comprehensive and complete, and to fully convey the concepts of the example implementation methods to those skilled in the art.

[0036] In addition, the described features, structures or characteristics may be combined in one or more embodiments in any suitable manner. In the following description, many specific details are provided to provide a full understanding of the embodiments of the present invention. However, those skilled in the art will appreciate that the technical solution of the present invention can be practiced without one or more of the specific details, or other methods, components, devices, steps, etc. may be adopted. In other cases, known methods, devices, implementations or operations are not shown or described in detail to avoid blurring various aspects of the present invention.

[0037] The block diagrams shown in the accompanying drawings are merely functional entities and do not necessarily correspond to physically independent entities. That is, these functional entities may be implemented in software form, or in one or more hardware modules or integrated circuits, or in different networks and / or processor devices and / or microcontroller devices.

[0038] The flowcharts shown in the accompanying drawings are only exemplary and do not necessarily include all the contents and operations / steps, nor must they be executed in the order described. For example, some operations / steps can be decomposed, and some operations / steps can be combined or partially combined, so the actual execution order may change according to actual conditions.

[0039] The present invention is described in detail below in conjunction with specific embodiments:

[0040] As attached Figure 1 As shown, Embodiment 1 of the present invention provides a method for dynamic monitoring of forest resources based on multi-temporal laser radar, and the specific steps of the method are as follows:

[0041] Step S1: using an airborne laser radar system, collecting data of a target forest at multiple different periods, including but not limited to performing two or more laser radar scans at different time points, to obtain multi-temporal point cloud data of the target forest, which is used to monitor changes in the target forest over a time series;

[0042] Step S2: Based on the echo characteristics of the laser radar, the point cloud of each period in the multi-temporal point cloud data is layered and processed with the echo order and echo intensity as constraints, and the point cloud is divided into a canopy layer point cloud and a ground layer point cloud to obtain canopy layer point cloud data and ground layer point cloud data. The echo characteristics include but are not limited to: the canopy layer is dominated by the first wave echo, and the ground layer is dominated by the last wave echo;

[0043] Step S3: Based on the canopy point cloud data, calculating the canopy change index of each period of the canopy point cloud, so as to characterize the dynamic change law of the target forest canopy in terms of density, height and coverage;

[0044] Step S4: Based on the crown change index, a crown change threshold is set, crown change analysis is performed on the crown layer point cloud data of each period, the area with significant crown changes is identified, and based on the combination pattern of the crown change index, the nature and cause of the significant crown change is judged, the crown change is classified, and the crown change analysis result is obtained;

[0045] Step S5: Based on the ground layer point cloud data, calculating the elevation distribution index of the ground layer point cloud in each period, so as to characterize the fluctuation and human disturbance of the target forest ground;

[0046] Step S6: Based on the elevation distribution index, a ground change threshold is set, and ground change analysis is performed on the ground layer point cloud data of each period, and areas where the ground elevation and point cloud density change significantly are identified and classified to obtain the spatial distribution pattern and possible causes of the ground change as the ground change analysis result;

[0047] Step S7: Combining the crown change analysis results and the ground change analysis results, extracting and mapping the resource change area data of the target forest, obtaining a forest resource change distribution map, and comprehensively describing the dynamic change process of the target forest resources.

[0048] Preferably, the data collection in step S1 is to input multi-temporal laser radar data, and the multi-temporal refers to performing two or more repeated data collections using an airborne laser radar system in the study area, i.e., the target forest area. Each data collection is called a period, and the data acquired in each period is a three-dimensional laser radar point cloud record of the study area at a specific time point. By comparing and analyzing the data of multiple periods, the change monitoring of the study area in the time series is realized;

[0049] Preferably, in step S1, in order to ensure the comparability of the laser radar data in different periods, each time the multi-temporal point cloud data is collected, the data consistency principle must be followed, and the data consistency principle includes but is not limited to:

[0050] (1) Using the same flight platform to maintain consistent flight attitude and route, such flight platform includes but is not limited to fixed-wing aircraft, helicopters and drones;

[0051] (2) Use the same type of LiDAR sensors, including but not limited to solid-state LiDAR and 3D imaging LiDAR, to keep the laser parameters such as wavelength, transmission power, and field of view consistent;

[0052] (3) Use the same scanning mode, including but not limited to line scanning, area scanning, and rotation scanning, to keep the point cloud distribution characteristics consistent;

[0053] (4) The same altitude and flight speed are used to ensure that the laser point clouds acquired at different times have similar point spacing and point density. The consistency of the above parameters is the prerequisite for ensuring the comparability of data from different periods.

[0054] Exemplarily, a drone is selected as a flight platform, equipped with a three-dimensional imaging lidar sensor, to carry out multi-temporal data collection on the target forest, and two data collection tasks are performed respectively to obtain high-precision laser point cloud data; after the flight, the onboard lidar sensor will synchronously record the aircraft's POS (a system used to track and monitor the position and attitude of the aircraft, Positioning and Orientation System, POS) data, including GPS time, flight attitude, track trajectory and other information; through differential positioning technology, the POS data is solved into high-precision three-dimensional coordinates and Euler angle data for subsequent laser point cloud solution; using the supporting lidar data processing software, the original lidar waveform data is solved into three-dimensional point cloud data, and each laser point contains Coordinates and reflection intensity attributes; at the same time, flight strip adjustment is performed according to the POS data, and multiple flight strips are spliced ​​to finally obtain full coverage seamless point cloud data of the target forest area, that is, the multi-temporal point cloud data.

[0055] It should be noted that after completing the multi-temporal airborne lidar data collection and preprocessing, a high-precision, high-density, full-coverage three-dimensional point cloud dataset of the target forest will be obtained in each period. These datasets are consistent in terms of spatial coordinate system, laser point attributes, distribution characteristics, etc., and the multi-temporal point cloud data will be obtained after integration.

[0056] Preferably, the layered processing in step S2 is based on the elevation distribution characteristics of the point cloud in the multi-temporal data in the vertical direction, and divides the point cloud into different layers, including but not limited to the canopy layer, shrub layer and ground layer; in this embodiment, the focus is on extracting the canopy layer point cloud and the ground layer point cloud, and by comparing the changes of the canopy layer point cloud and the ground layer point cloud in different periods, dynamic monitoring of forest resources is achieved.

[0057] Specifically, the core basis of point cloud layering processing is the echo characteristics of the laser radar. The laser radar transmits high-frequency laser pulses to the target, receives the reflected echo from the target, and calculates the propagation distance of the laser beam by measuring the time difference between transmission and reception, and then solves the three-dimensional coordinates of each echo point. When the laser beam penetrates the tree canopy, the upper layer of the tree canopy usually produces the first echo, and the lower layer of the tree canopy or the ground produces the last echo. Therefore, in a complete laser echo, the first echo mainly comes from the tree canopy layer, and the last echo mainly comes from the ground layer.

[0058] Specifically, in addition to the echo order, the canopy point cloud and the ground point cloud can be further distinguished based on the echo intensity. Since the canopy surface is rough and uneven, the laser beam is scattered and absorbed by leaves, branches, etc., and the echo intensity is relatively weak. The ground surface is relatively flat, and most of the laser beam is directly reflected back to the sensor. Therefore, the intensity of the ground echo is generally higher than that of the canopy echo. The combined use of echo order and echo intensity can more accurately extract the canopy point cloud and the ground point cloud.

[0059] Preferably, during the layered processing in step S2, a cutoff height threshold of the canopy layer needs to be set, and point clouds above the threshold are classified as canopy layer point clouds, and point clouds below the threshold are classified as ground layer point clouds; the selection of the threshold needs to comprehensively consider factors such as the tree height and terrain undulation of the target forest; in this embodiment, the cutoff height threshold is preferably 0.5 meters, that is, point clouds with an elevation greater than 0.5 meters are canopy layer point clouds, and point clouds with an elevation less than or equal to 0.5 meters are ground layer point clouds;

[0060] Specifically, after completing the above parameter settings, automatic stratification processing is performed through the lidar point cloud processing software, including: extracting the first echo point of each laser beam to obtain the first wave point cloud dominated by the canopy layer; extracting the last echo point of each laser beam, and performing ground filtering to obtain the last wave point cloud dominated by the ground layer; in the extraction process, taking 0.5 meters as the cutoff height threshold, the part of the first wave point cloud above 0.5 meters is automatically classified as the canopy layer, and the part of the last wave point cloud below 0.5 meters is classified as the ground layer; after the stratification processing is completed, the original point cloud data of each period is divided into two subsets: canopy layer point cloud data and ground layer point cloud data, among which the canopy layer point cloud data reflects the three-dimensional morphology and distribution of the forest canopy surface, including tree height, density, coverage and other attribute information; the ground layer point cloud depicts the ups and downs of the terrain and the coverage of the ground under the forest.

[0061] Preferably, in step S3, the tree crown change index includes but is not limited to: tree crown point cloud density, tree crown average height, calculated tree crown coverage and tree crown maximum height. The calculation of the tree crown change index specifically includes the following:

[0062] It should be noted that after completing the layered processing of multi-temporal LiDAR point cloud data, it is necessary to further calculate multiple statistical indicators of the canopy point cloud in each period to quantitatively describe the three-dimensional structural characteristics of the canopy. By comparing the changes in indicators in different periods, the dynamic changes in the density, height, coverage, etc. of the canopy can be revealed, providing a quantitative basis for forest resource monitoring;

[0063] Step S301, calculate the crown point cloud density, which is used to reflect the number of crown points per unit area. The higher the point cloud density, the denser the laser radar sampling of the crown surface and the more detailed the three-dimensional structure depiction. The calculation formula of the crown point cloud density is: ,in Indicates the density of tree crown point cloud, in points / square meter; Indicates the number of crown points, which can be obtained by statistically classifying the crown point cloud; It represents the area of ​​the target forest in square meters, which can be obtained by calculating the internal area of ​​the boundary vector. The crown point cloud density value of the corresponding period can be obtained by dividing the number of crown points in each period by the area of ​​the study area;

[0064] Step S302, calculate the average crown height to reflect the overall vertical elevation distribution of the crown. The larger the average crown height, the taller the trees are and the mature forest stands are. The smaller the average crown height, the younger or damaged forest stands are. The calculation formula for the average crown height is: ,in It represents the average height of the tree crown in meters; The elevation value of each crown point can be extracted from the point cloud attribute table; Represents the number of crown points; sum up the elevation values ​​of all crown points and divide by the number of points to obtain the average crown height value for each period;

[0065] Step S303, calculate the maximum crown height, that is, the maximum value of the crown point elevation, which reflects the height of the upper dominant trees in the forest stand. The calculation formula is: , that is, taking the maximum value of all crown point elevation values, represents the highest height of trees in each period; the change of the maximum crown height reflects the growth or felling of dominant trees;

[0066] Step S304, calculate the crown coverage, which means the proportion of the vertical projection area of ​​the crown to the total area of ​​the study area. The higher the crown coverage, the higher the forest density and the complete vegetation coverage. On the contrary, the forest is sparse and has more gaps. The calculation formula of the crown coverage is: ,in is the tree crown cover, expressed as a percentage; is the vertical projection area of ​​the tree crown, in square meters, which can be obtained by the convex hull or Shape extraction; represents the target forest area; dividing the tree crown vertical projection area by the target forest area to obtain the tree crown coverage in each period;

[0067] Preferably, after calculating the crown change index of the crown layer point cloud in each period in step S3, it is also necessary to calculate the change amount and change rate of the crown change index in different periods, so as to evaluate the structural dynamics of the crown within a certain time scale, wherein the change amount reflects the increase or decrease of the absolute value of the crown change index, and the change rate reflects the rate of change of the crown change index relative to the initial value;

[0068] For example, taking the tree crown point cloud density as an example, the tree crown point cloud density change calculation formula is: ,in is the density change of the crown point cloud, are the tree crown point cloud density of the latter period and the previous period respectively. The change can be positive or negative. A positive value indicates an increase in density, while a negative value indicates a decrease in density. The formula for calculating the change rate of the tree crown point cloud density is: ,in is the change rate of the crown point cloud density, expressed as a percentage; similarly, the calculation methods of the change amount and change rate of the average crown height and crown coverage are similar to those of the crown point cloud density, and are expressed as , , , )express.

[0069] Preferably, a raster data layer of the canopy point cloud data is generated to intuitively display the spatial change distribution of the canopy change index; illustratively, by using the "point to raster" tool of a comprehensive geographic spatial platform (such as ArcGIS), a raster data layer is generated based on the canopy point cloud data of each period; specifically, the canopy point cloud data is input into the tool, an appropriate interpolation algorithm (such as inverse distance weighted interpolation) is used, and a suitable spatial resolution (such as 1 meter) is set to obtain a raster data layer of indicators such as canopy point cloud density, average canopy height and canopy coverage; a "raster calculator" tool is used to perform raster algebraic operations to generate a raster of the change amount and change rate of the canopy change index in different time periods, to obtain a spatiotemporal change layer of the canopy structure, and to comprehensively display the dynamic changes of the canopy.

[0070] Preferably, based on the crown change index and the raster data layer, a quantitative assessment is performed on the dynamic processes of crown growth, loss, felling, etc. in the target forest to reveal the spatiotemporal differentiation of the three-dimensional structure of the crown.

[0071] Preferably, after setting the crown change threshold in step S4, threshold segmentation of the crown layer point cloud data is performed to divide the significance levels of the change amount grid and the change rate grid of the crown change index. The crown change threshold includes but is not limited to a point cloud density change amount threshold, a point cloud density change rate threshold, a maximum height change amount threshold, a coverage change amount threshold, a coverage change rate threshold, an average height change amount threshold and an average height change rate threshold, which are specifically as follows:

[0072] Step S401: For the tree crown point cloud density, define two tree crown change thresholds: point cloud density change threshold And the point cloud density change rate threshold , when area A in the target forest and When , there is a significant drop in the tree crown point cloud density in area A; for example, = Points / square meter, = , if the density change of the tree crown point cloud in area A in the tree crown point cloud data is Points / square meter, and the density change rate of the crown point cloud When the density decreases by more than 2 points and the decrease rate exceeds 20%, it is considered that a significant sudden drop in the crown point cloud density has occurred in the A area; the combined effect of the two crown change thresholds ensures that both the absolute and relative amounts of the sudden drop in density have reached a significant level;

[0073] Step S402: For the maximum height of the tree canopy, define a maximum height change threshold , when area A When , there is a significant drop in the maximum tree crown height in area A; for example, Meters, if the maximum height change of the tree crown in area A in the tree crown point cloud data is When the height of the tree canopy in the area A drops significantly, it is considered that the tree canopy maximum height drops significantly, indicating that the tall trees in the area have been significantly cut down or damaged.

[0074] Step S403: For the tree crown coverage, define a coverage change threshold and coverage change rate threshold , when the coverage change of area A And the coverage change rate When , there is a significant drop in canopy coverage in area A; for example, , , when the change in tree crown coverage in area A , and the crown coverage change rate When the coverage decreases by more than 5 percentage points and the rate of decrease exceeds 10%, it is considered that a significant drop in canopy coverage has occurred in area A, indicating that the forest canopy in the area has become sparse and gap-filled to a certain extent.

[0075] Step S404: For the average tree crown height, define an average height change threshold and the average height change rate threshold , when area A and , then the average canopy height in region A increased significantly; for example, 0.5 m, If the average height change of the tree crown in area A is 5%, meters, and the average crown height change rate When the average crown height increases by more than 0.5 meters and the growth rate exceeds 5%, it is judged that a significant increase in the average crown height has occurred in the A area, indicating that the trees in the area generally grow well and the crown develops upward as a whole; as above, in the crown change analysis, the positive growth of the change indicator is considered.

[0076] Preferably, the classification of crown changes in different periods in step S4 is achieved based on the setting of the crown change threshold, the crown change index of each grid unit in the grid data layer is judged and classified, the change amount and change rate of the crown change index are compared with the crown change threshold, the grid unit with significant change is found, and the type of change is judged in combination with the crown change index to obtain the crown change analysis result, which is specifically as follows;

[0077] Step S405: If area A satisfies the threshold conditions of a sudden drop in the crown point cloud density and a sudden drop in the crown coverage, and these two sudden drops are concentrated and contiguous in space, then area A is determined to be a suspected tree felling area, indicating that the trees in area A may have been artificially felled, resulting in a sudden drop in the crown point cloud density and the crown coverage, presenting a large area of ​​change patches;

[0078] Step S406: If area A satisfies the threshold conditions of a sudden drop in the maximum height of the crown and a sudden drop in the density of the crown point cloud, and these two sudden drops are concentrated and contiguous in space, then area A is determined to be a suspected tall tree felling area, indicating that the tall trees in area A may be selectively felled, resulting in a significant decrease in the height and density of the trees in the upper layer of the canopy;

[0079] Specifically, in addition to the identification of negative changes, there is also an analysis of positive changes, as follows:

[0080] Step S407: When area A satisfies both the conditions that the average crown height increases slowly and the point cloud density increases slightly, area A is determined to be a natural forest growth area, indicating that the trees in the area grow vigorously during this period, the crowns generally increase in height, and the density increases, reflecting the process of natural succession and recovery of the forest; when area A satisfies the conditions that the average height decreases slowly and the coverage decreases slowly, area A is determined to be a suspected area of ​​natural tree death or damage, indicating that the trees in the area may be affected by natural factors (such as drought, pests and diseases) or human factors (such as over-exploitation), resulting in a decrease in crown height and a sparse crown structure.

[0081] It should be noted that, as mentioned above, crown change analysis is based on crown change indicators, through the division of change thresholds and the discrimination of change patterns, to achieve qualitative identification and classification of dynamic changes in forest crowns. The crown change thresholds used in crown change analysis can be adjusted and optimized according to the actual situation of the target forest to meet the needs of different forest types and management objectives.

[0082] Preferably, the elevation distribution index in step S5, including but not limited to the average ground elevation, the lowest ground elevation, and the ground point cloud density, is obtained by performing a statistical analysis of the elevation distribution of the ground layer point cloud data of each period. The calculation of the elevation distribution index is specifically as follows:

[0083] Step S501: In the elevation range of 0-100 meters, with 1 meter as the group interval, the ground points are grouped according to the elevation value, the frequency of the points falling in each elevation interval is counted, and the elevation interval and the corresponding frequency are plotted into a histogram to intuitively display the distribution characteristics of the ground elevation; the frequency of the point reflects the distribution density of the ground point on the vertical elevation gradient, the abscissa of the histogram is the elevation interval, the ordinate is the frequency, the peak interval represents the dominant value of the ground point elevation, and the discrete degree of the histogram reflects the size of the ground elevation difference;

[0084] Step S502: Calculate the average ground elevation of the ground point cloud data in each period. The calculation formula is: ,in is the average ground elevation in meters; The elevation value of each ground point is extracted from the point cloud attribute table; is the total number of ground points; the average ground elevation of each period can be obtained by summing up the elevation values ​​of all ground points and dividing it by the total number of ground points. The average ground elevation reflects the vertical position of the entire ground and is an important parameter for characterizing the changes in ground undulation. The changes in the average elevation in different periods reflect the uplift or subsidence process of the ground.

[0085] Step S503: Calculate the lowest ground elevation of the ground point cloud data in each period. The calculation formula is: , that is, taking the minimum value of the elevation of all ground points; the lowest elevation of the point cloud reflects the lowest depression on the ground, corresponding to low terrain areas such as water bodies and valleys. The sudden change of the lowest elevation value (especially the appearance of zero value) often means that a negative elevation area has appeared on the ground, which may be related to factors such as artificial excavation and mining collapse;

[0086] Step S504: Calculate the ground point cloud density of the ground point cloud data in each period, and the calculation formula is: ,in is the ground point cloud density, in points / square meter, is the total number of ground points, is the target forest area in square meters; the total number of ground points is divided by the target forest area to obtain the ground point cloud density for each period; the ground point cloud density reflects the sampling density of the lidar on the ground, and indirectly reflects the density of the ground cover. A significant increase in the ground point cloud density indicates that there are traces of human activities on the ground, such as road construction, hardened surface, etc.;

[0087] Exemplarily, in this embodiment, the ground point cloud data can be batch processed by programming languages ​​such as Python and their scientific computing libraries (such as NumPy, SciPy), etc. Taking Python as an example, the histogram function of the extension library (NumPy) is used to count the elevation frequency distribution of ground points, set a suitable elevation range (such as 0-100 meters) and group distance (such as 1 meter), obtain the elevation histogram data, and use the histogram function of the drawing library to draw an intuitive frequency distribution image; specifically, the average function and minimum function of the extension library can also be used to quickly calculate the average elevation and minimum elevation of the ground point cloud in each period.

[0088] It should be noted that, through the calculation and analysis of the above-mentioned elevation distribution indicators, it is possible to quantitatively evaluate the spatiotemporal distribution characteristics of the target forest ground elevation and characterize the vertical structural changes of the topography.

[0089] Preferably, the ground change threshold in step S6 includes, but is not limited to: a minimum elevation change threshold, an average elevation change threshold, a ground point cloud density change threshold, and a ground point cloud density change rate threshold, which are specifically set as follows:

[0090] Step S601: For the lowest ground elevation, set the lowest elevation change threshold , which is used to determine whether there is a significant negative elevation anomaly in the lowest depression area on the ground. If the change in the lowest ground elevation in area A is , then the ground in area A has undergone a significant negative elevation change; for example, Meters, if the lowest ground elevation change in area A is , that is, when a zero-value area that is more than 1 meter lower than the previous period appears, it is considered that a significant negative elevation change has occurred in area A, indicating that the ground has obviously sunken or sunken, which may be related to factors such as surface collapse and mining settlement; the setting of the minimum elevation change threshold needs to consider the terrain undulation characteristics of the study area and the background conditions for collapse. In this embodiment, it is preferably 1 to 2 meters;

[0091] Step S602: For the average ground elevation, set the average elevation change threshold , used to determine whether the ground as a whole has a significant uplift change. If the average elevation change in area A is > , it is determined that there is a significant positive elevation change in area A; for example, suppose Meters, if the average elevation change in area A is >1 meter, that is, when the average elevation is locally raised by more than 1 meter, it is considered that a significant positive elevation change has occurred in area A; positive elevation change is often related to factors such as artificial accumulation and geological uplift. The setting of the average elevation change threshold needs to be combined with the topographic characteristics of the target forest and the intensity of human activities. In this embodiment, 0.5-1 meter is preferred;

[0092] Step S603: For the ground point cloud density, a ground point cloud density change threshold is set. And the ground point cloud density change rate threshold , which is used to determine whether the sampling density of the ground point cloud has increased significantly. If the density change of the ground point cloud in area A is > , ground point cloud density change rate > , then it is determined that the ground density in area A increases significantly; for example, let 5 points / square meter, is 50%, if the ground point cloud density change in area A >5, and the ground point cloud density change rate When the density of ground points increases by more than 5 points and the increase exceeds 50%, a significant sudden increase in density occurs in area A; the sudden increase in density indicates that the surface cover has changed significantly, including but not limited to vegetation felling and surface hardening, which enables the laser beam to penetrate the ground and generate more echo points; and The setting of needs to balance the land cover type of the study area and the technical parameters of the lidar. In this embodiment, it is preferred Set it between 3 and 5 points / square meter. Set it between 30% and 50%.

[0093] Preferably, the identification and classification of the areas where the ground elevation and point cloud density change significantly in step S6 is based on the ground change threshold, spatial analysis and pattern recognition are performed on the elevation distribution index of each period, the raster layer of each index in the elevation distribution index is compared with the corresponding ground change threshold, the change area that meets the threshold condition is extracted, and the combination relationship between the elevation distribution indexes is further analyzed to interpret the nature and cause of the change, so as to obtain the ground change analysis result, which specifically includes:

[0094] Step S604: if area A satisfies the following conditions simultaneously: the ground minimum elevation change is less than the minimum elevation change threshold, the ground point cloud density change is greater than the ground point cloud density change threshold, and the ground point cloud density change rate is greater than the ground point cloud density change rate threshold, that is, a zero value area appears in the ground minimum elevation and the ground point cloud density increases sharply, then area A is determined to be a newly added road area; road excavation will cause negative elevations to appear on the local ground, and the ground point cloud density will increase significantly due to road hardening;

[0095] Step S605: If area A satisfies the following conditions at the same time: the average elevation change is greater than the average elevation change threshold, the ground point cloud density change is greater than the ground point cloud density change threshold, and the ground point cloud density change rate is greater than the ground point cloud density change rate threshold, but the lowest ground elevation is not less than the lowest elevation change threshold, then area A is determined to be a newly added stockpiling area; stockpiling will cause the local ground to rise, forming a positive elevation anomaly, and will also increase the density of ground points, but will not change the elevation of the lowest depression;

[0096] Step S606: If the ground minimum elevation change in area A is not less than the minimum elevation change threshold, the average elevation change is not greater than the average elevation change threshold, the ground point cloud density change is not greater than the ground point cloud density change threshold, and the ground point cloud density change rate is not greater than the ground point cloud density change rate threshold, then area A is determined to be an area with no significant ground change; the topography in the area with no significant ground change remains relatively stable, has not been significantly disturbed by human beings or transformed by natural forces, and represents the normal background of the target forest;

[0097] Specifically, as mentioned above, ground change analysis is based on the calculation of ground change indicators, through the setting of ground change thresholds and the discrimination of change combination patterns, to achieve the identification and classification of areas with significant ground changes; the parameters selected for the ground change thresholds can be adjusted and optimized according to factors such as the geological and geomorphological background of the target forest and the intensity of human activities to meet the needs of different regions and application targets;

[0098] It should be noted that ground change analysis mainly targets topographic changes caused by human activities and geological processes. Some minor and gradual ground changes, such as natural weathering and soil erosion, may be difficult to accurately identify through threshold segmentation. In addition, some complex ground change situations, such as road excavation and material stacking, mining subsidence and vegetation restoration, may be difficult to fully distinguish with just a few simple change indicators.

[0099] Preferably, in step S7, extracting and mapping the resource change area data of the target forest is specifically as follows:

[0100] Step S701: spatially superimpose and perform logical operations on the change areas extracted from the crown change analysis and the change areas identified in the ground change analysis to extract the comprehensive change areas; in the superposition analysis, it is necessary to exclude known natural change factors, such as tree losses caused by natural disasters, forest pests and diseases, etc., to ensure that the extracted change areas mainly reflect the changes in forest resources caused by human activities;

[0101] Specifically, several typical comprehensive change zone types are defined as follows, and the corresponding extraction formulas and rules are given: Suspected logging area It refers to the area where the tree canopy has dropped sharply over a large area and new roads have appeared. The extraction formula is: , that is, the continuous crown drop area The spatial analysis function of GIS is used with the newly added road area to realize the fusion and reclassification of change information based on raster algebraic operations and vector overlay analysis methods. Taking ArcGIS software as an example, the tree crown sudden drop raster layer and the newly added road vector layer are first spatially superimposed, and the two layers are added using the "raster calculator" tool to derive the suspected felling area vector, and the "subsidence reduction analysis" is performed on the tree crown sudden drop area raster, that is, the known disaster-stricken area is subtracted from the tree crown sudden drop area in the raster calculator to derive the suspected damaged area vector, and finally the raster calculator is used to subtract the tree crown slow rise area, the suspected felling area and the suspected damaged area to derive the natural growth area vector, and the ground change vectors such as the newly added road and the stockpile yard are merged to derive the human interference area vector.

[0102] Step S702: After the comprehensive change zone extraction is completed, the change zone vectors are subsequently processed and mapped: the area of ​​the change zone vectors is calculated to obtain the area attribute of each change patch; the change zone vectors are attributed, and fields such as "change type" and "change time period" are added to the attribute table, and attributes are assigned according to the corresponding discrimination rules and time information; the vector layer with change attributes is superimposed on the administrative division map of the target forest, and the association between the change zone and the administrative district is realized through tools such as "spatial connection" or "intersection analysis", which is convenient for zoning statistics and mapping.

[0103] Step S703, using the Layout view of ArcGIS to perform thematic mapping of the comprehensive change zone vectors: design a color scheme according to the change type, including but not limited to: suspected logging areas are red, suspected damaged areas are yellow, natural growth areas are green, and human interference areas are purple, so as to intuitively express the spatial distribution of different change types; add administrative division base maps, change zone distribution maps, legends, compasses, scales and other mapping elements to build a complete thematic map; at the same time, load a summary statistical table of the change zone in the Layout view to display the area statistics of different change types of the target forest in each administrative area, so as to facilitate quantitative evaluation; finally, export the mapping results to general formats such as PDF and JPG to form a complete forest resource change distribution map.

[0104] It should be noted that, through the extraction and mapping of comprehensive change zones, the fusion and visualization of multi-source change information are realized, and the dynamic change pattern of forest resources in time and space is intuitively revealed; the comprehensive change zone map not only reflects the overall change trend of the forest ecosystem, such as tree growth, damage, felling, etc., but also reflects the interference intensity and spatial differentiation characteristics of human activities on the forest. This information is of great reference value for forest management departments to carry out resource monitoring, ecological protection, land planning and other work; since the extraction of comprehensive change zones involves multiple discrimination rules and threshold selection, there may be certain uncertainties and misjudgment risks in practical applications; in order to improve the reliability of comprehensive change zone extraction, in this embodiment, the change zone is verified and corrected by combining independent data sources such as high-resolution remote sensing images and field surveys, and the discrimination rules and threshold settings are optimized and adjusted according to regional characteristics and business needs to adapt to different application scenarios; in addition, in the process of extraction and mapping of change zones, it is also necessary to pay attention to factors such as quality control of data processing, matching accuracy of different data sources, and aesthetics of mapping to ensure the scientificity and practicality of the forest resource change distribution map.

[0105] In this embodiment, data of the target forest is collected at different periods to obtain multi-phase point cloud data for monitoring changes in the target forest in a time series; the point cloud of each period in the multi-phase point cloud data is layered with the echo order and the echo intensity as constraints to obtain canopy layer point cloud data and ground layer point cloud data, calculate the canopy change index and the ground change index, set the canopy change threshold and the ground change threshold, identify areas with significant changes in the canopy and the ground, classify the canopy and ground changes, obtain the canopy and ground change analysis results, and draw a distribution map of forest resource changes. A set of lidar point cloud change indicator systems including crown density, height, coverage, ground elevation, point density, etc. were proposed, which realized the quantitative analysis of forest change monitoring. Combined with change threshold segmentation, a variety of change threshold judgment criteria and change pattern recognition methods were designed, which greatly improved the efficiency and accuracy of automatic extraction of change information. The concept of comprehensive change zone was proposed, and a comprehensive extraction scheme integrating crown change and ground change was established. The dynamic change process of forest ecosystems was comprehensively characterized, and the automatic extraction, classification and mapping of forest crown and ground change information were realized, providing an important basis for forest management and ecological protection.

[0106] Those skilled in the art will readily appreciate other embodiments of the present invention after considering the specification and practicing the invention disclosed herein. This application is intended to cover any modification, use or adaptation of the present invention, which follows the general principles of the present invention and includes common knowledge or customary techniques in the art that are not disclosed by the present invention. The specification and examples are intended to be exemplary only, and the true scope and spirit of the present invention are indicated by the claims. It should be understood that the present invention is not limited to the precise structure described above and shown in the accompanying drawings, and various modifications and changes can be made without departing from the scope thereof. The scope of the present invention is limited only by the appended claims.

Claims

1. A method for dynamic monitoring of forest resources based on multi-temporal laser radar, characterized in that: include: Collecting data of the target forest at multiple different periods, including performing two or more laser radar scans at different time points, to obtain multi-temporal point cloud data of the target forest; Taking the laser radar echo order and echo intensity as constraints, performing layered processing on the point cloud of each period in the multi-temporal point cloud data, dividing the point cloud into a canopy layer point cloud and a ground layer point cloud, and obtaining canopy layer point cloud data and ground layer point cloud data; Based on the canopy point cloud data, a canopy change index of each period of the canopy point cloud is calculated, and based on the canopy change index, a canopy change threshold is set, and a canopy change analysis is performed on the canopy point cloud data of each period to identify areas where the canopy has changed significantly, and based on the combination pattern of the canopy change index, the nature and cause of the significant change in the canopy is judged to achieve classification of the canopy change and obtain a canopy change analysis result; Based on the ground layer point cloud data, calculating the elevation distribution index of the ground layer point cloud in each period, so as to characterize the fluctuation change and human disturbance status of the target forest ground; Based on the elevation distribution index, a ground change threshold is set, and a ground change analysis is performed on the ground layer point cloud data of each period. Areas with significant changes in ground elevation and ground point cloud density are identified and classified, and the spatial distribution pattern and possible causes of ground changes are obtained as the ground change analysis results.

2. The method for dynamic monitoring of forest resources based on multi-temporal laser radar according to claim 1 is characterized in that: The crown change indicators include: crown point cloud density, average crown height, crown coverage and maximum crown height; after calculating the crown change index of each period of the crown layer point cloud, the change amount and change rate of the crown change index in different periods are also calculated to evaluate the structural dynamics of the crown within a certain time scale. The change amount reflects the increase or decrease in the absolute value of the crown change index, and the change rate reflects the rate of change of the crown change index relative to the initial value.

3. The method for dynamic monitoring of forest resources based on multi-temporal laser radar according to claim 2 is characterized in that: Based on the canopy point cloud data, a raster data layer is generated to intuitively display the spatial change distribution of the canopy change index, and a change amount raster and a change rate raster of the canopy change index in different time periods are generated to obtain a spatiotemporal change layer of the canopy structure, which comprehensively displays the dynamic changes of the canopy layer; After setting the crown change threshold, threshold segmentation of the crown layer point cloud data is performed to divide the significance levels of the change amount grid and the change rate grid of the crown change index. The crown change threshold includes a point cloud density change amount threshold, a point cloud density change rate threshold, a maximum height change amount threshold, a coverage change amount threshold, a coverage change rate threshold, an average height change amount threshold, and an average height change rate threshold.

4. The method for dynamic monitoring of forest resources based on multi-temporal laser radar according to claim 3 is characterized in that: The setting of the point cloud density change amount threshold and the point cloud density change rate threshold is specifically as follows: for the tree crown point cloud density, two tree crown change thresholds are defined: point cloud density change amount threshold And the point cloud density change rate threshold , when the density change of the crown point cloud in area A of the target forest And the crown point cloud density change rate When , it is determined that a significant drop in the tree crown point cloud density has occurred in area A; The setting of the coverage change threshold and the coverage change rate threshold is specifically as follows: for the tree crown coverage, the coverage change threshold is defined as and coverage change rate threshold , when the coverage change of area A And the coverage change rate When , there was a significant drop in tree crown coverage in area A.

5. The method for dynamic monitoring of forest resources based on multi-temporal laser radar according to claim 4 is characterized in that: If area A in the target forest meets the threshold conditions of sudden drop in crown point cloud density and sudden drop in crown coverage at the same time, that is, and , and , and when these two sudden drops are concentrated and continuous in space, area A is judged to be a suspected deforestation area.

6. The method for dynamic monitoring of forest resources based on multi-temporal laser radar according to claim 5 is characterized in that: The elevation distribution indicators include: average ground elevation, minimum ground elevation, and ground point cloud density; The ground change thresholds include: a minimum elevation change threshold, an average elevation change threshold, a ground point cloud density change threshold, and a ground point cloud density change rate threshold; The setting of the minimum elevation change threshold is specifically as follows: for the lowest ground elevation, the minimum elevation change threshold is set , which is used to determine whether there is a significant negative elevation anomaly in the lowest depression area on the ground. If the change in the lowest ground elevation in area A is , it is determined that significant negative ground elevation changes have occurred in area A.

7. The method for dynamic monitoring of forest resources based on multi-temporal laser radar according to claim 6 is characterized in that: The setting of the average elevation change threshold is as follows: for the average ground elevation, the average elevation change threshold is set , used to determine whether the ground as a whole has a significant uplift change. If the average elevation change in area A is > , it is determined that a significant positive ground elevation change has occurred in area A.

8. The method for dynamic monitoring of forest resources based on multi-temporal laser radar according to claim 6 is characterized in that: The setting of the ground point cloud density change amount threshold and the ground point cloud density change rate threshold is specifically as follows: for the ground point cloud density, the ground point cloud density change amount threshold is set. And the ground point cloud density change rate threshold , which is used to determine whether the sampling density of the ground point cloud has increased significantly. If the density change of the ground point cloud in area A is > , ground point cloud density change rate > , it is determined that the ground density in area A has increased significantly.

9. The method for dynamic monitoring of forest resources based on multi-temporal laser radar according to claim 6 is characterized in that: If area A satisfies the following conditions simultaneously: the ground minimum elevation change is less than the minimum elevation change threshold, the ground point cloud density change is greater than the ground point cloud density change threshold, and the ground point cloud density change rate is greater than the ground point cloud density change rate threshold, that is, the ground minimum elevation has a zero value area and the ground point cloud density increases suddenly, then area A is determined to be a newly added road area; If area A simultaneously satisfies the following conditions: the average elevation change is greater than the average elevation change threshold, the ground point cloud density change is greater than the ground point cloud density change threshold, and the ground point cloud density change rate is greater than the ground point cloud density change rate threshold, but the lowest ground elevation is not less than the lowest elevation change threshold, then area A is determined to be a newly added stockpile area.

10. The method for dynamic monitoring of forest resources based on multi-temporal laser radar according to claim 9, characterized in that: Combining the crown change analysis results and the ground change analysis results, extracting and mapping the resource change area of ​​the target forest, and comprehensively describing the dynamic change process of the target forest resources; The resource change area includes the suspected tree felling area, the suspected tall tree felling area, the suspected tree natural death or damage area, the newly added road area and the newly added stockpile area.

Images