LabVIEW development platform provides interfaces for general 3D model files like VRML, STL, ASE, and many scene control methods, including background color, light source, perspective controller, automatic projection, display mode, and many other basic settings. Besides, setting options for texture, rotation, scale, translation are also available. Through selective and flexible use of these function controls, finally we can construct and render the desired scene. Modeling methods for 3D model files in the lidar system can be various and flexible. For instance, the terrain modeling can be finished with terrain plugin in 3dsmax, which can generate most of our required positions (different latitudes and longitudes) on earth by connecting the GIS. Trees and planes can be modeled in most 3D modeling tools like 3dsmax, AutoCAD, Solidworks. The above-mentioned 3D models are exemplified and shown in Fig.  3. Our target is to integrate different 3D models, configure appropriate setting options, and generate the desired scene for lidar simulation according to our actual requirements.

Fig.  3.

Figure Option ViewDownloadNew Window

Fig.  3. (a) 3D models of aircraft. It is designed to represent the airborne laser platform, which is located above the terrain and moves with the laser scanning. (b) 3D model of terrain. It is acquired from GIS system via a certain longitude and latitude, which contains geographic distribution data around the Earth. (c) 3D model of vegetation. It is designed to construct the forested areas representing the non-ground objects above the terrain. Appropriate location, height, branch parameters, density, and many other parameters for single tree need to be determined.

Simulation environment modeling for the lidar system is so crucial that it is directly related to the proximity between the simulation and the site trial. In other words, the authenticity of the simulation environment models reflects the practicality of our system.

This subsection is devoted to building four models: the laser pulse model, atmospheric transport model, target interaction model, and the receiver model. Tomas^[11] has made an overall mathematical model and also a detailed conclusion. Here, all of the above models are summarized as follows to illustrate the application in our system, and the entire structure of environment modeling is presented in Fig.  4.

Fig.  4.

Figure Option ViewDownloadNew Window

Fig.  4. Structure of environment modeling for simulation. The input parameters for the module are classified as two levels. The 1st level input parameters are the most directly related to each model, while the 2nd level parameters are the more detailed explanations of 1st level parameters. To generate the final simulation model, some additional inputs, like elevation difference, temperature, and track angle are taken into account.

(i) Laser pulse modelling is to describe the time propagation and irradiance distribution with the mathematical model. The former is described as^{[15, 16]}

where T_1/2 is the full width at half maximum, while the latter is

(ii) The target of atmospheric transport modelling is to describe the effects of beam broadening, intensity fluctuations, and aerial attenuation on the laser beam during the process of transmission.

The rule of beam broadening is

where , Λ ₁ = λ R/π ω ², and k = 2π /λ . The fundamental parameters in the equation are structure constant of atmospheric turbulence C_n, laser wavelength λ , distance from the beams to the target R, and beam diameter at a specific distance ω .

The probability density function of the scintillation is

where , , , and γ _target = (ρ _l/r_eff)^7/3. Here, r_eff = min (r_target, ω , υ _fov/2· R) represents the radius of an effective averaging spot. The main parameters in the above equations include the geometric radius of the target r_target, the laser beam radius at target ω , the instantaneous field of view υ _fov, the detected signal value S and its time average S_av, etc.

The aerial attenuation results from the absorption and reflection effect by the particles in the air, which can be described as

where

The parameter q is a variable relating to visibility distance V_M.

(iii) When the laser beam hits a diffuse surface of the target, the laser will be reflected into many directions, which results from effects of both diffuse and specular refection. Furthermore, the produced beams are independent and randomly phased, thus making dark and bright spots, which are called speckle cells, after they have interfered. The refection expression is

where A and B are correlation factors between diffuse and specular refection.

The speckle directly affects the waveform, which is different from the conventional additive noise, but a kind of multiplicative one. The probability density function of speckle values follows:

where M= D_rec/D_sp, with D_sp = λ /ϕ , and D_rec being the radius of optical input at the receiver.

(iv) Receiver model is modelled according to the noise equivalent power (NEP) from the detector and the amplifier, which can be expressed as

Since both of the expressions are representations of additive noise, the total NEP can be described as

The main parameters involved are explained as follows: e_c is the electronic charge; B is the bandwidth; I_ds is the dark surface current; I_db is the dark bulk current; I_b is the background illumination current; I_s is the signal current; M is the detector internal gain; F is the excess noise factor; ℜ is the responsivity; k is the Boltzmann’ s constant; T is the temperature in degrees Kelvin; N is the noise factor for the electronics; R_L is the load resistor.

In order to make our system available for site trial, a simple prototype of the lidar system is designed. As shown in Fig.  5, the system consists of four main parts: a laser transceiver system, a sensor array, a data acquisition system, and a data upper-controlling and processing system.

Fig.  5.

	Figure Option ViewDownloadNew Window
	Fig.  5. Diagram of our lidar system for site trial.

The working principle can be described as follows. The receiver of the laser transceiver system emits a laser pulse at a constant frequency, while a sample of it is extracted to trigger and drive the data acquisition system after having been converted into an electrical signal. The echo signals are collected on the receiver side, and input into the sensor array. The output signals are analog electrical signals with a certain gain, and they are input into multi-channels of the data acquisition system. Multi-ADCs are controlled by an FPGA chip to work, simultaneously, after the trigger signal has come. The sampled data of multi-channels are aligned with a high accuracy, and transmitted into SSD. The upper monitor on the PC side can acquire data from SSD and convert them into decimal data, and then construct the full-waveform graph with data stored.

In the prototype system, a 532-nm laser is used, with an output frequency of 10  kHz. The sensor array we adopt is composed of four-unit APDs, each of which works in linear mode and generates an output of an electrical signal with a certain gain and with a certain noise suppression. The data acquisition system provides four interfaces for data input, one interface for the trigger signal input, and one interface for a high-frequency signal input as the sampling rate control. The highest sampling the system can reach is up to 1  Gs· ps. The data upper-controlling and processing system is located on the PC side with an Ethernet cable connected to a data acquisition system, which is developed based on the virtual instrument technology. The upper monitor and full-waveform display in the site trial is shown in Fig.  6. Since the frequency of the laser emitting is 10  kHz, echo signals appear periodically. When the object with special geometric characteristics is projected, multiple individuals return and overlap within the current period.

Fig.  6.

	Figure Option ViewDownloadNew Window
	Fig.  6. Software interface of the upper monitor over the Ethernet and the full waveform display. It is developed based on LabVIEW. The main setting parameters include an IP address and port number of server, data length, data format, sending mode, and time interval, etc.

Full-waveform data processing system consists of three function modules, namely signal enhancement & de-noising, smoothing & deburring, and fitting & decomposition. The flow chart of full waveform data processing is shown in Fig.  7.

Fig.  7.

	Figure Option ViewDownloadNew Window
	Fig.  7. Flow chart of full waveform data processing.

3.4.1. Signal enhancement & de-noising

As mentioned above, noise sources in the lidar system mainly come from atmospheric effects, target interaction, and receiver internal influence, which can be divided into two kinds, namely Gaussian noise and speckle noise. There are various effective methods for Gaussian noise; for instance, the wavelet denoising method has good time– frequency characteristics, and wavelet basis selection is flexible.^{[17, 18]} Besides, wavelet denoising can deal with the edges, burr, and some other details of signals, to some extent. However, there is barely an effective processing method on speckle noise, which has been a bottleneck in the application of 3D imaging of the laser radar system. Owing to the multiplicative characteristic of speckle noise, a homomorphic filtering method is first considered. The principle is as follows: a noisy signal is first imposed on logarithmic transformation, thus converting the multiplicative noise into the additive; then frequency domain filtering is utilized to eliminate the ineffective part in the valid signal; finally the output signal is converted back to the original linear region by exponential transformation. The biggest problem of this method used in our system is that it requires a high performance in additive noise preprocessing, which is hardly realized in remote experiments and the low SNR condition.

Another method for dealing with a non-stationary random signal is the independent component correlation algorithm (ICA) method. Blanco et al.^[19] applied it to speckle signals, and proposed a speckle ICA (SICA) method which takes into account the multiplicative model, and finally it is proved to be clearly improved on the basis of a standard ICA method for the speckle mixture condition. The limitation of this method is that it has a requirement for the number of noisy signal samples, and also has a higher algorithmic complexity.

Pulse accumulation has been widely used as an effective method for a non-stationary random signal in a radar system. Especially, on condition that the echo signals are weak and even when they are submerged in noise, this method is effective for improving SNR, to extract useful signals. Actually, pulse accumulation has good performances in noise suppression for both linear Gaussian noise and also speckle noise. Assuming S to be a useful signal, and N a random noise, the noisy signal within a cycle can be expressed as

After noisy signals in M cycles are added up, then

Then the SNR is

Since noises within different cycles are statistically independent of each other, namely

Thus

The above equation means that the effect of speckle noise, which is involved in the envelope of a useful signal, will become almost M times as great as the original effect after pulse accumulation. To sum up, we adopt pulse accumulation as an approach to signal enhancement and denoise preprocessing. Owing to the fact that the number of samples for accumulation is limited by the pulse frequency and the platform speed, the wavelet filter is utilized for further denoising. However, there is something more crucial that the improvement of SNR and the integrity of details of the waveform should be taken into account, simultaneously. On the basis of the principle, an appropriate wavelet base, filter length, and decomposition level are selected. The processing results of the above-mentioned de-noising methods in our system are shown in Fig.  8, from which we can demonstrate their efficiency quantitatively.

Fig.  8.

Figure Option ViewDownloadNew Window

Fig.  8. Diagram of waveforms after enhancement and de-noising. (a) The original echo waveform presents an SNR of poor quality, which is about − 5.95  dB; (b) after accumulation within eight cycles, the SNR is improved by about eight times, larger than before, namely 3.20  dB; (c) the result after wavelet filtering, which is significantly enhanced without missing the details of the waveform, the SNR of which is 9.38  dB.

3.4.2. Smoothing

The purpose of the smoothing module is to deal with the detail problems of noisy signals, like local noise, burr, random fluctuation, etc. In this context, optimized signals are obtained, which are smoother and of higher SNR. There are some common smoothing methods, like polynomial smoothing, moving average, Gaussian smoothing, Vondrak smoothing, etc.

As mentioned above, one important principle for the noisy signal processing is to implement on the premise of ensuring the integrity of the wanted signal. To select an optimal approach as the smoothing module for our lidar system, various smoothing methods are compared. The Vondrak data smoothing method is proposed by an astronomer from Czech named Vondrak.^[20] This smoothing method works well on condition that both the change regulation of observation data and functional form for fitting are unknown. The basic criterion of this method is

where

with x(t_i) being the value to be smoothed. Here we define a parameter ε named the smoothing factor as the degree of smoothing, then ε = 1/λ ², which means that the degree of smoothing increases as ε decreases. Figure  9 shows the effect of Vondrak smoothing on the echo data. As shown in the figure, the original data have many burrs. To smooth the waveform, we utilize the Vondrak smoothing methods of three values of smoothing factor ε = 100, 10, 1 (marked with a dashed line, dotted– dash line, and dotted line respectively). The results show that the smoothing result under ε = 100 remains good integrity while burrs are eliminated effectively.

Fig.  9.

	Figure Option ViewDownloadNew Window
	Fig.  9. Vondrak smoothing results under different factors. Original data are marked with a solid line, which is smoothed with the three smoothing factor ε = 100, 10, 1 (marked with a dashed line, dotted– dash line, dotted line respectively).

3.4.3. Fitting & decomposition

In order to identify the component of the full-waveform, the fitting & decomposition method is utilized to acquire the parameters like amplitude, position, and half-width of the different component. One effective and highly precise approach to full-waveform lidar decomposition is one nonlinear fitting method based on Levenberg– Marquardt (LM), which combines the local convergence characteristics of the Gauss– Newton iterative algorithm, also with global characteristics of the gradient descent, thus having faster iteration and higher accuracy. Specifically, Hofton et al. proposed one composition method, which combines the nonnegative least-squares method (LSM), and the LM method. The implementation procedure is described as follows.^[21]

However, most of the decomposition methods still regard the component of the full-waveform as the Gaussian form. This means that the full-waveform is expressed as

Actually, the time propagation of a typical laser pulse is described as

where τ = T_1/2/3.5, with T_1/2 representing the full width at half maximum. So there is obvious deviation because of using the Gaussian component fitting method. To eliminate the deviation, the fitting template needs to be amended according to Eq.  (18). Then the new expression of a full-waveform is

The fitting results with both the Gaussian model and laser model are shown in Fig.  10. In order to compare and evaluate the quality of fitting, we propose a criterion called the coefficient of determination (or goodness of fit): a greater coefficient of determination means a better fitting performance of the observed data and a more accurate result of decomposition. The coefficient of determination is defined as

where TSS means the total square sum of the observed data, expressed as ESS represents the error square sum; RSS denotes the residual square sum. More specifically,

where is the fitting data.

Fig.  10.

	Figure Option ViewDownloadNew Window
	Fig.  10. Fitting results based on (a) Gaussian template and (b) laser template on simulated echo waveform. Two specific ranges are marked in the figure, within which a comparison of goodness of fit between them is obvious.

In Fig.  10(a), the value of R² is 0.9934, while it is 0.9935 in Fig.  10(b). It seems that there is little difference in fitting performance between both of the fitting methods. Actually, it is inappropriate to calculate R² directly through Eq.  (22). Within the time range from 226 to 246, we find that the fitting result with the laser model is much better than with the Gaussian fitting, just under the same condition with other ranges around the pulse peaks. Contrarily, fitting with the Gaussian model performs better within those ranges around the tail (end) of the pulse, for instance, the ranges between 303 to 319 shown in Fig.  10. The reason is that the pulse tails are distorted more under the effect of certain intensity noise than the pulse peaks. Therefore, the criterion needs to be amended as follows:

In Eq.  (23), we introduce a function W, which consists of a plurality of rectangular window functions, that is

where K is the number of rectangular window functions; R_{M ^j} (n − n_j) means the j-th rectangular window function, of which the width is M ^j, and the initial position is n_j. Conventionally, R_M (n) is defined as

Then the crucial point is how to determine the parameters K, M ^j, and n_j. Here, we propose one method and its steps are as follows.

Step 1 Initialize parameters: K = N, i = 1, j = 1.

Step 2 Calculate the starting frequency and cutoff frequency of the specified threshold (e.g. 3  dB) of whole sub-pulses.

Step 3 Let .

Step 4 If , i → i + 1, K = K − 1, repeat Step  4.

Step 5 If , , i → i + 1. If i ≥ N, the step ends; otherwise repeat Step  3.

As can be seen from the above, the updated method is essentially an approach to points’ selection and evaluation. With the new evaluation criteria, R̃ ² can be obtained from Fig.  10, which is 0.9984 in Fig.  10(a) and 0.9988 in Fig.  10(b). It shows that there is a maximum deviation of two points in the decomposition result after Gaussian fitting, which means 0.3-m deviation when the sampling rate is 1  Gs· ps. Compared with the Gaussian model, the fitting result with the laser model has better goodness of fit and thus higher accuracy.

In this section, we are to use the depth and intensity information from the data processing module to generate a digital elevation model (DEM) of the terrain. The implementation process is shown in Fig.  11. Two-dimensional (2D) images, including both a range image and an intensity image, are produced from the decomposition data. Then data filtering is utilized to finish the clustering, segmentation on point cloud, which is classified as ground data and non-ground data. Data filtering algorithms for lidar have been refined in recent years. The basic principles of different algorithms are unified, i.e., a particular screening method is selected to eliminate the non-ground point and wrong point from the point cloud according to the difference in elevation information among different types. Specifically, a typical filtering algorithm based on triangular irregular networks is proposed by Axelsson.^[22] The basic principle of the algorithm is described as follows. A sparse TIN is created from seed points, which are selected from those points with relatively low elevations within a sufficiently large area. Through an iterative method, the TIN network is constantly densified. The rule for densification is that the points from a new triangle mesh are added into the network once it satisfies the requirements for threshold constraint. The algorithm handles the surface with discontinuities in dense city areas. The commercial software TerraSolid adopts the algorithm as its filtering approach.

Fig.  11.

	Figure Option ViewDownloadNew Window
	Fig.  11. Flow chart of 3D reconstruction to build the DEM for Terrain object.