The aim of this project is to understand the advantages and pitfalls of an optimal estimation inversion retrieval scheme as described by Rodgers (1976,2000). As said by Clive Rodgers himself:
The problem that will be examined here is as follows: given a measurement or series of measurements of thermal radiation emitted by the atmosphere, the intensity and spectral distribution of which depend on the state of the atmosphere in a known way, deduce the best estimate of the state of the atmosphere.
There are two distinct aspects to this problem which are not always clearly separated; they may be described as the 'inverse' problem and the 'estimation' problem. The inverse problem is the matter of inverting a known equation which expresses radiation as a function of the atmospheric state, so as to express atmopsheric state in terms of radiation. This is usually an 'ill-posed' problem; i.e., it has no mathematically unique solution.
We therefore have an estimation problem, that is, to find the appropriate criteria which determine the best solution from all the possible ones which are consistent which the observations.
A simple retrieval scheme will be build and used to investigate the effect of the different parameters that determine the outcome of the final retrieval.
A full treatment on optimal estimation is beyond the scope of this course, but students would benefit from reading or skimming Daniel Jacob's introductory notes on the subject.
Use Equations (7.38) and (7.39) from Remote Sensing of the Lower Atmosphere: An Introduction (1994) to calculate 10 equally spaced temperature weighting functions W between 6 and 24 km. Use an atmospheric scale height H (for pressure) of 7 km. These weighting functions, while idealized, capture the general behavior of weighting functions calculated more exactly for realistic gas absorption cross sections.
Calculate synthetic radiances y from the tropical temperature profile using the linear forward model y = Wx where x represents the true temperature profile. Do not add instrument noise (yet).
You will now attempt to retrieve the temperature profile from the synthetic radiances.
Plot the weighting functions, the true temperature profile, the a priori temperature profile, and the retrived temperature profile. Also plot the error in the temperature profile (retrieved-true) as well as the a posteriori error profile (the square root of the diagonal elements of the posterior error covariance matrix). How does the retrieved profile compare to the true profile? How does it compare to the a priori temperature profile?
Verify the results by:
a. Showing that the retrieved profile leads to the observed (synthetic)
radiances.
b. The Weighting Function matrix W tells us how the radiances
y respond to changes in the true temperature profile x. That is,
Δy = WΔx. Determine the changes in the radiances for
a change (e.g., 10 K) in the temperature profile at 0, 20, and 50 km (separately).
How is the a posteriori error profile related to the sensitivity of
the radiances to changes in the temperature profile?
(Hint: You do not need to re-run the retrieval to answer this.)
a. Next, add random noise (such that values are normally distributed about zero with a standard
deviation of 5 Ky) to your synthetic radiances, and re-run the retrieval. Be sure
to tell the retrieval program that the assumed error is now 5 Ky (this goes into vector "y_err"
and specifically represents the assumed standard deviation of the error of each channel), and discuss how (and at which
altitudes) the retrieval is most affected.
b. What happens if you tell the retrieval program that the assumed error in
y is still very low, like 0.25K, even though it isn't? That is, use the same noisy values for
the measured radiances y as in step 6a above, but change y_err back to the small value assumed
in step 4. Why do you think this happened?
Complete steps 2-6 once again (you may exclude steps 5 and 6b this time), this time using either the mid-latitude winter temperature profile or mid-latitude summer profile to compute y (your choice) instead of the tropical temperature profile. Use a different a priori profile as well (your choice as well on the type of profile to use).
Rodgers, C.D., Inverse Methods for Atmospheric Sounding: Theory and Practice. World Scientific, 2000.
Rodgers, C.D., Retrieval of atmospheric temperature and composition from remote measurements of thermal radiation. Reviews of Geophysics and Space Sciences, 14, 609 - 624, 1976.
Rodgers, C.D., Characterization and error analysis of profiles retrieved from remote sounding measurements. Journal of Geophysical Research, 95, 5587 - 5595, 1990.
Stephens, G.L., Remote Sensing of the Lower Atmopshere. Oxford University Press, 1994.
To make life a little easier it is assumed in the above that the Planck function is linear in temperature (as is the case for microwave observations) and that there is no surface contribution to the observed radiances. Therefore, the radiances can be calculated as a matrix multiplication between the weighting functions and the temperature profile: y = W x. In this case, the radiances are actually brightness temperatures, and will have units of Kelvin.
The files tro.dat, mlw.dat, and mls.dat contain the necessary temperature profiles for this analysis. For this project, altitude and temperature are the only variables needed. All temperature profiles should have 61 levels, from an altitude of 0 km to 60 km (this goes for the isothermal profile in addition to the provided profiles). The information in the datasets are based on McClatchey et al. (1972) and are organized according to altitude (Z), atmospheric pressure, temperature, air density, water vapor density, and ozone density.
If using IDL to do these calculations, be sure that the input vectors "y" and the like all (1,M) column vectors (where M is the number of elements). This is because IDL orders its matrices (# of columns, # of rows) which is backwards from Fortran and how we generally think of matrices.
