A forest canopy tree height extraction method and system based on laser radar signal fusion

Abstract

The present application provides a kind of forest canopy height extraction method and system of laser radar signal fusion.The present application obtains high-resolution laser radar signal, and the topographic features of high-resolution laser radar signal are obtained by feature extraction;Obtain the topographic feature parameters of the subcomponent corresponding to the high-resolution channel band, combine the least square theory and two optimized paths A and B to carry out iterative cycle iterative calculation residual, obtain the center position of the nth subcomponent;Different waveform corresponding ground elevation is obtained by distance solving formula, atmospheric delay solving, error correction is sequentially obtained, and the corrected laser radar distance length is obtained, the geodetic height corresponding to laser radar footprint is obtained by laser radar footprint solving process, so as to obtain the geodetic height corresponding to forest ground surface.The present application effectively fuses the information of multiple channels, improves the spatial resolution of low channel information, and improves the capture ability of full-waveform laser radar ground surface information.

Description

A method and system for extracting forest canopy tree height by fusion of lidar signals Technical Field

[0001] The technology used in this study belongs to the field of environmental monitoring, and in particular, it relates to a method and system for extracting forest canopy tree height by fusion of lidar signals. Background Technology

[0002] Forests are the largest carbon sink on land, and forest canopy vertical structure information plays a crucial role in carbon cycling, carbon accounting, and the assessment of forest restoration activities. Currently, national forest monitoring inventories and publicly available laser altimeter satellite data are the main sources of canopy vertical structure information. Forest canopy vertical structure information refers to the physical parameters of forests (i.e., tree height, diameter at breast height, and other variables). Forest biomass quantification models are established based on relevant forest canopy vertical structure information. Forest canopy vertical structure information is an integration of data from multiple remote sensing technologies within a region, including ground-based lidar, unmanned aerial vehicles (UAVs), airborne lidar, multi-angle optical satellites, and space-based lidar. Ground-based and airborne lidar can provide high-precision data on forest canopy vertical structure attributes, offering high accuracy and rich spatial detail. However, they are characterized by high cost and small coverage area. UAV photogrammetry offers high temporal resolution, low cost, and flexible observation. However, it requires multiple observations at different times. Multi-angle optical satellites have the advantage of large-area global observation. However, compared to lidar, their accuracy is relatively low and affected by cloud cover. Space-based lidar is characterized by its large range, low cost, and high accuracy. Furthermore, by combining optical remote sensing and lidar data, machine learning methods can be used to map the vertical structure information of the canopy in lidar-unobserved areas. In summary, lidar is an important technology that can accurately capture canopy vertical structure information, featuring a large range, low cost, and high precision.

[0003] Currently, the internationally recognized spaceborne lidar satellites for measuring forest canopy vertical structure information are GEDI and ICESat-2. The GEDI satellite is carried by the International Space Station and is limited by its orbit (51.6°N, 51.6°S). Northern forests are located in the mid-to-high latitudes of the Northern Hemisphere, with a carbon storage as high as (272±23) Pg C, ranking second among all terrestrial ecosystems. Therefore, due to the limitations of its space observation range, the GEDI satellite cannot capture forest canopy vertical structure information in some cold-climate forest regions. GEDI temporarily ceased data collection on March 17, 2023. While the ICESat-2 satellite is designed to quantify ice sheet dynamics, it also provides information on terrestrial vegetation. The ICESat-2 lidar system is a single-photon counting lidar. Compared to full-waveform lidar (i.e., ICESat-1 and GEDI), single-photon lidar has limited measurement accuracy in complex terrain or densely vegetated areas.

[0004] Long-term, continuous forest canopy structure parameters can better reflect the spatial variability and dynamic changes of forest ecosystems. However, the GEDI satellite ceased operation on March 17, 2023. Since the GEDI satellite stopped collecting data, the forest vertical structure information collected by the full-waveform lidar satellite cannot be updated. On April 16, 2022, China launched the world's first atmospheric environment monitoring satellite (DQ-1) equipped with an integral path differential absorption (IPDA) lidar sensor. Its advanced scientific objective is to detect atmospheric aerosol and carbon dioxide concentrations. The DQ-1 satellite is equipped with four-wavelength solid-state lasers at 532.245, 1064.490, 1572.024, and 1572.085 nm, and the data received from the satellite is in full waveform. The spatial vertical accuracy of the 1064.490 nm channel of the DQ-1 satellite is 0.6 meters, and the spatial vertical accuracy of the 1572 nm (i.e., 1572.024 and 1572.085 nm) signal channel is 3 meters, with a footprint diameter of approximately 70 meters. Therefore, the 1064 and 1572 nm channels of the DQ-1 satellite launched after GEDI can invert global ground elevation, thereby further completing the extraction of the vertical structure of the forest canopy.

[0005] The technical problem this invention addresses lies in the multi-channel information fusion of the DQ-1 satellite. This method, combined with the hardware of the DQ-1 satellite, constructs a multi-channel fusion theory (MBFA). By fusing the 1064nm channel, the range resolution of the 1572nm channel is enhanced by approximately five times. The ground elevation and forest canopy products generated by the MBFA theory can be used to assist in generating forest canopy structure parameters for future plans. By combining the forest canopy structure parameters of interest, forest dynamics can be monitored, serving carbon sink and carbon measurement. The products generated by this algorithm ensure the continuity of forest canopy structure monitoring globally and fill the gap in northern forest canopy structure monitoring based on a full-waveform mechanism. Summary of the Invention

[0006] To address the aforementioned technical problems, this invention proposes a method and system for extracting forest canopy tree height through lidar signal fusion. This technical solution can be applied to multi-band lidar satellites to enhance low-resolution signal channel information.

[0007] The technical solution of this invention is a method for extracting forest canopy tree height by fusion of lidar signals, comprising the following steps:

[0008] Step 1: Acquire high-resolution lidar signals. Then, normalize the high-resolution lidar signals, resample the data, match the multi-band signals, and extract the waveform equation features to obtain the terrain features of the high-resolution lidar signals.

[0009] Step 2: Obtain the terrain feature parameters of the corresponding sub-components of the high-resolution channel band, and perform iterative calculation of the residual by combining the least squares theory and two optimized paths A and B. When the residual is less than the specified threshold requirement, the parameter d under different waveforms of the low-resolution band can be obtained, which is regarded as the center position of the nth sub-component.

[0010] Step 3: Obtain the surface elevation corresponding to different waveforms through the distance calculation formula, obtain the corrected lidar distance length through atmospheric delay calculation and error correction, and obtain the geodetic elevation corresponding to the lidar footprint through the lidar footprint calculation process, thereby obtaining the geodetic elevation corresponding to the forest surface.

[0011] Preferably, the data normalization described in step 1 is defined as:

[0012]

[0013] To ensure consistent terrain feature scales in the fused multi-band information, λ a and λ b The echo intensities of different bands are uniformly set to [0,1] according to their respective intensities.

[0014] in, This represents the amplitude value at the nth sampling point of the waveform in band i. It is the minimum value. It is the maximum value. This indicates its normalized amplitude;

[0015] Where i corresponds to the waveform data, that is: and

[0016] in, Indicates λ a On band under the band, Indicates λ a Off-band under the band, Indicates λ b On band under the band, Indicates λ b Off-band under the band;

[0017] The waveform equation described in step 1 is defined as follows:

[0018]

[0019]

[0020]

[0021] in, and Let n represent the nth component of the rising edge duration and the nth component of the falling edge duration of the channel i signal, respectively. This indicates the center position of the nth sub-component in the signal of band i. The constant representing the nth sub-component in band i signal. This represents the scaling factor of the nth sub-component in band i signal. This represents the background noise voltage of the nth sub-component in band i signal. The waveform function representing the nth sub-component in band i signal. This represents the first bias parameter of the nth sub-component in band i signal. This represents the second bias parameter of the nth sub-component in band i signal.

[0022] Based on and Using the waveform information as input and the characteristic parameter sets of different sub-signals in different bands as output, the signal waveform is decomposed using the least squares method to obtain the center position parameter of the nth sub-component in the solved band i signal. Solve for the first bias parameter of the nth sub-component in the i-th band signal. Solve for the second bias parameter of the nth sub-component in the i-th band signal. Solve for the waveform scaling factor of the nth sub-component in the i-th band signal. Solve for the waveform noise voltage of the nth sub-component in the i-th band signal. Solve for the waveform constant of the nth sub-component in the i-th band signal. and This is referred to as the waveform parameter related to terrain features. The solved parameter d is considered as the center position of the nth sub-component and used for subsequent calculation of the lidar column length.

[0023] Preferably, step 2 involves obtaining the terrain feature parameters of the sub-components corresponding to the high-resolution channel band, as follows:

[0024] Combining the method in step 1, the high-resolution channel λ of the DQ-1 satellite is... a Using the on and off bands under the band as input and the terrain feature parameter set as output, the terrain feature parameters of the corresponding sub-components of the high-resolution channel band are obtained.

[0025] Step 2 involves iteratively calculating the residual using least squares theory and two optimized paths A and B, as detailed below:

[0026] Based on the high-resolution channel λ of the DQ-1 satellite a The topographic feature parameter set taken from the on and off bands under the band is used as prior information and the low-resolution channel λ observed by the satellite. b The on-band and off-band information under the band is used as the observation information input, and combined with least squares theory and two optimized paths A and B, iterative loops are performed to process the low-resolution channel λ. b Topographic feature parameter sets of different sub-signals in the on and off bands of the band. and As output, isomorphic parameters related to terrain features are obtained;

[0027] The iterative calculation of residuals described in step 2, wherein the formula for calculating the residuals in each iteration is defined as:

[0028]

[0029]

[0030] in, R represents the observed low-resolution signal. i Represents the residual between i-band observation data and reconstructed data, RF i This represents low-resolution waveform data reconstructed based on terrain feature information extracted from high temporal resolution.

[0031] The optimized path A described in step 2 is defined by the following formula:

[0032]

[0033] HS=J T *J+λ*I,

[0034] Δ p =pinv(HS)*(J T *R),

[0035] p k+1 =p k +Δ p ,

[0036] Where, p k Δ represents the model parameters at the k-th iteration. p p represents the direction and step size of the parameter update. k+1This represents the model parameters after the k-th iteration update, where J represents the Jacobian matrix. λ is the first-order partial derivative of the parameter, m corresponds to the number of parameters, n corresponds to the total amount of data in channel i, λ is the damping coefficient, usually initially set to 0.01, I is the identity matrix, and HS is the Hessian matrix. R represents the observed data.

[0037] The optimized path B described in step 2 is defined by the following formula:

[0038]

[0039] Where N represents the total number of data points in channel i, R i NBL represents the model residual of channel i after the current iteration, and NBL represents the noise baseline being evaluated.

[0040] As a preferred embodiment, step 3 involves obtaining the ground elevation corresponding to different waveforms using the distance calculation formula, as detailed below:

[0041] λ a on-band and λ band under the band a The terrain feature parameter set information of the off channel under the band is used to calculate the surface elevation corresponding to different waveforms by inputting the feature parameter set corresponding to different waveforms into the distance solution formula.

[0042] The distance calculation described in step 3 is defined by the following formula:

[0043] L = cd / 2f,

[0044] Where L represents the spatial distance of the waveform position transition, c represents the speed of light propagation, d represents the data position of the waveform center, and f represents the DQ-1 hardware λ. b Channel acquisition frequency;

[0045] The atmospheric delay solution described in step 3 is defined as follows:

[0046]

[0047] Where N is the refractive index difference, r s It is the starting layer for the surface elevation of the observation station, r a It is the height of the uppermost atmosphere. It is the path integral along the zenith direction.

[0048] The error correction described in step 3 is defined as follows:

[0049]

[0050]

[0051] in, P represents the median error corresponding to the lidar signal. n The weights corresponding to the lidar signals are represented by L; L represents the lidar signal measurement value after processing the above data; and ML represents the corrected lidar distance length.

[0052] The process of solving the lidar footprint described in step 3 is as follows:

[0053]

[0054] Where ρ is the distance value of the lidar signal transmission, α, β, γ represent the deflection angle of the lidar beam relative to the satellite coordinates, and X I Y I and Z I These are the coordinates of a ground point in the satellite platform's reference coordinate system.

[0055]

[0056] Among them, X N Y N and Z N These are the coordinates of a ground point in the navigation coordinate system. This represents the transformation matrix between ground points in the satellite platform reference coordinate system and ground points in the satellite navigation coordinate system.

[0057]

[0058] Among them, H heading P represents the yaw angle. pitch R represents the pitch angle. roll Indicates the roll angle. This represents the transformation matrix between ground points in the satellite platform reference coordinate system and ground points in the satellite navigation coordinate system.

[0059]

[0060] in, Represents the satellite position vector, X gps ,Y gps and Z gps This indicates the GPS position coordinates of the satellite platform in the WGS84 coordinate system.

[0061]

[0062] in, Represents the satellite velocity vector, X v ,Y v and Z v This indicates the speed of the satellite platform in different directions.

[0063]

[0064]

[0065]

[0066] in, Represents the x-axis direction vector. This represents the y-axis direction vector. Represents the z-axis direction vector

[0067]

[0068] in, It represents the transformation matrix between the navigation coordinate system and the geodetic coordinate system for ground points.

[0069]

[0070] Where X, Y, and Z are the coordinates of the lidar footprint on the Earth's surface in a spatial rectangular coordinate system.

[0071]

[0072]

[0073] Where a is the semi-major axis of the reference ellipsoid, b is the semi-minor axis of the reference ellipsoid, e′ is the first eccentricity of the ellipsoid corresponding to the geodetic ellipsoid coordinate system, and c is the flattening of the reference ellipsoid.

[0074]

[0075] Where B is the longitude of the earth.

[0076]

[0077] Where L is the latitude of the Earth.

[0078]

[0079] Where H is the geodetic elevation.

[0080] The technical solution of this invention is a forest canopy tree height extraction system based on lidar signal fusion, comprising:

[0081] The terrain feature extraction module is used to acquire high-resolution lidar signals. It obtains the terrain features of the high-resolution lidar signals by performing data normalization, data resampling, multi-band signal matching, and waveform equation feature extraction.

[0082] The residual iterative calculation module is used to obtain the terrain feature parameters of the corresponding sub-components of the high-resolution channel band. It combines the least squares theory and two optimized paths A and B to perform iterative calculation of the residual. When the residual is less than the specified threshold requirement, the parameter d under different waveforms of the low-resolution band can be obtained after the solution. It is regarded as the center position of the nth sub-component.

[0083] The geodetic elevation calculation module is used to obtain the geodetic elevation corresponding to different waveforms through the distance solution formula. Then, the corrected lidar distance length is obtained through atmospheric delay solution and error correction. Through the lidar footprint solution process, the geodetic elevation corresponding to the lidar footprint is obtained, thereby obtaining the geodetic elevation corresponding to the forest surface.

[0084] The beneficial effects of this invention are as follows: This invention provides a technology for multi-band data fusion and forest canopy height extraction from a full-waveform lidar. This inversion method, combined with data from the high-resolution and low-resolution channels of the DQ-1 satellite, can effectively fuse information from multiple channels, thereby improving the spatial resolution of the low-channel information and enhancing the surface information acquisition capability of the full-waveform lidar. This provides support for the monitoring of global forest canopy height by the Chinese DQ-1 satellite. Attached Figure Description

[0085] Figure 1: Flowchart of the method according to an embodiment of the present invention. Detailed Implementation

[0086] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0087] An embodiment of the method of the present invention is a differential absorption lidar inversion method for carbon dioxide concentration based on a spectral energy model, and a method and system for extracting forest canopy tree height by lidar signal fusion, comprising the following steps:

[0088] Step 1: Acquire high-resolution lidar signals. Then, normalize the high-resolution lidar signals, resample the data, match the multi-band signals, and extract the waveform equation features to obtain the terrain features of the high-resolution lidar signals.

[0089] The data normalization described in step 1 is defined as follows:

[0090]

[0091] To ensure consistent terrain feature scales in the fused multi-band information, λ a and λ b The echo intensities of different bands are uniformly set to [0,1] according to their respective intensities.

[0092] in, This represents the amplitude value at the nth sampling point of the waveform in band i. It is the minimum value. It is the maximum value. This indicates its normalized amplitude;

[0093] Where i corresponds to the waveform data, that is: and

[0094] in, Indicates λ a = 1572nm band on band, Indicates λ a Off-band under the band, Indicates λ a =On band at 1064nm Indicates λ b Off-band under the band;

[0095] The waveform equation described in step 1 is defined as follows:

[0096]

[0097]

[0098]

[0099] in, and Let n represent the nth component of the rising edge duration and the nth component of the falling edge duration of the channel i signal, respectively. This indicates the center position of the nth sub-component in the signal of band i. The constant representing the nth sub-component in band i signal. This represents the scaling factor of the nth sub-component in band i signal. This represents the background noise voltage of the nth sub-component in band i signal. The waveform function representing the nth sub-component in band i signal. This represents the first bias parameter of the nth sub-component in band i signal. This represents the second bias parameter of the nth sub-component in band i signal.

[0100] Based on and Using the waveform information as input and the characteristic parameter sets of different sub-signals in different bands as output, the signal waveform is decomposed using the least squares method to obtain the center position parameter of the nth sub-component in the solved band i signal. Solve for the first bias parameter of the nth sub-component in the i-th band signal. Solve for the second bias parameter of the nth sub-component in the i-th band signal. Solve for the waveform scaling factor of the nth sub-component in the i-th band signal. Solve for the waveform noise voltage of the nth sub-component in the i-th band signal. Solve for the waveform constant of the nth sub-component in the i-th band signal. and This is referred to as the waveform parameter related to terrain features. The solved parameter d is considered as the center position of the nth sub-component and used for subsequent calculation of the lidar column length.

[0101] Step 2: Obtain the terrain feature parameters of the corresponding sub-components of the high-resolution channel band, and perform iterative calculation of the residuals by combining the least squares theory and two optimized paths A and B. The solved parameter d is regarded as the center position of the nth sub-component and used for subsequent calculation of the length of the lidar column.

[0102] Step 2 involves obtaining the terrain feature parameters of the corresponding sub-components of the high-resolution channel band, as detailed below:

[0103] Combining the method in step 1, the high-resolution channel λ of the DQ-1 satellite is... a Using the on and off bands under the band as input and the terrain feature parameter set as output, the terrain feature parameters of the corresponding sub-components of the high-resolution channel band are obtained.

[0104] Step 2 involves iteratively calculating the residual using least squares theory and two optimized paths A and B, as detailed below:

[0105] Based on the high-resolution channel λ of the DQ-1 satellite a The topographic feature parameter set taken from the on and off bands under the band is used as prior information and the low-resolution channel λ observed by the satellite. b The on-band and off-band information under the band is used as the observation information input, and combined with least squares theory and two optimized paths A and B, iterative loops are performed to process the low-resolution channel λ. b Topographic feature parameter sets of different sub-signals in the on and off bands of the band. and As output, isomorphic parameters related to terrain features are obtained;

[0106] The iterative calculation of residuals described in step 2, wherein the formula for calculating the residuals in each iteration is defined as:

[0107]

[0108]

[0109] in, R represents the observed low-resolution signal. i Represents the residual between i-band observation data and reconstructed data, RF i This represents low-resolution waveform data reconstructed based on terrain feature information extracted from high temporal resolution.

[0110] The optimized path A described in step 2 is defined by the following formula:

[0111]

[0112] HS=J T *J+λ*I,

[0113] Δ p =pinv(HS)*(J T *R),

[0114] p k+1 =p k +Δ p ,

[0115] Where, p k Δ represents the model parameters at the k-th iteration. p p represents the direction and step size of the parameter update. k+1 This represents the model parameters after the k-th iteration update, where J represents the Jacobian matrix. λ is the first-order partial derivative of the parameter, m corresponds to the number of parameters, n corresponds to the total amount of data in channel i, λ is the damping coefficient, usually initially set to 0.01, I is the identity matrix, and HS is the Hessian matrix. R represents the observed data.

[0116] The optimized path B described in step 2 is defined by the following formula:

[0117]

[0118] Where N represents the total number of data points in channel i, R i NBL represents the model residual of channel i after the current iteration, and NBL represents the noise baseline being evaluated.

[0119] Step 3: Obtain the surface elevation corresponding to different waveforms through the distance calculation formula, obtain the corrected lidar distance length through atmospheric delay calculation and error correction, and obtain the geodetic elevation corresponding to the lidar footprint through the lidar footprint calculation process, thereby obtaining the geodetic elevation corresponding to the forest surface.

[0120] Step 3 describes obtaining the surface elevation corresponding to different waveforms using the distance calculation formula, as follows:

[0121] λ a on-band and λ band under the band a The terrain feature parameter set information of the off channel under the band is used to calculate the surface elevation corresponding to different waveforms by inputting the feature parameter set corresponding to different waveforms into the distance solution formula.

[0122] The distance calculation described in step 3 is defined by the following formula:

[0123] L = cd / 2f,

[0124] Where L represents the spatial distance of the waveform position transition, c represents the speed of light propagation, d represents the data position of the waveform center, and f represents the DQ-1 hardware λ. b Channel acquisition frequency;

[0125] The atmospheric delay solution described in step 3 is defined as follows:

[0126]

[0127] Where N is the refractive index difference, r s It is the starting layer for the surface elevation of the observation station, r a It is the height of the uppermost atmosphere. It is the path integral along the zenith direction.

[0128] The error correction described in step 3 is defined as follows:

[0129]

[0130]

[0131] in, P represents the median error corresponding to the lidar signal. n The weights corresponding to the lidar signals are represented by L; L represents the lidar signal measurement value after processing the above data; and ML represents the optimized high-precision ranging value.

[0132] The process of solving the lidar footprint described in step 3 is as follows:

[0133]

[0134] Where ρ is the distance value of the lidar signal transmission, α, β, γ represent the deflection angle of the lidar beam relative to the satellite coordinates, and X I Y I and Z I These are the coordinates of a ground point in the satellite platform's reference coordinate system.

[0135]

[0136] Among them, X N Y N and Z N These are the coordinates of a ground point in the navigation coordinate system. This represents the transformation matrix between ground points in the satellite platform reference coordinate system and ground points in the satellite navigation coordinate system.

[0137]

[0138] Among them, H heading P represents the yaw angle. pitch R represents the pitch angle. roll Indicates the roll angle. This represents the transformation matrix between ground points in the satellite platform reference coordinate system and ground points in the satellite navigation coordinate system.

[0139]

[0140] in, Represents the satellite position vector, X gps ,Y gps and Z gpsThis indicates the GPS position coordinates of the satellite platform in the WGS84 coordinate system.

[0141]

[0142] in, Represents the satellite velocity vector, X v ,Y v and Z v This indicates the speed of the satellite platform in different directions.

[0143]

[0144]

[0145]

[0146] in, Represents the x-axis direction vector. This represents the y-axis direction vector. Represents the z-axis direction vector

[0147]

[0148] in, It represents the transformation matrix between the navigation coordinate system and the geodetic coordinate system for ground points.

[0149]

[0150] Where X, Y, and Z are the coordinates of the lidar footprint on the Earth's surface in a spatial rectangular coordinate system.

[0151]

[0152]

[0153] Where a is the semi-major axis of the reference ellipsoid, b is the semi-minor axis of the reference ellipsoid, e′ is the first eccentricity of the ellipsoid corresponding to the geodetic ellipsoid coordinate system, and c is the flattening of the reference ellipsoid.

[0154]

[0155] Where B is the longitude of the earth.

[0156]

[0157] Where L is the latitude of the Earth.

[0158]

[0159] Where H is the geodetic elevation.

[0160] The present invention provides a system for extracting forest canopy tree height through lidar signal fusion, comprising:

[0161] The terrain feature extraction module is used to acquire high-resolution lidar signals. It obtains the terrain features of the high-resolution lidar signals by performing data normalization, data resampling, multi-band signal matching, and waveform equation feature extraction.

[0162] The residual iterative calculation module is used to obtain the terrain feature parameters of the corresponding sub-components of the high-resolution channel band. It combines the least squares theory and two optimized paths A and B to perform iterative calculation of the residual. When the residual is less than the specified threshold requirement, the parameter d under different waveforms of the low-resolution band can be obtained after the solution. It is regarded as the center position of the nth sub-component.

[0163] The geodetic elevation calculation module is used to obtain the geodetic elevation corresponding to different waveforms through the distance solution formula. Then, the corrected lidar distance length is obtained through atmospheric delay solution and error correction. Through the lidar footprint solution process, the geodetic elevation corresponding to the lidar footprint is obtained, thereby obtaining the geodetic elevation corresponding to the forest surface.

[0164] The terrain feature extraction module, residual iterative calculation module, and geodetic elevation calculation module are all deployed on the server.

[0165] It should be understood that any parts not described in detail in this specification belong to the prior art.

[0166] It should be understood that the above description of the embodiments is quite detailed, but it should not be considered as a limitation on the scope of protection of this invention. Those skilled in the art can make substitutions or modifications under the guidance of this invention without departing from the scope of protection of the claims of this invention, and all such substitutions or modifications fall within the scope of protection of this invention. The scope of protection of this invention should be determined by the appended claims.

Claims

1. A method for extracting forest canopy tree height by fusion of lidar signals, characterized in that, Includes the following steps: Step 1: Acquire high-resolution lidar signals. Then, normalize the high-resolution lidar signals, resample the data, match the multi-band signals, and extract the waveform equation features to obtain the terrain features of the high-resolution lidar signals. Step 2: Obtain the terrain feature parameters of the corresponding sub-components of the high-resolution channel band. Combine the least squares theory and two optimized paths A and B to iteratively calculate the residual. When the residual is less than a specified threshold, the parameters of the low-resolution band under different waveforms can be obtained. The nth sub-component is considered as the center position; Step 3: Obtain the surface elevation corresponding to different waveforms through the distance solution formula, obtain the corrected lidar distance length through atmospheric delay solution and error correction, and obtain the geodetic elevation corresponding to the lidar footprint through the lidar footprint solution process, thereby obtaining the geodetic elevation corresponding to the forest surface; wherein, the waveform equation mentioned in Step 1 is defined as: = 1 / 4(1 / 1 / ) = 1 / 4(1 / +1 / )in, and They represent channels respectively. The nth component of the rising edge duration in the signal, channel The nth component of the duration of the falling edge of the signal. Indicates band The center position of the nth sub-component in the signal. Indicates band The constant of the nth sub-component in the signal, Indicates band The scaling factor of the nth sub-component in the signal. Indicates band The background noise voltage of the nth sub-component in the signal. Indicates band The waveform function of the nth sub-component in the signal. Indicates band The first bias parameter of the nth sub-component in the signal. Indicates band The second bias parameter of the nth sub-component in the signal. Indicates band The amplitude value of the waveform at the nth sampling point; based on and Using the waveform information as input and the characteristic parameter sets of different sub-signals in different bands as output, the signal waveform is decomposed using the least squares method to obtain the solved band. The center position parameter of the nth sub-component in the signal Solving for the post-band The first bias parameter of the nth sub-component in the signal Solving for the post-band The second bias parameter of the nth sub-component in the signal Solving for the post-band Waveform scaling factor of the nth sub-component in the signal Solving for the post-band Waveform noise voltage of the nth sub-component in the signal Solving for the post-band Waveform constant of the nth sub-component in the signal ； 、 、 、 、 and These are referred to as waveform parameters related to terrain features; the solved parameters The center position of the nth sub-component is considered as the location of the laser radar column length calculation step 2, which combines least squares theory and two optimized paths A and B to iteratively calculate the residual. Specifically, the residual is calculated based on the high-resolution channel of the DQ-1 satellite. The terrain feature parameter sets acquired from the on and off bands under the band are used as prior information and low-resolution channels observed by satellite. The on-band and off-band information under the band is used as the observation information input, and combined with least squares theory and two optimized paths A and B, iterative loops are performed to process the low-resolution channel. Topographic feature parameter sets of different sub-signals in the on and off bands of the band. 、 、 、 、 and As output, waveform parameters related to terrain features are obtained.

2. The method for extracting forest canopy tree height by fusion of lidar signals according to claim 1, characterized in that: The data normalization mentioned in step 1 is defined as follows: To ensure consistent terrain feature scale after fusion of multi-band information, and The echo intensities of different bands are uniformly set to [0,1] according to their respective intensities; where, Indicates band The amplitude value of the waveform at the nth sampling point. It is the minimum value. It is the maximum value. This represents its normalized amplitude; where Corresponding waveform data, i.e.: 、 、 and ;in, express On band under the band, express Off-band under the band, express On band under the band, express Off-band under the band.

3. The method for extracting forest canopy tree height by fusion of lidar signals according to claim 2, characterized in that: Step 2 involves obtaining the terrain feature parameters of the sub-components corresponding to the high-resolution channel bands, specifically as follows: Combining the method in Step 1, the high-resolution channel of the DQ-1 satellite... Using the on and off bands of the band as input and the terrain feature parameter set as output, the terrain feature parameters of the corresponding sub-components of the high-resolution channel band are obtained.

4. The method for extracting forest canopy tree height by fusion of lidar signals according to claim 3, characterized in that: The optimized path A described in step 2 is defined by the following formula: in, Indicates the first Model parameters at the next iteration Indicates the direction and step size of the parameter update. Indicates the first The model parameters are updated in the next iteration. This represents the Jacobian matrix. It is the first-order partial derivative of the parameter. The corresponding number of parameters correspond Total data volume of the channel This is the damping coefficient, which is usually initially set to 0.

01. It is the identity matrix. It is a Hessian matrix; The observed data; the optimized path B described in step 2, is defined by the formula: ,in, Representative Channel Total number of data points Represents the channel after the current iteration The model residuals This represents the noise baseline being evaluated.

5. The method for extracting forest canopy tree height by fusion of lidar signals according to claim 4, characterized in that: Step 3 describes obtaining the surface elevation corresponding to different waveforms using the distance calculation formula, as follows: On band under the band, The terrain feature parameter set information of the off channel under the band is used to obtain the distance solution formula for the surface elevation corresponding to different waveforms, as described in step 3. The formula is defined as follows: in, Represents the spatial distance of waveform position transition. Represents the speed of light propagation. The data position representing the center of the waveform. This represents the DQ-1 hardware. Channel acquisition frequency.

6. The method for extracting forest canopy tree height by fusion of lidar signals according to claim 5, characterized in that: The atmospheric delay solution described in step 3 is defined as follows: Where N is the refractive index difference, It is the starting layer for the surface elevation of the observation station. It is the height of the uppermost atmosphere. It is the path integral along the zenith direction; the error correction mentioned in step 3 is defined as follows: ， in, This represents the median error corresponding to the lidar signal; This indicates the weight corresponding to the lidar signal; This indicates the corrected lidar range length.

7. The method for extracting forest canopy tree height by fusion of lidar signals according to claim 6, characterized in that: The process of solving the lidar footprint described in step 3 is as follows: in, This is the distance value of the lidar signal transmission. ， ， This represents the angle of deviation of the lidar beam relative to the satellite coordinates. 、 and These are the coordinates of a ground point in the satellite platform's reference coordinate system. in, 、 and These are the coordinates of a ground point in the navigation coordinate system. This represents the transformation matrix between ground points in the satellite platform reference coordinate system and ground points in the satellite navigation coordinate system. in, Indicates the yaw angle. Indicates pitch angle, Indicates the roll angle. This represents the transformation matrix between ground points in the satellite platform reference coordinate system and ground points in the satellite navigation coordinate system. = in, Represents the satellite position vector. , and This indicates the GPS position coordinates of the satellite platform in the WGS84 coordinate system; = in, Represents the satellite velocity vector. , and Indicates the speed of the satellite platform in different directions; in, express Axial direction vector, express Axial direction vector, express Axial direction vector = in, It is a transformation matrix that represents the transformation between the navigation coordinate system and the geodetic coordinate system of the ground points; in, 、 and The coordinates of the lidar footprint on the Earth's surface in a spatial rectangular coordinate system. ,in, It is the semi-major axis of the reference ellipsoid It is the minor semi-axis of the reference ellipsoid. Let be the first eccentricity of the ellipsoid corresponding to the geodetic ellipsoid coordinate system. It is the flattening of the reference ellipsoid; in, It is the longitude of the earth; in, It is the latitude of the earth; in, It is the elevation of the earth.

8. A forest canopy tree height extraction system based on lidar signal fusion, characterized in that, The forest canopy tree height extraction method for performing lidar signal fusion according to any one of claims 1-7 comprises: a terrain feature extraction module for acquiring high-resolution lidar signals, and obtaining the terrain features of the high-resolution lidar signals through data normalization, data resampling, multi-band signal matching, and waveform equation feature extraction; and a residual iterative calculation module for acquiring the terrain feature parameters of the corresponding sub-components of the high-resolution channel bands, and performing iterative calculation of the residuals using least squares theory and two optimized paths A and B. When the residuals are less than a specified threshold, the parameters of different waveforms in the low-resolution bands can be obtained after the solution. It is regarded as the center position of the nth sub-component; the geodetic elevation calculation module is used to obtain the geodetic elevation corresponding to different waveforms through the distance solution formula, and obtain the corrected lidar distance length through atmospheric delay solution and error correction in sequence. Through the lidar footprint solution process, the geodetic elevation corresponding to the lidar footprint is obtained, thereby obtaining the geodetic elevation corresponding to the forest surface.

Images