Atmospheric Boundary Layer height estimation from aerosol lidar: a new approach based on morphological image processing techniques

The Atmospheric Boundary Layer (ABL) represents the lowermost part of the atmosphere directly in contact with the Earth surface. The estimation of its depth is of crucial importance in meteorology and for anthropogenic pollution studies. The ABL Height (ABLH) measurements are usually far from being adequate, both spatially and temporally. Thus, different remote sensing sources can be of great help in growing both the spatial and temporal ABLH measurement capabilities. To this aim, aerosol backscatter profiles are widely used as proxy to retrieve the ABLH. Hence, the scientific community is making 5 remarkable efforts in developing automatic ABLH retrieval algorithms applied to lidar observations. In this paper, we propose a ABLH estimation algorithm based on image processing techniques applied to the composite image of the total attenuated backscatter coefficient. A pre-processing step is applied to the composite total backscatter image based on morphological filters to properly set-up and adjust the image to detect edges. As final step, the detected edges are post-processed through both mathematical morphology and an object-based analysis. The performance of the proposed approach is assessed on real 10 data acquired by two different lidar systems, deployed in Potenza (Italy) and Evora (Portugal), belonging to the EARLINET network. The proposed approach has shown higher performance than the benchmark consisting of some state-of-the-art ABLH estimation methods.

cross-polarization channel measurement is available for the lidar instrument, the ABLH can be retrieved from the changes in the lidar depolarization ratio as showed in (Bravo-Aranda et al., 2017). (Hooper and Eloranta, 1986) retrieved the ABLH identifying the maximum variance of the backscattered signal, while in other several studies the ABLH is retrieved applying the wavelet covariance transform (WCT) to the lidar signal (Cohn and Angevine, 2000;Davis et al., 2000;Talianu et al., 2006; 60 Compton et al., 2013). (Comerón et al., 2013) put in evidence that a strong link exists between the wavelet transform and the gradient method. Other hybrid methods combine the wavelet approach with the image processing techniques (Morille et al., 2007;Lewis et al., 2013) or pair the measured aerosol-backscatter lidar profile with models, e.g. multi-wavelength numerical model (Dionisi et al., 2018); stability dependent model of ABLH temporal variation (Su et al., 2020). Recently, sophisticated algorithms that include very advanced techniques were proposed: (De Bruine et al., 2017) developed an algorithm based on 65 tracking (pathfinder) that employs graph theory to track the diurnal evolution of the ABLH. (Toledo et al., 2014) applied for the first time the Cluster Analysis (CA) to lidar measurements to obtain the ABLH. Another recent technology takes advantage of an Ensemble Kalman Filter (EKF) to trace ABLH evolution from ground-based lidar observations as shown in (Tomás et al., 2010). A more detailed review with more information about advantages and drawbacks of the reported algorithms can be found in (Dang et al., 2019). 70 In this manuscript we present the development and validation of a new algorithm to continuously retrieve the ABLH from measurements taken at different permanent observational sites deployed in the frame of the European Aerosol Lidar Network (EARLINET) included in the ACTRIS research infrastructure (Pappalardo et al., 2014). The new implemented algorithm is using a fully image-based methodology (2-D): instead of analyzing the lidar observations profile by profile, the retrieval takes into account also the temporal dimension (at very high resolution). The algorithm framework needs as input a statistically 75 significant temporal collection of the total attenuated backscatter vertically-resolved profiles. The retrieval consists in applying to the composite image during the pre-processing phase the morphological operators and the edge detectors. Instead, during the post-processing phase, after applying the morphological filters again, the significant edges are extracted through an objectbased analysis that is proven to be particularly successful in determining the ABLH. This proof-of-concept is a useful starting point to develop a central common strategy to produce high-quality and reliable ABLH retrieval in the frame of the Global 80 Atmosphere Watch (GAW) Aerosol Lidar Observation Network (GALION; (Bösenberg et al., 2008)) project of the WMO, which has as main objective to harmonize all the existing lidar and ceilometer networks. Indeed, the proposed approach is able, after a proper parameter tuning phase, to work with a huge variety of data addressing the task of the ABLH estimation even for data acquired by simpler networks or less advanced systems.
The manuscript is organized as follows. Sect. 2 is devoted to the presentation of widely used state-of-the-art methods to 85 estimate the ABLH. In Sect. 3, the proposed image processing based approach is described. Afterwards, the experimental results are shown in Sect. 4. Finally, conclusions are drawn in Sect. 5.
There are several methods to determine the ABLH from lidar observations that are based on the assumption that aerosol is trapped within the ABL. Those methods find the height where the aerosol concentration abruptly decreases. This happens because the aerosol particles within the ABL can be used as proxy to study the boundary-layer vertical structure and time variability. In fact, aerosols uplifted after sunrise by convective mixing can act as efficient tracers for the atmospheric portion over which mixing occurs. Aerosols can also be dispersed out of the ABL during strong convective events or temporary breaks 100 of the capping temperature inversion. Thus, elastic backscattered signals from aerosol particles measured by lidar systems can be used to determine the height and the internal structure of the ABL and, when possible, the residual layer and aerosol layers within and aloft the ABL (Di Girolamo et al., 1999).
For lidar systems, typically the detected backscattered light is much higher within the ABL than in the free troposphere due to the higher abundance of particles. The lidar equation is defined as: where λ 0 and λ are the emitted and received lidar wavelength (laser wavelength), z is the vertical height, O(z) is the overlap function, A refers to a system function (area and typical configuration), β mol and β par are the backscatter coefficients for molecular and particle components, respectively, T mol and T par indicate the atmospheric transmissivity and P bgd is the background signal. Often it is preferred to use the corrected signal for the square of the quota, named range corrected signal (RCS), 110 defined as: In order to determine an estimate of the height of ABLH, we directly apply the derivative method exploiting the logarithm of the quantity (2). In this case, the lidar elastic backscatter signal, P (λ, z), is used. Thus, we have that the ABLH, H(λ), can be obtained as follows: 115 The minimum of the quantity in (3) identifies the transitions between different layers, and the absolute minimum identifies the height of the boundary layer, because the largest variation of the lidar signal is considered corresponding to the largest 4 https://doi.org/10.5194/acp-2020-857 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.
variation of the aerosol concentration. H(λ) represents the transition from the stable layer to a neutral or unstable condition above. The method has been applied to the maximum vertical spatial and temporal resolutions. 120

Wavelet Covariance Transform for ABLH Estimation
The Wavelet Covariance Transform (WCT) is defined as (Brooks, 2003) where s(z) is the range-corrected lidar backscatter signal, z is the measurement height, z b and z t are the lower and upper limits of the lidar return signal profile, respectively, a is the dilation factor, b the translation, and h is defined as the Haar wavelet 125 function, i.e.
In order to have a more efficient implementation and to gain more insights about WCT, (4) can be seen as convolution. Thus, starting from (4) we have: * stands for the convolution operator and z b and z t are set to −∞ and +∞, respectively, without harming the generality. It is worth to be pointed out that the Haar function h(x) in (5) can be rewritten using the derivative and a triangular window. Hence, we have that: where is the triangular window.
analysis for linear time-invariant (LTI) systems. The convolution in time can be seen in the transformed domain as (Bracewell and Bracewell, 1986): where is the forward Fourier transform. Considering the derivative and the scaling properties of the Fourier transform (Bracewell and Bracewell, 1986) and that the dilation a ≥ 0, after simple 145 algebra, starting from (7) and (8) we have that: where i is the imaginary unit and is the normalized sine cardinal function.

150
Now if we consider the modulus of C(f ), i.e., |C(f )|, in (10) Thus, it is easy to see that if f → ±∞, |W (f )| → 0 and if f = 0, |W (f )| = 0. Indeed, the Fourier transform of the triangular window in (8) can be seen as a low-pass filter with a cut-off depending on a (i.e., the factor in (13) represented by the 155 sinc 2 function). Instead, the derivative in (8) leads to a multiplication by |f | in (13) generating a band-pass filter with cut-off frequencies depending on a, which rules the selection of a portion of the whole frequency spectrum. Hence, we have from (11) that higher frequencies of s(b) are passed with respect to the low-pass filter defined by the triangular windows. In this sense, the WCT method with the Haar wavelet function can be considered as a particular gradient-based method.
In all the equations above, we neglected the dependence on a of the WCT c. Indeed, the dilation a is set to a fixed value.

160
A rule of thumb can be found in (Baars et al., 2008). However, in this paper, in order to have a high performance method for our benchmark, we set a to its optimal value (which is sensor-dependent) obtained via a grid search approach. Furthermore, following the indications in (Baars et al., 2008), we normalize the range-corrected signal by its maximum value found below a given height (usually set around 1000 m). The normalization guarantees the applicability of a unique threshold (set to 0.05 as suggested in (Baars et al., 2008)) on c(b) in order to find the ABL even at very different backscatter conditions in rather clean 165 or very polluted air.

The Proposed Morphological Image Processing Approach (MIPA)
The composite image on which the proposed algorithm is applied consists in the temporal sequence of the range-corrected vertically-resolved acquired lidar profiles. The proposed MIPA algorithm has no prior knowledge on ABLH and the image processing approach relies on: i) a block that reduces the vertical spatial resolution to reach a working spatial resolution around 20m if the spatial resolution of the system is finer; ii) a pre-processing step applied to the daily lidar data exploiting mathematical morphology; iii) Canny's edge detection (Canny, 1986) applied to the pre-processed data; iv) a post-processing starting from the detected edges and based on both mathematical morphology and an object-based analysis to get the final outcome. In the following sections, the basics of morphological operators are presented first. Then, the proposed MIPA algorithm will be detailed.

Basics of Morphological Operators
An image I: E ⊆ Z 2 → V ⊆ Z is analyzed by morphological operators via the so-called structuring element (SE), here denoted as B (Soille, 2003), which can be defined through its spatial support N B (x) (i.e., the neighborhood with respect to the position x ∈ E in which B is centered) and by its values. For flat SEs (i.e., SEs with unitary values), the only free parameters are the origin and N B . We will focus on these SEs in this work.

180
Erosion, ε B [I], and dilation, δ B [I], are the two basic operators defined, for each x ∈ I, as follows: where S and S are the infimum and supremum values in the set S, respectively.

185
The erosion (respectively dilation) application has as filtering effect that suppresses bright (respectively dark) regions smaller than B and the enlargement of dark (respectively bright) ones. For bright and dark regions we mean that the local contrast in a certain region has intensity values greater or lower with respect to the surrounding ones, respectively. Erosion and dilation operators can be recast into minimum and maximum operators on B, respectively, if I is a binary image.
We also introduce for convenience the opening and closing that correspond to the two possible sequential compositions of 190 erosion and dilation. In particular, the opening is defined as: whereB denoting the SE obtained by reflecting B with respect to its origin. Instead, the closing is given by A closing removes dark regions smaller than B, whereas an opening suppresses bright ones. For further details about mor-195 phological operators, the interesting readers can refer to the related literature (Soille, 2003).
A number of morphological operators can be obtained by properly combining the above-mentioned operators. Two instances, which are of great interest in this work, are represented by the residuals of the application of erosion and dilation, usually called internal gradient and external gradient, respectively (Soille, 2003). In particular, the internal gradient, ρ − B [I], is defined as 200 and the external gradient, ρ + B [I], is given by These two gradients are also often called Half-Gradients (HGs).
In the remaining of this section, the four modules that constitutes the proposed approach will be detailed.

205
This first block starts from a matrix I: E ⊆ Z 2 → V ⊆ Z that is the daily sequence of the attenuated backscatter profiles. These latter form the columns of I. Thus, the number of rows is related to the maximum height and the spatial resolution of the system, instead the number of columns is about its temporal resolution. The downsampling with a factor R, which aims to reduce the bins' spatial resolution, is implemented by a low-pass filter along each column of I plus decimation with a factor R.
In particular, a moving-average filter is simply applied as low-pass filtering. The support (i.e., the length of the sliding window) 210 of the filter is R, again. Thus, R is the unique tuning parameter that is selected in order to have a spatial resolution not finer than 20 m. Hence, R can be directly calculated from the system's spatial resolution bringing the initial product to the target spatial resolution. This operation is performed to have a sharper edge defining the ABL. The outcome is an image I D used as input in the pre-processing step. It is worth to be remarked that if we work with data having a spatial resolution coarser than 20 m, R is set to one, thus skipping this step implying that I D is equal to I.

Pre-processing Based on Mathematical Morphology
In this work, we propose the use of a low-pass filter based on HGs to pre-process I D . Before coming into the details of the used operator, let us define the complement of a generic operator Ψ, i.e., Ψ, as Ψ = id − Ψ, where id is the identity operator. It is worth to be remarked that when discontinuities are present both the HGs assume positive values constituting an approximation of the norm of the signal gradient (Soille, 2003). Thus, the difference of the two (internal and external) HGs represents a detail 220 extraction operator, since it reproduces the variations of the function with respect to the local mean (Soille, 2003). This operator exploiting a SE B pre , denoted as Ψ HG,Bpre , is defined as follows: 8 https://doi.org/10.5194/acp-2020-857 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.
in which the factor 0.5 is applied to preserve the property of approximating the image gradient norm. The corresponding 225 low-pass filter is simply given by the complementary operator of Ψ HG,Bpre , i.e.
which corresponds to the semi-sum of the dilation and erosion operators. This operator is used in the pre-processing phase 230 applying it to I D and fixing B pre to a line SE in the horizontal direction (i.e., along the time direction) with a length l pre . It enables us to smooth the lidar image along the horizontal axis (where the dynamic of the ABL is quite slow), directionally reducing the noise and preserving the vertical edges that will be of crucial importance for the next step. The resulting image after pre-processing the image I D is indicated as I pre and it represents the input of the edge detection block.

235
This processing step can be implemented in several ways. Every edge detector can be exploited to extract a first estimation of the ABL starting from I pre . Thus, the proposed approach is flexible and this block can be changed to possibly improve the results.
Approaches, already discussed in this paper, as Wavelet Covariance Transform (WCT) or gradient-based could be adopted.
However, the analysis of the performance varying these strategies in the proposed framework is out-of-the-scope of this paper.
Thus, we employed Canny's edge detector (the classical version available in commercial software as MATLAB) (Canny, 1986) 240 to get the first estimation of the edge map, denoted as E. The detected edges in E are indicated with 1, instead, the rest of the map (background) is labeled as 0. All the bins labeled as 1 in the edge map are potential candidates to represent the ABL.

Post-processing
After applying Canny's edge detector to the pre-processed data, the edge map E is analyzed by using further signal processing.
In particular, morphological filters are exploited first. Then, a final post-processing phase relied upon an object-based analysis 245 is performed. The two steps will be detailed in the following.

Post-processing Based on Morphological Operators
The post-processing based on morphological filters is applied to the edge map E. This processing step is about removing unrealistic edges (i.e., edges that vary too fast with respect to the dynamic of the ABL). Thus, we apply a series of directional Low-Pass (LP) morphological filters. In particular, the used filters are obtained by sequentially applying an opening and a 250 closing operator using the same structuring element B post , i.e.
where B post is a line SE with a length l post and an angle θ. The application to E of these directional filters varying θ from θ min to θ max and combining the outputs with a maximum operator provides the output of this post-processing step, indicated as E post . 255 3.5.2 Object-based Post-processing The object-based post-processing is applied to the edge map E post . The detected edges are indicated with 1, instead, the rest of the map (background) is labeled as 0. The first layer in E post consists of the first (starting from the ground) bins detected as "edge" (i.e., labeled as 1) analyzing the edge map for each profile (i.e., in the vertical direction).
We work on these edges extracting objects. The main concept is the use of the connectivity in an edge map, i.e., the way 260 in which the bins labeled as "edge" (which assume value 1 in the edge map) are spatially related to their neighbors. A bin declared as "edge" is said "8-connected" if exists at least a bin belonging to its 8-neighborhood, i.e., the adjacent bins in vertical, horizontal, and diagonal directions, declared as "edge", as well. All the bins that are "8-connected" to each other form an object.
Thus, several objects, clustering the bins declared as "edge", are collected and analyzed. In particular, an analysis about the 265 spatial variability of these objects is performed. Indeed, if the absolute Euclidean distance between the mean of the heights for each extracted object (using the connectivity procedure explained above) and the related mean calculated on the objects in its neighborhood exceeds a threshold δ post , this object is removed from the solution. Finally, the outcome, i.e. the estimated ABL denoted as E out , is obtained by linearly interpolating the remaining objects in the edge map.

270
Algorithm 1 (MIPA) summarizes the sequence of the adopted signal processing steps in order to provide to the readers a complete overview of the proposed approach.
Algorithm 1 The proposed ABL estimation algorithm.
-Reduce the profile resolution of I by a factor R to get ID -Pre-process ID by low-pass filtering using HGs, see (24), as described in Sect. 3.3, to get Ipre -Estimate the edges of Ipre using Canny's edge detector obtaining the edge map E -Post-process the edge map E, as described in Sect. 3.5, using directional morphological filters as in (25)  Even if MUSA and PAOLI operate at the same wavelengths, their technical characteristics are very different: laser source, telescope, detection and acquisition system, space and time resolution, full overlap region. As a consequence, applying MIPA algorithm on both systems is a good benchmark to evaluate the algorithm performances.
MIPA algorithm uses as input the composite plot of the vertically-resolved attenuated backscattering coefficient at 1064 nm.

295
High spatial and temporal resolution is needed to increase sensibility in using the directional morphological filter while long

305
The total attenuated backscatter can be easily expressed in terms of measured lidar signals taking into account (1) and (2): where K is a calibration constant determined by imposing that β att (λ, z) = β mol (λ, z)T 2 mol (λ, z) in an aerosol-free atmospheric region. It is important to note that morphological filter techniques rely only on the correlation among adjacent lidar range bin and not on their absolute numeric values. As a consequence, the ABLHs computed by the proposed MIPA algorithm are rather insensitive to the accuracy of the calibration constant K. Actually, the morphological algorithm can be applied to the range corrected signal (P rcs ) time series instead of the total attenuated backscatter one, providing the same results in terms of ABLH. This is in general not true for the traditional ABLH retrieval algorithms (derivative, WCT) where proper thresholds on absolute signal value need to be defined. Tab. 1 summarizes all the input parameters exploited by the compared approaches for the two lidar systems considered. These have been defined after a tuning phase on different lidar scenarios. It is worth to be 315 remarked that the derivative approach requires no parameter to be set working at the higher resolutions possible.
All the input data sets considered in the study have been previously pre-processed at high resolution by using the EARLINET SCC. In particular, starting from raw lidar data, several instrumental corrections (like for example trigger delay correction, dead time correction, analog and photoncounting signal glueing, etc.) and all the required raw signal handlings (like atmospheric background subtraction, range correction) have been applied. More details on the pre-processing procedure implemented in 320 the EARLINET SCC are described in (D'Amico et al., 2016).
To asses the performance of all the considered ABLH retrievals, a reference needs to be set. As the most rigorous definition of the PBL is the one based on thermodynamic effects, we assume as reference for the ABLH the values retrieved from atmospheric temperature and pressure profiles of the co-located radiosondes and the ECMWF model (if available) related to positions closest to the site of the lidar system under study. The lidar retrieval uses the aerosols as tracers identifying the height 325 of the maximum vertical gradient indicating the drop in aerosol concentration. This corresponds to that part of the atmosphere where convection and turbulence drops and we get the ABLH to be used as reference for the approaches estimating the ABL from lidar data (Summa et al., 2013).
The sounding data available for all the measurement dates considered in this study are quite few and except from one case they are not co-located with the lidar station. In particular, there are no sounding data available for Potenza for any of We have used the atmospheric temperature and pressure profiles provided by the Numerical Weather Prediction (NWP) model over the two measurement sites for assessing the performance. Even though sounding data can be available sometimes (as reported in this work when available for the scenarios under test), these data are not enough to guarantee the calcula-tion of a reliable reference for the considered ABL retrievals due to the very low temporal resolution of these acquisitions.
To our knowledge, NWP is the best alternative to the use of sounding data. In particular, to calculate the ABLH reference points for all the Potenza cases, we have used the high-resolution NWP provided by the ECMWF Integrated Forecast System  Often, the aerosols in the free troposphere tend to intrude in the ABL making the separation between the ABL and the upper 365 atmospheric region not so clear. This is clearly visible around noon on 10 July and 11 July 2012. The general idea behind the retrieval of ABLH from lidar measurements is to use the aerosol as tracers of the ABL and assuming that the ABLH is given by the top of the first aerosol layer (starting from the surface). According to this assumption, we can clearly see that the ABLH is maximum around noon (when the solar convention is at its maximum) for all the three days and it reaches its minimum during nighttime. The black line in Fig. 1 shows the ABLH as retrieved by MIPA algorithm. In general, the expected temporal 370 evolution of the ABLH is well captured and the intrusions of the upper aerosol layers in the ABL do not seem to affect the outcomes, thus still obtaining reasonable ABLH estimates.

EARLINET 72h Exercise
It is important to highlight that the assumption to retrieve the ABLH as the top of the first detected aerosol layer is valid only if the ABL is above the full lidar overlap height. During night-time observations, this condition may be not always verified and in this case the algorithm detects as first edge the base of the layer next to the ABL top. This is clearly visible during the 375 nighttime period in Fig. 1 where the real ABL is too low to be detected with the MUSA system and the retrieved ABLH is typically overestimated (Di Girolamo et al., 1999).
As explained in Sect. 3, the MIPA algorithm is composed by four different modules: an edge detector based on Canny filter (see Sect. 3.4), a module for reducing the vertical resolution (see Sect. 3.2), the pre-processing module described in Sect. 3.3 Fig. 2 shows the role played by each of these modules in retrieving the ABLH for the Potenza 72h dataset and also the reference ABLH retrieved using the ECMWF forecasts (brown circles). The ABLHs shown in Fig. 2 were calculated considering the different steps of the MIPA algorithm, i.e. applying: the edge detector module on full resolution data (curve a in Fig. 2), the vertical resolution reduction module and the edge detector (curve b), the vertical resolution reduction module, the preprocessing module and the edge detector (curve c) and finally the whole algorithm (curve d). It is evident that the edge detector 385 is not sufficient to ensure a stable retrieval even if it is applied to a reduced resolution dataset. Consequently, the application of the post-processing module is crucial in the processing. The absolute differences of all the curves shown in Fig. 2 with respect to the reference are plotted in Fig. 3. In generating Fig. 3 (and all other figures showing the absolute difference with respect to the reference), the reference data have been interpolated at the same resolution of the lidar data assuming that the ABL is slowly varying. 390 Fig. 4 shows the comparison between the ABLH retrieved by the MIPA algorithm and the corresponding one obtained by using more traditional ABLH detection approaches. Specifically, we applied to the same dataset different algorithms for the ABLH estimation: the MIPA algorithm 1, the derivative method described in Sect. 2.1 and the procedure based on WCT as in Sect. 2.2. The ABLH retrieved using the ECMWF forecasts is reported as brown circles in Fig. 4. The agreement among the ABLH obtained by applying the considered algorithms on lidar data is satisfactory. The MIPA algorithm provides the smallest 395 discrepancies with respect to the reference. As expected, all the algorithms overestimate the ABLH in nighttime conditions where the real ABLH as retrieved by ECMWF forecasts data is clearly below the full overlap of the MUSA lidar, which is about 300 m.
The absolute difference of all the considered algorithms with respect to the reference is plotted in Fig. 5 and a statistical analysis is reported in Tab. 2 for all the considered algorithms. Both Fig. 5 and Tab. 2 confirm the better performance of the 400 MIPA algorithm with respect to the other approaches into the proposed benchmark. EARLINET station) are also reported as green triangles. The very good agreement with the corresponding ABLH retrieved by using the GDAS forecasts proves the good accuracy of the selected reference.
Tab. 3 sums up the statistical analysis of the absolute differences with respect to the reference for all the considered algorithms for the 72h Evora dataset. As for Potenza, the better performance is the one corresponding to the proposed ABLH casts and the outcomes of the MIPA approach, instead, during nighttime, the MIPA retrievals may be hindered by PAOLI's overlap for low ABLHs.

Other Potenza Cases
In this section we describe three additional case studies on which the proposed MIPA algorithm has been evaluated. The three cases refer all to MUSA lidar observations during three longer continuous observations ever made over Potenza with MUSA.

420
The first case study refers to April 21, 2010 when the Potenza EARLINET station performed a quite long record of lidar measurement concurrently with the eruption of the Icelandic volcano Eyjafjallajökull occurred in April-May 2010. During the eruption almost all the EARLINET stations promptly started continuous measurements, whenever weather conditions allowed it and the corresponding outcomes are described by (Pappalardo et al., 2013). Fig. 9 shows the comparison between the ABLH retrieved by the MIPA algorithm and the corresponding ABLH obtained 425 by using the other considered algorithms. ABLHs retrieved from ECMWF runs are assumed as reference (brown circles). The absolute difference of all the considered algorithms with respect to the reference is plotted in Fig. 10. For this case, from 08:30UTC to 17:00UTC, all the considered lidar retrieval algorithms give ABLH values systematically below the reference ones. This behavior can be explained looking at Fig. 11, where the total attenuated backscatter time series at 1064 nm measured by the MUSA system on April 21, 2010 is shown. Moreover, the black line shows the ABLH retrieved by MIPA and the white 430 line shows the reference ABLH. Before the 8:30UTC the real ABLH is below the MUSA overlap and, consequently, as already mentioned earlier, the proposed algorithm overestimates the ABLH capturing the edge of the next layer. Starting from the 08:30UTC, the solar activity initiates and the ABLH starts to increase above altitudes exceeding the lidar overlap. Thus, from this point on, the ABLH retrieved by using lidar measurements should agree with the reference one. However, in this particular case, the lidar observations show that two separated aerosol layers are present below the reference ABLH (about 435 3 km asl) indicating probably a not well mixed ABL. The bottom layer (below about 1.5-2.0 km asl) is probably composed by local particles while the upper one contains dust of mixed dust. Clearly, in conditions like this, the ABLH retrieved by lidar measurements (independently from the specific algorithm) will underestimate the real ABLH capturing the top of the first aerosol layer and not the top of the ABL.
The other two cases refer to classic nighttime observations taken on Nov 20, 2014 and July 13, 2015. Figs. 12 and 14 report 440 the comparison of PBHL retrieved by all the considered algorithms for the selected cases, respectively. Additionally, Fig. 12 shows also the ABLH retrieved by using temperature and pressure profiles taken from a correlative radiosounding launched at CIAO observatory (green triangle). The absolute differences for these two cases are shown in Figs. 13 and 15, respectively.
The agreement among all the considered algorithms is in general good for all the three cases. For all the cases in the dataset, MIPA algorithm shows the highest accuracy with respect to the reference. This is confirmed also by the mean and median 445 values of the absolute difference with respect to the reference summarized in Tab. 4 (minimum values of both these parameters are always obtained by using the proposed algorithm).
The estimation of the ABLH is of crucial importance both for meteorological and air-pollution related applications. In this work, we proposed a new algorithm to continuously retrieve the ABLH. This approach leverages on the use of a fully image-450 based methodology (instead of analyzing the lidar observations profile by profile). The retrieval consists in applying to the image, during the pre-processing phase, morphological operators. Afterwards, an edge detection is considered. Finally, during the post-processing phase, the significant edges are extracted through a further filtering phase based on mathematical morphology and an object-based analysis. This approach has been compared with a proper benchmark consisting of state-of-the-art ABLH estimation methods, i.e., a gradient-based approach and a WCT-based method. For the latter, the filtering capabili-455 ties of the approach were pointed out. Different datasets acquired by two lidar systems located in two separated EARLINET permanent observational sites have been considered to assess the performance.
The results, relying upon several statistical indexes, put in evidence that the proposed approach is more accurate than the compared approaches belonging to the benchmark. In particular, we observed an improvement of the accuracy of about 30% (on average) with respect to the closest state-of-the-art approach (i.e., the WCT). Moreover, the outcomes obtained by the 460 MIPA are more stable than the other benchmarking methods. This can be easily pointed out by having a look at the results depicted in this paper and it has also been corroborated by calculating measures of dispersion (e.g., the standard deviation) in the statistical analysis. The last concluding remark is about the computation analysis. Despite of the proposed approach seems quite complex, it leverages on the use of very efficient filters based on mathematical morphology. The running times on large datasets (72 hours) show excellent performance from this point of view requiring just few seconds for the execution 465 of the whole signal processing chain. The bottleneck of the system turns out to be the object analysis phase. However, the computation times can be considered comparable with the other approaches proposed into the benchmark.
Finally, it is worth to be stressed that the MIPA approach does not depend on the absolute calibrated values but rather only on the correlation among adjacent lidar range bins, which is of crucial importance in order to form the image. This interesting feature could be very useful for future developments. Indeed, the MIPA approach could be easily adapted to address the task 470 of the estimation of the ABLH using other widely available and continuously acquired data, such as, ceilometer data.    Figure 2. ABLH retrieved during the 72h EARLINET exercise (July 9-12, 2012) for Potenza. Gray, red, yellow and black lines show the ABLH retrieved from 1064 nm high-resolution total attenuated backscatter time series by applying the edge detector module on full resolution data (curve a), the vertical resolution reduction module and the edge detector (curve b), the vertical resolution reduction module, the pre-processing module and the edge detector (curve c) and the whole MIPA procedure (curve d). Brown circles represent the reference ABLHs as retrieved from atmospheric temperature and pressure profiles provided by ECMWF forecasts. Table 2. Statistical analysis of the absolute differences with respect to the reference of the ABLH retrieved by applying MIPA, WCT and derivative approaches on Potenza 72h high-resolution time series of the total attenuated backscatter at 1064 nm. The mean (∆mean), the median (∆mean), the standard deviation (∆ std ), the standard error (∆ste), the minimum and the maximum of the absolute differences are given in meters. N is the number of points on which the statistics are made. The reference is assumed to be the ABLH calculated from the co-located atmospheric temperature and pressure profiles provided by ECMWF forecasts.  Figure 3. Absolute differences between the ABLHs retrieved from Potenza 72h high-resolution total attenuated backscatter time series at 1064 nm and corresponding ABLHs derived from atmospheric temperature and pressure profiles as provided by ECMWF forecasts. ABLHs retrieved by using ECMWF forecasts data are assumed to be the reference. The original time resolution of ECMWF forecasts is 1 hour. They have been interpolated at the lidar data time resolution (1 minute). Gray, red, yellow and black lines show the absolute differences between the curves a), b), c) and d) shown in Fig. 2 and the reference, respectively. Table 3. Statistical analysis of the absolute differences with respect to the reference of the ABLH retrieved by applying MIPA, WCT and derivative approaches on Evora 72h high-resolution time series of the total attenuated backscatter at 1064 nm. The mean (∆mean), the median (∆mean), the standard deviation (∆ std ), the standard error (∆ste), the minimum and the maximum of the absolute differences are given in meters. N is the number of points on which the statistics are made. The reference is assumed to be the ABLH calculated from the co-located atmospheric temperature and pressure profiles provided by GDAS forecasts.  Brown circles represent the reference ABLH as retrieved from atmospheric temperature and pressure profiles provided by ECMWF forecasts. ABLHs retrieved by using ECMWF forecasts data are assumed to be the reference.       Figure 15. Absolute differences between the ABLHs retrieved from high-resolution total attenuated backscatter time series at 1064 nm measured over Potenza on July 13, 2015 and the corresponding ABLHs derived from atmospheric temperature and pressure profiles as provided by ECMWF forecasts. ABLHs retrieved by using ECMWF forecasts data are assumed to be the reference. The original time resolution of ECMWF forecasts data is 1 hour. They have been interpolated at the lidar data time resolution (1 minute). Black, red and yellow lines show the absolute differences with respect to the reference of the ABLHs retrieved by MIPA, derivative and WCT algorithms, respectively.