Linear Estimation of Particle Bulk Parameters from Multi-wavelength Lidar Measurements

An algorithm for linear estimation of aerosol bulk properties such as particle volume, effective radius and complex refractive index from multiwavelength lidar measurements is presented. The approach uses the fact that the total aerosol concentration can well be approximated as a linear combination of aerosol characteristics measured by multi-wavelength lidar. Therefore, the aerosol concentration can be estimated from lidar measurements without the need to derive the size distribution, which entails more sophisticated procedures. The definition of the coefficients required for the linear estimates is based on an expansion of the particle size distribution in terms of the measurement kernels. Once the coefficients are established, the approach permits fast retrieval of aerosol bulk properties when compared with the full regularization technique. In addition, the straightforward estimation of bulk properties stabilizes the inversion making it more resistant to noise in the optical data. Numerical tests demonstrate that for data sets containing three aerosol backscattering and two extinction coefficients (so called 3β + 2α) the uncertainties in the retrieval of particle volume and surface area are below 45 % when input data random uncertainties are below 20 %. Moreover, using linear estimates allows reliable retrievals even when the number of input data is reduced. To evaluate the approach, the results obtained using this technique are compared with those based on the previously developed full inversion scheme that relies on the regularization procedure. Both techniques were applied to the data measured by multiwavelength lidar at NASA/GSFC. The results obtained with both methods using the same observations are in good agreement. At the same time, the high speed of the retrieval using linear estimates makes the method preferable for generating aerosol information from extended lidar observations. To demonstrate the efficiency of the method, an extended time series of observations acquired in Turkey in May 2010 was processed using the linear estimates technique permitting, for what we believe to be the first time, temporal-height distributions of particle parameters.


Introduction
Theoretical and experimental studies of the last decade have demonstrated that multiwavelength (MW) Raman lidars based on a tripled Nd:YAG laser are able to provide the height distribution of particle physical parameters, such as radius, concentration and complex refractive index (Ansmann and Müller, 2005). Moreover, up to a certain limit such systems can reproduce the main features of the particle size distribution in the 0.075-10 µm radii range. To invert the aerosol extinction α and backscattering β coefficients measured at multiple wavelengths to particle parameters, numerous possible approaches have been considered but for routine processing of lidar measurements the inversion with regularization is now the most commonly used (see Müller et al., 1999;Veselovskii et al., 2002Veselovskii et al., , 2004Veselovskii et al., , 2009Kolgotin and Müller, 2008, and references therein). In order to adequately address the fundamental non-uniqueness of the 1136 I. Veselovskii et al.: Linear estimation of particle bulk parameters lidar data interpretation, a family of solutions is generated in the framework of this approach. Specifically, a series of solutions is generated using different initial guesses, different aerosol assumptions and different settings of a priori constraints. Each single solution is obtained using the regularization technique. Then the individual solutions corresponding to the smallest residuals are averaged and the result of the averaging is taken as the best estimate of the aerosol properties.
This approach has demonstrated the possibility of providing rather adequate retrievals of aerosol properties. However it is quite time-consuming, a fact that becomes an issue when large volumes of data need to be analyzed, as for example from an air-or space-borne lidar system. Installation of MW lidars on air or space-borne platforms poses another problem: the retrieval algorithm should be more tolerant to noise in the input data since reasonable averaging times are likely to be smaller for moving lidar systems. And finally, in the regularization approach described in (Müller et al., 1999;Veselovskii et al., 2002) at least five input optical data (three backscatterings and two extinctions, so called 3β + 2α) are needed to retrieve the particle size distribution (PSD), but in many applications it would be highly desirable to decrease the number of optical channels. So, the development of the approach permitting the reliable estimation of particle bulk properties such as volume, surface density and effective radius from a reduced number of optical channels would be an important improvement. One way to assess this possibility is to attempt to approximate the bulk properties by a linear combination of the input optical data (extinction and backscattering). The corresponding weight coefficients can be determined by expanding the PSD in terms of the measurement kernels (Twomey, 1977). Thomason and Osborn (1992) used this approach to estimate aerosol mass with the multiwavelength SAGE II extinction kernels. The interpretation of lidar measurements using the linear estimate techniques was explored in early studies by Chaikovskii and Shcherbakov (1985). The potential of this approach for treating elastic-Raman multiwavelength lidar measurements was studied by Donovan and Carswell (1997) under assumption of known refractive index. The technique was further explored in recent publications (De Graaf et al., 2009 where different aerosol models were used to invert optical data without prior information about the particle refractive index. In this paper we propose a modified technique, which here and below we refer to as "linear estimation" (LE). In difference with the mentioned above approaches the complex refractive index is derived as a part of retrieval procedure together with the bulk aerosol characteristics. In addition, in order to improve stability of solution we provide not a single solution but a family of solutions closely reproducing the measurements. To validate LE, we apply it and the full inversion algorithm (Veselovskii et al., 2009) to the same data and compare the results. Finally, we apply LE to an extended series of lidar measurements to evaluate heighttemporal variations of the particle bulk parameters.

Algorithm description
The aerosol extinction (α) and backscattering coefficients (β) are related to the particle volume size distribution v(r) via integral equations as follows: Index l labels the type of optical data (i = α, β) and wavelengths λ k ; K l (m, r) are the volume kernels (VK) depending on the complex refractive index m = m R − i · m I and particle radius r ∈[r min ,r max ]. To get kernels K l (m, r) in our study the Mie computations are used, thus the particles are assumed to be spherical. This approach can also be generalized to treat the particles of irregular shape, by using the kernels corresponding the ensemble of randomly oriented spheroids (Mishchenko et al., 2000;Dubovik et al., 2006). In vectormatrix form Eq. (1) can be rewritten as: Here v is the column vector with elements v k corresponding to the particle volume inside radii interval [r k , r k+1 ] and K is the matrix containing the discretized kernels as rows. The volume distribution v(r) can be expanded in terms of the kernels of Eq. (1), as prescribed by Twomey (1977). Such an expansion assumes that vector v corresponding to v(r) can be presented as a combination of the matrix K rows, i.e.
where v g is the projection of the volume distribution on the measurement kernels, while v ⊥ is the residual -the part of volume distribution orthogonal to these kernels (Kv ⊥ = 0) and x j are the weight coefficients of expansion. Using Eq.
(2) can be rewritten as: Then, the vector x of the expansion coefficients can be found as: and the volume distribution projected on the kernels is The residual then can be written as: Equations (6)-(7) can now be used to evaluate the linear estimations of the unmeasured aerosol characteristics. If there is one or several aerosol characteristics p i (i = 1, . . . , N p ) that are not measured but are needed to be estimated using measurements g, the dependence of p i on the size distribution can be described as: Here the elements p i of vector p are the unmeasured aerosol characteristics, and P is the matrix of the corresponding coefficients. Taking into account Eq. (3) the vector p can be expressed as: Here p g represents the vector of projections of characteristics p i on the measured set g and p ⊥ represents the vector of characteristics p i on the null-space v ⊥ . In other words, p g can be estimated from g, while the measurements provide no information about p ⊥ . Using Eqs. (6) and (7) the matrices of coefficients F and D ⊥ can be expressed as The elements of matrix F can be computed and stored in the look-up tables making computations of p g very fast. The residual term p ⊥ cannot be measured with the available set of observations g, but can be estimated from numerical modeling for typical situations. The situation is particularly favorable when particle bulk property p (for example, volume, surface, number density) needs to be estimated (Donovan and Carswell, 1997). In this case the matrix P contains the weight coefficients for different integral properties as rows. For example, for volume (i = 1) P 1k = 1, for surface (i = 2) P 2k = 3 r k and for number density (i = 3) P 3k = 3 4π r 3 k . In such case, the existence of the zero space v ⊥ does not have much importance and the residual p ⊥ is generally expected to be small, because the observations g are known to be strongly sensitive to aerosol total concentrations, while being less sensitive to the details of the size distribution. Thus, the projection p g can be estimated quickly from the observations g without calculating the full size distribution v, i.e. without performing a full inversion of Eq. (2). This is a significant advantage of the using the linear estimates p g as compared to the more conventional approach which yields PSD (e.g. Müller et al., 1999;Veselovskii et al., 2002Veselovskii et al., , 2004Veselovskii et al., , 2009Kolgotin and Müller, 2008). Indeed, the calculation of p g is fast and defining the coefficients F is straightforward. Following Eq. (10) the calculation of F involves the inversion of matrix KK T . In principle this operation can be ambiguous if KK T is ill-conditioned. However, in particular cases when only a very few measured characteristics g i are used (for example, in our case we use maximum 5 different observations), the matrix KK T has small dimension (maximum 5 × 5) and in the case such as here that each of the five measured characteristics g i is quite different, is well-posed and can be inverted exactly. By contrast, the conventional approaches that provide the size distribution v and must solve Eq. (2) may face significant difficulties. For example, the least square solution can of Eq. (2) is: Here the matrix K T K has to be inverted. This matrix has dimension N v × N v , where N v is the dimension of size distribution v. Therefore, matrix K T K has significantly larger dimension than matrix KK T . For example, the aerosol retrievals from sun-photometer observations discussed by Dubovik and King (2000) use N v = 22. In such situations K T K is known to be ill-conditioned and the inversion of this matrix becomes ambiguous. Therefore in many practical applications different types of constraints can be used to achieve unique and stable solution of Eq.
(2). For example, it can be constrained based on smoothness of the solution as suggested by Phillips (1962) and Twomey (1977). However, the use of smoothness or other a priori constraints may require rather sophisticated developments (see discussion by Dubovik, 2004). We should mention also that in most of algorithms for lidar data inversion (e.g. Müller et al., 1999;Veselovskii et al., 2002Veselovskii et al., , 2004Veselovskii et al., , 2009Kolgotin and Müller, 2008), the size distribution is described by a much smaller number of parameters than 22 (typically N v = 5-7). This reduces the difficulties associated with the inversion of matrix K T K, however this also may lead to the introduction of additional errors in the algorithm, since some features of the size distribution are neglected. In this respect the coefficients F can be calculated using the detailed size distribution with very large N v , since matrix KK T has the same dimension as number of rows of K. Correspondingly, coefficients F can be always found accurately. In our computations we normally use N v = 100 radii logarithmically distributed inside the inversion interval. The measured optical data g * contain the error g Thus the uncertainty of particle parameters estimation is: Then, if the measurement errors g are random, unbiased (i.e. g = 0) and have covariance matrix g T g = C g , the corresponding covariance matrix of retrieval uncertainties p can be written as www.atmos-meas-tech.net/5/1135/2012/

I. Veselovskii et al.: Linear estimation of particle bulk parameters
Here, term C random p represents the contribution of the random measurement errors g to C p and C systematic p represents the non-random part of the errors appearing due to existence of v ⊥ which is orthogonal to the kernels K. These terms can be expressed as: The first term can be calculated using the covariance matrix of the measurements C g and the second term can be estimated using a priori estimates of v. For example, our modeling experiments in the next section show that p ⊥ have generally rather small values and, therefore, the accuracy of the retrieval of p performed using the "projection" p g (see Eq. 9) is acceptable in the majority of cases. All equations given above are written with the assumption that the matrix K in Eq. (2) and, therefore, matrix F in Eq. (10) are known accurately. However, this is not really the case since K depends on the complex refractive index. Therefore, if the actual value of the complex refractive index m(λ) is not known and we use an estimatem(λ), instead of Eq. (10) we have: or it can be written in the same manner as Eq. (9) Correspondingly, if we use an estimate of the complex refractive index, the estimate p * g = Fg * , should be replaced bỹ p * g =Fg * . The corresponding uncertainties of the retrieval can be estimated as: where F = F −F. The covariance matrix of uncertainties is: Thus, if we compare Eqs. (15) and (21), the only difference is that Eq. (21) uses matrix of coefficientsF instead of F. It is reasonable to expect that the magnitudes of elements ofF are close to those of F since the optical characteristics generally do not exhibit very high sensitivity to variations of m(λ). Therefore one can expect that the components of C random p given by Eq. (21) will have magnitudes close to those given by Eq. (15).
By contrast C  (20) still retains the second term that represents the systematic bias: Thus, the use of estimatem(λ) that is different from the actual value of complex refractive index m(λ) always leads to the appearance of a systematic error term C systematic p . Moreover, it is rather clear that C systematic p may easily dominate over C random p in Eq. (23). This becomes rather obvious when using g = Kv in Eq. (23): Here the first term contains C g = g T g -the covariance matrix of errors g of g and the second term contains the matrix gg T . Since the magnitude of errors g are generally much smaller than the magnitudes of measurements g, the elements of C g are much smaller than gg T . Therefore, even if F = F −F is not very significant the magnitude of the second term of Eq. (23) is likely to remain considerable.
Thus, if the sensitivity of g * to m(λ) is high, the errors due to wrongly chosen m(λ) can be much higher than errors due to measurement uncertainties and existence of null-space. At the same time if the sensitivity of g * to m(λ) is high, and set of observations (g * ) T = g * 1 ; g * 2 ; ...; g * N 0 is quite representative we can attempt to estimate m(λ) from available observations. The input optical data (backscatters and extinctions) are themselves the particle properties and each g * j can be recalculated back from the rest of N 0 -1data using Eq. (9), as suggested in (De Graaf et al., 2009. By doing so for each optical data, we get N 0 estimates ofg j that we compare with the observations g * j . It should be mentioned that we can not make these estimates using all N 0 optical data, because in that caseg j and g * j coincide. If the sensitivity of g * to m(λ) is high the magnitudes of errors g j =g j −g * j should strongly depend on the assumed value of m(λ). Correspondingly, if we obtained the set of estimatesg j (m) using different assumed values of m(λ) we can attempt to estimate m(λ) by searching for the smallest errors g j =g j − g * j , for example by searching for the minimum of the following quadratic form: If the sensitivity of g * to m(λ) is high this form should have a well-defined minimum and m(λ) can be estimated using the available measurements g * . In our study we assume that the errors of the measurements are the same for all channels, and refractive index is spectrally independent inside the spectral range that is considered. Then the refractive index is found from the minimum of discrepancy: Since there is no a priori knowledge about the particle size distribution and refractive index we find the discrepancy ρ for all predefined values of r min , r max , lying in the interval 0.075 µm-10 µm, and for the set of values m R and m I from respective intervals 1.35-1.65 and 0.00-0.03 just as we did in our regularization algorithm (Veselovskii et al., 2002). Normally the total number of predefined combinations does not exceed N T = 3000. Based on our previous experience (Veselovskii et al., 2002) we prefer to average the solutions near the minimum of discrepancy rather than take a single solution. Such an averaging procedure stabilizes the inversion.
To choose the averaging interval the solutions are ranged in accordance with their discrepancy from minimum ρ min to maximum ρ max , normally 1 % of solutions are averaged. The retrieval for each vertical bin was done independently. The estimation of the refractive index from the minimization of ρ in Eq. (26) is illustrated by Fig. 1. The discrepancy ρ and the uncertainty of the volume retrieval ε V are given for different assumptions of the value of m R . Synthetic input data were generated assuming a log-normal aerosol distribution dV (r) d ln r with modal radius r 0 =1 µm and variance 0.4; the model refractive index is m = 1.5-i0.005. Both ρ and ε V have minima at m R = 1.5 corresponding to the "true" value of m R , thus the minimization of the discrepancy minimizes the uncertainty of the particle volume estimation. Similar plots can be provided also for the imaginary part of the refractive index. The presence of noise g in the input data naturally increases the minimum value of discrepancy ρ that can be achieved and thus the uncertainty of the refractive index retrieval. The actual increase of the retrieval uncertainly also depends on the particle size distribution and specific realization of errors g in the optical data. To evaluate the corresponding uncertainties of the estimation of particle parameters numerical simulations using different types of PSD and different input errors can be performed as we will illustrate in the next section.
Thus the main difference of described in this section algorithm from the approach presented previously by Donovan and Carswell (1997);De Graaf et al. (2009, is that we consider not a single solution but a family of linear solutions corresponding different inversion intervals r min , r max and different complex refractive indices. The average of solutions in the vicinity of the minimum of discrepancy (Eq. 26) is considered as most probable estimate of particle parameters.

Estimation of retrieval uncertainties
Numerical simulation is used here to test the algorithm and to estimate the retrieval uncertainties. In these simulations we used synthetic input optical data assuming a log-normal aerosol distribution dV (r) d ln r with r 0 = 0.2, and 2 µm, which are typical values for the fine and the coarse mode particles (e.g. see Dubovik et al., 2002), and the variance in all cases is 0.4. As discussed in the previous section, one of the principal questions in the application of the LE technique is the estimation of the residual v ⊥ in Eq. (3). This residual will depend on the PSD and the refractive index. The computations performed demonstrate that for all values of m considered here, the residual v ⊥ is below 4 % and 15 % for r 0 = 0.2, and Fig. 2. Cumulative probability of uncertainty ε V of volume density retrieval from 3β + 2α data with input errors ε = 10 %, 20 %, 30 %. Simulation was performed for r 0 = 0.2 µm using volume kernels.
2 µm respectively. So the existence of a null-space does not present a serious limitation to the LE technique for typical atmospheric aerosols.
To evaluate the effect of input uncertainties, the random errors in the range of [0, ±ε] were added to the data and from these distorted optical data, the particle parameters were retrieved. We assume that the uncertainties in all measurement channels are equivalent so that all the diagonal elements of the error covariance matrix C g = g T g are the same. The procedure was repeated 1000 times allowing robust statistics to be gathered. The retrieval uncertainties are presented in the form of probability distributions such as shown in Fig. 2 where a typical cumulative probability of volume density uncertainty ε V is shown. For every value of ε V the plot gives the probability that the retrieval uncertainty is below this value. For example, from the plots in Fig. 2 we can conclude that in 90 % of the cases the spread in the values of the volume estimation is below ε V ≈ 20 %, 35 %, 50 % for input errors ε = 10 %, 20 %, 30 % respectively. We take these values to represent the uncertainty in the retrieval. Thus the uncertainty rises approximately linearly with ε and the method can provide reasonable estimations even for 30 % input errors.
The results shown in Fig. 2 were obtained using volume kernels (VK) of Eq. (1), but Eq. (1) can be written also using other types of kernels corresponding to the number dN dr , surface dS dr or volume density size distribution dV d ln r in logarithmic space. All these kernels (henceforth referred to as NK, SK and VLK ) can also be used in retrievals. Donovan and Carswell (1997), reported that in their approach for the retrieval of surface density the surface kernels were preferable, while for volume retrieval the volume kernels were better suited. In our study, we also tested different types of kernels. We did not notice a significant difference between these kernels for small particles, but for particles in the coarse mode the volume kernels provided slightly better estimations of all parameters. The difference with the results of (Donovan and Carswell, 1997) may be due to the optimization of inversion intervals in our algorithm. All retrievals presented below were obtained with the volume kernels.
The simulation results are summarized in Table 1 showing the uncertainty of volume (ε V ), surface (ε S ), number (ε N ) density, effective radius (ε R eff ), and real part of refractive index (ε m R ) retrieval (taken at 90 % probability level) for input random uncertainties of ε = 0, 10 %, 20 %. The effective radius was estimated from the ratio of volume and surface density: r eff = 3 V S . The results are given for r 0 = 0.2 µm and 2 µm to separately characterize uncertainties for small and big particles. In the absence of input errors, the uncertainties of the retrieval are due to the null-space and the unknown value of the refractive index as follows from Eq. (22). Minimization of discrepancy (Eq. 26) keeps the uncertainty of the volume estimation below 5 % for small particle sizes characteristic of the fine mode and below 15 % for particles with sizes more consistent with the coarse mode particles.
From the results shown in Table 1, several conclusions can be made. First of all, the retrieval is stable for both small and big particles and an uncertainty of volume estimation below 45 % can be obtained even for 20 % input errors. For small particles the uncertainties of surface, volume and effective radius estimation are close, while for big particles the surface density is the most stable parameter in retrieval where the corresponding uncertainty of ε S is less than 30 % even for 20 % input errors. The most unstable parameter in the retrieval is the number density where the corresponding uncertainty for particles with r 0 = 2 µm is above 100 % for 20 % input errors. The real part of particle refractive index can be retrieved more accurately for big particles, where the corresponding uncertainty of the real part is below ±0.04 for ε = 20 %, while for small particles this uncertainty increases to ±0.07.
In our retrievals we considered the full data set 3β + 2α and the reduced one 3β + 1α, where extinction at 532 nm was removed. The important finding is that the uncertainties in the estimates of particle parameters from 3β + 1α data in most cases did not exceed the corresponding values for 3β + 2α, thus it seems apparent that the number of input data can be decreased when only particle bulk properties are desired. Evaluation of retrieval uncertainties for different combinations of the optical data and different particles characteristics is in our plans but beyond the scope of present paper.
The imaginary part of the refractive index is one of the most difficult parameters to estimate from multi-wavelength lidar as the kernels are not very sensitive to changes in the value of m I . As already mentioned, the inverse problem (1) is strongly underdetermined, so the solution depends on the constraints used, in particular on the range of refractive index values considered during the minimization of the discrepancy (Eq. 26). To evaluate the influence of the range of m I considered on the retrievals, we performed simulations for three intervals: 0 < m I < 0.01, 0 < m I < 0.02 and 0 < m I < 0.03 assuming 10 % errors in input data and model value m = 1.5-i0.005. Computations performed for particles with r 0 = 2 µm show that the uncertainties in the estimate of m I for these intervals are 50 %, 100 % and 140 % respectively. Hence, for reasonable estimation of the imaginary part of the refractive index it is very desirable to have a priori information about the aerosol type to constrain the range of m I that is considered.
To illustrate the influence of m I on the estimation of other parameters, Fig. 3 shows the cumulative probability plots for the volume retrieval of particles with r 0 = 0.2 µm and 2 µm using the three ranges of m I mentioned above. For particles in the coarse mode the uncertainty in ε V increases from 25 % to 35 % when maximal value of m I rises from 0.01 to 0.03. For particles in the fine mode the retrieval is essentially insensitive to the range of m I considered. Thus, in spite of the ambiguity in the retrieval of the imaginary part, the uncertainty in m I has little influence on the estimation of the other parameters.

Comparison with regularization retrievals
To validate the approach described in Sect. 2, the linear estimation (LE) and regularization (Veselovskii et al., 2002) algorithms were applied to the same experimental data obtained by multiwavelength Raman lidar at NASA/GSFC in Greenbelt, MD during August-September 2006 (Veselovskii et al., 2009). The lidar is based on a tripled Nd:YAG laser and provided three particle backscattering and two extinction coefficients. The retrieval of particle microphysical parameters from these 3β + 2α data using inversion with regularization was discussed in our earlier publication, where good agreement between AERONET and lidar observations was reported (Veselovskii et al., 2009). Figure 4 shows the vertical profiles of aerosol backscattering and extinction coefficients measured at 355, 532 and 1064 nm wavelengths on 15 August 2006. The backscattering shows a maximum at 1250 m and a secondary maximum at 1900 m. The particle size distribution for this day was represented mainly by the fine mode and the uncertainties of the optical data measurements were estimated to be below 10 % (Veselovskii et al., 2009). The vertical profiles of volume density, effective radius and real part of refractive index obtained with the regularization and LE approaches are shown in Fig. 5 where the results obtained with both techniques are similar. The volume density profile has a similar shape as the particle backscattering, meaning that the particle radius and refractive index do not change significantly with height. The retrieved effective radius, as shown in Fig. 5b, is about 0.22 ± 0.055 µm for all heights, the uncertainty of retrieval is estimated from the results of numerical simulations summarized in Table 1. It should be noted that the vertical profile of the effective radius obtained with LE oscillates less than the profile obtained with regularization, suggesting a more stable inversion. The refractive indices retrieved with both techniques agree reasonably well. The real part of the refractive index slightly rises with height from 1.37 ± 0.05 to 1.43 ± 0.05, the imaginary part m I is below 0.005 for all heights.
As discussed in the previous section, the number of input optical data can be reduced when only bulk particle properties are desired. To test this claim, we also performed the inversion using the reduced set of optical data given by 3β + 1α, where extinction at 532 nm is removed. The corresponding results are also shown in Fig. 5. The inversion using either the full (3β + 2α) or reduced (3β + 1α) data sets leads to similar results, supporting the conclusions made from the numerical simulations. This comparison of the regularization and LE approaches illustrates that the LE technique can provide trustworthy estimations of particle parameters. At the same time, the high speed of the retrieval using linear estimates makes the method preferable for generating bulk aerosol information from long-term series of lidar observations.

Inversion of long-term series of multiwavelength lidar observations
To test the retrieval of time-sequences of particle parameters we used data from the multiwavelength Raman lidar at TUBITAK Research Center located in the vicinity of Istanbul, Turkey. The lidar is based on a frequency-tripled Quantel Brilliant B Nd:YAG laser with 10 Hz repetition rate. The pulse energies at λ = 355, 532 and 1064 nm are 200, 250 and Vertical profiles of (a) particle volume density, (b) effective radius, (c) real part of refractive index retrieved with LE approach from 3β + 2α and 3β + 1α data and with regularization approach from 3β + 2α data.
300 mJ, respectively. The backscattered light is collected by a 40-cm aperture Newtonian telescope inclined so that the elevation angle is 40 degrees to the horizontal. The outputs of the detectors are recorded at 7.5 m range resolution using a Licel data acquisition system that incorporates both analog and photon-counting electronics. The system is able to monitor backscattering at 355, 532, 1064 nm, Raman nitrogen signals at 387, 608 nm and Raman water vapor signal at 408 nm. A polarizing beamsplitter cube in the 355 nm channel allows simultaneous monitoring of co-and cross-polarized components of backscattered radiation. The particle depolarization ratio was calculated from the ratio of co-and cross polarized components of the particle backscattering coefficients. For the calibration of depolarization measurements the molecular depolarization ratio in an aerosol-free region was used.
In each file 3000 laser pulses were accumulated, thus the temporal resolution of the measurements is 5 min. The measurements were performed during May 2010 when weak ash layers from the Eyjafjallajökull volcanic eruption periodically reached Turkey. The temporal evolution of the particle depolarization ratio δ p at 355 nm during the night of 20-21 May is shown in Fig. 6. The highly depolarizing volcanic layer with maximum particle depolarization ratio of approximately δ p = 20 ± 5 % appears at 22:00 UTC and is observed for about a period of approximately four hours at 2-3 km heights. During the same time but for altitudes below 2 km, the particle depolarization ratio did not exceed 5 %. The depolarization ratio in the ash layer is lower than the values of δ p ∼ 40 % that were observed over Northern Europe (Ansmann et al., 2010), implying that the ash may have been mixed with more locally produced aerosols. The analysis of the meteorological situation for this day and the optical properties of ash particles is given in Papayannis et al. (2011). The same figure shows the temporal variation of the aerosol extinction at 355 nm. The extinction is calculated from the Raman nitrogen signal (Ansmann et al., 1992). To decrease the uncertainty we reduced the height resolution up to 200 m and introduced a 3 point sliding average in the temporal domain. We estimate the uncertainty of the particle extinction and backscattering α 355 , β 355 , β 532 , β 1064 calculation at the heights of interest (below 2.5 km) to be less than 10 %. Comparing Fig. 6a and b we conclude that the optical depth of the ash layer is quite low (below 0.05). TheÅngström exponent calculated from the extinction coefficients at 355 and 532 nm was about 1.8 in 0.8-2 km height range, meaning that the particles were relatively small.
To retrieve the time-sequences of particle parameters, we used the 3β + 1α data set, because the uncertainty of extinction at 532 nm was too high for the chosen temporal resolution. It also allowed us to test the ability of a reduced dataset to provide useful time series results. Figure 7a shows the time-height distribution of the particle volume density, which is similar to the time-height extinction distribution in Fig. 6b. The region of enhanced volume density is contoured and shown also in Fig. 6a in order to illustrate that it coincides well with the region of enhanced particle extinction. This implies that the particle size and refractive index did not vary significantly in this region. In the color maps in Fig. 7b, c showing effective radius and the real part of refractive index, the regions where the particle extinction is low are removed, because no reliable retrieval could be performed there. The particle effective radius is about 0.22 µm in 0.8-2 km height range and it does not vary significantly over the night. Some increase of effective radius is observed near the ash plume. The retrievals inside the ash layer should be taken with care because ash particles are of irregular shape and the retrieval is based on Mie kernels, assuming a spherical particle shape, that introduces significant uncertainties. In particular the real part of refractive index is significantly underestimated with this approach . For the treatment of non-spherical particles, the kernels corresponding to randomly oriented spheroids (Dubovik et al., 2006) can be implemented as previously shown , but that effort goes beyond the scope of present paper.
The real part of refractive index shown in Fig. 7c varies in the range of 1.39-1.45. The marked region is characterized by a value m R ≈ 1.4, indicating that the aerosol contains a significant amount of water. At low altitudes after midnight some enhancement of m R up to 1.45 is observed. As mentioned, above 2 km the real part of the refractive index can be underestimated due to particle non-sphericity. The imaginary part of refractive index was estimated as 0.006 ± 0.003. The enhancement of m I up to 0.01 was observed inside the ash layer, but again, for accurate quantification of m I in this layer the spheroidal model should be used. Thus obtained results look reasonable and demonstrate that the use of linear estimates makes possible the fast retrieval of time-height distributions of particle parameters from extended lidar observations. The inversion of optical data to the aerosol parameters shown in Fig. 7 took approximately 5 min using a standard laptop computer illustrating the potential of the LE technique for processing large volumes of the MW lidar data.

Conclusions
An algorithm for the linear estimation of aerosol bulk properties such as particle volume and complex refractive index from multiwavelength lidar measurements is presented. The particle concentration is estimated from a linear combination of aerosol backscattering and extinction coefficients measured by multi-wavelength lidar while avoiding the retrieval of the particle size distribution. This approach is shown to both increase the speed and stability of the inversion. The definition of the coefficients required for the linear estimates is based on an expansion of the particle size distribution in terms of the measurement kernels. Once the coefficients for the linear estimates are established, the approach allows very fast retrieval of aerosol bulk properties. In addition, the straightforward estimation of bulk properties stabilizes the inversion making it more resistant to noise in the optical data: the retrieval does not fail even for input random uncertainties as large as 30 %. The uncertainties of the retrieval derived from numerical simulations are close to the values reported previously for the full inversion scheme that was used to derive the entire family of solutions using the regularization procedure (Veselovskii et al., 2002(Veselovskii et al., , 2004. The application of both techniques to the same lidar measurements did not reveal significant differences in the results of the two retrieval approaches.
An important finding of this study is that it is feasible to reduce the number of input optical characteristics and still retrieve useful bulk aerosol properties. A comparison of inversions using 3β + 2α and 3β + 1α data demonstrates that excluding particle extinction at 532 nm does not significantly degrade the retrieval. At the same time, removing extinction at 355 nm enhances uncertainties of retrieval The high speed of the retrieval using linear estimates makes the method preferable for generating aerosol information from long-term series of lidar observations. To demonstrate the efficiency of the method long-term series of aerosol physical properties derived from lidar observations performed in Turkey in May 2010 were processed. As a result, the multi-wavelength lidar data, for the first time, were inverted into time-height distributions of particle parameters. We should mention though that the algorithm studied here should not be considered as a replacement for the full inversion (regularization) approach, because in many applications 3β + 2α data exist and the retrieval of PSD is critical. However, if the data set is reduced, the current work demonstrates clearly that useful physical information may still be retrievable.