Edge Detection Method for Determining Boundary Layer Height Based on Doppler Lidar

: The top of the boundary layer, referred to as the planetary boundary layer height (BLH), is an important physical parameter in atmospheric numerical models, which has a critical role in atmospheric simulation, air pollution prevention, and climate prediction. The traditional methods for determining BLHs using Doppler lidar vertical velocity variance ( σ 2 w ) can be classiﬁed into the variance and peak methods, which depend on atmospheric conditions due to their use of a single threshold, hence limiting their ability to estimate diurnal BLHs. Edge detection (ED) was later introduced in BLH estimation due to its ability to identify the 2D gradient of an image. A key step in ED is automatically identifying the edge of BLHs based on the peaks of the proﬁle, hence avoiding the inﬂuence of extreme atmospheric conditions. Two cases in the diurnal cycle on 4 March 2019 and 8 July 2019 reveal that ED outperforms both the variance and peak methods in nighttime and extreme atmospheric conditions. The retrieved BLHs from 2018 to 2020 were compared with radiosonde (RS) measurements for the same time at the neutral, stable, and convective boundary layers. The correlation coefﬁcient (R: 0.4 vs. 0.05, 0.14; 0.26 vs. − 0.10, − 0.16; 0.35 vs. 0.01, 0.16) and root mean square error (RMSE (km): 0.58 vs. 0.82, 0.90; 0.37 vs. 1.01, 0.50; 0.66 vs. 0.98, 0.82) obtained by the ED method were higher and lower than those obtained by the variance and peak methods, respectively. The mean absolute error (MAE) of the ED method under the NBL, SBL, and CBL conditions are lower than the variance and peak methods (MAE (km): 0.44, 0.14, 0.50 vs. 0.62, 0.34, 0.64; 0.59, 0.75, 0.74), respectively. The mean relative error (MRE) of the ED method is lower than the variance and peak methods under the NBL condition (MRE: − 8.88% vs. − 18.39%, 13.91%). Under the SBL, the MRE of the ED method is lower than the variance method and higher than the peak method ( − 38.64%, vs. − 152.23%; 14.02%). Under the CBL, the MRE of the ED method is lower than the variance method and higher than the peak method ( − 15.07% vs. 2.24%; 5.64%). In addition, the comparison between ED and wavelet covariance transform (WCT) method and RS measurements showed that the ED method has a similar performance with the WCT method and is even better. In the long-term analysis, the hourly and monthly BLHs in the diurnal and annual cycles, respectively, as obtained by ED, were highly consistent with the RS measurements and obtained the lowest standard error. In the annual cycle, the retrieved BLHs in summer and autumn were higher than those retrieved in spring and winter. and the results show that ED was more accurate than the variance and peak methods under extreme and routine atmospheric conditions. Therefore, ED can estimate the BLHs in the diurnal cycle and perform better than the other two methods. The BLHs retrieved using the variance, peak, WCT, and ED methods from 2018 to 2020 were compared with the RS measurements under NBL, SBL, and CBL. The results show that ED was more consistent and had lower RMSE with the RS measurements than the variance and peak methods under NBL, SBL, and CBL. The lower MAE and MRE indicated the results of the ED method with reliability, and a number of the BLHs derived from the ED method are closer to the reference line than those obtained by the variance and peak methods. In addition, the ED method has similar results with the WCT method and is even better. In the long-term analysis, the hourly BLHs retrieved by ED were consistent with theory, whereas those retrieved using the variance and peak methods were not. The monthly BLHs retrieved using the peak and ED methods were more consistent with the RS measurements than those obtained using the variance method. The BLHs recorded in summer and autumn were also higher than those in spring and winter. In general, the ED method can estimate the BLHs in the diurnal cycle under extreme and routine atmospheric conditions and has a better performance than the variance, peak, and WCT methods. But the performance of these methods are unsatisfactory at present. Direct detection of turbulence changes can obtain accurate BLH theoretically, but the low SNR of DL may be one of the reasons that reduces the accuracy of the BLHs obtained from this remote sensing equipment. The poor performance of the ED method motivates us to improve or/and develop more practical approaches moving forward.


Introduction
The planetary boundary layer (PBL) is an atmospheric layer located 100 m to 3000 m above ground [1]. PBL is easily influenced by air pollution, meteorological factors, atmospheric flow, and climate [2,3], and its structure can be divided into the convective boundary layer (CBL), stable boundary layer (SBL), and neutral boundary layer (NBL)

Data
The following sections present the σ 2 w , signal-to-noise ratio (SNR), noise of DL wind and PBL types (NBL, SBL, and CBL), and the BLHs retrieved using RS at the SGP station of the ARM facility.
2.1.1. Vertical Velocity Variance of Doppler Lidar Wind (σ 2 w ) The ARM facility operates DL systems around the globe. The instruments of DL employ an eye-safe solid-state laser transmitter that operates at 1.548 µm. all fiber-coupled components were incorporated in these instruments, and a master oscillator power amplifier architecture was used to generate the coherent light with low pulse energy (<100 µJ) and high pulse repetition frequency (15 kHz). DL provides measurements of radial velocity, attenuated aerosol backscatter, and SNR. DL systems typically operate on a fixed scan schedule, performs plan-position-indicator scans several times an hour, and conducts vertical pointing at other times of the day. They measure clear air vertical velocity profiles in the lower troposphere with a temporal resolution of about 1 s and a range resolution of 30 m [14]. The σ 2 w measurements in clear air are captured below the lowest cloud base, with the range resolution of 30 m and the temporal intervals of 10 min. The corresponding noise and median SNR values are provided to control the quality of data. The corrected variances should be filtered based on SNR, noise information, or both [19].
In this study, SNR and noise were set to 0.0072 and 1, respectively, to filter the σ 2 w from the SGP station [21]. The data were captured in local time and above the mean sea level. The geographical location of the station (36 • 36 26.3592" N, 97 • 29 15.5148" W) is shown in Figure 1.

Data
The following sections present the 2 w σ , signal-to-noise ratio (SNR), noise of DL w and PBL types (NBL, SBL, and CBL), and the BLHs retrieved using RS at the SGP sta of the ARM facility.

Vertical Velocity Variance of Doppler Lidar Wind (
The ARM facility operates DL systems around the globe. The instruments of DL ploy an eye-safe solid-state laser transmitter that operates at 1.548 µm. all fiber-coup components were incorporated in these instruments, and a master oscillator power plifier architecture was used to generate the coherent light with low pulse energy (< µJ) and high pulse repetition frequency (15 kHz). DL provides measurements of ra velocity, attenuated aerosol backscatter, and SNR. DL systems typically operate on a fi scan schedule, performs plan-position-indicator scans several times an hour, and c ducts vertical pointing at other times of the day. They measure clear air vertical velo profiles in the lower troposphere with a temporal resolution of about 1 s and a range olution of 30 m [14]. The In this study, SNR and noise were set to 0.0072 and 1, respectively, to filter the from the SGP station [21]. The data were captured in local time and above the mean level. The geographical location of the station (36°36′26.3592″ N, 97°29′15.5148″ W shown in Figure 1.   [27]. The ARM facility provides the PBL types and BLHs for RS measurements. The PBL types were computed following the method of Liu and Liang [4], whereas the BLHs were derived using the bulk Richardson (BRN) method. The PBL types can be classified into the three major regimes, namely, CBL, NBL, and SBL [1,4]. Near-surface potential temperature gradients were used to identify the boundary layer regime. The following criteria were used to identify these regimes: where θ is potential temperature, and its subscript number denotes for the data level index assuming surface air at level = 1. θ 5 and θ 2 represent the fifth and second levels of the near-surface thermal, respectively. The values of δ s on land and on ocean/ice surface were 1.0 K and 0.2 K, respectively. BRN was used to estimate the BLH [9,28]. The BRN at height z AGL was defined as where g denotes acceleration due to gravity, θ v0 and θ vz denote the virtual potential temperatures at the surface and height z, respectively, and u z and v z denote the wind speed components at height z. The critical threshold of 0.25 was used in the BLH estimation.

The Variance Method
The variance method defines turbulence as significant when a given threshold is exceeded. Previous studies show that BLH can be estimated based on a specified threshold by comparing with RS measurements. A threshold of 0.04 m 2 s −2 was initially set to estimate the BLH based on the TexAQS-2006 dataset [8]. BLHs from London, UK were defined at a height where σ 2 w > 0.1 m 2 s −2 [29], and this same threshold was used to estimate the BLHs in an urban environment in Beijing, China [20]. The SGP station was located near an urban environment. Therefore, a threshold of 0.1 m 2 s −2 was used in estimating BLH [12].

The Peak Method
The peak method was used to automatically detect BLH by using the DL σ 2 w . This method sets an automatic threshold based on the peaks for the estimation. A point in a discrete profile was defined as a peak when its value exceeds that of its neighbors. The peak method involves the following steps [24]: Finding thresholds based on peak values. First, given an original threshold δ 0 w , the peak method filters the tiny scale effects, and only the data above this threshold were used to find the peaks. This method was also used to calculate the background level δ b w , which refers to the mean velocity variance on the part of the profile below δ 0 w . This level can be computed as The peaks {p 1 , . . . , p m } above δ 0 w were then detected. The thresholds based on these peaks can be calculated as The BLH was equal to the highest point of threshold τ p and lower than the position of threshold δ 0 w . In this study, the δ 0 w was used to filter the effect of high cloudy layers; therefore, the height of the δ 0 w should be higher than the BLH. BLHs are referred to the height where the σ 2 w is equal to 0.1 m 2 s −2 . Therefore, δ 0 w = 0.05 m 2 s −2 was empirically set. In a profile where the values of σ 2 w are large or small, δ 0 w may be too small to filter the effect of high layer or too large to estimate the BLHs.

The WCT Method
The WCT method is widely used to estimate the BLHs based on data from aerosol lidar. In the WCT method, the Haar wavelet function was defined as [18] h elsewhere.
where r is the altitude, a is the spatial extent or dilation of the function, and b is the center of the Haar function. WCT of Haar function defined as where σ 2 w (r) is the profile about σ 2 w , r t and r b are the lower and the higher limit of the σ 2 w profile, respectively. When W f (a, b) reaches the local maximum with a coherent scale of a, value b is usually considered the PBLH.

The ED Method
ED is widely used in the field of image processing [30]. This method considers the gradient of a 2D change to accurately determine the abrupt edge of an image. For the time series images of DL, the BLHs were defined as the positions where a sudden change in wind speed is detected. Therefore, the image ED method provides another way of estimating BLH. The ED method involves three main steps. First, SNR and noise were used to control the quality of primary data shown in Figure 2a. Second, the edge of BLHs was detected by searching the edges of the turbulence region, which is located between the minimum and maximum peaks of σ 2 w . To get the minimum and maximum peaks, the number of peaks should be more than 2. The data in the turbulence region were set to 1 and others were set to 0; the binary image was shown in Figure 2b. Then, the difference between the adjacent data (back-front) was determined to obtain an edge with a value of either −1 or 1; to show a binary image of the edge, the difference result was calculated as the absolute value shown in Figure 2c. Third, the BLHs were estimated based on a theory of the diurnal cycle of PBL shown in Figure 2d. In general, the BLH is the position where the value of edge is equal to −1. However, the BL has the diurnal variation. At the daytime, the BLH was defined as the position of the weakened turbulence, to avoid the estimation of abnormal BLH, the information of the former BLH was used to limit the height of BLH. The BLH was the height with the least difference from the previous one. At nighttime, the PBL exists on two layers, the lower SBL, and the higher residual layer. Therefore, BLH was the position where the lowest value equals −1. Unfortunately, there are likely to be some abnormal BLH values in the middle, or cases where BLH cannot be calculated. When the BLH changes more than 300 m between a time resolution, it is regarded as abnormal. Therefore, the difference between BLH and the before and after values was computed. When the difference is greater than 300 m, the average value of the before and after BLH is used to replace the current BLH. In addition, this measure was used to fill the BLH that cannot be estimated.
The flowchart of ED is depicted in Figure 3. The aforementioned steps are described in detail, as follows: was the position where the lowest value equals −1. Unfortunately, there are likely to be some abnormal BLH values in the middle, or cases where BLH cannot be calculated. When the BLH changes more than 300 m between a time resolution, it is regarded as abnormal. Therefore, the difference between BLH and the before and after values was computed. When the difference is greater than 300 m, the average value of the before and after BLH is used to replace the current BLH. In addition, this measure was used to fill the BLH that cannot be estimated. The flowchart of ED is depicted in Figure 3. The aforementioned steps are described in detail, as follows: Data preprocessing. To reduce the influence of data noise, controlling the quality of data is necessary. The SNR threshold (SNR = 0.0072) [21] and noise threshold (noise = 1) [19] were set to filter the noise and obtain valid 2 w σ data.
Edge detection. The peak value of the data obtained in step 1 was obtained, and then the number of peaks (Np) was calculated. If Np exceeds 2, then the data exceeding the minimum peak value and below the maximum peak value were set to 1, whereas the others were set to 0 to generate a binary image. The difference between the adjacent data (back-front) was determined to obtain an edge with a value of either -1 or 1. Meanwhile, if Np is less than 2, then the data cannot be estimated by BLH.
BLH estimation. Given a time processing mechanism, at daytime (9:00-15:00), the height of the edge (value = -1) with the smallest difference from the BLH for the previous period was treated as the BLH. At nighttime, a residual layer was set above the boundary layer instead of connecting a free atmosphere between these layers. Therefore, the upper edge (value = -1) of the lower layer was treated as the BLH. To eliminate the influence of a single BLH value, the differences 1 D and 2 D were calculated as If 1 D and 2 D both exceeded 300 m, then the BLH was replaced by the average value of the former and latter BLHs. The BLH values were then outputted.

Results and Discussion
The 2 w σ from 2018 to 2020 were used to evaluate the performance of ED. The BLHs computed using the variance, peak, and ED methods were compared with the RS measurements in three ways. Including the two cases in the diurnal cycle analysis, the BLHs derived using the variance, peak, and ED methods were compared with the RS measure- Data preprocessing. To reduce the influence of data noise, controlling the quality of data is necessary. The SNR threshold (SNR = 0.0072) [21] and noise threshold (noise = 1) [19] were set to filter the noise and obtain valid σ 2 w data. Edge detection. The peak value of the data obtained in step 1 was obtained, and then the number of peaks (Np) was calculated. If Np exceeds 2, then the data exceeding the minimum peak value and below the maximum peak value were set to 1, whereas the others were set to 0 to generate a binary image. The difference between the adjacent data (back-front) was determined to obtain an edge with a value of either −1 or 1. Meanwhile, if Np is less than 2, then the data cannot be estimated by BLH.
BLH estimation. Given a time processing mechanism, at daytime (9:00-15:00), the height of the edge (value = −1) with the smallest difference from the BLH for the previous period was treated as the BLH. At nighttime, a residual layer was set above the boundary layer instead of connecting a free atmosphere between these layers. Therefore, the upper edge (value = −1) of the lower layer was treated as the BLH. To eliminate the influence of a single BLH value, the differences D 1 and D 2 were calculated as If D 1 and D 2 both exceeded 300 m, then the BLH was replaced by the average value of the former and latter BLHs. The BLH values were then outputted.

Results and Discussion
The σ 2 w from 2018 to 2020 were used to evaluate the performance of ED. The BLHs computed using the variance, peak, and ED methods were compared with the RS measurements in three ways. Including the two cases in the diurnal cycle analysis, the BLHs derived using the variance, peak, and ED methods were compared with the RS measurements under the NBL, SBL, and CBL, and a long-term analysis of the hourly and monthly BLHs for the diurnal and annual cycles was conducted. The BLHs retrieved via RS measurements from the SGP station of the ARM facility were treated as ground truths.

Analysis of Diurnal Cases
To evaluate the performance of ED, the BLHs derived using the variance, peak, and ED methods were compared with the RS measurements under extreme and routine atmospheric conditions. Two cases were considered to analyze the performance of ED.

Extreme Atmospheric Condition
The σ 2 w recorded on 4 March 2019 was used to estimate the BLHs via the variance, peak, and ED methods, and the estimated BLHs were then compared with the corresponding RS measurements shown in Figure 4. The red triangle represents the BLHs obtained using RS, and the color gradient represents the DL σ 2 w . The purple x, red cross, and black circle represent the BLHs derived from the variance, peak, and ED methods, respectively. Most of the BLHs obtained using the variance method were overestimated and had a chaotic distribution. Meanwhile, many of the BLHs obtained using the peak method were empty, and a few BLHs were overestimated. The BLHs obtained using ED were stable and were highly consistent with the RS measurements.
The variance and peak methods were easily affected by the value of σ 2 w . Therefore, these methods tend to underestimate or overestimate the BLHs due to the large variances in σ 2 w . By contrast, ED can steadily estimate the BLHs and is not affected by the variances in σ 2 w , and the BLHs recorded in the diurnal cycle were consistent with the RS measurements. The variance and peak methods performed poorly probably due to extreme atmospheric conditions, and using a single threshold cannot adapt well to atmospheric variations, thereby resulting in inaccurate estimations of BLHs. By contrast, ED estimates the BLHs based on the variation in data and can avoid the influence of atmospheric conditions. using RS, and the color gradient represents the DL 2 w σ . The purple x, red cross, and black circle represent the BLHs derived from the variance, peak, and ED methods, respectively. Most of the BLHs obtained using the variance method were overestimated and had a chaotic distribution. Meanwhile, many of the BLHs obtained using the peak method were empty, and a few BLHs were overestimated. The BLHs obtained using ED were stable and were highly consistent with the RS measurements. The variance and peak methods were easily affected by the value of 2 w σ . Therefore, these methods tend to underestimate or overestimate the BLHs due to the large variances in 2 w σ . By contrast, ED can steadily estimate the BLHs and is not affected by the variances in 2 w σ , and the BLHs recorded in the diurnal cycle were consistent with the RS measurements. The variance and peak methods performed poorly probably due to extreme atmospheric conditions, and using a single threshold cannot adapt well to atmospheric variations, thereby resulting in inaccurate estimations of BLHs. By contrast, ED estimates the BLHs based on the variation in data and can avoid the influence of atmospheric conditions.

Routine Atmospheric Condition
The 2 w σ recorded on 8 July 2019 was used to estimate the BLHs via the variance, peak, and ED methods and was compared with the corresponding RS measurements shown in Figure 5. The red triangle represents the BLHs obtained using RS, and the color gradient represents the 2 w σ . The purple x, red cross, and black circle represent the BLHs derived from the variance, peak, and ED methods, respectively. At nighttime (00:00-06:00 and 18:00-00:00), the BLHs obtained using the variance method were higher than the RS measurements. The peak method cannot calculate the BLHs. By contrast, ED can provide stable and continuous BLH measurements with slight differences from the RS measurements.

Routine Atmospheric Condition
The σ 2 w recorded on 8 July 2019 was used to estimate the BLHs via the variance, peak, and ED methods and was compared with the corresponding RS measurements shown in Figure 5. The red triangle represents the BLHs obtained using RS, and the color gradient represents the σ 2 w . The purple x, red cross, and black circle represent the BLHs derived from the variance, peak, and ED methods, respectively. At nighttime (00:00-06:00 and 18:00-00:00), the BLHs obtained using the variance method were higher than the RS measurements. The peak method cannot calculate the BLHs. By contrast, ED can provide stable and continuous BLH measurements with slight differences from the RS measurements.  The variance and peak methods performed poorly compared with ED at nighttime. The BLHs obtained using the variance and peak methods were consistent with the RS measurements from 06:00 to 12:00. But a few retrieved BLHs obtained using the variance and peak methods are higher and empty. In addition, the BLHs obtained by the variance and peak methods largely differed from the RS measurements retrieved between 12:00 and 18:00. By contrast, ED performs steadily and is highly consistent with RS measurements at daytime. The values of 2 w σ were much larger than the threshold, which may lead to the underestimation of BLHs on the nearest surface. The turbulence at daytime was stronger than that at nighttime. Therefore, the values of  The variance and peak methods performed poorly compared with ED at nighttime. The BLHs obtained using the variance and peak methods were consistent with the RS measurements from 06:00 to 12:00. But a few retrieved BLHs obtained using the variance and peak methods are higher and empty. In addition, the BLHs obtained by the variance and peak methods largely differed from the RS measurements retrieved between 12:00 and 18:00. By contrast, ED performs steadily and is highly consistent with RS measurements at daytime. The values of σ 2 w were much larger than the threshold, which may lead to the underestimation of BLHs on the nearest surface. The turbulence at daytime was stronger than that at nighttime. Therefore, the values of σ 2 w were relatively higher at daytime and were lower at nighttime. A single threshold of σ 2 w could not estimate the BLHs in the diurnal cycle. The peak method was unable to calculate the BLHs due to the lack of obviousness of the contour profiles. However, on the basis of the peak value of σ 2 w , ED extracted the atmospheric turbulence movement area to obtain the edge of BLHs and then selected the corresponding BLHs according to the diurnal change characteristics. Therefore, ED is not easily affected by the value of σ 2 w and returns better estimations of BLHs than the variance and peak methods.

Statistical Comparison
The BLHs obtained using the variance, peak, and ED methods were then compared with the RS measurements under NBL, SBL, and CBL. These BLHs were matched with the RS measurements for the same moment. The BLH measurements obtained using RS were selected as the middle values of the detection period. The valid samples of the peak method were selected to compute the R and RMSE between the peak and RS measurements.
Data from 2018 to 2020 were used to verify the consistency of the three methods with RS under NBL, SBL, and CBL, respectively, as shown in Figure 6. A redder color corresponds to denser dots and more BLHs. The black dotted line represents the reference line with a slope of 1. On the NBL, SBL, and CBL, the R obtained by ED was larger than those obtained by the variance and peak methods (R: 0. Under the SBL, the MRE of the ED method is lower than the variance method and higher than the peak method (−38.64%, vs. −152.23%; 14.02%), respectively. Under the CBL, the MRE of the ED method is lower than the variance method and higher than the peak method (−15.07% vs. 2.24%; 5.64%), respectively. In addition, more BLHs (points at which color change from deep red to shallow blue) derived from the ED method are close to the reference line than those derived from the variance and peak methods. The lower MRE of the ED method makes its result reliable. In conclusion, the result of the ED method is more consistent with the variance and peak methods under the same conditions. Although, the R of these methods were all low.
Huge differences were observed between the BLHs obtained by the variance and peak methods and the RS measurements, and most of the BLHs obtained by the traditional methods were either overestimated or underestimated. Few of the BLHs obtained using ED were either overestimated or underestimated, but the majority approached the reference line. The temporal resolution of the BLHs based on DL was 10 min, and the detection time for the RS measurements was about 30 min from launch to detection. After its launching, the balloon may change its location due to atmospheric winds. Therefore, the time and position of the retrieved BLHs and RS measurements might not match, which may explain the weak R and large RMSE. In addition, the σ 2 w in clear air was retrieved below the lowest cloud base, which may explain why the above methods underestimated the BLHs occasionally. Any method has its drawback. The ED method cannot estimate the BLHs when peaks of a profile are below two, which is determined by the situation of the data. Though we will use the mean substitution for adjacent data, sometimes, some errors can occur. The ED method considers the time dimension of the BLHs; it can fill the missing of uncalculated BLHs and eliminate a single abnormal BLH. The continued abnormal BLHs cannot be identified, and the missing long-term BLHs will be replaced by one value. Then, the BLHs may lose reliability. More importantly, we found that the SNR data filtered more information of σ 2 w (such as the white area in Figures 4 and 5). Therefore, improving the SNR of the data may improve the results. −152.23%; 14.02%), respectively. Under the CBL, the MRE of the ED method is lower than the variance method and higher than the peak method (−15.07% vs. 2.24%; 5.64%), respectively. In addition, more BLHs (points at which color change from deep red to shallow blue) derived from the ED method are close to the reference line than those derived from the variance and peak methods. The lower MRE of the ED method makes its result reliable. In conclusion, the result of the ED method is more consistent with the variance and peak methods under the same conditions. Although, the R of these methods were all low. Figure 6. BLHs from 2018 to 2020 obtained using the variance, peak, WCT, and ED methods compared with the RS measurements.
Huge differences were observed between the BLHs obtained by the variance and peak methods and the RS measurements, and most of the BLHs obtained by the traditional methods were either overestimated or underestimated. Few of the BLHs obtained using ED were either overestimated or underestimated, but the majority approached the reference line. The temporal resolution of the BLHs based on DL was 10 min, and the detection Figure 6. BLHs from 2018 to 2020 obtained using the variance, peak, WCT, and ED methods compared with the RS measurements.
We used the WCT methods to estimate the BLHs for the DL dataset from 2018 to 2020 and compared those results with the RS measurements under the NBL, SBL, and CBL conditions, as shown in Figure 6g in clear air that are captured below the lowest cloud base, the clouds' effect was eliminated. The results of the ED method and the WCT method are similar, and the performance of the ED method is even better, combined with the analysis of the above factors. In addition, many BLHs derived from the WCT method were significantly underestimated under the NBL condition and overestimated under the SBL condition.
Compared with those methods proposed to estimate the BLHs based on the aerosol lidar, the ED method has a poor performance. Dang, Yang [31] improved the WCT method with a given top limiter altitude to reduce the interference of RL and cloud layer on BLHs estimation. The method was well performed and showed good results (R = 0.79) validated by radiosonde. Li, Chang [32] proposed the combining of the whale optimization algorithm and the top limit method to estimate the BLHs based on micro-pulse lidar. The R between this method, WCT, and gradient methods and RS measurements are 0.88, 0.82, and 0.78, respectively. Sleeman, Halem [33] proposed a deep learning method by employing deep neural networks for denoising and image segmentation to estimate the BLH, and it can infer BLHs in the presence of dense clouds. The results for this method are similar to the WCT method compared with the RS measurements. Among these methods, the WCT method was selected as the reference to estimated BLH from DL. Unfortunately, the results retrieved from the WCT method are not very well in this study, possibly because of the lower SNR and vertical resolution of the DL data. It is necessary to improve the quality of the DL signal for accurate BLHs. Comparatively speaking, the ED method also has better performance with the WCT method under the NBL, SBL, and CBL conditions. The poor performance of the ED method motivates us to improve or/and develop more practical approaches.

Long-Term Analysis
In this section, the BLHs from 2018 to 2020 derived using the variance, peak, and ED methods were used to analyze the long-term variations in the BLHs. Therefore, hourly average BLHs were used for the diurnal analysis, and monthly average BLHs were used for the annual analysis. The DL data for November and December 2020 were limited.

Diurnal Cycle
The hourly average BLHs are shown in Figure 7. The blue, orange, and yellow point lines represent the mean BLHs per hour in 2019 retrieved using the variance, peak, and ED methods, respectively. The purple asterisk represents the RS measurements for each hourly BLH. The hourly average BLHs are shown in Figure 7. The blue, orange, and yellow point lines represent the mean BLHs per hour in 2019 retrieved using the variance, peak, and ED methods, respectively. The purple asterisk represents the RS measurements for each hourly BLH. The hourly BLHs retrieved using the variance method were higher and lower than those obtained by RS measurements at nighttime (from 00:00 to 06:00 and from 18:00 to 24:00). BLHs obtained using the variance method were higher than RS measurements at 12:00. The BLHs remained stable from 07:00 to 14:00 and increased from 14:00 to 17:00. The BLHs obtained using the peak method slowly increased from 06:00 to 12:00 and decreased from 12:00 to 18:00. The BLHs retrieved using ED increased from 06:00 to 16:00, decreased from 16:00 to 18:00, and remained stable from 18:00 to 24:00. The difference between the peak and ED methods and RS measurements are both small.
The BLHs generally increased at sunrise and decreased at sunset. Therefore, the daytime BLHs should be higher than the nighttime BLHs. In the long-term analysis for the The hourly BLHs retrieved using the variance method were higher and lower than those obtained by RS measurements at nighttime (from 00:00 to 06:00 and from 18:00 to 24:00). BLHs obtained using the variance method were higher than RS measurements at 12:00. The BLHs remained stable from 07:00 to 14:00 and increased from 14:00 to 17:00. The BLHs obtained using the peak method slowly increased from 06:00 to 12:00 and decreased from 12:00 to 18:00. The BLHs retrieved using ED increased from 06:00 to 16:00, decreased from 16:00 to 18:00, and remained stable from 18:00 to 24:00. The difference between the peak and ED methods and RS measurements are both small.
The BLHs generally increased at sunrise and decreased at sunset. Therefore, the daytime BLHs should be higher than the nighttime BLHs. In the long-term analysis for the diurnal cycle, the BLHs retrieved by ED were more consistent with theory compared with those retrieved using the variance and peak methods. Moreover, the standard error of the ED method was smaller than those of the variance and peak methods, thereby suggesting that the BLHs retrieved by ED were close to the mean values.

Annual Cycle
The monthly average BLHs derived using the variance, peak, and ED methods and the RS measurements are shown in Figure 8. The blue, orange, yellow, and purple point lines represent the monthly mean BLHs from 2018 to 2020, retrieved using the variance, peak, and ED methods and RS measurements in the annual cycle, respectively. The monthly BLHs for the variances greatly differed from the RS measurements, particularly from June 2018 to June 2019. The monthly BLHs retrieved using the ED and peak methods were consistent with the RS measurements. However, there are differences between April and May 2018, and 2019. The monthly BLHs for November and December 2020 were limited. The BLHs from 2018 to 2020 retrieved using ED showed consistent trends with the RS measurements, although there is a big difference for these monthly BLHs in October 2020. The annual cycle BLHs showed that the BLHs in summer and autumn were higher than those in spring and winter. The BLHs increased slowly from January to May and stably decreased from June to December.

Conclusions
This paper introduced ED for estimating BLHs based on DL 2 w σ . The two cases on 4 March 2019 and 8 July 2019 were used to estimate the BLHs via the variance, peak, and ED methods in the diurnal cycle. The retrieved BLHs were compared with the RS measurements, and the results show that ED was more accurate than the variance and peak methods under extreme and routine atmospheric conditions. Therefore, ED can estimate the BLHs in the diurnal cycle and perform better than the other two methods. The BLHs retrieved using the variance, peak, WCT, and ED methods from 2018 to 2020 were compared with the RS measurements under NBL, SBL, and CBL. The results show that ED was more consistent and had lower RMSE with the RS measurements than the variance and peak methods under NBL, SBL, and CBL. The lower MAE and MRE indicated the results of the ED method with reliability, and a number of the BLHs derived from the ED method are closer to the reference line than those obtained by the variance and peak methods. In addition, the ED method has similar results with the WCT method and is even better. In the long-term analysis, the hourly BLHs retrieved by ED were consistent with theory, whereas those retrieved using the variance and peak methods were not. The monthly BLHs retrieved using the peak and ED methods were more consistent with the

Conclusions
This paper introduced ED for estimating BLHs based on DL σ 2 w . The two cases on 4 March 2019 and 8 July 2019 were used to estimate the BLHs via the variance, peak, and ED methods in the diurnal cycle. The retrieved BLHs were compared with the RS measurements, and the results show that ED was more accurate than the variance and peak methods under extreme and routine atmospheric conditions. Therefore, ED can estimate the BLHs in the diurnal cycle and perform better than the other two methods. The BLHs retrieved using the variance, peak, WCT, and ED methods from 2018 to 2020 were compared with the RS measurements under NBL, SBL, and CBL. The results show that ED was more consistent and had lower RMSE with the RS measurements than the variance and peak methods under NBL, SBL, and CBL. The lower MAE and MRE indicated the results of the ED method with reliability, and a number of the BLHs derived from the ED method are closer to the reference line than those obtained by the variance and peak methods. In addition, the ED method has similar results with the WCT method and is even better. In the long-term analysis, the hourly BLHs retrieved by ED were consistent with theory, whereas those retrieved using the variance and peak methods were not. The monthly BLHs retrieved using the peak and ED methods were more consistent with the RS measurements than those obtained using the variance method. The BLHs recorded in summer and autumn were also higher than those in spring and winter.
In general, the ED method can estimate the BLHs in the diurnal cycle under extreme and routine atmospheric conditions and has a better performance than the variance, peak, and WCT methods. But the performance of these methods are unsatisfactory at present. Direct detection of turbulence changes can obtain accurate BLH theoretically, but the low SNR of DL may be one of the reasons that reduces the accuracy of the BLHs obtained from this remote sensing equipment. The poor performance of the ED method motivates us to improve or/and develop more practical approaches moving forward.