Retrieval of temperature from a multiple-channel Rayleigh-scatter lidar using an optimal estimation method.

Retrieval of temperature from a multiple-channel Rayleigh-scatter lidar using an optimal estimation method R. J. Sica1,2,3,* and A. Haefele1 1

Federal Office of Meteorology and Climatology MeteoSwiss, Payerne, Switzerland

2

Laboratory of Environmental Fluid Mechanics and Hydrology, École Polytechnique Fédérale de Lausanne, Switzerland 3

Department of Physics and Astronomy, The University of Western Ontario, London, Canada *Corresponding author: [email protected] Received 21 August 2014; revised 23 December 2014; accepted 29 December 2014; posted 6 January 2015 (Doc. ID 221397); published 3 March 2015

The measurement of temperature in the middle atmosphere with Rayleigh-scatter lidars is an important technique for assessing atmospheric change. Current retrieval schemes for this temperature have several shortcomings, which can be overcome by using an optimal estimation method (OEM). Forward models are presented that completely characterize the measurement and allow the simultaneous retrieval of temperature, dead time, and background. The method allows a full uncertainty budget to be obtained on a per profile basis that includes, in addition to the statistical uncertainties, the smoothing error and uncertainties due to Rayleigh extinction, ozone absorption, lidar constant, nonlinearity in the counting system, variation of the Rayleigh-scatter cross section with altitude, pressure, acceleration due to gravity, and the variation of mean molecular mass with altitude. The vertical resolution of the temperature profile is found at each height, and a quantitative determination is made of the maximum height to which the retrieval is valid. A single temperature profile can be retrieved from measurements with multiple channels that cover different height ranges, vertical resolutions, and even different detection methods. The OEM employed is shown to give robust estimates of temperature, which are consistent with previous methods, while requiring minimal computational time. This demonstrated success of lidar temperature retrievals using an OEM opens new possibilities in atmospheric science for measurement integration between active and passive remote sensing instruments. © 2015 Optical Society of America OCIS codes: (010.3640) Lidar; (120.6780) Temperature; (120.0280) Remote sensing and sensors. http://dx.doi.org/10.1364/AO.54.001872

1. Introduction

Accurate determination of temperature change is of great importance in atmospheric science. To assess change requires that measurements be fully characterized for random and systematic uncertainties, as these changes can be less than 0.5% of the mean temperature and occur over decadal periods. One of the most important techniques for obtaining temperature in the middle atmosphere uses backscatter 1559-128X/15/081872-18$15.00/0 © 2015 Optical Society of America 1872

APPLIED OPTICS / Vol. 54, No. 8 / 10 March 2015

measurements obtained from Rayleigh-scatter lidars. This study introduces a methodology that allows, for the first time, complete characterization of Rayleigh-scatter lidar temperature retrievals to be determined on a profile-by-profile basis. Temperature retrieval from Rayleigh-scatter lidar measurements is a well-established technique in weather and climate studies of the middle atmosphere, both from the surface [1,2] and from space [3]. To summarize, a Rayleigh-scatter lidar (hence, lidar; other types of lidar will be described if required) is typically monostatic, meaning the laser and telescope are colocated so the lidar measures

return photocounts, which are directly backscattered into the telescope. The entire atmosphere can be assumed to obey the ideal gas law, and, above approximately 25 to 30 km altitude, the atmosphere is sufficiently optically thin to visible radiation that its transmission with height is essentially constant. Hence, the return photocounts are directly proportional to the atmospheric density. Hauchecorne and Chanin [4] recognized that, if one assumes hydrostatic equilibrium (HSEQ), the photocount profile can then be integrated to obtain a temperature profile, provided a pressure is assumed at the top of the atmosphere to “seed” the integration. The method of Hauchecorne and Chanin (hereafter, the HC method) or a variation thereof is used by systems which measure Rayleigh scattering height profiles to retrieve temperatures. At about the time that the HC method was developed, atmospheric scientists began applying, in particular with space-based passive measurements, inverse methods to retrieve profiles of composition and temperature. With these methods, a forward model has been developed, which captures all relevant physical and instrumental properties of the measurement. The forward model is then inverted to allow certain model parameters, which are unknown to be estimated within a space of states defined based on prior knowledge. One of the most robust and popular of the inversion methods is the optimal estimation method (OEM), which is explained in detail by Rodgers [5]. The OEM appeals to Bayes’ theorem, which yields the most likely state and which, under the assumption of Gaussian errors, is most consistent with the measurement and prior knowledge. A cost function is computed, and the state parameters vary until minimum cost is achieved. While OEM has become a standard tool in atmospheric remote sensing using passive detectors, it has seldom been utilized for lidar measurements. Khanna et al. [6] [henceforth, the grid search (GS) method] applied an inversion approach to the retrieval of temperature from lidar measurements. Their method has several advantages over the HC method, including the ability to integrate the temperature profile from the bottom instead of the top and relaxing the assumption of constant temperature between atmospheric layers. However, their method still relies on the assumption of HSEQ, and the only uncertainties considered are statistical and those due to the assumed seed pressure. Recently, Povey et al. [7] applied an OEM to lidar retrievals of tropospheric aerosols. Their paper gives an excellent overview of basic formalisms of the OEM technique. They show results consistent with traditional analysis techniques as well as improvements in the determination of uncertainties. In this paper, we apply the OEM to the retrieval of temperatures from lidar measurements, which allows complete inventory of the systematic uncertainties associated with the temperature retrieval to be computed on

a profile to profile basis. The temperature can be retrieved with or without requiring the assumption of HSEQ. The method quantitatively determines both the vertical resolution of the retrieval as a function of height and maximum height to which the retrieved profile is independent of the a priori temperature profile, eliminating the need to arbitrarily cut off the top of a temperature profile due to choice of seed pressure. 2. Temperature Retrieval Using Rayleigh-Scatter Lidar Measurements

In this section, the basic assumptions and corrections required to retrieve temperatures from lidar backscatter profiles are discussed. A. Lidar Equation

All retrieval methods for temperatures from Rayleigh-scatter lidars use the same form of the lidar equation to equate the measured photocount profiles to number density. This method is generally considered to be valid down to altitudes above the stratospheric aerosol layer, nominally 25 km. The top height of the retrieval is principally limited by the power-aperture product of the lidar system. The lidar equation in this region can be written in the following compact form: N t z ψz ·

nz B; z2

(1)

where N t z are the true photocounts as a function of altitude, z, and nz is the atmospheric number density as a function of altitude. The units of N and B are counts per height range bin per time and, if required, can be converted to a count rate (Hz). The instrument function, ψz, depends on the following quantities: 1. Efficiency of the system, accounting for detector efficiency and nonlinearities as well as losses in the optics. 2. Atmospheric transmission, assumed constant in the height range of the retrievals. 3. Number of transmitted photons. 4. Rayleigh scatter cross section, which becomes height dependent in the upper mesosphere. 5. Area of receiving telescope. B. Using the OEM to Determine the Dead-Time Correction

At high count rates, the lidar detector response becomes nonlinear. The nonlinearity is commonly modelled by the following equation, which relates the observed count rate, N o , to the true count rate: N o N t e−N t γ ;

(2)

where γ is the counting system dead time. In traditional analysis methods, the lidar signals are corrected for nonlinearity prior to being processed [8,9]. The dead time is determined by comparing two lidar signals in an altitude region, where one 10 March 2015 / Vol. 54, No. 8 / APPLIED OPTICS

1873

signal is linear while the other is not. The two signals can be an analog and a photon count signal or two photon count signals with different attenuations (high and low range). It will be shown that, for a suitable measurement configuration, the OEM method can derive the system dead time using Eq. (2). C.

Equation of State

All three methods use the same equation of state for the atmospheric gas, the ideal gas law. The ideal gas law allows the number density in the lidar equation to be related to temperature and pressure. The number density of the atmospheric gas, even at the surface, is sufficiently small enough for the gas to be considered ideal. D.

Hydrostatic Equilibrium

Hydrostatic equilibrium (HSEQ) is the balance of the upward pressure gradient force with the downward force of gravity on the gas: dp −gzρz; dz

(3)

Effect of Absorption by Ozone and Air

For temperature retrieval in the upper stratosphere and above, it is reasonable to assume that the atmospheric transmission is constant. Between 25 and 35 km, there are small effects due to absorption of ozone and, to a lesser extent, air. Ozone absorption in the Chappuis band can cause on the order of 0.5 K changes in temperature below 35 km [12]. In this study, the ozone model of McPeters et al. [13] is used for the ozone density, and the cross section is from Griggs [14]. Operationally, using the best ozone profile available would be recommended (for instance, using colocated ozone lidar or microwave radiometer, ozonesonde, or a model). The Rayleigh extinction is calculated using the cross section, as given by Nicolet [15], which, along with the a priori density profile, is used to calculate the extinction. For a Rayleigh-scatter lidar transmitting at 532 nm, the uncertainty introduced by using the a priori density is less than 0.05 K. If desired, the 1874


F.

Lidar Instrument Function

The variation of the lidar instrument function with height is small enough that it is often assumed constant, particularly for measurements in the upper stratosphere and lower mesosphere, the retrieval range of many systems. Assuming the atmospheric transmission is constant above 35 km, this assumption is reasonable. Below this height, there is a small change in transmission due to ozone and air. Above 80 km, the mean molecular mass and, hence, the Rayleigh-scatter cross section begin varying with height due to an increase in the relative amount of atomic oxygen, and these variations can affect Rayleigh-scatter temperature retrievals [16]. The form of the instrument function used in the OEM retrievals is ψz C · σ R z · e−2·τO3 z · e−2·τR z :

where the left-hand side of the equation is the pressure gradient and the right-hand side is the gravitational force on an air parcel. HSEQ is the critical assumption used by current methods of determining temperature from the measured photocount profile. The assumption of HSEQ is reasonable over much of the atmosphere, particularly the troposphere and lower stratosphere, when sufficiently long time periods and large spatial scales are considered. However, the upper stratosphere, mesosphere, and thermosphere typically support waves of long to short spatial scales that, over short time and spatial scales, are clearly not in HSEQ. Examples of these differences on individual nights and in a climatological sense have been shown by Sica and colleagues [10,11]. E.

measured lidar density could be used for the correction.

(4)

Here, C is the lidar constant, which includes the efficiency of the system, the number of transmitted photons, and the area of the receiving telescope; σ R z, the Rayleigh-scatter cross section; τO3 z, the ozone optical depth at the wavelength of the transmitter; and, τR z, the optical depth due to Rayleigh extinction. The ozone absorption is found from the ozone density, nO3 z, and cross section, σ O3 z: Z τO3 z

z

nO3 z0 · σ O3 z0 dz0 :

0

(5)

The optical depth due to Rayleigh extinction depends on the total atmospheric density and Rayleighscatter cross section, σ R z: Z τR z

z 0

nz0 · σ R z0 dz0 :

(6)

The lidar constant is the constant of proportionality between the measured photocounts and number density. The corrected photocounts are normalized in the clear signal region, here defined as a region of the count profile where the photomultipliers are linear and the signal-to-noise ratio is high (e.g., 55 to 60 km). The total number of counts in the clear signal region is then normalized to the integrated CIRA-86 number density [9]. The Rayleigh-scatter cross section is assumed constant below 80 km. Above 80 km, the cross-section value is calculated using the cross sections for N2, O2 , and O given by Argall [16], while the ratio of each species with respect to the total number density is taken from the US Standard Atmosphere [17]. This variation of the cross section with height becomes significant in the upper mesosphere and above. The systematic uncertainties introduced in the retrieval by this form of the instrument function will be discussed in the next section.

3. Current Lidar Temperature Retrieval Methods

Most Rayleigh-scatter temperature retrievals use a variation of the HC method. In this section, we review that method, the inversion approach of the GS method, and introduce the OEM applied to the temperature retrieval problem. A.

HC Method

The HC method employs numerical integration of the lidar equation, using the ideal gas law and the assumption of HSEQ to relate the retrieved photocounts to temperature. The integration requires the assumption that, within each retrieval layer, the temperature must be constant. The numerics are robust and, even when applied to measurements with extremely low signal-to-noise ratio, can retrieve reasonable temperature estimates. Statistical uncertainties due to the Poisson noise associated with the photocount measurements are fully treated. The method does not address systematic uncertainties. The intrinsic systematic uncertainty associated with this method is related to the integral, which correlates the pressure at the seeding altitude with the retrieved pressure at every height below. Hauchecorne and Chanin [4] showed that this uncertainty decreases proportionally to the density at the top height to the density at the each height below. Thus, in practice, when using the HC retrieval, the top 10 to 15 km (on the order of two scale heights) of the temperature profile should be discarded if the seeding pressure is uncertain, which is typically the case if it is chosen from an atmospheric model. Examples of this convergence with decreasing altitude are shown in Khanna et al. [6]. Other systematic uncertainties include uncertainties in the correction of the photocount profiles at both high signal rates (pulse pile up) and low signal rates (background correction), changes in atmospheric transmission due to Rayleigh and ozone extinction, variations in the mean molecular mass with altitude, uncertainties in the determination of the pressure, lidar constant, and the height dependence of the acceleration due to gravity. Traditionally researchers have attempted to correct their photocount profiles for some of these effects. For example, simultaneous measurements of vibrational Raman scattering can be used to correct for aerosols [18,19]. Keckhut et al. [20] give a detailed discussion and estimation of systematic uncertainties due to pulse pile up, background, geometric effects, and aerosols. These and other lidar groups have made careful measurements and estimates of various systematic uncertainties, primarily those due to aerosols and the photomultipliers (including their associated counting electronics). B.

GS Method

The GS method was the first attempt to use a statistical inversion approach to retrieve temperature. A grid search method was used to find a photocount profile that minimized the cost function relative to

the measured profile, though the cost function used did not include a regularization term and, thus, was not sensitive to departures of the solution from the a priori state. The forward model equation used was Z p0 ψ 1 M ztop gz0 0 No 2 dz B; (7) · exp − R z Tz0 z Tz where the instrument function ψ is height independent, p0 is the pressure at the reference altitude, M is the height-independent mean molecular mass, R is the gas constant, and Tz the temperature profile. Note that using this equation introduces two further unknowns into the analysis, the mean molecular mass, and the acceleration due to gravity. This forward model assumes HSEQ. The retrieved photocount profile can then be integrated in a manner similar to the HC method with two important differences. First, the temperature within a layer can vary with altitude. Second, the profile can be integrated from the bottom to the top of the retrieval grid, which has the advantage of minimizing seed pressure effects at all heights, thus, greatly reducing the seed pressure uncertainty at the greatest heights. Khanna et al. [6] showed that the systematic uncertainty associated with this integration from the bottom allowed the entire retrieval grid to be used without having to discard the top of the profile. The GS method does not calculate the uncertainties of the retrieved model parameters. The statistical uncertainty is found by a Monte Carlo simulation. Khanna et al. [6] showed that the statistical uncertainties found in this manner are fully consistent with the HC method; however, the computational time for the Monte Carlo simulation is long compared to the HC method. 4. Application of the OEM to Rayleigh-Scatter Lidar

While the OEM is well known in the atmospheric satellite and passive optical measurement community, it may not be familiar to the lidar specialist. Here, the specific application of OEM to Rayleigh-lidar temperature retrieval is discussed. A. OEM

The OEM offers significant advantages over the HC and GS methods: 1. Forward model does not require the assumption of HSEQ. 2. Allows direct determination of photomultiplier dead time to be found using the full nonlinear correction equation for a paralyzable system. 3. Allows physically and mathematically justified determination of the top of the retrieved profile. 4. Full uncertainty budget can be retrieved, including all systematic and statistical uncertainties on a profile by profile basis. 5. No post- or pre-filtering of the retrieval is required, and the height resolution of the retrieved temperature profile is the FWHM of the averaging kernels, which are computed at each height. 10 March 2015 / Vol. 54, No. 8 / APPLIED OPTICS

1875

6. Detector and sky background (either constant or height dependent) are retrieved along with temperature. 7. Can be applied at any required height or time resolution, e.g., nightly averaged profiles that are co-added in height or individual profiles at high temporal-spatial resolution. 8. Method is extremely fast computationally. These advantages will be explored in detail in the following sections. B.

The gain is the sensitivity of the retrieval to the measurement, or equivalently, to the measurement error: Gy

∂ˆx : ∂y

(11)

The averaging kernel gives the sensitivity of the retrieved state to the true state and can be shown to be equal to the product of the gain and the Jacobian:

Brief Description of the OEM

To apply the OEM, a forward model is required. The forward model contains all the physical and instrumental factors, which together describe the measurements. The general form of the forward model F is y Fx; b ϵ;

(8)

where y is the measurement vector, which is specified on the measurement grid. The forward model depends on the state vector, x, and the model parameters, b. The state vector and the model parameters describe the atmosphere and the instrument measuring it. The state vector contains all the parameters that are to be retrieved, while the model parameters contain all of the other parameters needed to fully model the measurement but which are not retrieved. Specific states include the true state (typically unknown except for synthetically generated states), xt , the a priori state, xa , and the retrieved state, xˆ . The measurement noise is given by ϵ. Finding the retrieved state from the measurements is made possible by Bayes’ theorem, which is the basis for the manipulation of conditional probabilities. In this case, the a priori state, xa , consists of a temperature profile, background, and, for the two-channel case, dead time. A probability is assigned reflecting the certainty of this state. There is also a measurement and its associated variance. Bayes’ theorem allows the most likely state to be found consistent with the a priori knowledge, the performed measurement, and the associated uncertainties. The optimum estimate for the retrieved state is found by minimizing a cost function formed from the measurement, y, and its covariance, Sy , the forward model, the retrieved state model parameters, and the a priori covariance, Sa , using Bayes’ theorem: cos t y − Fˆx; bT S−1 x; b y y − Fˆ x − xa : ˆx − xa T S−1 a ˆ

(9)

The retrieval is characterized by three matrices. The Jacobian gives the sensitivity of the forward model to the state: ∂F Kx : ∂x 1876


(10)

A

∂ˆx Gy Kx : ∂x

(12)

For the lidar temperature retrieval, the problem is nonlinear; hence, the retrieved state is found iteratively using the Marquardt–Levenberg method to find a solution that minimizes the cost function. Note that, in the practical implementation of the method, the vertical grid of the state vector is different from the measurement grid. Consequently, the forward model performs interpolations between the two grids, which are not explicitly shown in the following equations for the sake of brevity. All interpolations used are linear interpolations. Linear interpolation is sufficient, as the measurement and retrieval grid spacing is small relative to the scale height of the atmosphere. When pressure interpolation is required, the interpolation is performed on the logarithm of the pressure in order to ensure that the interpolation does not introduce fluctuations into the retrieval. The uncertainty budget is determined from the measurement, smoothing, and model parameter covariance matrices [5]. The retrieval’s covariance due to measurement noise, Sm , is Sm Gy Sy GTy :

(13)

The retrieval’s covariance due to smoothing, Ss , is Ss A − In Se A − In T :

(14)

Here, In is the unit matrix of order n and Se is the ensemble of states about the mean state: Se Kb Sb KTb Sy ;

(15)

which depends on the model parameter covariance, Sb , and the measurement covariance. Model parameters can have systematic and random uncertainties. The retrieval’s covariance due to the forward model parameters, SF , is SF Gy Kb Sb KTb GTy :

(16)

The total covariance, Stotal , is then Stotal Sm Ss SF :

(17)

N o z N t e−N t γ :

C. Forward Models for Rayleigh-Scatter Temperature Retrieval

After extensive experimentation, two forward models were chosen for further study. In both models, the retrieved parameters are temperature, dead time, and background (which includes photomultipler shot noise and sky background and, for some systems, depends on height). If only a single channel is available, dead time cannot be retrieved, but dead time is then treated as a forward model parameter. The lidar equation (LE) model is based solely on the lidar equation, with the number density replaced by pressure and temperature via the ideal gas law. For a linear system, the low-gain digital or analog Rayleigh [low-level Rayleigh (LLR)] channel forward model is N o z

ψz pz B: z2 kTz

(18)

For a channel with a height-dependent background, B would then be a function of altitude, and the additional coefficients would be included in the forward model. For the high-gain Rayleigh [high-level Rayleigh (HLR)] channel, the true count profile is transformed into the observed height profile using the dead-time correction in Eq. (2): N t z

ψz pz B; z2 kTz

N o z N t e−N t γ :

(19)

(20)

Unlike the GS method, the LE forward model has a height-dependent instrument function ψz. The model parameters here are the instrument function, the dead time, and a pressure profile. It is not possible for the instrument function and pressure to be included with temperature and background as retrieval parameters, as the Jacobians for the instrument function would not be linearly independent. The LE model does not assume HSEQ; thus, the pressure profile must be specified at all points on the retrieval grid. The second forward model is the OEM analog of the HC method, where HSEQ is assumed and combined with the lidar equation to give a forward model similar to that given for the GS method: Nt

ψz phseq z B; z2 kTz

(21)

where Z 1 ztop Mz0 gz0 0 phseq z p0 exp − dz : R z Tz0

(22)

As for the LE model, the true counts are related to the observed counts by

(23)

This second forward model, hereafter called the HSEQ forward model, has the same heightdependent instrument function as the LE forward model. The assumption of HSEQ introduces two new model parameters: the height-dependent mean molecular mass, Mz, and the acceleration due to gravity. In the lidar equation forward model, the pressure profile is fixed during the retrieval, but, in the HSEQ model, only p0 is fixed, so the height dependence of the pressure profile can change as the retrieved temperature profile changes. In the case of the HSEQ model, the resulting pressure profile may not be the true pressure profile, unless the atmosphere is in HSEQ. The advantage of the HSEQ forward model is that pressure only needs to be known at a single point to define a pressure profile at all points. However, two additional forward model parameters are now required: the mean molecular mass and the acceleration due to gravity. Given that the HSEQ forward model and the HC method contain similar physical assumptions, temperature retrieved from the same measurements should be similar. This assertion is shown to be true in Section 5. D.

Implementation of the OEM

For pilot studies of the effectiveness of retrieving various combinations of parameters, forward models, and uncertainties, a linear OEM retrieval was written and tested on synthetic lidar measurements with and without noise, as described in [5]. Results from the pilot study will not be shown, but, for reasonable perturbations (including adding realistic atmospheric waves) and a priori states, a robust temperature retrieval was obtained, giving confidence in the method and forward models used. Parameters that could (and could not) be retrieved were also investigated, and the best combination for use with a single data channel is to retrieve the temperature and instrumental background, leaving all other quantities as forward model parameters (e.g., pressure, instrument function, gravity). This study uses the Qpack software package developed by Eriksson et al. [21]. Qpack is a free MATLAB package designed for general OEM solutions. Qpack is part of the radiometer and microwave community ARTS model [22], which has been well tested for both ground-based and satellite retrieval. Qpack is used with the Marquardt–Levenberg method, which is also used because it is a nonlinear method and is often able to retrieve parameters even when the a priori is considerably far from the true solution, e.g., the US Standard Atmosphere temperature profile compared to an individual temperature profile on a given night. As discussed in Eriksson et al. [21], there is typically some covariance between the various quantities, and the covariance matrices should have off-diagonal elements. This correlation is provided by 10 March 2015 / Vol. 54, No. 8 / APPLIED OPTICS

1877

using correlation functions with a correlation length appropriate to the physical situation. For this study, we have used tent functions, with an a priori correlation length of 3000 m for temperature. The measurement covariance matrix, Sy , is a diagonal matrix. The solutions did not show much sensitivity to any reasonable choice of the correlation length. The Marquardt–Levenberg method then finds an iterative solution that minimizes the cost function to order unity in a few (typically four to six) iterations. 5. Results and Validation of the OEM Temperature Retrieval

In this section, the OEM retrieval is tested in several ways. These tests convincingly demonstrate the robustness of the method, its advantages compared to other methods, and the validity of the temperatures retrieved. The OEM retrieval is tested with both forward models using synthetic measurements and then with Rayleigh-scatter lidar measurements. The measurements used were obtained in the spring and summer of 2012 by The University of Western Ontario’s Purple Crow Lidar (PCL) at its new location, the Environmental Sciences Western Field Station (43.1°N, 279°W, 275 m altitude). This facility is about 20 km northeast of the PCL’s original location. The PCL is a large power-aperture monostatic lidar, which is described in detail in [9], with some significant upgrades for these measurements. The transmitter is now a Litron Nd:YAG laser that outputs 1000 mJ per pulse at 532 nm with a repetition rate of 30 Hz. The photocount measurements were obtained from three channels: a Hamamatsu R7400 photomultiplier, outputted in analog and digital mode to a Licel transient recorder and a Hamamatsu R5600 photomultiplier, outputted in digital mode to a SRS 430 multichannel scalar averager. Both counters accumulated returns over the same 1800 laser shots (∼60 s), at 7.5 m height resolution for the Licel analog (LLRa) and HHR channel and at 24 m for the LLR digital (LLRd) channel. The LLRd detector is attenuated with a neutral density filter, chosen to keep the count rate linear from the time when the system’s mechanical chopper fully opens (∼22 km altitude). The retrievals are performed between 30 and 120 km over integration periods of a few hours. This range is used to allow measurements to range from an extremely high signal-to-noise ratio up to altitudes in which the photocounts can be assumed to be due only to background light. The measurements were then coadded to 255 m for the Licel channel and 264 m for the LLR (i.e., the measurement grid). The retrieval grid has a vertical resolution of 1020 m. The clear signal region is 30 to 35 km for the LLR channel and 55 to 60 km for the HLR channel. The OEM requires an estimate of the covariances associated with the measurement and retrieval parameters. Since the photocounts obey Poisson counting statistics, the measurement variances 1878


are equal to the count rate for the LLRd and HLR channels, with no correlation between the photocount measurements assumed. The a priori temperature profile used for all retrievals is the US Standard Atmosphere, with a variance of 35 K2 . This choice of variance is consistent with the discussion of variability in the US Standard Model description as well as the various data sets used in the MSIS-91 extended model [23]. The background variance is estimated from the measurement profile by determining the standard deviation of the mean of the photocount returns above 115 km. The pressure profile is assigned a standard deviation of 1%, the lidar constant 1%, the O3 density 10%, the Chappuis band cross section 2%, the dead time 0.2%, acceleration due to gravity 0.1%, and the Rayleigh scattering cross section 0.2%. Correlation was introduced in the covariance matrices for temperature and the height-dependent forward model parameters using a tent function with a length of 3 km, as described in Section 4.D. These choices and their implications to the retrieval will be further discussed as follows. A. Retrievals Using Synthetic Measurements

Initial testing of the OEM was accomplished with both forward models by using the forward models themselves to generate synthetic photocount profiles with system parameters estimated for the PCL HLR channel with an integration time of 6.5 h (similar to the measurements on 24 May 2012). Poisson noise and a constant background were then added to the photocount profiles. For these tests, the true temperature profile and pressure values were taken from the CIRA-86 model for May at middle latitudes [24]. Neither forward model requires the model atmosphere to be in the HSEQ. The LE forward model does not require the assumption of the HSEQ. The HSEQ forward model only uses the atmospheric model to pick a pressure at one point (p0 , typically the top of the grid) on the measurement grid and then uses Eq. (22) to compute a pressure profile from the atmospheric model temperatures and p0 . 1. HSEQ Forward Model Figure 1 shows the temperature Jacobians, normalized by the number density. The Jacobians are normalized in the figure to allow them to be plotted on a linear scale. The main feature of the Jacobians is the sharp peaks at the level where the temperature has been perturbed. The Jacobians also show a second feature, which grows with decreasing altitude. Figure 2 shows the Jacobians for 45 km on the retrieval grid for the HSEQ and the LE forward models. Both Jacobians are identical from the top of the data grid down to where they begin to decrease rapidly with height. The LE Jacobian is symmetric and sharply peaked around 45 km. The HSEQ Jacobian is about twice as wide at half maximum and has a tail that increases with decreasing height. The use of a HSEQ pressure assumption causes the pressure

120

110

110

100

100

90

90

Altitude (km)

Altitude (km)

120

80 70

80 70

60

60

50

50

40

40

30 -8

-6

-4

-2

30

0

Jacobian / number density (x 10

-19

0

0.5

1

Temperature Averaging Kernels

)

Fig. 1. Temperature Jacobians, normalized to the number density (left panel) and averaging kernels (right panel) for a synthetic measurement using the HSEQ forward model. The horizontal dashed line on the right panel is the height above which the a priori temperature profile makes a significant contribution to the retrieval. This height is determined from the trace of the averaging kernel matrix.

at all heights to be correlated with the temperatures at all heights above, which is a consequence of the integral in phseq z. The tail in the HSEQ Jacobian below the peak is caused by changes in pressure below the peak of the Jacobian as the temperature changes because, if temperature below a given height changes, it causes a change in the pressure at all heights above it. The averaging kernels (Fig. 1) contain two important pieces of information not available in the HC or GS method. First, the trace of the averaging kernel matrix is the number of degrees of freedom of the retrieval. For the temperature retrievals shown here, and, in subsequent sections, there are about 90 retrieval levels extending above 120 km, with degrees of freedom in the mid-70 s, meaning the

retrieval contains mostly information from the measurement and not the a priori up to heights of 90 to 100 km. The OEM method eliminates any subjective judgements about the effect of seed pressure on the retrieval: the degrees of freedom directly give the cutoff, and the uncertainty budget (see below) gives the systematic uncertainty due to seed pressure. Second, the FWHM of an individual averaging kernel is the vertical resolution of the measurements at that point, which only can be equal to or greater than the retrieval grid spacing (Fig. 3). Often, lidar measurements are reported on a measurement grid and then digitally filtered, without the effect of the filter on the vertical resolution being specified. The plot of vertical resolution shows that, over most of the retrieval’s useful range, the vertical resolution of the retrieval is identical to the resolution of the retrieval grid. Figure 4 shows the percent difference between the synthetic measurement and the final iteration of the forward model. The agreement between the measurements and model is within the Poisson photocount uncertainty and has zero bias. The difference between the true and retrieved background (not shown) is about 0.5% and is within one standard deviation of the uncertainty of the retrieved background. The retrieved temperature profile is shown in Fig. 5. Below 100 km, the agreement is within the total uncertainty (Fig. 6). This figure reveals one of the major advantages of the OEM technique: the ability to compute a full uncertainty budget of random and systematic uncertainties. The following systematic uncertainties are considered: 1. Uncertainty of the Rayleigh cross section on Rayleigh extinction.

120 110

50

100

46

90

Altitude (m)

48

Altitude (m)

44 42 40

80 70 60

38 36

50

34

40

32 30

30 10 3 -2

-1.5

-1

-0.5

Temperature Jacobian (cts/K x 10 4)

0

Fig. 2. Representative temperature Jacobian at 45 km altitude for the HSEQ forward model (blue curve) and the LE forward model (red curve).

10 4

Vertical Resolution (m)

Fig. 3. Vertical resolution of the retrieval for a synthetic measurement using the HSEQ forward model. The horizontal dashed line is the height above which the a priori temperature profile makes a significant contribution to the retrieval. 10 March 2015 / Vol. 54, No. 8 / APPLIED OPTICS

1879

120 110 100

Altitude (m)

90 80 70 60 50 40 30 -10

-5

0

5

10

Percent Difference

Fig. 4. Residuals between the HSEQ forward model and synthetic measurements.

2. Uncertainty in the total density used to calculate the Rayleigh extinction. 3. Uncertainty in the ozone absorption cross section of the Chappuis band. 4. Uncertainty in the ozone concentration in the height range of the retrieval. 5. Uncertainty in the knowledge of the lidar constant. 6. Counting system dead time correction (HLR only). 7. Variation of Rayleigh-scattering cross section with height. 8. Pressure (either at one point for HSEQ forward model or a profile for LE forward model). 9. Acceleration due to gravity (HSEQ forward model only). 10. Variation of mean molecular mass with altitude (HSEQ forward model only). 120 110 Retrieved A priori True

100

Altitude (km)

90 80 70 60 50 40 30 150

200

250

300

350

Temperature (K)

Fig. 5. Temperature retrieval from synthetic measurements using the HSEQ forward model. The horizontal dotted line is the height above which the a priori temperature profile makes a significant contribution to the retrieval. 1880


Fig. 6. Uncertainty budget for the HSEQ forward model synthetic temperature retrieval. Left panel: Statistical (blue), smoothing (red), pressure (yellow circle), lidar constant (magenta x), ozone density (green +), ozone cross section (blue *), density profile for Rayleigh extinction (red square), Rayleigh extinction cross section (blue diamond), dead time (orange triangle), variation of Rayleighscatter cross section with height (yellow triangle), variation of mean molecular mass with height (purple triangle), gravity (green triangle), and total uncertainty (black line). The horizontal dashed line is the height above which the a priori temperature profile makes a significant contribution to the retrieval. Right panel: Expanded view of the uncertainties in the stratosphere.

The OEM also computes the smoothing error, which estimates the uncertainty due to the smoothing by the averaging kernels on the retrieved state, causing it to miss features that may be present in the true state. For the lidar retrievals, the smoothing error becomes important in the lower thermosphere, as the true count rate drops to the order of a few times the background. At the greatest heights, the variation of the mean molecular mass with altitude also begins to contribute. The random uncertainties are less than 0.8 K in the stratosphere. The total uncertainty is due primarily to the assumption of a seed pressure (here at 120 km), the uncertainty in the lidar constant, and knowledge of the local gravitational field. Small contributions result from uncertainties in the counting system dead time and ozone density. Random uncertainty is dominant until around 90 km, where the smoothing error begins to rapidly increase. The statistical uncertainties in the OEM method are approximately the same as the statistical uncertainties determined from the HC method. Hence, the statistical uncertainties in the OEM retrieval are similar to the HC method, as opposed to the GS method, where the statistical uncertainties are more uniform with height compared to the HC method [6]. Small contributions to the uncertainty budget arise from variations in the mean molecular mass at the greatest heights, the lidar constant, the seed pressure, and, at the lower heights, Rayleigh extinction and the dead time correction.

120 110 100 90

Altitude (m)

2. LE Forward Model Figure 7 shows the Jacobians and averaging kernels for the LE forward model. Without the assumption of the HSEQ, the Jacobians are more localized, with no contributions from lower heights as seen in the HSEQ case above. The vertical resolution is about 15% better at 100 km than for the HSEQ model (Fig. 8). The model fit is excellent, with small residual differences between the forward model and measurements on the order of the Poisson photocount noise (Fig. 9), leading to a successful temperature retrieval (Fig. 10).

80 70 60 50 40 30 -10

120

120

-5

0

5

10

110

110

100

100

90

90

Altitude (km)

Altitude (km)

Percent Difference

80 70

Fig. 9. Residuals between the LE forward model and the synthetic measurements.

120 80

110 70

Retrieved A priori True

100

60

60

50

50

40

40

30 -8

-6

-4

-2

0

30

Jacobian / number density (x 10-19 )

0

0.5

1

Temperature Averaging Kernels

Fig. 7. Temperature Jacobians, normalized to the number density (left panel) and averaging kernels (right panel) for a synthetic measurement using the LE forward model. The horizontal dashed line on the right panel is the height above which the a priori temperature profile makes a significant contribution to the retrieval. This height is determined from the trace of the averaging kernel matrix.

120

Altitude (km)

90 80 70 60 50 40 30 150

200

250

300

350

Temperature (K)

Fig. 10. Temperature retrieval from synthetic measurements using the LE forward model. The horizontal dotted line is the height above which the a priori temperature profile makes a significant contribution to the retrieval.

110 100

Altitude (m)

90 80 70 60 50 40 30

0

500

1000

1500

2000

2500

3000

3500

4000


Fig. 8. Vertical resolution of the retrieval for a synthetic measurement using the LE forward model. The horizontal dashed line is the height above which the a priori temperature profile makes a significant contribution to the retrieval.

The systematic uncertainties are significantly different between the two forward models (Fig. 11). For the LE forward model, the largest systematic uncertainty is due to the pressure profile and lidar constant, contributing between 4 to 7 K, with the maximum contribution around 45 km altitude. The contribution from the lidar constant is significant, with a similar shape to pressure. These differences are due to the assumption of HSEQ, which causes two important differences in the sensitivity of the forward models to their parameters. For pressure uncertainty, recall that, in the HSEQ forward model, the pressure is only assumed at one point (p0 ). The pressure profile is computed relative to this point, and p0 ’s contribution is primarily in the seeding of the pressure at the initial heights. The uncertainty due to p0 then primarily affects the pressures near 10 March 2015 / Vol. 54, No. 8 / APPLIED OPTICS

1881

110

48 46

100

44

Altitude (km)

Altitude (km)

90 80 70 Statistical Smoothing Pressure Constant [O ] 3 O3 Ray Exct Ray Exct Deadtime (z) Ray Total

60 50 40 30

0

10

20

Uncertainty (K)

42 40 38 36 34 32

30

30

0

0.2

0.4

Uncertainty (K)

Fig. 11. Uncertainty budget for the LE forward model synthetic temperature retrieval. Left panel: Statistical (blue), smoothing (red), pressure (yellow circle), lidar constant (magenta x), ozone density (green +), ozone cross section (blue *), density profile for Rayleigh extinction (red square), Rayleigh extinction cross section (blue diamond), dead time (orange triangle), and the variation of Rayleigh-scatter cross section with height (yellow triangle) and total uncertainty (black line). The horizontal dashed line is the height above which the a priori temperature profile makes a significant contribution to the retrieval. Right panel: Expanded view of the uncertainties in the stratosphere.

the seed height, with the initial pressure becoming less important farther from the initial point. This behavior is similar to the role of the constant of integration term in the HC method. In the LE forward model, pressure, and its associated uncertainty, is specified at every point on the retrieval grid, causing the pressure to contribute uncertainty locally at every height. The difference in uncertainty for the lidar constant can be seen by combining the HSEQ Eq. (3) and the ideal gas law to find a expression for Tz [25]. The resulting expression for the temperature profile depends on the relative density; hence, Tz is less sensitive to the choice of the lidar constant. The LE forward model depends on the absolute, not relative, density; thus, the choice of lidar constant affects the density at all heights. The uncertainty due to Rayleigh scatter crosssection variation is small. For this forward model, there is no uncertainty due to the acceleration due to gravity or the mean molecular mass, as these quantities are only used to compute the pressure profile when HSEQ is assumed. The largest uncertainty in the stratospheric temperature for the LE forward model is due to the uncertainty in the ozone cross section. By comparison, the uncertainty due to knowledge of the ozone density profile is about 30% of this amount. Small but significant contributions arise from the Rayleigh scatter cross section, dead time, number density, and, in the upper stratosphere, the statistical uncertainty of the measurement. 1882


3. Response of the OEM to Rapid Temperature Fluctuations To test the sensitivity of the OEM to rapid changes in temperature, an oscillation was introduced into the true profile that grew as a sinusoid in height with a growth rate of twice the scale height. This type of oscillation is what would be expected for a monochromatic gravity wave conserving kinetic energy density as it propagates upward. The perturbation is introduced at the lowest heights with an amplitude of 0.5 K and allowed to grow until the associated lapse rate exceeds the dry adiabatic lapse rate, after which the amplitude is fixed (“saturates”). Figure 12 shows the temperature profiles retrieved for an oscillation with a vertical wavelength of 10 km. The agreement with the true temperature profile shows that the OEM technique is capable of capturing extremely sharp features in the temperature profile. A less realistic case of a single wave with a vertical wavelength of 3 km is shown in Fig. 13. The OEM retrieval using the HSEQ forward model performs as well with these extreme perturbations, as it does in the cases without added oscillations discussed previously. The LE forward model performs equally as well as the HSEQ forward model for the same added oscillations. B. Temperature Retrievals for 24 May 2012

The results in the previous section give confidence in the OEM and forward models developed to fit the lidar measurements. In this section, these models are applied to PCL measurements on the night of 24 May 2012 and compared to retrievals using the HC method on the same lidar measurement.

120 120 110 110 100 100 Retrieved A priori True

90 80 70

90 Altitude (km)

50

Altitude (km)

120

80 70

60

60

50

50

40

40

30

150

200 250 300 Temperature (K)

350

30 −40

−20 0 Temperature (K)

20

Fig. 12. Temperature retrieval using the HSEQ forward model for a synthetic measurement, which includes a gravity-wave-like disturbance with a 10 km vertical wavelength. Left panel: Temperature profile. Right panel: Temperature difference between the retrieved and the true profiles (blue curve) and between the retrieved and a priori profiles (red curve). The horizontal dotted line is the height above which the a priori temperature profile makes a significant contribution to the retrieval.

120

120

120

50

110

110

110

48

100

100

70

80 70

80 70

60

60

50

50

60

40

40

50

30

30 −30

150

200 250 300 Temperature (K)

350

44

90

−20

−10 0 10 Temperature (K)

20

Fig. 13. Temperature retrieval using the HSEQ forward model for a synthetic measurement, which includes a gravity-wave-like disturbance with a 3 km vertical wavelength. Left panel: Temperature profile. Right panel: Temperature difference between the retrieved and the true profiles (blue curve) and between the retrieved and a priori profiles (red curve). The horizontal dotted line is the height above which the a priori temperature profile makes a significant contribution to the retrieval.

1. HSEQ Forward Model Figure 14 shows a comparison between the temperature retrieved by the OEM using the HSEQ forward model and the temperature calculated from the HC method. The a priori background is taken as 626 ct∕bin∕ min , and the retrieved background is 608 2.5 ct∕bin∕ min . Figure 15 shows the uncertainties associated with the retrieval. The first point to note is that the systematic uncertainties are largest until the greatest heights of the retrieval, where the statistical uncertainty and smoothing error account for most of the total uncertainty. Since the top 10 to 15 km would be removed from the profile, 120 110 100

Altitude (km)

80

90

Altitude (km)

Retrieved A priori True

Altitude (km)

Altitude (km)

90

46

100

Statistical Smoothing Pressure Constant [O 3 ] O 3 Ray Exct Ray Exct Ray DeadTime

40 38 36 34

Ray(z)

40 30

42

MMM Gravity Total

0

10

20

Uncertainty (K)

32 30

30

0

0.5

1

Uncertainty (K)

Fig. 15. Uncertainty budget for the HSEQ forward model temperature retrieval on 24 May 2012. Left panel: Statistical (blue), smoothing (red), pressure (yellow circle), lidar constant (magenta x), ozone density (green +), ozone cross section (blue *), density profile for Rayleigh extinction (red square), Rayleigh extinction cross section (blue diamond), dead time (orange triangle), variation of Rayleigh-scatter cross section with height (yellow triangle), variation of mean molecular mass with height (purple triangle), gravity (green triangle), and total uncertainty (black line). The horizontal dashed line is the height above which the a priori temperature profile makes a significant contribution to the retrieval. Right panel: Expanded view of the uncertainties in the stratosphere.

the HC retrieval would be reported starting around 85 km, where most of the uncertainty is statistical. However, with the OEM retrieval, a useful profile, with full uncertainty budget, is obtained up to 100 km. In addition to random uncertainties, variation in the mean molecular mass, knowledge of the lidar constant, p0 and the smoothing error all contribute to an overall uncertainty at 100 km of 15 K. In the stratosphere, the uncertainty increases to about 1 K at the bottom of the profile, with the largest contributions due to the assumed pressure at 120 km, the lidar constant, the gravitational profile, and the dead time.

Altitude (km)

90 80

Retrieved A priori HC

70 60 50 40 30 150

200

250

300

350

Temperature (K)

Fig. 14. Single channel HSEQ temperature retrieval from Purple Crow Lidar measurements on 24 May 2012 (red curve). The blue curve uses the same photocount profile to retrieve temperature using the HC method. The a priori temperature profile (cyan curve) is the US standard atmosphere. The horizontal dotted line is the height above which the a priori temperature profile makes a significant contribution to the retrieval.

2. LE Forward Model The LE forward model differs from the HSEQ forward model because a pressure profile, as opposed to a pressure seed height, is required. Temperature retrievals are performed using the LE forward model with different pressure profiles. Since the pressure profile is a model parameter, and not a retrieval parameter, utilizing the LE forward model requires the best available pressure profile to be useful, as pressure is inversely proportional to temperature. Few measurements of atmospheric pressure profiles exist, particularly on a routine basis. The Aura microwave limb sounder (MLS) provides excellent global coverage and retrieves temperature profiles as a function of pressure as well as geopotential heights, which are computed using the assumption of HSEQ [26]. Thus, pressure can be equated to 10 March 2015 / Vol. 54, No. 8 / APPLIED OPTICS

1883

altitude. Using this pressure profile measurement for an overpass during the measurement period on 24 May 2012 gives a temperature retrieval that is much hotter in the stratosphere and much cooler in the upper mesosphere than the HC retrieval. Similarly, using both CIRA-86 and the US Standard Atmosphere pressure profiles yielded temperature profiles considerably different than those computed from the HC method. To demonstrate that the LE model temperature differences were purely due to differences in the pressure profile, temperature was retrieved using a hybrid pressure approach. In the HSEQ forward model, the pressure profile is updated at each iteration in order to be consistent with the retrieved temperature profile at that iteration. The final pressure profile found from HSEQ is, thus, consistent with the retrieved temperature. When this pressure profile is used as the true pressure profile in the LE forward model, the retrieved temperature profile is equivalent to the HSEQ retrieval, and, like the HSEQ retrieval, is similar to the HC method. The reason for the large differences in temperature between the LE retrieval using the various pressure profiles and the HSEQ forward model is due to the large differences in pressure between the profiles (Fig. 16). The MLS pressure is less than the HSEQ pressure at all heights, with the difference between them almost a factor of 2 in the lower thermosphere. Both pressure profiles are significantly lower than the CIRA and US Standard models. The pressure in the stratosphere and mesosphere is greater in the CIRA model than the US Standard model. The HSEQ pressure is about 30% less than the CIRA pressure in the mesopause region. C. Dual Channel Retrieval of Temperature, Dead Time, and Background

The advantage of a two-channel system for large power-aperture product lidars is that the range over 120 CIRA HSEQ MLS

110 100

Altitude (km)

90 80 70 60 50 40 30 -60

-50

-40

-30

-20

-10

0

10

20

Pressure Difference to US Standard (%)

Fig. 16. Pressure percentage difference on 24 May 2012 for the CIRA-86 model (blue curve), pressure calculated from the temperature assuming HSEQ (red curve), and the MLS pressure profile (green curve) relative to the US standard model pressure. 1884


which temperatures can be retrieved is increased. However, even if the system is engineered so each channel is linear, some scheme must be used to combine the measurements into a single temperature profile. The most common approaches are to either form a common photocount profile from the two channels (often called “gluing”) or to retrieve temperature profiles from each channel and combine these into a single profile (often called “merging”). Each method suffers in practice from the same issue. For many nights, the slopes of the photocount profiles or the shape of the temperature profiles in the region of overlap are similar. However, for many other nights, the slopes are significantly different and gluing or merging introduces systematic uncertainty into the retrieved temperatures. Also, a quantitative criteria for what region to choose for overlap between the two channels has not been established. The gluing method also has the challenge of dealing with combining the statistical uncertainties, which are different for each channel. In an analog/digital system, there is the challenge of converting the raw digitized analog signal into virtual photocounts while maintaining appropriate statistical uncertainties. The two-channel OEM retrieval does not require solutions to any of these challenges. The measurement vector now contains the measurements from both channels, which are not required to be in the same units or even have the same spatial-temporal resolution. The retrieved temperature profile is a single profile, which minimizes the cost function by accounting for contributions from both measurements. The background for each channel is retrieved, in addition to the dead time for the HLR channel. Retrievals have also been successfully tested using analog and digital channels, using the native units for the analog signal. Extensive testing was performed using synthetic measurements to gain experience in the best combination of height ranges and overlap. The following retrievals use the appropriate HSEQ forward model. The lowest altitude for the HLR channel was chosen as 37.5 km. At this height, the nonlinearity of the detector is modest, e.g., ∼5%. Measurements were used up to 122 km, similar to the single channel case. The LLRd channel measurements are used from a lower height of 25 km to an upper height of 110 km. The normalization of the density profile is 55 to 60 km for the HLR channel and 30 to 35 km for the LLRd channel. The HLR measurements are coadded to 30 m, while the LLRd measurements are used at their native 24 m vertical resolution (Fig. 17). The retrieval grid is the same as the single channel case, 1020 m. The temporal co-adding for both channels is 6.5 h of data records consisting of 1800 laser shots each (1 min). The a priori estimate of the dead time is the manufacturer’s specification of 4 ns. Figure 18 shows the Jacobians for the two-channel retrieval on 24 May 2012. The sensitivity of the forward model at the lowest heights to temperature is larger for the HLR channel, due to its higher count

120 110

110

110

100

100

90

90

100

80 70 60 50 40 30 10 -4

10 -3

10 -2

10 -1

10 0

10 1

90

Altitude (m)

Altitude (km)

Altitude (km)

120

LLR HLR

120

80 70

60

50

50

40

40

30

30 0

Count Rate (MHz)

0.5

1

10

T Averaging Kernels

rate. The averaging kernels and vertical resolution are also similar to the one-channel case, except above 95 km (Fig. 19). In the 95 to 100 km region, the twochannel retrieval’s resolution increases more rapidly with altitude and is about 30% greater than the onechannel retrieval at 100 km. Figure 20 shows the differences between the forward model and the measurements. The LLRd measurements contribute from 25 km up to 37.5 km. Above this height, the measurement covariance weights the retrieval primarily to the HLR measurements. The handoff from the LLRd residual to the HLR residual is seamless, and the residuals are around 0.1% or less. At the greater heights, the LLRd residuals have a bias, but not the HRL residuals, which are determining the temperature profile at these heights. The HLR residuals become large above 100 km, as the measurements no longer

70

60

10 2

Fig. 17. HLR and LLR digital channel photocount measurements from the night of 24 May 2012. The count rates are typical for the PCL Rayleigh-scatter system.

80

3

10

4


Fig. 19. Averaging kernels (left panel) and vertical resolution for the two channel retrieval on 24 May 2012 using the HSEQ forward model. The horizontal dashed line is the height above which the a priori temperature profile makes a significant contribution to the retrieval.

control the retrieval and the model falls back to the a priori state. Figure 21 shows the retrieved temperature profile. The temperature profile retrieved from the twochannel measurements agrees to within the 1σ level with the temperature profiles calculated from the HC method. The profile simultaneously retrieved from both channels also agrees well with the HC retrieval for the LLRd channel. The retrieved backgrounds are within 0.21% and 1.9% of the a priori’s estimated from the measurements, and the retrieved dead time is 3.90 ns. The uncertainty budget is also similar above 40 km to the single-channel case (Fig. 22). The uncertainties below 37.5 km, where the LLRd makes its largest contribution are consistent with the 120

120

Low Digital Channel

120

110

110

100

100

90

90

Low Digital High Digital

High Digital Channel

110 100

80 70 60

Altitude (m)

Altitude (km)

Altitude (km)

90

80 70

80 70 60

60

50 50

50

40

40

30

30

40 30 -15 -1

-0.5

0

Jacobian / number density (x 10

-19

-6

-4

-2

) Jacobian / number density (x 10

-19

-10

-5

0

5

10

15

Percent Difference

0

)

Fig. 18. Temperature Jacobians, normalized to the number density for the LLR and HLR channels using the HSEQ forward model on 24 May 2012.

Fig. 20. Residuals between the two-channel photocount measurements with the HSEQ forward model. The LLR measurements contribute primarily below 37.5 km, the HLR measurements above 37.5 km. 10 March 2015 / Vol. 54, No. 8 / APPLIED OPTICS

1885

120 110 100 Retrieved A priori HC

80 70 60 50 40 30 150

200

250

300

350

Temperature (K)

Fig. 21. Two-channel HSEQ temperature retrieval from Purple Crow Lidar measurements on 24 May 2012 (red curve). The blue and green curve uses the HLR and LLR photocount profile, respectively, to retrieve temperature using the HC method. The a priori temperature profile (cyan curve) is the US Standard Atmosphere. The horizontal dotted line is the height above which the a priori temperature profile makes a significant contribution to the retrieval.

single-channel results, where the HLR channel was used down to 31 km. Since the retrieval now goes down to 25 km altitude, the uncertainty due to ozone density profile increases rapidly. In this situation, the additional ozone at the lower heights contributes 50

120 110

45

100

Altitude (km)

80 70 Statistical Smoothing Pressure Constant [O 3 ] O 3 Ray Exct Ray Exct

60 50 40

40

35

30

Ray

45

45

Ray(z)

MMM Gravity Total

30 0

10

20

Uncertainty (K)

30

25

0

0.2

0.4

Uncertainty (K)

Fig. 22. Uncertainty budget for the HSEQ forward model two-channel temperature retrieval on 24 May 2012. Left panel: statistical (blue), smoothing (red), pressure (yellow circle), lidar constant (magenta x), ozone density (green +), ozone cross section (blue *), density profile for Rayleigh extinction (red square), Rayleigh extinction cross section (blue diamond), variation of Rayleigh-scatter cross section with height (orange triangle), variation of mean molecular mass with height (yellow triangle), gravity (purple triangle), and total uncertainty (black line). The horizontal dashed line is the height above which the a priori temperature profile makes a significant contribution to the retrieval. Right panel: Expanded view of the uncertainties in the stratosphere. 1886

To illustrate the OEM’s robustness, Fig. 24 shows the comparison between the retrieval using the HSEQ forward model with two channels (LLRd and HLR) and the HC method on nine nights. For this comparison, the temperature profiles retrieved by the HC method are cut off at the height where the statistical uncertainty is above 10%. As the HC method temperatures are seeded at a signal-to-background level of 2, this choice ensures the top 10 to 15 km has been removed from the profile. The OEM retrieved temperatures are shown up to the greatest height given by the degrees of freedom of the retrieval on each night. The agreement is within the combined uncertainties, with each small bump in the HC method also occurring in the OEM retrieval.


Altitude (km)

Altitude (km)

90

D. Comparison of the OEM to the HC Method on Other Nights

Altitude (km)

Altitude (km)

90

to uncertainty in the stratosphere that exceeds the uncertainty in the knowledge of the ozone absorption cross section. The high signal-to-noise ratio measurements of the HLR channel at 37.5 km cause a jump downward in the statistical uncertainty along with increasing the systematic uncertainties due to the model parameters. This increase in uncertainty is due to the larger value of the temperature Jacobians at the lowest heights for the HLR channel compared to the LLRd channel Jacobians at the same height. The temperature retrieval above 37.5 km is predominantly due to the HLR channel measurements. To highlight the advantage of using the OEM instead of merging temperatures, consider the stratospheric temperature on two nights (Fig. 23). On 24 May 2012, the temperatures in the region of overlap between the two channels processed individually by the HC method have similar slopes, and it is straightforward to merge the temperatures. The OEM retrieval is also smooth in this region. Contrast this night to 7 August 2012, where the OEM retrieval is again smooth, but the temperature profiles retrieved by the HC method have different slopes and significant differences in magnitude below 40 km.

40

35

30

230

240

250

260

Temperature (K)

270

40

35

30

230

240

250

260

270

Temperature (K)

Fig. 23. Comparison of the OEM two-channel retrieval in the stratosphere (red curve) to single-channel profiles calculated from the HC method for the HLR and LLRd channels (blue and green curves) on 24 May 2012 and 07 August 2012, highlighting the differences in retrieved temperature profiles in the merging region. The cyan curve in both figures is the a priori temperature profile.

100 80 60 40 20 150 200 250 100 80 60 40 20 150 200 250

Temperature (K)

Temperature (K) 20120604

Temperature (K)

20120706

Temperature (K)

20120714

Temperature (K)

20120801

Temperature (K)

Altitude (km)

100 80 60 40 20 150 200 250

Altitude (km)

100 80 60 40 20 150 200 250

100 80 60 40 20 150 200 250

Altitude (km)

Altitude (km) Altitude (km)

100 80 60 40 20 150 200 250

20120526

20120525

100 80 60 40 20 150 200 250

Altitude (km)

Altitude (km) Altitude (km)

100 80 60 40 20 150 200 250

Altitude (km)

20120524

100 80 60 40 20 150 200 250

Temperature (K)

20120807

Temperature (K)

20120823

Temperature (K)

Fig. 24. Temperature retrieval (red curve) using two digital channels for nine nights of measurements. The blue curve shows the HC method retrieval for the same measurements using the HLR channel; the green curve is for the LLR channel. The temperatures from the HC method have been cut off at the height where the variance becomes greater than 10%. The shaded areas indicate the 1σ uncertainty including smoothing and statistical uncertainty for the OEM retrieval and the statistical uncertainties for the HC method.

The insensitivity of the retrieval to the choice of an a priori temperature profile is also demonstrated by these results, as the US Standard Atmosphere temperature profile was used as the a priori temperature profile for each night. The retrieved dead time was also consistent from night to night. The mean dead time was 4.02 0.066 ns. 6. Discussion

The comparisons in the previous sections show that both forward models, with and without assuming HSEQ, provide robust OEM retrievals, which yield temperature profiles, photocount background, and dead time. These retrievals compare favorably with expected values. The OEM retrieval also gives a full random and systematic uncertainty budget on a profile by profile basis. Furthermore, the OEM’s averaging kernels allow the determination of both the vertical resolution as a function of height and the maximum height to which the retrieval is mostly independent of the a priori profile. The HSEQ forward model is mathematically equivalent to the HC method, and this assertion is made clear by the nearly one-to-one correspondence in temperatures determined by both methods. The LE forward model requires that the true pressure profile be known to a high precision. The examples shown assumed the pressure profile is known to 1%, which results in systematic uncertainties due to pressure of about 4 K, with uncertainties due to knowledge of the lidar constant of the same order of magnitude. For more useful retrievals, it is recommended that the pressure profile be known through

the upper stratosphere, mesosphere, and lower thermosphere to less than 1%. This requirement may not be possible with the current generation of atmospheric instruments. In the past, individual investigators have addressed the question of data corrections, signalto-noise ratio, and seed pressure by attempting to correct their measurements before a temperature is retrieved. In an OEM framework, all significant instrumental and atmospheric corrections are built into the forward models, and level 0 measurements are used in the analysis. A temperature profile is retrieved along with a quantitative determination of the maximum retrieval height to which the retrieval is valid. The OEM determines both the vertical resolution of the retrieval as well as the full random and systematic budget for uncertainties. It is computationally efficient, performing the full nonlinear calculation in seconds on a laptop computer using a interpreted fourth-generation programming language. 7. Summary

This discussion of the application of OEMs to the retrieval of Rayleigh-scatter lidar temperatures can be summarized as follows. 1. The forward models presented for use with the OEM provide robust estimations of temperature. 2. OEM temperature retrievals using a forward model, which assumes HSEQ, are equivalent to the temperature profiles calculated using the HC method. 3. The OEM provides a cut-off height for retrievals. Below this height, the entire temperature profile can be used. 4. Vertical resolution is determined quantitatively at each height. 5. The OEM provides a complete uncertainty budget, including random uncertainties, smoothing error, and uncertainty due to model parameters including pressure, instrument parameters, and composition changes for each temperature profile retrieved. 6. If multiple channel measurements are available, a single temperature profile can be retrieved that is consistent with both measurements simultaneously. In this case, the dead time may be also retrieved in addition to the backgrounds and temperature profile. For the case of an analog and digital channel, no conversion of the analog signal into current or virtual photocounts is required. 7. OEMs are fast computationally and, thus, practical for routine data analysis. 8. Conclusions

OEM retrievals using the HSEQ forward model offer significant advantage over the HC and GS methods and are physically equivalent to traditional methods. The LE model allows the assumption of HSEQ to be dropped but requires high-quality estimates of the pressure as a function of height to produce useful 10 March 2015 / Vol. 54, No. 8 / APPLIED OPTICS

1887

results. Most if not all measurements of pressure available in the middle atmosphere at some point in their processing appeal to the assumption of HSEQ, in an atmosphere known to often not be in HSEQ. We encourage the atmospheric community to develop better methods for determining pressure profiles that do not depend on the assumption of HSEQ. When these pressure measurements are available, the LE forward model will be able to retrieve temperature without the assumption of HSEQ. The uncertainty budgets given in the study are heuristic and not meant to be definitive. The model parameter variances have been made as a “best guess” based on current best practice. The assessment and assignment of appropriate covariances to the model parameters should be a priority of the lidar measurement community. Much of this work is being performed by the International Space Science Institute’s Critical Assessment and Standardized Reporting of Vertical Filtering and Error Propagation in the Data Processing Algorithms of the NDACC Lidars working group (LeBlanc et al., http://www.issibern.ch/teams/ndacc/), but specific community recommendations for the OEM method on the uncertainties and best choices for model parameters would be beneficial. This study, as well as the recent application of OEM techniques to the retrieval of extinction and backscatter from an aerosol lidar by Povey et al. [7], establish OEM as an important new tool for retrieval of geophysical parameters from active remote sensing instruments. These studies are the first step in a new growth area, which can be called “measurement assimilation.” In the last two decades, dedicated effort has been made by meteorological services in the area of data assimilation, wherein forecast model predictions are improved by incorporating all available measurement information into the model to improve the forecast. In measurement assimilation, consider a station with numerous individual instruments, many of which individually determine or require the same atmospheric parameter. For instance, temperature information may come individually from microwave sounders, radiosondes, lidar, etc. The goal of measurement assimilation is to develop retrieval systems that analyze groups of measurements simultaneously, in a framework such as is described in this study. Forward models could be developed for various instruments, and joint parameters such as temperature could be retrieved. Parameters retrieved in this manner would be the best estimates based on all available information rather than a single instrument, and the weighting of individual instruments could change as a function of time, meteorological conditions, etc. The benefit for the user would be the ability to obtain on a continuous basis the best available measurement of a parameter from a station as opposed to multiple estimates at varying spatialtemporal resolutions. 1888


This work greatly benefited from scientific discussions with our colleagues. P. S. Argall, B. Calpini, G. Martucci, D. Ruffieux, and V. Simeonov. We especially thank T. LeBlanc for encouraging us to explore improvements to the quality of lidar temperature retrievals through his International Space Sciences Institute lidar algorithms initiative. We would like to thank the Federal Office of Meteorology and Climatology, MétéoSuisse, for its support of this project. R. J. S. would like to further thank the École Polytechnique Fédérale de Lausanne and The University of Western Ontario for supporting his research as well as the National Science and Engineering Research Council of Canada. The MLS measurements used in this study were acquired as part of the activities of NASA’s Science Mission Directorate and are archived and distributed by the Goddard Earth Sciences (GES) Data and Information Services Center (DISC). We also acknowledge the Aura MLS mission scientists and associated NASA personnel for the production of the data used in this study. References 1. A. Hauchecorne, M.-L. Chanin, and P. Keckhut, “Climatology and trends of the middle atmospheric temperature (33–87 km) as seen by Rayleigh lidar over the south of France,” J. Geophys. Res. Atmos. 96, 15297–15309 (1991). 2. P. S. Argall and R. J. Sica, “A comparison of Rayleigh and sodium lidar temperature climatologies,” Annales Geophysicae 25, 27–35 (2007). 3. R. T. Clancy, D. W. Rusch, and M. T. Callan, “Temperature minima in the average thermal structure of the middle mesosphere (70–80 km) from analysis of 40 to 92 km SME global temperature profiles,” J. Geophys. Res. Atmos. 99, 19001– 19020 (1994). 4. A. Hauchecorne and M. Chanin, “Density and temperature profiles obtained by lidar between 35 and 70 km,” Geophys. Res. Lett. 7, 565–568 (1980). 5. C. D. Rodgers, Inverse Methods for Atmospheric Sounding: Theory and Practice, Theory and Practice (World Scientific, 2011). 6. J. Khanna, J. Bandoro, R. J. Sica, and C. T. McElroy, “New technique for retrieval of atmospheric temperature profiles from Rayleigh-scatter lidar measurements using nonlinear inversion,” Appl. Opt. 51, 7945–7952 (2012). 7. A. C. Povey, R. G. Grainger, D. M. Peters, and J. L. Agnew, “Retrieval of aerosol backscatter, extinction, and lidar ratio from Raman lidar with optimal estimation,” Atmos. Meas. Tech. 7, 757–776 (2014). 8. D. P. Donovan, J. A. Whiteway, and A. I. Carswell, “Correction for nonlinear photon-counting effects in lidar systems,” Appl. Opt. 32, 6742–6753 (1993). 9. 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). 10. R. Sica and M. Thorsley, “Measurements of superadiabatic lapse rates in the middle atmosphere,” Geophys. Res. Lett. 23, 2797–2800 (1996). 11. R. J. Sica and P. S. Argall, “Seasonal and nightly variations of gravity-wave energy density in the middle atmosphere measured by the Purple Crow Lidar,” Ann. Geophys. 25, 2139–2145 (2007). 12. R. J. Sica, Z. A. Zylawy, and P. S. Argall, “Ozone corrections for Rayleigh-scatter temperature determinations in the middle atmosphere,” J. Atmos. Ocean. Technol. 18, 1223–1228 (2001). 13. R. D. McPeters, G. J. Labow, and J. A. Logan, “Ozone climatological profiles for satellite retrieval algorithms,” J. Geophys. Res. Atmos. 112, D05308 (2007).

14. M. Griggs, “Absorption coefficients of ozone in the ultraviolet and visible regions,” J. Chem. Phys. 49, 857–859 (1968). 15. M. Nicolet, “On the molecular scattering in the terrestrial atmosphere: an empirical formula for its calculation in the homosphere,” Planet. Space Sci. 32, 1467–1468 (1984). 16. P. S. Argall, “Upper altitude limit for Rayleigh lidar,” Annales Geophysicae 25, 19–25 (2007). 17. Committee on Extension to the Standard Atmosphere, “US standard atmosphere, 1976,” Tech. Rep. NOAA-S/T 76-1562 (National Oceanic and Atmospheric Administration, 1976). 18. M. Alpers, R. Eixmann, C. Fricke-Begemann, M. Gerding, and J. Hoffner, “Temperature lidar measurements from 1 to 105 km altitude using resonance, Rayleigh, and rotational Raman scattering,” Atmos. Chem. Phys. 4, 793–800 (2004). 19. M. Gerding, J. Höffner, J. Lautenbach, M. Rauthe, and F. J. Lübken, “Seasonal variation of nocturnal temperatures between 1 and 105 km altitude at 54° N observed by lidar,” Atmos. Chem. Phys. 8, 7465–7482 (2008). 20. P. Keckhut, A. Hauchecorne, and M. L. Chanin, “A critical review of the database acquired for the long-term surveillance of the middle atmosphere by the French Rayleigh lidars,” J. Atmos. Ocean. Technol. 10, 850–867 (1993).

21. P. Eriksson, C. Jiménez, and S. A. Buehler, “Qpack, a general tool for instrument simulation and retrieval work,” J. Quant. Spectrosc. Radiat. Transfer 91, 47–64 (2005). 22. S. A. Buehler, P. Eriksson, T. Kuhn, A. von Engeln, and C. Verdes, “ARTS, the atmospheric radiative transfer simulator,” J. Quant. Spectrosc. Radiat. Transfer 91, 65–93 (2005). 23. A. E. Hedin, “Extension of the MSIS thermosphere model into the middle and lower atmosphere,” J. Geophys. Res. Atmos. 96, 1159–1172 (1991). 24. J. J. Barnett and M. Corney, “Middle atmosphere reference model from satellite data,” in International Council of Scientific Unions Middle Atmosphere Program, Handbook for MAP (SCOSTEP, 1985), Vol. 16, pp. 47–85. 25. C. S. Gardner, D. C. Senft, T. J. Beatty, R. E. Bills, and C. Hostetler, “Rayleigh and sodium lidar techniques for measuring middle atmosphere density, temperature and wind perturbations and their spectra,” in World Ionosphere/ Thermosphere Study Handbook (SCOSTEP, 1989), Vol. 2, pp. 1–40. 26. M. J. Schwartz, A. Lambert, G. L. Manney, W. G. Read, N. J. Livesey, C. D. Boone, W. H. Daffer, M. G. Mlynczak, S. Pawson, J. M. Russell, M. L. Santee, K. A. Walker, and D. L. Wu, “Validation of the aura microwave limb sounder temperature and geopotential height measurements,” J. Geophys. Res. 113, 1–23 (2008).

10 March 2015 / Vol. 54, No. 8 / APPLIED OPTICS

1889

Stratospheric temperature measurement with scanning Fabry-Perot interferometer for wind retrieval from mobile Rayleigh Doppler lidar.

Searching for an optimal AUC estimation method: a never-ending task?

Motion field estimation for a dynamic scene using a 3D LiDAR.

Method to describe stochastic dynamics using an optimal coordinate.

Assessing and correcting topographic effects on forest canopy height retrieval using airborne LiDAR data.

Estimation of aircraft wake vortex parameters from data measured with a 1.5-μm coherent Doppler lidar.

The optimal temperature for 24-hour heart preservation using continuous coronary perfusion method.

Estimation of biomass burning influence on air pollution around Beijing from an aerosol retrieval model.

Photon-counting chirped amplitude modulation lidar using a smart premixing method.

Successful Retrieval of Uncoiled Coronary Guidewire Using Simple Balloon Method.

Estimation of left ventricular volumes using a simplified method.

A Low-Complexity Method for Two-Dimensional Direction-of-Arrival Estimation Using an L-Shaped Array.

Estimation of the refractive index structure characteristic of air from coherent Doppler wind lidar data.

Fat water decomposition using globally optimal surface estimation (GOOSE) algorithm.

Impact of data model and point density on aboveground forest biomass estimation from airborne LiDAR.

First measurements of a carbon dioxide plume from an industrial source using a ground based mobile differential absorption lidar.

Temperature retrieval from Rayleigh-Brillouin scattering profiles measured in air.

A method for optimal depth of the nasopharyngeal temperature probe: the philtrum to tragus distance.

A Robust Gradient Based Method for Building Extraction from LiDAR and Photogrammetric Imagery.

A modified Gustafson method of age estimation from teeth.

Storage and retrieval of THz-bandwidth single photons using a room-temperature diamond quantum memory.

Cancer incidence estimation method: an Apulian experience.

An improved method of fetal weight estimation using ultrasound measurements of fetal abdominal circumference.

Simulating acoustic scattering from atmospheric temperature fluctuations using a k-space method.