Decomposition of small-footprint full waveform LiDAR data based on generalized Gaussian model and grouping LM optimization

The classical LM algorithm is a blend of the Gauss–Newton method and the gradient descent method [30]. It has been proved that gradient descent and Gauss–Newton iteration are complementary in the advantages they provide. Assuming that there is a set of observations $\,\boldsymbol{X}\left(\boldsymbol{X}\in {{R}^{n}}\right)$. A mathematical model $\,f\left(\,\boldsymbol{p}\right)$, where $\boldsymbol{p}\in {{R}^{n}}~$ is the parameter vector whose entries need to be estimated, is used to fit the data. A common strategy in optimization is to set an initial vector of p, and then update it by an iteration process which minimizes $\,\boldsymbol{\varepsilon }^{{T}}\boldsymbol{\varepsilon }\left(\boldsymbol{\varepsilon }=\boldsymbol{X}-\tilde{{\boldsymbol{X}}}\,\right)$, where $\tilde{{\boldsymbol{X}}}\,\in \boldsymbol{R}^{{\text{n}}}$ is the vector whose entries calculated from $f\left(\boldsymbol{p}\right)$ with the currently updated p. In optimization jargon we are solving the following minimization problem:

$$\begin{eqnarray}&&\boldsymbol{p}_{{\text{opt}}}=\arg \underset{\boldsymbol{p}}{{\min}}\,||\boldsymbol{X}-f\left(\,\boldsymbol{p}\right)||\end{eqnarray} \tag{ 3 }$$

LM is one of the algorithms that achieves optimized parameters by iteration. Set an initial parameter vector $\boldsymbol{p}_{{k}}$ as the current parameter vector and substitute $\boldsymbol{p}_{{k}}$ into the equation $\tilde{{\boldsymbol{X}}}\,=f\left(\,\boldsymbol{p}\right)$ to calculate the estimated values. Then use the observations to compute the residual ${{\varepsilon}^{T}}\varepsilon$. The iteration process is performed until the residual $\boldsymbol{\varepsilon }^{{T}}\boldsymbol{\varepsilon }$ is less than the preset threshold $~\boldsymbol{\varepsilon }_{{0}}$, at which point the current $~\boldsymbol{p}_{{k}}$ is the final parameter vector. During each iteration the first order Taylor series expansion of $f\left(\,\boldsymbol{p}\right)$ near $\boldsymbol{p}_{{k}}$ is calculated:

$$\begin{eqnarray}&&f\left(\,\boldsymbol{p}_{{k}}+\boldsymbol{\delta }_{{k}}\right)=f\left(\,\boldsymbol{p}_{{k}}\right)+\boldsymbol{J}_{{k}}\boldsymbol{\delta }_{{k}}\end{eqnarray} \tag{ 4 }$$

where $\boldsymbol{J}_{{k}}=\partial f\left(~\boldsymbol{p}_{{k}}\right)/\partial \boldsymbol{p}_{{k}}$ is the Jacobian matrix evaluated at the point of $\boldsymbol{p}_{{k}}$, $\boldsymbol{\delta }_{{k~}}$ is the vector that minimizes the following:

$$\begin{eqnarray}&&\|\boldsymbol{X}-f\left(\,\boldsymbol{p}_{{k}}+\boldsymbol{\delta }_{{k}}\right)\|=\underset{{{\delta}_{k}}}{{\min}}\,\|\boldsymbol{\varepsilon }_{{k}}-\boldsymbol{J}_{{k}}\boldsymbol{\delta }_{{k}}\|\end{eqnarray} \tag{ 5 }$$

The precision of the final results and the speed of convergence are influenced by the initial parameters in the conventional LM algorithm [31]. When the initial values are closely matched to the true ones, the Jacobian matrix is non-singular; hence the LM algorithm is convergent with a high speed. Otherwise, non-numerical elements may appear in the Jacobian matrix, leading to an unexpected interruption of the iteration. To demonstrate this effect, let

$$\begin{eqnarray}&&{{f}_{j}}(x)={{A}_{j}}\exp \left(-\frac{\text{abs}{{\left(x-{{\mu}_{j}}\right)}^{\alpha _{j}^{2}}}}{2\sigma _{j}^{2}}\right)\end{eqnarray} \tag{ 6 }$$

be a Gaussian component, where $\boldsymbol{p}=\left\{{{\alpha}_{j}},{{\mu}_{j}},{{\sigma}_{j}},{{A}_{j}}\right\}$ is the parameter vector, $x$ is the sampling time, abs($x-{{\mu}_{j}}$) is the absolute value of ($x-{{\mu}_{j}}$), and ${{A}_{j}}$ is the amplitude. The partial derivative with respect to ${{\alpha}_{j}}$ is $\frac{{{A}_{j}}\centerdot {{\alpha}_{j}}\centerdot \text{abs}{{\left({{\mu}_{j}}-x\right)}^{\alpha _{j}^{2}}}\centerdot \log \left(\text{abs}\left({{\mu}_{j}}-x\right)\right)}{\sigma _{j}^{2}\centerdot \exp \left(\frac{\text{abs}{{\left({{u}_{j}}-x\right)}^{\alpha _{j}^{2}}}}{2\sigma _{j}^{2}}\right)}$, which has domain $\left\{{{\mu}_{j}}\ne x\right\}$ by definition. However, in the conventional LM optimization, $x$ can be any value including $~{{\mu}_{j}}$. Therefore an element of $\log (0)$ will occur in the Jacobian matrix when ${{\mu}_{j}}$  =  $~x$, resulting in the failure of the algorithm. This is caused by the selection of mathematical model (the phenomenon will not occur if a conventional Gaussian function is adopted; that is, if the shape index is set to $\sqrt{2}$), and by the initialization of the parameters (when the initial value of ${{\mu}_{j}}$ is $x$). There are no specific treatments for this problem in the existing literature. A common method is to modify the initial values of the parameters, but a major drawback of such an idea is equal treatment to all variables and it is not expected to be effective. Furthermore, if the initial values deviate far from the true values, then the algorithm is unlikely to be convergent. As a matter of fact, compared with the other parameters, the temporal difference $\,{{\mu}_{j}}-x$ is the most likely factor that causes the non numerical values in the Jacobian matrix when $\,{{\mu}_{j}}$ is assigned to with an integer. If dividing the involved parameters into distinct groups and then firstly optimize parameters in the group containing $\,{{\mu}_{j}}$ while treat parameters in other groups as constants, then the optimized $\,{{\mu}_{j}}$ is unlikely to be an integer. Iterate the process until convergence criteria are satisfied. This is what grouping LM works.

Bearing the abovementioned problem in mind, this paper proposes an improved LM algorithm: the grouping LM optimization. To employ this algorithm the full waveform parameters are grouped firstly by the definition domains of the elements in the Jacobian matrix. Then, the parameters in each group are sequentially optimized by the LM algorithm according to the number of constraints of each domain, from small to large, rather than optimizing all of the parameters at one time. When the parameters in one group are computed, the remaining parameters are treated as constants. Once the parameters in the first group have been optimized, then the remaining constants in the next group are optimized. Any non-numerical Jacobian entries can be effectively overcome in this way. To demonstrate, equation (6) in section 2.1 is revisited. It should be noted that a full waveform is not a mathematically continuous curve, but is rather a broken line in the 2D coordinate system spanned by the sampling time versus the amplitude. The parameter $x$ in equation (6) is an integer representing the sampling time. If the group including $~{{\mu}_{j}}$ is optimized at first, then the optimized value of ${{\mu}_{j}}$ should be a non-integer. Therefore, in the following process where other parameters such as ${{\alpha}_{j}}$, ${{\sigma}_{j}}$, etc. are to be optimized, the situation where ${{\mu}_{j}}$  =  $~x$ will not occur since an integer cannot be equal to a non-integer; hence the non-numerical elements in the Jacobian matrix are avoided. The following describes the detailed process of the grouping LM by taking a full waveform with two Gaussian components for illustration. The mathematical model for this particular situation is expressed as:

$$\begin{eqnarray}&&{{y}_{j}}\left(\,\boldsymbol{p}\right)=\underset{i=1}{\overset{k}{\sum}}\,{{A}_{{{j}_{i}}}}\exp \left(-\frac{\text{abs}{{\left(x-{{\mu}_{ji}}\right)}^{\alpha _{{{j}_{i}}}^{2}}}}{2\sigma _{{{j}_{i}}}^{2}}\right)\end{eqnarray} \tag{ 7 }$$

where $i$ denotes the number of components and p denotes the parameter vector that needs to be optimized such that $\boldsymbol{p}=\left\{{{\alpha}_{ji}},{{\mu}_{ji}},{{\sigma}_{ji}},{{A}_{ji}}\right\}_{i=1}^{k}$ (in the current case k  =  2). ${{A}_{ji}}$ is the amplitude of the full waveform, ${{\mu}_{ji}}$ denotes the location of the component peak, ${{\sigma}_{ji}}$ denotes the half width of the pulse, and ${{\alpha}_{ji}}$ is the shape parameter which can fine-tune the shape of the waveform. The objective function for optimization is:

$$\begin{eqnarray}&&\arg \underset{\boldsymbol{p}}{{\min}}\,\|x-y\left(\boldsymbol{p}\right)\|\end{eqnarray} \tag{ 8 }$$

The optimization process under the umbrella of the grouping LM can be described in detail as follows:

• (1)  
  Initialization: Determine the initial values of $\left\{\alpha _{ji}^{0},\mu _{ji}^{0},\sigma _{ji}^{0},A_{ji}^{0}\right\}_{i=1}^{2}$, where the superscript 0 denotes the number of iterations. Suppose that the raw waveform data have been pre-processed, and the peak amplitudes and locations $A_{ji}^{0}$ and $\mu _{ji}^{0}$ are determined by local maxima as follows. For a given set of pre-processed waveform data, assign the maximum amplitude and its corresponding time position plus half a sampling step to $A_{ji}^{0}$ and $\mu _{ji}^{0}$, respectively. For each peak, scan the waveform along both sides to find two points with amplitudes equal to half that of the peak. Calculate the distance along the $x$ axis between these two points, and divide the distance by two to determine $\,\sigma _{ji}^{0}$. We assume that the initial shape parameter of each echo is 1.5.
• (2)  
  Grouping: the parameters are divided into two groups: $\left\{\mu _{ji}^{0},A_{ji}^{0}\right\}_{i=1}^{2}$ and $\left\{\alpha _{ji}^{0},\sigma _{ji}^{0}\right\}_{i=1}^{2}$.
• (3)  
  Optimizing: First, the parameters $\left\{\alpha _{ji}^{0},\sigma _{ji}^{0}\right\}_{i=1}^{2}$ are treated as constants and the conventional LM algorithm is used to optimize the parameters $\left\{\mu _{ji}^{0},A_{ji}^{0}\right\}_{i=1}^{2}$ resulting in $\left\{\mu _{ji}^{\text{opt}},\,A_{ji}^{\text{opt}}\right\}_{i=1}^{2}$. Now keep $\left\{\mu _{ji}^{\text{opt}},A_{ji}^{\text{opt}}\right\}_{i=1}^{2}$ as constants to optimize $\left\{\alpha _{ji}^{0},\sigma _{ji}^{0}\right\}_{i=1}^{2}$ resulting in $\left\{\alpha _{ji}^{\text{opt}},\sigma _{ji}^{\text{opt}}\right\}_{i=1}^{2}$.
• (4)  
  Iterating: repeat step (3) with updated parameters $\left\{\mu _{ji}^{\text{opt}},A_{ji}^{\text{opt}}\right\}_{i=1}^{2}$ and $\left\{\alpha _{ji}^{\text{opt}},\sigma _{ji}^{\text{opt}}\right\}_{i=1}^{2}$, until the stop criterion mentioned in section 2.1 is satisfied.

Both simulation and real data were used to test our algorithm. 10 000 simulation data were generated in accordance to formula (2) where we assumed $\,{{A}_{i}}=1$, therefore, they are noise-free data suitable for theoretically test the proposed algorithm; 5000 simulated waveforms were generated from two Gaussian functions superpositioned, while the other 5000 waveforms were the superposition of three. The real data were acquired by the RIEGL LMS-Q560 over a mountainous area full of vegetation. The full waveform data are contained in two separate files with suffixes LWF and LGC. A ^*.LWF file records the raw full waveform data and a ^*.LGC file records attributes related to the raw waveform. The average density of the point cloud provided by the system is about 0.8 pts m⁻². The dataset contains 5487 240 waveforms. Each waveform consists of 58 samples from the emitted signal and the same number of samples from the echo.

The Gaussian decomposition of the waveform data is performed according to the steps described in section 2.2. The point cloud is generated by a set of formulae (9) [32] after parameters of the Gaussian model have been estimated by the grouping LM.

$$\begin{eqnarray}&&\left\{\begin{array}{c} {{E}_{i}}={{E}_{0}}+\text{d}E\ast \left(\text{WFOFFSET}+i\right)~ \\ {{N}_{i}}={{N}_{0}}+\text{d}N\ast \left(\text{WFOFFSET}+i\right)~ \\ {{H}_{i}}={{H}_{0}}+\text{d}H\ast \left(\text{WFOFFSET}+i\right)~ \end{array}\right. \end{eqnarray} \tag{ 9 }$$

where $i$ denotes the $i\text{th}$ sample of a waveform with $\left\{{{E}_{i}},{{N}_{i}},{{H}_{i}}\right\}$ denoting its geodetic coordinates. $\left\{{{E}_{0}},{{N}_{0}},{{H}_{0}}\right\}$ denote the geodetic coordinates of the first sample of the emission pulse. $\left\{\text{d}E,\text{d}N,\text{d}H\right\}$ denote the coordinate differentials of the waveform sampling unit. $\text{WFOFFSET}$ is the offset between the first sample of the emitted waveform and the received one.

Precision is evaluated by comparing the point clouds generated by the proposed algorithm with those provided by the system:

$$\begin{eqnarray}&&\left\{\begin{array}{c} {{\sigma}_{x}}=\sqrt{\underset{i=1}{\overset{n}{\sum}}\, \Delta x_{i}^{2}/n} \\ {{\sigma}_{y}}=\sqrt{\underset{i=1}{\overset{n}{\sum}}\, \Delta y_{i}^{2}/n} \\ {{\sigma}_{z}}=\sqrt{\underset{i=1}{\overset{n}{\sum}}\, \Delta z_{i}^{2}/n} \end{array}\right. \end{eqnarray} \tag{ 10 }$$

where $\Delta x={{X}_{1}}-{{X}_{2}}, \Delta y={{Y}_{1}}-{{Y}_{2}}, \Delta z={{Z}_{1}}-{{Z}_{2}}$; $\left\{{{X}_{1}},{{Y}_{1}},{{Z}_{1}}\right\}$ denote the 3D coordinates of the point cloud generated by the proposed method. $\left\{{{X}_{2}},{{Y}_{2}},{{Z}_{2}}\right\}\,$ denote the 3D coordinates of the point cloud provided by the system. ${{\sigma}_{x}}$, ${{\sigma}_{y}}$ and ${{\sigma}_{z\,}}$ represent the root mean square errors (RMSE) of deviation along the $X$, $Y$ and $Z$ axis respectively. $n$ is the number of data points used in the experiment.

3.3.1. Results of the simulation data.

We fitted the waveforms by using conventional LM and grouping LM and calculated the RMSE for the parameters of echo (peak) location, half width and shape index, respectively. Table 1 shows the comparison results. Generally speaking, parameters estimated by LM grouping are more precise than estimated by conventional LM in terms of their RMSEs. More specifically, the pulse width and the shape index estimated by the grouping LM are closer to the true values. Both algorithms can estimate comparable echo locations, though the grouping LM obviously outperforms the conventional LM. We can conclude that the grouping LM algorithm is superior to the conventional algorithm in all the aspects of waveform fitting.

Table 1. Comparison of the RMSEs of the parameters estimated by the conventional LM and the grouping LM.

Return echo	Algorithm	RMSE
		Echo location $\mu$	Half width $\sigma$	Shape parameter $\alpha$
The 1st return	Conventional LM	0.021	0.12	0.036
	Grouping LM	0.019	0.07	0.003
The 2nd return	Conventional LM	0.29	0.15	0.041
	Grouping LM	0.11	0.09	0.007
The 3rd return	Conventional LM	0.22	0.19	0.049
	Grouping LM	0.10	0.08	0.005

3.3.2. Results of real data.

We picked up 487 240 waveforms from the whole real dataset to generate a point cloud by the grouping LM and the conventional LM algorithm respectively in order to evaluate the precision of the point cloud generated by the proposed method. During the iteration of the conventional LM, no discrete points are generated when non- numerical elements occurred in the Jacobian matrix. This is why the grouping LM generates more points than that of the conventional LM, as listed in table 2.

Table 2. Comparing the number of points provided by the system with those generated by the conventional LM and the grouping LM.

Return echo	The system	The conventional LM	The grouping LM
The 1st return	445 010	441 727	391 777
The 2nd return	104 116	110 559	140 638
The 3rd return	28 921	29 095	52 739
The 4th return	347	466	762
The 5th return	29	41	64
The 6th return	0	3	12
Total number of generated points	578 423	581 891	585 992

As shown in the last row in table 2, both of the total number of points generated by conventional LM and the grouping LM are larger than those provided by the system, while the grouping LM has the largest number. It can be discovered by simple calculation that the number of first echoes of grouping LM is about 12% less than the other two, while the number of second echoes of grouping LM is about 26% and 21% more than that of the system and the conventional LM, respectively. More obvious variations occur with the number of third and fourth echoes: the grouping LM generates double the third and fourth echo rates of the others. These facts indicate that it is highly likely the grouping LM can generate more multiple echoes in an area of vegetation, meaning that it is the preferred decomposition algorithm in forest regions.

Deviations were calculated by subtracting the coordinates provided by the system from the corresponding values calculated by both the conventional LM and the grouping LM. As shown in figure 1, the horizontal axis indicates the serial number of the points and the vertical axis indicates the deviation. Around 12 000 points were picked up from the total dataset to plot in figure 1. The first column of figure 1 shows the deviations between the systematic point cloud and those generated by the LM algorithm, while the second column depicts the deviation between the systematic point cloud and those generated by the grouping LM algorithm.

Zoom In Zoom Out Reset image size

It is obvious from the figures that the deviation along the $X$ axis is generally within  ±0.05 m, while a small portion lies between  ±0.05 m and  ±0.4 m. Similarly, the deviation along the $Y$ axis is generally within  ±0.1 m, while a small portion lies between  ±0.1 m and  ±0.6 m. The deviation along the $Z$ axis is generally between −0.3 m and +0.2 m with the majority being less than zero. This shows that the elevation of the generated point cloud is slightly lower than that provided by the system. Though the deviations by the conventional LM and the grouping LM have the same trend, it is highly likely that deviations found with the grouping LM are slightly smaller than those of the conventional LM. This conclusion is further supported by the RMSE of deviations between the grouping LM algorithm and the conventional LM algorithm which were calculated using equation (10), and which are listed in table 3.

Table 3. Comparing the RMSE of deviations along different axis.

Algorithm	$X$ axis (m)	$Y$ axis (m)	$Z$ axis (m)
The conventional LM	0.0353	0.0726	0.2835
The grouping LM	0.0255	0.0521	0.2116

In order to further determine which algorithm can extract a digital elevation model (DEM) with a better accuracy, we collected 154 check points from the test area as ground truth data to assess the accuracies of the three datasets: point cloud provided by the LiDAR system and generated from the conventional LM and grouping LM. All the check points were collected by real time kinematic (RTK) technique, a technique commonly used for accuracy checking in real LiDAR projects [33]. The RMSEs of the elevations in the three datasets with respect to the check points were calculated, as shown in table 4; moreover, the error distributions of the three DEM datasets were demonstrated as frequency histograms, as shown in figure 2. Table 4 and figure 2 show that the grouping LM algorithm has the smallest RMSE and the most concentrated error distribution toward zero, indicating that the DEM generated by the grouping LM algorithm has a better accuracy than the others. Therefore, we conclude that the grouping LM outperforms the conventional LM in terms of DEM accuracy.

Table 4. Comparing the RMSEs of the elevations in the three datasets with respect to checkpoints.

Generation method	Number of checkpoints	RMSE (m)
The Grouping LM	154	0.28
The conventional LM	154	0.31
The system	154	0.32

Zoom In Zoom Out Reset image size

In order to evaluate the weak-pulse detectability both of the conventional LM and the grouping LM, we assessed the section profiles of the point cloud generated by the two algorithms, as shown in figure 3. The red and the green points highlight those points with low amplitudes which the system failed to detect. Obviously, more points are generated by the grouping LM algorithm than by the conventional LM algorithm. Because weak pulses usually occur at positions of potential ground points beneath forest canopies, it implies that the grouping LM can recover more ground points in vegetative areas, implying its potential usefulness for DEM generation from full waveform data. Furthermore, more points from canopies can be derived by the grouping algorithm as well, as shown in figure 3(c). These strengths would be of significant value in areas with a lower laser energy penetration rate, such as forests with heavy canopy, since in such areas the determination of weak returns from the ground is especially useful for DEM generation and for subsequent estimation of tree heights [34].

Zoom In Zoom Out Reset image size