## Abstract

The conventional method of calculating atmospheric temperature profiles using Rayleigh-scattering lidar measurements has limitations that necessitate abandoning temperatures retrieved at the greatest heights, due to the assumption of a pressure value required to initialize the integration at the highest altitude. An inversion approach is used to develop an alternative way of retrieving nightly atmospheric temperature profiles from the lidar measurements. Measurements obtained by the Purple Crow lidar facility located near The University of Western Ontario are used to develop and test this new technique. Our results show temperatures can be reliably retrieved at all heights where measurements with adequate signal-to-noise ratio exist. A Monte Carlo technique was developed to provide accurate estimates of both the systematic and random uncertainties for the retrieved nightly average temperature profile. An advantage of this new method is the ability to seed the temperature integration from the lowest rather than the greatest height, where the variability of the pressure is smaller than in the mesosphere or lower thermosphere and may in practice be routinely measured by a radiosonde, rather than requiring a rocket or satellite-borne measurement. Thus, this new technique extends the altitude range of existing Rayleigh-scatter lidars 10–15 km, producing the equivalent of four times the power-aperture product.

© 2012 Optical Society of America

## 1. Introduction

The most frequently employed temperature analysis for Rayleigh-scatter lidar measurements is that of Haunchecorne and Chanin [1] (henceforth referred to as the CH method). In the CH method, the atmosphere is considered an ideal gas in hydrostatic equilibrium. Since the lidar equation relates the returned photocounts to the density of the atmosphere, the density can be related to temperature. Since hydrostatic equilibrium is a statement of the rate of change of atmospheric pressure with height and not a relation between absolute pressure values and altitude, the CH method requires an assumption of a reference pressure (or temperature) commonly referred to as the seed condition and chosen to be at the greatest height where sufficient signal-to-noise ratio (SNR) measurements are available. This choice is a major limitation in the conventional technique of temperature retrieval, as an incorrect choice of the seed pressure at the top can introduce significant errors (e.g., much larger than the random uncertainties due to photon counting statistics) near the top of the altitude range in the retrieved temperature profile. Another assumption of this method is an isothermal atmosphere between retrieval points to further simplify the analysis procedure and produce an analytical solution for temperatures. With improvements in computational capabilities it has now become possible to solve such problems numerically, without the requirement of these simplifying assumptions. This paper describes a solution to this limitation using a numerical technique of mathematical inversion with grid search optimization. There are two key advantages of this new method. First, a bottom seed condition can be used, allowing a more accurate guess of the seed pressure. Second, the choice of seed pressure in the inversion approach extends the lidar’s height range by allowing the uncertainties due to choice of seed pressure to be much less than the statistical uncertainties of the measurement at the greatest heights.

## 2. Background

The CH method of atmospheric temperature retrieval from Rayleigh-scatter lidar measurements is based on two assumptions, the first that the atmosphere behaves as an ideal gas and the second that the atmosphere can be considered static and in hydrostatic equilibrium at the spatial-temporal resolution at which the measurements are obtained.

It can be shown using these assumptions that the temperature, ${T}_{i}$ at height ${z}_{i}$ is

The second term in Eq. (1) is independent of the seed pressure value. Thus, the systematic uncertainty from the estimation of ${P}_{0}$ affects only the first term of the right-hand side. The seed pressure value can be represented as a sum of the true signal ${P}_{\mathrm{true}}$ and the error $\mathrm{\Delta}P$ as ${P}_{0}={P}_{\text{true}}+\mathrm{\Delta}P$. Then Eq. (1) becomes

whereIt is evident that the uncertainty in calculated temperatures is directly proportional to $\mathrm{\Delta}P$ and because of its inverse relationship with atmospheric density, $\mathrm{\Delta}T$ increases as the density decreases with altitude. Also, because the seed pressure is chosen at the top of the atmosphere, a large uncertainty is expected due to the lack of routine pressure measurements at these heights compounded by the large geophysical variability that exists in the upper atmosphere due to atmospheric waves.

Equation (2) also shows why integration from the bottom is not possible in the CH method. Since the pressure decreases as one moves up in altitude, the cumulative difference ${P}_{0}-\sum _{j=1}^{i}\rho ({z}_{j})g({z}_{j})\mathrm{\Delta}z$, instead of a summation, may become negative because a difference between two large and approximately equal numbers is performed. This difference can make the pressure estimates negative and hence the retrieved temperatures start to diverge rapidly from their real value. This divergence will happen even if the seed pressure value is only slightly erroneous, making the retrieved temperatures highly sensitive to the choice of ${P}_{0,\text{bottom}}$. The new method presented next will be based on the same physical principles as the CH method but offer solutions to its limitations.

## 3. Inversion Method

The method presented here is an inversion approach. This method is applied to retrieve atmospheric temperatures from the coupled form of the system equations. The theory of mathematical inversion is based upon the relationship between the state of the system or system parameters and the observed variables of the system. This relationship can be expressed by the forward model of the system given by

where, ${\mathbf{d}}_{N\times 1}$ is the vector of observed data and ${\mathbf{m}}_{M\times 1}$ is the vector of unknown model parameters. $\mathbf{F}$ represents the forward model (letters in bold face represent matrices or vectors).The inversion approach, like other optimization schemes, is an iterative approach to approximate the true value of parameters by starting from an initial value for them and updating them with an improved value, step by step by following a predefined optimization criterion. If ${\mathbf{m}}_{0}$ is an initial guessed value of the parameter vector $\mathbf{m}$ and ${\mathbf{m}}_{0}+\delta \mathbf{m}$ the true value that is sought, $\mathbf{F}({\mathbf{m}}_{0})$, becomes an estimate of $\mathbf{d}$ and $\mathbf{F}({\mathbf{m}}_{0}+\delta \mathbf{m})$, the true value ${\mathbf{d}}_{\text{true}}$. The crux of the inversion method is to minimize the difference between the estimated and true value of $\mathbf{m}$ by minimizing the difference between the estimated and observed value of $\mathbf{d}$. The condition on convergence of the guessed solution to the true one is implemented by minimizing the objective function (${\chi}^{2}$), i.e., a weighted least-squares minimization (Bevington and Robinson Chap. 4 [2]). Following the method of maximum likelihood, the ${\chi}^{2}$ difference is defined assuming that the data are sampled from a Gaussian distribution and is expressed as follows:

The objective function is defined for a data vector $\mathbf{X}$ of size $N$ and mean vector $\overline{\mathbf{X}}$, where ${x}_{j}$ is an element of the data vector, $\overline{{x}_{j}}$ is the expected value of that variable, and ${\sigma}_{j}^{2}$ is the data variance for the $j$th variable.

To apply this method for Rayleigh-scatter temperatures requires a forward model (discussed in the next section). The model is initialized (in this case with a temperature profile) and the model parameters varied until the difference between the model and the measurements is a minimum.

A reasonable choice for the convergence condition is that the change in ${\chi}^{2}$ per degree of freedom (or for each parameter) per minimization step is less than or equal to 0.1%. If the parameters are correlated to one other, the true definition of the ${\chi}^{2}$ function for correlated data can be used, which is

Equation (8) is a general expression for ${\chi}^{2}$, the square root of which is called the Mahalanobis distance [3]. This expression is used as the objective function in the current study as the system parameters (temperatures) have correlation among themselves due to geophysical variation. The next section will develop how the grid search method of inversion can be applied to the temperature retrieval from Rayleigh-scatter lidar measurements.

## 4. Implementation of the Inversion Approach

In the CH method, the decoupled system equations are used to obtain an approximate analytical form. The more sophisticated forward model in the inversion approach allows the assumption of constant temperature between layers to be removed. This removal is achieved by integrating the equation by using a trapezoidal-rule integration instead of a simple rectangular integration in the following manner:

From Eq. (13), $\rho ({z}_{i})$ can be written in terms of the lidar return signal at the altitude ${z}_{i}$, which modifies Eq. (12) to

Equation (14) has three unknowns: $T(z)$, ${P}_{0},$ and $C$. As explained below, the constant $C$ is obtained in the usual manner, by normalizing a nightly co-added counts profile to an altitude where the density is reasonably well known. This density value is obtained either from some other available measurement or a model at the seeding altitude ${z}_{0}$. This choice introduces an additional uncertainty in the retrieval process. One way of avoiding this uncertainty is to seed with temperature instead of pressure. Using Eq. (13), Eq. (14) can be written in terms of the noise-subtracted lidar signal and the density at the seeding altitude ${z}_{0}$. Next the ideal gas law can be utilized to write the expression in terms of a seed temperature, ${T}_{0}$, in the following manner:

*et al.*[6] have previously shown Eq. (15) is an equivalent form of the forward model, with the advantage of being independent of the constant $C$. Hence, one source of uncertainty can be removed by using Eq. (15) as a forward model instead of using Eq. (14). Note that the retrieved temperatures in this form of the forward model are dependent only on a relative density (counts) profile and hence, no assumption of a density at the seeding altitude has to be made.

A seed temperature value is expected to introduce the same amount of systematic uncertainty, as does the combination of seed pressure and seed density if both seed pressure and density are obtained from a model, which uses the ideal gas law to calculate temperature, pressure, or density from measured values of the other two quantities. Such a situation exists for the current Purple Crow Lidar (PCL) system; thus, it will not matter which of the two equations, Eq. (14) or Eq. (15), is used for the retrieval. The temperature retrievals shown in the following sections use Eq. (14) as a forward model.

The PCL’s Rayleigh-scatter system, as configured for this study, collects collects 1200 laser shots per minute in 24 m altitude range bins to form single minute count profiles throughout the night [7]. For a nightly average temperature profile, the single minute scans are co-added together. To obtain a value of the normalization constant $C$, the nightly co-added profile, after background correction, is scaled to a density model. In this study the CIRA 1990 atmospheric model [8] was used as a standard density model. Normalization of the density profiles is done by integrating the model and measurements between 45 and 60 km. The covariance matrix, as previously described, is constructed from the night’s single minute scans and then used with the inversion approach on the nightly co-added measurements to retrieve a single average temperature profile for a night.

In this study the seed pressure value ${P}_{0}$ was also obtained from the CIRA model [8]. However since the starting altitude is in the range of radiosonde measurements($\approx 35\text{\hspace{0.17em}}\mathrm{km}$), a radiosonde-derived density can be used if available for a more accurate seed pressure.

The output from the CH method was used as an initial temperature profile for the inversion. Atmospheric models can also be used for an initial profile, but convergence to the same solution takes more iterations. In the present analysis the step size $\mathrm{\Delta}{T}_{i}$ was chosen to be 2 K. This choice was made after experimenting with other possible values for $\mathrm{\Delta}{T}_{i}$ by performing the search process and selecting the optimum value. At each iteration the step sizes were divided by the current number of iterations.

All of the temperatures retrieved for this study were spatially and temporally co-added to get nightly averaged temperature profiles. Co-adding counts increases the SNR, which helps extend the upper altitude limit of integration. Co-adding was done over a 500 m vertical resolution. In the CH method, using PCL data, it is customary to start integration, after co-adding, from the altitude where the SNR becomes approximately equal to 2 to separate data from noise [9].

Last, the “Condition for Convergence” is set so that the ${\chi}^{2}$ per degree of freedom, between two successive iterations, changes only by 0.1% or less, i.e., if at the $t$th iteration, for all the parameter values $\mathbf{T}$, $({\chi}^{2}(t-1)-{\chi}^{2}(t))/{\chi}^{2}(t)\times 100\%\le 0.1\%$ then the search process stops with the current value of $\mathbf{T}$ as the final value.

## 5. Error Analysis

Since the Grid Search method is numerical in nature, analytical uncertainties cannot be found for the random uncertainties as in the CH method and thus Monte Carlo techniques are used. A detailed description of the propagation of uncertainties using the Monte Carlo method is published by the Joint Committee for Guides in Meteorology [10,11]. The Monte Carlo technique employed for this analysis takes the final converged upon count profile, ${N}_{\text{model}}$, to which noise is added and the inversion approach is repeated on the new noisy count profile to obtain a corresponding temperature profile. These Monte Carlo runs were repeated a set number of times, and the error in the temperature profile is the standard deviation of all the temperature profiles.

Over a night’s measurements the distribution of total uncertainties for single-minute lidar measurements tend to a Gaussian distribution for similar geophysical conditions. The assumption of a Gaussian distribution is adequate here for the random uncertainties as the mean count rate even at the greater heights is sufficiently large that a Poisson distribution is indistinguishable from a Gaussian distribution [2]. This condition provides a means to model the uncertainty budget by using a Gaussian random number generator and then multiplying the output by the uncertainty in the converged upon lidar returned count profile ${N}_{\text{model}}$.

For the systematic uncertainties, beginning with ${N}_{\text{model}}$, the variance, ${\sigma}_{N}^{2}$, is taken appropriately for Poisson counting statistics. Next any systematic corrections that were previously applied to the raw measurements to retrieve ${N}_{\text{exp}}$ must also be taken into account to correct the variance ${\sigma}_{N}^{2}$. Such corrections could include the sky background, detector nonlinearity, atmospheric ozone absorption, or any other parameter specific to the lidar system. For example, the background counts are typically determined from the mean counts in an altitude region selected as the background range. To model the correct uncertainty, the mean of the background, ${N}_{BG}$, is added back to the variance such that

Similar steps can be taken for other systematic uncertainties applied in order to obtain an accurate estimate for the uncertainty of ${N}_{\text{model}}$ for uncorrelated uncertainties. With the necessary corrections to the variance completed, the noisy count profiles needed for the Monte Carlo error analysis can be generated by multiplying the error at a given altitude by a random normalized Gaussian distribution generator and adding this quantity to ${N}_{\text{model}}$ such that

where ${N}_{MC}$ is the generated Monte Carlo count profile, and $\lambda $ is the normalized random Gaussian number generator. If a large number of noisy counts profiles are generated by adding Gaussian noise to the counts profile using the aforementioned method, the known analytic uncertainty and the standard deviation from all the Monte Carlo generated profiles can be compared. It was found that the reported uncertainties for both were in agreement for the entire profile. In addition, with the developed Monte Carlo method, it is possible to isolate systematic and random uncertainties and model the effect on the final uncertainty in the temperature profile. For the PCL, it was found that 50 runs were sufficient for the Monte Carlo error analysis to converge.## 6. Results

The primary disadvantage of the CH method is that the top 10–15 km of the temperature retrieval has to be discarded due to the uncertainty associated with the choice of seed pressure, which propagates downward. With the inversion approach, a bottom seed pressure can be used, and due to the numerical optimization technique, an incorrect choice of seed pressure does not impact the retrieved temperature profile as much as the CH method. This assertion was confirmed by taking some representative nights of PCL measurements, then applying both methods and changing the seed pressure used by equal proportions and viewing the effect on the retrieved temperature profiles. The temperature profile retrieved with the initial seed pressure without any variations can be taken as the true profile for the night, while the other profiles model incorrect choices of pressure and the effect on the retrieved temperatures. Figures 1 and 2 show the results for the CH method and the inversion approach respectively for the night of 1 September 2005. The initial seed pressure used in each case was 0.00021 hPa for the top-down CH method, and 2.92 hPa for the bottom-up inversion approach. For both cases the altitude bins were co-added in bins of 504 m resolution, and then smoothed with a low-pass 3s and 5s filter [12].

Inspection of Figs. 1 and 2 shows that the temperature profiles found from the CH method are rapidly divergent towards the top of the profile due to the variation in top seed pressure. One can see that a 10% error in the estimate of the top seed pressure, which is not uncommon at such heights in the lower thermosphere, can result in a difference of almost 25 K at the top and it is not until below 75 km that the uncertainty due to the guess is negligible and the profiles converge. The results shows why it is necessary to remove the top 10–15 km of a nightly temperature profile so that the uncertainty due to the seed pressure is within the statistical uncertainty for temperatures. For the inversion approach, the maximum difference due to the choice of seed pressure is less than 5 K. The rippling effects seen are a numerical effect of the Grid Search method where the temperatures are perturbed to minimize the objective function. It is noteworthy that the largest difference at the top for the inversion approach is smaller than the difference for the CH method even after discarding the top 10 km. It can be concluded that temperature retrieval due to variations in the selection of a bottom pressure for the inversion approach has a smaller effect than the variation of seed pressure in the CH method. The removal of the top scale height of temperature retrievals is not necessary for the inversion approach, as the differences due to an incorrect seed pressure uncertainty are not as prominent. This ability to keep the retrievals at the greatest heights is the primary advantage of using an inversion approach.

To validate the temperature uncertainties found using the Monte Carlo model described in the previous section, the results from this model were compared on several nights to the analytic results given by the CH method for the statistical uncertainties. Figure 3 shows a comparison between the methods for the measurements from Sept. 1, 2005, shown in the previous figures. The statistical temperature uncertainties retrieved from the Monte Carlo model using the Grid Search method are in very favorable agreement with the CH method, giving us confidence in our technique.

Results for the nightly temperature retrieval with their respective uncertainties using the inversion approach are presented for three selected nights: 22 August 1995, 4 August 2000, and 5 May 2006 (Figs. 4–6). The red portion on each profile shows the extended temperatures that the Grid Search method can obtain over the CH method, where the top portion has to be discarded. The errors were retrieved through the Monte Carlo technique with a total of 50 runs for each night and compare well with the CH method.

## 7. Conclusions

Using mathematical inversion techniques allows Rayleigh-scatter temperature profiles retrieved from lidar measurements to be extended up to heights where previously retrievals were not possible due to the limitations in the CH method in regard to the top-down pressure integration. Gaining 10 km at the top of a lidar temperature profile requires increasing the power-aperture product of the lidar by a factor of 4. For systems such as the lidars in the Network for the Detection of Atmospheric Climate Change (NDACC) that typically already employ high-power laser transmitters, this change would require doubling the telescope diameter, which is typically prohibitively expensive. Hence, applying an inversion approach is equivalent to a significant hardware upgrade to an existing lidar, with the added benefit that past as well as future measurements enjoy the increase in effective power-aperture product.

Implementing an inversion approach for temperature retrievals offers two significant advantages compared to the CH method. First, the bottom upwards integration allows the choice of a seed pressure in the stratosphere rather than the mesosphere or the lower thermosphere. The geophysical variability at all spatial-temporal scales is much smaller in the stratosphere than the upper atmosphere. Hence, the systematic uncertainty due to the choice of seed pressure is much smaller. In addition, an inversion approach allows this uncertainty to be propagated over all heights in the model, contributing less to the systematic uncertainty near the top of the retrieval compared to the CH method.

Second, an inversion approach allows for a more thorough, albeit more complicated, uncertainty analysis. The use of a Monte Carlo technique allows systematic as well as random uncertainties to be estimated.

Our next efforts will be twofold. We are working on how best to apply this technique when higher temporal resolution scans are required during a measurements period. We also will be reprocessing our Rayleigh-scatter temperature climatology (originally given by Argall and Sica [9]) to extend the climatology to significantly greater heights.

We would like to acknowledge support of this research by the Natural Sciences and Engineering Research Council of Canada. We would also like to thank the reviewers for their insightful and helpful comments.

## References

**1. **A. Hauchecorne and M. L. Chanin, “Density and temperature profiles obtained by lidar between 35 and 70 km,” Geophys. Res. Lett. **7**, 565–568 (1980). [CrossRef]

**2. **P. R. Bevington and D. K. Robinson, *Data Reduction and Error Analysis for the Physical Sciences*, 2nd ed. (McGraw-Hill, 1992).

**3. **P. C. Mahalanobis, “On the generalized distance in statistics,” Proc. Natl. Inst. Sci. India **2**, 49–55 (1936).

**4. **V. A. Kovalev and W. E. Eichinger, *Elastic Lidar: Theory, Practice and Analysis Method* (Wiley, 2004).

**5. **P. B. Russell and B. M. Morley, “Orbiting lidar simulations. 2: density, temperature, aerosol, and cloud measurements by a wavelength combining technique,” Appl. Opt. **21**, 1554–1563 (1982). [CrossRef]

**6. **J. P. Thayer, N. B. Nielsen, R. E. Warren, C. J. Heinselman, and J. Sohn, “Rayleigh lidar system for middle atmosphere research in the arctic,” Opt. Eng. **36**, 2045–2061 (1997). [CrossRef]

**7. **R. J. Sica, S. Sargoytchev, P. S. Argall, E. F. Borra, L. Girard, C. T. Sparrow, and S. Flatt, “Lidar measurements taken with a large-aperture liquid mirror. 1. Rayleigh-scatter system,” Appl. Opt. **34**, 6925–6936 (1995). [CrossRef]

**8. **E. L. Fleming, S. Chandra, J. J. Barnett, and M. Corney, “Zonal mean temperature, pressure, zonal wind and geopotential height as functions of latitude,” Adv. Space Res. **12**, 11–59 (1990).

**9. **P. S. Argall and R. J. Sica, “A comparison of Rayleigh and sodium lidar temperature climatologies,” Ann. Geophys. **25**, 27–35 (2007). [CrossRef]

**10. **JCGM, “Evaluation of measurement data: guide to the expression of uncertainty in measurement,” Tech. Rep. (Joint Committee for Guides in Meteorology, 2008).

**11. **JCGM, “Evaluation of measurement data,” supplement 1 to the “Guide to the expression of uncertainty in measurement propagation of distributions using a Monte Marlo method,” Tech. Rep. (Joint Committee for Guides in Meteorology, 2008).

**12. **R. W. Hamming, *Digital Filters* (Prentice-Hall, 1977).