ANGELINI_corrected:Layout 6 ANNALS OF GEOPHYSICS, 57, 2, 2014, A0218; doi:10.4401/ag-6408 A0218 Some remarks about lidar data preprocessing and different implementations of the gradient method for determining the aerosol layers Federico Angelini1,2,*, Gian Paolo Gobbi2 1 ENEA, Laboratorio Diagnostiche e Metrologia (UPATPRAD-DIM), Frascati (Rome), Italy 2 CNR, Istituto di Scienze dell’Atmosfera e del Clima (ISAC), Rome, Italy ABSTRACT The determination of atmospheric aerosol layers from lidar returns is possible through automated algorithms. This product is useful, for ex- ample, in monitoring the Boundary Layer Height (BLH) as well as vol- canic plumes. Aerosol layers are usually detected using the gradient method, i.e. by finding the inflection points of the range-corrected backscattered signal. These points can be either calculated as the minima of the numerical derivative of the signal, or by zero-order Digital Wavelet Transforms. Since for low signal-to-noise ratios the numerical derivative is very prone to noise-induced fluctuations, a moving average is often per- formed. We demonstrate why this procedure should be avoided in the gra- dient calculation. Finally, an alternative approach to the Digital Wavelet Transform is proposed, giving the same results but lowering the compu- tational times by about one order of magnitude. 1. Introduction During the last years both research lidar and ceilometers networks have been set up worldwide [Flentje et al. 2010b] to monitor the vertical structure of aerosol and clouds. More than 1000 ceilometers have been operating in Europe in 2012, producing a very huge amount of data. Although the raw signal contour plot for visual inspection is an important output of these systems, some quantities are calculated by auto- mated algorithms (cloud bottom, cloud coverage). Among them, a very important quantity is the Bound- ary Layer Height (BLH), often also called Mixing Layer Height or Mixing Height, used as a first approximation of atmospheric dispersion capability. On the other hand, the retrieval of the BLH is not straightforward since there are not direct methods for its measurement. A thorough discussion about the prob- lems related to the determination of the BLH can be found in Seibert et al. [2000]. Elastic lidars allow the de- termination of the BLH under some assumptions that will be discussed later on. At the basis of this method, however, there is the assumption that the aerosol is pro- duced at ground and dispersed by turbulence within the PBL. The BLH is one of the most important quantities related to the Atmospheric Boundary Layer, having many theoretical and practical applications. This quan- tity, in fact, is important in air pollution modelling, in aviation applications and evaluation of model per- formances. In spite of its importance, the definition of MLH is not straightforward. According to the COST710 Action final report, “the practical determination of the MLH, and sometimes even its definition, is not trivial, and there are many associated practical and theoretical problems”. For our purposes, the following definition, from Seibert et al. [2000] may be assumed: “The mixing height is the height of the layer adjacent to the ground over which pollutants or any constituents emitted within this layer or entrained into it become vertically dis- persed by convection or mechanical turbulence within a time scale of about an hour”. Because of the importance of the BLH, many pa- rameterization have been developed and implemented in Numerical Weather Prediction models [Arya 1981]. In fact, the turbulence in the PBL is a sub-grid process that cannot be explicitly resolved by the NWP, and must be somehow parameterized. The experimental determination of the BLH is then important for the evaluation of the performances of the different para- meterization employed in the numerical weather fore- cast models [Bei et al. 2010, Hu et al. 2010]. There are a few ways to identify the mixing layer Article history Received September 20, 2013; accepted March 25, 2014. Subject classification: Data processing, Volcano monitoring, Volcanic effects, Atmosphere: composition and structure, Instruments and techniques. from aerosol elastic backscatter signals. Some employ the temporal variance of the signal [Hennemuth and Lammert 2006] or a threshold on the backscatter coef- ficient [Melfi 1985], but the most used criterion is to search for inflection points in the backscatter coefficient profile. This method is referred to as the gradient method, proposed initially by Endlich et al. [1979], and since then widely used to detect the MLH. In fact, the other methods are not suitable for automated detection by low power systems. Moreover, while the variance method performs well under convective regimes, the detection of stable layers becomes hard, as witnessed by the fact that most works intentionally refer only to the unstable boundary layer [Hennemuth and Lam- mert 2006, Lammert and Bosenberg 2006, Pal et al. 2010]. As a consequence, the gradient method is the most used for MLH detection, and much effort is done to make it more efficient. Using the gradient method is also possible to infer the thickness of the entrainment zone [Flamant et al. 1997]. Unfortunately, however, the main drawback of this method is given by the possible presence of multiple aerosol layers, since strong gradi- ents are found at the top of each layer. Recent works [Haeffelin et al. 2011] have demonstrated, through an in depth comparison among different algorithms, that the major uncertainty comes from the attribution of the BLH to a layer than the determination of the layers it- self. The best way to automatically attribute the BLH to the right layer is still under investigation [Angelini et al. 2009, Di Giuseppe et al. 2011, Haeffelin et al. 2011, Pal et al. 2013]. It is worth to sketch quickly that the results of the gradient method differ from a system to another also depending on the wavelength employed. In fact, the relative contribution of aerosols and molecules strongly depend on the laser wavelength, as discussed for exam- ple in Haeffelin et al. [2011]. In any case, the importance of the detection of all the aerosol layers is a key step in every gradient-based algorithm and should then be per- formed at better and in the most efficient way. Furthermore, the determination of aerosol layers is also useful in monitoring volcanic plumes, although because of their sporadic nature such observations are still less systematic than estimations of the BLH. Nev- ertheless, the development of ceilometer networks worldwide can make this technique very powerful in monitoring the status of volcanoes [Sassen et al. 2007, Flentje et al. 2010a]. 2. Discussion of method and results Strictly speaking, the gradient analysis should be performed on the backscatter coefficient, which repre- sents the volume backscatter cross section of the total aerosol amount. However, this requires a calibration of the signal over the attenuated molecular backscattering profile. This is not convenient for automated retrievals and, moreover, the low power of the ceilometers cou- pled to the typical use of wavelengths in the infrared re- gion does not allow this easily. Anyway, the logarithm of the Range-Corrected Signal (RCS) is a good alterna- tive, being roughly proportional to the total backscatter coefficient, apart from the attenuation due to extinction. The lidar equation can be written in the form: Where B represents the background signal and C the lidar constant, which takes into account a variety of factors (power of the emitted pulses, efficiency of the collection and detection system), and bm, ba repre- sent the backscatter coefficient of molecules and aerosols, respectively, and T the atmospheric transmis- sivity given by: , being a the extinction coefficient, and where i represents either the molecules or the aerosols. The term O(z) should be introduced to take into account the fact that the overlap between the lidar beam and the telescope field of view is not usually uniform in space, but is in- stead a function of the distance. O(z) is typically zero in the close field, and reaches 1 at long distances when the laser enters completely in the receiver’s field of view, and the full overlap is achieved. This means that almost all li- dars are blind at very short distances (from tens to hun- dreds of meters), and hence the detection of very low layers is sometimes impossible. This represents often an insurmountable problem when shallow PBLs have to be detected, like during nighttime and under very stable at- mospheric conditions in general. Anyway, these prob- lems lie outside the purposes of this paper, and O(z) can be considered as a constant. In practice, we will consider only the lidar equation after the full overlap. The range-corrected signal is usually calculated as: , so that the lidar equation becomes: Since Tm is a very smooth function, (at 855 nm, for standard atmosphere, it equals 0.987 between ground and 4 km), no influence on gradients is expected and it can be neglected. Also Ta may be generally neglected. In fact, while the backscatter coefficient may show rapid ( ) ( ) ( ) ( ) ( ) .RCS z C z z T z T zm a m a 2 2b b= +6 @ ( ) ( ) ( ) ( ) ( ) ( ) BS z O z z C z z T z T zm a m a2 2 2b b= + +6 @ 2- ( ) i r dr z 0 ai ( )T z e= # ( ) ( ) RCS z z S z B 2= - ANGELINI AND GOBBI 2 3 variations in case of superimposed layers, the transmis- sivity depends on the integral of the extinction coeffi- cient, and usually it varies smoothly. Since we are concerned in detecting the gradients of the signal, such approach actually works well. As an example, we can consider a typical optical thickness of 0.1 at 870 nm: the transmissivity is 0.82 and hence, neglecting the correc- tion, this leads to an underestimation of about 20% at the top of the layer. Moreover, T is a monotonic func- tion and therefore neglecting it cannot introduce large errors in the detection of the inflection points. Hence, we can neglect the transmission coefficients, and we can write, assuming an exponential distribution of the mo- lecular density, and if C is expressed in opportune units: where N0 is the Loschmidt constant at standard condi- tions (i.e. the molecular density in air), H the scale height for molecular density (~ 8 km), and rep- resents the Rayleigh backscatter volume cross section. For aerosol free atmosphere, ba=0 and the equa- tion above becomes: Since ~5 10-33 m2 sr-1 at 855 nm and N0 ~2.7 10 25 m-3, the first addend on the right becomes the of the order of −15, while within the boundary layer z/H spans from 0 to 0.5. Hence, if expressed in RCS, the lidar attenuated molecular backscattering profile is approximately a straight line, with a small angular coefficient that makes ln(RCS) a good choice for gradient detection. Another interesting feature of such quantity is that the slope of the attenuated molecular backscattering profile does not depend on the wavelength of the system. On the contrary, the aerosol contribution of course does. In presence of aerosols (ba >0 and dependent on wave- length) the RCS is still greater than in the case of pure molecular echo, at least as long as the aerosol extinc- tion can be neglected. Though this representation of the lidar equation is not convenient to retrieve the aerosol backscatter and attenuation, it is useful for the application of the gra- dient method. In spite of the definition, however, it should be noted that the determination of the inflection points by gradient method can be determined in a variety of ways, i.e. the direct numerical differentiation, the Dig- ital Wavelet Transform (DWT) or even through best fit- ting procedures with Gaussian or error functions. For systems having a limited signal-to-noise ratio, the direct differentiation of the lnRCS may introduce unrealistic gradients. For this reason, spatial and tem- poral averages are introduced to improve the signal-to- noise ratio. 3. Signal averaging Both temporal and spatial averaging of the lidar signal are widely used to improve the signal-to-noise ratio. The problems related to this procedure is the lost of information of small-scale events, such turbulence or microphysics. Orlanski’s classification of atmos- pheric processes [Orlanski 1975] showed a continuous cascade of spatial and temporal scaled from planetary waves to turbulence. The processes that involve the PBL are quoted in Table 1. From Table 1 it is clear that different averaging will point out some patterns and will completely mask the processes at smaller scales. Hence, suitable averaging windows must be carefully exploited to highlight the -( ) ( ),ln ln lnRCS d d N e Ca H z 0b v X = + + r ` j d dv X r ( ) ( ) .ln ln lnRCS z d d N z H C0 v X = - + r ^ h ` j ON THE APPLICATION OF GRADIENT METHOD TO LIDAR DATA Process Temporal scale Spatial scale Scale definition Nocturnal low-level jet Inertial waves 1d - 2h 20 - 200 km Meso-b Urban effects Internal Gravity Waves 10 h - 10 min 2 - 20 km Meso-c Deep convection Short gravity waves 1h - 2 min 0.2 - 2 km Micro-a Dust devils Thermals Wakes 10 - 1 min 20 - 200 m Micro-b Plumes Roughness Turbulence 1 minute - 1 s < 20 m Micro-c Table 1. Scales of atmospheric processes involved in the PBL. Adapted from Orlanski [1975]. d dv X r processes investigated. The average can be performed either using moving average or a block average. Ad- vantages and drawbacks of these approaches will be now discussed. The running average may be employed to improve the signal-to-noise ratio of the lidar data without re- ducing the resolution. This is true only apparently, since the running average may introduce artifacts such spu- rious maxima and saturation plateau around peaks (Fig- ure 1). Of course, the points averaged may be weighted by suitable filters (i.e. binomial, Gaussian), but in this case the effectiveness of the smoothing in removing noise lowers considerably. The other procedure to increase the signal-to noise ratio consists in lowering the temporal resolution per- forming a block average over the profiles. Since the standard deviation of the background signal (SSD) can be thought as noise in lidar measurements, to under- stand the effect of temporal averaging on the spatial standard deviation of the signal between 6 and 7 km has been calculated for several temporal averaging in- tervals, spanning from 1 to 240 min. It has been found that with a 30 min-average, the SSD falls closer than 10% to its minimum value, and after 60 min, falls below 5%. These results are shown in Figure 2. After 200 min the SSD raises slightly, likely because of changes in the overall status of the atmosphere. In the spatial dimension, the signal is often smooth- ed with range-dependent windows (from 50 to 200 m) and then the direct differentiation is performed on smoothed data. Anyway, it is possible to demonstrate that the mov- ing average smoothing is completely not useful for gra- dient calculation, since the resulting data are as noisy as in the case of no smoothing. To demonstrate that, let us consider a function f (k) the ±n-point window centered around the k-th point of the discrete domain of f. The smoothed function f ’ is obtained as an average over the considered window: At the point k+1, the smoothed function is ob- tained as: Hence, the difference f ’(k) − f (k+1) becomes: The gradient is then expressed by the difference of two points, far 2n+1 points one from the other. The first effect of this operation is that only one point is used for the gradient calculation, strongly limiting the effectiveness of the aveaging. Furthermore, the points used lay as far as the window width, causing the values of the gradient to become meaningless. This is visible in Figure 3, where the red solid curve shows the gradi- ent of the function smoothed over a moving average, compared to the DWT and a third algorithm that will be described hereafter. 4. Wavelet algorithms Algorithms based on the use wavelets are used to overcome this problem, being them especially effective for noisy data. These methods employ the convolution ( ) .n f k j2 1 1 1j n n + + +=-f =( )k 1+’ / ( ) ( 1) ... ( ) ( 1 ) ... ( ) ( 1 ) ( ) ( ) . n f k n f k n f k n f k n f k n f k n n f k n f k n 2 1 1 2 1 1 1 1 + - + - + + + + - - + - - - + - + + = + - - - + + = ^ ^ h h ( ) .n f k j2 1 1 j n n + +=-f =( )k’ / ANGELINI AND GOBBI 4 Figure 1. The effect of running average on close peaks. Around 5 a spurious peak (‘ghost’) appears although the original signal showed a minimum in this zone. These effect may affect both temporal and spatial averages, for example when small clouds are detected. Figure 2. Effects of time averaging on the noise. At 30 (60) min. the 90% (95%) of the lowest standard deviation on the signal is found. The SSD is the average of seven different profiles obtained through- out the day, for different days. 5 product of the function and a Haar wavelet (also known as D2, the 0-th order Daubechies wavelet) and have been abundantly described in the literature [Brooks 2003, de Haij et al. 2007, Morille et al. 2007]. Anyway a short description of the WCT will be given. The Haar wavelet is defined by the equation: It is characterized by a step-like appearance (Fig- ure 4), and the convolution of such function with the signal shows maxima in correspondence of inflection points. The advantage of this technique is that, being based on an integral rather than a derivative, it is less sensitive to noise, and hence it is widely used to ana- lyze low power lidar profiles. The dilation of the wavelet (i.e. the size of its non null part) can be varied within the characteristic scales of the signal variations, so that the average of these multiple wavelets allow the so-called “multi scale approach” which is useful for bet- ter discriminate the noise-induced inflection points. The dilation of the specific algorithm implemented spans from 90 to 360 m. The main drawback of this method is that no determination is possible in that re- gion closer to the ground than the half dilation of the largest wavelet. For our algorithm, in particular, no layer detection is possible below 180 m, which is defi- nitely too high for most nocturnal layers. Anyway, it must be noted that this problem usually arises in paral- lel to the lack of overlap in the close field and in many cases it does not extend really the blind zone of a lidar. The DWT is then compared with a threshold, and all the maxima exceeding this threshold are stored. The same is done for the numerical derivative. In order to take into account a temporal consistency of the layers, three consecutive profiles are grouped (representing hence a lapse of time of 45 min), all their inflection points are joined together and sorted by increasing height. Then, the consecutive points whose distance ex- ceeds a season-dependent threshold are considered be- longing to different layers. Such clustering allows to detect the mean height and the thickness of each layer. This latter is intended as twice the standard deviation of a k 01 1-, ( ) , , H k a a a k otherwise 1 1 0 02 $= - Z [ \ ]] ]] ON THE APPLICATION OF GRADIENT METHOD TO LIDAR DATA Figure 3. Upper panel: Signal (green) and smoothed signals using different window sizes. Bottom panel: The gradient of the signal as cal- culated bt DWT, Derivative on Moving Average and Derivative on Block Average. In all cases, the multi-scale approach was the same, span- ning from 9 (15) to 24 points of dilation with an increment of 3 each time. Each resulting curve is the average of the six profiles obtained at different dilation, as explained for example in de Haij et al. [2007]. The vertical lines represent the top of the aerosol layers detected when the functions in the bottom panel exceed a threshold at its maxima. The threshold in this case is set at 0.05. Figure 4. The effect of the block average over the windows W1 and W2 and the DWT on the function f. the heights of the points belonging to the layer. Since clusters containing only one point are excluded, this method lets one exclude some isolated spurious layers. The average values of the DWT or the first derivative maxima relative to all the points of each layer are also computed, since those may be included among the pa- rameters to be considered for the attribution. Anyway, the computation of the DWT is quite time consuming, since the product of the wavelet for the function (even though only the window of finite values of the wavelet can be considered) must be calculated for all the points of the function. This requires a big amount of opera- tions, although this is not usually prohibitive for mod- ern computers. However, if a lot of profiles must be analyzed and a multiscale approach is adopted, a quicker algorithm is certainly desirable. Here we will illustrate a different algorithm to cal- culate the gradient, which is as robust as the DWT, but decreases the computation time of about one order of magnitude. It will be referred to as DBA, since it is based on the Derivative of Box Averaged data. Let us now consider the same window sizes used for the moving average and the multiscale approach in the DWT. If the function f is now averaged over disjoint win- dows (the spatial resolution is reduced by a factor of 2n+1), the difference of the function between two con- tiguous points becomes similar to the Wavelet trans- form. The scheme is presented in Figure 4: if the function f is averaged over the two windows W1 and W2 and then the difference is computed, the result is the same of the convolution of the wavelet function H defined with the same dilation. After the calculation of the gradient with the box average, an interpolation is necessary to retrieve the same resolution of the signal. This is also necessary to compare gradients at different dilations. Both linear in- terpolation and nearest neighbor approximation have been used, and the results do not change appreciably. At the end, as usually done with DWT, the average of the different windowed gradients is calculated. This multiscale approach helps refining the calculation of the gradient maxima, and almost the same results of the DWT are obtained. To evaluate the usefulness of the method, the al- gorithm has been implemented in Matlab environ- ment, and the computation times of both techniques have been considered. For the dilations used (15, 18, 21, 24 points) the ratio between the average computation times was about 10. This means that ‘jumping’ the moving window, especially at larger scales, does not worsen the results of the wavelet transform. The computation times of the DWT and the DBA for 8 whole days of measurements, consisting in 96 dif- ferent ceilometers profiles per day have been calculated. The ratio of the mean calculation times is about 10. If the smallest dilations are discarded, this ratio raises again, since the larger the windows the less calculations are required for the block average. 5. Conclusions The lidar Range-Corrected Signal (RCS) is widely employed in automated analysis for the detection of the atmospheric Boundary Layer Height (BLH), and, more in general, for detection of aerosol layers, in- cluding clouds. In this paper it has been demonstrated that ln(RCS) is a convenient quantity for the application of the gradient method, since in absence of aerosol it behaves as a straight line with small slope. The unhelp- fulness of moving average has been demonstrated in gradient detection, since the cancellation of most terms makes the derivative extremely noisy and poorly mean- ingful. Finally, an algorithm based on the Derivative of Block-Averaged data (DBA) for aerosol layer detection is compared with the classical Digital Wavelet Trans- form. This algorithm is built starting from the strong analogy between the derivative of block-averaged data and the Haar wavelet transform, but allows to save about the 90% of operations when detecting the in- flection points: this aspect is helpful in the automated processing of large datasets. The DBA algorithm can be applied to other kinds of analysis, wherever saving computational time is often mandatory. ANGELINI AND GOBBI 6 Figure 5. Computation times, in seconds, for multi scale DWT and DBA, using linear interpolation (green) and ‘nearest neighbor ap- proximation’ (red). The ratio between the mean computation times is 10.91 for the first case and 10.56 for the latter. This is dependent on the dilations used, since the time save for larger windows is higher than for small windows. 7 References Angelini, F., F. Barnaba, T.C. Landi, L. Caporaso and G.P. Gobbi (2009). Study of atmospheric aerosols and mixing layer by LIDAR, Radiat. Prot. Dosim., 137 (3-4), 275-279; doi:10.1093/rpd/ncp219. Arya, S.P.S. (1981). Parameterizing the height of the sta- ble atmospheric boundary layer, J. Appl. Meteorol., 20, 1192-1202. Bei, N., W. Lei, M. Zavala and V. Molina (2010). Ozone predictabilities due to meteorological uncertainties in the Mexico City basin using ensemble forecasts, Atmos. Chem. Phys., 10, 6295-6309; doi:10.5194/acp- 10-6295-2010. Brooks, I. (2003). Finding boundary layer top: application of a wavelet covariance transform to lidar backscat- ter profiles, J. Atmos. Ocean. Tech., 20, 1092-1105. de Haij, M., W. Wauben and H. Klein Baltink (2007). Continuous mixing layer height determination using the LD-40 ceilometer: a feasibility study, KNMI sci- entific report, Royal Netherlands Meteorological In- stitute, 102 pp. Di Giuseppe, F., A. Riccio, L. Caporaso, G. Bonafè, G.P. Gobbi and F. Angelini (2011). Automatic detection of atmospheric boundary layer height using ceilome- ter backscatter data assisted by a boundary layer model, Q. J. Roy. Meteor. Soc., 138 (664), 649-663; doi:10.1002/qj.964. Endlich, R.M., F. Ludwig and E.E. Uthe (1979). An au- tomated method for determining the mixing layer depth from lidar observations, Atmos. Environ., 13, 1051-1056. Flamant, C., J. Pelon, P. Flamant and P. Durand (1997). Lidar determination of the entrainment zone thickness at the top of the unstable marine atmos- pheric boundary layer, Bound.-Lay. Meteorol., 83, 247-284. Flentje, H., H. Claude, T. Elste, S. Gilge, U. Kohler, C. Plass-Dulmer, W. Steinbrecht, W. Thomas, A. Werner and W. Fricke (2010a). The Eyjafjallajokull eruption in April 2010 - detection of volcanic plume using in- situ measurements, ozone sondes and lidar-ceilome- ter profiles, Atmos. Chem. Phys., 10, 10085-10092. Flentje, H., B. Heese, J. Reichardt and W. Thomas (2010b). Aerosol profiling using the ceilometers net- work of the German Meteorological Service, Atmos. Meas. Tech. Discuss., 3, 3643-3673. Haeffelin, M., F. Angelini, Y. Morille, G. Martucci, S. Frey, G.P. Gobbi, S. Lolli, C.D. O’Dowd, L. Sauvage, I. Xueref-Rémy and B. Wastine (2011). Evaluation of mixing height retrievals from automatic profiling li- dars and ceilometers in view of future integrated networks in Europe, Bound.-Lay. Meteorol., 143 (1), 49-75. Hennemuth, B., and A. Lammert (2006). Determination of atmospheric boundary layer height from ra- diosonde and lidar backscatter, Bound.-Lay. Meteo- rol., 120, 181-200. Hu, X., F. Zhang and J.W. Nielsen Gammon (2010). En- semble based simultaneous state and parameter es- timation for treatment of mesoscale model error: A real data study, Geophys. Res. Lett., 37, L08802; doi:10.1029 /2010GL043017. Lammert, A., and J. Bosenberg (2006). Determination of the convective boundary-layer height with laser remote sensing, Bound.-Lay. Meteorol., 119, 159-170. Melfi, S.H., J.D. Spinhirne, S.H. Chou, S.P. Palm (1985). Lidar observations of vertically organized convec- tion in the planetary boundary layer over the ocean, J. Clim. Appl. Meteorol., 24, 806-821. Morille, Y., M. Haeffelin, P. Drobinski and J. Pelon (2007). STRAT: an automated algorithm to retrieve the vertical structure of the atmosphere from sin- gle-channel lidar data, J. Atmos. Ocean. Tech., 24, 761-775. Orlanski, I. (1975). A rational subdivision of scales for atmospheric processes, B. Am. Metereol. Soc., 56 (5), 527-530. Pal, S., A. Behrendt and V. Wulfmeyer (2010). Elastic- backscatter-lidar-based characterization of the con- vective boundary layer and investigation of related statistics, Annales Geophysicae, 28, 825-847. Pal, S., M. Haeffelin and E. Batchvarova (2013). Ex- ploring a geophysical process-based attribution tech- nique for the determination of the atmospheric boundary layer depth using aerosol lidar and near- surface meteorological measurements, J. Geophys. Res.-Atmos., 118, 1-19; doi:10.1002/jgrd.50710. Sassen, K., J. Zhu, P. Webley, K. Dean and P. Cobb (2007). Volcanic ash plume identification using po- larization lidar: Augustine eruption, Alaska, Geo- phys. Res. Lett., 34, L08803; doi:10.1029/2006GL0 27237. Seibert, P., F. Beyrich, S.E. Gryning, S. Joffre, A. Ras- mussen and P. Tercier (2000). Review and intercom- parison of operational methods for the determination of the mixing height, Atmos. Environ., 34, 1001- 1027. *Corresponding author: Federico Angelini, ENEA, Laboratorio Diagnostiche e Metrologia (UPATPRAD-DIM), Frascati (Rome), Italy; email: federico.angelini@enea.it. © 2014 by the Istituto Nazionale di Geofisica e Vulcanologia. All rights reserved. ON THE APPLICATION OF GRADIENT METHOD TO LIDAR DATA << /ASCII85EncodePages false /AllowTransparency false /AutoPositionEPSFiles false /AutoRotatePages /None /Binding /Left /CalGrayProfile (Dot Gain 20%) /CalRGBProfile (sRGB IEC61966-2.1) /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2) /sRGBProfile (sRGB IEC61966-2.1) /CannotEmbedFontPolicy /Warning /CompatibilityLevel 1.3 /CompressObjects /Tags /CompressPages true /ConvertImagesToIndexed true /PassThroughJPEGImages true /CreateJobTicket false /DefaultRenderingIntent /Default /DetectBlends true /DetectCurves 0.1000 /ColorConversionStrategy /LeaveColorUnchanged /DoThumbnails false /EmbedAllFonts true /EmbedOpenType false /ParseICCProfilesInComments true /EmbedJobOptions true /DSCReportingLevel 0 /EmitDSCWarnings false /EndPage -1 /ImageMemory 1048576 /LockDistillerParams true /MaxSubsetPct 100 /Optimize false /OPM 1 /ParseDSCComments true /ParseDSCCommentsForDocInfo true /PreserveCopyPage true /PreserveDICMYKValues true /PreserveEPSInfo true /PreserveFlatness true /PreserveHalftoneInfo false /PreserveOPIComments false /PreserveOverprintSettings true /StartPage 1 /SubsetFonts true /TransferFunctionInfo /Apply /UCRandBGInfo /Preserve /UsePrologue false /ColorSettingsFile (None) /AlwaysEmbed [ true /AndaleMono /Apple-Chancery /Arial-Black /Arial-BoldItalicMT /Arial-BoldMT /Arial-ItalicMT /ArialMT /CapitalsRegular /Charcoal /Chicago /ComicSansMS /ComicSansMS-Bold /Courier /Courier-Bold /CourierNewPS-BoldItalicMT /CourierNewPS-BoldMT /CourierNewPS-ItalicMT /CourierNewPSMT /GadgetRegular /Geneva /Georgia /Georgia-Bold /Georgia-BoldItalic /Georgia-Italic /Helvetica /Helvetica-Bold /HelveticaInserat-Roman /HoeflerText-Black /HoeflerText-BlackItalic /HoeflerText-Italic /HoeflerText-Ornaments /HoeflerText-Regular /Impact /Monaco /NewYork /Palatino-Bold /Palatino-BoldItalic /Palatino-Italic /Palatino-Roman /SandRegular /Skia-Regular /Symbol /TechnoRegular /TextileRegular /Times-Bold /Times-BoldItalic /Times-Italic /Times-Roman /TimesNewRomanPS-BoldItalicMT /TimesNewRomanPS-BoldMT /TimesNewRomanPS-ItalicMT /TimesNewRomanPSMT /Trebuchet-BoldItalic /TrebuchetMS /TrebuchetMS-Bold /TrebuchetMS-Italic /Verdana /Verdana-Bold /Verdana-BoldItalic /Verdana-Italic /Webdings ] /NeverEmbed [ true ] /AntiAliasColorImages false /CropColorImages true /ColorImageMinResolution 150 /ColorImageMinResolutionPolicy /OK /DownsampleColorImages true /ColorImageDownsampleType /Bicubic /ColorImageResolution 300 /ColorImageDepth -1 /ColorImageMinDownsampleDepth 1 /ColorImageDownsampleThreshold 1.10000 /EncodeColorImages true /ColorImageFilter /DCTEncode /AutoFilterColorImages true /ColorImageAutoFilterStrategy /JPEG /ColorACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /ColorImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 150 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 300 /GrayImageDepth -1 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.10000 /EncodeGrayImages true /GrayImageFilter /DCTEncode /AutoFilterGrayImages true /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /GrayImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 1200 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.08250 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects false /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile (None) /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName (http://www.color.org) /PDFXTrapped /Unknown /CreateJDFFile false /SyntheticBoldness 1.000000 /Description << /ENU (Use these settings to create PDF documents with higher image resolution for high quality pre-press printing. The PDF documents can be opened with Acrobat and Reader 5.0 and later. These settings require font embedding.) /JPN /FRA /DEU /PTB /DAN /NLD /ESP /SUO /NOR /SVE /KOR /CHS /CHT /ITA >> >> setdistillerparams << /HWResolution [2400 2400] /PageSize [595.000 842.000] >> setpagedevice