Forest canopy height extraction method based on space-borne photon counting lidar

Abstract

The present application relates to a forest canopy height extraction method based on a spaceborne photon counting lidar, and belongs to the technical field of photon point cloud data processing.The present application comprises the following steps: firstly, accurately dividing signal and noise data of photon point cloud data; then, based on the separated signal data, fitting a ground curve by using an iterative filtering algorithm of local terrain, and removing pseudo-ground points by using mutation detection of a gradient descent method and decomposition and reconstruction functions of EMD, so as to accurately identify ground point and non-ground point data; finally, extracting canopy top points from non-ground points by using a dynamic binning statistical method, and calculating canopy height according to the height difference between the canopy top points and the ground points.The present application realizes accurate segmentation of signals and noises of different terrain scenes and determines a noise threshold; further removes pseudo-ground points and adjusts canopy top points under the premise of ensuring that ground signals are not distorted, and greatly improves tree height extraction precision.

Description

TECHNICAL FIELD

[0001] The present application relates to a forest canopy height extraction method based on a spaceborne photon counting lidar, and belongs to the technical field of photon point cloud data processing. BACKGROUND

[0002] In recent years, the greenhouse effect and global warming are increasingly severe, and the forest ecosystem has received widespread attention from the international community due to its stable carbon sequestration capacity. Therefore, the carbon distribution and carbon sink estimation of the forest ecosystem has substantial significance for studying the process of global climate change. Among them, the forest canopy height, as an important part of the vertical structure parameters of the forest, plays an indispensable role in the estimation of forest carbon storage. As a new type of spaceborne multi-beam micro-pulse photon counting lidar, ICESat-2 / ATLAS can provide narrow-band elevation profile point clouds, and through the accurate identification of canopy and ground photons, higher precision forest vertical structure parameters can be extracted [7] . ATLAS reduces the spot diameter and improves the along-track sampling density, and to some extent, reduces the influence of terrain compared with ICESat-1 full waveform data. However, the pulse signal of ICESat-2 photon data itself is weak, and it is easily disturbed by solar background noise and other interference when returning the signal, forming dense and uniform scene noise. Therefore, how to accurately denoise the photon point cloud data and extract higher precision forest canopy height under different terrain conditions has become the focus and difficulty of the research.

[0003] Many of the current photon point cloud data processing methods are based on ICESat-2 airborne simulation test data, such as Multiple Altimeter Beam Experimental Lidar (MABEL), Slope Imaging Multi-polarization Photon-counting Lidar (SIMPLE), Sigma Space and MATLAS data. Xia Shaobo et al. used local distance statistics and least squares curve fitting to denoise and filter photon point cloud data, but this method is suitable for high signal-to-noise ratio data (i.e. less noisy scenes), and it is difficult to remove noise signals close to the ground but below the ground, resulting in corresponding filtering errors affecting the accuracy of tree height estimation. Tang's voxel-based spatial filtering algorithm (reference: Tang H, Swatantran A, Barrett T, et al. Voxel-Based Spatial Filtering Method for Canopy Height Retrieval from Airborne Single-Photon Lidar [J]. Remote Sensing, 2016, 8(9)), also only adapts to high signal-to-noise ratio photon data. Popescu's window-based clustering statistical filtering denoising algorithm (reference: Popescu SC, Zhou T, Nelson R, et al. Photon counting LiDAR: An adaptive ground and canopy height retrieval algorithm for ICESat-2 data [J]. Remote Sensing of Environment, 2018, 208:154-170), which has made great progress in processing low signal-to-noise ratio data, uses dynamic binning statistics to identify ground points and canopy points when extracting canopy height; but in the case of high canopy coverage and sparse ground return signals, it is easy to misclassify false ground points and cause filtering errors.To solve the above problems, Nie Sheng proposed the DBSCAN density clustering denoising algorithm with buffer and mirror strategy and elliptical search (reference: Nie Sheng, Wang Cheng, Xi Xiaohuan, Luo Shezhou, Li Guoyuan, Tian Jinyan, Wang Hongtao. Estimating the vegetation canopy height using micro-pulse photon-counting LiDAR data. [J]. Optics express, 2018, 26(10)), combined with the progressive irregular triangle network encryption method to iteratively identify ground points and percentile statistics to obtain the crown point. This method effectively solves the edge effect and filtering error problem, but needs to be improved in removing the noise inside the crown layer. Xiaoxiao Zhu gradually filters out the noise signal by height statistical histogram coarse denoising, multi-directional elliptical search clustering statistical fine denoising and interval statistical final denoising (reference: Xiaoxiao Zhu, Sheng Nie, Cheng Wang, Xiaohuan Xi, Zhenyue Hu. (2018) A Ground Elevation and Vegetation Height Retrieval Algorithm Using Micro-Pulse Photon-Counting Lidar Data. Remote Sensing 10:12, pages 1962); After denoising, according to the local density statistics and EMD decomposition and reconstruction, the ground points and non-ground points are classified, and finally the crown top is identified from the non-ground points by combining the percentile statistics method. This algorithm effectively removes noise and filters out pseudo-ground points, improving the accuracy of crown height extraction. At the same time, based on the spatial characteristics of the point cloud in different directions by elliptical search statistics, it is found that the profile point cloud is more concentrated in the horizontal direction, and then Zhu Xiaoxiao proposed an improved OPTICS that can also extract high-precision signals (reference: Zhu X, Nie S, Wang C, et al. A Noise Removal Algorithm Based on OPTICS for Photon-Counting LiDAR Data [J]. IEEE Geoscience and Remote Sensing Letters, 2020).Qin Lei improved the DBSCAN algorithm and proposed an outlier detection method (reference: Qin Lei, Xing Yanqiu, Huang Jiapeng, et al. ICESat-2 airborne experimental photon cloud data adaptive denoising and classification algorithm [J]. Remote Sensing, 2020, 24(12): 1476-1487), which realizes the denoising of the crown, intracoronary and underground layers. However, the algorithm still has many uncertain parameter values, and the denoising lacks adaptability for different signal-to-noise ratios, which may cause errors in photon point cloud type recognition. In the above profile point cloud crown height extraction research, most scholars obtain ground photon points and crown top photon points after photon classification, and fit the corresponding terrain surface and crown layer surface, and then define the forest height by the height difference. In addition, some studies also directly use the crown top photon to subtract the corresponding terrain surface to realize the extraction of vegetation height. However, the ICESat-2 / ATLAS light spot in the forest coverage area only returns a few signal photons, which are not only scattered, but also have uncertain vertical positions in the crown layer. Therefore, the received crown top photon may not be located at the true tree crown top, which leads to the forest height extracted by the above method being usually lower than the true height. Existing research shows that the vegetation height based on photon data has a high correlation with the 95% percentile height of conventional airborne LiDAR point cloud. This view has good reference significance for higher precision crown top photon recognition.

[0004] Due to the limited MABEL airborne experimental data and high collection cost, and the great difference between the data and ICESat-2 data, the ICESat-2 currently uses the point cloud processing algorithm of DRAGANN and surface finding (Surface Finding) to process the original photon point cloud. The specific theory can be referred to the Algorithm Theoretical Basis Document (ATBD) of ICESat-2 / ATL08 product data basic theory algorithm. On this basis, the ATL08 product data has been publicly released, but the algorithm is not mature enough in terms of parameter adaptation, and the details still need to be further improved. Directly applying the MABEL airborne experimental data processing algorithm to the spaceborne ICESat-2 / ATLAS data still has many problems. Therefore, the present application improves the point cloud processing algorithm in the ATBD, so that it has strong adaptability and can still obtain high-precision forest parameters under complex terrain conditions. SUMMARY

[0005] In view of the problems existing in the prior art mentioned above, the present application provides a forest crown height extraction method based on spaceborne photon counting lidar, which can accurately remove photon point cloud noise and extract high-precision forest crown height.

[0006] The technical scheme of the present application is: a forest canopy height extraction method based on a spaceborne photon counting laser radar, first, accurate division of signal and noise data of photon point cloud data is performed; then based on the separated signal data, an iterative filtering algorithm of local terrain is used to fit the ground curve, and the gradient descent method is used to detect the mutation and the decomposition and reconstruction function of EMD to remove the pseudo-ground points; finally, a dynamic binning statistical method is used to extract the canopy points from the non-ground points, and the canopy height is calculated according to the height difference between the canopy points and the ground points.

[0007] As a further scheme of the present application, the specific steps of the method are as follows:

[0008] Step1, data preprocessing: coordinate projection conversion and segmentation processing are performed on the original point cloud to form an elevation profile point cloud with a suitable span;

[0009] Step2, point cloud denoising and classification: noise, ground, canopy and canopy top photons are identified through point cloud precise denoising and filtering classification;

[0010] Step3, canopy height extraction: the canopy height is calculated by calculating the height difference between the canopy points and the ground points.

[0011] As a further scheme of the present application, the Step1 includes:

[0012] Step1.1, the original three-dimensional point cloud data based on ICESat-2 / ATL03 is converted into elevation profile point cloud data through coordinate projection;

[0013] Step1.2, then the strip profile is segmented according to a span of 10km, each segmented unit is called a Seg-Km processing segment; at the same time, in order to ensure the accuracy of the photon point coordinates matching when the canopy height extracted by the present application is compared with the ICESat-2 / ATL08 product, the corresponding ATL08 data is also processed in Seg-Km; in addition, in order to solve the problem of edge effect exhibited by the denoising method based on local statistical parameters, the statistical parameter here is photon density, which is usually processed in the denoising of the Seg-Km processing segment in combination with the moving window method, and the binning span is set to 3.4km, at the same time, along the horizontal orbit, 200m buffer zones are set before and after each bin.

[0014] As a further scheme of the present application, in the Step2, the point cloud denoising includes the following specific steps:

[0015] Step2.1.1, debugging filter parameters and point-by-point local density statistics: after data testing, it is found that the filter parameter P is positively correlated with the signal-to-noise ratio SNR, and then a relationship model is constructed, such as formula (1), to realize the self-adaptation of the filter parameter P; then a K-D tree index structure is constructed for all photons in the profile point cloud, and the photon density is calculated point by point after searching the radius combined with formula (2), to form a local density histogram;

[0016] P = Int (2.96273 * SNR + 2.24325) (1)

[0017] In the formula, P is the number of neighborhood points, SNR is the rough estimated signal-to-noise ratio, by counting the elevation frequency histogram of the profile point cloud, then taking the median of the frequency as the threshold, the frequency below the threshold is noise, otherwise it is signal, and the mean value is obtained and divided to obtain SNR, Int is the integer function;

[0018]

[0019] In the formula, R is the adaptive search radius, P is the number of neighborhood points, N Total is the total number of photons under each moving window in Seg_Km processing section, wherein the span of the moving window is 200m+3.4km+200m, and Int is the integer function;

[0020] Step2.1.2, Gaussian decomposition and parameter optimization, double-peak Gaussian parameter screening: the histogram curve is filtered and decomposed by Gaussian low-pass filter, then the parameters of the Gaussian components are optimized by Levenberg-Marquardt L-M, and the signal Gaussian and noise Gaussian are finally screened by the following conditions: the leftmost peak is marked as noise Gaussian, and other Gaussian components are sorted according to α i ×σ i , and the first one that satisfies μ i -μ1≥2σ1 is marked as signal Gaussian;

[0021] Step2.1.3, accurate extraction of noise threshold: after obtaining the noise Gaussian and signal Gaussian parameters, the double-peak feature is reconstructed by formula (3) and the only wave trough point is extracted as the signal-to-noise threshold by using second-order derivative;

[0022]

[0023] In the formula, (α i ,μ i ,σ i ) are the optimized noise Gaussian and signal Gaussian parameters, i.e. amplitude, mean value and standard deviation.

[0024] As a further aspect of the present invention, Step 2, the point cloud classification includes the following specific steps:

[0025] Step 2.2.1: Define an iterative filter synthesized by the median filtering algorithm and the moving average algorithm, and since it is based on local terrain conditions, the moving window span of the filter will be set.

[0026] Step 2.2.2: Detrending processing is performed on the Seg-Km processing segment. First, the original signal after median filtering is denoised under the elevation threshold. Then, the retained signal is subjected to conformal segmentation, cubic Hermite interpolation and iterative filtering to extract the terrain trend line. Finally, the difference between the denoised profile point cloud and the terrain trend line is calculated to achieve detrending.

[0027] Step 2.2.3: Using all detrended data as the initial value and combining iterative filters to continuously update the fitting point input, the upper and lower bounds of potential ground points are obtained. The ground points within the interval are refined by iterative filters. At the same time, the ground points are interpolated to the original resolution by linear interpolation to obtain the ground curve. EMD-DISPO filters and gradient methods are used to further remove pseudo ground photons in different cases in the ground curve.

[0028] Step 2.2.4: Photons with percentiles in the range of [96, 99] are identified as top canopy photons using the binning statistical method; the binning span of the profile point cloud is usually set to 20m.

[0029] As a further aspect of the present invention, in Step 2.2.3, a method combining a DISPO filter and an Empirical Mode Decomposition (EMD) spatiotemporal scale filter is used to remove pseudo-ground points, as shown in formula (4).

[0030] The EMD algorithm decomposes the data in the time domain, and then uses time-scale denoising or threshold denoising to denoise and reconstruct the decomposition results. In this process, time-scale denoising is prone to causing distortion of the original signal, and threshold denoising is difficult to solve the problem of adaptive threshold calculation. The introduction of the DISPO algorithm first filters the obvious noise in the decomposed signal and then reconstructs it, which can solve the above problems. The original signal and reconstructed signal functions are shown below.

[0031]

[0032] In the formula, H(x) and H1(x) represent the original signal and the reconstructed signal, respectively, and imf i1 / i2(x) represents a single-component IMF function, r(x) represents the trend term, and N represents the total number of IMF components formed during the EMD decomposition of the ground curve. In the equation H(x), i1 represents the label of each IMF component. When reconstructing the ground curve in EMD, imf1(x) is discarded, and the remaining IMF component indices in [2, N / 2] are filtered by DISPO, while those in [N / 2, N] are not filtered. In the equation H1(x), i2 represents the label of each IMF component that needs to be filtered by DISPO, while j represents the label of each IMF component that does not need to be filtered.

[0033] As a further aspect of the present invention, in Step 2.2.3, when there are two land types, forest land and non-forest land, in the Seg_Km data processing segment, the density of ground photons under the forest canopy and the density of ground noise photons under the non-forest land are similar. When fitting the ground curve, it is easy to fit the ground noise photons under the non-forest land to the ground, forming abrupt changes in terrain. EMD-DISPO is not suitable for processing similar abrupt change signals. The gradient descent method is used to detect the amplitude and location of the ground point abrupt change interval, and the corresponding translation sequence is calculated based on the gradient information and the difference is made with the original signal to obtain a more accurate ground curve.

[0034] As a further aspect of the present invention, Step 3, the extraction of canopy height includes the following specific steps:

[0035] Step 3.1: The top of the canopy is slightly lower than the interior of the canopy. Therefore, before calculating the canopy height, the photons extracted from the top of the canopy in Step 2.2.4 need to be further filtered, a KD Tree spatial index is constructed for them, the number of neighborhood points within a 15m local range is counted, and the photon point cloud that deviates from the range [3,10] is re-marked as noise, and the rest are the final canopy top photons;

[0036] Step 3.2: Using the final crown photon extracted in Step 3.1 as the interpolation point, the conformal piecewise cubic Hermite interpolation algorithm is used to interpolate the ground curve obtained in Step 2.2.3, thereby obtaining the ground photon at the corresponding position in the ground curve;

[0037] Step 3.3: Finally, calculate the height difference between the canopy photons and the interpolated ground photons to extract the canopy height.

[0038] The beneficial effects of this invention are:

[0039] 1. This invention improves the point cloud processing algorithms of DRAGANN and Surface Finding in ICESat-2 / ATBD, realizes parameter adaptation in the process of photon point cloud processing of elevation profile, and improves the accuracy of forest canopy height extraction.

[0040] 2. In the process of preprocessing profile point cloud data to remove photon point cloud noise, this invention does not use a fixed value for the filter parameter P. Instead, it constructs a correlation model between parameter P and signal-to-noise ratio to achieve adaptive P values ​​under different geographical conditions, thereby improving the universality of the denoising algorithm.

[0041] 3. The TS-SCABR filter proposed in this invention not only has adaptive search capability, but also can accurately segment signals and noise in different terrain scenarios and determine the noise threshold.

[0042] 4. This invention further filters out false ground points and adjusts crown vertex identification while ensuring that the ground signal is not distorted, which greatly improves the tree height extraction accuracy. The tree height extraction accuracy is increased from 85.2% (ICESat-2 / ATL08 tree height extraction accuracy) to 96.8%, which is of great significance for the high-precision inversion of forest height on a large scale in the future. Attached Figure Description

[0043] Figure 1 This is a schematic diagram illustrating the working principle of the ICESat-2 platform of the present invention;

[0044] Figure 2 The geographical distribution of the test area in this invention;

[0045] Figure 3 This is a flowchart of the experimental process in this invention;

[0046] Figure 4 This is a schematic diagram of the ICESat-2 / ATL03 photon point cloud segmentation processing rules in this invention;

[0047] Figure 5 This is a schematic diagram of the threshold extraction results of the original DRAGANN under high signal-to-noise ratio in this invention;

[0048] Figure 6 This is a schematic diagram of the threshold extraction results of TS-SCABR under high signal-to-noise ratio in this invention;

[0049] Figure 7 This is a schematic diagram of the original DRAGANN threshold extraction results under low signal-to-noise ratio in this invention;

[0050] Figure 8 This is a schematic diagram of the threshold extraction results of TS-SCABR under low signal-to-noise ratio in this invention;

[0051] Figure 9 This is a schematic diagram of the point cloud denoising results under low signal-to-noise ratio in this invention;

[0052] Figure 10 This is a schematic diagram of the point cloud denoising results under high signal-to-noise ratio in this invention;

[0053] Figure 11 This invention uses gradient descent to assist in removing terrain abrupt changes in the translation sequence.

[0054] Figure 12 This is a schematic diagram of the original iterative moving surface filtering effect in this invention;

[0055] Figure 13 This is a schematic diagram of point cloud classification results under high signal-to-noise ratio and low signal-to-noise ratio conditions in this invention.

[0056] Figure 14 This is a schematic diagram illustrating the accuracy verification of the estimated tree height in this invention. Detailed Implementation

[0057] Example 1: This example uses ICESat-2 as a successor to ICESat-1, equipped with an Advanced Topographic Laser Altimeter System (ATLAS). (Reference) Figure 1 It is known that ATLAS employs Micro-Pulse Photon Counting Lidar (MPCL) technology. Guided by the cross-track channel, it emits six laser pulses at a repetition frequency of 10 kHz, with a duration of 1 ns. Along the Reference Ground Track (RGT), it acquires three sets of 17m diameter strong and weak light spot stripes, each set with an inter-track spacing of approximately 3.3 km and a spot spacing of approximately 0.7m along the track direction. This results in six ground tracks acquired with the RGT as the center line, with an inter-track spacing of approximately 90m within each set. ATLAS emits laser pulses towards the target, which are then scattered and reflected by ground objects and received by a satellite receiver. The system then counts and processes the echo photon events to extract the true distance data from background noise and dark counts. Currently, I-SIPS (ICESat-2 Science Investigator-led Processing System) provides 21 data products (ATL00 to ATL21) across four levels, from Level-0 to Level-3. The main research of this invention is based on the Level-2 data product ATL03. ATL08 is a Level-3A product generated based on ATL03, which is the global land vegetation height, and will be used in this invention for comparison and inversion of canopy height. All data can be downloaded free of charge from the Ice and Snow Data Center (https: / / icesat-2.gsfc.nasa.gov / icesat-2-data / );

[0058] like Figures 1-14As shown, the method for extracting forest canopy height based on spaceborne photon counting lidar includes the following specific steps:

[0059] Step 1: Data preprocessing: Perform coordinate projection transformation and segmentation on the original point cloud to form an elevation profile point cloud with a suitable span.

[0060] The specific steps of Step 1 are as follows:

[0061] Step 1.1: Based on ATL03, the original 3D point cloud (i.e., horizontal distance along the rail, vertical distance along the rail, and 3D elevation) data is transformed into elevation profile point cloud data through coordinate projection.

[0062] Step 1.2: The strip profile is then segmented according to a 10km span, with each segment unit called a Seg-Km processing segment. Simultaneously, to ensure accurate matching of photon point coordinates between the extracted canopy height and the ICESat-2 / ATL08 products during precision comparison, the corresponding ATL08 data is also processed using Seg-Km. Furthermore, to address the edge effect issues exhibited by density clustering denoising methods based on local statistical parameters (here, photon density), a moving window method is typically incorporated into the denoising of the Seg-Km processing segments. The span of each bin is set to 3.4km, and a 200m buffer zone is established before and after each bin along the transverse track.

[0063] Step 2, Point Cloud Denoising and Classification: Noise, ground, canopy and canopy apex photons are identified through fine point cloud denoising and filtering classification;

[0064] This invention improves the photon point cloud processing algorithms of DRAGANN and Surface Finding in ICESat-2 / ATBD, realizes parameter adaptation in the process of elevation profile photon point cloud processing, and improves the accuracy of forest canopy height extraction.

[0065] The noise removal process for the original ICESat-2 / AIL08 data product used the Differential Regressive and Gaussian Adaptive Nearest Neighbor (DRAGANN) algorithm. This algorithm first constructs a KD Tree spatial index, introduces the number of neighboring points P to calculate the search radius R, and then performs a point-by-point range search to generate a local density histogram of the signal. Typically, in photon point cloud data, the local density of noisy photons is lower than that of signal photons. Especially after airborne simulation experiments at MABEL, it was shown that a high and narrow noise peak forms on the left side of the histogram, while a low and wide signal peak forms on the right side.

[0066]

[0067] In the formula, R is the adaptive search radius, P is the number of neighborhood points (empirically set to 20), and N... Total This represents the total number of photons within the moving window in the Seg-Km processing segment.

[0068] After generating the histogram, the Savitzky-Golay filtering algorithm is used to smooth it, generating a smooth curve with an approximate "bimodal" distribution. Then, the noise peak and signal peak are identified by bimodal Gaussian fitting (the formula below), and the denoising threshold is identified by the bimodal threshold segmentation method.

[0069]

[0070] In the formula, (α) i ,μ i ,σ i The parameters are the optimized noise Gaussian and signal Gaussian parameters, namely amplitude, mean, and standard deviation.

[0071] In processing ICESat-2 / ATL03 data, this invention discovered the following problems with the DRAGANN algorithm: 1) The DRAGANN algorithm uses a fixed P-value to construct a local density histogram, making it easy to confuse noise and signal peaks, and difficult to accurately determine the denoising threshold; 2) Under different spatiotemporal conditions, photon point cloud signal-to-noise ratios often differ, which further leads to different overall shapes of the local density histogram. Therefore, in addition to the previously summarized histogram shape of "high and narrow noise peaks, low and wide signal peaks," this invention newly summarizes a histogram shape of "low and narrow noise peaks, wide and high signal peaks." In the latter case, due to the complex conditions of the ground features, most signals exhibit a "multi-peak" distribution. If the bimodal Gaussian fitting in the DRAGANN algorithm is used, it is difficult to accurately fit the noise peaks and signal peaks.

[0072] To address the aforementioned problems, this invention proposes a threshold segmentation based on spatial clustering and bimodal reconstruction (TS-SCABR) denoising algorithm based on the original theory of DRAGANN. The denoising algorithm includes the following specific steps:

[0073] Step 2.1.1, Adjusting filter parameters and performing local density statistics point by point: Adaptively obtain the number of neighborhood points P; Considering that the adaptability of parameter P under different natural conditions (i.e., terrain undulation and solar radiation) greatly affects the accuracy of local density statistics, this invention uses 450 samples for testing. The results show that parameter P is positively correlated with signal-to-noise ratio (SNR), and a relationship model is constructed as shown in formula (3) to achieve adaptive filter parameter P. Then, a KD Tree index structure is constructed for all photons in the profile point cloud. After calculating the search radius in combination with formula (4), the photon density is statistically analyzed point by point based on this structure to form a local density histogram.

[0074] P=Int(2.96273*SNR+2.24325) (3)

[0075] Where P is the number of neighborhood points, SNR is a rough estimate of the signal-to-noise ratio, and SNR is obtained by statistically analyzing the elevation frequency histogram of the profile point cloud, then using the median frequency as a threshold. Frequencies below this threshold are considered noise, and those above are considered signals. The mean of each frequency is calculated and then divided to obtain the SNR. Int is the floor function. This invention aims to achieve dynamic P values ​​for different regions by constructing a correlation model between the parameter P and the signal-to-noise ratio, thereby improving the universality of the denoising algorithm.

[0076]

[0077] In the formula, R is the adaptive search radius, P is the number of neighborhood points, and N is the number of neighboring points. Total The Seg_Km function processes the total number of photons within each moving window in the segment, where the span of the moving window is 200m + 3.4km + 200m, and Int is the floor function.

[0078] Step 2.1.2, Gaussian Decomposition and Parameter Optimization, Bimodal Gaussian Parameter Selection: The histogram curve is filtered using a Gaussian low-pass filter and then subjected to Gaussian decomposition. The Gaussian components are then optimized using the Levenberg-Marquardt LM algorithm. The signal Gaussian and noise Gaussian components are finally selected based on the following criteria: the leftmost peak is marked as noise Gaussian, and other Gaussian components are selected according to α. i ×σ i Sort and the first one satisfies μ i -μ1≥2σ1 is marked as a Gaussian signal;

[0079] Step 2.1.3 Accurate extraction of noise threshold: After obtaining the noise Gaussian and signal Gaussian parameters, the bimodal features are reconstructed by combining formula (5) and the second derivative is used to extract the unique valley point as the signal-to-noise threshold.

[0080]

[0081] In the formula, (α)i ,μ i ,σ i The parameters are the optimized noise Gaussian and signal Gaussian parameters, namely amplitude, mean, and standard deviation.

[0082] After achieving fine denoising of photon point clouds, the ICESat-2 basic theoretical algorithm document uses the basic principle of surface detection for point cloud classification.

[0083] Step 2, point cloud classification includes the following specific steps:

[0084] Step 2.2.1: Define an iterative filter synthesized by the median filtering algorithm and the moving average algorithm, and since it is based on local terrain conditions, the moving window span of the filter will be set.

[0085] Step 2.2.2: Detrending processing is performed on the Seg-Km segment. First, the original signal after median filtering is denoised at a specific elevation threshold. Then, the retained signal is subjected to conformal segmentation, cubic Hermite interpolation, and iterative filtering to extract the terrain trend line. Finally, the difference between the denoised profile point cloud and the terrain trend line is calculated to achieve detrending.

[0086] Step 2.2.3: Using all detrended data as initial values ​​and combining iterative filters, the fitted point input is continuously updated to obtain the upper and lower bounds of potential ground points. Ground points within the interval are refined using iterative filters, and simultaneously, linear interpolation is used to interpolate the ground points to the original resolution to obtain the ground curve. Furthermore, this invention employs an EMD-DISPO filter and gradient method to further eliminate pseudo-ground photons in different situations within the ground curve.

[0087] Step 2.2.4: Photons with percentiles in the range of [96, 99] are identified as top canopy photons using the binning statistical method; the binning span of the profile point cloud is usually set to 20m.

[0088] As a further aspect of the present invention, in areas with high vegetation coverage, laser pulses may not be able to penetrate the vegetation canopy. In this case, only a small number of signal photons can be returned from the ground surface, while there are also a lot of noise signals near the ground surface. If the noise signals are identified as ground photons, it will greatly affect the fitting accuracy of the ground curve; in Step 2.2.3, a method combining a Digital smoothing polynomial D (ISPO) filter and an Empirical Mode Decomposition (EMD) spatiotemporal scale filter is used to remove pseudo ground points, as shown in formula (6);

[0089] The EMD algorithm decomposes the data in the time domain, and then uses time-scale denoising or threshold denoising to denoise and reconstruct the decomposition results. In this process, time-scale denoising is prone to distorting the original signal, and threshold denoising is difficult to solve the problem of adaptive threshold calculation. Introducing the DISPO algorithm to filter the obvious noise in the decomposed signal before reconstruction can solve the above problems. The original signal and reconstructed signal functions are shown below;

[0090]

[0091] In the formula, H(x) and H1(x) represent the original signal and the reconstructed signal, respectively, and imf i1 / i2 (x) represents a single-component IMF function, r(x) represents the trend term, and N represents the total number of IMF components formed during the EMD decomposition of the ground curve. In the equation H(x), i1 represents the label of each IMF component. When reconstructing the ground curve in EMD, imf1(x) is discarded, and the remaining IMF component indices in [2, N / 2] are filtered by DISPO, while those in [N / 2, N] are not filtered. In the equation H1(x), i2 represents the label of each IMF component that needs to be filtered by DISPO, while j represents the label of each IMF component that does not need to be filtered.

[0092] As a further aspect of the present invention, in Step 2.2.3, when there are two land types, forest land and non-forest land, in the Seg_Km data processing segment, the density of ground photons under the forest canopy and the density of ground noise photons under the non-forest land are similar. When fitting the ground curve, it is easy to fit the ground noise photons under the non-forest land as the ground, forming abrupt changes in terrain. EMD-DISPO is not suitable for processing similar abrupt change signals. The gradient method is used to detect the amplitude and location of the ground point abrupt change interval, and the corresponding translation sequence is calculated based on the gradient information and the difference is made with the original signal to obtain a more accurate ground curve.

[0093] Step 3, Canopy Height Extraction: The canopy height is calculated by calculating the elevation difference between the canopy apex and the ground point.

[0094] Step 3, the extraction of canopy height, includes the following specific steps:

[0095] Step 3.1: The top of the canopy is usually slightly lower than the interior of the canopy. Therefore, before calculating the canopy height, the photons extracted from the top of the canopy in Step 2.2.4 need to be further filtered, a KD Tree spatial index is constructed for them, the number of neighborhood points within a 15m local range is counted, and the photon point cloud that deviates from the range [3,10] is re-marked as noise, and the rest are the final top photons.

[0096] Step 3.2: Using the final crown photon extracted in Step 3.1 as the interpolation point, the conformal piecewise cubic Hermite interpolation algorithm is used to interpolate the ground curve obtained in Step 2.2.3, thereby obtaining the ground photon at the corresponding position in the ground curve;

[0097] Step 3.3: Finally, calculate the height difference between the canopy photons and the interpolated ground photons to extract the canopy height.

[0098] This invention utilizes measured tree height data under different terrain conditions to verify the canopy height extraction method of ICESat-2. Part of the measured data comes from the Canopy Height Model (CHM) product provided by the Airborne Observation Platform (AOP) operated by the National Ecological Observatory Network (NEON) of the US National Science Foundation, with a spatial resolution of 1 m. CHM samples are distributed across five NEON sites in Alabama, Massachusetts, Georgia, and California, USA, where the terrain is relatively gentle (Table 1). The other part of the measured data consists of 25 sample plots in Shangri-La City, Yunnan Province, China. This data was primarily obtained through manual surveying, ground-based 3D laser scanning, and UAV lidar. Shangri-La is located in the eastern part of the Hengduan Mountains on the southeastern edge of the Qinghai-Tibet Plateau, with an elevation difference of approximately 4000 m. Data from this region is used to verify the accuracy of the algorithm in calculating tree height in complex terrain areas (Table 1). The specific distribution of the verification sample plots is as follows: Figure 2 As shown.

[0099] Table 1. Relevant surface information of the test site

[0100]

[0101]

[0102] Note: The slope parameters mentioned above are based on the slope range and mean obtained from surface analysis of the 30m SRTM DEM. The vegetation cover is calculated based on the GlobeLand30 (http: / / www.globallandcover.com) primary land use classification product.

[0103] (1) To further illustrate the effects of the present invention, this embodiment presents a comparative analysis of threshold segmentation accuracy, the details of which are as follows:

[0104] To achieve precise denoising with DRAGANN, two core conditions are theoretically required: 1) the local density histogram must have significant segmentation in its shape; and 2) accurate identification of signal peaks and noise peaks must be guaranteed under different signal-to-noise ratios. This invention therefore explores which of these conditions is easily missing during DRAGANN processing in different signal-to-noise ratio scenarios, and proposes corresponding improvement methods for each missing condition. Figure 5 -a is an example of a local density histogram in high signal-to-noise ratio scenarios, characterized by "low and narrow noise peaks and wide and high signal peaks." It already has good signal-to-noise segmentation conditions, but after Savitzky-Golay smoothing, one of its principal components, the noise peak, is filtered out, which in turn causes the denoising threshold to shift towards the signal peak direction. Figure 5 -b). The reason is that the smoothing strength of Savitzky-Golay is often affected by the polynomial order, and the shape characteristics of such histograms make the polynomial order setting sensitive. Considering the difficulty in controlling its smoothing strength, this invention first uses Gaussian low-pass filtering for smoothing, then performs Gaussian decomposition on the smoothing result and optimizes the Gaussian components using the LM algorithm. Finally, under specific conditions, Gaussian components that can represent signal peaks and noise peaks are selected and subjected to bimodal Gaussian refitting, thereby extracting an accurate denoising threshold ( Figure 6 ). Figure 7 -a represents an example of a local density histogram in low signal-to-noise ratio scenarios, characterized by "high and narrow noise peaks and low and wide signal peaks." In this case, the extracted denoising threshold either filters out signal photons or retains some noise photons. Figure 7 -b). This situation is often caused by an improper value of P=20. Therefore, this invention adjusts the value of P using formula (1). Under the same data but different values ​​of P, Figure 8 -a indicates better segmentation conditions. At the same time, it complements... Figure 6 The improved method in the paper applies similar processing to the local density histogram, thus extracting an accurate denoising threshold. Figure 8 -b). In summary, after achieving dynamic acquisition of parameter P and adjusting the signal-to-noise peak identification method, the algorithm of this invention is more accurate and universal in extracting the denoising threshold, which is of great significance for eliminating noise in the scene as much as possible.

[0105] (2) To further illustrate the effects of the present invention, this embodiment presents a comparative analysis of point cloud denoising results, the details of which are as follows:

[0106] Due to solar background noise, the signal-to-noise ratio (SNR) of photon point clouds acquired by ICESat-2 during the day is low, while the SNR of data acquired at night is high. Noise is widely distributed across different locations within the data set, with noise near but below the surface, noise clusters above the canopy, and noise within the forest being the most difficult to remove. Analysis of the ICESat-2 data revealed an issue with signal point filtering in the ALT08 data. Figure 9 -a and Figure 10 -a), where signal filtering is most severe in high signal-to-noise ratio data. Signal filtering is related to the failure to correctly establish the segmentation threshold between noise peaks and signal peaks. The TS-SCABR algorithm effectively solves the point cloud denoising problem in different scenarios, greatly preserving signal data while eliminating noise ( Figure 9 -b and Figure 10 -b); In this invention, the point cloud denoising aspect is an improvement on the official algorithm DNRGANN of the ICESat-2 / ATL08 data product. To facilitate the differentiation of the expression, the improved denoising algorithm is renamed Threshold Segmentation based on Spatial Clustering and Bimodal Reconstruction (TS-SCABR).

[0107] (3) To further illustrate the effects of the present invention, this embodiment presents a comparative analysis of point cloud classification results, the details of which are as follows:

[0108] In scenarios with significant terrain undulations, such as mountainous terrain, the profile point cloud simultaneously exhibits both data trends resulting from changes in photon elevation and terrain trends. Under the influence of terrain trends, the profile point cloud is prone to baseline shift, leading to substantial errors in identifying ground curves using iterative filters. To remove terrain trends from the profile point cloud and focus on analyzing the fluctuation characteristics of data trends, detrending processing is necessary. Furthermore, the terrain trend line fitted during the detrending process can eliminate some unremoved noise clusters in the atmosphere. This type of noise has density characteristics similar to signal photons and is easily identified as signal photons using clustering algorithms. Previous researchers have proposed distance-based outlier detection methods to remove this noise, but these require distance statistics on the signal, making the denoising process complex. This invention uses a terrain trend line and then defines suitable upper and lower bounds as non-noise point intervals for further denoising. The upper limit threshold is typically set at 150m, and the lower limit threshold is set at twice the standard deviation of the Seg-Km processed data.

[0109] After detrending the profile point cloud, an iterative filter was used to effectively identify the upper and lower boundaries of the ground points, and the ground curve was extracted after EMD-DISPO filtering. When a significant abrupt change occurs in the ground points ( Figure 11 The gradient descent method is used to detect the amplitude and location of abrupt changes. Then, a corresponding translation sequence is generated about the horizontal axis (H=0m). The original signal is then subtracted from the translation sequence, and EMD-DISPO filtering is used to eliminate signal distortion, thus completing the de-mutation operation. This de-mutation process eliminates false ground points, greatly improving the accuracy of the ground curve.

[0110] The crown ridge curve fitting method in ICESat-2 / ATBD is similar to that of the ground curve fitting method. However, the crown ridge curve itself has more significant fluctuations than the ground curve. Therefore, even when using iterative moving surface filtering for ground and crown ridge fitting, the ground curve is closer to reality, while the crown ridge curve is often distorted and cannot better fit the actual crown ridge. Figure 12 This situation is particularly prominent in environments with complex land cover types. Therefore, this invention employs dynamic binning statistics to identify crown vertices. Specifically, canopy photons at the [96, 99] percentile are extracted as crown vertices based on binning with a 20m span. Xiaoxiao Zhu et al. (reference: Xiaoxiao Zhu, Sheng Nie, Cheng Wang, Xiaohuan Xi, Zhenyue Hu. (2018) A Ground Elevation and Vegetation Height Retrieval Algorithm Using Micro-Pulse Photon-Counting Lidar Data. Remote Sensing 10:12, pages 1962) also empirically concluded that when statistically analyzing crown vertices at percentiles, daytime and nighttime elevations in the [96, 100] and [99, 100] intervals should be classified as canopy noise, while only the [96, 99] interval can be classified as crown vertices. Furthermore, the detrending processing performed to adapt to the iterative moving surface filtering method for identifying ground curves still has shortcomings. Specifically, areas with significant changes in terrain slope (mountain and valley points) often lead to local widening or overlap of canopy points, thus reshaping the density of local canopy points. Therefore, to reduce the impact of pseudo-canopy vertices on canopy height extraction, neighborhood statistics are also used for further screening of canopy vertices. After a series of algorithm improvements, this invention can achieve accurate point cloud classification under different signal-to-noise ratios. Figure 13 -a and Figure 13 -b) is of great significance for the accurate extraction of the height of the high canopy.

[0111] (4) To further illustrate the effects of the present invention, this embodiment conducted a comparative analysis of the accuracy of canopy height extraction, the details of which are as follows:

[0112] This invention selected tree height data from 150 sampling points to verify the extracted canopy height. For example... Figure 14 -a reflects the linear correlation between the measured data and the ICESat-2 / ATL08 data products, with a root mean square error (RMSE) of 4.018m and a coefficient of determination (R²). 2 It is 0.85212. For example... Figure 14 -b reflects the linear correlation between the measured data and the estimation by the algorithm of this invention, with a root mean square error (RMSE) of 2.013m and a coefficient of determination (R²). 2 The value is 0.96834. Compare the evaluation coefficients RMSE and R... 2 It can be seen that the algorithm in this invention further improves the accuracy of canopy height extraction.

[0113] In summary, based on the official ICESat-2 documentation (ATL08_ATBD), this invention proposes a forest canopy height extraction method with high adaptability, accuracy, and systematicity. After comparing the method with measured data under different terrain conditions and the ATL08 product, the following conclusions can be drawn:

[0114] 1) For clustering denoising algorithms similar to DRAGANN, the universality of their core clustering statistical parameters is particularly important. DRAGANN's original reference value of P=20 only satisfies point cloud denoising conditions in a few high signal-to-noise ratio (SNR) scenarios, but at low SNR, it easily leads to a lack of suitable SNR segmentation conditions in the local density histogram (signal-to-noise peak confusion and crossover). This invention establishes a quantitative relationship model between P and SNR through data experiments at different SNRs, and uses this model to calculate clustering statistical parameters for specific scenarios. Furthermore, while ensuring that the local density histogram of the profile point cloud meets the signal-to-noise segmentation conditions, a systematic mechanism is also needed to accurately locate the position of the signal-to-noise peak. This invention applies RBF filtering and Gaussian decomposition before bimodal Gaussian fitting. By decomposing first and then reconstructing, it effectively avoids the problems that were originally prone to occur: ① In high signal-to-noise ratio scenarios, the noise peak features are not obvious and are filtered out, resulting in a shift in the position of the bimodal fitting; ② In scenarios with complex ground cover types, there are large differences between signal photons, forcing the signal peaks to show a significant "multi-peak" distribution. If the multiple signal peaks are not properly screened, the bimodal Gaussian fitting will tend to fit the relatively obvious peak, resulting in a large error in the denoising threshold at the trough of the fitted bimodal peak.

[0115] 2) Regarding ground curve recognition: ① In scenarios with high vegetation cover, where fewer ground photons return, EMD-DISPO was used to remove pseudo-ground points (canopy photons or noise located near but below the ground); ② In the intersection of forest and non-forest areas, where abrupt changes in surface recognition are prone to occur, a gradient method was used for targeted processing. For canopy recognition, canopy photons at the [96, 99] percentile were statistically analyzed using binning, and the results showed that this recognition method is feasible.

[0116] 3) The extracted forest canopy heights under different terrains were comprehensively verified. Compared with the ATL08 vegetation canopy height product, the estimation accuracy was improved from 85.2% to 96.8%. This shows that the original algorithm still has a lot of room for optimization, while the improved algorithm of this invention has good universality under different undulating terrains and is of great significance for subsequent large-scale forest canopy height inversion.

[0117] The specific embodiments of the present invention have been described in detail above with reference to the accompanying drawings. However, the present invention is not limited to the above embodiments. Within the scope of knowledge possessed by those skilled in the art, various changes can be made without departing from the spirit of the present invention.

Claims

1. A method for extracting forest canopy height based on spaceborne photon counting lidar, characterized in that: First, the signal and noise data of the photon point cloud data are accurately separated. Then, based on the separated signal data, the ground curve is fitted using an iterative filtering algorithm for local terrain. Pseudo-ground points are removed by abrupt change detection using gradient descent and decomposition and reconstruction functions of EMD. Finally, the crown apex is extracted from non-ground points using a dynamic binning statistical method, and the crown height is calculated based on the elevation difference between the crown apex and the ground point. The specific steps of the method are as follows: Step 1: Data preprocessing: Perform coordinate projection transformation and segmentation on the original point cloud to form an elevation profile point cloud with a suitable span. Step 2, Point Cloud Denoising and Classification: Noise, ground, canopy and canopy top photons are identified through point cloud denoising and filtering classification; Step 3, Canopy Height Extraction: The canopy height is calculated by calculating the elevation difference between the canopy apex and the ground point; Step 2, point cloud denoising, includes the following specific steps: Step 2.1.1 Debug filter parameters and perform local density statistics point by point: After data testing, it was found that the filter parameter P is positively correlated with the signal-to-noise ratio (SNR). Then, a relationship model was constructed, formula (1), to achieve adaptive filter parameter P; then, a KD Tree index structure was constructed for all photons in the profile point cloud. After calculating the search radius using formula (2), the photon density was statistically analyzed point by point based on this structure to form a local density histogram. （1）； In the formula, P is the number of neighborhood points, SNR is a rough estimate of the signal-to-noise ratio, and SNR is obtained by statistically analyzing the elevation frequency histogram of the profile point cloud, and then using the median frequency as a threshold. If the frequency is below this threshold, it is considered noise, and otherwise it is considered signal. The mean of each is calculated and divided to obtain SNR, and Int is the floor function. （2）； In the formula, R is the adaptive search radius, P is the number of neighborhood points, and N is the number of neighboring points. Total The Seg_Km function processes the total number of photons within each moving window in the segment, where the span of the moving window is 200m + 3.4km + 200m, and Int is the floor function. Step 2.1.2, Gaussian Decomposition and Parameter Optimization, Bimodal Gaussian Parameter Selection: The histogram curve is filtered using a Gaussian low-pass filter and then subjected to Gaussian decomposition. The Gaussian components are then optimized using the Levenberg-Marquardt LM algorithm. The signal Gaussian and noise Gaussian components are finally selected based on the following criteria: the highest peak on the left is marked as noise Gaussian, and other Gaussian components are selected according to... Sort and the first one satisfies The label is a Gaussian signal; Step 2.1.3 Accurate extraction of noise threshold: After obtaining the noise Gaussian and signal Gaussian parameters, the bimodal features are reconstructed by combining formula (3) and the second derivative is used to extract the unique valley point as the signal-to-noise threshold. （3）； In the formula, The parameters are, in order, the optimized noise Gaussian and signal Gaussian parameters: amplitude, mean, and standard deviation.

2. The method for extracting forest canopy height based on spaceborne photon counting lidar according to claim 1, characterized in that: Step 1 includes: Step 1.1: The original 3D point cloud data based on ICESat-2 / ATL03 is transformed into elevation profile point cloud data through coordinate projection. Step 1.2: Then, the strip profile is segmented according to a 10km span, and each segment is called a Seg-Km processing segment. At the same time, in order to ensure that the extracted canopy height matches the photon point coordinates accurately when compared with the ICESat-2 / ATL08 product, the corresponding ATL08 data is also processed using Seg-Km. In addition, in order to solve the problem of edge effect in the denoising method based on local statistical parameters, the statistical parameter here is the photon density. Usually, the moving window method is combined in the denoising of the Seg-Km processing segment. The span of the box is set to 3.4km, and a 200m buffer is set at the front and back positions of each box along the transverse track.

3. The method for extracting forest canopy height based on spaceborne photon counting lidar according to claim 1, characterized in that: Step 2, point cloud classification includes the following specific steps: Step 2.2.1: Define an iterative filter synthesized by the median filtering algorithm and the moving average algorithm, and since it is based on local terrain conditions, the moving window span of the filter will be set. Step 2.2.2: Detrending processing is performed on the Seg-Km processing segment. First, the original signal after median filtering is denoised under the elevation threshold. Then, the retained signal is subjected to conformal segmentation, cubic Hermite interpolation and iterative filtering to extract the terrain trend line. Finally, the difference between the denoised profile point cloud and the terrain trend line is calculated to achieve detrending. Step 2.2.3: Using all detrended data as the initial value and combining iterative filters to continuously update the fitting point input, the upper and lower bounds of potential ground points are obtained. The ground points within the interval are refined by iterative filters. At the same time, the ground points are interpolated to the original resolution by linear interpolation to obtain the ground curve. EMD-DISPO filters and gradient methods are used to further remove pseudo ground photons in different cases in the ground curve. Step 2.2.4: Photons with percentiles in the range of [96, 99] are identified as top canopy photons using the binning statistical method; the binning span of the profile point cloud is usually set to 20m.

4. The method for extracting forest canopy height based on spaceborne photon counting lidar according to claim 3, characterized in that: In Step 2.2.3, a method combining a DISPO filter and an Empirical Mode Decomposition (EMD) spatiotemporal scale filter is used to remove pseudo-ground points, as shown in formula (4). The EMD algorithm decomposes the data in the time domain, and then uses time-scale denoising or threshold denoising to denoise and reconstruct the decomposition results. In this process, time-scale denoising is prone to causing distortion of the original signal, and threshold denoising is difficult to solve the problem of adaptive threshold calculation. The introduction of the DISPO algorithm first filters the obvious noise in the decomposed signal and then reconstructs it, which can solve the above problems. The original signal and reconstructed signal functions are shown below. （4）； In the formula, H(x) and H1(x) represent the original signal and the reconstructed signal, respectively, and imf i1 / i2 (x) represents a single-component IMF function, r(x) represents the trend term, and N represents the total number of IMF components formed during the EMD decomposition of the ground curve. In the equation H(x), i1 represents the label of each IMF component. When reconstructing the ground curve in EMD, imf1(x) is discarded, and the remaining IMF component indices in [2, N / 2] are filtered by DISPO, while those in [N / 2, N] are not filtered. In the equation H1(x), i2 represents the label of each IMF component that needs to be filtered by DISPO, while j represents the label of each IMF component that does not need to be filtered.

5. The method for extracting forest canopy height based on spaceborne photon counting lidar according to claim 3, characterized in that: In Step 2.2.3, when there are two land types, forest land and non-forest land, in the Seg_Km data processing segment, the density of ground photons under the forest canopy and the density of ground noise photons under the non-forest land are similar. When fitting the ground curve, it is easy to fit the ground noise photons under the non-forest land to the ground, resulting in abrupt changes in terrain. EMD-DISPO is not suitable for processing similar abrupt change signals; the gradient descent method is used to detect the amplitude and location of ground point abrupt change intervals, and the corresponding translation sequence is calculated based on the gradient information and subtracted from the original signal to obtain a more accurate ground curve.

6. The method for extracting forest canopy height based on spaceborne photon counting lidar according to claim 1, characterized in that: Step 3, the extraction of canopy height, includes the following specific steps: Step 3.1: The top of the canopy is slightly lower than the interior of the canopy. Therefore, before calculating the canopy height, the extracted photons from the top of the canopy need to be further filtered. A KD Tree spatial index is constructed for them, and the number of neighborhood points within a 15m local range is counted. Photon point clouds that deviate from the range [3,10] are re-marked as noise, and the rest are the final top photons. Step 3.2: Using the final crown photon extracted in Step 3.1 as the interpolation point, the obtained ground curve is interpolated using the conformal piecewise cubic Hermite interpolation algorithm to obtain the ground photon at the corresponding position in the ground curve. Step 3.3: Finally, calculate the height difference between the canopy photons and the interpolated ground photons to extract the canopy height.

Images