Summer Research Fellowship Programme of India's Science Academies

To study cloud extinction using Lidar technique

Pranjal Dave

Msc. Mathematics, B.E. Electrical and Electronics Engineering, Birla Institute of Technology and Science, Pilani, Pilani Campus 333031

Dr. Y. Bhavani Kumar

Scientist/Engineer - SF, Project Leader, LIDAR project, National Atmospheric Research Laboratory, Department of Space, Tirupati 517502


LIDAR (Light Detection and Ranging) is a form of active remote sensing technique that uses the reflection of light pulses from objects for observations. It is similar in principle to RADAR (Radio Detection and Ranging) but in RADAR radio waves are used as source signal while in LIDAR light waves are used. LIDAR has a wide range of applications from Driver-less cars to atmospheric studies. The LIDAR system basically consists of a transmitter and a receiver. The transmitter sends light pulses of some fixed wavelengths into the atmosphere. The receiver consists of an apparatus to detect the returning pulses for acquiring data. In the present study, we use an atmospheric LIDAR to observe and measure the various properties of the atmosphere, particularly focusing on clouds. LIDAR cloud measurements form an important data-set for atmospheric studies. They can be used for detecting cloud base, for visibility studies and to study the boundary layer dynamics. As any light signal passes through the atmosphere, it gets attenuated (scattered) due to the various molecules and particles present in the atmosphere. The main aim of this study is to calculate the optical depth of cloud layers that advect in the atmosphere. For this purpose, one has to determine the range resolved light attenuation through the cloud layer. The loss of optical energy due to scattering in the atmosphere is known as extinction. The study of LIDAR cloud measurement data provides a way to estimate the attenuation coefficient. In this report, we have used the Klett (1981) inversion method to determine the extinction coefficient from lidar signals.

Keywords: Lidar equation, range-corrected signal, Klett's inversion method, slope method, molecular backscatter coefficient, optical depth


The abbreviations used in this report are explained below:

LIDARLight Detection and Ranging
RCSRange Corrected Signal
RTI Range- time Intensity 



Remote Sensing (RS) is a branch of modern science that acquires data from an object without making physical contact with the object. RS uses wide-band of Electromagnetic (EM) radiation that also covers light spectrum. RS is of two types: Active and Passive. In the active RS technique, the source of EM radiation is human-made. The sensor itself emits EM beams and then records response from the target, whereas in the case of passive RS method, the source of the EM radiation is natural like sunlight, moonlight, etc. In this type, the rays from the Sun hit the object of interest and then the sensor records the response from the object. LIDAR (Light Detection and Ranging) is an example of the active RS technique. It is similar to RADAR but it uses light waves instead of radio waves. In this technique, lidar sends pulses of light to the object to be detected and records time of flight (TOF) from the reflection. In the present study, we used a ground based lidar to study the range resolved loss of optical energy due to clouds. For this, we employed an elastic backscatter lidar (EBL) system that operates at the harmonic wavelengths of Nd:YAG laser (​Yellapragada, et al, 2016​; ​Vishnu, et al, 2017​). The EBL system is an indigenously developed ground based lidar instrument that operates at IR (1064 nm) and visible (532 nm) wavelengths even in daylight conditions. In this report, the EBL system signal returns that correspond to the Nd:YAG laser wavelengths 1064 and 532 nm were used for estimation of light extinction and backscattering coefficient. The block diagram of lidar setup ( ​Vishnu, et al, 2017​ ) used in the report is shown in Figure 1. The laser unit of EBL is a compact flash-lamp pumped system that sends light pulses correspond to the colors of IR and Green simultaneously at repetition rate of 20 per second. The backscattered light returns are collected using a 150 mm telescope. The collected light returns are converted into electrical signals using a PMT or APD depending upon wavelength of interest. Transient recorder (TR) unit is employed in the receiver system to record the range resolved TOF signals. A computer system uses Ethernet bus to obtain signals from TR unit.

The EBL system provides us a way to visualize clouds as they will be responsible for most of the reflections and hence a peak in intensity of the reflected signal indicates the height from ground at which a cloud is present. We can also find the attenuation of the light signals due to the atmosphere and also the optical depth of the atmosphere. This parameter plays a vital role in atmosphere science and aviation meteorology which allows us to find out the visibility, optical depth of the atmosphere etc

Screenshot (19).png
    Block diagram of EBL system used in the study ( Vishnu, et al, 2017 )

    This provides us with a way to visualize the clouds as they will be responsible for most of the reflections and hence a peak in intensity of the reflected signal indicates the height from ground at which a cloud is present. We can also find the attenuation of the light signals due to the atmosphere and also the optical depth of the atmosphere.These factors are of crucial importance as they are used to then find the visibility, energy loss due to the atmosphere etc.

    LIDAR equation

    LIDAR equation forms the basis of the equations used for the atmospheric studies. The various equations used in this study can be derived from the LIDAR equation. In basic terms, LIDAR equation expresses the relation between the received signal and the different variables related to the measuring instrument, atmospheric conditions and the properties of the beam sent. The equation can be expressed as:

    P(r)  =  Pocτ2Aβ(r)r2exp[20rσ(r)  dr]P(r)\;=\;P_o\frac{c\tau}2A\frac{\beta(r)}{r^2}exp\left[-2\int_0^r\sigma(r')\;dr'\right]

    where ,

    P(r) is the instantaneous received power at time t

    𝛔(r) is the attenuation coefficient

    𝜷(r) is the backscatter coefficient

    A is the effective receiver area

    Po is the transmitted power

    ͳ is the pulse duration

    r is the range

    As the required coefficients are present on the right hand side of the equation, we use inversion of the LIDAR equation to get their expressions. Some of the quantities like A, Po etc. are controllable by the experimentalist and are thus combined into a single constant.

    Statement of the Problem

    "To calculate various parameters of the atmosphere using LIDAR data "

    The following needs to be addressed:

    • Extracting the raw data from the text files.
    • Pre-processing the data to remove the noise and applying range correction.
    • Calculating molecular backscatter and total backscatter.
    • Using the Klett method to obtain the sigma values and thus obtain the LIDAR ratio.
    • Finding the optical density.

    Objectives of the Research

    Overall objective

    • To process data by removing background noise and applying range correction
    • To obtain range-time intensity profile
    • To calculate backscatter and extinction coefficient of the atmosphere
    • To find the optical density of the clouds


    The study concerns with the observation and study of the atmosphere and clouds over a tropical rural site Gadanki in Andhra Pradesh. The data was recorded using multi-wavelength EBL system. Two wavelengths were used in the study namely 532nm and 1064nm (​Vishnu, et al, 2017​). Only the atmospheric conditions concerning clouds were focused upon and factors like wind speed, humidity etc. were not used. The cloudless and partially cloudy atmospheres were considered and phenomenon like rains was ignored for the study.

    There are various methods available to obtain the extinction coefficient from the LIDAR equation. In the present study, we have focused on two main methods, one describes slope method and the other uses Klett (1981) approach. The other methods have their advantages as well as disadvantages but these two methods are used due to the availability of vast literature about these methods and were known to provide reliable results.


    The method followed in this study is taken from the research paper by James D. Klett titled ‘Stable analytical inversion solution for processing lidar returns’ (​James D. Klett, 1981​). In this paper, Klett discusses an analytical method to calculate attenuation and backscatter coefficients. The paper begins with a review of slope method which is used to calculate the attenuation 𝛔 directly from the plot of S(r). Here S(r) is defined as:

    S(r)  =    ln  [r2  P(r)]S(r)\;=\;\;\ln\;\lbrack r^2\;P(r)\rbrack

    In this method, the atmosphere is considered to be homogeneous so that d𝜷/dr = 0, giving us a direct relation between 𝛔 and S(r).

    σ=  12dSdr\sigma=-\;\frac12\frac{dS}{dr}

    So the slope of the least squares straight line fit to the curve S = S(r), is used to calculate 𝛔 from the above formula. The 𝛔  obtained from the above method is not accurate as in the case of clouds, the atmosphere cannot be considered homogeneous. To make the slope method more accurate, we also have to consider the variation of 𝜷 with the range (i.e. d𝜷/dr ≠ 0). 

    A relation between 𝜷 and 𝛔 is assumed. As mentioned by Klett in his paper, they can be related by a power relation

    𝜷 = constant 𝛔k

    Here k is a constant which is known to vary between 0.67-1. This reduces the LIDAR equation into a usable form as described in the methodology.

    A report titled "On the inversion of lidar equation" (​Evans, et al​,1984) was also referred. The report discusses in detail the derivation of the LIDAR equation from the basic laws of radiation. It also discusses the inversion of LIDAR equation is great depth and then goes to discuss the different ways to obtain the coefficients from the equations. A new method is also suggested by the author, but it was not considered for this study as it is not as well established as Klett's method. Finally, the author provides with the MATLAB code for the method described in the report. The code was referred while writing the code for this study.

    The paper Carnuth, et al, 1986 was also consulted as it is in line with the objectives of this study. The paper discusses using Klett's inversion method to find cloud extinction profiles which is similar to the process followed in this study to find the extinction coefficient. The figures presented in the paper were also consulted in this study.



    Inversion of LIDAR equation

    As mentioned before, the inversion of LIDAR equation is an essential step to obtain the atmospheric coefficients. After substituting S(r), as defined above, into the equation, we obtain

    S    So  =  lnββo2rorσ  drS\;-\;S_o\;=\;\ln\frac\beta{\beta_o}-2\int_{r_o}^r\sigma\;dr'

    here βo = β(ro) . Now we differentiate both sides to obtain the corresponding differential equation

    dSdr  =  1βdβdr  2σ\frac{dS}{dr}\;=\;\frac1\beta\frac{d\beta}{dr}-\;2\sigma

    Using the power relation (​​Sec. LIDAR equation​​), we eliminate β from the equation to get the equation only in terms of σ.

    dSdr=kσdσdr2σ\frac{dS}{dr}=\frac k\sigma\frac{d\sigma}{dr}-2\sigma

    This transforms the equation into a non-linear ordinary differential equation but the form of this equation is the well known Bernoulli or homogeneous Ricatti equation and solution for this has been known for hundreds of years. We obtain the following form of solution

    σ=exp(SSo/k)σo1    2krorexp(SSo/k)dr\sigma=\frac{exp(S-S_o/k)}{\sigma_o^{-1}\;-\;{\displaystyle\frac2k}\int_{r_o}^rexp(S-S_o/k)dr'}

    But this form of the equation does not represent the Klett's form. Here due to the fact that two relatively large terms are getting subtracted in the denominator, the solution becomes unstable. Klett suggested that instead of ro , rm (the maximum distance from the lidar system) should be taken as the reference point. This change makes the solution stable and is known as the Klett's inversion method. The final form used in this study has been taken from the paper ​​Uri Dayan, 1991​​. Here the constant k has been assumed to be 1.

    σ=exp(SSm)σm1  +  2rrmexp(SSm)dr\sigma=\frac{exp(S-S_m)}{\sigma_m^{-1}\;+\;{\displaystyle2}\int_r^{r_m}exp(S-S_m)dr'}

    here σm refers to the reference value taken for sigma. In the present study, we use slope method to get an initial approximation of sigma values.


    Data processing

    The LIDAR scan data was provided in the form of .txt files. The files contains data about the various parameters like range scale, number of data values etc. used for recording the data. It also contains the data recorded by the various channels of the LIDAR system. It was converted to excel format using Origin. The data was then exported to MATLAB for further processing. 

    First, the data were corrected for the background noise present in it. For this, the data values after a specific range were considered as noise and they were summed together. The sum was then divided by the number of bins that were present in the given range. Thus the noise per bin was obtained. This was subtracted from each data point as the background noise. The result is the noise corrected data. This data is then processed further to obtain the range-corrected data. The squares of range values were multiplied to corresponding data points to get the range-corrected data. This was used for further operations. This represents [r2P(r)], we then take ln of the data to obtain S(r). In ​​Fig 2​​, the first panel represents the raw data, the second panel represents the noise corrected data, the third panel represents the range corrected data and the fourth panel represents the natural logarithm of the range corrected data.

      Sequence of stages showing the LIDAR data processing

      Obtaining the molecular backscatter coefficient (𝜷)

      Backscatter coefficient in simple terms refers to the number of light waves that scatter by 180 degrees i.e. in the direction of the incoming light. This coefficient can be broken down into backscatter due to molecules and due to aerosols. The coefficient gives the backscatter per length unit and per solid angle (steradian).

      𝜷total = 𝜷aero + 𝜷mol

      The molecular backscatter coefficient varies with the height from the surface of the earth and can be calculated using temperature and pressure at the given heights.’betamol_standard’ function of the R Atmosphere library of R was used to calculate the molecular backscatter coefficient for the given range. The data was then imported into MATLAB for further processing. Figures​​ Fig 3a & 3b​​ depicting the variation of 𝜷 for 532nm and 1064nm respectively, with range are shown below.

      Altitude distribution for molecular backscatter cross section at 532 nm
        Altitude distribution for molecular backscatter cross section at 1064 nm
          Profiles of molecular backscatter coefficient at 532nm and 1064nm

          Algorithm to calculate extinction coefficient

          As in the case of backscatter coefficient, the extinction coefficient can be expressed as the sum of two components, namely the aerosol extinction coefficient and the molecular extinction coefficient.

          σtotal =σaero + σmol

          The formula described above ( ​​Sec. 3.1​​ ) gives the total extinction coefficient. Form this value, the molecular extinction coefficient has to subtracted to calculate the aerosol extinction coefficient.

          To calculate total extinction coefficient

          After obtaining S(r) from the data (for both 532 nm and 1064 nm), the equation for σ(r) was implemented in MATLAB. The following algorithm was used:

          • Calculate reference sigma (σm) by using slope method. The ' \ ' operator was to used find the slope of the best fit line on the plot of S(r) vs range. The values of S(r) were only taken upto a specific range as they represented the signal of interest. The slope was multiplied by -0.5 to obtain σm.
          • Set Sm and S in the equation. The value of S was set as the column values of the S(r) matrix. S m represents the value of the particular column at mth row. Here m is the reference range taken.
          • Evaluate exp(S-Sm) and hence to calculate the integral term. The expression exp(S-Sm) can be easily evaluated in MATLAB using the exponential function. For finding the integral term, the 'trapz' function of MATLAB was used. Input to it was the array of elements obtained after evaluating the previous expression. A loop on the row numbers was used to implement the definite integral.
          • Calculate σ. After finding all the terms, they were into the equation to calculate the extinction coefficient.

          To calculate molecular extinction coefficient

          The molecular extinction coefficient can be evaluated using the molecular backscatter coefficient for the particular wavelength. They have the following relationship (Thorin, et al, 2006 )

          σmolβmol = 8π3

          Thus the obtained molecular extinction coefficient can be subtracted from total extinction coefficient to obtain the aerosol extinction coefficient.

          Calculating Optical density

          Optical density of the atmosphere can be calculated using the values of the extinction coefficient. Two types of optical density can be calculated:

          1. Cloud optical density

          2. Aerosol optical density

          Firstly, using the plot of RCS vs range,the top and the bottom cloud boundaries are determined. Then, on the basis of which optical density is to be calculated, the extinction coefficient values are integrated within the specific ranges. If the range is taken from ground to the bottom cloud boundary, we obtain the aerosol optical density. If the range is taken from the bottom cloud boundary to the top cloud boundary, we obtain the cloud optical density

          CBCT σtol dr = σ1(30) + σ2(30)+ ......

          Here σ1, σ2 etc. refer to the values of the extinction coefficient at each point in the range CB to CT and 30 refers to the range scale as each data point is 30m apart from the other.


          Range-time Intensity (RTI) profile

          RTI plot refers to a 3-dimensional representation of the lidar data with range as the y axis, time as the x axis and the intensity being represented using various colors. It allows for a better understanding of the data, and also helps in visualizing the clouds and phenomenon such as rains etc. The following RTI plot was plotted using 1 hr of LIDAR data. The range increment in the RTI plot shown in ​​Fig 4​​ is 3.75 meters.

            RTI plot indicating cloud and rain periods over lidar sit

            Waterfall plot

            Waterfall plot shown in ​​Fig 5​​ is similar to the RTI plot, however, it gives additional information on discrimination between clouds and rain events. In addition to intensity being represented by colors, it is also represented on the z axis, so it is a 3-D plot. The waterfall plot below is of the same data as the RTI plot, but the time and range intervals are different. Here also y axis represents range and x axis represents time. The intensity is represented on the z axis.(​​Fig 5​​). ​​Fig 6​​ is a zoomed part of ​​Fig 5​​ shows clearly the rainfall event. The temporal and spatial variation of rain is depicted in Fig.6.

            Waterfall plot.jpg
              Waterfall plot
                A zoomed view of the above plot to show the peaks and troughs

                Sigma Values of different time scales

                Here sigma refers to the total extinction coefficient. The values obtained were plotted in different ways to ascertain their validity. The following four figures denote the sigma values calculated by taking time averages of the data. As the data was provided in intervals of 30 sec, to obtain 1 min, 2 min etc. data, the required number of files were averaged together. ​​Fig 7a​​, ​​Fig 7b​​, ​​Fig 7c​​ and ​​Fig 7d​​ represent the various time intervals (1 min, 2 min, 4 min and 8 min) used to generate the values of the extinction coefficient respectively. The extinction profiles were plotted over 500 receiver bins which represent a range corresponding to 15km. The lidar was operated with 30 bin range resolution at a time sampling of 1 sec.

                1 min values
                  2 min values
                    4 min values
                      8 min values
                        Various time averaged sigma values

                        Backscatter Coefficient

                        The backscatter coefficients were calculated for 532nm. Further research on these coefficients and the working of the LIDAR system can be found in ​Yellapragada, et al, 2016​. The following figures (​​​Fig 8a​​,​​Fig 8b​​ and ​​Fig 8c​​) compares the backscatter coefficients for 10 minutes intervals, by plotting total backscatter with molecular backscatter coefficient. They appear to be overlapping after a point. It can be explained by the fact that in lower altitudes the backscatter is mainly due to aerosol and hence,the total backscatter is more than molecular backscatter coefficient. But as the altitude increases, molecular backscattering starts dominating and hence both of the plots start overlapping.

                        10:00 - 10:10 am
                          10:10- 10:20 am
                            10:20 - 10:30 am
                              Total backscatter coefficient with molecular backscatter coefficient

                              Extinction Coefficient

                              The extinction was calculated using the method given above. The altitude profiles of extinction were calculated for 532 nm wavelength, using the lidar data. These profiles are shown in figures ​​Fig 9a​​,​ ​Fig 9b​​ and ​​Fig 9c​. As in the case of backscatter coefficient, the extinction coefficients are also plotted in intervals of 10 minutes (​​Fig 9a​​,​ ​Fig 9b​​ and ​​Fig 9c​​) with molecular extinction coefficient plot superimposed on the total extinction coefficient plot. Similarly, the plots come close after certain range due to the fact that molecular extinction coefficient becomes dominant on higher altitudes.

                              10:00 - 10:10 am
                                10:10 - 10:20 am
                                  10:20 - 10:30 am
                                    Total extinction coefficient with molecular extinction coefficient

                                    LIDAR ratio

                                    LIDAR ratio refers to the ratio of the extinction to backscatter coefficient. It represents the efficiency of the scatterer with the light interaction process which generates attenuation to the light passing through the scatterers. We have calculated the above ratio utilizing backscatter and extinction coefficients that are calculated form the lidar. The ratio should remain constant upto the cloud boundary layer as only aerosol is present till there. After that the ratio varies due to the clouds present. ​​Fig 10c​​ can be obtained by dividing the data of ​​Fig 9a​​ by the data of ​​Fig 8a​​. Similar process can be carried out to obtain ​​Fig 10c​​ and ​​Fig 10c​​.

                                    10:00 - 10:10 am
                                      10:10 - 10:20 am
                                        10:20 - 10:30 am
                                          LIDAR ratio

                                          Optical Density

                                          The following figure (​​Fig 11a​​) denotes the variation of optical density of clouds throughout the day for 1st April, 2019 over lidar site. At around 3:30 pm we can see a peak in optical density. It can be explained by observing the RTI plot(​​Fig 11b​​) of the same. Around the same time we can observe a cloud coming at a height of around 4 km from the ground, thus increasing the optical density.

                                          Variation of optical density with time
                                            RTI plot
                                              Figures showing increase in optical density due to clouds


                                              In this study, we verified that Klett's inversion method is reliable for calculating the backscatter and extinction coefficients. The LIDAR ratio was found to be around 16 for the aerosol and it varied for the clouds.

                                              It is also observed that the changes in optical density corresponed to the changes in the cloud layer. It affirms that fact that the extinction coefficient values correspond to the real world phenomenon and hence are valid.

                                              It can also be established from the results that the aerosol layer extends upto 4 km from the ground level. It is shown by the fact that the LIDAR ratio remains constant till 4 km and also by the fact that the plots of extinction coefficient and backscatter coefficient start overlapping after 4 km.


                                              I am thankful to Dr.Y Bhavani Kumar for providing me with this opportunity to work under his guidance. I am also very grateful to him for his constant support and his constructive feedback throughout the project.

                                              I am also very thankful to Mrs.Sandhya for her constant aid in this project by providing the data and helping in any software related problem.

                                              I am also grateful to the Indian Academy of Sciences for providing me with this wonderful opportunity to work with scientists and learn from them. I am also very thankful to them for providing me with scholarship for this internship.

                                              I am also thankful to my university, Birla Institute of Technology and Science, Pilani for allowing me and assisting me in pursuing the internship at NARL, Tirupati.


                                              • Yellapragada, Bhavani Kumar (2016). Lidar research activities and observations at NARL site, Gadanki, India.

                                              • Vishnu, R. and Kumar, Y. Bhavani and Sinha, P. R. and Rao, T. Narayana and James Jebaseelan Samuel, E. and Kumar, P. (2017). Comparison of mixing layer heights determined using LiDAR, radiosonde, and numerical weather prediction model at a rural site in southern India. 38,

                                              • James D. Klett, 1981, Stable analytical inversion solution for processing lidar returns, Applied Optics, vol. 20, no. 2, pp. 211

                                              • Evans, B. T. N., 1984: On the inversion of the lidar equation. DREV Rep. 4343/84, Defense Research Establishment Valcartier, 1–50. [Available from DREV, P.O. Box 8800, Courcelette, PQ G0A 1R0, Canada.]. 

                                              • Carnuth, Walter and Reiter, Reinhold (1986). Cloud extinction profile measurements by lidar using Klett’s inversion method. 25,

                                              • Ram Hashmonay, Ariel Cohen, Uri Dayan, 1991, Lidar Observation of the Atmospheric Boundary Layer in Jerusalem, Journal of Applied Meteorology, vol. 30, no. 8, pp. 1228-1236

                                              • Thorin, E. (2006). Design of Algorithms to Extract Atmospheric Aerosol Extinction from Raman Lidar Data (Dissertation). Institutionen för fysik, kemi och biologi. Retrieved from http://urn.kb.se/resolve?urn=urn:nbn:se:liu:diva-7912

                                              Written, reviewed, revised, proofed and published with