Full waveform data decomposition method for satellite-borne large-light-spot laser radar

文档序号：1612729 发布日期：2020-01-10 浏览：16次中文

阅读说明：本技术 一种星载大光斑激光雷达全波形数据分解方法 (Full waveform data decomposition method for satellite-borne large-light-spot laser radar ) 是由唐新明刘诏高小明薛玉彩于 2019-11-01 设计创作，主要内容包括：本发明公开了一种星载大光斑激光雷达全波形数据分解方法,该方法包括以下步骤：步骤S1,对星载大光斑激光雷达全波形数据包括的每个原始回波波形数据进行预处理,得到滤波后的回波波形；步骤S2,检测所述滤波后的回波波形中的明显峰值点和与所述明显峰值点对应的拐点,并估计有效峰值点对应的高斯分量的参数的初始值；步骤S3,检测有效明显峰值点两侧的叠加高斯分量,并估计叠加高斯分量的参数初始值；步骤S4,对滤波后的回波波形中检测到的高斯分量的总个数进行更新并对所述高斯分量的参数进行最优化估计。本发明结合拐点匹配和迭代优化算法,改进高斯分量的初始检测策略,得到了较好的回波波形的高斯分量初始估计。(The invention discloses a full waveform data decomposition method of a satellite-borne large-spot laser radar, which comprises the following steps of: step S1, preprocessing each original echo waveform data included in the satellite-borne large-spot laser radar full waveform data to obtain a filtered echo waveform; step S2, detecting an obvious peak point in the filtered echo waveform and an inflection point corresponding to the obvious peak point, and estimating an initial value of a parameter of a Gaussian component corresponding to an effective peak point; step S3, detecting the superposition Gaussian components at both sides of the effective obvious peak point, and estimating the parameter initial value of the superposition Gaussian components; and step S4, updating the total number of the Gaussian components detected in the filtered echo waveform and carrying out optimization estimation on the parameters of the Gaussian components. The initial detection strategy of the Gaussian component is improved by combining the inflection point matching and the iterative optimization algorithm, and the better initial estimation of the Gaussian component of the echo waveform is obtained.)

1. A full waveform data decomposition method for a satellite-borne large-spot laser radar is characterized by comprising the following steps of:

step S1, preprocessing each original echo waveform data included in the satellite-borne large-light-spot laser radar full waveform data to judge whether each original echo waveform is an effective original echo waveform, and performing one-dimensional Gaussian filtering processing on the effective original echo waveform to obtain a filtered echo waveform;

step S2, detecting an obvious peak point and an inflection point corresponding to the obvious peak point in the filtered echo waveform, judging whether the detected obvious peak point is effective or not based on the detected obvious peak point and the inflection point corresponding to the obvious peak point, and estimating an initial value of a parameter of a Gaussian component corresponding to the effective peak point;

step S3, detecting superposition Gaussian components at two sides of the effective obvious peak point, and estimating parameter initial values of the superposition Gaussian components;

and step S4, updating the total number of gaussian components detected in the filtered echo waveform and performing an optimized estimation on the parameters of the gaussian components.

2. The method according to claim 1, wherein step S1 comprises the following sub-steps:

s1.1, estimating the background noise level of each original echo waveform in the satellite-borne large-spot laser radar full waveform data;

step S1.2, setting a background noise threshold value based on the background noise level, comparing the maximum amplitude value of the sampling point of the original echo waveform with the background noise threshold value to judge whether the original echo waveform is an effective original echo waveform, and only reserving the effective original echo waveform;

and S1.3, performing one-dimensional Gaussian filtering processing on the effective original echo waveform to obtain a filtered echo waveform.

And S1.4, calculating the maximum amplitude value of the sampling point of the filtered echo waveform, comparing the maximum amplitude value with a background noise threshold of the echo waveform, if the maximum amplitude value is less than or equal to the background noise threshold, retaining the original effective echo waveform before filtering, and otherwise, retaining the filtered echo waveform.

3. The method according to claim 2, characterized in that in step S1.1, the background noise level of the original echo waveform comprises an average value μ of the background noise of the original echo waveform_nAnd standard deviation of background noise delta_nWherein:

mean value of background noise mu_nThe formula of (c) is shown as follows:

wherein, N is the sum of the sampling points in a selected original echo waveform, N (preferably, N is 20) is the sampling points respectively selected at the beginning part and the ending part of the sampling point of the original echo waveform, and V_iIs the amplitude of the ith sample point.

Standard deviation delta of background noise_nThe formula of (c) is shown as follows:

wherein, N is the sum of the sampling points in a selected original echo waveform, and N (preferably, N is 20) is the sum of the sampling points in the original echo waveformNumber of sampling points, V, respectively selected at the end_iIs the amplitude of the ith sample point.

4. A method according to claim 2, characterized in that in step S1.1, the threshold (th) for background noise is set according to the following formula:

th＝μ_n+m·δ_n

where th is the background noise threshold, μ_nIs the average value of the background noise, δ_nM is a constant value (preferably, m is 4.5) as the standard deviation of the background noise.

If the maximum amplitude value of the sampling point in the original echo waveform is greater than the background noise threshold (th), the original echo waveform can be judged as an effective original echo waveform, otherwise, the original echo waveform is judged as a noise echo waveform, and only the effective original echo waveform is reserved.

5. The method according to claim 2, characterized in that step S1.3 comprises in particular the following sub-steps:

step S1.3.1, determining the width of a one-dimensional Gaussian filter and the window size of the one-dimensional Gaussian filter, wherein the width (sigma) of the one-dimensional Gaussian filter is selected as the full width at half maximum (FWHM) of the emission waveform of the satellite-borne laser radar; the size of the one-dimensional gaussian filter window W is calculated as:

W＝1+2·ceil(3·σ)

ceil (x) is a function that functions to return the smallest integer greater than or equal to a specified expression, W is the window size of the one-dimensional gaussian filter, and σ is the width of the one-dimensional gaussian filter;

step S1.3.2, calculating a gaussian convolution kernel of the one-dimensional gaussian filter, specifically as follows:

C＝floor(W/2)+1

where K (i) is the Gaussian convolution kernel of the one-dimensional Gaussian filter, W is the window size of the one-dimensional Gaussian filter, C is the center position of the Gaussian template, σ is the width of the one-dimensional Gaussian filter, floor (x) is a function whose function is to take the largest integer no greater than x.

The values in the convolution kernel of gaussian are normalized, and the specific calculation is shown as the following formula:

wherein, k (i) is a gaussian convolution kernel of the one-dimensional gaussian filter, W is a window size of the one-dimensional gaussian filter, C is a center position of the gaussian template, and σ is a width of the one-dimensional gaussian filter.

And S1.3.3, performing one-dimensional Gaussian filtering on all sampling points in the effective original echo waveform based on the determined width and window size of the one-dimensional Gaussian filter and the Gaussian convolution kernel.

Specifically, one-dimensional gaussian filtering is performed on all sampling points in the effective original echo waveform according to the following formula:

wherein y (i) is an amplitude value of an ith sampling point of the effective original echo waveform before filtering, y (i) is an amplitude value of an ith sampling point of the filtered echo waveform, L is floor (W/2), W is a window size of the one-dimensional gaussian filter, K (L + j +1) is a gaussian convolution kernel of the one-dimensional gaussian filter, and N is the total number of sampling points in the effective original echo waveform.

6. The method according to claim 1, wherein step S2 comprises the following sub-steps:

s2.1, traversing the filtered echo waveform through a neighborhood window to detect an obvious peak point in the filtered echo waveform;

s2.2, calculating second-order derivatives of the sampling points of the filtered echo waveform, and detecting inflection points on the left side and the right side of the obvious peak point according to the criterion that whether the second-order derivatives of two adjacent sampling points have opposite signs;

s2.3, judging whether the obvious peak point is effective or not based on the detected obvious peak point and an inflection point corresponding to the obvious peak point;

and S2.4, estimating the initial value of the parameter of the Gaussian component corresponding to the determined effective obvious peak value.

7. The method according to claim 6, characterized in that step S2.3 comprises in particular the following sub-steps:

step S2.3.1, if the left inflection point and the right inflection point of the obvious peak point do not exist, judging the obvious peak point as an invalid peak point; if there is only a left inflection point, performing steps S2.3.2-S2.3.4, if there is only a right inflection point, performing steps S2.3.5-S2.3.7, and if there are both a left inflection point and a right inflection point, performing steps S2.3.8-S2.3.10;

step S2.3.2, calculating the average value of the time coordinates of all left inflection points of the obvious peak point;

wherein, mu^LIs the average of the time coordinates of all left inflection points of the distinct peak point, m is the total number of left inflection points of the distinct peak point,the time coordinate value of the k-th inflection point on the left side of the obvious peak point.

Step S2.3.3, calculating the time distance from the left inflection point to the obvious peak point based on the average value of the time coordinates of all the left inflection points and the time coordinate value of the obvious peak point;

d_L＝|μ^L-X(i)|

wherein d is_LThe time distance from the left inflection point to the obvious peak point, X (i) is the obvious peak point Y (i) pairCorresponding time coordinate, mu^LThe average of the time coordinates of all left inflection points of the apparent peak point.

Step S2.3.4, judging whether the time distance from the left inflection point to the obvious peak point is more than or equal to half of the full width at half maximum of the transmitted waveform of the satellite-borne laser radar, if so, determining that the obvious peak point is an effective peak point, otherwise, determining that the obvious peak point is an invalid peak point;

wherein d is_LAnd the FWHM is the full width at half maximum of the transmitted waveform of the satellite-borne laser radar.

Step S2.3.5, calculating the average value of the time coordinates of all right inflection points of the obvious peak point;

wherein, mu^RIs the average of the time coordinates of all right inflection points of the distinct peak point, n is the total number of right inflection points of the distinct peak point,

Step S2.3.6, calculating the time distance from the right inflection point to the obvious peak point based on the average value of the time coordinates of all the right inflection points and the time coordinate value of the obvious peak point;

d_R＝|μ^R-X(i)|

wherein d is_RThe time distance from the right inflection point to the obvious peak point, X (i) is the corresponding time coordinate of the obvious peak point Y (i), mu^RThe average of the time coordinates of all right-hand inflection points of the apparent peak point.

Step S2.3.7, judging whether the time distance from the right inflection point to the obvious peak point is more than or equal to half of the full width at half maximum of the transmitted waveform of the satellite-borne laser radar, if so, determining that the obvious peak point is an effective peak point, otherwise, determining that the obvious peak point is an invalid peak point;

wherein d is_RAnd the FWHM is the full width at half maximum of a transmitting waveform of the satellite-borne laser radar.

Step S2.3.8, calculating the average value of the time coordinates of all left inflection points and all right inflection points of the significant peak point respectively, specifically calculating the following formula:

wherein, mu^LAverage of the time coordinates of all left-hand inflection points of the significant peak point, μ^RIs an average value of time coordinates of all right inflection points of the apparent peak point, m is a total number of left inflection points of the apparent peak point, n is a total number of right inflection points of the apparent peak point,

At step S2.3.9, the time distances from the left inflection point and the right inflection point to the distinct peak point are calculated, respectively. Specifically, the time distances from the left inflection point and the right inflection point to the obvious peak point are calculated according to the following formula:

wherein d is_LIs left inflection point to brightTime distance of peak point, d_RThe time distance from the right inflection point to the obvious peak point, X (i) is the corresponding time coordinate of the obvious peak point Y (i), mu^LAverage of the time coordinates of all left-hand inflection points of the significant peak point, μ^RThe average of the time coordinates of all right-hand inflection points of the apparent peak point.

Step S2.3.10, respectively determining whether the time distance from the left inflection point and the right inflection point to the distinct peak point is greater than or equal to half of the full width at half maximum of the transmitted waveform of the laser radar, if the time distance from at least one of the left inflection point and the right inflection point to the distinct peak point is greater than or equal to half of the full width at half maximum of the transmitted waveform of the laser radar, determining that the distinct peak point is an effective peak point, otherwise, determining that the distinct peak point is an ineffective peak point.

8. The method according to claim 1, wherein step S3 comprises the following sub-steps:

step S3.1, setting an interval [ X (i) -D, X (i) + D ], wherein X (i) is a time coordinate value corresponding to the Gaussian component corresponding to the effective obvious peak point, and D is an estimated initial value of a standard deviation of the Gaussian component corresponding to the effective obvious peak point;

step S3.2, based on the set interval, detecting the superimposed Gaussian component on the left side of the effective obvious peak point according to the following formula:

Wherein j1 and k1 are respectively corresponding positions of two adjacent inflection points detected outside the interval set in the monotonic interval on the left side of the effective significant peak point in X time coordinate data, X (j1) and X (k1) are respectively time coordinate values of two adjacent inflection points detected outside the interval set in the monotonic interval on the left side of the effective significant peak point, y (j1) and y (k1) are respectively amplitude values of the two adjacent inflection points detected in the monotonic area on the left side of the effective significant peak point in an effective original echo waveform before filtering, and FWHM represents the full width at half maximum of a transmitting waveform of the satellite-borne laser radar; if the condition of the above formula is met, detecting a superimposed Gaussian component in two adjacent inflection point regions with corresponding sequence numbers j1 and k1 in the X time coordinate data on the left side of the effective obvious peak point;

step S3.3, based on the set interval, detecting the superimposed Gaussian component on the right side of the effective obvious peak point according to the following formula:

Wherein j2 and k2 are respectively the corresponding positions of two adjacent inflection points detected outside the set interval of the monotonic interval on the right side of the effective significant peak point in X time coordinate data, X (j2) and X (k2) are respectively the time coordinate values of two adjacent inflection points detected outside the set interval of the monotonic interval on the right side of the effective significant peak point, y (j2) and y (j2) are respectively the amplitude values of the two adjacent inflection points detected in the monotonic area on the right side of the effective significant peak point in the effective original echo waveform before filtering, and FWHM represents the full width at half maximum of the transmitted waveform; if the condition of the above equation is satisfied, a superimposed gaussian component is detected in two adjacent inflection regions with corresponding sequence numbers j2, k2 in the X-time coordinate data to the right of the significant peak point.

S3.4, estimating an initial value of the parameter of the superposition Gaussian component; step S3.4 specifically includes the following substeps:

step S3.4.1, estimating the parameter initial value of the Gaussian component superimposed on the left side of the effective obvious peak point; estimating the amplitude A of the left-side superimposed Gaussian component as the amplitude value of the inflection point which is positioned on the right side and is closer to the effective obvious peak point in the first two adjacent inflection points positioned outside the set interval on the left side of the effective obvious peak point; estimating a value of the central time coordinate μ of the superimposed gaussian component to a value of X (k 1); estimating the standard deviation sigma of the superimposed Gaussian component to be half of the distance between two first adjacent inflection points which are positioned on the left side of the effective obvious peak point and are positioned outside the set interval;

step S3.4.2, estimating the parameter initial value of the superimposed Gaussian component on the right side of the effective obvious peak point; estimating the amplitude A of the right-side superimposed Gaussian component as the amplitude value of the inflection point which is positioned on the left side and is closer to the effective obvious peak point in the first two adjacent inflection points positioned outside the set interval on the right side of the obvious peak point; estimating the value of the central time coordinate mu of the superimposed Gaussian component as the value of X (j 2); estimating the standard deviation sigma of the superimposed Gaussian component to be half of the distance between two first adjacent inflection points which are positioned at the right side of the effective obvious peak point and are outside the set interval;

and S3.5, calculating the number of the Gaussian components detected in the filtered echo waveform.

9. The method according to claim 1, wherein step S3 comprises the following sub-steps:

and S4.1, based on an LM algorithm, carrying out iterative optimization on initial values of parameters of Gaussian components (including the Gaussian component corresponding to the effective peak point and the superposed Gaussian component) detected in the filtered echo waveform. Wherein, the step S4.1 specifically comprises the following substeps:

step S4.1.1, an error function and a target function between the fitting function and the observation value selected in the effective original echo waveform (namely the actual observation waveform) are established.

And S4.1.2, iteratively optimizing parameters of all gaussian components (including the gaussian component corresponding to the effective peak point and the superimposed gaussian component) detected in the filtered echo waveform based on the objective function and the error function.

Step S4.2, combining threshold th of fitting residual error RMSE_RMSEAnd further performing iterative optimization on parameter values of Gaussian components contained in the filtered echo waveform. Wherein, the step S4.2 specifically comprises the following substeps:

step S4.2.1 calculates a residual RMSE between the amplitude values of the sampling points in the filtered echo waveform and the fitting values of the amplitude values of the sampling points obtained after the LM algorithm is iteratively optimized.

And step S4.2.2, finding the Gaussian components of the missed ground object target by fitting the threshold of the residual RMSE, updating the number of the Gaussian components and carrying out initial value estimation on the found missed Gaussian components.

Step S4.2.3, further iterative optimization of parameters of the detected Gaussian components in the filtered echo waveform until the fitting residual RMSE is less than or equal to the fitting residual RMSE threshold th_RMSEAnd when so, stopping iteration.

10. A computer program, characterized in that the method of claims 1-9 is performed by executing the computer program.

Technical Field

The invention relates to the technical field of satellite-borne laser height measurement data processing, in particular to a full waveform data decomposition method of a satellite-borne large-spot laser radar.

Background

The full-waveform laser radar technology can record complete waveform information of ground object echo pulses, and physical attributes of the ground objects including laser foot points can be obtained by analyzing backscattering waveforms to reflect vertical distribution characteristics of the ground objects. With the development of the laser radar technology, the satellite-borne full-waveform laser radar technology becomes an effective active earth observation means. In 1999, NASA designed a set of vehicle-mounted large-spot laser altimeter system LVIS (laser targeting Imaging sensor), which can digitally sample laser emission pulses and echo pulses to acquire topographic data of sub-tree crowns and top of canopy. NASA launched the first Ice, Cloud, and land Elevation measurement Satellite ICESat (Ice, Cloud, and land Elevation Satellite) used primarily for polar Ice measurements in 2003, which carried the earth Laser altimetry system glas (geoscience Laser Altimeter system). The GLAS can record full waveform information of return pulses of detected ground objects and can be used for measuring the elevation of the ground surface, extracting the height of the ground objects and the like.

For the full-waveform lidar technology, the processing of full-waveform data is a key step, and the processing methods of full-waveform data mainly include three methods, namely a threshold method, a deconvolution method and a waveform decomposition method. Compared with the former two methods, a waveform decomposition method based on gaussian decomposition is often adopted to process the received waveform.

Disclosure of Invention

Aiming at the problems, the invention provides a method for decomposing full waveform data of a satellite-borne large-spot laser radar, which solves the problem of decomposing the full waveform data of the satellite-borne laser radar under different ground surface coverage types, particularly the problems of detecting and decomposing superposed waveforms, complex multi-peak waveforms and weak echo waveforms in received waveform data, so as to improve the extraction precision of the waveform characteristic parameters of ground objects and further improve the extraction precision of satellite laser elevation control points and the inversion precision of the ground object parameters. Accordingly, it is seen that the present invention is directed to a solution that substantially obviates one or more of the problems due to limitations and disadvantages of the related art.

The invention provides a full waveform data decomposition method of a satellite-borne large-spot laser radar, which is characterized by comprising the following steps of:

step S3, detecting superposition Gaussian components at two sides of the effective obvious peak point, and estimating parameter initial values of the superposition Gaussian components;

step S4, updating the number of gaussian components included in the filtered echo waveform, and performing an optimized estimation on the parameters of the gaussian components.

Preferably, step S1 specifically includes the following sub-steps:

s1.1, estimating the background noise level of each original echo waveform in the satellite-borne large-spot laser radar full waveform data;

and S1.3, performing one-dimensional Gaussian filtering processing on the effective original echo waveform to obtain a filtered echo waveform.

Preferably, in step S1.1, the background noise level of the original echo waveform comprises an average value μ of the background noise of the original echo waveform_nAnd standard deviation of background noise delta_nWherein:

mean value of background noise mu_nThe formula of (c) is shown as follows:

Standard deviation delta of background noise_nThe formula of (c) is shown as follows:

Preferably, in step S1.1, the threshold (th) for background noise is set according to the following formula:

th＝μ_n+m·δ_n

where th is the background noise threshold, μ_nIs the average value of the background noise, δ_nM is a constant value (preferably, m is 4.5) as the standard deviation of the background noise.

Preferably, step S1.3 specifically comprises the following sub-steps:

W＝1+2·ceil(3·σ)

step S1.3.2, calculating a gaussian convolution kernel of the one-dimensional gaussian filter, specifically as follows:

C＝floor(W/2)+1

The values in the convolution kernel of gaussian are normalized, and the specific calculation is shown as the following formula:

Specifically, one-dimensional gaussian filtering is performed on all sampling points in the effective original echo waveform according to the following formula:

Preferably, step S2 specifically includes the following sub-steps:

s2.1, traversing the filtered echo waveform through a neighborhood window to detect an obvious peak point in the filtered echo waveform;

s2.3, judging whether the obvious peak point is effective or not based on the detected obvious peak point and an inflection point corresponding to the obvious peak point;

and S2.4, estimating the initial value of the parameter of the Gaussian component corresponding to the determined effective obvious peak value.

Preferably, step S2.3 specifically comprises the following sub-steps:

step S2.3.2, calculating the average value of the time coordinates of all left inflection points of the obvious peak point;

wherein, mu^LIs the average of the time coordinates of all left inflection points of the distinct peak point, m is the total number of left inflection points of the distinct peak point,

the time coordinate value of the k-th inflection point on the left side of the obvious peak point.

d_L＝|μ^L-X(i)|

wherein d is_LThe time distance from the left inflection point to the obvious peak point, X (i) is the corresponding time coordinate of the obvious peak point Y (i), mu^LThe average of the time coordinates of all left inflection points of the apparent peak point.

wherein d is_LAnd the FWHM is the full width at half maximum of the transmitted waveform of the satellite-borne laser radar.

Step S2.3.5, calculating the average value of the time coordinates of all right inflection points of the obvious peak point;

wherein, mu^RIs the average of the time coordinates of all right inflection points of the distinct peak point, n is the total number of right inflection points of the distinct peak point,

the time coordinate value of the k-th inflection point to the right of the sharp peak point.

d_R＝|μ^R-X(i)|

wherein d is_RAnd the FWHM is the full width at half maximum of a transmitting waveform of the satellite-borne laser radar.

the time coordinate value of the k-th inflection point to the left of the obvious peak point,the time coordinate value of the k-th inflection point to the right of the sharp peak point.

wherein d is_LDistance in time from left inflection point to distinct peak point, d_RThe time distance from the right inflection point to the obvious peak point, X (i) is the corresponding time coordinate of the obvious peak point Y (i), mu^LAverage of the time coordinates of all left-hand inflection points of the significant peak point, μ^RThe average of the time coordinates of all right-hand inflection points of the apparent peak point.

Preferably, step S3 specifically includes the following sub-steps:

step S3.2, based on the set interval, detecting the superimposed Gaussian component on the left side of the effective obvious peak point according to the following formula:

and y (j1 < y (k1)

step S3.3, based on the set interval, detecting the superimposed Gaussian component on the right side of the effective obvious peak point according to the following formula:

and y (j2 > y (2) (26)

S3.4, estimating an initial value of the parameter of the superposition Gaussian component; step S3.4 specifically includes the following substeps:

and S3.5, calculating the number of Gaussian components contained in the filtered echo waveform.

Preferably, step S4 specifically includes the following sub-steps:

and S4.1.1, establishing an error function and a target function between a fitting function and the observation value selected from the effective original echo waveform (namely the actual observation waveform).

Step S4.2, combining threshold th of fitting residual error RMSE_RMSEAnd performing further iterative optimization on the parameter values of the detected Gaussian components in the filtered echo waveform. Wherein, the step S4.2 specifically comprises the following substeps:

step S4.2.1 is to calculate the residual error RMSE between the amplitude values of the sampling points in the effective original echo waveform and the fitting values of the amplitude values of the sampling points obtained after the LM algorithm is iteratively optimized.

The invention also provides a computer program characterized in that the method of claims 1-9 is performed by executing the computer program.

Due to the adoption of the technical scheme, compared with the prior art, the invention has the following advantages and positive effects:

(1) the initial detection strategy of the Gaussian component is improved by combining the inflection point matching and the iterative optimization algorithm, and the initial detection strategy comprises detecting the Gaussian component corresponding to the obvious peak point and superposing the Gaussian component, so that the better initial estimation of the Gaussian component of the echo waveform is obtained. On the basis, an optimization estimation strategy of Gaussian component parameters is improved, in the parameter iteration optimization process of the Gaussian components, a threshold value of fitting residual errors RMSE is set, the missing Gaussian components (including the Gaussian components in the complex multi-peak waveforms) are continuously searched in an iteration mode, constraint conditions for eliminating invalid Gaussian components are set, and finally the decomposition result of the Gaussian components is more reliable, effective and accurate.

(2) For some original echo waveforms, before gaussian filtering, the maximum amplitude value of an original echo waveform sampling point is greater than the calculated background noise threshold, and after gaussian filtering, the maximum amplitude value of the echo waveform sampling point is less than or equal to the background noise threshold. This makes the filtered echo waveform unusable for subsequent detection of significant peak and inflection points. For the situation, the original effective echo waveform before filtering is reserved and used as the filtered echo waveform to be used for subsequent detection of the obvious peak point and the inflection point, so that the effective utilization rate of the original echo waveform is improved.

(3) The invention removes noise from the original echo waveform through Gaussian filtering to reduce the influence of background noise, detects the obvious peak point and inflection point of the obtained filtered echo waveform and matches the inflection point, can detect the Gaussian component of the ground object target which is easy to distinguish, mainly corresponds to the ground covering areas such as bare land, farmland areas, building areas and the like, and detects the superposed Gaussian component of the ground object target through matching the inflection points at the two sides of the obvious peak point, which corresponds to the areas of the ground object target which are close to each other in the vertical direction, mainly corresponds to the ground covering areas such as the building areas, vegetation areas and the like with similar vertical height, and completes parameter initial value estimation of the Gaussian component detected in the filtered echo waveform through the detection. And then, carrying out parameter iterative optimization on the Gaussian components through a nonlinear iterative optimization algorithm, in the process, searching the Gaussian components of the missed ground object target through fitting a threshold value of residual RMSE, and removing invalid Gaussian components through an additional constraint condition. Based on the process, in the area with a complex ground object structure, the method can also finish the detection and decomposition of the Gaussian component of the ground object target, such as a forest vegetation area and the like. Therefore, the method provided by the invention has certain applicability and reliability under different surface and ground object coverage conditions, such as urban building areas, farmland areas, forest vegetation areas and the like.

(4) The method comprises the steps of preprocessing an original echo waveform by Gaussian filtering and the like, detecting an obvious peak point and an inflection point of the obtained filtered echo waveform and matching the inflection point, finishing detection of Gaussian components in the filtered echo waveform and parameter initial value estimation of corresponding Gaussian components, performing optimization estimation of Gaussian component parameters by a nonlinear iterative optimization algorithm, searching for the Gaussian components of a leaked ground object target by fitting a threshold value of residual RMSE in the optimization process, and removing invalid Gaussian components by combining an additional constraint condition. By the means, the difficult problems in the original echo waveform data decomposition, such as the decomposition problems of the original echo waveform containing the superimposed Gaussian component, the original echo waveform containing complex multiple peaks and the original echo waveform containing the weak echo component, can be effectively solved.

Drawings

Fig. 1 is a flowchart of a full waveform data decomposition method of a satellite-borne large-spot lidar according to an embodiment of the invention.

Fig. 2(a) is a transmission waveform of the GLAS altimeter system according to an embodiment of the present invention.

Fig. 2(b) shows an original echo waveform corresponding to a transmitted waveform, which has been compressed during reception, according to an embodiment of the present invention.

Fig. 2(c) is a diagram illustrating the original echo waveform after decompression and the background noise estimation according to the embodiment of the invention.

FIG. 2(d) is a diagram illustrating the result of Gaussian filtering according to the embodiment of the invention.

Fig. 2(e) is a schematic diagram illustrating a waveform decomposition result corresponding to an original echo waveform after decompression according to an embodiment of the present invention.

Fig. 3(a) is a waveform decomposition result corresponding to the decompressed original echo waveform, wherein the decomposition result of the superimposed original echo waveform is shown.

Fig. 3(b) is a waveform decomposition result corresponding to the decompressed original echo waveform according to the embodiment of the present invention, wherein the decomposition result of the original echo waveform containing a weak echo component is shown.

Fig. 3(c) is a waveform decomposition result corresponding to the decompressed original echo waveform according to the embodiment of the present invention, wherein the decomposition result of the original echo waveform of the complex multiple peaks is shown.

FIG. 4(a) is a graph of R of a large number of waveform decomposition results of Beijing plain test area according to an embodiment of the present invention²And (6) counting the histogram.

FIG. 4(b) is a graph showing R of a result of decomposing a large number of waveforms in a Daxing-AnLing test area of Heilongjiang province according to an embodiment of the present invention²And (6) counting the histogram.

Detailed Description

The present invention is described more fully hereinafter with reference to the accompanying drawings, in which exemplary embodiments of the invention are shown.

The embodiment of the invention provides a full waveform data decomposition method of a satellite-borne large-spot laser radar in combination with an inflection point matching and iterative optimization algorithm, which mainly comprises the following steps of:

The satellite-borne large-spot laser radar full-waveform data can be downloaded and obtained from foreign official websites, for example, GLAS data can be downloaded and obtained from national snow and ice data center (NSDIC) websites of America, GLA01 data in the GLAS data provided by the websites comprise satellite-borne full-waveform laser radar data, and a GLA14 product comprises result parameters after the full-waveform data is processed. The full waveform data of the satellite-borne large-spot laser radar comprises a certain amount of original echo waveforms to be processed, and the original echo waveforms are influenced by background noise in the receiving process, so that each original echo waveform in the full waveform data of the satellite-borne large-spot laser radar needs to be preprocessed respectively. The pretreatment operation specifically comprises the following steps: and judging whether each original echo waveform is an effective original echo waveform, and performing one-dimensional Gaussian filtering processing on the effective original echo waveform to obtain a filtered echo waveform.

Specifically, step S1 specifically includes the following sub-steps:

and S1.1, estimating the background noise level of each original echo waveform in the satellite-borne large-spot laser radar full waveform data.

The original echo waveform is composed of N discrete sampling points, wherein the sampling points of the beginning part and the ending part are usually generated by random fluctuation of background noise, namely belong to a background noise signal, in order to estimate the background noise level, the traditional method generally selects the sampling points of the beginning part in the original echo waveform to estimate the background noise, and the improvement of the invention is that each N sampling points (N is a natural number) of the beginning part and the ending part are respectively selected from the sampling points of the original echo waveform. The n sampling points of the selected beginning part and the ending part are combined and then are uniformly counted, and the background noise level can be estimated more fully and reasonably by the selection mode. Preferably, n is 20, that is, 2n (40) sampling points are selected in total, and the average value (mu) is calculated_n) And standard deviation (delta)_n) And will mu_nAnd delta_nAs the mean and standard deviation, respectively, of the background noise.

In particular, the background noise level of the original echo waveform comprises an average μ of the background noise of the original echo waveform_nAnd standard deviation of background noise delta_nWherein:

mean value of background noise mu_nIs calculated byAs shown in the following formula:

wherein N is the sum of the number of sampling points in a selected original echo waveform, N is the number of sampling points respectively selected at the beginning part and the ending part in the sampling point of the original echo waveform, and V_iIs the amplitude of the ith sample point.

Standard deviation delta of background noise_nThe formula of (c) is shown as follows:

From the above calculations, a good estimate of the mean and standard deviation of the background noise, μ_nAnd delta_nThese two parameter values may represent the level of background noise.

And S1.2, setting a background noise threshold based on the background noise level, comparing the maximum amplitude value of the sampling point of the original echo waveform with the background noise threshold to judge whether the original echo waveform is an effective original echo waveform, and only retaining the effective original echo waveform.

In this step, a background noise threshold is set for judging an original echo waveform as a whole as an effective original echo waveform or a noise echo waveform.

The threshold (th) for background noise is set according to the following formula:

th＝μ_n+m·δ_n(3)

where th is the background noise threshold, μ_nIs the average value of the background noise, δ_nM is a constant value (empirical value, according to the experimental sum) for the standard deviation of the background noiseActual demand set). In order to reduce the influence of background noise on the detection and subsequent processing of the echo waveform, it is preferable that m be 4.5.

If the maximum amplitude value of the sampling point in the original echo waveform is larger than the background noise threshold (th), the original echo waveform can be judged as an effective original echo waveform, otherwise, the original echo waveform is judged as a noise echo waveform, and only the effective original echo waveform is reserved. It should be noted that: since each original echo waveform is independent, each echo waveform individually corresponds to a background noise threshold.

According to the preferred embodiment of the present invention, in consideration of the difference between the software and hardware of the satellite laser load, for each original echo waveform, if the echo waveform has been compressed based on the redundancy of the reduction data when it is acquired, the original echo waveform data needs to be decompressed first before being preprocessed, and the original echo waveform data is restored to a complete echo signal through the decompression.

For example, the echo waveform data of GLAS is obtained after compression processing, so decompression processing needs to be performed first to recover a complete echo signal. Specifically, the original echo waveform data is decompressed according to the following formula:

wherein p, q and L are Npq type compression constants, i is the serial number of the sampling point in the original echo waveform, N_gatesAnd t (i) is the sequence number of the ith sampling point in the original echo waveform before compression.

Before data compression, each original echo waveform of the GLAS comprises N (the system given value is 1000) sampling points, and according to the Npq type compression principle, for 1-N sampling points, a plurality of continuous sampling points are replaced by the average value of the sampling points until the number of the sampling points of the echo waveform reaches N_gates(system setpoint 544), if 1000 sample points are not sufficient to fill the entire echo waveform, the remaining echo waveform is filteredThe remaining sample points are filled with 0.

For the original echo waveform after compression processing, from 1 to L-1 sampling points, the amplitude value of each sampling point is obtained by averaging continuous p sampling points of the original echo waveform before compression, from L to N_gatesAnd the amplitude value of each sampling point is obtained by averaging the continuous q sampling points of the original echo waveform before compression.

Let x (i) be the time coordinate corresponding to each original echo waveform data of GLAS before compression, where i is 1,2, …, N, and the corresponding amplitude value is y (i), where i is 1,2, …, N. Setting the time coordinate corresponding to the original echo waveform data obtained by compression processing as X_Npq(i)Where i is 1,2, …, N_gatesAnd the corresponding amplitude value is Y_Npq(i)Where i is 1,2, …, N_gates. Specifically, according to the above equation (4), the original echo waveform data recovered by the decompression process is:

X[t(i)]＝X_Npq(i) where i is 1,2, …, N_gates。

Y[t(i)]＝Y_Npq(i) Where i is 1,2, …, N_gates。

Based on the above assignments, for assignments of other sampling points of the original echo waveform before compression, when i is 1,2, …, N_gatesWhen the temperature of the water is higher than the set temperature,

if i is less than or equal to L-1, then:

if L is less than or equal to i is less than or equal to N_gatesThen, there are:

wherein p, q and L are Npq type compression constants, i is the serial number of the sampling point in the original echo waveform, N_gates(the system given value is 544) is the total number of sampling points of the original echo waveform, and t (i) is the ith sample in the original echo waveformNumber of dots before compression, X [ t (i)]Is the time coordinate corresponding to the t (i) th sampling point of the original echo waveform before compression, Y [ t (i)]The amplitude value corresponding to the t (i) th sampling point of the original echo waveform before compression. And according to the decompression process, restoring and obtaining each original echo waveform before compression.

Fig. 2(a) shows a transmission waveform of a GLAS altimeter system according to an embodiment of the present invention, and fig. 2(b) shows an original echo waveform corresponding to the transmission waveform according to the embodiment of the present invention, which is compressed during reception. Fig. 2(c) is a schematic diagram of an original echo waveform after decompression and background noise estimation according to an embodiment of the present invention, specifically, fig. 2(c) is an original echo waveform obtained by decompressing the original echo waveform of fig. 2(b), calculating a background noise average value and a standard deviation of the waveform, and estimating a noise threshold, where a black solid line is an estimated background noise threshold.

And S1.3, performing one-dimensional Gaussian filtering processing on the effective original echo waveform to obtain a filtered echo waveform.

Each of the original valid echo waveforms may appear to jitter in amplitude due to the influence of background noise, and if the valid original echo waveforms are directly used for the detection of the significant peak and inflection points in the subsequent step S2, a large number of significant peak and inflection points having low amplitude and low pulse width are obtained, and these significant peak and inflection points are substantially caused by the background noise, thereby affecting the detection result of the gaussian component included in the valid original echo waveform. Therefore, it is necessary to smooth the effective echo waveform by performing one-dimensional gaussian filtering processing on the effective echo waveform, so as to reduce the influence of background noise on the effective echo waveform.

The gaussian filtering is a low-pass filter, and the principle of the one-dimensional gaussian filtering is to perform convolution processing on each effective original echo waveform and a gaussian function, suppress high-frequency background noise, and realize smoothing of the effective original echo waveform, so as to obtain the echo waveform after the gaussian filtering.

The one-dimensional gaussian filter function g (x, σ) is given by:

where σ is the filter width, the parameter is variable, and t is the filter argument on the abscissa axis.

Let f (t) represent the signal to be processed, the one-dimensional gaussian filtered signal s (t, σ) is expressed as:

s(t,σ)＝f(t)*g(t,σ) (6)

where, is the convolution operator, σ is the filter width, and t is the position of the signal to be processed.

Step S1.3 specifically includes the following substeps:

step S1.3.1, determine the width of the one-dimensional gaussian filter and the window size of the one-dimensional gaussian filter.

For the satellite-borne large-spot laser radar, since the orbit of the satellite is high and is influenced by the atmosphere and the earth surface, the background noise has a large influence on the received echo waveform, so that the signal-to-noise ratio of the echo signal of the satellite-borne laser radar is low, and therefore, a larger filter width needs to be selected to reduce the influence of the background noise, for this reason, preferably, the width (i.e., σ) of the one-dimensional gaussian filter is selected to be the full width at half maximum (FWHM) of the transmission waveform of the satellite-borne laser radar. The selection avoids the poor denoising effect caused by the small width of the one-dimensional Gaussian filter and the insufficient smoothness degree of the echo waveform, and also avoids the loss of the distance precision of the echo waveform caused by the excessive smoothness of the original echo waveform due to the excessive width of the filter. For example, for a GLAS altimetry system, where the full width at half maximum (FWHM) value of the transmit waveform is 4ns, the one-dimensional Gaussian filter width is selected to be 4 ns.

For the gaussian function, after the one-dimensional gaussian filter width σ is determined, the window size of the one-dimensional gaussian filter needs to be determined. The Gaussian function is a bell-shaped curve, the bell-shaped curve being in the interval (mu)_n-3σ,μ_nArea in the range of +3 σ) accounts for 99.7% of the total area under the curve, where μ_nIs the mean of a gaussian function, and therefore a one-dimensional gaussian filter with a radius of 3 sigma is chosenA window of the one-dimensional gaussian filter can thus almost contain the effective range of values of the gaussian function.

The size of the one-dimensional gaussian filter window W is calculated as:

W＝1+2·ceil(3·σ) (7)

where ceil (x) is a function that functions to return a minimum integer greater than or equal to a specified expression, W is the window size of the one-dimensional gaussian filter, and σ is the width of the one-dimensional gaussian filter, where preferably the width σ of the one-dimensional gaussian filter is selected to be the full width at half maximum (FWHM) of the transmit waveform.

Step S1.3.2, calculate the gaussian convolution kernel of the one-dimensional gaussian filter.

The one-dimensional gaussian filter width and the filter window size are determined through the step S1.3.1, in this step, the gaussian template center position C is calculated by using the two parameters, and a gaussian convolution kernel k (i) of the one-dimensional gaussian filter is established through the filter function g (x, σ) in the formula (5), where i ═ 1,2, …, W, and the convolution kernel of gaussian is a discrete approximation to continuous gaussian. The specific calculation is shown as the following formula:

C＝floor(W/2)+1 (8)

After the computation is completed, the values in the convolution kernel of the gaussian are normalized, and the purpose of this is to ensure that the sum of the values in the convolution kernel of the gaussian must be 1.

The specific calculation is shown as the following formula:

Through the above steps S1.3.1 and S1.3.2, the values of the gaussian convolution kernel required for the one-dimensional gaussian filtering are obtained, and these values are used as input in the one-dimensional gaussian filtering of this step.

Specifically, one-dimensional gaussian filtering is performed on all sampling points in the effective original echo waveform according to the following formula:

wherein y (i) is an amplitude value of an ith sampling point of an effective original echo waveform before filtering, y (i) is an amplitude value of an ith sampling point of the filtered echo waveform, L is floor (W/2), W is a window size (floor (x)) of a one-dimensional gaussian filter, and the function is to take a maximum integer not greater than x, and L is obtained by rounding half of the window size of the filter, so that L can be regarded as the radius of the filter window, K (L + j +1) is a gaussian convolution kernel of one-dimensional gaussian filtering, and N is the total number of sampling points in the effective original echo waveform.

Fig. 2(d) shows the filtered echo waveform obtained by performing one-dimensional gaussian filtering on the effective original echo waveform, which can be seen to achieve the effect of smoothing the echo waveform and reduce the influence of background noise.

The echo waveform after one-dimensional gaussian filtering obtained through the above steps is used for the detection of the significant peak point and the inflection point in step S2, so that the number of the detected significant peak points and inflection points having low amplitude and low pulse width can be reduced, and the number of the detected significant peak points and inflection points can be reduced as a whole, thereby improving the accuracy of the result of the following step S2.

In a preferred embodiment of the present invention, after step S1.3, the following step S1.4 is further included:

step S1.4, calculating the maximum amplitude value of the sampling point of the filtered echo waveform, comparing the maximum amplitude value with the background noise threshold of the echo waveform, if the maximum amplitude value is less than or equal to the background noise threshold, retaining the original effective echo waveform before filtering and using the original effective echo waveform as the filtered echo waveform (namely, not performing filtering processing on the original effective echo waveform, but directly using the original effective echo waveform before filtering as the filtered echo waveform), otherwise, retaining the filtered echo waveform.

It is noted that for some echo waveforms, the maximum amplitude value of the original echo waveform sampling point is greater than the noise threshold before gaussian filtering, and the maximum amplitude value of the echo waveform sampling point is equal to or less than the noise threshold after gaussian filtering. Since each original echo waveform is independent, each echo waveform is individually corresponding to a background noise threshold. For the above situation, detection of an obvious peak point and an inflection point in the subsequent steps cannot be performed using the echo waveform after gaussian filtering. For this purpose, the following is specifically dealt with:

the maximum amplitude value in the sample points of the filtered echo waveform data y (i) (where i ═ 1,2, …, N) is calculated by the max () function, where max () represents the maximum value of the set of data calculated, i.e. max (y) is found.

(1) If max (y) ≦ th, the valid original echo waveform data y (i) (where i ═ 1,2, …, N) before gaussian filtering is retained for detection of significant peak points and detection of inflection points on the left and right sides of significant peak points in step S2 described below.

(2) If max (y) > th, the filtered echo waveform data y (i) (where i ═ 1,2, …, N) is retained for detection of an apparent peak point and detection of inflection points on the left and right sides of the apparent peak point in step S2 described below. In most cases, the original echo waveform is the case.

By performing step S1.4, the utilization of the effective original echo waveform can be further improved, which is also one of the improvement points of the present invention.

Step S2, detecting an obvious peak point and an inflection point corresponding to the obvious peak point in the filtered echo waveform, determining whether the detected obvious peak point is valid based on the detected obvious peak point and the inflection point corresponding to the obvious peak point, and estimating an initial value of a parameter of a gaussian component corresponding to the valid peak point.

Since the more easily detected sharp peak in the filtered echo waveform is a point corresponding to a gaussian component of a potential clutter object. Therefore, in this step, the significant peak point of the filtered echo waveform is first detected, and on this basis, the inflection points corresponding to the left and right sides of the significant peak point are detected according to the principle that the second derivative of the sampling point in the filtered echo waveform has different signs at both sides of the inflection point. Generally, one distinct peak point matches with the left and right inflection points, and considering the complexity of the echo waveform, the inflection points should be present on at least one side of the two sides of the distinct peak point. In order to avoid the influence of background noise, the detection of an obvious peak point and an inflection point is only carried out on the sampling points of which the amplitude values are above a noise threshold.

The Gaussian components contained in the effective original echo waveform can be reflected by detecting the obvious peak point and the inflection points on the left side and the right side of the filtered echo waveform, so that the number of the Gaussian components in the original echo waveform can be preliminarily judged according to the detection result of the filtered echo waveform, and the initial values of the parameters such as the amplitude, the central position, the standard deviation and the like of the Gaussian components contained in the effective original echo waveform are estimated by using the obvious peak point and the matched inflection point detected in the filtered echo waveform.

Step S2 specifically includes the following substeps:

and S2.1, traversing the filtered echo waveform through a neighborhood window to detect an obvious peak point in the filtered echo waveform.

Step S2.1 is specifically as follows: for the filtered echo waveform y (i) (where i is 1,2, …, N), using a local maximum method, selecting 2r +1 sampling points as a neighborhood window size (where r is a radius of the selected neighborhood window, and preferably, r is 2), traversing the filtered echo waveform with the neighborhood window, if the amplitude values of 2r +1 sampling points of the filtered echo waveform covered by the neighborhood window are all greater than a background noise threshold of the echo waveform (it should be noted that, since each original echo waveform is independent, each echo waveform individually corresponds to a background noise threshold), the ith sampling point of the filtered echo waveform is located at the center of the neighborhood window, and the amplitude value y (i) of the ith sampling point is the largest among the amplitude values of the 2r +1 sampling points and is on the left side of the ith sampling point, starting from the (i-1) th sampling point, sequentially leftwards, and gradually reducing the amplitude value corresponding to the corresponding sampling point; and on the right side of the ith sampling point, the amplitude value of the (i +1) th sampling point is less than or equal to the amplitude value of the ith sampling point, and the amplitude values corresponding to the corresponding sampling points are gradually reduced from the (i +1) th sampling point to the right in sequence. If the above conditions are satisfied, the ith sampling point is a local maximum point, that is, the ith sampling point is determined to be an obvious peak point.

Specifically, the apparent peak point is determined according to the following formula:

wherein i is 3,4, …, N-2(12)

Wherein, y (i) is an amplitude value of an ith sampling point of the filtered echo waveform, r is a radius of a neighborhood window, th is a background noise threshold, and N is a total number of sampling points of the filtered echo waveform. And if the sampling point i of the filtered echo waveform meets the formula (12), determining the sampling point i as an obvious peak value point.

And S2.2, calculating second-order derivatives of the sampling points of the filtered echo waveform, and detecting inflection points on the left side and the right side of the obvious peak point according to the criterion that whether the second-order derivatives of two adjacent sampling points have opposite signs.

After a plurality of obvious peak points are detected, each obvious peak point corresponds to a gaussian component of a potential ground object target, and for the gaussian component, generally speaking, one obvious peak point and two inflection points, namely, the obvious peak point and the two inflection points are matched with one gaussian component. In consideration of the complexity of the echo waveform, inflection points should be present on at least one side of both sides of the significant peak point.

The improvement of the invention is that: and performing inflection point detection on the left side and the right side of each obvious peak point, calculating second-order derivatives of sampling points of the filtered echo waveform, and detecting inflection points in the filtered effective echo waveform according to the criterion whether the second-order derivatives of two adjacent sampling points have opposite signs. The general inflection point detection method includes: according to the principle that the concave-convex property of curves on two sides of an inflection point is opposite, for a sampling point to be detected, neighborhood windows (such as 6 sampling points) with a certain size are selected in the left side area and the right side area of the sampling point, the sampling points are divided into two parts, each part is provided with 3 sampling points, in the 3 sampling points, the sum of the amplitude values of the sampling points at two ends is subtracted from twice of the amplitude value of the sampling point in the middle, then the calculation result values obtained by the two parts are multiplied, and if the multiplication result is smaller than 0, the inflection point is found. Compared with a common inflection point detection method, the inflection point detection method is more accurate and effective.

Specifically, inflection points on the left and right sides of each obvious peak point are detected according to the following formula:

y '(i). Y' (i +1) < 0 and Y (i) > th and Y (i +1) > th (14)

Wherein, Y (i) is an amplitude value of an ith sampling point of the gaussian-filtered echo waveform, Y ″ (i) is a second derivative value of the ith sampling point of the gaussian-filtered echo waveform, and th is a background noise threshold; if the second derivatives Y "(i) and Y" (i +1) of two adjacent sampling points Y (i) and Y (i +1) on both sides of the significant peak point, which are greater than the background noise threshold, satisfy the above equation (14), then the sampling point corresponding to Y (i) is an inflection point.

The above equation (13) is a five-point numerical differential equation of the second derivative of the discrete points, which is used to calculate the second derivative value of the sampling points, wherein the constant values are all specified by the five-point numerical differential equation of the second derivative of the discrete points.

The above equation (14) is a discrimination formula of inflection points, and is used to detect inflection points on both sides of each significant peak point. And in the detection process, only the sampling point of the filtered echo waveform which is greater than the background noise threshold th is considered, and the sampling point of the filtered echo waveform which is less than or equal to the background noise threshold th is not considered. It should be noted that: since each original echo waveform is independent, each echo waveform individually corresponds to a background noise threshold. The inflection point is calculated by using the formula (14) according to the principle that the second derivative of the Gaussian component has opposite signs on both sides of the inflection point.

If the second derivatives Y "(i) and Y" (i +1) of the adjacent two sampling points Y (i) and Y (i +1) greater than the background noise threshold in the monotonous region on the left of the obvious peak point satisfy the above equation (14), then Y (i) corresponds to the sampling point being an inflection point and the corresponding time coordinate value being x (i). Similarly, if the second derivatives Y "(i) and Y" (i +1) of two sampling points corresponding to adjacent Y (i) and Y (i +1) in the monotonous region on the right side of the sharp peak point satisfy the above equation (14), the sampling point corresponding to Y (i) is an inflection point whose time coordinate value is x (i). And if the amplitude value corresponding to one inflection point is Y (i), the time coordinate value corresponding to the inflection point is X (i), and all inflection points on two sides of the obvious peak point are found in the mode.

And (4) completing the detection of all obvious peak points and inflection points by utilizing the steps. For each gaussian component significant peak point, there may be more than one inflection point detected on both sides, and assuming that m inflection points are detected in the monotonous region on the left side of the significant peak point, the time values of the inflection points on the abscissa are: x (i1), X (i2), …, X (im) (i1, i2, …, im is the serial number of the detected inflection point in the X data), which can be expressed as:

and n inflection points are detected in a monotonous area on the right side of the obvious peak point, and the time values of the inflection points on the abscissa are respectively as follows: x (j1), X (j2), …, X (jn) (j1, j2,…, jn is the sequence number of the detected inflection point in the X data), which can be expressed as:

and S2.3, judging whether the obvious peak point is effective or not based on the detected obvious peak point and an inflection point corresponding to the obvious peak point.

After the above steps are completed, all the obvious peak points and the corresponding inflection points are detected as the input data of the step. Since an obvious peak point does not necessarily correspond to the gaussian component of the ground object target, that is, the obvious peak point may be invalid, the validity of the obvious peak point is determined.

The method specifically comprises the following substeps:

and if the left inflection point and the right inflection point on the two sides of the obvious peak point do not exist, judging that the obvious peak point does not correspond to the Gaussian component, judging the obvious peak point to be an invalid obvious peak point, and discarding the invalid obvious peak point. The principle of making the above judgment is that for the gaussian component of the ground object target, an obvious peak point should be included, and the left and right sides of the obvious peak point respectively correspond to a left inflection point and a right inflection point. In practical situations, under the influence of background noise, a plurality of inflection points may appear on the left side and the right side of the distinct peak point, and if the left inflection point and the right inflection point on both sides of the distinct peak point do not exist in the filtered echo waveform, the distinct peak point is likely to be generated by the background noise.

Step S2.3.2, calculating the average value of the time coordinates of all left inflection points of the obvious peak point;

wherein, mu^LIs the average of the time coordinates of all left inflection points of the distinct peak point, m is the total number of left inflection points of the distinct peak point,

the time coordinate value of the k-th inflection point on the left side of the obvious peak point.

d_L＝|μ^L-X(i)|(16)

wherein d is_LAnd the FWHM is the full width at half maximum of the transmitted waveform of the satellite-borne laser radar.

Step S2.3.5, calculating the average value of the time coordinates of all right inflection points of the obvious peak point;

wherein, mu^RAll right of the significant peak pointThe average value of the time coordinates of the side inflection points, n is the total number of right side inflection points of the obvious peak point,the time coordinate value of the k-th inflection point to the right of the sharp peak point.

d_R＝|μ^R-X(i)| (19)

wherein d is_RAnd the FWHM is the full width at half maximum of a transmitting waveform of the satellite-borne laser radar.

And S2.3.8, respectively calculating the average value of the time coordinates of all left inflection points and all right inflection points of the obvious peak point.

If the number of the detected inflection points in the monotonous area on the left side of the obvious peak point is more than or equal to 1, the obvious peak point has a left inflection point, otherwise, the obvious peak point does not have a left inflection point. If the obvious peak point has a left inflection point, taking the average value of the time coordinates of all inflection points on the left side of the obvious peak point as the time coordinate of the left inflection point of the obvious peak point on the abscissa, and recording the time coordinate as mu^L. Similarly, if the monotonic region on the right side of the obvious peak point is detectedIf the number of the inflection points is more than or equal to 1, the obvious peak point exists at the right inflection point, otherwise, the obvious peak point does not exist at the right inflection point. If the obvious peak point exists at the right inflection point, the average value of the time coordinates of all inflection points at the right side of the obvious peak point is taken as the time coordinate of the right inflection point of the obvious peak point on the abscissa and is marked as mu^R. The following formula is specifically calculated:

the time coordinate value of the k-th inflection point to the left of the obvious peak point,

the time coordinate value of the k-th inflection point to the right of the sharp peak point.

At step S2.3.9, the time distances from the left inflection point and the right inflection point to the distinct peak point are calculated, respectively.

Specifically, the time distances from the left inflection point and the right inflection point to the obvious peak point are calculated according to the following formula:

If the peak point is judged to be a valid obvious peak point, the peak point corresponds to a Gaussian component of a potential surface feature target. According to the preferred embodiment of the present invention, it should be noted that, for the filtered echo waveform and satisfying max (y) > th, if the steps of detecting the significant peak point and detecting the inflection points at the left and right sides of the significant peak point are performed, the corresponding gaussian component is not found, and for this case, the method of the present invention is as follows:

calculating the maximum amplitude value of the sampling point of the filtered echo waveform, namely max (Y), and considering that the position of the maximum amplitude value corresponds to a Gaussian component, namely, calculating the amplitude value, the time coordinate and the amplitude value of the sampling point corresponding to the maximum amplitude value

The (estimated value of standard deviation) value is used as an initial parameter of the gaussian component (where FWHM is the full width at half maximum of the emission waveform). This is due to the fact that the numerical relationship between the standard deviation σ and the FWHM of the theoretical gaussian component is

And enters the optimized estimation of the gaussian component of the subsequent step.

In addition, it should be noted that the correspondence between the transmit waveform and the echo waveform is as follows: the echo waveform is formed after response between the transmitting waveform and the ground object targets, the echo waveform is the superposition sum of Gaussian components corresponding to the ground object targets, and each Gaussian component in the echo waveform reflects the characteristics of the corresponding ground object target. In general, the full width at half maximum of each gaussian component in the echo waveform is equal to or greater than the full width at half maximum of the transmit waveform.

And S2.4, estimating the initial value of the parameter of the Gaussian component corresponding to the determined effective obvious peak value.

After the above steps, the effective significant peak point is determined and used as the input data of the step. For each determined valid significant peak point, there is a gaussian component corresponding to the clutter object in a potential valid raw echo waveform. And (5) assuming that the corresponding sampling point of Y (i) is an obvious peak value point, estimating the initial value of the parameter of the corresponding Gaussian component.

Since the original echo waveform of a single ground object target is approximately gaussian distributed, it can be expressed mathematically through a gaussian function, while there may be multiple ground objects on the transmission path of the satellite laser pulse, and the echo waveform received by the instrument is the sum of the echoes of multiple ground object targets, for this purpose, each received original echo waveform is expressed through a gaussian mixture model, as shown in the following formula:

theoretically, the original echo waveform can be regarded as the superposition of a plurality of gaussian functions, which can be expressed by a gaussian mixture model, and therefore, the received original echo waveform is expressed as the superposition of k gaussian components

Wherein y is an approximate waveform of the received original echo waveform, t is time coordinate data corresponding to a sampling point of the original echo waveform received in a discretization form, k is the number of Gaussian components in the original echo waveform, A_i、μ_i、σ_iRespectively, the amplitude, the central position (i.e. time coordinate) and the standard deviation of the ith Gaussian component, wherein epsilon is the background noise of the original echo waveform, and the average value mu of the background noise can be used_nAnd (4) directly estimating. On the basis of establishing a Gaussian mixture model, the full waveform data decomposition method provided by the invention is utilized to complete the fitting of the waveform and the decomposition of the Gaussian component.

For the Gaussian mixture model in equation (23), k and A_i、μ_i、σ_i(i＝1，2, …, k) are final parameters required to be obtained by the original echo waveform processing, and the parameter values are unknown, and since the gaussian mixture model is a nonlinear model, the parameters need to be iterated and optimized continuously by using a nonlinear least square method in subsequent steps, and finally the optimal values of the parameters are obtained. Before the iteration and optimization of the parameters, the initial value of the number k of the gaussian components contained in the received original echo waveform needs to be estimated, and the parameter a of each gaussian component_i、μ_i、σ_iAnd (i ═ 1,2, …, k), and these initial values are used as initial input parameters.

Based on the above principle, the gaussian component corresponding to the effective distinct peak point is one of gaussian components contained in the original echo waveform, and the mathematical expression of the gaussian component is as follows:

wherein the parameter A of the Gaussian component_i、μ_i、σ_iAre unknown and therefore the initial values of these parameters need to be estimated first.

The method comprises the following specific steps:

amplitude A of the Gaussian component corresponding to the significant peak point_iIs estimated as the amplitude value y (i) of the effective significant peak point;

the central position mu of the Gaussian component corresponding to the effective obvious peak point_iThe initial value of the sampling point is estimated as the time coordinate value of the ith sampling point corresponding to the effective obvious peak point in the filtered echo waveform;

standard deviation sigma of Gaussian component corresponding to effective obvious peak point_iThe initial value of D is estimated as follows:

(1) if it is not

And the right turning point is not present, D ═ D_L；

(2) If it is not

And the left inflection point does not exist, D ═ D_R；

(3) If it is notAnd isThen D ═ D_L；

(4) If it is not

And is

Then D ═ D_R；

(5) If it is not

And isFurther judge this if d_L≥d_RWhen D is D_ROtherwise D ═ D_L。

The above judgment principle is: and obtaining the effective obvious peak point according to the judgment of the obvious peak point in the steps. The position of each effective peak point corresponds to the Gaussian component of a ground object target contained in an effective original echo waveform, and the relationship between the emission waveform and the echo waveform indicates that the distance d between the left inflection point or the right inflection point of the Gaussian component corresponding to the effective peak point and the effective peak point is the distance d between the left inflection point or the right inflection point of the Gaussian component corresponding to the effective peak point and the effective peak point_LAnd d_RShould be numerically equal to or greater than half FWHM, but the distance value d takes into account the influence of background noise_LAnd d_RAt least 1 should be equal to or greater than half the FWHM.

By detecting the obvious peak point and the inflection point in the echo waveform after Gaussian filtering in the steps, the Gaussian component of the ground object target which is easy to distinguish can be detected, and the ground object target corresponds to ground surface coverage areas such as bare land, farmland areas, building areas and the like. And then finishing the initial judgment of the Gaussian component contained in the corresponding original echo waveform and the estimation of the parameter initial value of the corresponding Gaussian component.

And step S3, detecting the superposition Gaussian components at both sides of the effective obvious peak point, and estimating the parameter initial value of the superposition Gaussian component.

Through the step S2, effective significant peak points included in the echo waveform after gaussian filtering are detected, initial values of parameters such as amplitudes, center positions, standard deviations, and the like of gaussian components corresponding to all the effective significant peak points are estimated, and these results are used as input for the superimposed gaussian component detection in this step.

On both sides of the distinct peak point, there is a possible superimposed gaussian component, which needs to be further detected in the echo waveform. The superimposed gaussian component is generated because two ground object targets are close to each other, for example, when the distance between the two ground object targets is less than 1.5 meters, in the corresponding original echo waveforms, the gaussian components corresponding to the two ground object targets are superimposed together, and a superimposed echo waveform with only one obvious peak point may be generated. The echo waveform containing the superimposed gaussian component mainly corresponds to the ground surface coverage areas such as building areas and vegetation areas with similar vertical heights.

In this step, for the gaussian components corresponding to the effective significant peak points, intervals are set

[X(i)-D,X(i)+D](24)

Wherein, x (i) is a time coordinate value corresponding to the gaussian component corresponding to the effective distinct peak point, and D is an initial value of the standard deviation of the estimated gaussian component corresponding to the effective distinct peak point.

The interval in the above formula is considered to be on both sides of the peak point of the gaussian component, and the potential superimposed gaussian component should be distributed outside the interval, so as to avoid detecting false superimposed gaussian component. Therefore, the superimposed gaussian component is sought on both sides of the peak point of the gaussian component and outside the interval of the above equation. It is theorized that the effective superimposed gaussian component should be located at a distance of at least 1 pulse width (FWHM) from the adjacent gaussian component, and 2 · D may be approximately equal to FWHM. Therefore, for the superimposed gaussian components present to the left of the apparent peak point, the inflection points should be present at positions to the left of x (i) -D, and for the superimposed gaussian components present to the right of the apparent peak point, the inflection points should be present at positions to the right of x (i) + D.

The invention provides an improved detection scheme of a superposition Gaussian component, which judges whether a superposition Gaussian component exists or not by detecting two adjacent inflection points at two sides of the effective obvious peak point and increasing comparison of corresponding amplitude values of the two adjacent inflection points in an original echo waveform. And the superposed Gaussian components detected at two sides of the effective obvious peak point belong to the corresponding effective original echo waveform.

Step S3 specifically includes the following substeps:

and S3.1, setting an interval [ X (i) -D, X (i) + D ], wherein X (i) is a time coordinate value corresponding to the Gaussian component corresponding to the effective obvious peak point, and D is an estimated initial value of a standard deviation of the Gaussian component corresponding to the effective obvious peak point.

And S3.2, detecting a superimposed Gaussian component on the left side of the effective obvious peak point based on the set interval.

In the step S2, the effective significant peak point in the echo waveform after gaussian filtering is determined, and initial detection and parameter initial value estimation are performed on the gaussian component belonging to the effective original echo waveform corresponding to the effective significant peak point. And then, in the step, for the Gaussian component corresponding to each effective obvious peak point, detecting the superposed Gaussian component in the monotonous area on the left side of the effective obvious peak point. The effective original echo waveform data y (i) (where i is 1,2, …, N) before filtering, the filtered echo waveform data y (i) (where i is 1,2, …, N) and the corresponding time coordinate data x (i) (where i is 1,2, …, N), and the gaussian component and the corresponding parameter corresponding to the preliminarily obtained effective significant peak point are used as input data.

And (3) detecting inflection points in the sampling points of the monotonous area on the left side of the effective obvious peak point and outside the interval of [ X (i) -D, X (i) + D ] according to an inflection point detection method combining the formulas (13) and (14), wherein if two first adjacent inflection points are detected, the corresponding sequence numbers of the two adjacent inflection points in X time coordinate data are respectively recorded as j1 and k1, and the corresponding time coordinate values are respectively X (j1) and X (k 1).

In order to improve the accuracy of the superimposed gaussian component, it is specified that the distance between X (j1) and X (k1) is greater than half the full width at half maximum of the transmitted waveform, and the amplitude value y (j1) of the left-hand inflection point of the two adjacent inflection points in the original echo waveform should be smaller than the amplitude value y (k1) of the right-hand inflection point in the original echo waveform.

Specifically, the superimposed gaussian component to the left of the significant peak point is detected according to the following formula:

and y (j1 < y (k1) (25)

The amplitude value in the effective original echo waveform before filtering is selected to be added for judgment because the amplitude value of the sampling point of the effective original echo waveform after filtering is compared with the amplitude value of the echo waveform after filtering, the amplitude value of the sampling point is not lost due to filtering processing and can be used as a reference for auxiliary judgment, so that whether superposition Gaussian components exist on the left side of the effective obvious peak point can be judged more effectively. If the condition of the above formula is satisfied, a superimposed gaussian component is detected in two adjacent inflection point regions with corresponding sequence numbers j1 and k1 in the X time coordinate data on the left side of the effective significant peak point, and the time coordinate value of the inflection point closer to the effective significant peak point located on the right side of the two inflection points is taken as the center position of the superimposed gaussian component.

And S3.3, detecting the superposition Gaussian component on the right side of the effective obvious peak point based on the set interval.

After the detection of the superimposed gaussian component on the left side of the significant peak point in the above step S2 is completed, the superimposed gaussian component is detected in the monotonous region on the right side of the significant peak point in a manner similar to that in the previous step S2.2.1. The effective original echo waveform data y (i) (where i is 1,2, …, N) before filtering, the filtered echo waveform data y (i) (where i is 1,2, …, N) and the corresponding time coordinate data x (i) (where i is 1,2, …, N), and the gaussian component and the corresponding parameter corresponding to the preliminarily obtained effective significant peak point are used as input data.

Specifically, in the sampling points in the monotonous region on the right side of the significant peak point, and outside the interval of [ X (i) -D, X (i) + D ], inflection point detection is performed according to the inflection point detection method combining equations (13) and (14), if the first two adjacent inflection points occur, and the corresponding positions of the two adjacent inflection points in the X time coordinate data are respectively recorded as j2 and k2, the corresponding time coordinate values are respectively X (j2) and X (k 2). If the following conditions are satisfied:

and y (j2 > y (k2) (26)

Wherein j2 and k2 are respectively corresponding positions of two adjacent inflection points detected outside the interval set in the monotonic interval on the right side of the effective significant peak point in X time coordinate data, X (j2) and X (k2) are respectively time coordinate values of two adjacent inflection points detected outside the interval set in the monotonic interval on the right side of the effective significant peak point, y (j2) and y (j2) are respectively amplitude values of the two adjacent inflection points detected in the monotonic area on the right side of the effective significant peak point in the effective original echo waveform before filtering, and FWHM represents the full width at half maximum of the transmitted waveform. The amplitude value in the effective original echo waveform before filtering is selected to be added for judgment because the amplitude value of the effective original echo waveform after filtering is compared, the amplitude value of the effective original echo waveform is not lost due to filtering processing, and the judgment can be assisted, so that whether superposition Gaussian components exist on the right side of the effective obvious peak point or not can be judged more effectively. If the condition of the above formula is satisfied, a superimposed gaussian component is found in two adjacent inflection point regions with corresponding sequence numbers j2 and k2 in the X time coordinate data on the right side of the effective significant peak point, and the time coordinate value of the inflection point closer to the effective significant peak point located on the left side of the two inflection points is taken as the center position of the superimposed gaussian component.

And S3.4, estimating an initial value of the parameter of the superposition Gaussian component.

And after the detection of the superposed Gaussian components on the left side and the right side of all the effective obvious peak points is finished, estimating the initial value of the parameter of the superposed Gaussian components.

Step S3.4 specifically includes the following substeps:

and S3.4.1, estimating the initial value of the parameter of the Gaussian component superimposed on the left side of the effective obvious peak point.

For the gaussian component superimposed on the left side of the significant peak point, two adjacent inflection points are detected, and the time coordinate values of the two adjacent inflection points are X (j1) and X (k1), then:

estimating the amplitude a of the left-side superimposed gaussian component as an amplitude value of an inflection point located at a right position and closer to the effective significant peak point, of first two adjacent inflection points located outside the set interval on the left side of the effective significant peak point, that is, a ═ Y (k1), where k1 is a sequence number in the X-time coordinate data of the inflection point located at the right position and closer to the effective significant peak point, of the first two adjacent inflection points, and Y (k1) is an amplitude value of the inflection point located at the right position and closer to the effective significant peak point, of the first two adjacent inflection points;

estimating a central time coordinate μ value of the superimposed gaussian component to be a value of X (k1), that is, μ ═ X (k1), where k1 is a sequence number in X time coordinate data of an inflection point located closer to the effective significant peak point in a right position of first two adjacent inflection points located outside the set section on the left side of the effective significant peak point, and X (k1) is a time coordinate value corresponding to an inflection point located closer to the effective significant peak point in the right position of the first two adjacent inflection points;

estimating the standard deviation σ of the superimposed gaussian component as half of the distance between two first adjacent inflection points outside the set interval on the left side of the effective significant peak point, i.e., | X (j1) -X (k1) |/2, where j1 and k1 are the sequence numbers of the two first adjacent inflection points in the X time coordinate data, and X (j1) and X (k1) are the time coordinate values corresponding to the inflection points located at the left position and the right position in the two first adjacent inflection points, respectively.

And S3.4.2, estimating the initial value of the parameter of the Gaussian component superimposed on the right side of the effective obvious peak point.

For the superimposed gaussian components to the right of the significant peak, two adjacent inflection points are detected, and the time coordinate values of the two adjacent inflection points are X (j2) and X (k2), then:

estimating the amplitude a of the right-side superimposed gaussian component as an amplitude value of an inflection point located closer to the effective significant peak point at the left position out of the first two adjacent inflection points located outside the set interval on the right side of the significant peak point, that is, a ═ Y (j2), where j2 is a serial number in the X-time coordinate data of the inflection point located closer to the effective significant peak point at the left position out of the first two adjacent inflection points, and Y (j2) is an amplitude value of the inflection point located closer to the effective significant peak point at the left position out of the first two adjacent inflection points;

estimating a central time coordinate μ of the superimposed gaussian component to be a value of X (j2), that is, μ ═ X (j2), where j2 is a sequence number in X time coordinate data of an inflection point located closer to the effective significant peak point at a left position among two first adjacent inflection points located outside the set interval on the right side of the effective significant peak point, and X (j2) is a time coordinate value corresponding to an inflection point located closer to the effective significant peak point at a left position among the two first adjacent inflection points;

the standard deviation σ of the superimposed gaussian component is estimated to be half of the distance between two first adjacent inflection points outside the set interval on the right side of the effective significant peak point, i.e., | X (j2) -X (k2) |/2, where j2 and k2 are the serial numbers of the first adjacent inflection points in the X time coordinate data, and X (j2) and X (k2) are the corresponding time coordinate values of the first adjacent inflection points.

And S3.5, calculating the number of the Gaussian components detected in the filtered echo waveform.

Specifically, the sum of the number of the detected effective significant peak points and the number of superimposed gaussian components on both sides of the effective significant peak points is calculated, so as to obtain the number of gaussian components contained in the filtered echo waveform.

Table 1 shows the initial value estimation results of the gaussian component parameters such as amplitude, center position, standard deviation, etc. of the selected original echo waveform, and as seen from table 1, through the above step S2, for one selected original echo waveform case, 9 initial gaussian components are detected, and the parameters of the gaussian components are shown in table 1.

Table 1 estimation results of initial parameter values of gaussian components included in the corresponding original echo waveform in fig. 1

After the above steps are completed, it is necessary to optimally estimate the initial values of the parameters of all gaussian components detected in the filtered echo waveform obtained in steps S2 and S3 in step S4. For the superposed echo waveform, the complex multi-peak echo waveform and the weak echo waveform, it is difficult to obtain all gaussian components that should be included in the original echo waveform through the above steps, that is, the gaussian components corresponding to the ground object target are leaked. Therefore, in the optimized estimation of the gaussian components (including the gaussian component corresponding to the effective peak point and the superimposed gaussian component) detected in the filtered echo waveform of step S4, it is necessary to find the gaussian component of the feature target that cannot be detected in steps S2 and S3, that is, the missing gaussian component, by waveform fitting the residual threshold value.

And step S4, updating the total number of gaussian components detected in the filtered echo waveform and performing an optimized estimation on the parameters of the gaussian components.

After the steps S2 and S3 are completed, obtaining gaussian components and superimposed gaussian components corresponding to all the effective distinct peak points, calculating the total number of the gaussian components (including the gaussian components corresponding to the effective peak points and the superimposed gaussian components) detected in the steps S2 and S3, and estimating a of each gaussian component_i、μ_i、σ_i(where i ═ 1,2, …, k) and the like, and the above results, as well as valid raw echo waveform data y (i) (where i ═ 1,2, …, N) and corresponding time coordinate data x (i) (where i ═ 1,2, …, N) are used as inputs in this step.

The detection results obtained in steps S2 and S3 represent initial results obtained in the effective original echo waveform, and therefore, according to the results obtained in steps S2 and S3, the final result of the number of gaussian components contained in the effective original echo waveform and the optimization result of the parameters of the gaussian components are obtained by updating the total number of gaussian components (including the gaussian component corresponding to the effective peak point and the superimposed gaussian component mentioned above) and optimizing the parameters of the gaussian components (including the gaussian component corresponding to the effective peak point and the superimposed gaussian component mentioned above) and approximating the fitted echo waveform to the effective original echo waveform in this step.

The improvement of the invention is that after the initial judgment of the ground object target Gaussian component and the detection of the superposition Gaussian component and the estimation of the parameter initial value of the corresponding Gaussian component are carried out in the filtered echo waveform, all the parameter initial values of the Gaussian components are combined into a whole to be used as input data, in the iterative optimization process of the parameters of the Gaussian components (including the Gaussian component corresponding to the effective peak point and the superposition Gaussian component), the threshold value of a waveform fitting residual error RMSE is set to continuously search the Gaussian components which are contained in the effective original echo waveform and are missed due to the undetected residual errors in the steps S2 and S3, including the Gaussian components in the complex multi-peak echo waveform, and the constraint condition is set to reject the ineffective Gaussian components until the termination condition of the echo waveform decomposition is met, and finally the Gaussian components detected in the filtered echo waveform (including the Gaussian components corresponding to the effective peak point) are finished Component and the superimposed gaussian component described above and the gaussian component of the missed surface feature target) to obtain a final result of the total number of gaussian components included in the effective original echo waveform and an optimized result of the parameters of the gaussian components. Step S4 specifically includes the following substeps:

In the step, an LM parameter iterative optimization algorithm is required, wherein the LM algorithm is provided by Levenberg and Marquardt, and the method introduces damping on the basis of a Gauss-Newton method, so that a damping Gauss-Newton method is formed. The LM algorithm firstly establishes an error function f (x) between a fitting function and an observed value, wherein x is an object parameter needing to be solved, and then establishes an object function F (x):

and minimizing the result value of the objective function through nonlinear iterative optimization, thereby obtaining the optimal result of the objective parameter.

According to the principle of the LM algorithm, there are:

(J^TJ+μI)h_lm＝-g (27)

h_lm＝-(J^TJ+μI)^-1g (28)

wherein J is the Jacobian matrix of the error function f (x), g ═ J^Tf (x), mu is damping coefficient (mu > 0), h_lmFor the update amount of the target parameter x to be solved, I is AND matrix J^TJ is the same number of rows and columns of the identity matrix.

The specific process of LM algorithm iteration is as follows:

(1) setting an initial value of mu, setting a coefficient beta, setting a maximum iteration number m, and estimating an initial value x of a target parameter x₀Calculating Jacobian matrix J and matrix J^TJ. The number of iterations that have been performed is set to j, and the initial j is 0.

(2) According to the initial value x of the target parameter x₀Is represented by x_jCalculating an error function f (x)_j) To calculate the objective function F (x)_j) Value e of_jAnd calculating g ═ J^TThe value of f.

(3) From the above equations (27) and (28), x is obtained_jIncrement h of_lmLet x_j+1＝x_j+h_lmX is to be_j+1Substituting the error function f (x)_j+1) Calculating an objective function F (x)_j+1) Value e of_j+1。

(4) If e_j+1＜e_jThen the value of the target parameter x, i.e. x, is updated_j+1＝x_j+h_lmIf mu is equal to mu/beta; if e_j+1≥e_jThen the value of the target parameter x, i.e. x, is not updated_j+1＝x_jWhen μ is μ β, μ is set to μ β.

According to the iteration process, if the process (1) - (4) is iteration 1 time, the executed iteration time is updated to j ═ j +1, and when the iteration time j of the LM algorithm reaches the maximum iteration time limit, the whole iteration process is ended, and finally the optimal result x of the target parameter x to be obtained is obtained_m。

The above is the iterative optimization principle of the LM algorithm, and the following is the specific method steps. Step S4.1 specifically includes the following substeps:

According to the iteration principle of the gaussian component parameters, because the original echo waveform is an actual observed waveform to be subjected to approximation by the gaussian mixture model, in order to avoid the influence of background noise on the parameter iteration of the gaussian component, the amplitude value of a sampling point of which the amplitude value is greater than the background noise threshold value in the effective original echo waveform is taken as an observed value, and the amplitude value participates in the following iterative optimization process for fitting the effective original echo waveform, wherein n observed values meeting the condition in the effective original echo waveform are provided, and the corresponding serial number of the observed values in the time coordinate data X is p_iWhere i is 1,2, …, n.

Establishing an error function between a fitting function and the observed value, wherein a specific formula is as follows:

setting the target parameter x to be obtained as (A)₁,μ₁,σ₁,A₂,μ₂,σ₂,…,A_k,μ_k,σ_k) Then, then

Where i is 1,2, …, n (30)

Wherein the content of the first and second substances,

taking the time coordinate data of the sampling points with the amplitude values larger than the background noise threshold value in the effective original echo waveform as a fitting function of an independent variable f_i(x) As an error function, representing the p-th of the effective original echo waveform_iThe error between the amplitude value of the sampling point and the fitted amplitude value, where x is the target parameter to be found, i.e. x ═ a₁,μ₁,σ₁,A₂,μ₂,σ₂,…,A_k,μ_k,σ_k) (ii) a n is the number of sampling points selected from the effective original echo waveform as observed values, y (p)_i) For the p-th in the effective original echo waveform_iAmplitude value of each sample point, p_iSelecting sampling points from the effective original echo waveform as sequence numbers corresponding to observed values, wherein i is 1,2, …, n, X (p)_i) For the p-th in the effective original echo waveform_iTime coordinate value of each sampling point, k is the number of Gaussian components, A_j、μ_j、σ_jRespectively representing the amplitude, center position, standard deviation, mu, of the jth Gaussian component_nIs the average of the background noise.

The p-th wave form of the effective original echo wave form_iTime coordinate value X (p) of each sampling point_i) (where i ═ 1,2, …, n), the p-th echo in the effective original echo waveform_iAmplitude value y (p) of each sampling point_i) And the average value mu of the background noise is substituted into the above equation (30) to form an error function containing the target parameter to be obtained.

According to the fitting function and the error function, an objective function F (x) is established, and the specific formula is as follows:

wherein F (x) is an objective function, f_i(x) Representing selected p-th of the effective original echo waveform participating in fitting as an error function_iError between amplitude values of the sampling points and the fitted amplitude values, where x ═ a₁,μ₁,σ₁,A₂,μ₂,σ₂,…,A_k,μ_k,σ_k) And n is the number of sampling points selected from the effective original echo waveform as observed values.

The meaning of equation (31) above is to calculate half the sum of the squared errors between the fitted echo waveform and the observed values of the valid raw echo waveform. The objective of using the LM algorithm is to minimize the value of the objective function f (x).

K gaussian components and a of each gaussian component obtained according to the above step S2_i、μ_i、σ_i(where i is 1,2, …, k), and so on, in conjunction with the fitting function containing the target parameter to be found established in step s4.1.1 above

Error function f_i(x) And an objective function F (X), time coordinate data X of sampling points in the effective original echo waveform and n observation values (the corresponding serial number in the time coordinate data X is p) selected in the step S4.1.1 and participating in the fitting process_iWhere i ═ 1,2, …, n). These results are used as input in this step.

The method specifically comprises the following substeps:

step S4.1.2.1, setting the initial values of the damping coefficient mu and coefficient beta, the maximum iteration times and the initial values of the parameters of k Gaussian components to be solved according to the LM algorithm principle and the iteration process.

According to the LM algorithm principle and the specific iteration process, before the LM algorithm is adopted, an initial value of mu is set, generally, the value of mu is 0.01, an initial value of a coefficient beta is set, generally, the value of beta is 10, the maximum iteration number m is set, and generally, m is 30. Let r be the number of iterations performed, then r is 0 before the LM algorithm iteration process begins.

At step S4.1.2.2, a Jacobian matrix of the fitting function and the error function is calculated.

According to the iterative optimization principle of the LM algorithm and the above equations (29) and (30), the error function f containing the target parameter to be solved is established through the above step S4.1.1_i(x) (where i is 1,2, …, n), calculating the fitting function and the error function f_i(x) Jacobian matrix J (x)₀)，J(x₀) The calculation method of (c) is as follows:

wherein the content of the first and second substances,

is the initial estimated value of the parameter x of k Gaussian components to be solved, k is the initial number of the Gaussian components, f₁(x)，f₂(x)，…，f_n(x) Is an error function representing the error function between the amplitude values of the N sample points selected from the valid raw echo waveform (also called observed waveform) participating in the fitting process and the fitted amplitude values, where i ═ 1,2, …, N,

wherein the content of the first and second substances,

as a function of the error f₁(x) For amplitude parameter A in target parameter₁The partial derivative is calculated to obtain a value, where x ═ x₀，

As a function of the error f₁(x) For the central position (time coordinate) parameter mu in the target parameter₁The partial derivative is calculated to obtain a value, where x ═ x₀，

As a function of the error f₁(x) For standard deviation parameter sigma in target parameter₁The partial derivative is calculated to obtain a value, where x ═ x₀Similarly, other partial derivatives have similar meanings.

The above-mentioned initial estimation value x of parameter x of k Gaussian components to be solved₀Detecting a of k gaussian components contained in the effective original echo waveform in the step S2_i、μ_i、σ_i(where i is 1,2, …, k), or the like, i.e. estimating the initial values of the parameters, i.e. the initial values of the parametersWhere k is the initial number of gaussian components obtained in step S2,

is amplitude A_iIs determined by the first time estimate of (a),as a central position (time coordinate) mu_iIs determined by the first time estimate of (a),

is the standard deviation sigma_iIs estimated.

Initial estimated value x₀Substituting into the above equation (32), calculating Jacobian matrix J (x) of error function between the fitting function and the observed values of the selected sampling points in the effective original echo waveform₀). Then, J (x) is calculated₀)^TJ(x₀) Wherein the superscript T represents the Jacobian matrix J (x)₀) And (5) calculating transposition.

Step S4.1.2.3 calculates the values of the error function and the target function between the fitted function established at step s4.1.1 and the observed values of the selected sample points in the effective original echo waveform.

First estimate of x

And the time coordinate data X of the sampling points in the effective original echo waveform and the n observation values (the corresponding serial number in the time coordinate data X is p) selected in the step S4.1.1 and participating in the fitting process_iWhere i is 1,2, …, n) is substituted into the above equation (30), an error function f between the fitting function and the observed values of the selected sample points in the effective original echo waveform is calculated_i(x₀) Wherein i is 1,2, …, n. Then, the obtained error function f_i(x) Is substituted into the above equation (31), the objective function F (x) in the above step s4.1.1 is calculated₀) A value of (i), i.e

Is provided with

The above calculation formula can be uniformly expressed as:

wherein F (x)_r) Representing A of k Gaussian components to be solved_i、μ_i、σ_i(where i is 1,2, …, k) and the like at the beginning of the optimization of the r-th iteration of the LM algorithm

Corresponding objective function value, using e_rIs represented by e_r＝F(x_r)。

At step S4.1.2.4, the increments of the parameter estimates for the k Gaussian components to be solved are calculated.

According to the principle of the above formula (28), A of k Gaussian components to be solved is obtained_i、μ_i、σ_i(where i is 1,2, …, k) and so on_lmI.e. h_lm＝(ΔA₁,Δμ₁,Δσ₁,ΔA₂,Δμ₂,Δσ₂,…,ΔA_k,Δμ_k,Δσ_k) The specific calculation formula is as follows:

h_lm＝-(J(x_r)^TJ(x_r)+μI)^-1*J(x_r)^T*[f₁(x_r)，f₂(x_r)，…，f_n(x_r)]^T(34)

wherein h is_lmFor A of k Gaussian components to be solved_i、μ_i、σ_i(where i is 1,2, …, k), J (x)_r) For obtaining A of k Gaussian components to be solved at the beginning of the r-th iterative optimization of the LM algorithm_i、μ_i、σ_i(where i is 1,2, …K) estimated value x of the same parameter_rCorresponding Jacobian matrix, calculation mode and primary estimation value x₀The corresponding Jacobian matrices are similar, except for the parameter estimates.

f₁(x_r)，f₂(x_r)，…，f_n(x_r) And representing the error between the amplitude value of the selected sampling point participating in fitting in the effective original echo waveform and the fitting amplitude value by using an error function. Let x_r+1＝x_r+h_lmX is to be_r+1Substituting into the above equation (30) to calculate the error function f_i(x_r+1) (where i is 1,2, …, n), after which the objective function F (x) is calculated_r+1) A value of (i), i.e

Let e_r+1＝F(x_r+1)。

Step S4.1.2.5 determines whether to update A for k Gaussian components to be solved_i、μ_i、σ_i(where i is 1,2, …, k) and the like.

If e_r+1＜e_rThen the estimated value of the target parameter x, i.e. x, is updated_r+1＝x_r+h_lmAfter being updated

According to the iteration principle of the LM algorithm, mu is equal to mu/beta; if e_r+1≥e_rThen the target parameter A is not updated_i、μ_i、σ_i(where i is 1,2, …, k), or the like, i.e. x_r+1＝x_r，

According to the iterative principle of the LM algorithm, μ ═ μ · β.

According to the above iterative process, if the process from step S4.1.2.1 to step S4.1.2.5 is that the LM algorithm iterates 1 time, the number of iterations executed is updated to be r ═ r +1, and when the number of iterations executed by the LM algorithm r is equal to the maximum number of iterations m, generally, m is 30, and the whole iterative process is ended.According to the above iterative process of LM algorithm, obtain the optimized result of target parameter x, namely

As optimized amplitude values

Center position (time coordinate) value

And standard deviation value

Obtaining A of k Gaussian components_i、μ_i、σ_i(where i is 1,2, …, k) and obtaining the approximate waveform of the effective original echo waveform

Can be expressed as:

wherein, mu_nIs the average value of background noise, k is the number of Gaussian components,

andrespectively, the amplitude, the central position (i.e. time coordinate) and the standard deviation parameter of the ith Gaussian component after the iterative optimization of the LM algorithm.

Step S4.2, combining threshold th of fitting residual error RMSE_RMSEAnd performing further iterative optimization on parameter values of the detected Gaussian components in the filtered echo waveform.

Through the step S4.1, the optimized amplitude values of k Gaussian components are obtained

Center position (i.e. time coordinate) valueStandard deviation valueAnd fitted amplitude values of a fitted waveform of the effective raw echo waveformWhere i is 1,2, …, N, these results are used as input in this step. Since there is a possibility that the gaussian component of the ground object target is not detected and is leaked out after the initial determination of the gaussian component detected in the filtered echo waveform and the detection of the superimposed gaussian component detected in the above steps S2 and S3 are performed, it is necessary to set the threshold th of the fitting residual RMSE by this step after step S4.1_RMSEAnd judging whether the condition of the Gaussian component of the ground object target is leaked, if the condition of the Gaussian component of the ground object target is leaked, updating the total number and the parameter values of the Gaussian components detected in the filtered echo waveform, and further performing iterative optimization on the parameter values of the Gaussian components detected in the filtered echo waveform.

Step S4.2 specifically includes the following substeps:

After the LM algorithm iterative optimization process in step S4.1, calculating a residual RMSE between the amplitude values of the sampling points in the effective original echo waveform and the fitting values of the amplitude values of the sampling points obtained after the LM algorithm iterative optimization, the calculation is as follows:

wherein n represents the number of sampling points in the effective original echo waveform data participating in fitting, only the sampling points which are more than or equal to the background noise threshold value in the effective original echo waveform are selected for fitting approximation, y (i) is the amplitude value of the ith sampling point in the effective original echo waveform,and obtaining the fitting amplitude value of the ith effective original echo waveform sampling point after LM algorithm iterative optimization.

After the steps are carried out, setting a fitting residual error threshold th_RMSEFor determining the gaussian component missing in step S4.1. Specifically, for a satellite-borne laser altimetry system, the signal-to-noise ratio of a laser echo signal is low and is susceptible to noise, and a fitting residual error threshold th is set_RMSEThe threshold value is expressed by the following formula:

th_RMSE＝4.5·δ_n(37)

wherein, delta_nRepresenting the standard deviation of the background noise.

According to the fitting residual error RMSE obtained in the step S4.1 and the set fitting residual error RMSE threshold th_RMSEAnd determining whether effective Gaussian components are missed or not, updating the number of the Gaussian components, and optimizing the Gaussian components again. The method comprises the following specific steps:

if residual RMSE > th_RMSEThen, finding a new Gaussian component and calculating

Finding out the position of the maximum value of the error in the effective original echo waveform, marking as ii, adding a new Gaussian component at the position, and respectively setting the amplitude, the central position and the standard deviation of the Gaussian component as A_ii，μ_ii，σ_iiFor the 3 piecesParameters are assigned to initial values:

A_ii＝y(ii)，μ_ii＝X(ii)，

wherein A is_iiFor the initial amplitude value of the found new Gaussian component, y (ii) for the amplitude value corresponding to the ii-th position of the found new Gaussian component in the sampling point of the effective original echo waveform, μ_iiAn initial value of the center position of the found new Gaussian component, X (ii) a time coordinate, σ, corresponding to the ii-th position of the found new Gaussian component in the sampling point of the effective original echo waveform_iiThe initial value of the standard deviation of the new gaussian component found.

And updating the number of the detected Gaussian components in the filtered echo waveform, namely increasing the total number of the Gaussian components by 1, namely k is k + 1. The principle of estimating the initial value of the found new gaussian component is that the initial value of the amplitude of the new gaussian component can be estimated by using the amplitude value corresponding to the position of the new gaussian component in the effective original echo waveform, the initial value of the central position of the new gaussian component can be estimated by using the time coordinate corresponding to the position of the new gaussian component in the effective original echo waveform, and the initial value of the standard deviation of the new gaussian component can be obtained by dividing the full width at half maximum of the transmitting waveform of the satellite-borne laser radar by the full width at half maximum of the transmitting waveform of the satellite-borne laser radar

The estimation is performed because the numerical relationship between the standard deviation σ and FWHM of the theoretical gaussian component is

The foregoing steps S4.2.1 and S4.2.2 determine the number of the missing gaussian components, obtain initial values of parameters of the gaussian components, update the number of the gaussian components detected in the filtered echo waveform, and use the updated number of the gaussian components detected in the filtered echo waveform and the updated parameter values of the gaussian components as input for the present step. This gaussian component is added to the gaussian mixture model of equation (23), and the iterative optimization process of the LM algorithm of step S4.1.2 described above is performed again.

Through the steps S4.2.1 and S4.2.2, the number of gaussian components and parameters of gaussian components after iterative optimization are obtained. Then, the following judgment is made:

(1) if RMSE > th_RMSEAnd adding a new Gaussian component and updating the total number and parameters of the Gaussian components. Then, the parameter iteration process of the gaussian components in sub-step S4.1.2 in step S4.1 is repeated, and after that, the above steps S4.2.1, S4.2.2 and S4.2.3 are continuously performed, and additional constraint conditions are added in the process of optimizing the gaussian component parameters, so as to remove invalid gaussian components, and update the total number of gaussian components and the parameters of the gaussian components. Until RMSE is less than or equal to th_RMSEAnd finally, stopping the whole process of the parameter iterative optimization of all Gaussian components contained in the effective original echo waveform.

(2) If RMSE is less than or equal to th_RMSEThen, the following additional constraint conditions are added for eliminating invalid gaussian components and updating the number of gaussian components and parameters of the gaussian components. Then, the parameter iteration process of the gaussian component in sub-step S4.1.2 in step S4.1 is performed again until all invalid gaussian components are removed, and the parameter iteration process of the entire gaussian component is terminated.

Additional constraints:

in the process of updating the number of gaussian components and iteratively optimizing the parameters of the gaussian components in step S4.2, additional constraint conditions are set to remove invalid gaussian components, and the conditions are as follows:

(1) minimum distance between gaussian components > minimum distance threshold in ns;

(2) minimum amplitude of gaussian component > noise threshold th;

(3) minimum pulse width of gaussian component > full width at half maximum (FWHM) of emission waveform;

(4) the number n of Gaussian components is less than or equal to the maximum number threshold.

The reason for setting the above condition is that the minimum distance threshold of the condition (1) is set because of influence of noise, and if the distance between gaussian components is set too small, there is a case where a signal of an unreal feature is erroneously determined as a valid gaussian component, and therefore the condition of (1) needs to be set to improve the resolution accuracy of the gaussian component, and preferably, the minimum distance threshold is 10ns, corresponding to a distance of 1.5 meters. The condition (2) is set to avoid the influence of background noise on the detection of the gaussian component. The setting condition (3) is to avoid the influence of invalid gaussian components, and the setting condition (4) is a limiting condition that is provided in consideration of the characteristics of the vertical distribution of the ground features on the satellite laser transmission path and the limitations of the characteristics of the echo waveform signal itself, and the threshold value of the maximum number of gaussian components is usually 6. If the Gaussian component does not meet any condition, the echo waveform is judged to be invalid and eliminated, and then the fitting of the echo waveform and the optimization of the parameters are carried out again. Through the further iterative optimization step of the parameter values of the Gaussian components, the method can also finish the detection and decomposition of the Gaussian components of the ground object target in the area with a complex ground object structure, such as a forest vegetation area and the like. Through the steps S1, S2, S3 and S4, the invention finally completes the decomposition of the satellite-borne large-spot laser radar full waveform data, and Gaussian components of each ground object target contained in the effective original echo waveform are all decomposed from the effective original echo waveform, so that the amplitude A of each Gaussian component in each effective original echo waveform is extracted and obtained_iCentral position μ_iStandard deviation σ_i(where i ═ 1,2, …, k) and the like, and the number of gaussian components.

Fig. 2(e) is a schematic diagram of a waveform decomposition result corresponding to the decompressed original echo waveform according to the embodiment of the present invention, in which the final original echo waveform decomposition result obtained through step S3 is shown. Through the optimized estimation of Gaussian component parameters in the echo waveform in the step, 6 Gaussian components are obtained through decomposition, and the number of the Gaussian components and parameters such as the amplitude, the central position and the standard deviation of each Gaussian component are obtained. The decomposition result values of the original echo waveform corresponding to fig. 2(e) are shown in table 2.

Table 2 decomposition results of the original echo waveform corresponding to fig. 2

Taking GLAS data as an example, fig. 3(a) - (c) show typical examples of decomposition results of 3 kinds of original echo waveforms which are difficult to process, and show 3 kinds of echo waveforms which are difficult to process, such as a superposition echo waveform, a weak echo waveform, and a complex multi-peak echo waveform. Fig. 3(a) is one of typical examples of waveform decomposition results corresponding to an original echo waveform after decompression according to the embodiment of the present invention, which represents decomposition results of an original echo waveform that is superimposed. Fig. 3(b) shows a second typical example of the waveform decomposition result corresponding to the decompressed original echo waveform according to the embodiment of the present invention, which represents the decomposition result of the original echo waveform containing weak echo components. Fig. 3(c) shows a third exemplary waveform decomposition result corresponding to the decompressed original echo waveform according to the embodiment of the present invention, which represents the decomposition result of the original echo waveform with multiple peaks in a complex pattern. Therefore, the echo waveforms are effectively decomposed by the method, and the Gaussian component and the parameter thereof contained in each original echo waveform are obtained.

In order to verify the full waveform data decomposition method of the satellite-borne large-spot laser radar, echo waveform data of the GLAS are used as input data, a Beijing plain area and a great-fun and safe-land area of Heilongjiang province are respectively selected as test areas, and a large number of full waveform data decomposition tests are carried out. Tables 3 and 4 show statistics of the number of gaussian components of the waveform decomposition results in the beijing plain region and the black longjiang great khingan region, respectively. The NSIDC (national Snow and IceCorter) algorithm represents the official waveform data decomposition method of GLAS.

Compared with the NSIDC method, the number of original echo waveforms with the number of Gaussian components of 1 is reduced by 444, and the number of original echo waveforms with the number of Gaussian components of 2, 3,4, 5 and 6 is increased in the Beijing plain area. The above results show that: the waveform decomposition method provided by the invention can obtain more surface feature waveform components in plain areas. In the great-Xingan mountain area of Heilongjiang, the number of original echo waveforms with the number of Gaussian components of 1 is reduced by 309, the number of original echo waveforms with the number of Gaussian components of more than 4 is reduced, and the number of original echo waveforms with the number of Gaussian components of 2 and 3 is increased. The above results show that: the waveform decomposition method provided by the invention can represent the vertical distribution condition of the ground objects in the laser spot by less Gaussian components in the forest vegetation area while ensuring the fitting accuracy of Gaussian component decomposition.

TABLE 3 Gaussian component count of waveform decomposition results in Beijing plain area

TABLE 4 Gaussian component count for waveform decomposition results in the region of great Kyoho

FIGS. 4(a) - (b) show R of the results of waveform decomposition corresponding to the experiments²And (6) counting the histogram. FIG. 4(a) is a graph of R of a large number of waveform decomposition results according to an embodiment of the present invention²Statistical histograms are represented in the graph by the Beijing plain test area. FIG. 4(b) is a graph of R of a large number of waveform decomposition results according to an embodiment of the present invention²And (5) counting the histogram, wherein the graph represents a test area of Daxing-an mountain in Heilongjiang province.

It can be seen that the waveform fit determines the coefficient R²Mainly distributed at 0.95, it is shown that the fitting result of the waveform decomposition method of the present invention has high accuracy.

The above description is only a preferred embodiment of the present invention, and for those skilled in the art, the present invention should not be limited by the description of the present invention, which should be interpreted as a limitation.

45页详细技术资料下载

上一篇：一种医用注射器针头装配设备

下一篇：光接收模块

Full waveform data decomposition method for satellite-borne large-light-spot laser radar

相关技术

网友询问留言