# Observation Selection, Total Variation, and L-Curve Methods for LiDAR Data Denoising

Authors: Jie Zhang; Pak Wai Chan; Michael K. Ng
Journal: Advances in Meteorology
 (2022)

---

## Abstract

In this article, we propose a light detection and ranging (LiDAR) data denoising scheme for wind profile observation as a part of quality control procedure for wind velocity monitoring and windshear detection. The proposed denoising scheme consists of several components. (i) It selects LiDAR observations according to their SNR values so that serious noisy data can be removed. (ii) A polar-based total variation smoothing term is employed to regularize LiDAR observations. (iii) The regularization parameters are automatically determined to balance the data-fitting term and the total variation smoothing term. Numerical results for LiDAR data collected at the Hong Kong International Airport are reported to demonstrate that the denoising performance of the proposed method is better than that of the testing LiDAR data denoising schemes in the literature.


---

## Body

## 1. Introduction

Light detection and ranging (LiDAR) technique [1, 2] is a remote sensing tool that plays a significant role in environmental monitoring sciences. It is widely used in meteorological data observing. Generally, LiDAR emits a beam of light to the observation region. Some of this light would be backscattered towards the LiDAR receiver since it interacts with the medium or particles under observation [3]. The backscattered light captured by the LiDAR receiver is used to determine the characteristics of the observation area, e.g., the velocity of wind. Due to the impact of measurement environments and some other reasons, there would be some observation errors and very noisy observations in the observational LiDAR data as the range of observation increases [4]. It can have a serious effect in different LiDAR data applications such as windshear detection. Therefore, it is indispensable to develop an effective denoising method as a part of quality control for LiDAR observational data to improve the data quality and remove bad observations.Several denoising methods have been developed to improve the quality of LiDAR data. There are mainly two different types of methods for LiDAR data denoising. One is for the time-varying but location-fixed LiDAR data. For example, a stationary wavelet domain spatial filtering-based denoising method was proposed by Yin et al. [5]. The method can effectively remove noise and detect the sudden change of LiDAR signal. Hassanpour [6, 7] proposed a singular value decomposition-based Savitzky-Golay approach for signal denoising. Similarly, Azadbakht et al. [8] employed the Savitzky-Golay method for full-waveform LiDAR data denoising. The second one is for the time-varying and distance-varying LiDAR data. For instance, Wu et al. [9] proposed a biorthogonal discrete wavelet transform (DWT) with a distance-dependent threshold algorithm to do the line-of-sight wind velocity denoising. In [10], Wu et al. studied the empirical mode decomposition (EMD) method [11] to analyze LiDAR data. Liu et al. applied the EMD method to a Doppler wind LiDAR acquisition system and got much better denoising results than the original method in [12]. In [13], Zhang et al. combined the EMD method with the Savitzky-Golay filtering algorithm, which can retain more features of LiDAR signal. Li et al. proposed a LiDAR denoising method based on ensemble empirical mode decomposition in [14]. This method can overcome the mode mixing phenomenon that occurs with the EMD method. Also Tian et al. improved the EMD method for range and frequency analysis in LiDAR data and proposed an automatic EMD denoising method in [15]. In [16], the EMD-CIIT method [17] was applied to reduce the noise of large-scale LiDAR data.However, the abovementioned methods cannot remove the bad observations in LiDAR data. To address this issue, several methods including bad data removal were used in the LiDAR data applications. The most commonly used methods are based on the signal-to-noise ratio (SNR). For example, Baranov et al. [18] filtered the bad observations in LiDAR data by the threshold level of SNR and then applied a smoothing algorithm to improve the data quality. In [19], Newsom et al. removed the bad observations based on an SNR threshold level. The Hong Kong Observatory (HKO) currently used a SNR-based quality control method to preprocess the LiDAR data collected at the Hong Kong International Airport (HKIA) for windshear detection (refer to Algorithm 3 for more detailed description). However, the SNR threshold level the above methods used is really empirical that requires numerous tests based on the observational data. Moreover, we remark that all the abovementioned methods do not consider LiDAR data at the other azimuth angles together in the data processing procedure.In this article, we aim to investigate a mathematical denoising scheme that can not only reduce the noise level for a whole scan but also remove the bad observations in LiDAR data adaptively based on the corresponding SNR values. By analysis of the LiDAR data collected at the Hong Kong International Airport (see Section2 for more details), we propose a LiDAR data denoising method based on the minimization of an objective function containing (i) the data-fitting term between the observed LiDAR data and the denoised data; (ii) the polar-based total variation regularization term that is used to smooth LiDAR observations; and (iii) the weighting term of LiDAR observations that is employed to control whether LiDAR observations are used in the denoising procedure based on their SNR values and neighborhood observations. In the optimization scheme, we also propose an L-curve selection method for several regularization parameters to be used for balancing the contribution of the above three terms in the objective function. The whole scan data are used in the proposed scheme by the polar total variation method instead of the only data at the same azimuth angle used by the other methods given in [9, 16, 18, 19]. The SNR-based weighting term makes it possible to remove the very bad observations more flexibly instead of just relying on a fixed empirical threshold. Numerical results demonstrate the usefulness of the proposed denoising scheme compared with testing LiDAR data denoising schemes in the literature.This paper is organized as follows. In Section2, we introduce our proposed denoising scheme and the parameter selection method. In Section 3, results and discussions are presented. Finally, some concluding remarks are given in Section 4.


## 2. Methods

In this section, the information about LiDAR data we investigate in this paper is given in Subsection2.1. Next, we introduce the proposed model in Subsection 2.2. Also, the algorithm and parameter selection method are given in Subsections 2.3 and 2.4, respectively.


### 2.1. Data Sets

The Hong Kong International Airport (HKIA) is located at the place lying to the north of Lantau Island that is quite mountainous with heights ranging from 300 m to 900 m. Due to the complex terrain near the airport, it is necessary to collect LiDAR data of wind velocities and observe any windshear phenomena appearing over the flight paths of the airport. To provide timely windshear alerting, the Hong Kong Observatory devised a Doppler LiDAR system (see [20, 21] for more details). Due to the highly cluttered environment around HKIA such as vehicles, derricks, barges, windmills, and cable cars, there are lots of noise and bad observations in the observational LiDAR data. For example, we show in Figure 1(a) the LiDAR radial velocity data of conical scan, where the radius and the polar angle of the scan refer to the slant range and the azimuth angle of LiDAR beam, respectively. In Figure 1(b), we show the signal-to-noise ratio (SNR) of the measured wind velocities corresponding to Figure 1(a). It is obvious that there are noise and some outliers whose SNR values are from −10 to −60 in the observational data. Therefore, it is significant to develop an efficient data denoising method to remove bad observations and improve the quality of observational LiDAR data.Figure 1
An example of the LiDAR data on 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25. (a) Wind velocities and (b) signal-to-noise ratio of LiDAR observations. The observations were collected at the Hong Kong International Airport under LiDAR system by Hong Kong Observatory [21, 22].
(a)(b)The set of data used in this study was collected at HKIA from 1 March to 31 March 2015, including the 147 windshear cases that reported in the pilot report and several nonwindshear cases. Each scan of windshear cases tookyi∈ℝn about 25 seconds (see one example in Figure 1).


### 2.2. The Proposed Model

Let be the LiDAR data observed at azimuth angleθi. For simplicity, we assume that there are m azimuth angles (θ1<θ2<⋯<θm) to be recorded in between 0° and 359°, and also there are n range values to be recorded, i.e.,(1)yi=yi,1,yi,2,…,yi,n,1≤i≤m.Note that the locations of such range values are not necessary to be uniform, and the distance between the LiDAR centre and the observed valueyi,j is equal to rj. We are interested to compute the denoising LiDAR data as follows:(2)xi=xi,1,xi,2,…,xi,n,1≤i≤m,according to the given SNR values:(3)si=si,1,si,2,…,si,n,1≤i≤m,from the observed LiDAR data. When yi,j is a missing LiDAR observation, the corresponding si,j can be set to be −∞. In total, there are nm observations and SNR values and nm unknowns covered in the conical scan in a two-dimensional plane. In the proposed minimization model, there are three components in the objective function.(i)
The first term is the data-fitting term between the observed LiDAR datayi and the denoised data xi (1≤i≤m). In order to determine whether yi,j at the azimuth angle θi and range value rj to be used in the model, we incorporate a nonnegative weight wi,j for yi,j. When the value of wi,j is equal to zero, the LiDAR observation yi,j is not used. The resulting data-fitting term is given by(4)12∑i=1m∑j=1nwi,jyi,j−xi,j2.It is clear that when the LiDAR data observation is removed (wi,j=0), the corresponding data-fitting term does not appear.(ii)
In the model, we would like to smooth LiDAR observations over all possible range and azimuth values together. In image processing, total variation regularization techniques based on rectangular pixel-based form were shown to be a very useful denoising method (see, for instance, [23–25]). Here, the polar-based total variation regularization term is proposed and employed as follows:(5)TVxi,j=xi,j−xi,j−1dxi,j,xi,j−1+xi,j−xi−1,jdxi,j,xi−1,j,wheredxi,j,xi,j−1 refers to the distance between the two observed range values at the same azimuth angle, i.e.,(6)dxi,j,xi,j−1=rj−rj−1,anddxi,j,xi,j−1 refers to the distance between the same observed range value at the two adjacent azimuth angles θi and θi−1, i.e.,(7)dxi,j,xi−1,j=2rjsinθi−θi−12.In (7), the first and the second terms are the one-dimensional total variation regularization. However, we apply the one-dimensional total variation regularization to the polar-based LiDAR data observations by using their actual distances in the formula. The combined total variation regularization term is given as follows:(8)∑i=1m∑j=1nwi,jTVxi,j.We see from (7) and (8) that the denoised values xi,j are coupled into the regularization formula, and therefore the denoising procedure is adapted to the whole polar-based LiDAR data observations instead of observations in range values only such as the methods we introduce in Section 1.(ii)
In the model, we should include SNR valuessi,j to determine whether the LiDAR observation yi,j is used in the denoising procedure. Here, we propose the following term in the model:(9)λ1∑i=1m∑j=1nexp−si,jwi,j2+λ2∑i=1m∑j=1nexpsi,jwi,j−12,where λ1 and λ2 are the positive numbers to control the balance between the two terms in (9). We note in (9) that when si,j is too small (a negative number), the LiDAR observation is very noisy, and therefore exp−si,j is large and expsi,j is small. The optimization process drives wi,j to be zero. Similarly, when si,j is not small (the LiDAR observation is not noisy), exp−si,j is small and expsi,j is large and therefore the optimization process drives wi,j to be one. Finally, we also incorporate nonnegativity constraint on wi,j in the optimization model:(10)wi,j≥0,1≤i≤m,1≤j≤n.For the ease of presentation, we denotex=xi,j and w=wi,j. According to (4), (6), (9), and (10), the combined optimization model is given as follows:(11)x,w=arg minx,w12∑i=1m∑j=1nwi,jyi,j−xi,j2+α∑i=1m∑j=1nwi,jTVxi,j+χw+λ1∑i=1m∑j=1nexp−si,jwi,j2+λ2∑i=1m∑j=1nexpsi,jwi,j−12,where α is a positive number to balance the contribution of the total variation regularization term and the other terms and χWw is the characteristic function of the nonnegativity constraint set W=w:wi,j≥0,1≤i≤m,1≤j≤n(12)χWw=0,w∈W,+∞,otherwise.


### 2.3. The Algorithm

In this subsection, we will introduce the algorithm for the proposed model. Since bothdxi,j,xi−1,j and dxi,j,xi−1,j are fixed in LiDAR observations, for simplicity, we use τi,j and μi,j to represent them, respectively, in the following discussion. Let pi,j=xi,j−xi,j−1,qi,j=xi,j−xi−1,j. We can rewrite (11) as follows:(13)x,w,p,q=argminx,w,p,q12∑i=1m∑j=1nwi,jyi,j−xi,j2+α∑i=1m∑j=1nwi,jpi,jτi,j+α∑i=1m∑j=1nwi,jqi,jμi,j+χWw+λ1∑i=1m∑j=1nexp−si,jwi,j2+λ2∑i=1m∑j=1nexpsi,jwi,j−12,s.t.pi,j=xi,j−xi,j−1,qi,j=xi,j−xi−1,j.The augmented Lagrangian of the above-constrained optimization problem is given as follows:(14)Lγ1,γ2x,w,p,q;Λ,Γ=12∑i=1m∑j=1nwi,jyi,j−xi,j2+α∑i=1m∑j=1nwi,jpi,jτi,j+α∑i=1m∑j=1nwi,jqi,jμi,j+χWw+λ1∑i=1m∑j=1nexp−si,jwi,j2+λ2∑i=1m∑j=1nexpsi,jwi,j−12+∑i=1m∑j=1nΛi,jpi,j−xi,j−xi,j−1+∑i=1m∑j=1nΓi,jqi,j−xi,j−xi−1,j+γ12∑i=1m∑j=1npi,j−xi,j−xi,j−12+γ22∑i=1m∑j=1nqi,j−xi,j−xi−1,j2,where γ1,γ2>0 are penalization parameters and ΛandΓ are Lagrange multipliers.With an initial guess ofx0,w0,p0,q0;Λ0,Γ0, the iterations of the alternating direction method of multipliers [26] are given as follows:(15)xk+1=argminxLγ1,γ2x,wk,pk,qk;Λk,Γk,wk+1=argminwLγ1,γ2xk+1,w,pk,qk;Λk,Γk,pk+1=argminpLγ1,γ2xk+1,wk+1,p,qk;Λk,Γk,qk+1=argminqLγ1,γ2xk+1,wk+1,pk+1,q;Λk,Γk,Λi,jk+1=Λi,jk+γ1pi,jk+1−xi,jk+1−xi,j−1k+1,Γi,jk+1=Γi,jk+γ2qi,jk+1−xi,jk+1−xi−1,jk+1.Next, we demonstrate how to solve the abovementioned subproblems (15) with respect to x,w,p,q.Let us consider thex-subproblem in (15):(16)xk+1=argminx12∑i=1m∑j=1nwi,jyi,j−xi,j2+γ12∑i=1m∑j=1npi,jk+Λi,jkγ1−xi,j−xi,j−12+γ22∑i=1m∑j=1nqi,jk+Γi,jkγ2−xi,j−xi−1,j2.Noting that the data-fitting term will vanish if there is no observation datayi,j. The optimality condition of x-subproblem is given by(17)wi,jk+2γ1+2γ2xi,j−γ1xi,j+1+xi,j−1−γ2xi−1,j+xi+1,j=wi,jkyi,j+γ1pi,jk−pi,j+1k+γ2qi,jk−qi+1,jk+Λi,jk−Λi,j+1k+Γi,jk−Γi+1,jk,ifyi,jexists,2γ1+γ2xi,j−γ1xi,j+1+xi,j−1−γ2xi−1,j+xi+1,j=γ1pi,jk−pi,j+1k+γ2qi,jk−qi+1,jk+Λi,jk−Λi,j+1k+Γi,jk−Γi+1,jk,otherwise.We can solve the above set of equations by using the conjugate gradient (CG) method [27, 28].Forw-subproblem, we have(18)wk+1=argminw≥012∑i=1m∑j=1nwi,jyi,j−xi,jk+12+α∑i=1m∑j=1nwi,jpi,jkτi,j+α∑i=1m∑j=1nwi,jqi,jkμi,j+λ1∑i=1m∑j=1nexp−si,jwi,j2+λ2∑i=1m∑j=1nexpsi,jwi,j−12.We note thatwi,j are decoupled, and we only need to deal with scalar optimization problems. More precisely, if yi,j exists, the scalar optimization problem is given by(19)wi,jk+1=argminw≥012wi,jyi,j−xi,jk+12+αwi,jpi,jkτi,j+αwi,jqi,jkμi,j+λ1exp−si,jwi,j2+λ2expsi,jwi,j−12.Ifyi,j does not exist, it is equal to(20)wi,jk+1=argminw≥0αwi,jpi,jkτi,j+αwi,jqi,jkμi,j+λ1exp−si,jwi,j2+λ2expsi,jwi,j−12.Both the above two scalar optimization problems are convex. And one can readily get that(21)wk+1=maxw¯k+1,0,where(22)w¯i,jk+1=2λ2expsi,j−1/2yi,j−xi,jk+12−αpi,jk/τi,j−αpi,jk/μi,j2λ1exp−si,j+2λ2expsi,j,ifyi,jexists,2λ2expsi,j−αpi,jk/τi,j−αqi,jk/μi,j2λ1exp−si,j+2λ2expsi,j,otherwise.Forp-subproblem and q-subproblem, they are given by(23)pk+1=argminpα∑i=1m∑j=1nwi,jpi,jτi,j+γ12∑i=1m∑j=1npi,j+Λi,jkγ1−xi,jk+1−xi,j−1k+12,and(24)qk+1=argminqα∑i=1m∑j=1nwi,jqi,jμi,j+γ22∑i=1m∑j=1nqi,j+Γi,jkγ2−xi,jk+1−xi−1,jk+12,respectively. These two subproblems can be solved by using the soft-thresholding technique [29], and their solutions are given by(25)pi,jk+1=maxxi,jk+1−xi,j−1k+1−Λi,jkγ1−αwi,jk+1γ1τi,j,0signxi,jk+1−xi,j−1k+1−Λi,jkγ1,qi,jk+1=maxxi,jk+1−xi−1,jk+1−Γi,jkγ2−αwi,jk+1γ2μi,j,0signxi,jk+1−xi−1,jk+1−Γi,jkγ2,where signp≔o/p.Finally, the overall algorithm is summarized in Algorithm1.Algorithm 1:
Input: Noisy data y, SNR data s, and parameters α,λ1,λ2,γ1,γ2.Initialization:x0=y,w0=1,p0=0,q0=0,Λ0=0,Γ0=0,k=0.(1)
while Stopping criteria is not satisfied do(2)
Solvexk+1 by using CG for (17)(3)
Solvewk+1 by (21)(4)
Findpk+1 by (24)(5)
Findqk+1 by (25)(6)
Update the multipliers byΛi,jk+1=Λi,jk+γ1pi,jk+1−xi,jk+1−xi,j−1k+1 and Γi,jk+1=Γi,jk+γ2qi,jk+1−xi,jk+1−xi−1,jk+1(7)
k←k+1(8)
end while


### 2.4. The Calculation of Parameters

In the proposed denoising scheme, there are several parameters to be determined. In this subsection, we incorporate the L-curve method to determine the values of parameters. The L-curve is a tradeoff curve between two quantities that both need to be controlled and balanced [30, 31]. It is widely used in the engineering and applied mathematics field.First, the regularization parameterα in (14) is crucial in our proposed method since it controls the ratio between the data-fitting term and the total variation regularization term. When the value of α is large (small), the importance of the total variation regularization term is large (small) and the denoising results are more (less) regularized. Moreover, there are two parameters λ1 and λ2 in (14), and they can affect the values of wi,j according to the SNR values of LiDAR observations. When the ratio λ1/λ2 is large (i.e., λ1 is larger than λ2), the term λ1∑i=1m∑j=1mexp−si,jwi,j2 is more important than the term λ2∑i=1m∑j=1mexp−si,jwi,j−12. It follows that when LiDAR observations have small SNR values, there would be more zero values of wi,j in the optimization procedure. In contrast, when the ratio λ1/λ2 is small and LiDAR observations have large SNR values, there would be more values of wi,j to be one in the optimization procedure.In the L-curve method, we consider the following two quantities:(26)Dxα,λ1,λ2,wα,λ1,λ2≔−log10∑i=1m∑j=1nwi,jα,λ1,λ2yi,j−xi,jα,λ1,λ22∑i=1m∑j=1nwi,jα,λ1,λ2,and(27)Sxα,λ1,λ2,wα,λ1,λ2≔−log10∑i=1m∑j=1nwi,jα,λ1,λ2TVxi,jα,λ1,λ2∑i=1m∑j=1nwi,jα,λ1,λ2.The first quantityDxα,λ1,λ2,wα,λ1,λ2 refers to the degree of data fitting when x and w are computed by Algorithm 1 for fixed α, λ1, and λ2. Similarly, the second quantity Sxα,λ1,λ2,wα,λ1,λ2 refers to the degree of smoothness in the regularization for the computed x and w. In order to compare the degrees of data fitting and smoothness with respect to different numbers of denoising values wi,j across different values of λ1 and λ2, the weighted average of data fitting and smoothness are accounted in the two quantities. The main idea of the L-curve method is to select the suitable values of α, λ1, and λ2 such that both Dxα,λ1,λ2,wα,λ1,λ2 and Sxα,λ1,λ2,wα,λ1,λ2 are balanced.As an illustration, we apply Algorithm1 to the LiDAR observations in Figure 1 for different values of α, λ1, and λ2. More precisely, in Figure 2(a), we fix λ1 (100) and λ2 (0.5,5,50,500) and compute x and w for different values of α (0.01,0.05,0.1,0.3,0.5,0.8,1,5). In Figure 2(a), we plot Dxα,λ1,λ2,wα,λ1,λ2 and Sxα,λ1,λ2,wα,λ1,λ2. The lines of points of D and S with respect to α are generated for each fixed λ1 and λ2. In total, there are four lines of points in Figure 2(a). We observe in each line that when α is small, D is large and S is small, and thus the points appear in the left hand side of the line. When α is large, D is small and S is large, and thus the points appear in the right hand side of the line. We find that there are several corner points in the line, and they refer to several large rates of change of data fitting with respect to the change of smoothness. Here, we can pick up the corner point with the largest rate of change and employ the corresponding value of α for regularization. Next, we select the values of λ1 and λ2 by comparing the selected corner points of different lines. Here, we can pick up the selected corner point with the largest D (the data fitting is good) and the largest S (the smoothness is large). According to the plot in Figure 2(a), we select α∗=0.5, λ1∗=100,and λ2∗=500.Figure 2
L-curves of different parametersα, λ1, and λ2.
(a)(b)Based on the above idea, we summarize our L-curve method with automatic selection ofα, λ1, and λ2 in Algorithm 2.Algorithm 2: L-curve Method.
Input: The set of parameters values α, λ1, λ2.Initialization: Corner point c0=0,0(1)
while the set of parameter values is not empty do(2)
Select (λ1, λ2) from the parameter set(3)
Find the corner pointc∗ of the generated lines of points Dxα,λ1,λ2,wα,λ1,λ2,Sxα,λ1,λ2,wα,λ1,λ2 from a set of α values, and record the corresponding α˜.(4)
If the position ofc∗ is higher than the position of c0 in the D-S plot, set c0=c∗, α∗=α˜, λ1∗=λ1, λ2∗=λ2(5)
end whileFrom Figure2, it is obvious that the model is very sensitive for different values of the regularization parameter α due to its balance effect for data fitting and smoothing. However, from Figure 2(a), for fixed values of parameters α and λ1, the changes of D and S for different values of λ2 are small (less than 1), which means that the model is pretty robust for parameter λ2. To evaluate the robustness of parameter λ1, we plot the four L-curves with respect to α (0.01,0.05,0.1,0.3,0.5,0.8,1,5) for fixed λ2 in Figure 2(b), where λ2=100, λ1 varies in 0.5,5,50,500 for different curves. Similar to the results of parameter λ2, one can readily find that the proposed model is quite robust for parameter λ1.


## 2.1. Data Sets

The Hong Kong International Airport (HKIA) is located at the place lying to the north of Lantau Island that is quite mountainous with heights ranging from 300 m to 900 m. Due to the complex terrain near the airport, it is necessary to collect LiDAR data of wind velocities and observe any windshear phenomena appearing over the flight paths of the airport. To provide timely windshear alerting, the Hong Kong Observatory devised a Doppler LiDAR system (see [20, 21] for more details). Due to the highly cluttered environment around HKIA such as vehicles, derricks, barges, windmills, and cable cars, there are lots of noise and bad observations in the observational LiDAR data. For example, we show in Figure 1(a) the LiDAR radial velocity data of conical scan, where the radius and the polar angle of the scan refer to the slant range and the azimuth angle of LiDAR beam, respectively. In Figure 1(b), we show the signal-to-noise ratio (SNR) of the measured wind velocities corresponding to Figure 1(a). It is obvious that there are noise and some outliers whose SNR values are from −10 to −60 in the observational data. Therefore, it is significant to develop an efficient data denoising method to remove bad observations and improve the quality of observational LiDAR data.Figure 1
An example of the LiDAR data on 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25. (a) Wind velocities and (b) signal-to-noise ratio of LiDAR observations. The observations were collected at the Hong Kong International Airport under LiDAR system by Hong Kong Observatory [21, 22].
(a)(b)The set of data used in this study was collected at HKIA from 1 March to 31 March 2015, including the 147 windshear cases that reported in the pilot report and several nonwindshear cases. Each scan of windshear cases tookyi∈ℝn about 25 seconds (see one example in Figure 1).


## 2.2. The Proposed Model

Let be the LiDAR data observed at azimuth angleθi. For simplicity, we assume that there are m azimuth angles (θ1<θ2<⋯<θm) to be recorded in between 0° and 359°, and also there are n range values to be recorded, i.e.,(1)yi=yi,1,yi,2,…,yi,n,1≤i≤m.Note that the locations of such range values are not necessary to be uniform, and the distance between the LiDAR centre and the observed valueyi,j is equal to rj. We are interested to compute the denoising LiDAR data as follows:(2)xi=xi,1,xi,2,…,xi,n,1≤i≤m,according to the given SNR values:(3)si=si,1,si,2,…,si,n,1≤i≤m,from the observed LiDAR data. When yi,j is a missing LiDAR observation, the corresponding si,j can be set to be −∞. In total, there are nm observations and SNR values and nm unknowns covered in the conical scan in a two-dimensional plane. In the proposed minimization model, there are three components in the objective function.(i)
The first term is the data-fitting term between the observed LiDAR datayi and the denoised data xi (1≤i≤m). In order to determine whether yi,j at the azimuth angle θi and range value rj to be used in the model, we incorporate a nonnegative weight wi,j for yi,j. When the value of wi,j is equal to zero, the LiDAR observation yi,j is not used. The resulting data-fitting term is given by(4)12∑i=1m∑j=1nwi,jyi,j−xi,j2.It is clear that when the LiDAR data observation is removed (wi,j=0), the corresponding data-fitting term does not appear.(ii)
In the model, we would like to smooth LiDAR observations over all possible range and azimuth values together. In image processing, total variation regularization techniques based on rectangular pixel-based form were shown to be a very useful denoising method (see, for instance, [23–25]). Here, the polar-based total variation regularization term is proposed and employed as follows:(5)TVxi,j=xi,j−xi,j−1dxi,j,xi,j−1+xi,j−xi−1,jdxi,j,xi−1,j,wheredxi,j,xi,j−1 refers to the distance between the two observed range values at the same azimuth angle, i.e.,(6)dxi,j,xi,j−1=rj−rj−1,anddxi,j,xi,j−1 refers to the distance between the same observed range value at the two adjacent azimuth angles θi and θi−1, i.e.,(7)dxi,j,xi−1,j=2rjsinθi−θi−12.In (7), the first and the second terms are the one-dimensional total variation regularization. However, we apply the one-dimensional total variation regularization to the polar-based LiDAR data observations by using their actual distances in the formula. The combined total variation regularization term is given as follows:(8)∑i=1m∑j=1nwi,jTVxi,j.We see from (7) and (8) that the denoised values xi,j are coupled into the regularization formula, and therefore the denoising procedure is adapted to the whole polar-based LiDAR data observations instead of observations in range values only such as the methods we introduce in Section 1.(ii)
In the model, we should include SNR valuessi,j to determine whether the LiDAR observation yi,j is used in the denoising procedure. Here, we propose the following term in the model:(9)λ1∑i=1m∑j=1nexp−si,jwi,j2+λ2∑i=1m∑j=1nexpsi,jwi,j−12,where λ1 and λ2 are the positive numbers to control the balance between the two terms in (9). We note in (9) that when si,j is too small (a negative number), the LiDAR observation is very noisy, and therefore exp−si,j is large and expsi,j is small. The optimization process drives wi,j to be zero. Similarly, when si,j is not small (the LiDAR observation is not noisy), exp−si,j is small and expsi,j is large and therefore the optimization process drives wi,j to be one. Finally, we also incorporate nonnegativity constraint on wi,j in the optimization model:(10)wi,j≥0,1≤i≤m,1≤j≤n.For the ease of presentation, we denotex=xi,j and w=wi,j. According to (4), (6), (9), and (10), the combined optimization model is given as follows:(11)x,w=arg minx,w12∑i=1m∑j=1nwi,jyi,j−xi,j2+α∑i=1m∑j=1nwi,jTVxi,j+χw+λ1∑i=1m∑j=1nexp−si,jwi,j2+λ2∑i=1m∑j=1nexpsi,jwi,j−12,where α is a positive number to balance the contribution of the total variation regularization term and the other terms and χWw is the characteristic function of the nonnegativity constraint set W=w:wi,j≥0,1≤i≤m,1≤j≤n(12)χWw=0,w∈W,+∞,otherwise.


## 2.3. The Algorithm

In this subsection, we will introduce the algorithm for the proposed model. Since bothdxi,j,xi−1,j and dxi,j,xi−1,j are fixed in LiDAR observations, for simplicity, we use τi,j and μi,j to represent them, respectively, in the following discussion. Let pi,j=xi,j−xi,j−1,qi,j=xi,j−xi−1,j. We can rewrite (11) as follows:(13)x,w,p,q=argminx,w,p,q12∑i=1m∑j=1nwi,jyi,j−xi,j2+α∑i=1m∑j=1nwi,jpi,jτi,j+α∑i=1m∑j=1nwi,jqi,jμi,j+χWw+λ1∑i=1m∑j=1nexp−si,jwi,j2+λ2∑i=1m∑j=1nexpsi,jwi,j−12,s.t.pi,j=xi,j−xi,j−1,qi,j=xi,j−xi−1,j.The augmented Lagrangian of the above-constrained optimization problem is given as follows:(14)Lγ1,γ2x,w,p,q;Λ,Γ=12∑i=1m∑j=1nwi,jyi,j−xi,j2+α∑i=1m∑j=1nwi,jpi,jτi,j+α∑i=1m∑j=1nwi,jqi,jμi,j+χWw+λ1∑i=1m∑j=1nexp−si,jwi,j2+λ2∑i=1m∑j=1nexpsi,jwi,j−12+∑i=1m∑j=1nΛi,jpi,j−xi,j−xi,j−1+∑i=1m∑j=1nΓi,jqi,j−xi,j−xi−1,j+γ12∑i=1m∑j=1npi,j−xi,j−xi,j−12+γ22∑i=1m∑j=1nqi,j−xi,j−xi−1,j2,where γ1,γ2>0 are penalization parameters and ΛandΓ are Lagrange multipliers.With an initial guess ofx0,w0,p0,q0;Λ0,Γ0, the iterations of the alternating direction method of multipliers [26] are given as follows:(15)xk+1=argminxLγ1,γ2x,wk,pk,qk;Λk,Γk,wk+1=argminwLγ1,γ2xk+1,w,pk,qk;Λk,Γk,pk+1=argminpLγ1,γ2xk+1,wk+1,p,qk;Λk,Γk,qk+1=argminqLγ1,γ2xk+1,wk+1,pk+1,q;Λk,Γk,Λi,jk+1=Λi,jk+γ1pi,jk+1−xi,jk+1−xi,j−1k+1,Γi,jk+1=Γi,jk+γ2qi,jk+1−xi,jk+1−xi−1,jk+1.Next, we demonstrate how to solve the abovementioned subproblems (15) with respect to x,w,p,q.Let us consider thex-subproblem in (15):(16)xk+1=argminx12∑i=1m∑j=1nwi,jyi,j−xi,j2+γ12∑i=1m∑j=1npi,jk+Λi,jkγ1−xi,j−xi,j−12+γ22∑i=1m∑j=1nqi,jk+Γi,jkγ2−xi,j−xi−1,j2.Noting that the data-fitting term will vanish if there is no observation datayi,j. The optimality condition of x-subproblem is given by(17)wi,jk+2γ1+2γ2xi,j−γ1xi,j+1+xi,j−1−γ2xi−1,j+xi+1,j=wi,jkyi,j+γ1pi,jk−pi,j+1k+γ2qi,jk−qi+1,jk+Λi,jk−Λi,j+1k+Γi,jk−Γi+1,jk,ifyi,jexists,2γ1+γ2xi,j−γ1xi,j+1+xi,j−1−γ2xi−1,j+xi+1,j=γ1pi,jk−pi,j+1k+γ2qi,jk−qi+1,jk+Λi,jk−Λi,j+1k+Γi,jk−Γi+1,jk,otherwise.We can solve the above set of equations by using the conjugate gradient (CG) method [27, 28].Forw-subproblem, we have(18)wk+1=argminw≥012∑i=1m∑j=1nwi,jyi,j−xi,jk+12+α∑i=1m∑j=1nwi,jpi,jkτi,j+α∑i=1m∑j=1nwi,jqi,jkμi,j+λ1∑i=1m∑j=1nexp−si,jwi,j2+λ2∑i=1m∑j=1nexpsi,jwi,j−12.We note thatwi,j are decoupled, and we only need to deal with scalar optimization problems. More precisely, if yi,j exists, the scalar optimization problem is given by(19)wi,jk+1=argminw≥012wi,jyi,j−xi,jk+12+αwi,jpi,jkτi,j+αwi,jqi,jkμi,j+λ1exp−si,jwi,j2+λ2expsi,jwi,j−12.Ifyi,j does not exist, it is equal to(20)wi,jk+1=argminw≥0αwi,jpi,jkτi,j+αwi,jqi,jkμi,j+λ1exp−si,jwi,j2+λ2expsi,jwi,j−12.Both the above two scalar optimization problems are convex. And one can readily get that(21)wk+1=maxw¯k+1,0,where(22)w¯i,jk+1=2λ2expsi,j−1/2yi,j−xi,jk+12−αpi,jk/τi,j−αpi,jk/μi,j2λ1exp−si,j+2λ2expsi,j,ifyi,jexists,2λ2expsi,j−αpi,jk/τi,j−αqi,jk/μi,j2λ1exp−si,j+2λ2expsi,j,otherwise.Forp-subproblem and q-subproblem, they are given by(23)pk+1=argminpα∑i=1m∑j=1nwi,jpi,jτi,j+γ12∑i=1m∑j=1npi,j+Λi,jkγ1−xi,jk+1−xi,j−1k+12,and(24)qk+1=argminqα∑i=1m∑j=1nwi,jqi,jμi,j+γ22∑i=1m∑j=1nqi,j+Γi,jkγ2−xi,jk+1−xi−1,jk+12,respectively. These two subproblems can be solved by using the soft-thresholding technique [29], and their solutions are given by(25)pi,jk+1=maxxi,jk+1−xi,j−1k+1−Λi,jkγ1−αwi,jk+1γ1τi,j,0signxi,jk+1−xi,j−1k+1−Λi,jkγ1,qi,jk+1=maxxi,jk+1−xi−1,jk+1−Γi,jkγ2−αwi,jk+1γ2μi,j,0signxi,jk+1−xi−1,jk+1−Γi,jkγ2,where signp≔o/p.Finally, the overall algorithm is summarized in Algorithm1.Algorithm 1:
Input: Noisy data y, SNR data s, and parameters α,λ1,λ2,γ1,γ2.Initialization:x0=y,w0=1,p0=0,q0=0,Λ0=0,Γ0=0,k=0.(1)
while Stopping criteria is not satisfied do(2)
Solvexk+1 by using CG for (17)(3)
Solvewk+1 by (21)(4)
Findpk+1 by (24)(5)
Findqk+1 by (25)(6)
Update the multipliers byΛi,jk+1=Λi,jk+γ1pi,jk+1−xi,jk+1−xi,j−1k+1 and Γi,jk+1=Γi,jk+γ2qi,jk+1−xi,jk+1−xi−1,jk+1(7)
k←k+1(8)
end while


## 2.4. The Calculation of Parameters

In the proposed denoising scheme, there are several parameters to be determined. In this subsection, we incorporate the L-curve method to determine the values of parameters. The L-curve is a tradeoff curve between two quantities that both need to be controlled and balanced [30, 31]. It is widely used in the engineering and applied mathematics field.First, the regularization parameterα in (14) is crucial in our proposed method since it controls the ratio between the data-fitting term and the total variation regularization term. When the value of α is large (small), the importance of the total variation regularization term is large (small) and the denoising results are more (less) regularized. Moreover, there are two parameters λ1 and λ2 in (14), and they can affect the values of wi,j according to the SNR values of LiDAR observations. When the ratio λ1/λ2 is large (i.e., λ1 is larger than λ2), the term λ1∑i=1m∑j=1mexp−si,jwi,j2 is more important than the term λ2∑i=1m∑j=1mexp−si,jwi,j−12. It follows that when LiDAR observations have small SNR values, there would be more zero values of wi,j in the optimization procedure. In contrast, when the ratio λ1/λ2 is small and LiDAR observations have large SNR values, there would be more values of wi,j to be one in the optimization procedure.In the L-curve method, we consider the following two quantities:(26)Dxα,λ1,λ2,wα,λ1,λ2≔−log10∑i=1m∑j=1nwi,jα,λ1,λ2yi,j−xi,jα,λ1,λ22∑i=1m∑j=1nwi,jα,λ1,λ2,and(27)Sxα,λ1,λ2,wα,λ1,λ2≔−log10∑i=1m∑j=1nwi,jα,λ1,λ2TVxi,jα,λ1,λ2∑i=1m∑j=1nwi,jα,λ1,λ2.The first quantityDxα,λ1,λ2,wα,λ1,λ2 refers to the degree of data fitting when x and w are computed by Algorithm 1 for fixed α, λ1, and λ2. Similarly, the second quantity Sxα,λ1,λ2,wα,λ1,λ2 refers to the degree of smoothness in the regularization for the computed x and w. In order to compare the degrees of data fitting and smoothness with respect to different numbers of denoising values wi,j across different values of λ1 and λ2, the weighted average of data fitting and smoothness are accounted in the two quantities. The main idea of the L-curve method is to select the suitable values of α, λ1, and λ2 such that both Dxα,λ1,λ2,wα,λ1,λ2 and Sxα,λ1,λ2,wα,λ1,λ2 are balanced.As an illustration, we apply Algorithm1 to the LiDAR observations in Figure 1 for different values of α, λ1, and λ2. More precisely, in Figure 2(a), we fix λ1 (100) and λ2 (0.5,5,50,500) and compute x and w for different values of α (0.01,0.05,0.1,0.3,0.5,0.8,1,5). In Figure 2(a), we plot Dxα,λ1,λ2,wα,λ1,λ2 and Sxα,λ1,λ2,wα,λ1,λ2. The lines of points of D and S with respect to α are generated for each fixed λ1 and λ2. In total, there are four lines of points in Figure 2(a). We observe in each line that when α is small, D is large and S is small, and thus the points appear in the left hand side of the line. When α is large, D is small and S is large, and thus the points appear in the right hand side of the line. We find that there are several corner points in the line, and they refer to several large rates of change of data fitting with respect to the change of smoothness. Here, we can pick up the corner point with the largest rate of change and employ the corresponding value of α for regularization. Next, we select the values of λ1 and λ2 by comparing the selected corner points of different lines. Here, we can pick up the selected corner point with the largest D (the data fitting is good) and the largest S (the smoothness is large). According to the plot in Figure 2(a), we select α∗=0.5, λ1∗=100,and λ2∗=500.Figure 2
L-curves of different parametersα, λ1, and λ2.
(a)(b)Based on the above idea, we summarize our L-curve method with automatic selection ofα, λ1, and λ2 in Algorithm 2.Algorithm 2: L-curve Method.
Input: The set of parameters values α, λ1, λ2.Initialization: Corner point c0=0,0(1)
while the set of parameter values is not empty do(2)
Select (λ1, λ2) from the parameter set(3)
Find the corner pointc∗ of the generated lines of points Dxα,λ1,λ2,wα,λ1,λ2,Sxα,λ1,λ2,wα,λ1,λ2 from a set of α values, and record the corresponding α˜.(4)
If the position ofc∗ is higher than the position of c0 in the D-S plot, set c0=c∗, α∗=α˜, λ1∗=λ1, λ2∗=λ2(5)
end whileFrom Figure2, it is obvious that the model is very sensitive for different values of the regularization parameter α due to its balance effect for data fitting and smoothing. However, from Figure 2(a), for fixed values of parameters α and λ1, the changes of D and S for different values of λ2 are small (less than 1), which means that the model is pretty robust for parameter λ2. To evaluate the robustness of parameter λ1, we plot the four L-curves with respect to α (0.01,0.05,0.1,0.3,0.5,0.8,1,5) for fixed λ2 in Figure 2(b), where λ2=100, λ1 varies in 0.5,5,50,500 for different curves. Similar to the results of parameter λ2, one can readily find that the proposed model is quite robust for parameter λ1.


## 3. Results and Discussion

In this section, to illustrate the effectiveness of the proposed scheme (denoted as the “TV” method), we will show some denoising results by this scheme and compare them with the results generated by other denoising methods for the LiDAR data collected at HKIA. Clearly, the LiDAR data we study in this paper are time-varying and distance-varying, so that we compare the results by the DWT method [9] and the EMD-CIIT method [16]. They are denoted as “DWT” and “EMD-CIIT.” Moreover, we consider the method used by Baranov et al. [18] (denoted as “Baranov” in the rest of this article) and Newsom et al. [19] (denoted as “Newsom” in the rest of this article) as well as the method used by Hong Kong Observatory (denoted as the HKO method in the rest of this article) (see Algorithm 3 for more detailed description). Since Baranov et al. did not give an exact SNR threshold level in [18], we set the SNR threshold level as −5, which is used by Hong Kong Observatory for the LiDAR data observed at HKIA. The effects of averaging smoothing algorithm and median smoothing algorithm are shown to be about the same in [18], so that we apply a five-point-based averaging smoothing algorithm to the data selected by the SNR threshold level. Refer to Algorithm 4 for more detailed description.Algorithm 3: The Denoising Method used by Hong Kong Observatory.
Input: LiDAR observations yi,j and their SNR values si,j for 1≤i≤m and 1≤j≤n(1)
i=1 (for each azimuth angle)(2)
forj=1 to n (along each range) do(3)
Ifsi,j<=−5, then yi,j=NA(4)
end for(5)
xi,1=yi,1(6)
forj=2 to ndo(7)
Ifxi,j−1HKO−yi,j>16, then xi,jHKO=xi,j−1HKO×0.95+yi,j×0.05(8)
Ifxi,j−1HKO=NA and yi,j+1=NA, then xi,jHKO=NA(9)
end for(10)
Seti=i+1 and Goto Step 2Algorithm 4: The Denoising Method used by Baranov et al.
Input: LiDAR observations yi,j and their SNR values si,j for 1≤i≤m and 1≤j≤n(1)
Seti=1 (for each azimuth angle)(2)
forj=1 to n (along each range) do(3)
Ifsi,j<=−5, then yi,j=NA(4)
end for(5)
xi,1Baranov=yi,1(6)
xi,2Baranov=yi,1+yi,2/2(7)
forj=3 to n−2do(8)
xi,jBaranov=yi,j−2+yi,j−1+yi,j+yi,j+1+yi,j+2/5(9)
end for(10)
xi,n−1Baranov=yi,n−3+yi,n−2+yi,n−1+yi,n/4(11)
xi,nBaranov=yi,n−2+yi,n−1+yi,n/3(12)
Seti=i+1 and Goto Step 2


### 3.1. Comparison Results with the “DWT” and “EMD-CIIT” Methods

Figures3 and 4 show the results of the “DWT” method, the “EMD-CIIT” method, and the proposed TV method for LiDAR data collected on 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 with a fixed azimuth angle 224.0° and 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39 with a fixed azimuth angle 314.0°, respectively. According to the record from Hong Kong International Airport, these two wind profiles from LiDAR data correspond to two windshear cases. According to Figures 3 and 4, one can readily find that(1)
The “DWT” method and the “EMD-CIIT” method cannot remove very noisy observations, but our proposed scheme can remove noisy observations with very low SNR values effectively.For example, there is a high peak appearing at the range of around 8200 m in the LiDAR data observation in Figure3. The DWT denoising method (xDWT) and the EMD-CIIT denoising method (xEMD−CIIT) cannot remove this very noisy observation. We see from Figure 3(d) that the corresponding SNR value is about −8 at the range of around 8200 m. Indeed, it is quite acceptable to remove such noisy observation compared with neighborhood observations. Similarly, in Figure 4, there is a rock bottom appearing at the range of around 5000 m in the LiDAR data observation. From Figure 4(d), one can readily get that the corresponding SNR value is about −7 at the range of around 5000 m, which indicates the high noise intensity at this slant range. We see from Figures 4(a) and 4(b) that the “DWT” method and the “EMD-CIIT” method cannot remove such noisy observations. However, the proposed scheme (Figures 3 and 4(c)), according to the corresponding SNR values, filters both the high peak and rock bottom as well as the noisy data around them out.(2)
The method can over-smooth the input LiDAR data, and some LiDAR observations may be lost. As comparison, the denoising results of the proposed scheme are smooth but not distorted.Figure 3
The denoising results by the “DWT” (a), “EMD-CIIT” (b), and our proposed scheme (c) for the LiDAR data collected on 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 with a fixed azimuth angle (224.0°). The SNR values of the corresponding observed LiDAR data are shown in (d).
(a)(b)(c)(d)Figure 4
The denoising results by the “DWT” (a), “EMD-CIIT” (b), and our proposed scheme (c) for the LiDAR data collected on 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39 with a fixed azimuth angle (314.0°). The SNR values of the corresponding observed LiDAR data are shown in the (d).
(a)(b)(c)(d)


### 3.2. Comparison Results with the “Baranov” Method

Figures5(a) and 5(c) show the results by the “Baranov” method and the proposed scheme for LiDAR data collected on 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39 with a fixed azimuth angle 314° and 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 with a fixed azimuth angle 277°, respectively. According to Figures 5(a) and 5(c), we can easily find that(1
)The Baranov method might restore few missing points as well as some data points that are removed by the SNR threshold level by the averaging smooth algorithm. However, the values of them seem to be quite the same as the previous observed one, which play a tiny role in the following data applications such like windshear detection.For example, in Figure5(a), there is a rock bottom appearing at the range of around 5000 m in the LiDAR data observation. From Figure 5(b), the corresponding SNR value is about −7, which indicates the data point is removed by the SNR threshold level. However, this data point is restored by the averaging smoothing algorithm in the results by the “Baranov” method (xBaranov). It is obvious that the velocity values of data points between the range of 4800 m and 5800 m are quite the same in xBaranov. Similarly, there are some missing values at the range of around 9000 m and 10800 m. Some missing data points are restored with almost same values in xBaranov.(2)
Due to the effect of the averaging smoothing algorithm, the “Baranov” method might mis-shift the velocities of some data points.For example, there are some mis-shifted data points in the result by the “Baranov” method at the range of around 4500 m in Figure5(a). From 5(b), the SNR values around 4500 m are around 0, which are still acceptable. Similarly, the SNR values of data points at the range of around 5500 m are also acceptable, but the velocity values are mis-shifted.(3)
In comparison, the results by the proposed method (xTV) remove the bad observations properly and provide less useless information.Figure 5
The denoising results by the proposed methodxTV and the “Baranov” method xBaranov. (a) The LiDAR data were on 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39 with a fixed azimuth angle 314° and the denoising results. (b) The SNR values of the LiDAR data in the 1st row. (c) The LiDAR data were on 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 with a fixed azimuth angle 277° and the denoising results. (d) The SNR values of the LiDAR data in the 3rd row.
(a)(b)(c)(d)


### 3.3. Comparison Results with HKO Method

In Figures6(a) and 6(c), we consider LiDAR data collected on 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39 with a fixed azimuth angle 124° and 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 with a fixed azimuth angle 139°, respectively. The denoising results xHKO of LiDAR data by the HKO method are shown in Figures 6(a) and 6(c). We have the following observations.(1)
There are some mis-shifted results at the range values around 4000m-4200 m in Figure6(a) and around 4700m-5000 m in Figure 6(c). For reference, the SNR values corresponding to LiDAR data in Figures 6(a) and 6(c) are given in Figures 6(b) and 6(d), respectively. We note in Figure 6(b) that the SNR values at 4000m-4200 m are 2–10. Similarly, the SNR values at 4700 m−5000 m are 13–18 in Figure 6(d). In these two cases, SNR values can still be in an acceptable range. However, the method by Hong Kong observatory cannot provide a reasonable denoising results in these two range values, and the detected wind profiles may be affected.(2)
In comparison, we show the denoising resultsxTV of LiDAR data by the proposed method in Figures 6(a) and 6(c). We find that the proposed method fits the LiDAR data quite well with as much smoothness as possible based on the SNR values of the given observations.Figure 6
The denoising results by the proposed methodxTV and the HKO method xHKO. (a) The LiDAR data were collected on 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39 with a fixed azimuth angle 124° and the denoising results. (b) The SNR values of the LiDAR data in the 1st row. (c) The LiDAR data were collected on 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 with a fixed azimuth angle 139° and the denoising results. (d) The SNR values of the LiDAR data in the 3rd row.
(a)(b)(c)(d)


### 3.4. Comparison Results with the “Newsom” Method

Note that the “Newsom” method just removes the bad observations based on the SNR values (the SNR threshold level is set to be 0.02 in [19]) without any other operations for LiDAR data, so that there are no value differences between the data after processing and the observed data. For LiDAR data in Figure 3, all the observed data after the range 4400 m are removed because their SNR values are less than 0.02. It seems that it is not reasonable to remove all these observed data points. In Figure 4, the observed data points near the rock bottom are removed. In Figures 5 and 6, the processed data by the “Newsom” method are about the same as the original observed data.


### 3.5. Comparison Results in Conical Scans

In Figures7 and 8, we plot the conical scans of LiDAR data collected on 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39 and 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 and the results by the “DWT” method, the “EMD-CIIT” method, the“Newsom” method, the “Baranov” method, the HKO method, and the proposed method. The computational time for these methods is also given. From these figures, we can readily find that(1)
There are some spikes in the denoising results of the “DWT” method and the “EMD-CIIT” method. For example, in Figures7(e) and 7(f), there are some spikes in the range of azimuth angles from 206° to 250°.(2)
The over-smoothness of the “EMD-CIIT” method is also overcome by the proposed scheme.(3)
There are some outliers that provide little information of the wind profile in the denoising results generated by the “Baranov” method. In comparison, the denoising results of our proposed method are much clear.We note that Figures7 and 8 correspond to two windshear cases in March 2015. Although the existence of outliers does not have an effect on visual judgement for windshear analysis, it can have a strong effect on machine learning methods in windshear detection. For instance, the wind profile features suggested by the machine learning method in [32] are calculated based on the maximal difference in wind velocities within a range of azimuth angles, and it is clear that the outliers can change the resulting wind profile features.(4)
Since the SNR threshold level (SNR>0.02) used by “Newsom” method is high, there are only few data points retined in the results. And the proposed method, without any fixed SNR threshold level, retained more data points which contain useful information of wind profiles.(5)
Comparing with the results by HKO method, the results of our proposed scheme is much smoother.(6)
The proposed scheme is faster than the EMD-CIIT method, but it is slower than the other methods.Figure 7
(a) The LiDAR data and (b) SNR values collected on 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39. The denoising results by (c) the “DWT” method, (d) the “EMD-CIIT” method, (e) the “Newsom” method, (f) the “Baranov” method, (g) the HKO method, and (h) the proposed scheme (denoted by TV).
(a)(b)(c)(d)(e)(f)(g)(h)Figure 8
(a) The LiDAR data and (b) SNR values collected on 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25. The denoising results by (c) the “DWT” method, (d) the “EMD-CIIT” method, (e) the “Newsom” method, (f) the “Baranov” method, (g) the HKO method, and (h) the proposed scheme (denoted by TV).
(a)(b)(c)(d)(e)(f)(g)(h)For example, there are some fluctuations in the range of azimuth angles from257° to 359° in Figure 8(g), but the data are much smoother in Figure 8(h).To further evaluate the performance of the proposed method, we show a nonwindshear case in Figure9. In the figure, we plot the conical scans of the LiDAR data and the denoising results by the abovementioned six methods. The LiDAR data were collected on 4 March 2015 from 00 : 00 : 16 to 00 : 00 : 39. Similarly, the denoising results given by the proposed scheme do not contain any outliers and retain more data points with good smooth effect.Figure 9
(a) The LiDAR data and (b) SNR values collected on 4 March 2015 from 00 : 00 : 16 to 00 : 00 : 39. The denoising results by (c) “DWT” method, (d) “EMD-CIIT” method, (e) “Newsom” method, (f) “Baranov” method, (g) HKO method, and (h) the proposed scheme (denoted by TV).
(a)(b)(c)(d)(e)(f)(g)(h)


### 3.6. Comparison Results in Noisy Observations Removal

Finally, we would like to demonstrate the capability of the proposed method in removing noisy observations compared with the “Newsom” method,the “Baranov” method, as well as the HKO method. As a matter of fact, when SNR values are too small, the LiDAR observations may not be suitable for data analysis usage. In Table1, we show the noisy observation removal results for a set of LiDAR observations for windshear cases at Hong Kong International Airport. The set of observations was collected from 1 March to 31 March 2015. There were 147 windshear cases reported in the pilot report. Each scan of windshear cases took about 25 seconds (see the two examples in Figures 7 and 8). We see from Table 1 that(1)
Since the SNR threshold level of the“Newsom” method is 0.02, all the noisy observations are removed in the results of the “Newsom” method, but the number of the retained data points are really small.(2)
The number∑i=1m∑j=1nwi,j of denoising values by the proposed method is about the same as the number of original observations yi,j and the number of denoising values xi,jBaranov by the “Baranov” method and the number of denoising values xi,jHKO by HKO method in the slant range values from 359 m to 3509 m. Note that the SNR values of LiDAR observations in the slant range values from 359 m to 3509 m are usually high, which means LiDAR observations are quite accurate.(3)
There are usually more noisy observations in large slant range values. We see from Table1 that there are only 281927 LiDAR observations in the slant ranges from 6659 m to 9809 m and 793282 LiDAR observations in the slant ranges from 3509 m to 6659 m compared with 1033268 LiDAR observations in the slant ranges from 359 m to 3509 m.(4)
The number∑i=1m∑j=1nwi,j of denoising values by the proposed method is less than the number of original observations yi,j and the number of denoising values xi,jBaranov by the “Baranov” method and the number of denoising values xi,jHKO by the HKO method when there are many noisy observations in the slant ranges from 6659 m to 9809 m. However, the proposed method preserves more data points with high SNR values in the slant ranges from 6659 m to 9809 m.Table 1
The comparison of noisy observations removal by the “Newsom” method, “Baranov” method, HKO method, and proposed method.
RangeDataSNR>=−5SNR>=−3SNR>=0359 m–3509 m (total number: 1036067)LiDAR observations103326810328851031484Denoising values by the “Newsom” method103147210314721031472Denoising values by the “Baranov” method103326810328851031484Denoising values by the HKO method103321210328741031484Denosing values by the proposed method1032987103279310313943509 m–6659 m (total number: 836510)LiDAR observations793282772883678138Denoising values by the “Newsom” method677112677112677112Denoising values by the “Baranov” method793282772883678138Denoising values by the HKO method791440772670678130Denosing values by the proposed method7867727726326780436659 m–9809 m (total number: 455697)LiDAR observations28192718442577073Denoising values by the “Newsom” method766457664576645Denoising values by the “Baranov” method28192718442577073Denoising values by the HKO method26303318302677018Denoising values by the proposed method24771318435377070


## 3.1. Comparison Results with the “DWT” and “EMD-CIIT” Methods

Figures3 and 4 show the results of the “DWT” method, the “EMD-CIIT” method, and the proposed TV method for LiDAR data collected on 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 with a fixed azimuth angle 224.0° and 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39 with a fixed azimuth angle 314.0°, respectively. According to the record from Hong Kong International Airport, these two wind profiles from LiDAR data correspond to two windshear cases. According to Figures 3 and 4, one can readily find that(1)
The “DWT” method and the “EMD-CIIT” method cannot remove very noisy observations, but our proposed scheme can remove noisy observations with very low SNR values effectively.For example, there is a high peak appearing at the range of around 8200 m in the LiDAR data observation in Figure3. The DWT denoising method (xDWT) and the EMD-CIIT denoising method (xEMD−CIIT) cannot remove this very noisy observation. We see from Figure 3(d) that the corresponding SNR value is about −8 at the range of around 8200 m. Indeed, it is quite acceptable to remove such noisy observation compared with neighborhood observations. Similarly, in Figure 4, there is a rock bottom appearing at the range of around 5000 m in the LiDAR data observation. From Figure 4(d), one can readily get that the corresponding SNR value is about −7 at the range of around 5000 m, which indicates the high noise intensity at this slant range. We see from Figures 4(a) and 4(b) that the “DWT” method and the “EMD-CIIT” method cannot remove such noisy observations. However, the proposed scheme (Figures 3 and 4(c)), according to the corresponding SNR values, filters both the high peak and rock bottom as well as the noisy data around them out.(2)
The method can over-smooth the input LiDAR data, and some LiDAR observations may be lost. As comparison, the denoising results of the proposed scheme are smooth but not distorted.Figure 3
The denoising results by the “DWT” (a), “EMD-CIIT” (b), and our proposed scheme (c) for the LiDAR data collected on 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 with a fixed azimuth angle (224.0°). The SNR values of the corresponding observed LiDAR data are shown in (d).
(a)(b)(c)(d)Figure 4
The denoising results by the “DWT” (a), “EMD-CIIT” (b), and our proposed scheme (c) for the LiDAR data collected on 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39 with a fixed azimuth angle (314.0°). The SNR values of the corresponding observed LiDAR data are shown in the (d).
(a)(b)(c)(d)


## 3.2. Comparison Results with the “Baranov” Method

Figures5(a) and 5(c) show the results by the “Baranov” method and the proposed scheme for LiDAR data collected on 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39 with a fixed azimuth angle 314° and 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 with a fixed azimuth angle 277°, respectively. According to Figures 5(a) and 5(c), we can easily find that(1
)The Baranov method might restore few missing points as well as some data points that are removed by the SNR threshold level by the averaging smooth algorithm. However, the values of them seem to be quite the same as the previous observed one, which play a tiny role in the following data applications such like windshear detection.For example, in Figure5(a), there is a rock bottom appearing at the range of around 5000 m in the LiDAR data observation. From Figure 5(b), the corresponding SNR value is about −7, which indicates the data point is removed by the SNR threshold level. However, this data point is restored by the averaging smoothing algorithm in the results by the “Baranov” method (xBaranov). It is obvious that the velocity values of data points between the range of 4800 m and 5800 m are quite the same in xBaranov. Similarly, there are some missing values at the range of around 9000 m and 10800 m. Some missing data points are restored with almost same values in xBaranov.(2)
Due to the effect of the averaging smoothing algorithm, the “Baranov” method might mis-shift the velocities of some data points.For example, there are some mis-shifted data points in the result by the “Baranov” method at the range of around 4500 m in Figure5(a). From 5(b), the SNR values around 4500 m are around 0, which are still acceptable. Similarly, the SNR values of data points at the range of around 5500 m are also acceptable, but the velocity values are mis-shifted.(3)
In comparison, the results by the proposed method (xTV) remove the bad observations properly and provide less useless information.Figure 5
The denoising results by the proposed methodxTV and the “Baranov” method xBaranov. (a) The LiDAR data were on 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39 with a fixed azimuth angle 314° and the denoising results. (b) The SNR values of the LiDAR data in the 1st row. (c) The LiDAR data were on 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 with a fixed azimuth angle 277° and the denoising results. (d) The SNR values of the LiDAR data in the 3rd row.
(a)(b)(c)(d)


## 3.3. Comparison Results with HKO Method

In Figures6(a) and 6(c), we consider LiDAR data collected on 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39 with a fixed azimuth angle 124° and 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 with a fixed azimuth angle 139°, respectively. The denoising results xHKO of LiDAR data by the HKO method are shown in Figures 6(a) and 6(c). We have the following observations.(1)
There are some mis-shifted results at the range values around 4000m-4200 m in Figure6(a) and around 4700m-5000 m in Figure 6(c). For reference, the SNR values corresponding to LiDAR data in Figures 6(a) and 6(c) are given in Figures 6(b) and 6(d), respectively. We note in Figure 6(b) that the SNR values at 4000m-4200 m are 2–10. Similarly, the SNR values at 4700 m−5000 m are 13–18 in Figure 6(d). In these two cases, SNR values can still be in an acceptable range. However, the method by Hong Kong observatory cannot provide a reasonable denoising results in these two range values, and the detected wind profiles may be affected.(2)
In comparison, we show the denoising resultsxTV of LiDAR data by the proposed method in Figures 6(a) and 6(c). We find that the proposed method fits the LiDAR data quite well with as much smoothness as possible based on the SNR values of the given observations.Figure 6
The denoising results by the proposed methodxTV and the HKO method xHKO. (a) The LiDAR data were collected on 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39 with a fixed azimuth angle 124° and the denoising results. (b) The SNR values of the LiDAR data in the 1st row. (c) The LiDAR data were collected on 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 with a fixed azimuth angle 139° and the denoising results. (d) The SNR values of the LiDAR data in the 3rd row.
(a)(b)(c)(d)


## 3.4. Comparison Results with the “Newsom” Method

Note that the “Newsom” method just removes the bad observations based on the SNR values (the SNR threshold level is set to be 0.02 in [19]) without any other operations for LiDAR data, so that there are no value differences between the data after processing and the observed data. For LiDAR data in Figure 3, all the observed data after the range 4400 m are removed because their SNR values are less than 0.02. It seems that it is not reasonable to remove all these observed data points. In Figure 4, the observed data points near the rock bottom are removed. In Figures 5 and 6, the processed data by the “Newsom” method are about the same as the original observed data.


## 3.5. Comparison Results in Conical Scans

In Figures7 and 8, we plot the conical scans of LiDAR data collected on 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39 and 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 and the results by the “DWT” method, the “EMD-CIIT” method, the“Newsom” method, the “Baranov” method, the HKO method, and the proposed method. The computational time for these methods is also given. From these figures, we can readily find that(1)
There are some spikes in the denoising results of the “DWT” method and the “EMD-CIIT” method. For example, in Figures7(e) and 7(f), there are some spikes in the range of azimuth angles from 206° to 250°.(2)
The over-smoothness of the “EMD-CIIT” method is also overcome by the proposed scheme.(3)
There are some outliers that provide little information of the wind profile in the denoising results generated by the “Baranov” method. In comparison, the denoising results of our proposed method are much clear.We note that Figures7 and 8 correspond to two windshear cases in March 2015. Although the existence of outliers does not have an effect on visual judgement for windshear analysis, it can have a strong effect on machine learning methods in windshear detection. For instance, the wind profile features suggested by the machine learning method in [32] are calculated based on the maximal difference in wind velocities within a range of azimuth angles, and it is clear that the outliers can change the resulting wind profile features.(4)
Since the SNR threshold level (SNR>0.02) used by “Newsom” method is high, there are only few data points retined in the results. And the proposed method, without any fixed SNR threshold level, retained more data points which contain useful information of wind profiles.(5)
Comparing with the results by HKO method, the results of our proposed scheme is much smoother.(6)
The proposed scheme is faster than the EMD-CIIT method, but it is slower than the other methods.Figure 7
(a) The LiDAR data and (b) SNR values collected on 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39. The denoising results by (c) the “DWT” method, (d) the “EMD-CIIT” method, (e) the “Newsom” method, (f) the “Baranov” method, (g) the HKO method, and (h) the proposed scheme (denoted by TV).
(a)(b)(c)(d)(e)(f)(g)(h)Figure 8
(a) The LiDAR data and (b) SNR values collected on 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25. The denoising results by (c) the “DWT” method, (d) the “EMD-CIIT” method, (e) the “Newsom” method, (f) the “Baranov” method, (g) the HKO method, and (h) the proposed scheme (denoted by TV).
(a)(b)(c)(d)(e)(f)(g)(h)For example, there are some fluctuations in the range of azimuth angles from257° to 359° in Figure 8(g), but the data are much smoother in Figure 8(h).To further evaluate the performance of the proposed method, we show a nonwindshear case in Figure9. In the figure, we plot the conical scans of the LiDAR data and the denoising results by the abovementioned six methods. The LiDAR data were collected on 4 March 2015 from 00 : 00 : 16 to 00 : 00 : 39. Similarly, the denoising results given by the proposed scheme do not contain any outliers and retain more data points with good smooth effect.Figure 9
(a) The LiDAR data and (b) SNR values collected on 4 March 2015 from 00 : 00 : 16 to 00 : 00 : 39. The denoising results by (c) “DWT” method, (d) “EMD-CIIT” method, (e) “Newsom” method, (f) “Baranov” method, (g) HKO method, and (h) the proposed scheme (denoted by TV).
(a)(b)(c)(d)(e)(f)(g)(h)


## 3.6. Comparison Results in Noisy Observations Removal

Finally, we would like to demonstrate the capability of the proposed method in removing noisy observations compared with the “Newsom” method,the “Baranov” method, as well as the HKO method. As a matter of fact, when SNR values are too small, the LiDAR observations may not be suitable for data analysis usage. In Table1, we show the noisy observation removal results for a set of LiDAR observations for windshear cases at Hong Kong International Airport. The set of observations was collected from 1 March to 31 March 2015. There were 147 windshear cases reported in the pilot report. Each scan of windshear cases took about 25 seconds (see the two examples in Figures 7 and 8). We see from Table 1 that(1)
Since the SNR threshold level of the“Newsom” method is 0.02, all the noisy observations are removed in the results of the “Newsom” method, but the number of the retained data points are really small.(2)
The number∑i=1m∑j=1nwi,j of denoising values by the proposed method is about the same as the number of original observations yi,j and the number of denoising values xi,jBaranov by the “Baranov” method and the number of denoising values xi,jHKO by HKO method in the slant range values from 359 m to 3509 m. Note that the SNR values of LiDAR observations in the slant range values from 359 m to 3509 m are usually high, which means LiDAR observations are quite accurate.(3)
There are usually more noisy observations in large slant range values. We see from Table1 that there are only 281927 LiDAR observations in the slant ranges from 6659 m to 9809 m and 793282 LiDAR observations in the slant ranges from 3509 m to 6659 m compared with 1033268 LiDAR observations in the slant ranges from 359 m to 3509 m.(4)
The number∑i=1m∑j=1nwi,j of denoising values by the proposed method is less than the number of original observations yi,j and the number of denoising values xi,jBaranov by the “Baranov” method and the number of denoising values xi,jHKO by the HKO method when there are many noisy observations in the slant ranges from 6659 m to 9809 m. However, the proposed method preserves more data points with high SNR values in the slant ranges from 6659 m to 9809 m.Table 1
The comparison of noisy observations removal by the “Newsom” method, “Baranov” method, HKO method, and proposed method.
RangeDataSNR>=−5SNR>=−3SNR>=0359 m–3509 m (total number: 1036067)LiDAR observations103326810328851031484Denoising values by the “Newsom” method103147210314721031472Denoising values by the “Baranov” method103326810328851031484Denoising values by the HKO method103321210328741031484Denosing values by the proposed method1032987103279310313943509 m–6659 m (total number: 836510)LiDAR observations793282772883678138Denoising values by the “Newsom” method677112677112677112Denoising values by the “Baranov” method793282772883678138Denoising values by the HKO method791440772670678130Denosing values by the proposed method7867727726326780436659 m–9809 m (total number: 455697)LiDAR observations28192718442577073Denoising values by the “Newsom” method766457664576645Denoising values by the “Baranov” method28192718442577073Denoising values by the HKO method26303318302677018Denoising values by the proposed method24771318435377070


## 4. Conclusion

In this article, we have proposed a LiDAR data denoising scheme for wind profile observations to improve data quality and remove bad observations. By introducing a weight in each observation, the scheme can filter out very noisy observations according to their SNR values and neighborhood observations. Combining with the data-fitting term and the polar-based total variation smoothing term, a global-based denoising model is built. The alternating direction method of multipliers is applied to find a solution of the proposed model. To find suitable regularization parameters of the proposed scheme, we consider an L-curve based parameter selection method that can select parameters automatically via the balance between the data-fitting term and the polar-based total variation regularization term. By applying the proposed scheme to the LiDAR data collected at the Hong Kong International Airport, we find that (i) our proposed scheme performs better than the denoising methods such as DWT and EMD-CIIT, where they cannot handle very noisy observations and they can conduct a denoising procedure along each slant range; (ii) our proposed method can handle noisy observations quite well, and its performance is better than that of HKO method and “Newsom” method as well as “Baranov” method; (iii) our proposed scheme can balance the data-fitting and the smoothing regularization via suitable parameter selection. However, the proposed scheme is slower than the comparing methods except the EMD-CIIT method. As a future research work, we will investigate a much faster algorithm for the proposed model. And we also need to consider an online denoising scheme such that it can handle LiDAR data sets in a continuous-time manner and study three-dimensional LiDAR data observations for denoising.



---

*Source: 3792225-2022-01-22.xml*