Regional scale under-forest terrain generation method based on multi-mode satellite-borne laser radar

Abstract

The invention discloses a regional scale under-forest terrain generation method based on a multi-mode satellite-borne laser radar, and relates to the technical field of satellite-borne laser radar data application. According to the method, multi-criterion data screening is carried out on satellite-borne laser radar data ICESat-2 and GEDI through characteristic parameters, and a high-precision multi-mode satellite-borne laser radar light spot data set is obtained; secondly, soil types serve as distinguishing conditions, under-forest terrain models for different soil types are constructed through a random forest model, under-forest terrain area scale inversion is achieved from under-forest terrain data of the light spot scale, high-precision under-forest terrain of the area scale is obtained, and the method can be used for conducting under-forest terrain research of the area scale.

Description

Technical Field

[0001] The present invention belongs to the technical field of satellite-borne laser radar data application, and in particular relates to a method for generating regional-scale understory terrain based on multi-modal satellite-borne laser radar. Background Art

[0002] Digital elevation models (DEMs) are widely used in various fields such as topographic mapping, hydrological analysis, forestry research, and ecology. However, in the understory environment, due to the shielding or absorption of optical and radar electromagnetic waves by vegetation such as forests, the understory elevation information detected by traditional remote sensing methods cannot truly reflect the surface conditions. Therefore, how to accurately estimate the understory terrain elevation information has always been a difficult problem in the field of quantitative remote sensing research.

[0003] As an advanced active remote sensing technology, LiDAR (Light Detection and Ranging) has shown great advantages in the field of accurate estimation of forest structural parameters. Unlike airborne LiDAR data, which can only achieve high-precision data measurement of a small coverage area, satellite-borne LiDAR data can provide scientific data for accurate estimation of forest understory terrain in large areas or even on a global scale. In 2018, NASA successfully launched a new generation of satellite-borne LiDAR satellites, ICESat-2 / ATLAS (Ice, Cloud and land Elevation Satellite-2 / Advanced Topographic Laser Altimeter System) and GEDI (Global Ecosystem Dynamics Investigation). The two satellites can provide scientific data for accurate estimation of ground elevation and vegetation canopy height on a global scale. ICESat-2 is a photon counting satellite-borne LiDAR data, and GEDI is a linear satellite-borne LiDAR data. Compared with the previous generation ICESat / GLAS (Ice, Cloud and land Elevation Satellite / Geoscience LaserAltimetry System), ICESat-2 / ATLAS and GEDI have improved to varying degrees in sampling density, spot diameter, resolution, sensitivity and spot interval. The lasers of ICESat-2 / ATLAS and GEDI can emit different types of laser beams, capture laser beams that interact with vegetation, and accurately quantify forest canopy and understory terrain information by analyzing the returned laser beams, thereby providing reliable data for generating high-precision understory terrain. Therefore, ICESat-2 / ATLAS and GEDI are widely used in understory terrain generation.

[0004] However, although ICESat-2 / ATLAS and GEDI can provide high-precision altimetry data, both of them have challenges in data processing to varying degrees. Among them, due to the photon counting system of ICESat-2 / ATLAS, there is a challenge of lack of effective classification of noise photons and signal photons. Although NASA uses algorithms to extract signal photons and effectively divide ground photons and canopy photons in forest environments, and stores them in the ATL08 product, some ATLAS data still have low elevation accuracy due to interference from factors such as atmospheric scattering, solar background radiation and instrument noise, and cannot be directly used as elevation control points for understory terrain. GEDI uses a linear system of recording, which has the problem of difficulty in effectively identifying signal waveforms and forest characteristic waveform data. Although NASA uses algorithms to extract signal waveforms and waveform characteristic information in forest environments for effective division, and stores them in the GEDI L2A data product, some GEDI data still have low elevation accuracy due to interference from factors such as data transmission, atmospheric scattering and instrument noise, and cannot be directly used as elevation control points for understory terrain. At the same time, although collaborative multimodal spaceborne lidar can provide high-precision understory terrain information at the spot scale, it is impossible to estimate the understory terrain at the regional scale using only spaceborne lidar data. Summary of the invention

[0005] In view of the shortcomings of the existing technology, the present invention proposes a regional-scale forest understory terrain generation method based on multi-modal satellite-borne lidar, in order to generate high-precision digital elevation product data of forest understory terrain.

[0006] The regional scale understory terrain generation method based on multi-modal spaceborne lidar includes the following steps:

[0007] Step 1: Determine the forest understory terrain study area and its geographical location information, obtain data products from multimodal satellite-borne lidar data based on the geographical location information, and simultaneously obtain modeling data and classification data of the forest understory terrain study area;

[0008] The multimodal satellite-borne laser radar data is in the form of spots, including satellite-borne laser radar data GEDI and satellite-borne laser radar data ICESat-2 / ATLAS; the data products include GEDIL2A data obtained from satellite-borne laser radar data GEDI, global photon positioning data product ATL03 and land and vegetation height data product ATL08 obtained from satellite-borne laser radar data ICESat-2 / ATLAS;

[0009] Furthermore, the geographical location information is the latitude and longitude information of the forest terrain study area or the boundary vector information of the forest terrain study area; the modeling data includes the digital elevation product, canopy height data and forest coverage data product of the forest terrain study area; the classification data is soil type data;

[0010] Step 2: Extract characteristic parameters from the global photon positioning data product ATL03 and the land and vegetation height data product ATL08 of the forest topography study area, and associate the global photon positioning data product ATL03 with the land and vegetation height data product ATL08 to obtain the associated data product;

[0011] Step 2.1: Extract characteristic parameters from the global photon positioning data product ATL03;

[0012] The characteristic parameters in the global photon positioning data product ATL03 include the longitude lon_ph of the photon, the latitude lat_ph of the photon, the elevation information h_ph of the photon, the 20m segment identifier segment_id, the photon signal confidence signal_conf_ph, the geoid correction information geoid and the photon starting index ph_index_beg;

[0013] Step 2.2: Extract characteristic parameters from land and vegetation height data product ATL08;

[0014] The characteristic parameters in the land and vegetation height data product ATL08 include the longitude of the photon, the latitude of the photon, the best fitting terrain elevation within 20m h_te_best_fit_20m, the terrain elevation uncertainty h_te_uncertainty, the optical centroid height centroid_height, the classified photon cloud flag classed_PC_flag, the classified photon cloud index classed_PC_indx, the 20m segment identifier ph_segment_id, the segment ground cover type segment_landcover, the night flag night_flag, the skewness flag h_te_skew, the multipath scattering warning flag msw_flag, the cloud confidence flag cloud_flag_atm, the signal-to-noise ratio snr, the photon number flag n_seg_ph and the number of photons for the terrain height estimation of a specific segment n_te_photons;

[0015] Step 2.3: By matching the 20m segment identifier segment_id in the global photon positioning data product ATL03 with the 20m segment identifier ph_segment_id in the land and vegetation height data product ATL08, the global photon positioning data product ATL03 is associated with the land and vegetation height data product ATL08 to obtain an associated data product;

[0016] Specifically, firstly, the 20m segment identifier ph_segment_id in the land and vegetation height data product ATL08 is matched with the 20m segment identifier segment_id in the global photon positioning data product ATL03, so as to match the photons in each 20m data segment in the land and vegetation height data product ATL08 and the global photon positioning data product ATL03, and obtain the photon serial numbers in the corresponding 20m data segments. Secondly, the global photon serial number is calculated, and the classified photon cloud index class in ATL08 is ed_PC_indx, i.e., the photon serial number in the 20m data segment, plus the photon starting index ph_index_beg of the corresponding 20m data segment of the global photon positioning data product ATL03 and minus 1, obtains the global photon serial number of the photon in the global photon positioning data product ATL03; based on the global photon serial number, the associated data product is obtained, and the photon corresponding to each global photon serial number in the associated data product has both the characteristic parameters in the land and vegetation height data product ATL08 and the characteristic parameters in the global photon positioning data product ATL03.

[0017] Step 3: Extract the parameter information of waveform spot events and spot quality screening parameters in the GEDI L2A data of the understory terrain study area;

[0018] The parameter information of the waveform light spot event includes: latitude information of the light spot lat_lowestmode_aN, longitude information of the light spot lon_lowestmode_aN, elevation information of the ground under the forest of the light spot elev_lowestmode_aN and geoid correction information of the light spot mean_sea_surface;

[0019] The light spot quality screening parameters include: waveform quality assessment parameter quality_flag, waveform signal-to-noise ratio Sensitivity and satellite status parameter degrade_flag;

[0020] Step 4: Perform geoid correction information correction on the obtained linked data products and GEDI L2A data, so that the linked data products, GEDI L2A data and modeling data are in the same coordinate system, and then resample the linked data products and GEDI L2A data to obtain a multi-modal spaceborne lidar spot data set, which includes the resampled linked data products and the resampled GEDI L2A data;

[0021] Step 4.1: For GEDI L2A data, use the spot's geoid correction information mean_sea_surface to calibrate the spot's understory ground elevation information elev_lowestmode_aN to obtain the elevation correction value H GEDI ;

[0022] H GEDI =elev_lowestmode_aN-mean_sea_surface (1)

[0023] Step 4.2: For linked data products, perform elevation calibration based on the best fitting terrain elevation h_te_best_fit_20m within 20m using the geoid correction information geoid to obtain the elevation correction value H ICESat-2 ;

[0024] H ICESat-2 =h_te_best_fit_20m-geoid (2)

[0025] Step 4.3: Resample the linked data products and GEDI L2A data according to the resolution of the digital elevation products in the modeling data to obtain the resampled linked data products and GEDI L2A data;

[0026] Step 5: Perform multi-criteria data screening on the multi-modal spaceborne lidar spot dataset to obtain the multi-modal spaceborne lidar spot dataset screened under multi-constraint conditions; the multi-modal spaceborne lidar spot dataset screened under multi-constraint conditions includes the ICESat-2 / ATLAS spaceborne lidar spot dataset screened under multi-constraint conditions and the GEDI L2A spaceborne lidar spot dataset screened under multi-constraint conditions;

[0027] Step 5.1: Perform preliminary screening of photons in the resampled correlation data products;

[0028] Specifically: according to the classified photon cloud flag classed_PC_flag in the land and vegetation height data product ATL08, determine whether the photons in the resampled associated data product are ground photons, only retain ground photons with classed_PC_flag=1, and remove other photons; filter according to the photon signal confidence signal_conf_ph in the global photon positioning data product ATL03, and only retain high-confidence photons with signal_conf_ph as the set value;

[0029] Step 5.2: Perform secondary screening on the photons in the initially screened correlation data products;

[0030] Specifically, first, the terrain elevation uncertainty h_te_uncertainty is used for screening, and the screening condition is: retain the photons with h_te_uncertainty≠3.4028235E+38; then the optical centroid height centroid_height is used for screening, and the screening condition is: retain the photons with centroid_height≠3.4028235E+38;

[0031] Step 5.3: Filter the photons in the associated data products filtered in step 5.2 three times to obtain the ICESat-2 / ATLAS space-borne lidar spot dataset filtered under multiple constraints;

[0032] Specifically: use the segment ground cover type segment_landcover to filter, and only retain the forest terrain photons with segment_landcove as the set value; use the night flag night_flag to filter, and retain the night photons with night_flag=1; use the skewness flag h_te_skew to filter, and only retain the photons with |h_te_skew| less than the set threshold; use the multipath scattering warning flag msw_flag to filter, and only retain the photons with msw_flag in the set range; use the cloud confidence flag cloud_flag_atm to filter, and only retain the photons with cloud_flag_atm not greater than the set threshold; use the signal-to-noise ratio snr to filter, and only retain the photons with snr greater than the set threshold; use the photon number flag n_seg_ph to filter, and only retain the photons with n_seg_ph less than the set threshold; use the number of photons estimated by the terrain height of a specific segment n_te_photons to filter, and retain the photons in the 20m data segment where the proportion of terrain photons exceeds the set value;

[0033]

[0034] Where T is the ratio of terrain photons; N terrain is the number of terrain photons, N total is the total number of photons;

[0035] Step 5.4: Filter the photons in the resampled GEDI L2A data according to the spot quality screening parameters; obtain the GEDI L2A space-borne lidar spot data set after screening under multiple constraints;

[0036] Specifically, the photons with quality_flag as the set value are selected according to the waveform quality assessment parameter quality_flag; the photons with sensitivity greater than the set threshold are selected according to the waveform signal-to-noise ratio Sensitivity; the photons with degrade_flag as the set value are selected according to the satellite status parameter degrade_flag;

[0037] Step 6: At the locations with light spots, different soil types are used as differentiation conditions, modeling data is used as independent variables, and the multimodal satellite-borne lidar light spot dataset screened under multiple constraints is used as the dependent variable. The random forest model is used to construct understory terrain models for different soil types.

[0038] Step 6.1: Divide the multimodal satellite-borne lidar spot dataset screened under multiple constraints according to different soil types to obtain the multimodal satellite-borne lidar dataset under different soil types;

[0039] Step 6.2: Use the multimodal satellite-borne lidar datasets under different soil types as dependent variables, digital elevation products, forest canopy height and forest cover as independent variables, and use the random forest model to construct understory terrain models for different soil types;

[0040] Step 7: The coverage area of ​​the multimodal satellite-borne lidar spot data set screened under multiple constraints is used to implement regional scale extrapolation of the understory terrain using modeling data combined with understory terrain models of different soil types to obtain high-precision understory terrain data;

[0041] Step 7.1: The coverage area of ​​the multimodal satellite-borne lidar spot dataset screened under multiple constraints is used to invert the regional-scale forest terrain using forest terrain models of different soil types, combined with digital elevation products, forest canopy height and forest coverage, to generate a new regional-scale forest terrain;

[0042] Step 7.2: Expand the newly generated regional-scale understory terrain in ArcGIS, extract the regional-scale information based on the longitude and latitude information, form a data file, and obtain high-precision understory terrain data.

[0043] Compared with the prior art, the present invention has the following beneficial effects:

[0044] The present invention provides a method for generating regional-scale understory terrain based on multi-modal satellite-borne lidar. Multi-criteria data screening is performed on satellite-borne lidar data ICESat-2 and GEDI through characteristic parameters to obtain a high-precision multi-modal satellite-borne lidar spot data set. Secondly, soil type is used as a distinguishing condition, and understory terrain models for different soil types are constructed respectively through a random forest model. The regional-scale inversion of understory terrain is realized from the spot-scale understory terrain data, and a high-precision regional-scale understory terrain is obtained, which can be used to carry out regional-scale understory terrain research. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 It is a flow chart of a method for generating regional-scale understory terrain based on multi-modal spaceborne laser radar in an embodiment of the present invention;

[0046] Figure 2 This is a schematic diagram of SRTM DEM in an embodiment of the present invention;

[0047] Figure 3 This is a multi-criteria screening flow chart in an embodiment of the present invention;

[0048] Figure 4 is a regional scale forest topography of a certain soil type in an embodiment of the present invention;

[0049] Figure 5 A schematic diagram of the geographical location of the forest topography research area and the light spot distribution of its airborne data and satellite-borne laser radar data according to an embodiment of the present invention;

[0050] Among them, (a) is a schematic diagram of the geographical location of the understory terrain study area; (b) is a schematic diagram of the spot distribution of G-LiHT DTM; (c) is a schematic diagram of the spot distribution of ICESat-2; (d) is a schematic diagram of the spot distribution of GEDI; (e) is a schematic diagram of the spot distribution of the multi-modal satellite-borne lidar spot data set. DETAILED DESCRIPTION

[0051] The present invention is described in detail below with reference to the accompanying drawings and embodiments.

[0052] A regional-scale understory terrain generation method based on multi-modal spaceborne lidar, such as Figure 1 As shown, the following steps are included:

[0053] Step 1: Determine the forest understory terrain study area and its geographical location information, obtain data products from multimodal satellite-borne lidar data based on the geographical location information, and simultaneously obtain modeling data and classification data of the forest understory terrain study area;

[0054] The multi-modal satellite-borne laser radar data is in the form of spots, including satellite-borne laser radar data GEDI and satellite-borne laser radar data ICESat-2 / ATLAS;

[0055] The data products include GEDI L2A data obtained from space-borne lidar data GEDI, global photon positioning data product ATL03 and land and vegetation height data product ATL08 obtained from space-borne lidar data ICESat-2 / ATLAS;

[0056] The geographical location information is the longitude and latitude information of the forest topography research area or the boundary vector information of the forest topography research area;

[0057] The modeling data includes digital elevation products, canopy height data and forest coverage data products of the forest topography research area; the classification data is soil type data;

[0058] In the embodiment of the present invention, the selected understory terrain study area is located in Maryland, and the specific geographical coordinates are 38°-39°43'N and 75°-79°31'W. Maryland has rich vegetation types and a high forest coverage rate, with the highest tree height being 45m. Its altitude ranges from -46m to 1025m. The vegetation types in the understory terrain study area are mainly forests and herbaceous covers, among which the forests are mainly closed evergreen broad-leaved forests and closed deciduous broad-leaved forests. The terrain of the understory terrain study area is flat, with a slope within 10 degrees, which is conducive to ignoring the impact of terrain differences on the uncertainty of the accuracy of satellite-borne lidar.

[0059] Through the official website of NASA (https: / / www.earthdata.nasa.gov / ), the ICESat-2 satellite-borne lidar data and GEDI satellite-borne lidar data of Maryland in similar years were downloaded; the multimodal satellite-borne lidar data GEDI and ICESat-2 / ATLAS were searched and selected according to the longitude and latitude information of the understory terrain study area or the vector information of the boundary of the understory terrain study area, and finally the data products were obtained according to the geographical location, year, data quality and other information of the data products.

[0060] Through the GEE official website, download and obtain the digital elevation product of Maryland. The digital elevation product used in this embodiment is SRTM DEM, which is a global digital elevation product (DEM) obtained by collecting data and performing relevant analysis using a C-band interferometric synthetic aperture radar (InSAR) system in the Shuttle Radar Topography Mission (SRTM) conducted by the National Aeronautics and Space Administration (NASA) and the National Geospatial-Intelligence Agency (NGA). The version used in this embodiment is SRTM V3, which has a resolution of 30m and is based on the 1984 World Geodetic System (WGS84) and the 1996 Earth Gravitational Model (EGM96). This version of SRTM DEM can more accurately extract the terrain information of the United States. Figure 2 Shown is a schematic diagram of the SRTM DEM of Maryland, the understory terrain study area in this embodiment.

[0061] The canopy height data in this embodiment is GEDI V27 data, and the forest cover data product is GFCC30 data;

[0062] Note: The download years of the above data should be close to avoid the influence of time differences on the generated forest topography.

[0063] Step 2: Extract the characteristic parameters of the global photon positioning data product ATL03 and the land and vegetation height data product ATL08 in the forest terrain study area. Since the land and vegetation height data product ATL08 does not record the photon coordinates but only records the classification attributes, the global photon positioning data product ATL03 is associated with the land and vegetation height data product ATL08 to obtain the associated data product. In this way, not only the spatial distribution characteristics of photons can be grasped, but also the detailed classification attributes can be obtained at the same time;

[0064] Step 2.1: Extract characteristic parameters from the global photon positioning data product ATL03;

[0065] The characteristic parameters in the global photon positioning data product ATL03 include the longitude (lon_ph) of the photon, the latitude (lat_ph) of the photon, the elevation information (h_ph) of the photon, the 20m segment identifier (segment_id), the photon signal confidence (signal_conf_ph), the geoid correction information (geoid) and the photon start index (ph_index_beg);

[0066] Step 2.2: Extract characteristic parameters from land and vegetation height data product ATL08;

[0067] The characteristic parameters of the land and vegetation height data product ATL08 include the longitude of the photon, the latitude of the photon, the best fit terrain elevation within 20m (h_te_best_fit_20m), terrain elevation uncertainty (h_te_uncertainty), optical centroid height (centroid_height), classification attributes (classified photon cloud flag (classed_PC_flag), classified photon cloud index (classed_PC_indx), 20m segment identifier (ph_segm ent_id), segment ground cover type (segment_landcover), night flag (night_flag), skew flag (h_te_skew), multipath warning flag (msw_flag), cloud confidence flag (cloud_flag_atm)), signal-to-noise ratio (snr), number of photons flag (n_seg_ph) and number of photons estimated for the terrain height of a specific segment (n_te_photons), where no information is recorded in the longitude of the photons and the latitude of the photons;

[0068] Step 2.3: By matching the 20m segment identifier (segment_id) in the global photon positioning data product ATL03 with the 20m segment identifier (ph_segment_id) in the land and vegetation height data product ATL08, the global photon positioning data product ATL03 is associated with the land and vegetation height data product ATL08 to obtain an associated data product;

[0069] First, the 20m segment identifier (ph_segment_id) in the land and vegetation height data product ATL08 is matched with the 20m segment identifier (segment_id) in the global photon positioning data product ATL03 to match the photons in each 20m data segment in the land and vegetation height data product ATL08 and the global photon positioning data product ATL03, and the photon serial numbers in the corresponding 20m data segments are obtained. Secondly, the global photon serial number is calculated, and the classified photon cloud index (classed_PC_indx) in ATL08, that is, the photon in the 20m data segment, is converted into the global photon serial number. The sub-sequence number is added to the photon starting index (ph_index_beg) of the corresponding 20m data segment of the global photon positioning data product ATL03 and subtracted by 1 (because ph_index_beg starts from 1, and the calculation usually requires an index starting from 0) to obtain the global photon sequence number of the photon in the global photon positioning data product ATL03; according to the global photon sequence number, the associated data product is obtained, and the photon corresponding to each global photon sequence number in the associated data product has both the characteristic parameters in the land and vegetation height data product ATL08 and the characteristic parameters in the global photon positioning data product ATL03.

[0070] Step 3: Extract the parameter information of waveform spot events and spot quality screening parameters in the GEDI L2A data of the understory terrain study area;

[0071] The parameter information of the waveform light spot event includes: latitude information of the light spot (lat_lowestmode_aN), longitude information of the light spot (lon_lowestmode_aN), elevation information of the ground under the forest of the light spot (elev_lowestmode_aN) and geoid correction information of the light spot (mean_sea_surface);

[0072] The light spot quality screening parameters include: waveform quality assessment parameter (quality_flag), waveform signal-to-noise ratio (Sensitivity) and satellite status parameter (degrade_flag);

[0073] Step 4: Perform geoid correction information correction on the obtained linked data products and GEDI L2A data, so that the linked data products, GEDI L2A data and modeling data (digital elevation products, canopy height data, forest cover data) and classification data (soil type data) are all in the WGS84 coordinate system, ensuring that all data sets are at the same scale. Then resample the linked data products and GEDI L2A data to obtain a multimodal spaceborne lidar spot data set, which includes the resampled linked data products and the resampled GEDI L2A data, so that it has the same resolution as other data.

[0074] The specific steps are as follows:

[0075] Step 4.1: For GEDI L2A data, use the spot's geoid correction information (mean_sea_surface) to calibrate the spot's understory ground elevation information (elev_lowestmode_aN) to obtain the elevation correction value H GEDI ;

[0076] H GEDI =elev_lowestmode_aN-mean_sea_surface (4)

[0077] Step 4.2: For linked data products, perform elevation calibration based on the best fit terrain elevation (h_te_best_fit_20m) within 20m using the geoid correction information (geoid) to obtain the elevation correction value H ICESat-2 ;

[0078] H ICESat-2 =h_te_best_fit_20m-geoid (5)

[0079] Among them, the best fitting terrain elevation (h_te_best_fit_20m) and geoid correction information (geoid) within the range of 20m have been synchronously associated with the global photon positioning data product ATL03 and the land and vegetation height data product ATL08 in step 2 according to the ATLAS official data description document and latitude and longitude information. Therefore, the best fitting terrain elevation (h_te_best_fit_20m) in the global photon positioning data product ATL03 and the geoid correction information (geoid) in the land and vegetation height data product ATL08 corresponding to a photon can be obtained at the same time;

[0080] Step 4.3: In order to ensure that the linked data products, GEDI L2A data and modeling data are at the same scale, the linked data products and GEDI L2A data are resampled according to the resolution of the digital elevation products in the modeling data to obtain the resampled linked data products and GEDI L2A data;

[0081] There are many resampling methods. This embodiment takes the bilinear interpolation method as an example. The grid value is calculated by taking the distance from the sampling point to the surrounding four neighboring pixels and interpolating in the X and Y directions respectively. The X direction usually refers to the horizontal direction of the image or grid data, and the Y direction refers to the vertical direction of the image or grid data. First, a linear interpolation is performed in the X direction (or Y direction) to obtain two intermediate values, and then a second linear interpolation is performed in the Y direction (or X direction, depending on the direction of the first interpolation) to finally obtain the value of the point to be interpolated. This interpolation method has an averaged low-pass filtering effect and can produce a relatively coherent output image;

[0082] Assuming α=xi, β=yj, the formula of bilinear interpolation can be expressed as formula (6):

[0083] f(x,y)=(1-α)(1-β)f(i,j)+(1-α)βf(i,j+1)+α(1-β)f(i+1,j)+αβf(i+1,j+1) (6)

[0084] The coordinates of the sampling point are (x, y), and the coordinates of the four nearest neighboring points around it are (i, j), (i+1, j), (i, j+1), and (i+1, j+1), where i and j are integers, x and y can be integers or decimals, and f is the estimated value of the coordinate point.

[0085] Step 5: Figure 3 As shown in the figure, multi-criteria data screening is carried out on the multi-modal satellite-borne lidar spot dataset to reduce the influence of noise and environmental factors on data accuracy, and the multi-modal satellite-borne lidar spot dataset screened under multi-constraint conditions is obtained, thereby ensuring the understory terrain elevation data (elevation correction value H) provided by the screened multi-modal satellite-borne lidar spot dataset. GEDI and elevation correction value H ICESat-2 ) reliability; the multi-modal satellite-borne laser radar spot dataset screened under the multi-constraint conditions includes the ICESat-2 / ATLAS satellite-borne laser radar spot dataset screened under the multi-constraint conditions and the GEDI L2A satellite-borne laser radar spot dataset screened under the multi-constraint conditions;

[0086] Step 5.1: Perform a preliminary screening of the photons in the resampled linked data products to avoid the impact of noise photons on data accuracy and ensure that the photons finally screened can provide real surface information;

[0087] Specifically: according to the classified photon cloud flag (classed_PC_flag) in the land and vegetation height data product ATL08, determine whether the photons in the resampled associated data product are ground photons, and only retain the ground photons with classed_PC_flag=1, and remove other photons; according to the photon signal confidence (signal_conf_ph) in the global photon positioning data product ATL03, only retain the high-confidence photons with signal_conf_ph=4 to avoid the influence of noise photons on data accuracy;

[0088] Step 5.2: Perform secondary screening on the photons in the initially screened associated data products, and try to reduce the low-precision photons caused by factors such as atmosphere, terrain, and sunlight, thereby ensuring the elevation correction value H provided by the remaining photons ICESat-2 Reliability;

[0089] Specifically, first, the terrain elevation uncertainty (h_te_uncertainty) is used for screening, and the screening condition is: retain the photons with h_te_uncertainty≠3.4028235E+38; then the optical centroid height (centroid_height) is used for screening, and the screening condition is: retain the photons with centroid_height≠3.4028235E+38. If it is an invalid value (3.4028235E+38), it cannot reflect the variation range of the elevation measurement results in a specific area and provide meaningful terrain or vegetation height information, so the invalid value photons are deleted;

[0090] Step 5.3: Filter the photons in the associated data products filtered in step 5.2 three times to obtain the ICESat-2 / ATLAS space-borne lidar spot dataset filtered under the most constraints;

[0091] Specifically, the segment ground cover type (segment_landcover) is used for screening, and only the forest terrain photons with segment_landcove=20, 111, 112, 113, 114, 115, 116, 121, 122, 123, 124, 125, and 126 are retained; the night flag (night_flag) is used for screening, and the night photons with night_flag=1 are retained; the skewness flag (h_te_skew) is used for screening, and only the photons with |h_te_skew|<1 are retained to ensure that the elevation values ​​are evenly distributed and reduce the occurrence of outliers; the multipath scattering warning flag (msw_flag) is used for screening, and only the photons with msw_flag=[0,3] are retained to filter out unreliable data, thereby providing Improve the accuracy of data; use the cloud confidence flag (cloud_flag_atm) for screening, only retaining photons with cloud_flag_atm≤1, reducing the impact of attenuation and scattering of photon signals caused by clouds or aerosols on the measurement accuracy of ground elevation; use the signal-to-noise ratio (snr) for screening, only retaining photons with snr>1:3 to ensure the accuracy of elevation measurement; use the photon number flag (n_seg_ph) for screening, only retaining photons with n_seg_ph<500 to avoid introducing errors to ensure the accuracy of elevation measurement; use the number of photons estimated by terrain height in a specific segment (n_te_photons) for screening, retaining photons in the 20m data segment where the proportion of terrain photons exceeds 50%, that is, T>0.5, to optimize the quality of elevation data;

[0092]

[0093] Where T is the ratio of terrain photons; N terrain is the number of terrain photons, N total is the total number of photons;

[0094] Step 5.4: Filter the photons in the resampled GEDI L2A data according to the spot quality screening parameters, and obtain the GEDI L2A space-borne lidar spot data set after screening under multiple constraints;

[0095] Specifically, the photons with quality_flag=1 are screened according to the waveform quality assessment parameter (quality_flag); the photons with sensitivity>0.95 are screened according to the waveform signal-to-noise ratio (Sensitivity); the photons with degrade_flag=0 are screened according to the satellite status parameter (degrade_flag);

[0096] Step 6: At the positions with light spots, using different soil types as the differentiation condition, taking the modeling data as the independent variable, and taking the multi-modal spaceborne lidar light spot dataset screened under multiple constraints as the dependent variable, use the random forest model to construct an understory terrain model for different soil types;

[0097] The specific steps are as follows:

[0098] Step 6.1: Divide the multi-modal spaceborne lidar light spot dataset screened under multiple constraints according to different soil types to obtain multi-modal spaceborne lidar datasets under different soil types;

[0099] Step 6.2: Taking the multi-modal spaceborne lidar datasets under different soil types as the dependent variable, using the digital elevation product, forest canopy height, and forest coverage as the independent variables, use the random forest model to construct an understory terrain model for different soil types;

[0100] The random forest model belongs to a type of Bagging algorithm. Its essence is the idea of ensemble learning in machine learning. The main advantage is that by integrating multiple models, it can effectively handle the overfitting problem and improve the generalization ability of the model. This model is composed of multiple regression decision trees, and its final prediction result is the average of the prediction results of all regression decision trees. The main characteristics of the random forest lie in sample randomness and feature randomness, which are also the reasons why it is superior to a single decision tree. So-called sample randomness means randomly extracting part of the samples from the original dataset by sampling with replacement to construct a new subset. The samples are random and allow a certain degree of repetition, ensuring data diversity, thereby significantly improving the resolution of the feature space and forming a more accurate and smooth decision boundary. Feature randomness means that when choosing the optimal features for splitting, only a randomly selected part of the features are considered, rather than all features. The steps for establishing the random forest model are as follows:

[0101] S1: Use the bootstrap sampling method to randomly extract N samples from the training dataset with replacement to form multiple sub-datasets; the training dataset includes the multi-modal spaceborne lidar dataset, digital elevation product, forest canopy height, and forest coverage under a certain soil type;

[0102] S2: Randomly select m independent variables (m < M, where M is the number of independent variables in the training dataset) from each sub-dataset for splitting to construct a single regression decision tree;

[0103] S3: Repeat the above two steps to generate a large number of regression decision trees, and then obtain the random forest model. The final result of the random forest model is the average of the prediction results of all regression decision trees. The calculation formula is as shown in formula (8):

[0104]

[0105] Where N is the number of decision trees, F i is the predicted value of the i-th decision tree for the input sample;

[0106] In the random forest model, there are two key parameters that affect the running speed of the model and its generalization performance: the number of decision trees and the feature variables randomly selected when each node is split. The more decision trees there are, the better the performance of the model, but the speed of the model will be correspondingly slower; the independent variables randomly selected at each node are calculated based on the minimum mean square error, that is, the data set from any divided variable A to the split point b is D 1 , the data set from the split point b to any partitioned variable S is D 2 , the following calculation formula (9) is used to make D 1 and D 2 The variables and variable division points corresponding to the smallest mean square error of the sets and the smallest sum of the mean square errors of the two sets are:

[0107]

[0108] In the formula, x n is the characteristic variable; n is the true value of the sample; c 1 D 1 The sample output mean of the data set; c 2 D 2 The sample output mean of the data set; the characteristic variables include independent variables and dependent variables.

[0109] Step 7: The coverage area of ​​the multi-modal satellite-borne lidar spot data set screened under multiple constraints is used to achieve regional scale extrapolation of the understory terrain using modeling data combined with understory terrain models of different soil types to obtain high-precision understory terrain data;

[0110] Step 7.1: The multimodal satellite-borne lidar spot dataset covered the area after screening under multiple constraints, and the forest understory models with different soil types were combined with digital elevation products, forest canopy height and forest coverage to realize regional scale forest understory inversion, and generate new regional scale forest understory, such as Figure 4 As shown;

[0111] Step 7.2: Expand the newly generated regional-scale understory terrain in ArcGIS, extract the regional-scale information based on the longitude and latitude information, form a data file, and obtain high-precision understory terrain data to prepare data for the subsequent quantitative evaluation of the understory terrain.

[0112] In order to measure the accuracy of the generated forest terrain, the present invention uses the official free airborne LiDAR data (G-LiHT) product as the accuracy verification data to evaluate the elevation accuracy of the generated high-precision forest terrain data.

[0113] The verification data G-LiHT can be obtained by free downloading from the G-LiHT official website https: / / gliht.gsfc.nasa.gov / . Figure 5 Schematic diagram of the spot distribution of airborne data and satellite-borne lidar data in Maryland, the forest topography study area in an embodiment of the present invention, wherein (a) is a schematic diagram of the geographical location of the forest topography study area; (b) is a schematic diagram of the spot distribution of G-LiHT DTM; (c) is a schematic diagram of the spot distribution of ICESat-2; (d) is a schematic diagram of the spot distribution of GEDI; (e) is a schematic diagram of the spot distribution of a multi-modal satellite-borne lidar spot data set;

[0114] According to the high-precision understory terrain data formed in step 7-2, it is matched with the verification data set (G-LiHT) in the evaluation area through longitude and latitude to form a new data file, and the accuracy verification study of the regional scale understory terrain of the multi-modal satellite-borne lidar spot data set under multi-constraint conditions is carried out;

[0115] In addition, the regional-scale understory terrain inversion achieved by multiple machine learning models can be statistically analyzed, and they can be matched with verification data through information such as longitude and latitude. The efficiency of different machine learning models can be analyzed by calculating accuracy evaluation indicators, thereby ensuring that the accuracy of the regional-scale understory terrain generated by the multi-modal satellite-borne lidar spot data set under multi-constraint conditions reaches a high level.

[0116] The forest terrain elevation accuracy evaluation index in the embodiment of the present invention includes: determination coefficient R 2 , root mean square error RMSE, standard deviation SD and mean square error MSE.

[0117] Among them, the coefficient of determination R 2 The calculation of is as follows (10):

[0118]

[0119] The root mean square error (RMSE) is calculated as follows:

[0120]

[0121] The mean square error (MSE) is calculated as follows:

[0122]

[0123] The standard deviation (SD) is calculated as follows (13):

[0124]

[0125] Among them, y l Extract values ​​for the validation dataset, is the extracted value of the forest terrain at the regional scale, y is the mean extracted value of the validation dataset, and x l It is called the lth training point, and L is the number of instance training data sets in this study.

[0126] Table 1 Comparison of regional scale understory terrain accuracy between SRTM DEM and multi-modal satellite-borne lidar

[0127]

[0128] It can be concluded from Table 1 that the RMSE of the understory terrain accuracy generated by the multi-modal satellite-borne lidar spot dataset under multi-constraint conditions is improved from the original 9.847m to 0.869m, and the accuracy is improved by 91%. 2 It has also been improved from the original 0.996 to 1. Overall, the regional-scale understory terrain generated by the multi-modal satellite-borne lidar spot dataset under multiple constraints has achieved a fairly high accuracy.

Claims

1. A method for generating regional scale understory terrain based on multi-modal spaceborne lidar, characterized in that: The following steps are involved: Step 1: Determine the forest understory terrain study area and its geographical location information, obtain data products from multimodal satellite-borne lidar data based on the geographical location information, and simultaneously obtain modeling data and classification data of the forest understory terrain study area; The multimodal satellite laser radar data is in the form of spots, including satellite laser radar data GEDI and satellite laser radar data ICESat-2 / ATLAS; the data products include GEDI L2A data obtained from satellite laser radar data GEDI, global photon positioning data product ATL03 and land and vegetation height data product ATL08 obtained from satellite laser radar data ICESat-2 / ATLAS; Step 2: Extract characteristic parameters from the global photon positioning data product ATL03 and the land and vegetation height data product ATL08 of the forest topography study area, and associate the global photon positioning data product ATL03 with the land and vegetation height data product ATL08 to obtain the associated data product; Step 3: Extract the parameter information of waveform spot events and spot quality screening parameters in the GEDI L2A data of the understory terrain study area; Step 4: Perform geoid correction information correction on the obtained linked data products and GEDI L2A data, so that the linked data products, GEDI L2A data and modeling data are in the same coordinate system, and then resample the linked data products and GEDIL2A data to obtain a multi-modal spaceborne lidar spot data set, which includes the resampled linked data products and the resampled GEDI L2A data; Step 5: Perform multi-criteria data screening on the multi-modal spaceborne lidar spot dataset to obtain the multi-modal spaceborne lidar spot dataset screened under multi-constraint conditions; the multi-modal spaceborne lidar spot dataset screened under multi-constraint conditions includes the ICESat-2 / ATLAS spaceborne lidar spot dataset screened under multi-constraint conditions and the GEDI L2A spaceborne lidar spot dataset screened under multi-constraint conditions; Step 6: At the locations with light spots, different soil types are used as differentiation conditions, modeling data is used as independent variables, and the multimodal satellite-borne lidar light spot dataset screened under multiple constraints is used as the dependent variable. The random forest model is used to construct understory terrain models for different soil types. Step 7: The coverage area of ​​the multimodal satellite-borne lidar spot dataset screened under multiple constraints is used to achieve regional scale extrapolation of the understory terrain using modeling data combined with understory terrain models of different soil types to obtain high-precision understory terrain data.

2. The method for generating regional-scale understory terrain based on multi-modal spaceborne laser radar according to claim 1, characterized in that: The geographic location information in step 1 is the latitude and longitude information of the forest understory terrain study area or the boundary vector information of the forest understory terrain study area; the modeling data includes digital elevation products, canopy height data and forest coverage data products of the forest understory terrain study area; and the classification data is soil type data.

3. The method for generating regional-scale understory terrain based on multi-modal spaceborne laser radar according to claim 1, characterized in that: The step 2 specifically includes: Step 2.1: Extract characteristic parameters from the global photon positioning data product ATL03; The characteristic parameters in the global photon positioning data product ATL03 include the longitude lon_ph of the photon, the latitude lat_ph of the photon, the elevation information h_ph of the photon, the 20m segment identifier segment_id, the photon signal confidence signal_conf_ph, the geoid correction information geoid and the photon starting index ph_index_beg; Step 2.2: Extract characteristic parameters from land and vegetation height data product ATL08; The characteristic parameters in the land and vegetation height data product ATL08 include the longitude of the photon, the latitude of the photon, the best fitting terrain elevation within 20m h_te_best_fit_20m, the terrain elevation uncertainty h_te_uncertainty, the optical centroid height centroid_height, the classified photon cloud flag classed_PC_flag, the classified photon cloud index classed_PC_indx, the 20m segment identifier ph_segment_id, the segment ground cover type segment_landcover, the night flag night_flag, the skewness flag h_te_skew, the multipath scattering warning flag msw_flag, the cloud confidence flag cloud_flag_atm, the signal-to-noise ratio snr, the photon number flag n_seg_ph and the number of photons for the terrain height estimation of a specific segment n_te_photons; Step 2.3: By matching the 20m segment identifier segment_id in the global photon positioning data product ATL03 with the 20m segment identifier ph_segment_id in the land and vegetation height data product ATL08, the global photon positioning data product ATL03 is associated with the land and vegetation height data product ATL08 to obtain an associated data product; Specifically, firstly, the 20m segment identifier ph_segment_id in the land and vegetation height data product ATL08 is matched with the 20m segment identifier segment_id in the global photon positioning data product ATL03, so as to match the photons in each 20m data segment in the land and vegetation height data product ATL08 and the global photon positioning data product ATL03, and obtain the photon serial numbers in the corresponding 20m data segments. Secondly, the global photon serial number is calculated, and the classified photon cloud index class in ATL08 is ed_PC_indx, i.e., the photon serial number in the 20m data segment, plus the photon starting index ph_index_beg of the corresponding 20m data segment of the global photon positioning data product ATL03 and minus 1, obtains the global photon serial number of the photon in the global photon positioning data product ATL03; based on the global photon serial number, the associated data product is obtained, and the photon corresponding to each global photon serial number in the associated data product has both the characteristic parameters in the land and vegetation height data product ATL08 and the characteristic parameters in the global photon positioning data product ATL03.

4. The method for generating regional-scale understory terrain based on multi-modal satellite-borne laser radar according to claim 3 is characterized in that: The parameter information of the waveform light spot event in step 3 includes: the latitude information lat_lowestmode_aN of the light spot, the longitude information lon_lowestmode_aN of the light spot, the elevation information elev_lowestmode_aN of the ground under the forest of the light spot, and the geoid correction information mean_sea_surface of the light spot; The light spot quality screening parameters include: waveform quality assessment parameter quality_flag, waveform signal-to-noise ratio Sensitivity and satellite status parameter degrade_flag.

5. The method for generating regional-scale understory terrain based on multi-modal spaceborne laser radar according to claim 4, characterized in that: Step 4 specifically includes: Step 4.1: For GEDI L2A data, use the spot's geoid correction information mean_sea_surface to calibrate the spot's understory ground elevation information elev_lowestmode_aN to obtain the elevation correction value H GEDI ; H GEDI =elev_lowestmode_aN-mean_sea_surface (1) Step 4.2: For linked data products, perform elevation calibration based on the best fitting terrain elevation h_te_best_fit_20m within 20m using the geoid correction information geoid to obtain the elevation correction value H ICESat-2 ; H ICESat-2 =h_te_best_fit_20m-geoid (2) Step 4.3: Resample the linked data products and GEDI L2A data according to the resolution of the digital elevation products in the modeling data to obtain the resampled linked data products and GEDI L2A data.

6. The method for generating regional scale understory terrain based on multi-modal spaceborne laser radar according to claim 4, characterized in that: Step 5 specifically includes: Step 5.1: Perform preliminary screening of photons in the resampled correlation data products; Specifically: according to the classified photon cloud flag classed_PC_flag in the land and vegetation height data product ATL08, determine whether the photons in the resampled associated data product are ground photons, only retain ground photons with classed_PC_flag=1, and remove other photons; filter according to the photon signal confidence signal_conf_ph in the global photon positioning data product ATL03, and only retain high-confidence photons with signal_conf_ph as the set value; Step 5.2: Perform secondary screening on the photons in the initially screened correlation data products; Specifically, first, the terrain elevation uncertainty h_te_uncertainty is used for screening, and the screening condition is: retain the photons with h_te_uncertainty≠3.4028235E+38; then the optical centroid height centroid_height is used for screening, and the screening condition is: retain the photons with centroid_height≠3.4028235E+38; Step 5.3: Filter the photons in the associated data products filtered in step 5.2 three times to obtain the ICESat-2 / ATLAS space-borne lidar spot dataset filtered under multiple constraints; Specifically: use the segment ground cover type segment_landcover to filter, and only retain the forest terrain photons with segment_landcove as the set value; use the night flag night_flag to filter, and retain the night photons with night_flag=1; use the skewness flag h_te_skew to filter, and only retain the photons with |h_te_skew| less than the set threshold; use the multipath scattering warning flag msw_flag to filter, and only retain the photons with msw_flag in the set range; use the cloud confidence flag cloud_flag_atm to filter, and only retain the photons with cloud_flag_atm not greater than the set threshold; use the signal-to-noise ratio snr to filter, and only retain the photons with snr greater than the set threshold; use the photon number flag n_seg_ph to filter, and only retain the photons with n_seg_ph less than the set threshold; use the number of photons estimated by the terrain height of a specific segment n_te_photons to filter, and retain the photons in the 20m data segment where the proportion of terrain photons exceeds the set value; Where T is the ratio of terrain photons; N terrain is the number of terrain photons, N total is the total number of photons; Step 5.4: Filter the photons in the resampled GEDI L2A data according to the spot quality screening parameters; obtain the GEDI L2A space-borne lidar spot data set after screening under multiple constraints; Specifically: filter according to the waveform quality assessment parameter quality_flag, and retain the photons with quality_flag as the set value; filter according to the waveform signal-to-noise ratio Sensitivity, and retain the photons with sensitivity greater than the set threshold; filter according to the satellite status parameter degrade_flag, and retain the photons with degrade_flag as the set value.

7. The method for generating regional scale understory terrain based on multi-modal spaceborne laser radar according to claim 1, characterized in that: Step 6 specifically includes: Step 6.1: Divide the multimodal satellite-borne lidar spot dataset screened under multiple constraints according to different soil types to obtain the multimodal satellite-borne lidar dataset under different soil types; Step 6.2: Use the multimodal satellite-borne lidar dataset under different soil types as the dependent variable, digital elevation products, forest canopy height and forest cover as independent variables, and use the random forest model to construct the understory terrain model for different soil types.

8. The method for generating regional-scale understory terrain based on multi-modal spaceborne laser radar according to claim 2, characterized in that: Step 7 specifically includes: Step 7.1: The coverage area of ​​the multimodal satellite-borne lidar spot dataset screened under multiple constraints is used to invert the regional-scale forest terrain using forest terrain models of different soil types, combined with digital elevation products, forest canopy height and forest coverage, to generate a new regional-scale forest terrain; Step 7.2: Expand the newly generated regional-scale understory terrain in ArcGIS, extract the regional-scale information based on the longitude and latitude information, form a data file, and obtain high-precision understory terrain data.

Images