An optimal pulse width Gaussian decomposition method for spaceborne full-waveform lidar echoes

Abstract

The invention relates to the technical field of full-waveform radar echo data processing, and discloses an optimal pulse width Gaussian decomposition method for satellite-borne full-waveform laser radar echoes. The method adopts an iterative strategy, starts the decomposed waveform from the highest peak, identifies its peak value and peak position, determines the left and right inflection points of the peak, calculates the initial pulse width parameter according to the inflection point, then widens the numerical range of the pulse width parameter, takes RMSE as an optimization criterion, gradually approaches the pulse width parameter of the target waveform, and improves the pulse width parameter accuracy through repeated iterative calculations. After each parameter decomposition of a Gaussian component is completed, the component is removed from the waveform, and then iterative decomposition is continued with the next maximum peak as the target, until the peak value of the remaining Gaussian component peak is lower than zero, and thus the Gaussian decomposition process of the entire waveform is completed. The algorithm can extract waveform parameters more accurately, and can effectively improve the calculation accuracy of subsequent echo ground feature classification, tree height inversion, leaf area index, etc.

Description

Technical Field

[0001] The invention relates to the technical field of full-waveform radar echo data processing, and in particular to an optimal pulse width Gaussian decomposition method for satellite-borne full-waveform laser radar echoes. Background Art

[0002] In recent years, with the rapid development of remote sensing technology, space-borne laser radar (LiDAR) systems have become an important tool in the field of earth observation. It measures the distance of the surface and objects above it by emitting laser pulses to the surface and receiving reflected echoes. This technology can generate high-precision three-dimensional surface models, which is crucial for ecological monitoring, forest management, topographic mapping, urban planning and other fields. In particular, space-borne laser radar systems provide unique capabilities for large-scale earth observation because they can cover a wide geographical area.

[0003] Full-waveform lidar echo data records the entire process of the interaction between the laser pulse and the ground surface and its covering, and contains rich information on the reflection characteristics of the ground objects. However, the complexity of these data and the diversity of the information they contain make it a technical challenge to extract accurate and reliable information from them. Gaussian decomposition method, as a technology for processing full-waveform echo data, aims to improve the accuracy and efficiency of identifying ground features by decomposing complex echo signals into several Gaussian components.

[0004] The existing Gaussian decomposition methods for satellite-borne full-waveform lidar echoes still have the following shortcomings: (1) The decomposition is not universal. Most of the current decomposition methods are mainly aimed at waveforms of specific structures or targets. For satellite-borne lidar echoes with complex structures, it is difficult to effectively decompose the overlapping peaks or hidden peaks therein; (2) The decomposition accuracy needs to be improved, and the decomposed parameters often have obvious deviations; (3) Most of the current decomposition methods require parameter optimization based on algorithms such as LM to obtain better decomposition results, but the decomposition of some waveforms is prone to fall into local minima, resulting in parameter optimization failure. Summary of the invention

[0005] 1. Technical issues to be solved

[0006] In view of the shortcomings of the prior art, the present invention provides an optimal pulse width Gaussian decomposition method for satellite-borne full-waveform lidar echoes, which has the advantages of accurate decomposition and no need for parameter optimization, thereby solving the above-mentioned technical problems.

[0007] (II) Technical solution

[0008] To achieve the above object, the present invention provides the following technical solution: an optimal pulse width Gaussian decomposition method for a satellite-borne full-waveform laser radar echo, comprising the following steps:

[0009] S1. Extract GEDI 01B data, including satellite-borne full-waveform lidar echo data, background noise mean and background noise standard deviation;

[0010] S2. Processing the satellite-borne full-waveform laser radar echo data waveform, specifically including the following steps:

[0011] S2.1, remove background noise, the specific calculation expression is as follows:

[0012] p1(t)=p(t)-ε

[0013] Where p1(t) represents the waveform data after background noise removal, ε represents the background noise, p(t) represents the spaceborne full-waveform lidar echo data, and p(t) is a discrete data point [(t1, y1), (t2, y2),…, (t m ,y m )], (t1, y1), (t2, y2),…, (t m ,y m ) represent the coordinates of 1 to m discrete points, t1, t2, ..., t m Represent 1 to m different times, y1, y2, ..., y m represents 1 to m different amplitudes, [*] represents the discrete data matrix of full-waveform spaceborne lidar echo data;

[0014] S2.2, extracting the signal waveform of the waveform data p1(t) after removing the background noise in step S2.1, specifically comprising:

[0015] S2.2.1. Determine the local peak maximum value y max , with the local maximum peak position As the dividing line, the waveform is divided into two parts, left and right, and the peak and trough points of the waveform are determined by solving the first-order partial derivative of the waveform data;

[0016] S2.2.2. Find the left and right trough points of the local maximum signal waveform, and use the local maximum peak position To start, look for a valley point to the left and record it as the left valley point p1(μ l ), left valley point p1(μ l ) satisfies the trough value less than 0, and the left trough point p1(μ l )The first peak point on the left and the left valley point p1(μ l ) is greater than the background noise standard deviation δ, with the local maximum peak position To start, look for a valley point to the right and record it as the right valley point p1(μ r ), right valley point p1(μ r) satisfies the trough value less than 0, and the right trough point p1(μ r )The first peak point value on the right and the right trough point p1(μ r ) is greater than the standard deviation of the background noise δ;

[0017] S2.2.3, find the left and right nodes: take the left valley point p1(μ l ) The first one on the right that is greater than 0 is taken as the left node p1(tl), and the right valley point p1(μ r )The first node on the left that is greater than 0 is the right node p1(t r );

[0018] S2.2.4. Combine all extracted local signal waveforms to obtain a complete signal waveform p signal (t);

[0019] S3, waveform denoising, use SG filter to remove the complete signal waveform p signal (t) glitch noise;

[0020] S4, waveform compensation, convolving the filtered waveform with the Gaussian sharpening operator to restore the peak and valley values ​​​​distorted during the filtering process;

[0021] S5, waveform decomposition and calculation of optimal pulse width, including the following steps:

[0022] S5.1. Take the waveform after waveform compensation as the waveform to be decomposed, and find the maximum peak from the waveform to be decomposed to obtain the peak value A of the maximum Gaussian component i and peak position μ i ;

[0023] S5.2. Take the second-order derivative of the decomposed waveform and find the peak position μ i The inflection point closest to the left and right ends is recorded as the left inflection point t i,l and right turning point t i,r , and passes through the left turning point t i,l , right turning point t i,r and peak position μ i Calculate the initial pulse width σ of the Gaussian component i0 ;

[0024] S5.3. Calculate the optimal pulse width σ i :Comprising the following steps:

[0025] S5.3.1, with the initial pulse width σ i0 As the initial value, expand the pulse width search range, perform an iteration, and determine the pulse width σ of the nth iteration in scope;

[0026] S5.3.2, with the left turning point t i,land right turning point t i,r The root mean square error RMSE between the Gaussian component and the waveform to be decomposed inf Find the optimal pulse width σ for this round in* , where the optimal pulse width σ in* , specifically, the minimum root mean square error RMSE of the waveform to be decomposed inf The corresponding optimal pulse width σ in* ;

[0027] S5.3.3, and the optimal pulse width σ in* Re-initialize and expand the pulse width search range to repeat the iteration until the minimum RMSE in the current round is reached. inf The value is compared with the minimum RMSE of the previous round inf The decrease in value is less than the minimum RMSE in the previous round inf The iteration ends when the value is 10%, and the estimated optimal pulse width σ is obtained. i* ;

[0028] S5.3.4. Peak A i , peak position μi and the estimated optimal pulse width σ output in step S5.3.3 i* The method includes removing the waveform to be decomposed, and comprising the following steps for preventing the waveform to be decomposed from having a concave waveform below zero:

[0029] S5.3.4.1. Set the number of cycles nc and the pulse width parameter accuracy σ of the cycle iteration p0 , and set the judgment area d to μi-3σ i* to μi+3σ i* ;

[0030] S5.3.4.2. Reduce the estimated optimal pulse width σ i* The value of the pulse width parameter precision σ is reduced by one time each time the loop iteration p0 ;

[0031] S5.3.4.3. Determine the waveform to be decomposed minus the estimated optimal pulse width σ i* Is there a sinking waveform below zero in the range d after the Gaussian component of the value? If not, the optimal pulse width σ is set for this round. i* As the optimal pulse width σ i , if yes, return to step S5.3.4.2;

[0032] S5.4. After completing the waveform decomposition on the waveform to be decomposed, return to step S5.1 and remove the peak value A in step S5.3.4. i , peak position μi and the current optimal pulse width σ output in step S5.3.3 in* The waveform is re-decomposed until the peak value A i<0 or reaches the set number of returns, the decomposition of the entire waveform is completed;

[0033] S6, accuracy assessment, performing error analysis between the waveform to be decomposed and the waveform obtained by waveform decomposition in step S5.

[0034] As a preferred technical solution of the present invention, the specific expression for dividing the waveform into left and right parts in step S2.2.1 is:

[0035] Left part of the waveform data:

[0036]

[0037] The right part of the waveform data:

[0038]

[0039] Among them, t1, ..., tm represent 1 to m different times respectively, y1, ..., ym represent 1 to m different amplitudes, ymax represents the local maximum peak value, and tymax represents the local maximum peak position.

[0040] As a preferred technical solution of the present invention, the calculation expression for determining the peak point and the trough point of the waveform by solving the first-order partial derivative of the waveform data in step S2.2.1 is as follows:

[0041]

[0042] in, A represents the partial derivative calculation of the i-th Gaussian component of the waveform data with respect to t, i represents the peak value, μi represents the peak position, σ represents the pulse width, and e represents the natural logarithm.

[0043] As a preferred technical solution of the present invention, the SG filter in the waveform denoising in step S3 is specifically a fifth-order polynomial SG filter with a window length of 31.

[0044] As a preferred technical solution of the present invention, the specific calculation formula for waveform compensation in step S4 is as follows:

[0045]

[0046] Wherein, g(t) represents the Gaussian sharpening operator, fg(t) represents the Gaussian component of the Gaussian sharpening operator, Ag, μg, σg, Lg are the peak value, peak position, pulse width and half window width of the Gaussian sharpening operator respectively, and λ represents the Gaussian sharpening factor.

[0047] As a preferred technical solution of the present invention, the initial pulse width σ of the Gaussian component is calculated in step S5.2. i0 The expression is as follows:

[0048]

[0049] Among them, ,l Indicates the left turning point, ti ,r represents the right inflection point, μi represents the peak position, σ i0 Indicates the initial pulse width of the Gaussian component.

[0050] As a preferred technical solution of the present invention, the pulse width σ of the nth iteration in step S5.3.1 is in The range expression is as follows:

[0051]

[0052] Among them, σ i0 represents the initial pulse width, σ in represents the pulse width of the nth iteration, n represents the number of iterations, σ p Indicates the accuracy of pulse width, when σ in When it is less than 0, it does not participate in the calculation. i0 -5σ p ,σ i0 -4σ p ,…,σ i0 ,σ i0 +1σ p ,…,σ i0 +5σ p They represent the initial pulse width σ i0 With 5 times σ p The difference between the initial pulse width σ i0 With 4 times σ p The difference between, ..., the initial pulse width σ i0 With 1 times σ p The sum of the values, ..., initial pulse width σ i0 With 5 times σ p The sum between .

[0053] As a preferred technical solution of the present invention, in step S5.3.2, the left turning point ti ,l and right turning point ti ,r The root mean square error RMSE between the Gaussian component and the waveform to be decomposed inf The expression is as follows:

[0054]

[0055] in, Indicates the left turning point ti ,l To the right turning point ,r Sum the data of all points between inf Indicates the left turning point ti,l and right turning point ti ,r The root mean square error between the Gaussian component and the waveform to be decomposed, ti ,l Indicates the left turning point, ti ,r represents the right inflection point, yf(t) represents the waveform to be decomposed, and fi(t) represents the i-th Gaussian component of the waveform data.

[0056] As a preferred technical solution of the present invention, the accuracy evaluation in step S6 includes the following expression:

[0057]

[0058] Wherein, yf(t) represents the waveform to be decomposed, yfit(t) represents the waveform obtained by waveform decomposition in step S5, and represent the mean of the waveform to be decomposed and the waveform obtained by waveform decomposition in step S5, sf 2 and sfit 2 represents the variance of the waveform to be decomposed and the waveform obtained by waveform decomposition in step S5, s represents the covariance of the waveform to be decomposed and the waveform obtained by waveform decomposition in step S5, C1 and C2 represent constants, RMSE represents root mean square error, R 2 represents the coefficient of determination, SSIM represents the structural similarity index, N represents the number of data points, Indicates the sum of N data.

[0059] Compared with the prior art, the present invention provides an optimal pulse width Gaussian decomposition method for spaceborne full-waveform lidar echoes, which has the following beneficial effects:

[0060] 1. The present invention adopts an iterative strategy to start the decomposed waveform from the highest peak, identify its peak value and peak position, determine the left and right inflection points of the peak, and calculate the initial pulse width parameters based on the inflection points, then widen the numerical range of the pulse width parameters, and gradually approach the pulse width parameters of the target waveform with RMSE as the standard, and improve the accuracy of the pulse width parameters by repeated iterative calculations. After each parameter decomposition of a Gaussian component is completed, the component is removed from the waveform, and then iterative decomposition is continued with the next largest peak as the target, until the peak value of the remaining Gaussian component peak is lower than zero, thus completing the Gaussian decomposition process of the entire waveform. The algorithm lays a good foundation for subsequent research on waveform feature parameter extraction, echo feature classification, tree height inversion, leaf area index calculation, etc. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 It is a schematic diagram of the process of the present invention;

[0062] Figure 2It is a schematic diagram of the denoising result of the present invention;

[0063] Figure 3 It is a schematic diagram of the filtering effect and waveform compensation effect of the present invention;

[0064] Figure 4 This is a waveform decomposition diagram of GEDI 01B echo data of the present invention;

[0065] Figure 5 This is a waveform decomposition diagram of GEDI 01B echo data of the present invention;

[0066] Figure 6 It is a schematic diagram of the accuracy evaluation calculation results of the present invention;

[0067] Figure 7 It is a schematic diagram of the process of the present invention. DETAILED DESCRIPTION

[0068] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0069] In this example, Kunming City, Yunnan Province, China (102°10′10″E~103°40′6″E, 24°23′17″N~26°32′42″N), with a data period from December 2018 to June 2022, contains more than 8,000 sets of GEDI 01B echo data of three types of landforms: forests, buildings, and farmland. The data used is obtained from https: / / lpdaac.usgs.gov / tools / earthdata-search /

[0070] See also Figure 1-7 , an optimal pulse width Gaussian decomposition method for spaceborne full-waveform laser radar echoes, comprising the following steps:

[0071] S1. Extract GEDI 01B data, which includes the satellite-borne full-waveform lidar echo data, background noise mean and background noise standard deviation δ;

[0072] S2. Processing the satellite-borne full-waveform laser radar echo data waveform, specifically including the following steps:

[0073] S2.1, remove background noise, the specific calculation expression is as follows:

[0074] p1(t)=p(t)-8

[0075] Wherein, p1(t) represents the waveform data after background noise is removed, ε represents the background noise, p(t) represents the spaceborne full-waveform lidar echo data, and p(t) is a discrete data point [(t1,y1),(t2,y2),…,(tm,ym)], (t1,y1),(t2,y2),…,(tm,ym) represent 1 to m discrete point coordinates respectively, t1, t2,…,tm represent 1 to m different times respectively, y1, y2,…,ym represent 1 to m different amplitudes, [*] represents the discrete data matrix of the full-waveform spaceborne lidar echo data;

[0076] S2.2, extracting the signal waveform of the waveform data p1(t) after removing the background noise in step S2.1, specifically comprising:

[0077] S2.2.1. Determine the local peak maximum value ymax, taking the local maximum peak position As the dividing line, the waveform is divided into two parts, left and right, and the peak and trough points of the waveform are determined by solving the first-order partial derivative of the waveform data;

[0078] S2.2.2. Find the left and right trough points of the local maximum signal waveform, and use the local maximum peak position To start, look for a valley point to the left and record it as the left valley point p1(μ l ), left valley point p1(μ l The valley point value is less than 0, and the left valley point p1(μ l )The first peak point on the left and the left valley point p1(μ l The difference is greater than the background noise standard deviation δ, with the local maximum peak position To start, look for a valley point to the right and record it as the right valley point p1(μ r ), right valley point p1(μ r ) satisfies the trough value less than 0, and the right trough point p1(μ r )The first peak point value on the right and the right trough point p1(μ r ) is greater than the standard deviation of the background noise δ;

[0079] S2.2.3, find the left and right nodes: take the left valley point p1(μ l ) The first one on the right that is greater than 0 is taken as the left node p1(tl), and the right valley point p1(μ r )The first node on the left that is greater than 0 is the right node p1(t r );

[0080] S2.2.4, combining all extracted local signal waveforms to obtain a complete signal waveform psignal(t);

[0081] S3, waveform denoising, using SG filter to remove the burr noise on the complete signal waveform psignal(t);

[0082] S4, waveform compensation, convolving the filtered waveform with the Gaussian sharpening operator to restore the peak and valley values ​​​​distorted during the filtering process;

[0083] S5, waveform decomposition and calculation of optimal pulse width, including the following steps:

[0084] S5.1. Take the waveform after waveform compensation as the waveform to be decomposed, and find the maximum peak from the waveform to be decomposed to obtain the peak value A of the maximum Gaussian component i and peak position μi;

[0085] S5.2. Take the second-order derivative of the decomposed waveform and find the inflection point closest to the left and right ends of the peak μi, which is recorded as the left inflection point ti ,l and right turning point ti ,r and passes through the left turning point ti ,l , right turning point ti ,r and peak position μi to calculate the initial pulse width σ of the Gaussian component i0 ;

[0086] S5.3. Calculate the optimal pulse width σ i :Comprising the following steps:

[0087] S5.3.1, with the initial pulse width σ i0 As the initial value, expand the pulse width search range, perform an iteration, and determine the pulse width σ of the nth iteration in scope;

[0088] S5.3.2, with the left turning point ti ,l and right turning point ti ,r The root mean square error RMSE between the Gaussian component and the waveform to be decomposed inf Find the optimal pulse width σ for this round in* , where the optimal pulse width σ in* , specifically, the minimum root mean square error RMSE of the waveform to be decomposed inf The corresponding optimal pulse width σ in* ;

[0089] S5.3.3, and the optimal pulse width σ in* Re-initialize and expand the pulse width search range to repeat the iteration until the minimum RMSE in the current round is reached. inf The value is compared with the minimum RMSE of the previous round inf The decrease in value is less than the minimum RMSE in the previous round inf The iteration ends when the value is 10%, and the estimated optimal pulse width σ is obtained.i* , set the maximum number of iterations to 4;

[0090] S5.3.4. Peak A i , peak position μi and the estimated optimal pulse width σ output in step S5.3.3 i* The method includes removing the waveform to be decomposed, and comprising the following steps for preventing the waveform to be decomposed from having a concave waveform below zero:

[0091] S5.3.4.1. Set the number of cycles nc and the pulse width parameter accuracy σ of the cycle iteration p0 , and set the judgment area d to μi-3σ i* to μi+3σ i* ;

[0092] S5.3.4.2. Reduce the estimated optimal pulse width σ i* The value of the pulse width parameter is reduced by one time each time the cycle iteration is repeated. p0 ;

[0093] S5.3.4.3. Determine the waveform to be decomposed minus the estimated optimal pulse width σ i* Is there a sinking waveform below zero in the range d after the Gaussian component of the value? If not or the number of cycles reaches 40, the optimal pulse width σ will be estimated. i* As the optimal pulse width σ i , if yes, return to step S5.3.4.2;

[0094] S5.4. After completing the waveform decomposition on the waveform to be decomposed, return to step S5.1 and remove the peak value A in step S5.3.4. i , peak position μi and the current optimal pulse width σ output in step S5.3.3 in* The waveform is re-decomposed until the peak value A i <0 or reaches the set number of returns, the decomposition of the entire waveform is completed;

[0095] S6, accuracy assessment, performing error analysis between the waveform to be decomposed and the waveform obtained by waveform decomposition in step S5.

[0096] Furthermore, the specific expression for dividing the waveform into the left and right parts in step S2.2.1 is:

[0097] Left part of the waveform data:

[0098]

[0099] The right part of the waveform data:

[0100]

[0101] Among them, t1, ..., tm represent 1 to m different times, y1, ..., ym represent 1 to m different amplitudes, ymax represents the local peak maximum value, Indicates the local maximum peak position.

[0102] Furthermore, in step S2.2.1, the calculation expressions for determining the peak points and trough points of the waveform by solving the first-order partial derivative of the waveform data are as follows:

[0103]

[0104] in, A represents the partial derivative calculation of the i-th Gaussian component of the waveform data with respect to t, i represents the peak value, μi represents the peak position, σ represents the pulse width, and e represents the natural logarithm.

[0105] Furthermore, in step S3, the SG filter in the waveform denoising is specifically a fifth-order polynomial SG filter with a window length of 31. The reason for using the fifth-order polynomial SG filter with a window length of 31 is as follows: based on the experience of processing multiple groups of waveforms with different window lengths and orders, the fifth-order polynomial with a window length of 31 has a better processing effect on the GEDI signal waveform, which can effectively remove the burr noise on the waveform without causing obvious signal distortion. Using SG filters of other orders or other lengths to process GEDI echo data will either result in poor denoising effect or cause serious waveform distortion and loss of characteristic information.

[0106] Furthermore, the specific calculation formula for the waveform compensation in step S4 is as follows:

[0107]

[0108] Wherein, g(t) represents the Gaussian sharpening operator, fg(t) represents the Gaussian component of the Gaussian sharpening operator, Ag, μg, σg, and Lg are the peak value, peak position, pulse width, and half-window width of the Gaussian sharpening operator, respectively, λ represents the Gaussian sharpening factor, usually σg is equal to the pulse width of the transmitted pulse, and Lg is taken as 0.2 times the length of the signal waveform.

[0109] Furthermore, in step S5.2, the initial pulse width σ of the Gaussian component is calculated. i0 The expression is as follows:

[0110]

[0111] Among them, ,l Indicates the left turning point, ti ,r represents the right inflection point, μi represents the peak position, σ i0 Indicates the initial pulse width of the Gaussian component.

[0112] Furthermore, the pulse width σ of the nth iteration in step S5.3.1 is in The range expression is as follows:

[0113]

[0114] Among them, σ i0 represents the initial pulse width, σ in represents the pulse width of the nth iteration, n represents the number of iterations, σ p Indicates the accuracy of pulse width, when σ in When it is less than 0, it does not participate in the calculation. i0 -5σ p ,σ i0 -4σ p ,…,σ i0 ,σ i0 +1σ p ,…,σ i0 +5σ p They represent the initial pulse width σ i0 With 5 times σ p The difference between the initial pulse width σ i0 With 4 times σ p The difference between, ..., the initial pulse width σ i0 With 1 times σ p The sum of the values, ..., initial pulse width σ i0 With 5 times σ p If the initial pulse width σ is calculated for the first time, i0 is 4.5ns. In the first iteration, the accuracy of the extended pulse width range is calculated to be 1ns. The extended pulse widths are σ i1 =-0.5,0.5,1.5,…,4.5,…9.5(ns).

[0115] Furthermore, in step S5.3.2, the left turning point t i,l and right turning point t i,r The root mean square error RMSE between the Gaussian component and the waveform to be decomposed inf The expression is as follows:

[0116]

[0117] in, Indicates the left turning point t i,l To the right turning point t i,r Sum the data of all points between inf Indicates the left turning point t i,l and right turning point t i,r The root mean square error between the Gaussian component and the waveform to be decomposed, t i,l Indicates the left turning point, ti,r Indicates the right turning point, y f (t) represents the waveform to be decomposed, f i (t) represents the i-th Gaussian component of the waveform data.

[0118] Furthermore, the accuracy evaluation in step S6 includes the following expression:

[0119]

[0120] Among them, y f (t) represents the waveform to be decomposed, y fit (t) the waveform obtained by waveform decomposition in step S5, and represent the mean of the waveform to be decomposed and the waveform obtained by waveform decomposition in step S5, respectively, f 2 and fit 2 represents the variance of the waveform to be decomposed and the waveform obtained by waveform decomposition in step S5, s represents the covariance of the waveform to be decomposed and the waveform obtained by waveform decomposition in step S5, C1 and C2 represent constants, RMSE represents root mean square error, R 2 represents the coefficient of determination, SSIM represents the structural similarity index, N represents the number of data points, It means to sum N data, C1=6.5025, C2=58.5225, the result of decomposition accuracy index calculation is as follows Figure 6 As shown, Figure 5 It is a schematic diagram of the waveform decomposition of the GEDI 01B echo data of the present invention, wherein Figures a, d, g, and j are schematic diagrams of the decomposition of forest type echoes, Figures b, e, h, and k are schematic diagrams of the decomposition of building type echoes, and Figures c, f, i, and l are schematic diagrams of the decomposition of farmland type echoes.

[0121] Although embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions and variations may be made to the embodiments without departing from the principles and spirit of the present invention, and that the scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. An optimal pulse width Gaussian decomposition method for spaceborne full-waveform laser radar echoes, characterized in that: The following steps are involved: S1. Extract GEDI01B data, which includes satellite-borne full-waveform lidar echo data, background noise mean and background noise standard deviation; S2. Processing the satellite-borne full-waveform laser radar echo data waveform, specifically including the following steps: S2.1, remove background noise, the specific calculation expression is as follows: p1(t)=p(t)-ε Where p1(t) represents the waveform data after background noise removal, ε represents the background noise, p(t) represents the spaceborne full-waveform lidar echo data, and p(t) is a discrete data point [(t1, y1), (t2, y2), …, (t m ,y m )], (t1, y1), (t2, y2),…, (t m ,y m ) represent the coordinates of 1 to m discrete points, t1, t2, ..., t m Represent 1 to m different times, y1, y2, ..., y m represents 1 to m different amplitudes, [*] represents the discrete data matrix of full-waveform spaceborne lidar echo data; S2.2, extracting the signal waveform of the waveform data p1(t) after removing the background noise in step S2.1, specifically comprising: S2.2.

1. Determine the local peak maximum value y max , with the local maximum peak position As the dividing line, the waveform is divided into two parts, left and right, and the peak and trough points of the waveform are determined by solving the first-order partial derivative of the waveform data; S2.2.

2. Find the left and right trough points of the local maximum signal waveform, and use the local maximum peak position To start, look for a valley point to the left and record it as the left valley point p1(μ l ), left valley point p1(μ l ) satisfies the trough value less than 0, and the left trough point p1(μ l )The first peak point on the left and the left valley point p1(μ l ) is greater than the background noise standard deviation δ, with the local maximum peak position To start, look for a trough point to the right and record it as the right trough point p1(μ r ), right valley point p1(μ r ) satisfies the trough value less than 0, and the right trough point p1(μ r )The first peak point value on the right and the right trough point p1(μ r ) is greater than the standard deviation of the background noise δ; S2.2.3, find the left and right nodes: take the left valley point p1(μ l )The first one on the right that is greater than 0 is the left node p1(t l ), with the right trough point p1(μ r )The first node on the left that is greater than 0 is the right node p1(t r ); S2.2.

4. Combine all extracted local signal waveforms to obtain a complete signal waveform p signal (t); S3, waveform denoising, use SG filter to remove the complete signal waveform p signal (t) glitch noise; S4, waveform compensation, convolving the filtered waveform with the Gaussian sharpening operator to restore the peak and valley values ​​​​distorted during the filtering process; S5, waveform decomposition and calculation of optimal pulse width, including the following steps: S5.

1. Take the waveform after waveform compensation as the waveform to be decomposed, and find the maximum peak from the waveform to be decomposed to obtain the peak value A of the maximum Gaussian component i and peak position μ i ; S5.

2. Take the second-order derivative of the decomposed waveform and find the peak position μ i The inflection point closest to the left and right ends is recorded as the left inflection point t i，l and right turning point t i，r , and passes through the left turning point t i，l , right turning point t i，r and peak position μ i Calculate the initial pulse width σ of the Gaussian component i0 ; S5.

3. Calculate the optimal pulse width σ i :Comprising the following steps: S5.3.1, with the initial pulse width σ i0 As the initial value, expand the pulse width search range, perform an iteration, and determine the pulse width σ of the nth iteration in scope; S5.3.2, with the left turning point t i，l and right turning point t i，r The root mean square error RMSE between the Gaussian component and the waveform to be decomposed inf As the optimization criterion, find the optimal pulse width σ in this round in* , where the optimal pulse width σ in* , specifically, the minimum root mean square error RMSE of the waveform to be decomposed inf The corresponding optimal pulse width σ in* ; S5.3.3, and the optimal pulse width σ in* Re-initialize, expand the pulse width search range, and repeat the iteration until the minimum RMSE in the current round is obtained. inf The value is compared with the minimum RMSE of the previous round inf The decrease in value is less than the minimum RMSE in the previous round inf The iteration ends when the value is 10%, and the estimated optimal pulse width σ is obtained. i* ; S5.3.

4. Peak A i , peak position μ i and the estimated optimal pulse width σ output in step S5.3.3 i* The method includes removing the waveform to be decomposed, and comprising the following steps for preventing the waveform to be decomposed from having a concave waveform below zero: S5.3.4.

1. Set the number of cycles n c and the pulse width parameter accuracy of the loop iteration σ p0 , and set the judgment area d to μ i -3σ i* To μ i +3σ i* ; S5.3.4.

2. Reduce the estimated optimal pulse width σ i* The value of the pulse width parameter precision σ is reduced by one time each time the loop iteration p0 ; S5.3.4.

3. Determine the waveform to be decomposed minus the estimated optimal pulse width σ i* Is there a sinking waveform below zero in the range d after the Gaussian component of the value? If not, the optimal pulse width σ is estimated i* As the optimal pulse width σ i , if yes, return to step S5.3.4.2; S5.

4. After completing the waveform decomposition on the waveform to be decomposed, return to step S5.1 and remove the peak value A in step S5.3.

4. i , peak position μ i and the current optimal pulse width σ output in step S5.3.3 in* The waveform is re-decomposed until the peak value A i <0 or reaches the set number of returns, the decomposition of the entire waveform is completed; S6, accuracy assessment, performing error analysis between the waveform to be decomposed and the waveform obtained by waveform decomposition in step S5.

2. The optimal pulse width Gaussian decomposition method of a satellite-borne full-waveform laser radar echo according to claim 1 is characterized in that: The specific expression for dividing the waveform into the left and right parts in step S2.2.1 is: Left part of the waveform data: The right part of the waveform data: Among them, t1, ..., t m Represent 1 to m different times, y1, ..., y m Represents 1 to m different amplitudes, y max represents the local peak maximum, Indicates the local maximum peak position.

3. The optimal pulse width Gaussian decomposition method of a satellite-borne full-waveform laser radar echo according to claim 2 is characterized in that: The calculation expression for determining the peak and trough points of the waveform by solving the first-order partial derivative of the waveform data in step S2.2.1 is as follows: in, A represents the partial derivative calculation of the i-th Gaussian component of the waveform data with respect to t, i Indicates peak value, μ i represents the peak position, σ represents the pulse width, and e represents the natural logarithm.

4. The optimal pulse width Gaussian decomposition method of a satellite-borne full-waveform laser radar echo according to claim 1, characterized in that: The SG filter in the waveform denoising in step S3 is specifically a fifth-order polynomial SG filter with a window length of 31.

5. The optimal pulse width Gaussian decomposition method of a space-borne full-waveform laser radar echo according to claim 4, characterized in that: The specific calculation formula for the waveform compensation in step S4 is as follows: Among them, g(t) represents the Gaussian sharpening operator, f g (t) represents the Gaussian component of the Gaussian sharpening operator, A g , μ g , σ g , L g are the peak value, peak position, pulse width and half window width of the Gaussian sharpening operator respectively, and λ represents the Gaussian sharpening factor.

6. The optimal pulse width Gaussian decomposition method of a satellite-borne full-waveform laser radar echo according to claim 1, characterized in that: The initial pulse width σ of the Gaussian component is calculated in step S5.

2. i0 The expression is as follows: Among them, t i，l Indicates the left turning point, t i，r Indicates the right turning point, μ i represents the peak position, σ i0 Indicates the initial pulse width of the Gaussian component.

7. The optimal pulse width Gaussian decomposition method of a space-borne full-waveform laser radar echo according to claim 6, characterized in that: The pulse width σ of the nth iteration in step S5.3.1 in The range expression is as follows: Among them, σ i0 represents the initial pulse width, σ in represents the pulse width of the nth iteration, n represents the number of iterations, σ p Indicates the accuracy of pulse width, when σ in When it is less than 0, it does not participate in the calculation. i0 -5σ p , σ i0 -4σ p , ..., σ i0 , σ i0 +1σ p , ..., σ i0 +5σ p They represent the initial pulse width σ i0 With 5 times σ p The difference between the initial pulse width σ i0 With 4 times σ p The difference between, ..., the initial pulse width σ i0 With 1 times σ p The sum of the values, ..., initial pulse width σ i0 With 5 times σ p The sum between .

8. The optimal pulse width Gaussian decomposition method of a space-borne full-waveform laser radar echo according to claim 7, characterized in that: In step S5.3.2, the left turning point t i，l and right turning point t i，r The root mean square error RMSE between the Gaussian component and the waveform to be decomposed inf The expression is as follows: in, Indicates the left turning point t i，l To the right turning point t i，r Sum the data of all points between inf Indicates the left turning point t i，l and right turning point t i，r The root mean square error between the Gaussian component and the waveform to be decomposed, t i，l Indicates the left turning point, t i，r Indicates the right turning point, y f (t) represents the waveform to be decomposed, f i (t) represents the i-th Gaussian component of the waveform data.

9. The optimal pulse width Gaussian decomposition method of a space-borne full-waveform laser radar echo according to claim 1, characterized in that: The accuracy evaluation in step S6 includes the following expression: Among them, y f (t) represents the waveform to be decomposed, y fit (t) the waveform obtained by waveform decomposition in step S5, and represent the mean of the waveform to be decomposed and the waveform obtained by waveform decomposition in step S5, respectively, f 2 and fit 2 represents the variance of the waveform to be decomposed and the waveform obtained by waveform decomposition in step S5, s represents the covariance of the waveform to be decomposed and the waveform obtained by waveform decomposition in step S5, C1 and C2 represent constants, RMSE represents root mean square error, R 2 represents the coefficient of determination, SSIM represents the structural similarity index, N represents the number of data points, Indicates the sum of N data.

Images