CN117541679B

CN117541679B - Forest canopy height mapping method and system based on sample point individual representativeness

Info

Publication number: CN117541679B
Application number: CN202410014839.4A
Authority: CN
Inventors: 洪亮; 何丽; 史正涛
Original assignee: Yunnan Normal University
Current assignee: Yunnan Normal University
Priority date: 2024-01-05
Filing date: 2024-01-05
Publication date: 2024-04-02
Anticipated expiration: 2044-01-05
Also published as: CN117541679A

Abstract

The invention provides a method and a system for mapping forest canopy height based on sample point individual representativeness, and belongs to the field of forest environment monitoring. The method comprises the following steps: firstly, determining a geographic area for drawing, and collecting geographic environment factors highly related to a forest canopy to characterize the geographic environment conditions of footprint points and unknown points; secondly, measuring the representing degree of the footprint point to the unknown point by comparing the comprehensive similarity of the geographic environment conditions between the two points; and finally, estimating the forest canopy height of the unknown points according to the representing degree of all footprint points of the area to the unknown points and quantifying the space estimation uncertainty of the forest canopy height. The method breaks through the strict requirement of the existing method on the global representativeness of the footprint points, and makes up the limitation of uncertainty in speculation caused by the fact that the representativeness of the footprint points cannot be measured, and the obtained forest canopy height speculation result has better spatial continuity and higher drawing precision.

Description

Forest canopy height mapping method and system based on sample point individual representativeness

Technical Field

The invention belongs to the field of forest environment monitoring, and particularly relates to a method and a system for mapping forest canopy height based on sample point individual representativeness.

Background

Forest is the main regulator of earth climate, forest canopy is the main producer of forest ecosystem, and scholars usually adopt the mode of highly drawing forest canopy to carry out forest degradation monitoring and forest recovery effect evaluation.

At present, two main types of methods for mapping the height of a large-scale forest canopy at home and abroad are adopted: the first method is based on a regression analysis method, the footprint point forest canopy height estimated by satellite-borne LiDAR data is used as a ground truth value, a 'best fit' regression model between the forest canopy height and geographical environment covariates (climate, topography, forest canopy spectral reflectivity and the like) is constructed by utilizing a statistical theory, and the forest canopy height value at an unknown position is estimated by utilizing the regression model, so that space continuous forest canopy height drawing is realized. Although the method is simple in principle, easy to understand and calculate; however, the method needs a large number of footprint points with good global representativeness for the forest canopy height spatial change of the research area to construct a regression model meeting the stationarity assumption, and ignores the influence of individual representativeness differences of the footprint points of the area on the forest canopy height estimation of unknown points. The second method is based on spatial interpolation of the forest canopy height map, namely, a spatial autocorrelation relation between the forest canopy height difference and the geographical environment covariate difference (such as the forest canopy spectrum, climate, terrain and the like) is established from a footprint point set by using a statistical method, the spatial autocorrelation relation is used for determining the spatial interpolation weight of the footprint point to the forest canopy height of the unknown point, and then the weighted summation is carried out on the forest canopy height values of the unknown point adjacent to the footprint point to infer the forest canopy height of the unknown point. Compared with a forest canopy height mapping method based on a regression model, the method considers the spatial autocorrelation of the differences of geographic environment covariates of unknown points and adjacent footprint points, reduces errors caused by the spatial representativeness differences of the footprint points to the unknown points, but the method still needs a large number of footprint points with good global representativeness to construct the spatial autocorrelation relation, and the representativeness of the adjacent footprint points to the unknown points is limited, so that the spatial autocorrelation relation between the unknown points in a complex terrain area and the forest canopy heights of the footprint points is difficult to accurately express.

In summary, the forest canopy height mapping methods based on regression analysis and spatial interpolation all need a large number of footprint points with good global representativeness to construct a correlation relationship of 'forest canopy height-geographic environment covariates' meeting the regional stability assumption. However, the forest canopy height is comprehensively influenced by geographical environment covariates such as soil, topography, climate, illumination and the like, and the regional forest canopy height and the geographical environment covariates have a certain synergistic relationship. However, the relationship of 'forest canopy height-geographical environment covariates' in the complex terrain area has stronger local spatial heterogeneity, is influenced by data acquisition cost and traffic accessibility, is difficult to acquire a group of typical footprint points capable of comprehensively reflecting the spatial variation of the forest canopy height of the area, and is difficult to meet strict requirements of the existing method on the number and spatial distribution of the footprint points due to the limited footprint points, so that the precision of the obtained forest canopy height spatial distribution information is difficult to ensure by utilizing the existing forest canopy height drawing method based on the footprint points, and the method is not suitable for the accurate estimation of the forest canopy height in the complex terrain area.

Disclosure of Invention

In view of the foregoing drawbacks or deficiencies of the prior art, the present invention is directed to a method and system for mapping forest canopy heights based on sample point individual representatives, wherein each footprint point is considered to be a case containing a specific "forest canopy height-geographical environment covariate" relationship, and may represent a forest canopy height of an unknown point having similar geographical environment conditions thereto, based on the assumption that the geographical environments are more similar and the forest canopy height values are more similar. The present invention refers to the "individual representativeness" of such a single footprint point as a sample point, the degree of representativeness being measured by the similarity of the geographic environment between the two points. According to the method, the forest canopy height of the area which can be represented by the sample point is estimated by the geographic environment similarity between the unknown point and the sample point, so that the strict requirement of the existing method on the global representation of the footprint point of the area is broken through, the accuracy of the forest canopy height drawing of the area with complex terrain can be effectively improved, and the limitation that the existing method cannot measure the estimation uncertainty caused by the poor representation of the footprint point can be overcome.

In order to achieve the above purpose, the embodiment of the present invention adopts the following technical scheme:

in a first aspect, an embodiment of the present invention provides a method for mapping a height of a canopy of a forest based on individual representativeness of sample points, including the following steps:

step S1, determining a geographic area for drawing, and collecting data of the geographic area related to canopy height; the related data at least comprises ICESat-2 ATL08 data and Global Ecosystem dynamic survey (Global Ecosystem DynamicsInvestigation, GEDI) data;

s2, configuring geographic environment covariates;

step S3, constructing a forest canopy height consistency model, and carrying out GEDI and ICESat-2 ATL08 facula scale forest canopy height consistency processing;

s4, calculating the feature importance of each geographic environment covariates in the quantitative characterization area forest canopy height space change by using a random forest regression algorithm based on the star LiDAR footprint point forest canopy height and the geographic environment conditions of the star LiDAR footprint point forest canopy height;

s5, respectively calculating the geographic environment similarity between each footprint point and the point to be speculated in the region by utilizing the Gaussian similarity function and the feature importance of each geographic environment covariate, and quantitatively representing the representing degree of the footprint point to be speculated by the geographic environment similarity between the footprint point and the point to be speculated;

S6, quantifying the presumption uncertainty of the point to be presumed according to the representing degree of the point to be presumed of all footprint points of the area, and setting a presumption uncertainty threshold; when the presumption uncertainty is smaller than the threshold value, presuming the forest canopy height of the point to be presumed by taking the representative degree as the forest canopy height of the weight comprehensive footprint point; when the uncertainty of the estimation is greater than or equal to a threshold value, the forest canopy height of the point to be estimated is assigned a null value;

and S7, drawing a forest canopy height spatial distribution map and a presumption uncertainty distribution map of the geographic area according to the calculated forest canopy height values of all the points to be presumed.

As a preferred embodiment of the present invention, the drawing method further includes:

and S8, evaluating the precision of the drawing.

As a preferred embodiment of the present invention, the forest canopy height-related data further includes: airborne LiDAR data, field measured data, sentinel-2 data, SRTM-DEM data, climate data and soil moisture data.

As a preferred embodiment of the present invention, the geographical environment covariates include forest canopy characteristics, topography characteristics, climate characteristics, and soil moisture response index.

As a preferred embodiment of the present invention, the topographical features include: elevation, grade, slope direction; the climate characteristics include: annual average temperature, seasonal temperature, annual average precipitation, seasonal precipitation; the forest canopy features include: annual maximum normalized vegetation index NDVI _MAX Coefficient of annual NDVI variation NDVI _CV The average spectrum variation vegetation index SVVImean, the greenness index Tcgreeneness in the thysancap transformation.

As a preferred embodiment of the present invention, the step S3 of performing GEDI and ICESat-2 att 08 spot scale forest canopy height consistency processing includes:

step S31, GEDI and ICESat-2 ATL08 data screening;

step S32, extracting overlapped footprints of GEDI and ICESat-2 satellite-borne LiDAR data;

s33, constructing a consistency model;

step S34, obtaining GEDI and ICESat-2 ATL08 data after consistency test.

As a preferred embodiment of the present invention, the geographic environment similarity in step S5 includes: forest canopy feature similarity, terrain similarity, climate similarity and soil similarity.

As a preferred embodiment of the present invention, step S5 includes:

step S51, calculating the similarity of the unknown pixels and the footprint point pixels on each geographic environment covariates by using a Gaussian similarity function;

And S52, integrating the similarity of the geographical environment covariates by utilizing the depicting capacity of the geographical environment covariates on the synergistic relationship of the forest canopy height and the geographical environment covariates of different forest vegetation areas to obtain the geographical environment similarity of the unknown pixels and the footprint pixels.

As a preferred embodiment of the present invention, the precision evaluation in step S8 includes: GEDI verification data evaluation, field actual measurement data evaluation and airborne LiDAR evaluation.

In a second aspect, an embodiment of the present invention further provides a system for mapping a canopy height of a forest based on individual representativeness of sample points, where the system includes: the device comprises a data collection module, a covariate configuration module, a data preprocessing module, a feature importance calculation module, a similarity calculation module, a height value calculation module and a drawing module; wherein,

the data collection module is used for determining a cartographic geographic area and collecting data of the geographic area and canopy height; the related data at least comprises ICESat-2 ATL08 data and global ecological system dynamic investigation GEDI data;

the covariate configuration module is used for configuring geographic environment covariates;

the data preprocessing module is used for constructing a forest canopy height consistency model and carrying out GEDI and ICESat-2 ATL08 facula scale forest canopy height consistency processing;

The feature importance calculating module is used for calculating the feature importance of each geographic environment covariates when the quantitative characterization area forest canopy height space changes by using a random forest regression algorithm based on the star-borne LiDAR footprint point forest canopy height and the geographic environment conditions of the star-borne LiDAR footprint point forest canopy height in the ICESat-2 ATL08 data;

the similarity calculation module is used for calculating the geographical environment similarity between each footprint point and the point to be speculated in the area by utilizing the Gaussian similarity function and the feature importance of each geographical environment covariate, and quantitatively representing the representing degree of the footprint point to be speculated by using the geographical environment similarity between the footprint point and the point to be speculated;

the height value calculation module is used for quantifying the presumption uncertainty of the point to be presumed according to the representing degree of all footprint points of the area to be presumed, and setting a presumption uncertainty threshold; when the presumption uncertainty is smaller than the threshold value, presuming the forest canopy height of the point to be presumed by taking the representative degree as the forest canopy height of the weight comprehensive footprint point; when the uncertainty of the estimation is greater than or equal to a threshold value, the forest canopy height of the point to be estimated is assigned a null value;

the drawing module is used for drawing a forest canopy height spatial distribution map and a presumption uncertainty distribution map of the geographic area according to the calculated forest canopy height values of all points to be presumed.

The technical scheme provided by the embodiment of the invention has the following beneficial effects:

according to the method and the system for mapping the forest canopy height based on the sample point individual representativeness, namely the IPFCHM method, the forest canopy height mapping is carried out by utilizing the geographic environment similarity between unknown points and footprint points based on GEDI and ICESat-2 data and based on the third law of geography, mapping results have good consistency with GEDI, airborne LiDAR and field verification data, compared with the existing large-scale forest canopy height mapping products with the same spatial resolution, the forest canopy height spatial estimation precision and robustness of the IPFCHM method in a complex terrain area are obviously superior to those of regression analysis (RM) and spatial interpolation (Neuralnetwork guided interpolation, NNGI) methods, and the IPFCHM method is not influenced by the number of sample points and spatial distribution pattern. The IPFCHM method provided by the invention reduces the influence of the saturation effect of the optical remote sensing data on the estimation of the forest canopy height space in the complex geographic environment area, and has high drawing precision and accuracy, and the obtained forest canopy height drawing result has better space continuity.

Of course, it is not necessary for any one product or method of practicing the invention to achieve all of the advantages set forth above at the same time.

Drawings

In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for the description of the embodiments will be briefly described below, and it is apparent that the drawings in the following description are only some embodiments of the present invention, and other drawings may be obtained according to these drawings without inventive effort for a person of ordinary skill in the art.

FIG. 1 is a flow chart of a method for mapping a height of a canopy of a forest based on a representative individual sample point provided by an embodiment of the present invention;

FIG. 2 is a plot of a consistency analysis scatter plot of the computed forest canopy height percentile parameters, including GEDI RH100, RH95, RH90, RH85, RH80, RH75, RH70, RH65, and RH60, and on-board LiDAR data, for GEDI in an embodiment of the present invention;

FIG. 3 is a scatter plot of estimated forest canopy height versus GEDI RH65 for ICESat-2 ATL08 raw data (a) and after consistency verification (b) in an embodiment of the present invention;

FIG. 4 is a graph of the Yunnan forest canopy height map based on the IPFCHM method in an embodiment of the present invention;

FIG. 5 is a graph accuracy evaluation result of the Yunnan forest canopy height of the Yunnan province based on the IPFCHM method in the embodiment of the invention; (a) GEDI verification data, (b) airborne LiDAR data, (c) field measured data;

FIG. 6 is a graph of forest canopy height mapping suitability evaluation results of IPFCHM, RM and NNGI methods in different topography index intervals, (a) RMSE, (b) MAE, and (c) BIAS in an embodiment of the present invention;

FIG. 7 is an illustration of forest canopy height mapping suitability assessment of IPFCHM, RM and NNGI methods in different slope intervals, (a) RMSE, (b) MAE, (c) BIAS, in an embodiment of the invention;

FIG. 8 is an illustration of forest canopy height mapping suitability assessment of IPFCHM, RM and NNGI methods in different terrain relief intervals, (a) RMSE, (b) MAE, (c) BIAS, in an embodiment of the invention;

FIG. 9 is a graph showing the effect of a multi-source on-satellite LiDAR data consistency check and spatial distribution pattern on forest canopy height spatial estimation accuracy in an embodiment of the present invention; (a) The accuracy comparison with airborne Lidar data, and (b) the comparison with field actual measurement data;

FIG. 10 is a graph of comparison of IPFCHM forest canopy heights based on different combinations of on-board LiDAR data in an embodiment of the present invention, comprising: (a-d) use only the ICESat-2 ATLAS footprint point, (e-h) use only the consistency-verified ICESat-2 ATLAS footprint point, (i-l) use only the GEDI footprint point, (m-p) use both the GEDI and the ICESat-2 ATLAS footprint point, (q-t) use both the GEDI and the consistency-verified ICESat-2 ATLAS footprint point;

FIG. 11 is a graph showing the effect of the number of footprint points and the spatial distribution pattern of the LiDAR on the accuracy of the height mapping of the canopy of a forest in an embodiment of the present invention, including (a) the comparison result of the LiDAR data with the onboard LiDAR and (b) the comparison result of the LiDAR data with the field actual measurement data;

FIG. 12 is a graph comparing the accuracy of the canopy height map of a forest in the IPFCHM method according to an embodiment of the present invention, including Random Forest (RF), weighted Average (WA), least restriction factor (MO), fused least restriction factor (MTF) and weighted average (CMPWA);

FIG. 13 is a graph comparing the results of the forest canopy height map of Yunnan province based on the IPFCHM method with the products of the forest canopy height map based on RM (d) and NNGI method (c) in the example of the present invention.

Detailed Description

The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. It will be apparent that the described embodiments are only some, but not all, embodiments of the invention. The components of the embodiments of the present invention generally described and illustrated in the figures herein may be arranged and designed in a wide variety of different configurations. It should be noted that, in the case of no conflict, the embodiments of the present invention and features in the embodiments may also be combined with each other.

It should be noted that: like reference numerals and letters denote like items in the following figures, and thus once an item is defined in one figure, no further definition or explanation thereof is necessary in the following figures. In the description of the present invention, the terms "first," "second," "third," "fourth," and the like are used merely to distinguish between descriptions and are not to be construed as indicating or implying relative importance.

After finding the problems, the inventor of the application conducted intensive researches on the existing forest canopy height mapping method and system. According to the third law of geography, the more similar the combination of the elements of the geographic environment between two points is under the same forest type and the same forest age, the more similar the height value of the canopy of the forest is; that is, the higher the overall similarity between LiDAR footprint points and points to be speculated on the geographic environmental elements (including spatial and non-spatial elements), the more representative an individual is speculated on the high-altitude space of the forest canopy to be speculated. Therefore, in the vegetation type range of each forest, the comprehensive similarity degree of the unknown points and the footprint points on the geographic environment element combination is utilized to infer the forest canopy height value of the unknown points, so that the strict requirement of the existing method on the global representativeness of the sample points is broken through, and the precision and the accuracy of the forest canopy height drawing of the complex geographic environment area can be improved.

It should be noted that the above drawbacks and solutions of the prior art solutions are all results obtained by the inventor after practice and careful study, and thus the discovery process of the above problems and the solutions presented below by the embodiments of the present invention for the above problems should be all contributions of the inventor to the present invention during the process of the present invention.

After the above deep analysis, the application provides a method (Individual predictive Forest Canopy Height mapping, IPFCHM) and a system for mapping the height of the forest canopy based on sample point individual representativeness based on that the similar geographical environments are, and the similar height values of the forest canopy are. Aiming at a forest coverage area with extremely complex geographic environment, the IPFCHM method takes footprint point forest canopy height values extracted from GEDI and ICESat-2 ATLAS data subjected to consistency test as ground real samples, realizes a forest canopy height map by using the IPFCHM method based on Sentinel-2, topography, climate, soil moisture and other data, compares precision and robustness with the existing 30m resolution forest canopy height map with global and national dimensions manufactured by using regression and space autocorrelation methods, and shows that the IPFCHM method can effectively utilize limited LiDAR footprint points of the complex geographic environment area to carry out forest canopy height mapping and quantify uncertainty of space speculation, has higher mapping precision and better algorithm robustness, and the precision is not influenced by the quantity and space distribution of the footprint points.

As shown in fig. 1, the method for mapping the height of the canopy of the forest based on the individual representativeness of the sample points comprises the following steps:

step S1, determining a mapped geographic area, and collecting data of the geographic area related to canopy height.

In this step, the related data at least includes ICESat-2 ATL08 data, global ecosystem dynamic survey GEDI data, and may also include airborne LiDAR data, field actual measurement data, sentinel-2 data, SRTM-DEM data, climate data and soil moisture data.

And S2, configuring the geographical environment covariates.

In this step, the geographical environment covariates include forest canopy characteristics, terrain characteristics, climate characteristics, and soil moisture response index.

And S3, constructing a forest canopy height consistency model, and carrying out GEDI and ICESat-2 ATL08 footprint scale forest canopy height consistency processing.

In this step, the consistency processing includes:

step S31, GEDI and ICESat-2 ATL08 data screening;

step S32, extracting overlapping footprints of GEDI and ICESat-2 ATL 08;

s33, constructing a consistency model;

step S34, obtaining GEDI and ICESat-2 ATL08 data after consistency test.

And S4, calculating the feature importance of each geographic environment covariates in the quantitative characterization area forest canopy height space change by using a random forest regression algorithm based on the space-borne LiDAR footprint point forest canopy height and the geographic environment conditions of the space-borne LiDAR footprint point forest canopy height.

And S5, respectively calculating the geographical environment similarity between each footprint point and the point to be speculated in the region by utilizing the Gaussian similarity function and the feature importance of each geographical environment covariate, and quantitatively representing the representing degree of the footprint point to be speculated by the geographical environment similarity between the footprint point and the point to be speculated.

In this step, the geographical environment similarity includes similarity of 12 geographical environment covariates in four aspects, which are respectively: forest canopy feature similarity, terrain similarity, climate similarity and soil similarity. Although the single on-board LiDAR footprint point has limited representatives of the area "forest canopy height-geographic environment covariate" correlation, it still contains some valuable forest canopy height spatial presumption information. For example, each footprint point contains the forest canopy height and the geographical environmental "niche" features of its location. Therefore, the application assumes that the forest canopy heights of two spatial positions with similar geographic environments are similar, and considers that each footprint point can be used as a case containing a specific relationship of forest canopy height-geographic environment covariate, and can represent an unknown point similar to the geographic environment; this representation is referred to herein as "individual representation" of footprint points, and the degree of representation can be measured by the similarity of the geographic environments between the two points. For each geographic environment covariate, calculating the similarity of the unknown point and the footprint point on the single geographic environment covariate by using a Gaussian similarity function; for the geographic environment covariate combination, the method constructs a quantitative relation of forest canopy height-geographic environment covariate by using a random forest regression algorithm, calculates weights of geographic environment covariates, and calculates the geographic environment similarity of each unknown point and footprint point by using the similarity of single geographic environment covariates and the weights of the geographic environment covariates.

S6, quantifying the presumption uncertainty of the point to be presumed according to the representing degree of the point to be presumed of all footprint points of the area, and setting a presumption uncertainty threshold; when the presumption uncertainty is smaller than the threshold value, presuming the forest canopy height of the point to be presumed by taking the representative degree as the forest canopy height of the weight comprehensive footprint point; and when the uncertainty of the estimation is greater than or equal to a threshold value, assigning the forest canopy height of the point to be estimated as a null value.

And S8, evaluating the precision of the drawing result.

In this step, the precision evaluation includes: GEDI verification data evaluation, field actual measurement data evaluation and airborne LiDAR evaluation.

The embodiment of the invention also provides a system for mapping the height of the canopy of the forest based on the individual representativeness of the sample points, which comprises: the device comprises a data collection module, a covariate configuration module, a data preprocessing module, a feature importance calculation module, a similarity calculation module, a height value calculation module and a drawing module; wherein,

In addition, it should be noted that, in this embodiment, the system for mapping the height of the canopy of the forest based on the individual representativeness of the sample point corresponds to the method for mapping the height of the canopy of the forest based on the individual representativeness of the sample point, and the description and limitation of the method are also applicable to the system and are not repeated herein.

And the method for drawing the canopy height in the selected area is carried out by taking Yunnan province in southwest China as a research area and adopting the method for drawing the canopy height based on the sample point individual representativeness.

Step S1, determining a geographic area for drawing, and collecting data related to the geographic area and canopy height.

The combination part of China and three areas of southeast Asia and south Asia at the position of Yunnan province (97 degrees 31-106 degrees 11'E, 21 degrees 8-29 degrees 15' N) belongs to the mountain area province of a plateau, the topography is declined from northwest to southeast, the topography types are complex and various, and the topography types of basins, hills, mountains, plateaus and the like are all distributed. In addition, the whole Yunnan province belongs to subtropical plateau monsoon climate, but is influenced by complex topography and topography in the region and hydrothermal condition differences of latitudes, warp directions and vertical zones, a complex and various natural geographic environment is formed, a forest ecological system with various biological species and vegetation types is inoculated, and the feasibility of the traditional forest canopy height drawing method in the complex topography region is explored.

After selecting the geographic area of the drawing, the relevant data in the area is collected.

In this example, 88 plot-scale forest canopy height data from a field survey conducted at 5 months 2023 was used to verify the validity of the IPFCHM method. In the actual measurement, the height value of the forest canopy of each plot (0.09 ha) is measured by a three-tree method, namely, the average value of the height values of the forest canopy of 3 minimum trees, 3 average trees and 3 dominant trees in the plot, namely, the arbor with the chest diameter of more than or equal to 5cm, is obtained by measuring the height value of the forest canopy of each tree through a laser range finder. The central position of each sample plot is determined by GPS positioning, and mature pure forest with low growth speed and large-area coverage is selected as a field actual measurement sample plot in the embodiment in consideration of lower GPS positioning precision under the forest and inconsistent acquisition time with other remote sensing data (for example, acquired in 2020).

Collecting 4 regions of Yunnan province together 355.51 km ² The effectiveness of the IPFCH method is verified by taking the airborne LiDAR data of the aircraft as ground truth values, all the airborne LiDAR data utilize the aerial survey/remote sensing system D2-LiDAR210 and Riegl VQ580-II in month 12 of 2020, and the aerial flying area of each airborne LiDAR is 0.42km ² To 353.67km ² And not equal. All on-board LiDAR data is processed using LiDAR360 software, including outlier processing, filtering, and CHM extraction. First, the outlier processing aims at removing high-level gross errors (such as echoes caused by clouds and low-altitude flyers) and low-level gross errors (such as echoes caused by multipath effects and very low points generated by laser range finder errors), and in this embodiment, the outlier of the airborne LiDAR data is removed by checking whether a point is more than a mean+5std threshold from its nearest neighbor 10 LiDAR points. In a preferred embodiment, the abnormal value of the regional airborne LiDAR data is removed by adopting a manual identification method. Secondly, the purpose of filtering is mainly to distinguish ground points and non-ground points (namely forest vegetation points in the embodiment) from airborne LiDAR data, and the ground points and the forest vegetation points are extracted by adopting an improved progressive irregular triangular filtering algorithm. Finally, the surface points and forest vegetation points are interpolated into a Digital Terrain Model (DTM) and a Digital Surface Model (DSM) using Kriging and IDW methods, and the difference between the DSM and DTM is taken as its corresponding CHM value. To ensure consistent resolution, the present embodiment resamples CHM data for each region to 30m spatial resolution using the bilinear method and sets the pel value to 98 th percentile forest canopy height (ghm_rh 98).

And acquiring satellite-borne LiDAR data. In this embodiment, the satellite LiDAR data includes two types of satellite LiDAR data, GEDI and ICESat-2 ATL 08.

And S2, configuring the geographical environment covariates.

In this embodiment, according to the knowledge of the relationship between the height of the canopy of the forest and the geographical environment (Simard et al 2011), 12 geographical environment covariates are selected from four geographical environment elements including climate, topography, spectrum of the canopy of the forest and soil to characterize the geographical environment characteristics related to the spatial variation of the height of the canopy of the forest, as shown in table 1.

Table 1 geographical environment covariates

The 12 geographical environment covariates include a topographical feature, a climate feature, a canopy spectral feature, and a soil moisture response index, wherein the topographical feature comprises: elevation, grade, slope direction; the climate characteristics include: annual average temperature, seasonal temperature, annual average precipitation, seasonal precipitation; the forest canopy spectral features include: NDVI _MAX 、NDVI _CV SVVImean, TCgreeness. The topographical features (grade, elevation and slope direction) were extracted from 30m resolution SRTMDEM data, grade and slope direction calculated using terrain api and SRTMDEM in GEE; the climate characteristics are extracted based on month climate data WorldClim (version 2.1) with a resolution of 1km worldwide in 1970-2000, and 4 climate characteristics of annual average temperature, seasonal temperature, average precipitation and seasonal precipitation are calculated by using a coefficient of variation method. The canopy spectral features are calculated based on the Sentinel-2A/B data of the GEE platform. In order to reduce the influence of cloud pollution and atmospheric attenuation on the height mapping of the canopy of the forest, the annual NDVI maximum value, the annual NDVI variation coefficient value (Coefficient of variation, NDVIcv, used for evaluating the stability of the annual NDVI value), the average spectrum variation vegetation index (SVVImean, spectral variabilityvegetation index) and the greenness index in the thysancap transformation are calculated by utilizing Sentinel-2A/B SR data available all year round in Yunnan province of 2020, and 4 characteristics are totally calculated as shown in the formula (3-4). The soil moisture response index (SoilMoistureResponseIndex, SMRI) is based on the soil moisture data of 10-100cm of 1km day scale estimated based on meteorological stations in 2020, and the slope is used for calculating the energy of soil moisture retention in the depth of 100cm Force to characterize the ability of the soil to provide moisture required for forest growth.

Forest vegetation type adopts the special data of the forest vegetation classification of Yunnan province with 30m resolution in 2016, and the data is the forest vegetation classification based on the difference of the physical and environmental characteristics of the forest vegetation. In this embodiment, the forest vegetation type masks are performed by using the Yunnan forest vegetation classification data to eliminate the influence of the non-forest area on the height drawing of the forest canopy. To maintain consistency of the geographical environment factors, the geographical environment covariates employed were re-projected to WGS 84N and all geographical environment covariates were re-sampled to 30m using a bilinear interpolation method.

(1)

(2)

(3)

(4)

In the formulas (1) - (4), e is regional geographical environment configuration, and m is the number of geographical environment covariates; MS is the average seasonal temperature (precipitation),and->The standard deviation and the mean of the month temperature (precipitation) are respectively. />Andstandard deviation and mean of all NDVI image synthesis in 2020.Is the standard deviation of the blue, green, red, near infrared, SWIR1 and SWIR2 bands. />Is the standard deviation of the NIR, SWIR1, SWIR2 bands.

And S3, constructing a forest canopy height consistency model, and carrying out GEDI and ICESat-2 ATL08 footprint scale forest canopy height consistency processing. The method specifically comprises the following steps:

Step S31, GEDI and ICESat-2 ATL08 data screening.

The preprocessing process of the satellite-borne GEDI data is as follows:

GEDI was launched by the United states aerospace agency at day 12 and 5 of 2018, and was installed in an International space station on-board LiDAR system with a spatial coverage of 51.6S-51.6N, a footprint diameter of about 25 meters, a ground track spacing of about 60m, a ground inter-track distance of about 600m, and containing L1-L4 four-stage data products. Research has shown that compared with the GEDIL2A data of the first edition, the GEDIL2A data of the second edition is improved in positioning accuracy and selection of a forest canopy attribute optimal algorithm, and 100 relative forest canopy height indexes provided at each footprint point are more accurate. Thus, to match the on-board LiDAR data acquisition time, the present embodiment downloads data of GEDILevel-2 Averson 2 from Yunnan 1 month to 2021 month from the NASA Earth data website, and for each footprint point, extracts a series of RH indicators including RH60, RH65, RH70, RH75, RH80, RH85, RH90, RH95, RH98, RH100 from the GEDILevel-2 Averson 2 data from Yunnan. In this embodiment, the RH index of each footprint point is a value obtained by using a GEDI default algorithm. And comparing each RH index with a forest canopy height value obtained by inversion of the airborne LiDAR, and selecting a forest canopy height percentile parameter with highest consistency as the forest canopy height value of the GEDI footprint point.

Under different weather and vegetation conditions, there may be significant differences in the ability of different beams to penetrate the forest canopy, thereby affecting the accuracy of the GEDI forest canopy height estimation. Thus, the GEDI footprint points are further filtered using thresholding to preserve the best quality GEDI footprint points. The threshold analysis specifically operates as follows: firstly deleting invalid GEDI footprint points (quality_flag=1, grade_flag=0) by using a quality evaluation index to obtain GEDI footprint points of valid waveforms; secondly, considering the signal-to-noise ratio of the waveform, the GEDI footprint can penetrate the canopy to the earth's surface; further filtering the GEDI footprint points according to sensitivity greater than or equal to 0.95; finally, in order to reduce the influence of forest-free vegetation area data and cloud coverage on the height drawing of the forest canopy in Yunnan province, GEDI footprint points, wherein the estimated height of the forest canopy is lower than 3m, and the difference between the ground elevation estimated by GEDI and the TanDEM-X elevation carried by the GEDI is larger than 50m, are removed. After the thresholding, the remaining 2235413 (223 ten thousand) GEDI footprint points are screened out.

The preprocessing process of the on-board ICESat-2 ATLAS data is as follows:

ICESat-2 ATLAS was successfully transmitted in 2018 on 9 months and 15 days, and is a new generation of satellite-borne laser radar satellites transmitted after ICESat failure in the United states, the orbit height is about 500km, and the satellite-borne laser radar satellites are circulated once every 91 days. Unlike GEDI, ICESat-2 uses photon counting laser altimeter for data acquisition, each individual photon is marked with time and location, greatly expanding detection performance. The ICESat-2 ATLAS comprises three pairs of laser beams, each pair having a spacing of about 3.3km between them, each pair having an inter-pair spacing of about 90m, and comprising a high beam and a low beam of laser light, each ground track providing an overlapping footprint spaced about 0.7m along the track and about 17.5m in diameter. In this example, ICESat-2ATL08 data from Yunnan, 10 months of 2018 to 3 months of 2021 was downloaded from an ice and snow data center in the United states, and the downloaded data provided various canopy height, terrain elevation and RH related indicators along the 100m segment of the ground track. However, while the ICESat-2 ATLAS sensor improves spatial coverage, its high detection sensitivity makes it susceptible to solar background noise. First, to ensure accuracy of the ICESat-2ATL08 derived forest canopy height data, three weak beams and daytime footprint points are excluded because these footprint points are relatively low in accuracy in forest canopy height estimation; secondly, the uncertainty of forest canopy height estimation is larger (h_copy_uncertaitry=3.4028235e+37) and the footprint of canopy height estimation anomaly (h_copy >160 m) are further filtered; finally, the footprint points of ICESat-2 estimated floor elevation and SRTMDEM elevation differences greater than 50m and closed_flag_atm <2 are typically removed under the influence of cloud coverage. The remaining 1736774 (173 ten thousand) ICESat-2ATL08 footprint points were screened out.

Step S32, extracting the overlapped footprints of the GEDI and ICESat-2 satellite borne LiDAR data.

And step S33, constructing a consistency model.

In this embodiment, the ICESat-2ATL08 and GEDI data are obtained in approximately the same time, and the two data are fused, so that not only can the density of the forest canopy height sample points be increased, but also geographic position complementation can be realized, the randomness of the spatial distribution of the regional LiDAR footprint points can be increased, and a reliable data source is provided for forest canopy height mapping. However, the two kinds of satellite-borne LiDAR data have obvious differences in the aspects of forest detection mechanism, digital recording mode, data distribution mode and the like, so that the forest canopy height precision obtained by inversion of footprint dimensions of the two kinds of satellite-borne LiDAR data is inconsistent. Compared to the GEDI data, each pulse of ICESat-2ATL08 only receives a small amount of signal photons, resulting in lower accuracy of estimation of forest canopy height in dense forests. In order to integrate the data of the GEDI mode and the ICESat-2ATL08 mode to improve the drawing precision of the forest canopy height, a proposed forest canopy height consistency model is adopted to correct the precision of the forest canopy height inverted by the ICESat-2ATL08 data in a footprint scale. The model is divided into two key steps; firstly, acquiring overlapping footprint points of GEDI and ICESat-2ATL08 data; secondly, taking the height of the forest canopy of the footprint point estimated by GEDI as a response variable (GEDI RH 65), taking the characteristic parameter of ICESat-2ATL08 as an interpretation variable, and respectively establishing a forest canopy height consistency model of different forest vegetation type areas in the overlapping footprint point set of GEDI and ICESat-2ATL08 by combining a random forest and a stepwise linear regression algorithm as shown in a table 2; and the consistency model is used for correcting the forest canopy height value of all ICESat-2ATL08 footprint points in Yunnan province.

TABLE 2 ICESat-2 ATL08 characterization parameters

Step S34, obtaining GEDI and ICESat-2 ATL08 data after consistency test.

In this step, the geographic environment similarity measures unknown points(/>=1, 2,3, …, k) and on-board LiDAR footprint points(/>Similarity between =1, 2,3, …, k) at the geographical environment covariate scale and sample scale, from the environment vector, +.>And->Grid pixels respectively representing 30m resolution:

(5)

in the formula (5), the amino acid sequence of the compound,representing unknown point->And LiDAR footprint point->Geographical environmental similarity of->And->The two geographic positions of the unknown point and the sample point are respectively corresponding +.>Values of individual geographical environment covariates, +. >(/>=1, 2,3, …, m); e (-) and P (-) represent the similarity function of the unknown point and the footprint point on a single geographic environment variable, respectively, and the comprehensive similarity function of the geographic environments of the unknown point and the footprint point.

And S51, calculating the similarity of the unknown pixels and footprint point pixels in each geographic environment covariate. E (·) is used to describe unknown pointsAnd footprint point->In geographical environment covariates->As shown in equation (6), as a function of the change in similarity with the change in feature distance between two points. IPFCHM utilizes Gaussian similarity function (equations 6 and 7)And calculating the similarity of the single geographic environment covariates.

(6)

(7)/>

(8)

In the formulae (6) to (8),is unknown->And footprint point->In geographical environment covariates->The degree of similarity in terms of the degree of similarity,is the geographical environment covariates of the whole area +.>Standard deviation of>Is footprint point->And all unknown points->(=1, 2,3, …, k) in the geographical environment covariates +.>Square root of the average feature distance of (c).

The function P (-) is to establish a relation between the height of the forest canopy and the covariate of the geographical environment by utilizing the constraint structure and the characteristics of the covariate of the geographical environment to the height of the forest canopy, and establish a comprehensive function by utilizing the relation and E (-)A function of the similarity of the covariates of the geographic environments as shown in equation (9). The IPFCHM utilizes a random forest regression algorithm to construct a relation of forest canopy height-geographical environment covariates, calculates the feature importance of each geographical environment covariates, and utilizes the feature importance and the geographical environment covariate scale similarity to calculate unknown points ++ >And all footprint points in the area->Is a geographic environmental similarity of (c).

(9)

In the formula (9), the amino acid sequence of the compound,is unknown->And footprint point->Comprehensive similarity in geographical environment covariates, < ->Andrepresenting unknown points->And footprint point->Is a combination of all geographical environment covariates.

And S52, integrating the similarity of the covariates of the geographic environments to obtain the overall environmental similarity of the unknown pixels and the footprint pixels.

And (4) comprehensively representing the similarity of the geographic environment according to the relative importance of various geographic environment covariates in influencing the spatial distribution pattern of the height of the forest canopy. Therefore, the IPFCHM method adopts a random forest regression method to respectively determine the characteristic importance of each geographical environment covariate in each forest vegetation type area to the regional forest canopy height drawing, and adopts a weighted summation method to calculate the geographical environment similarity of each pixel to be presumed and the sample point according to the characteristic importance.

The forest canopy height spatial presumption uncertainty of each point to be presumed is inversely proportional to its similarity to the existing LiDAR footprint geographic environment. If the geographical environment similarity of unknown points to footprint points is poor, then the existing set of footprint points is insufficient to represent the "forest canopy height-geographical environment covariate" correlation of the location, and using these footprint points to predict the forest canopy height of the location would lead to a large uncertainty. Thus, the uncertainty estimate in the IPFCHM is an existing footprint point set to unknown point forest canopy height spatial speculation reliability assessment. The calculation formula is shown as formula (10):

(10)

wherein,refers to the point to be speculated +.>Is a forest canopy height spatial presumption uncertainty,refers to the LiDAR footprint points and the points to be speculated +.>Is a maximum value of the geographic environmental similarity.

If the uncertainty of the estimate of the unknown point is too high, specifying that none of the existing footprint points in the region represent well the "forest canopy height-geographical environment covariate" relationship for that location, using these footprint points to estimate the forest canopy height of the unknown point will result in a large uncertainty. Therefore, the IPFCHM selects a footprint point set with higher representativeness to the unknown point by setting an uncertainty threshold, takes the similarity of the geographic environment between the two points as a weight, synthesizes the forest canopy height values of all footprint points by using a weighted summation method, and further presumes the forest canopy height of the unknown point (formula 11).

(11)

In the formula (11), the amino acid sequence of the compound,is unknown->Forest canopy height estimate of +.>Is the sample dot->And unknown point->Is similar to the geographic environment of (1)Degree (f)>Forest canopy height value estimated as LiDAR footprint point, < >>Is a threshold for similarity segmentation of geographic environment, and the embodiment selects unknown points +.>And all footprint points->And the average value of the similarity of the comprehensive geographic environment is used as a segmentation threshold value.

And S8, evaluating the precision of the drawing.

To evaluate the effectiveness of the IPFCHM method in mapping the forest canopy height in a complex geographic environment area, the present embodiment uses three different scale validation data sets to quantitatively evaluate the accuracy of mapping the forest canopy height in Yunnan province based on the IPFCHM method, including 355.51 km ² Is a real-time measurement of forest canopy height data in the field of 88 sample plot scales (30 m x 30 m), and 10% of GEDI footprint scale forest canopy height data (223541). In order to maintain consistency, root Mean Square Error (RMSE), mean absolute error (meanabsolute error, MAE), mean BIAS (BIAS) three statistical parameters were calculated using each validation dataset and the forest canopy height mapping result of this embodiment, respectively, the calculation formula of which is shown in formula (12-14):

(12)

(13)

(14)

In the formulas (12-14), n is the number of verification samples,is a verification sample->Corresponding IPFCHM method generates Yunnan forest canopy height space presumption value, </i >>Forest canopy height values measured for footprint points, on-board LiDAR, and field plots were validated for GEDI.

The IPFCHM gets rid of the dependence on the regional average state 'forest canopy height-geographical environment covariate' correlation relation to satisfy the stationarity assumption, the forest canopy height spatial presumption precision and the robustness in the complex geographical environment region are obviously superior to those of the existing regression and spatial autocorrelation methods, and the IPFCHM method is not influenced by the number of sample points and the spatial distribution pattern. However, under different geographical environmental conditions, there is a certain difference in stress factors of forest canopy height, for example, in tropical rainforest areas, precipitation, temperature and illumination are major factors affecting forest canopy height in the area; in dry and hot valley areas, precipitation and soil are the main stress factors of forest canopy height in the area.

Comparing the above-mentioned Yunnan province forest canopy height map manufactured by the IPFCHM method of the present embodiment with the existing forest canopy height map with 30m resolution of the whole world and the whole country manufactured by regression and space autocorrelation methods, and comparing the precision, the robustness and the like of the map of the same area, the result is as follows:

Because the method for measuring forest canopy height by GEDI and ICESat-2 ATLAS is different from that of airborne LiDAR, in the embodiment, CHM obtained by inversion of airborne LiDAR data is used as a ground truth value, and the forest canopy height percentile parameter derived from the GEDI data is compared with the CHM generated by the airborne LiDAR data respectively, so that the forest canopy height percentile parameter with the highest consistency is selected as the forest canopy height value derived from the GEDI. As can be seen from fig. 2, as the forest canopy height percentile decreases, the consistency between GEDI and on-board LiDAR extraction CHM tends to increase and decrease. When the forest height percentile parameter is RH65, the RMSE, MAE, BIAS value between the two is minimum, 7.314m, 5.769m and-0.0165 m respectively. Thus, the present embodiment takes the GEDI RH65 as the forest canopy height value for the GEDI footprint point. In addition, FIG. 3 shows the difference in forest canopy height values obtained by inversion of the GEDI and ICESat-2 ATL08 data before and after consistency verification. As can be seen from fig. 3, the consistency test shows that the ICESat-2 atl08 forest canopy height value has good consistency with GEDIRH65 (r2=0.719, rmse= 3.171m, mae=2.011), and can be used together in the mapping process. Thus, the present embodiment uses GEDI RH65 and consistency checked ICESat-2 ATL08 data to represent the forest canopy height values generated by the on-board LiDAR footprint points; and resamples these footprint points to 30m resolution.

Forest canopy height spatial distribution map with 30m resolution in 2020 of Yunnan province is estimated based on IPFCHM method. Because of the high density of the satellite borne LiDAR footprint points (49.75/km 2) in Yunnan province, the height of the forest canopy in the whole Yunnan province is directly presumed to cause memory overflow. Therefore, according to the repetition rate of 10%, all the geographical environment elements in the whole Yunnan province are cut into 7224 image blocks of 400 multiplied by 400 by adopting a sliding window cutting method, for each image block, the forest canopy height of each pixel in the image block is estimated by adopting an IPFCHM method, and finally all the image blocks are spliced to obtain a forest canopy height spatial distribution map of the Yunnan province, fig. 7. From the whole area scale, the Yunnan forest canopy height is generally subjected to a normal distribution with an average value of 13.797m and a standard deviation of 2.902m (fig. 4). From the aspect of geographic distribution, the height spatial distribution pattern of the forest canopy in Yunnan province has obvious spatial heterogeneity under the comprehensive influence of hydrothermal condition difference and complex geographic environment, and the height of the forest canopy in the northwest, the Yunnan and the Yunnan areas is obviously higher than that of the regions of the northeast, the middle and the Yunnan areas.

And (3) precision verification:

the IPFCHM method-based forest canopy height spatial estimation result of Yunnan province has good consistency relation with three verification data sets (RMSE: 3.4205-5.4556 m, MAE: 2.9846-4.5034 m, BIAS: -0.0011-0.798 m), which shows that the IPFCHM method is suitable for forest canopy height mapping in a complex geographical environment area of Yunnan province (figure 5). As can be seen from fig. 5 (a), the forest canopy height (bias= -0.0011 m) was slightly underestimated in forests with a height greater than 25m, as compared to the geni verification data. In the evaluation results using the onboard LiDAR data (fig. 5 (b)), the mapping results based on the IPFCHM method slightly overestimate the measured forest canopy height (bias=0.017 m) in a forest with a height of less than 25 m. Comparing the results of this example with the field plot measurements (fig. 5 (c)), the results of this example overestimate the measured canopy height (bias= 0.7980 m) in forests with a forest canopy height of less than 15 m.

Comparing the drawing results of the embodiment with the drawing precision of the Chinese forest canopy height product manufactured based on the spatial autocorrelation idea and the global forest canopy height product manufactured based on the regression analysis method in the Yunnan province, the drawing precision of the embodiment is found to be remarkably improved (table 3). In order to quantitatively evaluate the drawing precision of the products, the embodiment evaluates the drawing precision of the two products in the forest canopy height of Yunnan province according to GEDI footprint point screening conditions and verification sample occupation ratio, and respectively verifies the drawing precision of the two methods in the Yunnan province by using the airborne LiDAR and field measured data of the embodiment, wherein the precision verification results are shown in Table 3. First, from the on-board LiDAR footprint point verification results, compared with the GEDI verification data of the NNGI method, the MAE, the RMSE and the BIAS of the IPFCHM method are respectively reduced by 1.4529m,1.8831m and 2.0414m; MAE, RMSE and BIAS of the method of this example were reduced by 0.0469m, 0.3984m and 0.6549m, respectively, compared to the GEDI validation data of RM method. Second, compared to the on-board LiDAR data of this example, the RMSE, MAE, and BIAS of the IPFCHM method were reduced by 71.02% and 28.22%,70.88% and 28.62%,6.0385m, and 2.9524m, respectively, compared to the NNGI and RM methods. Compared to the field measured samples of this example, RMSE, MAE and BIAS for the IPFCHM method were reduced by 82.66% and 57.91%,91.71% and 39.88%,3.2376m and 1.777m, respectively, compared to NNGI and RM methods. In conclusion, the IPFCHM method is obviously higher in the forest canopy height drawing precision in Yunnan province than the NNGI and RM methods.

Table 3 shows the comparison of drawing accuracy with other large-scale forest canopy height drawing products

In order to verify the applicability of the IPFCHM method in a complex geographical environment area, the embodiment compares the difference of forest canopy height spatial estimation precision of the IPFCHM, NNGI and RM methods in different forest type areas, gradient areas, topography relief areas and topography index areas respectively. Wherein the verification samples all use 10% GEDI data corresponding to each product.

The topography index comprehensively reflects the change characteristics of the elevation and the gradient in the research area, and can be used for representing the complexity of the geographic environment. Therefore, in this embodiment, the accuracy of the forest canopy height mapping of the IPFCHM, NNGI and RM methods in different topography index areas of the Yunnan province is evaluated by using the gendi verification data, and the evaluation result is shown in fig. 6. Overall, the applicability of the IPFCHM method to forest canopy height mapping of complex geographical environment areas is significantly higher than that of NNGI and RM methods. From the forest canopy height space presumption precision change trend, the forest canopy height drawing precision and the topography index are in negative correlation relation based on IPFCHM, NNGI and RM methods. However, since the NNGI and RM methods are both to apply the "forest canopy height-geographical environment covariate" correlation obtained from large-area-scale statistics to local-scale forest canopy height spatial estimation, the model has poor local applicability, so that the variation amplitude of the forest canopy height mapping accuracy in different topography index areas is large (the variation amplitudes of RMSE and MAE of NNGI are 2.4348 times and 2.1546 times of that of IPFCHM method respectively, and the variation amplitudes of RMSE and MAE of RM are 1.4669 times and 1.8117 times of that of IPFCHM respectively), which makes it difficult to meet the requirement of complex geographical environment forest canopy height mapping. The IPFCHM method adopts the similarity of the geographic environments of the unknown points and the footprint points to infer the height value of the forest canopy of the unknown points, gets rid of dependence on the assumption of the stability of the related relation of the regional forest canopy height-geographic environment covariates, and has better robustness on the forest canopy height drawing of the complicated geographic environment region.

Fig. 7 shows the difference in the accuracy of the IPFCHM, NNGI and RM methods for the forest canopy height mapping in different slope intervals in the Yunnan province. Firstly, from the difference of forest canopy height space presumption precision, the IPFCHM method has higher forest canopy height drawing precision in all gradient ranges of Yunnan province than NNGI and RM methods. Secondly, from the trend of the drawing precision of the canopy of the forest along with the change trend of the gradient, the drawing precision of the canopy height products of the three forests all show a descending trend along with the increase of the gradient, but the descending speeds of the three forests have larger difference. The drawing accuracy of IPFCHM was reduced by the smallest extent (rmse=2.2402 m, mae= 1.4969m, absolute bias= 1.8894 m), followed by NNGI (rmse= 5.6492, mae= 3.9366, absolute bias= 5.5609 m) and RM (rmse=6.8553 m, mae= 5.2202m, absolute bias= 5.497 m). Therefore, the gradient change has less influence on the IPFCHM method, and the forest canopy height spatial estimation robustness of the method is better than NNGI and RM in a complex geographic environment area.

The relief is the difference between the highest point elevation and the lowest point elevation in a particular area, the greater the relief, the greater the complexity of the terrain. In the embodiment, the mean value variable point method and the maximum height difference method are adopted to obtain the topographic relief degree of Yunnan province, the difference of the forest canopy height drawing precision of the IPFCHM, NNGI and RM methods in different topographic relief degree areas is evaluated based on GEDI data, and the verification result is shown in figure 8. First, from the perspective of forest canopy height mapping accuracy, the mapping accuracy of all terrain relief areas of the IPFCHM method is higher than that of NNGI and RM. Secondly, from the trend of drawing precision change, as the relief degree of the topography increases, the forest canopy height drawing precision based on the IPFCHM, NNGI and RM methods shows a trend of decreasing. However, because of the difference of the conditions assumed by the algorithm, the drawing precision variation amplitude of the complex geographic environment area is greatly different. The variation amplitude of the forest canopy height drawing precision based on the IPFCHM method with the increase of the complexity of the geographical environment is the smallest (rmse= 1.4424m, mae= 0.9778m, absolute bias= 1.2889), followed by RM (rmse=4.9723 m, mae=6.6903 m, bias= 2.9831 m) and NNGI (rmse=7.3206m, mae=5.4808 m, bias= 8.1986 m). In conclusion, the IPFCHM has good portability in different topographic relief areas, and has higher applicability to forest canopy height mapping in complex geographic environment areas.

In the drawing process, the influence of IPFCHM parameter setting on drawing precision is considered. FIG. 9 shows the results of accuracy verification of forest canopy height mapping with on-board LiDAR (FIG. 9 (a)) and field measured data (FIG. 9 (b)) based on a combination of different on-board LiDAR data. As can be seen from fig. 10, compared with the original ICESat-2 ATL08 data, the ATL08 data after consistency test significantly improves the forest canopy height drawing precision of the complex geographic environment region, and compared with the on-board and field verification data, the RMSE is reduced by 0.9m and 1.8501m respectively. Meanwhile, as can be seen from fig. 10, the forest canopy height mapping result based on the ATL08 data after consistency test better accords with the actual situation that the height of the valley forest canopy is greater than the mountain top. Secondly, compared with the forest canopy height spatial estimation using GEDI and ATL08 data alone, the forest canopy height mapping result integrating the two footprint points is smoother, which further proves the conclusion that the stripe effect of the forest canopy height spatial estimation can be reduced by integrating the GEDI and ATL08 data. However, the accuracy of the mapping was slightly lower than that of the forest canopy height mapping using GEDI data alone due to the difference in ATL08 and GEDI footprint scale forest canopy height inversion accuracy (FIG. 9), and there was a sudden change in the spatial distribution pattern of the forest canopy height (Zhu et al 2022). As can be seen from fig. 9 and 10, the fused GEDI and consistency-checked ATL08 data obtained the best forest canopy height mapping results in Yunnan province, and the RMSE values in the on-board and sample areas were reduced by 0.1544m and 0.1285m, respectively, compared to the forest canopy height mapping results of the fused GEDI and ATL08 data. The reason is probably that the ATL08 footprint point forest canopy height data after consistency test has good consistency (R2=0.79) with GEDI data, so that the influence of the difference of GEDI and ATL08 on the inversion precision of the forest canopy height on the large-scale forest canopy height drawing is greatly reduced.

In the drawing process, the influence of the number of the footprint points of the satellite LiDAR and the spatial distribution pattern of the number of the footprint points on the IPFCHM method is also considered. In this embodiment, the effect of the number of footprint points and the spatial distribution pattern thereof on the IPFCHM method is analyzed by comparing the height drawing precision of the forest canopy in the Yunnan province based on 10-100% of the satellite LiDAR data (fig. 11). As can be seen from fig. 11, the cloud south province forest canopy height mapping accuracy based on 30% footprint points is highest (mae= 4.4486m, rmse= 5.4112 m) compared to CHM data generated from on-board LiDAR data; forest canopy height mapping accuracy based on 100% footprint points is lowest (mae= 4.5034m, rmse= 5.4556 m); the amplitude of variation of MAE and RMSE values based on the 10% -100% footprint point forest canopy height mapping results were 4.476 (+ -0.0274) and 5.4334 (+ -0.0222), respectively. Compared with field measured data, the cloud-south forest canopy height mapping accuracy based on 100% of the satellite-borne LiDAR footprint points is highest (MAE= 2.9846m, RMSE= 3.4205 m), and the forest canopy height spatial estimation accuracy based on 10% of the satellite-borne LiDAR footprint points is lowest (MAE= 3.3725m, RMSE= 3.8174 m); the range of variation based on MAE, RMSE, BIAS values of 10% -100% is 3.1559 (+ -0.2127) and 3.619 (+ -0.1985) respectively. In conclusion, the number of the footprint points of the satellite LiDAR and the spatial distribution pattern of the footprint points have no obvious influence on the IPFCHM method.

In the drawing process, the influence of the comprehensive similarity measurement method of the geographic environment on the forest canopy height drawing result is also considered.

Fig. 12 shows the difference between the RF-based geographical environment similarity calculation method and the minimum limiting factor method, the weighted average method, and the forest canopy height mapping accuracy combining the weighted average method and the minimum limiting factor method in the IPFCHM method, and the mapping accuracy based on the RF method is found to be significantly improved. As can be seen from fig. 12 (a), the MAE of the RF-based geographical environment similarity calculation method is reduced by 0.1905m, 0.2904m, and 0.0467m, respectively, as compared to the minimum limiting factor method, the weighted average method, and the integrated minimum limiting factor method and the weighted average method; RMSE was reduced by 0.2394m, 0.4388m and 0.0064m, respectively; BIAS was reduced by 0.0069m, 0.3219m and 0.0881m, respectively. The significant improvement of the forest canopy height drawing precision proves the effectiveness of the IPFCHM method in the forest canopy height drawing of complex geographic environment areas. As can be seen from fig. 12 (b), the consistency of the yunnan province forest canopy height generated using RF-calculated geographical environment comprehensive similarity with the field measured sample is highest (rmse=2.985, mae= 3.421m, bias=0.798 m); the reason may be that the RF method considers the applicability of the covariates of the geographic environments to the spatial estimation of the forest canopy height of the forest type areas, reducing the influence of the complexity of the geographic environments on the spatial estimation of the forest canopy height of the areas.

Compared with 3 verification data, the forest canopy height drawing RMSE based on the IPFCHM method is 3.4205-5.4556 m, the MAE is 2.9846-4.5034 m, and the BIAS is-0.0011-0.798 m. These numbers approach the previous large scale forest canopy height mapping study (RMSE: 3.3-6.1 m) based on-board LiDAR. However, in view of the difference in spatial resolution (30 m,500m or 1 km), the present example compares the results of this example to the pixel scale of global and chinese forest canopy height map products made using RM and NNGI methods. On the pixel scale, the result of the embodiment is well matched with the existing large-scale forest canopy height space estimation product overall, the forest canopy height difference overall accords with normal distribution, the average values are 5.323m and 7.097m respectively, and the standard values are 4.321m and 5.135m respectively (fig. 13 (c) and fig. 13 (d)). From the overall spatial distribution pattern, the height of the forest canopy of the results of this example is slightly higher overall than the other two products, especially in the Xishuang Banna and Yangtze river areas. The reason for this phenomenon may be that NNGI and RM methods are "forest canopy height-geographical environment covariate" correlations extracted from all LiDAR footprint points in the study, and apply them to other areas, and the relationships obtained in this way are correlations of the average states in the area footprint point set, affected by sample acquisition cost and traffic accessibility, and it is difficult to meet the needs of forest canopy height mapping in complex geographical environment areas.

According to the technical scheme, the forest canopy height mapping method based on sample point individual representativeness provided by the embodiment of the invention is based on GEDI and ICESat-2 data and is based on the third law of geography, the geographical environment similarity between unknown points and footprint points is utilized for carrying out forest canopy height mapping, mapping results have good consistency with GEDI, airborne LiDAR and field verification data, compared with the existing large-scale forest canopy height mapping products with the same spatial resolution, the forest canopy height spatial estimation precision and robustness of the IPFCHM method in a complex terrain area are obviously superior to those of NNGI and RM methods, and the IPFCHM method is not influenced by the number of sample points and spatial distribution pattern. The IPFCHM method provided by the invention reduces the influence of the saturation effect of the optical remote sensing data on the estimation of the forest canopy height space in the complex geographic environment area, and has high drawing precision and accuracy, and the obtained forest canopy height drawing result has better space continuity.

The above description is only of the preferred embodiments of the present invention and the description of the technical principles applied is not intended to limit the scope of the invention as claimed, but merely represents the preferred embodiments of the present invention. It will be appreciated by persons skilled in the art that the scope of the invention referred to in the present invention is not limited to the specific combinations of the technical features described above, but also covers other technical features formed by any combination of the technical features described above or their equivalents without departing from the inventive concept. All other embodiments, which can be made by a person skilled in the art without making any inventive effort, are intended to be within the scope of the present invention.

Claims

1. A forest canopy height mapping method based on sample point individual representativeness is characterized by comprising the following steps:

step S1, determining a cartographic geographical area, and collecting data of the geographical area related to forest canopy height; the related data at least comprises ICESat-2 ATL08 data and global ecological system dynamic investigation GEDI data;

s2, configuring geographic environment covariates; the geographical environment covariates comprise forest canopy characteristics, terrain characteristics, climate characteristics and soil moisture response indexes;the topographical features include: elevation, grade, slope direction; the climate characteristics include: annual average temperature, seasonal temperature, annual average precipitation, seasonal precipitation; the forest canopy features include: annual maximum normalized vegetation index NDVI _MAX Coefficient of annual NDVI variation NDVI _CV The average spectrum variation vegetation index SVVImean, the greenness index TCgreeneess in the thysancap transformation; the soil moisture response index SMRIs are soil moisture data of 10-100cm in 1km daily scale estimated based on meteorological stations, and the capability of soil for keeping moisture in a depth of 100cm is calculated by adopting a slope so as to characterize the capability of the soil for providing moisture required by forest growth;

s4, calculating the feature importance of each geographic environment covariates in quantitative characterization area forest canopy height space change by using a random forest regression algorithm based on the ICESat-2 ATL08 data and the spaceborne LiDAR footprint point forest canopy height in the GEDI data and the geographic environment conditions of the ICESat-2 ATL08 data and the GEDI data;

2. A method of mapping forest canopy height based on sample point individual representatives as claimed in claim 1, further comprising:

and S8, evaluating the precision of the drawing.

3. A method of mapping forest canopy height based on sample point individual representatives as claimed in claim 1 or 2, wherein the forest canopy height related data further comprises: airborne LiDAR data, field measured data, sentinel-2 data, SRTM-DEM data, climate data and soil moisture data.

4. The method for mapping the height of the canopy of the forest based on the individual representativeness of the sample points according to claim 1 or 2, wherein the step S3 of performing the consistency processing of the height of the canopy of the forest of GEDI and the ICESat-2 att 08 facula scale forest comprises the following steps:

step S31, GEDI and ICESat-2 ATL08 data screening;

s33, constructing a consistency model;

step S34, obtaining GEDI and ICESat-2 ATL08 data after consistency test.

5. The method for mapping the height of the canopy of the forest based on the individual representativeness of the sample points according to claim 1 or 2, wherein the similarity of the geographic environment in step S5 comprises: forest canopy feature similarity, terrain similarity, climate similarity and soil similarity.

6. A method for mapping forest canopy height based on sample point individual representativeness as set forth in claim 1 or 2, wherein step S5 includes:

7. The method for mapping forest canopy height based on sample point individual representativeness as set forth in claim 2, wherein the precision evaluation in step S8 includes: GEDI verification data evaluation, field actual measurement data evaluation and airborne LiDAR evaluation.

8. A system for mapping forest canopy height based on sample point individual representativeness, the system comprising: the device comprises a data collection module, a covariate configuration module, a data preprocessing module, a feature importance calculation module, a similarity calculation module, a height value calculation module and a drawing module; wherein,

The data collection module is used for determining a cartographic geographic area and collecting data of the geographic area and forest canopy height; the related data at least comprises ICESat-2 ATL08 data and global ecological system dynamic investigation GEDI data;

the covariate configuration module is used for configuring geographic environment covariates; the geographical environment covariates comprise forest canopy characteristics, terrain characteristics, climate characteristics and soil moisture response indexes; the topographical features include: elevation, grade, slope direction; the climate characteristics include: annual average temperature, seasonal temperature, annual average precipitation, seasonal precipitation; the forest canopy features include: annual maximum normalized vegetation index NDVI _MAX Coefficient of annual NDVI variation NDVI _CV The average spectrum variation vegetation index SVVImean, the greenness index TCgreeneess in the thysancap transformation; the soil moisture response index SMRIs are soil moisture data of 10-100cm in 1km daily scale estimated based on meteorological stations, and the capability of soil for keeping moisture in a depth of 100cm is calculated by adopting a slope so as to characterize the capability of the soil for providing moisture required by forest growth;