Software for retrieval of aerosol particle size distribution from multiwavelength lidar signals

Computer Physics Communications 199 (2016) 53–60

Contents lists available at ScienceDirect

Computer Physics Communications journal homepage: www.elsevier.com/locate/cpc

Software for retrieval of aerosol particle size distribution from multiwavelength lidar signals✩ S. Sitarek a,∗ , T. Stacewicz b , M. Posyniak c a

Maksymilian Pluta Institute of Applied Optics, ul. Kamionkowska 18. 03-805 Warsaw, Poland

b

Institute of Experimental Physics, Faculty of Physics, University of Warsaw, ul. Pasteura 5, 02-093 Warsaw, Poland

c

Institute of Geophysics, Polish Academy of Sciences, ul. Księcia Janusza 64, 01-452 Warsaw, Poland

article

info

Article history: Received 2 February 2015 Accepted 15 August 2015 Available online 26 October 2015 Keywords: Lidar Remote sensing Atmospheric aerosol Particle size distribution

abstract Software to retrieve profiles of aerosol particle size distribution (APSD) from multiwavelength lidar signals is presented. The approach consists in direct fit of artificial signal generated using predefined distribution to the experimental signals. Combination of two lognormal functions with a few free parameters is applied for the predefined APSD. The minimization technique allows finding lognormal function parameters which provide the best fit. The approach was tested on the experimental signals registered at 1064, 532 and 355 nm. The software is designated for processing on PCs. The computation time was about several minutes. Program summary Program title: APSD_Software Catalogue identifier: AEXV_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEXV_v1_0.html Program obtainable from: CPC Program Library, Queen’s University, Belfast, N. Ireland Licensing provisions: Open Licence No. of lines in distributed program, including test data, etc.: 12813 No. of bytes in distributed program, including test data, etc.: 878099 Distribution format: tar.gz Programming language: Delphi 2010. Computer: PC. Operating system: Windows XP, Vista, 7, 8, 8.1. RAM: 37MB Classification: 13. Nature of problem: Aerosol Particle Size Distribution (APSD) is a very significant function for evaluation of atmospheric optical properties. Remote determination of APSD might be performed with multiwavelength lidar. Retrieval of profiles of APSD from multiwavelength lidar signals is an example of ill-posed problem in the atmospheric physics (in the sense of Jacques Hadamard). Solution method: The approach consists in direct fit of artificial signals to the experimental signals. The artificial signals are generated using predefined aerosol particle distribution, Combination of two lognormal functions with a few free parameters is applied for the predefined APSD. The minimization technique used in the software allows finding lognormal function parameters which provide the best fit.

✩ This paper and its associated computer program are available via the Computer Physics Communication homepage on ScienceDirect (http://www.sciencedirect.com/ science/journal/00104655). ∗ Corresponding author. E-mail addresses: [email protected] (S. Sitarek), [email protected] (T. Stacewicz), [email protected] (M. Posyniak).

http://dx.doi.org/10.1016/j.cpc.2015.08.024 0010-4655/© 2015 Elsevier B.V. All rights reserved.

54

S. Sitarek et al. / Computer Physics Communications 199 (2016) 53–60

Additional comments: The input signal should be located into the SIGNALS directory which is automatically created by the software. The results are presented in the main form of the software. They are also saved into the same directory from which it was loaded input file. They are saved as a BitMap as well as the ASCII files. The software returns two main sets of files: first one ‘‘APSD.bmp’’ and ‘‘APSD_Iet.txt’’, and the second one ‘‘EffectiveRadius.bmp’’ and ‘‘EffectiveRadius.txt’’. Running time: 195 s for an input file consisting 12 measurement points (altitudes). © 2015 Elsevier B.V. All rights reserved.

1. Introduction Aerosols have a substantial and complex impact on atmospheric processes. For many aerosol types their contribution to various phenomena (for example the Earth surface radiation budget [1]) is poorly constrained. Understanding their optical properties and the distribution in the atmosphere is crucial for assessing their role in local and global scale. Active remote sensing with lidar technique is an important tool for investigation of the atmospheric aerosols [2,3]. The opportunities provided by multiwavelength lidars (MWL) are especially useful in this case since they enable determination of range resolved Aerosol Particle Size Distribution (APSD) which is a very important function for evaluation of aerosol optical properties.

evaluated using the approach of Bodhaine et al. [6] assuming the standard atmosphere profile as the basis. Aerosol scattering coefficients αAλ (z ) and βAλ (z ) described by relevant integrals can be evaluated by one of several possible aerosol models [7–13]. Mie theory can be used in the case of spherical particles [14]. In Eqs. (3)–(4) r denotes particle radius, and QλE and QλB — extinction and backscattering efficiencies, respectively. The inversion of Eqs. (3) and (4) provides n(z , r ) — distance dependent aerosol particle size distribution (APSD) which is the matter of investigation. For each distance APSD can be derived in a predefined form. Recently more often log-normal function: Cj (z )

1

 

· · exp − nj (z , r ) = √ 2π · log σj (z ) r

log r − log Rj (z ) 2 · log2 σj (z )

2  , (5)

or a linear combination of log-normal modes: 2. Approach Active remote sensing with lidar consists in sending of laser pulses at several wavelengths λ (λ = 1, 2, . . . , Λ) to atmosphere [2,4]. The light that is backscattered at distance z reaches the lidar receiver, where wavelength is separated and digitized. In consequence it provides the signals S1 (z ), S2 (z ), . . . , SΛ (z ): Sλ (z ) =

Aλ βλ (z )

(z − z0 )2



z



exp −2

 αλ (y) dy .

(1)

z0

Here z0 corresponds to lidar position, Aλ are the wavelength dependent apparatus constants while αλ (z ) and βλ (z ) denote the spatial distribution of total atmospheric extinction and backscattering coefficients respectively. The exponential factor reflects the use of Lambert–Beer law for the radiation extinction while the factor of 2 in the exponent is caused due to double run of the laser beam to the point z and back. Range corrected form of the signals (1) is more useful for the analysis:



Lλ (z ) = Aλ βλ (z ) exp −2



z



αλ (y) dy .

(2)

z0

Solution of these lidar equations provides both αλ and βλ coefficients. Each coefficient can be expressed as a sum of constituents [5]:

αλ (z ) = αRλ (z ) +  αAλ (z ) ∞

= αRλ (z ) +

QλE (r ) π r 2 n (z , r ) dr ,

(3)

0

βλ (z ) = βRλ (z ) + β  Aλ (z ) ∞

= βRλ (z ) +

QλB (r ) π r 2 n (z , r ) dr .

(4)

0

The first constituent of (3) or (4) corresponds to Rayleigh scattering. For standard atmosphere αRλ (z ) and βRλ (z ) can be

n( z , r ) =

K 

nj (z , r ),

(6)

j =1

is applied. Here Rj denotes the modal radius, Cj — the concentration of aerosol in a given mode, and σj — the mode width. Junge function with two free parameters is also used for this purpose. Then the inversion is reduced to the determination of optimal function parameters. Various inversion techniques have been proposed. The review of early approaches to this problem was presented by Zuev and Naats [15]. Herman et al. [16] fitted Junge distribution to the signals, while Rajeev and Parameswaran [17] have shown the iterative method without any assumed shape. The approaches to multifrequency lidar signals inversion were also performed by Piskozub [18] and Zieliński [19]. Heintzenberg et al. [20] proposed randomized minimization to derive an assumed histogram distribution. Similarly Müller and Quenzel [21] used this technique to test APSD determination assuming Junge and lognormal functions. Müller et al. applied Tikhonov regularization for the inversion [22– 25]. Veselovskii et al. [26] utilized the eigenvalue analysis for this purpose. A detailed review of different approaches to APSD retrieval was done by Böckmann [27] whose hybrid method was later applied to the experimental data analysis and the refractive index and single-scattering albedo retrieving [28]. Another solution was proposed by Shifrin and Zolotov [29]. They assumed APSD function in a form of combination of several log-normal functions. Using these distributions αλ and βλ were calculated and their values were compared to those measured by lidar technique. Then the mean ordinates over the solutions were calculated then the APSD closest to the mean ordinates was taken as the most probable aerosol distribution [30]. Ligon et al. [31,32] used the Monte Carlo method in order to shorten the calculation time. Kaufman et al. [33] expanded the method to invert Moderate Resolution Imaging Spectroradiometer (MODIS) and Cloud-Aerosol Lidar and Infrared Pathfinder Satellite

S. Sitarek et al. / Computer Physics Communications 199 (2016) 53–60

(CALIPSO) data simultaneously. Following aerosol models used in the MODIS data inversion they assumed APSD consisting of two log-normal modes, which were selected among four fine modes and five coarse modes. First twenty possible combinations of these modes are calculated on the basis of lidar signals only. Then with those profiles twenty simulated MODIS data are prepared. APSD providing the data that fits the best the real MODIS measurements is chosen as the result. In the approaches mentioned above the solution of lidar equations (2) requires additional information since both αλ and βλ are unknown. Commonly, one uses the lidar ratio: qλ =

βλ . αλk

3. APSD retrieval The lidar returns which are provided by the apparatus are quantized in space with the interval ∆z resulting from the digitization rate and smoothing procedure. At a distance zl for each wavelength the signals can by expressed by:

 Lλ (zl ) = Aλ βλ (zl ) exp −∆z

l 

[αλ (zi−1 ) + αλ (zi )] .

(8)

The integral from Eq. (2) was approximated by use of the trapezoidal approach. For further analysis the ratio of the signals at distances zl and zl+1 = zl + ∆z is taken: Lλ (zl )

Lλ (zl )

∞

QλB (r )π r n (zl+1 , r )dr

= 0 ∞ 0

× exp −∆z

QλB (r )π r n (zl , r )dr ∞



QλB (r )π r n (zl , r )dr

∞

 + 0

QλE (r )π r n (zl+1 , r )dr



.

(10)

As it was mentioned above the efficiencies QλE and QλB can be calculated using Mie theory [14]. In our approach n(z , r ) in a predefined form of two mode lognormal function (5) with several free parameters {C1 , R1 , σ1 , C2 , R2 , σ2 } is used. The solution of Eq. (10) is performed with minimization technique. For this purpose we construct a cost function on the basis of Eq. (10):

χ (zl ) = 2

 Λ  Lλ (zl+1 ) λ=1

Lλ (zl )

−

βλ (zl+1 ) βλ (zl ) 2

× exp {−∆z [αλ (zl ) + αλ (zl+1 )]}

.

(11)

The APSD with the set of the parameters that corresponds to minimal value of the cost function is accepted as the solution. Such nmin (r , z ) fits the best to the experimental signals. Deriving APSD from the extinction and backscattering coefficients is a classic example of ill-posed mathematical problem [47]. The number of the APSD parameters is often greater than the number of equations that are inverted. Moreover since APSD is related with the backscatter and extinction parameters by means of the integrals shown in (3) and (4), their inversion might not lead to unambiguous solutions. Therefore the approach requires additional constraints such as smoothness and positivity of derived function. The constraints might also consist in consideration of the parameters values from particular intervals. Limits of the intervals should follow from physical premises of analyzed aerosol. When the smoothness and the positivity of the derived functions are not satisfied, the change of the intervals should be considered. 4. Software In order to retrieve Aerosol Particle Size Distribution from multiwavelength lidar signals software described below was prepared. The programming language was Delphi 2010. Schematic diagram of the algorithm was shown in Fig. 1. 4.1. After start



i=2

=

Lλ (zl+1 )

0

The exponent fulfils the relation: 0.67 ≤ k ≤ 1.0 [34–36]. The value of the lidar ratio is often guessed or assumed [37,38]. The doubt of such approach consists in the fact that one usually assumes the independence of the lidar ratio on the distance from the lidar (or the height in the atmosphere). Such problems are widely discussed in Kovalev and Eichinger’s book [4]. Moreover it is not clear whether the relation (7) is unique one. Analysis of large dataset of EARLINET lidar returns shows about 40 % variability of the ratio for aerosols in boundary layer [37], therefore the use of such assumption can lead to significant errors [39]. Additional doubts are related to precision of the scattering and/or extinction coefficients (αλ and βλ ), which are necessary for APSD inversion. These coefficients can be found from the lidar signal with the inversion technique by Klett [36] and Fernald [40]. However these methods require independent information on aerosol parameters at a certain reference point. When the vertical profiling is performed the reference point is typically selected at high altitudes, where the aerosol concentration is often negligible and the molecular lidar ratio βRλ /αRλ = 3/8π can be used [6,41,42]. Then lidar equations are solved backward, starting from the reference point. Another approach may be applied in the case of clear sky. An independent measurement of aerosol optical thickness using sun-photometers allows to deduce the lidar ratio [6,43,44], but under assumption that it is constant with height. Such solutions are not applicable in many experimental situations, e.g. in horizontal measurements. Some of these problems may be overcome by simultaneous use of aerosol and Raman lidar [2,45]. In our recent paper a different approach to the APSD retrieving was presented [46]. It consists in direct substitution of αλ and βλ coefficients in form ((3), (4)) to the Eqs. Eq. (2). As a result only one unknown function remains: the APSD. The method does not require knowledge of lidar ratio and the experimental estimates of αλ and βλ values are not necessary. The aerosol is assumed to consist of droplets of known refractive index and of a predefined distribution form. Such approach is also applied in this software.

Lλ (zl+1 )

This allows omitting the apparatus constants Aλ , which are usually unknown. The left-hand side of formula (9) describes the ratio of the experimental signals while the extinction and backscattering coefficients in form (3) and (4) are substituted to the right side of the equations. Since both coefficients depend on aerosol distribution in Eq. (9) only one unknown function (i.e. APSD) remains. Therefore:



(7)

55

βλ (zl+1 ) exp {−∆z [αλ (zl ) + αλ (zl+1 )]} . βλ (zl )

(9)

After a start of the software the user can see a main window divided vertically into two sections as shown in Fig. 2. On the left side groups of panels are present. The panels called Wavelength, Particle Radius, Refractive Index, Signals, Calculations, Messages correspond to preparation of the input data as well as to maintaining of the calculation. The messages about the calculation progress and possible errors are located there. The user has to use these panels in order of their appearance.

56

S. Sitarek et al. / Computer Physics Communications 199 (2016) 53–60

practical form:

αAλ (z ) =



Rmax

QλE (r ) π r 2 n (z , r ) dr ,

(12)

QλB (r ) π r 2 n (z , r ) dr .

(13)

Rmin

βAλ (z ) =



Rmax Rmin

Fig. 1. Diagram of the algorithm.

On the right side the following display charts are located: Signals, Refractive Index, Mie Coefficient, APSD, Effective Radius, 3D Plot. Each diagram occurs in a separate table. The diagrams are successively filled with data under the calculations. 4.2. Input data The user has to input several initial parameters enabling the calculations. The set of wavelengths of analyzed signals is necessary (panel Wavelength). Default values are 1064, 532 and 355 nm. The user can add and remove other wavelengths given in nanometers. These values are automatically sorted in descending order. Next panels, called Particle Radius, Refractive Index are also used to provide the relevant initial parameters. The information about range of radii of aerosol particles {Rmin , Rmax } as well as their step (dr ) should be introduced. The default values of radii span from Rmin = 1 nm to Rmax = 10000 nm with a step of dr = 1 nm. In practice the integration over a limited radii range takes place in expressions for αAλ (z ) and βAλ (z ) in spite of the integration over the radii within the infinite range. Due to that Eqs. (3), (4) take the

The inversion of MWL signals requires the knowledge of refraction index of aerosol particles. Default value of this parameter was not defined. The user can pre-select aerosol type from drop-down menu or to load a text file containing wavelength dependent refractive index. The software will automatically fit the index of refraction for previously selected wavelengths. The software allows using an absorbing aerosol with imaginary part of the refractive index. For this reason the data in the input file is arranged in three columns: wavelength [nm], real part of ref. ind, imaginary part of ref. ind. Complex refractive index can be generated for two component aerosol mixture when its percentage composition is given. Well known rule that the value of refractive index of the mixture is equal to weighted average of its components is used in this algorithm [48]. The aerosol of the three or more components is not supported by the software. The refractive index of such mixture can be loaded from the text file as well. 4.3. Preparation of measured signals Loading of the signal file is done using the panel Signals shown as P3 in Fig. 2. The signals must be given in a single text file with proper number of columns. The first column has to denote the altitude in meters. The altitude z0 = 0 m should be equivalent to the beginning of lidar returns. Others columns have to be the signals from lidar photodetector (i.e. in the form of S (z ) records (1)) with the order in accordance to the descending values of the wavelengths. As far as in Eqs. (9)–(11) the ratio of the experimental returns is used, then the signals might be expressed in arbitrary units. For the calculations the signals have to be used in the rangecorrected form ((2), (8)). Such returns can be just loaded from the text file in its final form prepared by the user or can be loaded as a raw data for a processing using our software. Determination of APSD from multiwavelength lidar returns requires good preparation of the signals. The raw signals usually cannot be used for future elaboration because they are affected by various noises and interferences which provide the artifacts in final data.

Fig. 2. Main window of the APSD Software. The numbers represent the main panels: P1 — Wavelength, P2 — Particle Radius and Refractive Index, P3 — Signals, P4 — Calculations, P5 — Messages; P6 — represents the set of Display Charts.

S. Sitarek et al. / Computer Physics Communications 199 (2016) 53–60

57

Our software performs the following operations: change the sign of the signals, offset removal, smoothing and generation of range-corrected signals. 4.3.1. Change signals sign Raw optical signals (for example from photomultipliers) have negative values. For the user convenience the software allows to change their sign for all registered optical channels. 4.3.2. Offset removing The software allows removing the systematic errors of the lidar signal, usually called as an offset. The offset might occur due to interferences by the discharges in laser circuits, etc. Such systematic errors follow from the electric interferences that are synchronous with laser pulses. They might change with time. Other type of the offset occurs due to interference of the photonic signal by light following from the sources different than the laser of the transmitter (sun, moon, etc.) as well as due to dark current of the photodetector. As far as the processed signals are averaged over large number of laser shots then such kind of the offset usually does not change with time within the lidar return. The user is responsible for the proper recognition of a type of the error and the selection of relevant procedure of signal improvement. The relative part of the software is initiated by Offset Removing command. It provides several options. The best approach for removing of time dependent offset is the use of zero-file method. It consists in registration of lidar signals for blinded output of the transmitter. In such case only the electric interferences as well as the laser light scattered inside the lidar system, and the light following from the external sources (sun moon, lamps) are registered. Under the processing, the signal from the zero-file is initially smoothed. Adjacent averaging or Gaussian smoothing might be applied (see Section 4.3.3). Both procedures might be used with a window of constant width or that one of noise dependent width. These methods are described in Sections 4.3.3 and 4.3.4. Then the resulting file is simply subtracted from the lidar signal using pixel to pixel operation. After such procedure a main offset can be efficiently removed from the lidar signal. For time independent offset (Oλ ) resulting mainly from dark currents of the photodetectors and the photonic background the elimination software with a set of simple mathematic approaches was proposed (Constant Offset). First one, so called parabolic procedure was presented in Fig. 3. It is based on conversion of off experimental lidar signal Sλ (z )affected by the offset [49]: off

Sλ (z ) = Sλ (z ) + Oλ ,

(14)

(Fig. 3(a)) to the range corrected form: off

off

Lλ (z ) = (z − z0 )2 · Sλ (z ) = Lλ (z ) + Oλ · (z − z0 )2 .

(15)

(Fig. 3(b)) In many cases for large distances the value of Lλ (z ) vanishes with z − z0 increase while the offset constituent rises quadratically. Simple fit of quadratic function to the signal at large ¯ λ value distances (z > zmin ) provides opportunity of mean O off determination, which then is subtracted from Sλ signal. In order off

¯ λ determination the function Lλ (−z ), to increase the precision of O which is reflected in respect to OY axis, is artificially generated off as well (Fig. 3(c)). While at z = z0 the value Lλ = 0 the fit off

off

¯ λ (z − z0 )2 function is performed to both Lλ (z ) and Lλ (−z ) of O signals. The user selects the fit range for |z | > |zmin | using Scroll Bar in the software. In order to avoid the affecting of the procedure by residual value of Lλ (z ) signal the value of zmin should be distant enough in respect to z0 . Another possibility of constant offset determination is the use of the pre-trigger signal. In many cases the digitizers of lidar

Fig. 3. Offset removing and signal smoothing. (a) example of lidar return (1000 points), (b) range corrected signal, (c) preparation of the signal and fit of offset parabola, (d) range corrected signal with offset removed, (e) range corrected signal after adjacent averaging with the window of 20 points, (f) range corrected signal after adjacent averaging with noise dependent window ranging from 10 to 100 points.

receivers provide the signal that contains the points corresponding to the time moment before the laser shot (z < z0 — so called ‘‘negative’’ time moments). Averaging of these values provides ¯ λ determination and its subtraction from the the opportunity of O lidar signal. In our software the range of the altitudes for these calculations is also selected by the Scroll Bars. 4.3.3. Signal-to-noise improvement by adjacent averaging As it was mentioned above the raw single-shot lidar returns are usually affected by uncertainty (∆S ) of about several tens of percents. In order to improve signal-to-noise ratio the signals should be averaged over n number of √ laser shots. That leads to diminution of signal uncertainty over n value [50]. Typically this operation is implemented by the hardware under the registration of the signals. The assumption has been done that the loaded signals were initially prepared by the user this way. So the software does not contain any kind of such operation. Further signal-to-noise improvement can be achieved due to adjacent averaging over selected number of signal points (p). Due to both denoising procedures the uncertainty of the signals which √ are taken to calculations can be reduced by n · p factor [50]. Our software allows smoothing with rectangular or Gaussian kernel of selectable width p (Fig. 3(e)). When the digitizers of high frequency (f ) are used the signals might consist of hundreds or even thousands of points corresponding to the spatial resolution ∆z = c /2f which reaches the value of single meters sometimes. Such huge number of points is useless in many cases. Moreover when the averaging over the p — points is performed then the effective spatial signal resolution

58

S. Sitarek et al. / Computer Physics Communications 199 (2016) 53–60

is p∆z. The software creates new bins of the signals by putting the average values just in the middle of the adjacent intervals. 4.3.4. Adjacent averaging with noise dependent window As it is seen in Fig. 3(d) the noise to signal ratio of range corrected lidar returns usually rises quickly with the distance increase. The smoothing by adjacent averaging with constant window is not sufficient there, especially when the signal spans over large range of the distances (Fig. 3(e)). For such cases the option of adjacent averaging over noise dependent window is supported by our software. In order to determine optimal local value of the window width the distance range is initially divided into intervals of equal length. For each interval the mean signal value (M ) and the standard deviation (D) are calculated. Then D to M ratio (D/M ) as a function of distance z is plot and eD/M (z ) function is fit to such dependence. The fit is performed with second order polynomial because mainly the quadratic rise of L(z ) error is observed with the distance increase (see Fig. 3(b)). The function eD/M (z ) is interpreted as measure of fluctuations and noises. Next, taking into account the root dependence between averaging window width and D/M ratio improvement [50], the window width pD/M is determined according to the formula: pD/M (z ) = a eD/M (z )



2

+ b,

(16)

where a and b are the constants selected by the user. The value of pD/M (z ) varies usually from 10 points close to the signal beginning to 100 points at the end of the return range. Such procedure leads to the efficient increase and stabilization of the signal-to-noise ratio along the distance (Fig. 3(f)). The price is a gradually decreasing spatial signal resolution which is equal to pD/M (z )∆z. Therefore the using of such procedure is not suggested at small signal ranges where the smoothing with constant window is sufficient. Combination of averaging with noise dependent window and Gaussian averaging leads to using Gaussian smoothing with noise dependent widths. Such option is also supported by our software. It is worth to point out that as far as the adjacent averaging with constant window can be used for any signal, but the smoothing with noise dependent window requires good offset elimination before the using of the procedure.

Preparation of matrices of synthetic signals coefficients is done by means of the panels Particle Radius and Refractive Index. It takes place automatically when the buttons Matrix Gener in the panel Calculations are pressed consecutively. The procedures do not need any service by the user. 4.5. APSD retrieval The minimization of Eq. (11), that consists in search of optimal

op op op {Rop 1 , C1 , R2 , C2 } values in the look-up tables for all possible combinations of {R1 , C1 , R2 , C2 }, values, requires the calculation

of the cost function for tens of millions of cases. Since such a computation is very time consuming the Monte Carlo method is applied. The look-up tables are randomly drawn by multiple times. About 2 % of all the cases (i.e. about 300 thousands) is probed in this way. The value of χR21 C1 R2 C2 λ for each probe is determined. The optimal distribution nmin (r , z ) is selected as the case with minimum of χ 2 . However the bimodal lognormal function found by the minimization procedure is not unique. Quite satisfactory results (i.e. the comparable values of χR21 C1 R2 C2 λ parameter) can be achieved for several different {R1 , C1 , R2 , C2 }, sets. Therefore a group of 300 of results with a smallest χ 2 is chosen. Their χR21 C1 R2 C2 λ (i) values (i = 1 ÷ 300) do not differ more than about several tens of percent. Bimodal lognormal function nmin (r , z , i) is generated for each of them. As a final result of APSD retrieval (nopt (r , z )) the weighted average of these functions is accepted. The values of χR−12C1 R2 C2 λ (i) are taken as weights for this averaging. Reasonability of such approach was confirmed in our experiment [46]. During this campaign the APSD retrieved from lidar signals was compared with size distribution measured by particle counter. Finally in the software the APSD function is used for calculation of effective radius of aerosol particles using the formula:



reff (z ) = 

r 3 nopt (r , z )dr r 2 nopt (r , z )dr

.

(17)

This calculation is started when the button Fit I et is pressed. The reveal of results occurs when Analyze N And Surf is switched. 5. Exemplary analysis

4.4. Preparation of matrices of synthetic signals coefficients Minimization of the cost function (11) is time consuming. Therefore it requires special preparation of synthetic signals coefficients αλ and βλ . As it was mentioned above we attempt to approximate APSD within bimodal lognormal distribution ((5), (6)). Each of the lognormal function is characterized by Rj , Cj and σj coefficients (j = 1, 2). In such case the look-up tables of αRj Cj σj λ and βRj Cj σj λ coefficients are prepared for each mode separately. For each mode the modal radii Rj and the amplitudes Cj are arranged in geometric progression with a common ratio of 1,05. The minimization procedure of (11) selects the best solutions among these coefficients. For small particle mode the modal radii R1 are in a range 5,5–237,7 nm and the amplitudes C1 change from 70 to 9205 cm−3 . For large particle mode we use R2 from 230 to 1619,2 nm and C2 from 0,15 to 52, 3 cm−3 [5]. As far as σ — parameters concerns, it was stated, that within the radii range of 100–5000 nm many types of APSD can be well approximated by combination of lognormal functions with the widths equal of 2 [11]. Therefore this value of σ1 and σ2 was selected for both modes. Such approximation provides the opportunity to reduce the look-up table size by one order of magnitude or more. Due to that the job time was reduced by two orders of magnitude.

The example of analysis presented in this section was performed with data that were collected on July 30th of 2008 during lidar campaign in Warsaw. On this day a mass of marine polar air advected from the northeast in the peripheries of the Scandinavian high-pressure system. According to balloon atmospheric sounding from 12 UTC nearby WMO 12374 Legionowo station there was a well-mixed boundary layer extending up to 1660 m. The measurement was performed with lidar within this layer. Construction of the lidar was described in our previous papers [46,51]. Pulsed Nd:YAG laser generating light pulses at 1064, 532 and 355 nm (10 Hz repetition rate, 6 ns FWHM pulse duration time, the energies 100, 60 and 40 mJ respectively) was used in the transmitter. In the optical receiver, a Newtonian telescope with a 400 mm diameter mirror with a focal length of 1200 mm was used. The return light pulses collected by the telescope were spectrally separated by a polychromator and registered in three channels corresponding to the consecutive wavelengths. The signals from the photomultipliers installed in each channel were digitized by 12bit, 50 MHz A/D converters. During the data acquisition, the lidar profiles were averaged over 150 pulses. The results are shown in Fig. 4. Dot lines correspond to the APSD retrieved from consecutively collected lidar signals within the time period of 6 min. 24 signals were collected and for each

S. Sitarek et al. / Computer Physics Communications 199 (2016) 53–60

Fig. 4. APSD retrieved from lidar signals smoothed over 30 m. Dot lines — individual cases, continuous line — average result.

Fig. 5. Effective radius calculated with APSD presented in Fig. 4. Continuous line — the average result.

of them the APSD was retrieved. While the conversion frequency of our ADconverter corresponds to 3 m spatial resolution, in this case the smoothing of the signals over the window of 30 m was applied. Huge spread of the retrieved APSD functions is seen. The values of particle size distributions that were achieved for various time moments differ about one order of magnitude. That disagrees with physical premises about atmospheric properties and with the information about ‘‘the well-mixed boundary layer’’ that was mentioned above. That also leads to spread of effective radii of aerosol particles from about 300 nm to 1700 nm (Fig. 5). Noises affecting the signals lead to unreasonable results when APSD retrieving. Necessity of using of high quality lidar signals for APSD retrieval was discussed in our previous paper [46]. In the next approach the smoothing of the returns with the window ten times longer (300 m) was used in order to improve the signal quality. The results are shown in Fig. 6. Spread of the APSD is much lower. As it is seen in the inset of the figures it is smaller than 30 %. Mean effective radius (Fig. 7) reaches the value of 1050 nm for the altitude of 900 m and 1025 nm at 1200 m with the mean error of ±30 nm. These data fit well to other data that were considered for similar atmospheric conditions [52]. Effective range — 1, 5 micrometer, the exceeding range data follow from extrapolation by two mode lognormal function found with the technique described above. Nevertheless the experiment showed that even in this range the measurements with particle counter and the APSD retrieval from lidar signal were in satisfactory agreement. 6. Conclusion Software to retrieve spatial profiles of aerosol particle size distribution from multiwavelength lidar signals was presented. The

59

Fig. 6. APSD retrieved from lidar signals smoothed over 300 m. Dot lines — individual cases, continuous line — average result.

Fig. 7. Effective radius calculated with APSD’s smoothed over 300 m. Continuous line — the average result.

software enables full preparation of raw lidar returns. Offset elimination and signal-to-noise improvement are implemented. Various techniques applied for these operations provide opportunity to select the best approach. The purified experimental returns are compared with artificial signals generated using predefined distributions. The minimization technique allows finding APSD. The method was tested within the experiments with maritime aerosol as well with the aerosol under cumulus cloud [46,53]. Acknowledgments This work was supported by the Polish Ministry of Science and Higher Education with the statutory funds for research (2013/14, BST 166700/OP), the Faculty of Physics of the University of Warsaw (Institute of Experimental Physics) and Ministry of Economy with the statutory funds for research (NF/604), the Maksymilian Pluta Institute of Applied Optics. Appendix A. Supplementary data Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.cpc.2015.08.024. References [1] K.P. Shine, P.M.D.F. Forster, Glob. Planet. Change 20 (1999) 205–225. [2] R.M. Measures, Laser Remote Sensing. Fundamentals and Applications, Wiley & Sons, New York, 1984. [3] C. Weitkamp, Lidar, Range-Resolved Optical Remote Sensing of the Atmosphere, Springer, Berlin, 2005. [4] V.A. Kovalev, W.E. Eichinger, Elastic Lidar, A. John Willey & Sons Inc. Publication, Hoboken, New Jersey, 2004.

60

S. Sitarek et al. / Computer Physics Communications 199 (2016) 53–60

[5] J.H. Seinfeld, S.N. Pandis, Atmospheric Chemistry and Physics, John Wiley & Sons, New York, 1997. [6] B.A. Bodhaine, N.B. Wood, E.G. Dutton, J.R. Slusser, J. Atmos. Ocean. Technol. 16 (1999) 1854–1861. [7] M.I. Mishchenko, L.D. Travis, J. Quant. Spectrosc. Radiat. Transfer 60 (1998) 309–324. [8] M. Kahnert, J. Quant. Spectrosc. Radiat. Transfer 85 (2004) 231–249. [9] FORTRAN routine available at http://www.giss.nasa.gov/∼crmim/t_matrix. html. [10] M.I. Mishchenko, L.D. Travis, R. Kahn, R. West, J. Geophys. Res. 102 (1997) 16831–16847. [11] M. Hess, P. Koepke, I. Schult, Bull. Am. Meteorol. Soc. 79 (1998) 831–844. [12] O.V. Kalashnikova, I.M. Sokolik, J. Quant. Spectrosc. Radiat. Transfer 87 (2004) 137–166. [13] T. Nousiainen, M. Kahnert, B. Veihelmann, J. Quant. Spectrosc. Radiat. Transfer 101 (2006) 471–487. [14] C.F. Bohren, D.R. Huffman, Absorption and Scattering of Light by Small Particles, John Wiley & Sons, New York, 1999. [15] V.E. Zuev, I.E. Naats, Inverse Problems of Lidar Sensing of the Atmosphere, Springer Verlag, Berlin–Heidelberg–New York, 1983, p. 260. [16] B.M. Herman, S.R. Browning, J.A. Reagan, J. Atmos. Sci. 28 (1981) 763–771. [17] K. Rajeev, K. Parameswaran, Appl. Opt. 37 (1998) 4690–4700. [18] J. Piskozub, Oceanologia 28 (1990) 69–76. [19] T. Zieliński, Physical Properties of the Near-water Aerosol Layer in Coastal Areas, Institute of Oceanology, Polish Academy of Sciences, Sopot, 2006. [20] J. Heintzenberg, H. Müller, H. Quenzel, E. Thomalla, Appl. Opt. 20 (1981) 1308–1315. [21] H. Müller, H. Quenzel, Appl. Opt. 24 (1985) 648–654. [22] D. Müller, U. Wandinger, A. Ansmann, Appl. Opt. 38 (1999) 2346–2357. [23] D. Müller, U. Wandinger, A. Ansmann, Appl. Opt. 38 (1999) 2358–2368. [24] I. Veselovskii, A. Kolgotin, V. Griaznov, D. Müller, U. Wandinger, D.N. Whiteman, Appl. Opt. 41 (2002) 3685–3699. [25] I. Veselovskii, A. Kolgotin, V. Griaznov, D. Müller, K. Franke, D.N. Whiteman, Appl. Opt. 43 (2004) 1180–1195. [26] I. Veselovskii, A. Kolgotin, D. Müller, D.N. Whiteman, Appl. Opt. 44 (2005) 5292–5303. [27] C. Böckmann, Appl. Opt. 40 (2001) 1329–1342. [28] C. Böckmann, I. Mironova, D. Müller, L. Schneidenbach, R. Nessler, J. Opt. Soc. Amer. 22 (2005) 518–528. [29] S.K. Shifrin, I. Zolotov, Appl. Opt. 36 (1997) 6047–6056.

[30] [31] [32] [33] [34] [35] [36] [37]

[38]

[39] [40] [41] [42] [43]

[44]

[45] [46] [47] [48] [49] [50] [51] [52] [53]

S.K. Shifrin, I. Zolotov, J. Atmos. Ocean. Technol. 20 (2003) 1411–1420. D.A. Ligon, J.B. Gillespie, T.W. Chen, Appl. Opt. 35 (1996) 4297–4303. D.A. Ligon, J.B. Gillespie, P. Pellegrino, Appl. Opt. 39 (2000) 4402–4410. Y.J. Kaufman, D. Tanre, J.-F. Leon, J. Pelon, IEEE Trans. Geosci. Remote Sens. 41 (2003) 1743–1754. J.A. Curcio, G.L. Knestric, J. Opt. Soc. Amer. 48 (1958) 686–689. O.D. Barteneva, Bull. Acad. Sci. USSR 12 (1960) 1237–1244. J.D. Klett, Appl. Opt. 20 (1981) 211–220. E. Landulfo, A. Papayannis, P. Artaxo, A.D.A. Castanho, A.Z. De Freitas, R.F. Souza, N.D. Vieira Jr., M.P. Jorge, O.R. Sanchez-Coyllo, D.S. Moreira, Atmos. Chem. Phys. 3 (2003) 1523–1539. Y. Iwasaka, T. Shibata, T. Nagatani, G.Y. Shi, Y.S. Kim, A. Matsuki, D. Trochkine, D. Zhang, M. Yamada, M. Nagatani, H. Nakata, Z. Shen, G. Li, B. Chen, K. Kawahira, J. Geophys. Res. 108 (2003) ACE 20-1-8. S. Twomey, H.B. Howell, Appl. Opt. 4 (1965) 501–506. F.G. Fernald, Appl. Opt. 23 (1984) 652–653. J.E. Hansen, L.D. Travis, Space Sci. Rev. 16 (1974) 527–610. C. Frohlich, G.E. Shaw, Appl. Opt. 19 (1980) 1773–1775. E.J. Welton, Measurements of Aerosol Optical Properties over the Ocean Using Sunphotometry and Lidar (Ph.D. dissertation), University of Miami, Coral Gable, Florida, 1998. G. Karasiński, A.E. Kardaś, K. Markowicz, S.P. Malinowski, T. Stacewicz, K. Stelmaszczyk, S. Chudzyński, W. Skubiszak, M. Posyniak, A.K. Jagodnicka, C. Hochhertz, L. Woeste, Eur. Phys. J. 144 (2007) 129–140. B.R. Herman, B. Gross, F. Moshary, S. Ahmed, Appl. Opt. 47 (2008) 1617–1627. A.K. Jagodnicka, T. Stacewicz, G. Karasiński, M. Posyniak, S.P. Malinowski, Appl. Opt. 48 (2009) B8–B16. C.E. Groetsch, Inverse Problems in the Mathematical Sciences, Vieweg, 1993. Y. Liu, P.H. Daum, Aerosol Sci. 39 (2008) 974–986. M. Posyniak, T. Stacewicz, M. Miernecki, A.K. Jagodnicka, S.P. Malinowski, Opt. Appl. 40 (2010) 601–608. H.W. Ott, Noise Reduction Techniques in Electronic Systems, WileyInterscience, New York, 1988. A.K. Jagodnicka, T. Stacewicz, M. Posyniak, S.P. Malinowski, Proc. SPIE 7479 (2009) 747903. S. Arabas, H. Pawlowska, Geosci. Model Dev. 4 (2011) 15–31. T. Stacewicz, M. Posyniak, S.P. Malinowski, S. Sitarek, Atmos. Res. 142 (2014) 32–39.