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.
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 Section 2 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 Section 2, 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 Subsection 2.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.
(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 took about 25 seconds (see one example in Figure 1).
2.2. The Proposed Model
Let be the LiDAR data observed at azimuth angle . For simplicity, we assume that there are m azimuth angles () to be recorded in between and , and also there are range values to be recorded, i.e.,
Note that the locations of such range values are not necessary to be uniform, and the distance between the LiDAR centre and the observed value is equal to . We are interested to compute the denoising LiDAR data as follows:according to the given SNR values:from the observed LiDAR data. When is a missing LiDAR observation, the corresponding can be set to be . In total, there are observations and SNR values and 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 data and the denoised data (). In order to determine whether at the azimuth angle and range value to be used in the model, we incorporate a nonnegative weight for . When the value of is equal to zero, the LiDAR observation is not used. The resulting data-fitting term is given by It is clear that when the LiDAR data observation is removed (), 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: where refers to the distance between the two observed range values at the same azimuth angle, i.e., and refers to the distance between the same observed range value at the two adjacent azimuth angles and , i.e., 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: We see from (7) and (8) that the denoised values 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 values to determine whether the LiDAR observation is used in the denoising procedure. Here, we propose the following term in the model:where and are the positive numbers to control the balance between the two terms in (9). We note in (9) that when is too small (a negative number), the LiDAR observation is very noisy, and therefore is large and is small. The optimization process drives to be zero. Similarly, when is not small (the LiDAR observation is not noisy), is small and is large and therefore the optimization process drives to be one. Finally, we also incorporate nonnegativity constraint on in the optimization model:
For the ease of presentation, we denote and . According to (4), (6), (9), and (10), the combined optimization model is given as follows:where is a positive number to balance the contribution of the total variation regularization term and the other terms and is the characteristic function of the nonnegativity constraint set
2.3. The Algorithm
In this subsection, we will introduce the algorithm for the proposed model. Since both and are fixed in LiDAR observations, for simplicity, we use and to represent them, respectively, in the following discussion. Let . We can rewrite (11) as follows:
The augmented Lagrangian of the above-constrained optimization problem is given as follows:where are penalization parameters and are Lagrange multipliers.
With an initial guess of , the iterations of the alternating direction method of multipliers [26] are given as follows:
Next, we demonstrate how to solve the abovementioned subproblems (15) with respect to .
Let us consider the -subproblem in (15):
Noting that the data-fitting term will vanish if there is no observation data . The optimality condition of -subproblem is given by
We can solve the above set of equations by using the conjugate gradient (CG) method [27, 28].
For -subproblem, we have
We note that are decoupled, and we only need to deal with scalar optimization problems. More precisely, if exists, the scalar optimization problem is given by
If does not exist, it is equal to
Both the above two scalar optimization problems are convex. And one can readily get thatwhere
For -subproblem and -subproblem, they are given byandrespectively. These two subproblems can be solved by using the soft-thresholding technique [29], and their solutions are given bywhere .
Finally, the overall algorithm is summarized in Algorithm 1.
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 and in (14), and they can affect the values of according to the SNR values of LiDAR observations. When the ratio is large (i.e., is larger than ), the term is more important than the term . It follows that when LiDAR observations have small SNR values, there would be more zero values of in the optimization procedure. In contrast, when the ratio is small and LiDAR observations have large SNR values, there would be more values of to be one in the optimization procedure.
In the L-curve method, we consider the following two quantities:and
The first quantity refers to the degree of data fitting when and are computed by Algorithm 1 for fixed , , and . Similarly, the second quantity refers to the degree of smoothness in the regularization for the computed and . In order to compare the degrees of data fitting and smoothness with respect to different numbers of denoising values across different values of and , 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 , , and such that both and are balanced.
As an illustration, we apply Algorithm 1 to the LiDAR observations in Figure 1 for different values of , , and . More precisely, in Figure 2(a), we fix (100) and () and compute and for different values of (). In Figure 2(a), we plot and . The lines of points of and with respect to are generated for each fixed and . In total, there are four lines of points in Figure 2(a). We observe in each line that when is small, is large and is small, and thus the points appear in the left hand side of the line. When is large, is small and 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 and by comparing the selected corner points of different lines. Here, we can pick up the selected corner point with the largest (the data fitting is good) and the largest (the smoothness is large). According to the plot in Figure 2(a), we select , ,and .
(a)
(b)
Based on the above idea, we summarize our L-curve method with automatic selection of , , and in Algorithm 2.
|
From Figure 2, 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 , the changes of and for different values of are small (less than 1), which means that the model is pretty robust for parameter . To evaluate the robustness of parameter , we plot the four L-curves with respect to () for fixed in Figure 2(b), where , varies in for different curves. Similar to the results of parameter , one can readily find that the proposed model is quite robust for parameter .
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 , 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.
|
| ||||||||||||||||||||||||||
3.1. Comparison Results with the “DWT” and “EMD-CIIT” Methods
Figures 3 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 and 4 March 2015 from 21 : 26 : 16 to 21 : 46 : 39 with a fixed azimuth angle , 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 Figure 3. The DWT denoising method () and the EMD-CIIT denoising method () 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 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
3.2. Comparison Results with the “Baranov” Method
Figures 5(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 and 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 with a fixed azimuth angle , 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 Figure 5(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 , 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 (). It is obvious that the velocity values of data points between the range of 4800 m and 5800 m are quite the same in . 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 .(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 Figure 5(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 () remove the bad observations properly and provide less useless information.
(a)
(b)
(c)
(d)
3.3. Comparison Results with HKO Method
In Figures 6(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 and 5 March 2015 from 04 : 16 : 02 to 04 : 16 : 25 with a fixed azimuth angle , respectively. The denoising results 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 Figure 6(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 results 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.
(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 Figures 7 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 Figures 7(e) and 7(f), there are some spikes in the range of azimuth angles from to .(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 Figures 7 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 ) 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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
For example, there are some fluctuations in the range of azimuth angles from to 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 Figure 9. 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.
(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 Table 1, 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 of denoising values by the proposed method is about the same as the number of original observations and the number of denoising values by the “Baranov” method and the number of denoising values 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 Table 1 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 of denoising values by the proposed method is less than the number of original observations and the number of denoising values by the “Baranov” method and the number of denoising values 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.
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.
Data Availability
The request for data used to support the findings of this study could be addressed to the Hong Kong Observatory. The data provision would be considered on a case-by-case basis.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The research was supported in part by HKRGC GRF 17201020 and 17300021 and UGC-RMGS 207300829.