Smoothing radio occultation bending angles above 40 km

The 'statistically optimal' approach to smoothing bending angles derived from radio occultation (RO) mea- surements is outlined. This combines a measured bending angle profile with an a priori or background estimate derived from climatology, in order to obtain the most probable bend- ing angle profile. However, the method is only optimal if the error statistics of both the measured and background profiles are known and applied accurately. In this work it is shown that correlations in the background estimate have a signifi- cant role in determining the degree of smoothing in the solu- tion. We find that smooth profiles, consistent with the mea- sured values, can be derived if the correlations are approxi- mated analytically with a Gaussian, assuming a scale length of 6km. In regions where the observed and background error levels are comparable, the solutions take the general shape from the background estimate, centred on the observation data. The effects of correlated observation errors are also considered. It is shown that the quality of the temperature retrievals can be significantly affected by the choice of cli- matology used for background estimate.


Introduction
Radio occultation (RO) measurements of the earth's atmosphere using the GPS satellite constellation are a relatively new source of meteorological data (Kursinski et al., 1996;Rocken et al., 1997), with potential applications in both operational numerical weather prediction (NWP) (for example, the validation of NWP models in the stratosphere) and climate research (Leroy, 1997). The technique, which has been widely used used in the study of planetary atmospheres (e. g., Fjeldbo et al., 1971), is based on measuring how radio waves are bent by refractive index gradients in an atmosphere. It Correspondence to: S. B. Healy (sbhealy@meto.gov.uk) can be shown (for example, Kursinski et al., 1997) that this information can be inverted with an Abel transform to give a vertical profile of refractive index and, subsequently, temperature.
A detailed description of the inversion of RO measurements has been given by a number of authors (for example, Kursinski et al., 1997) and it has been recognised that the temperature retrievals at heights above 30 km are sensitive to residual ionospheric noise on the measured bending angles. Calculations with a bending angle profile provided by the Jet Propulsion Laboratory (JPL), smoothed to a sampling rate consistent with the time required for a ray to traverse the local Fresnel diameter at the limb (Kursinski et al., 1997), indicate that an error in a single bending angle value of 10 µrads near 50 km can lead to temperature errors of ∼1.5 K at 30 km. Consequently, well known 'statistical optimization' techniques have been employed in an attempt to smooth the noise and improve the accuracy of the retrievals, as originally suggested by Sokolovskiy and Hunt (1996). In principle, these methods combine the observed bending angle profile with a smooth background estimate, derived from climatology in a statistically optimal way, to obtain the most probable bending angle profile. However, it is only optimal if the error statistics (magnitude and correlations) of both the background and observations are known and applied accurately. To date (Hocke, 1997;Gorbunov and Gurvich, 1998;Steiner et al., 1999), this approach has been employed by making the additional assumption that both the background and the observation errors are not vertically correlated. While this may be a reasonable approximation for the latter, it will not be the case for the background estimate derived from climatology. Furthermore, it will be demonstrated that the correlations in background errors can have a significant impact on both the degree of smoothing in the optimal bending angle profile and in the retrieved temperatures.
The purpose of this paper is to illustrate the limitations of the 'statistically optimal' approach when correlations in the background estimate are neglected. It is demonstrated that smooth profiles consistent with the observation data can be derived by introducing a simple Gaussian form to approximate the error correlations. In Sect. 2, the processing of RO data is briefly outlined and a description of the statistically optimal smoothing technique is given in Sect. 3. Results are presented in Sect. 4, followed by discussion and conclusions in Sects. 5 and 6 respectively.

Inversion of RO measurements
The theory of RO measurements, using GPS satellites, has been described by Kursinski et al. (1997); Rocken et al. (1997) and is only outlined here. The technique is based on measuring how radio waves transmitted by a GPS satellite are bent by refractive index gradients before being received at a low earth orbiting (LEO) satellite. If the atmosphere is spherically symmetric, the ray path will be determined by Bouguer's formula (Born and Wolf, 1986), where n is the refractive index, r is the radius value, φ is the angle between the ray path and the local radius vector and a is a constant for a ray path, known as the impact parameter. The bending angle, α, is then, where x = nr. The GPS satellite transmits at two frequencies (L1=1575.42 MHz and L2=1227.6 MHz) and the bending due to ionospheric plasma is removed or 'corrected', to a first order, by taking a linear combination of the L1 and L2 bending angle values (Vorob'ev and Krasil'nikova, 1994). For a spherically symmetric atmosphere, the variation of the corrected bending angle with impact parameter can be inverted with an Abel transform (Fjeldbo et al., 1971) to recover the refractive index profile, The refractive index is related to the total pressure, temperature and water vapour pressure P , T and P w through, where c 1 (= 7.76 × 10 −5 K/hPa) and c 2 (= 3.73 × 10 −1 K 2 / hPa) are known constants (Bean and Dutton, 1968). This can be rewritten as n = 1 + 10 −6 N , where N is the refractivity. For a dry atmosphere (P w = 0), the refractivity profile can be used to integrate the hydrostatic equation dP /dz = −ρg and a vertical temperature profile is derived from the ideal gas law, since P = ρRT = 10 −6 N T /c 1 . The calculation is initialised with an estimate for the temperature at the upper level boundary. This is typically 260 K at 50 km (Kursinski et al., 1997).

Statistically optimal smoothing
A number of authors (Rocken et al., 1997;Hocke, 1997;Gorbunov and Gurvich, 1998;Steiner et al., 1999) have applied a statistically optimal retrieval technique (also known as statistical regularization) to smooth the bending angle profiles above 40 km. The technique has a clear theoretical basis (Rodgers, 1976) and is commonly used in NWP, in the operational processing of satellite sounding data (e.g. Eyre et al., 1993) and, more generally, when solving the 'data assimilation problem' (Lorenc, 1986). In current context, the aim of this approach is to obtain the most probable bending angle profile,α, given the noisy observation, α o , and a smooth, a priori, or background estimate, α b , which is derived from climatology. Given unbiased errors with Gaussian probability distributions, it can be shown that the optimal solution minimizes the quadratic cost function J , where B and O are the background and observation error covariance matrices respectively. If the problem is linear, the solution that minimizes the cost function can be written as, This can also be solved in an alternative form using the ma- The method is a generalised least squares approach, which incorporates a priori information. The solution,α, represents the optimal, simultaneous fit, to within the expected errors of both the background and observation vectors. The errors are statistically characterised with the covariance matrices B and O. The diagonal elements of these matrices contain the expected error variances of the individual bending angle values; the off-diagonal elements specify the degree to which these errors are correlated with each other. The background error matrix in bending angle is theoretically given by mapping the climatological covariance matrix C of, for example, temperature on fixed pressure levels (commonly used in meteorology) into bending angle space with a linear matrix transformation B = KCK T . The matrix K contains the partial derivatives of the simulated bending angles with respect to the temperature values (i.e., K ij = ∂α i /∂T j ). In physical terms, this transformation is simply mapping how statistical errors in the climatology produce equivalent statistical errors in the simulated bending angle values. Observation errors can be estimated from simulation (Kursinski et al., 1997) and monitoring the measured values. It should be emphasised that the method requires both the background and observation errors (ie, the values minus the true bending angle values) must be statistically bias-free to ensure an unbiased solution. If either α b or α o are found to be biased, a bias correction must be included prior to the smoothing step. This is often necessary in the operational processing of satellite sounding products for NWP, as discussed by Eyre et al. (1993).
In this smoothing problem we generally assume that the background and observation vectors have very different error characteristics. The former will be large and correlated, as a result of correlations in the climatology errors and transforming into bending angle space. Therefore, the overall shape of α b should be broadly correct but it can be shifted, or offset, relative to the true profile. Consequently, on performing an eigenvector/value decomposition of B most of the uncertainty (the largest eigenvalues) in the background estimate will be found in the slowly varying, low frequency eigenmodes. The observation errors should be smaller than the background error throughout most of the atmosphere and they are assumed to be uncorrelated or high frequency (however, note that correlated observation errors are considered in Sect. 4.2). In physical terms, Eq. 6 represents a statistically optimal filter which damps the high frequency noise. When the background and observation errors are of comparable magnitude, the solution,α, will tend to take the general shape of bending angle profile from the background α b , but centred on observation data α o . See Rodgers (1976, Eqs. 46-55) for a more detailed discussion.
A significant difficulty with the 'statistically optimal' approach is deriving an accurate covariance matrix B. As noted, theoretically this requires mapping a climatological covariance matrix, C, (which should ideally extend up to around 100 km) into bending angle space but it is not always possible to find a suitable matrix C and it may be necessary to derive an approximate B. Given reasonable estimates for the variance values, σ 2 b (i), and a scale length, l, 'artificial' covariance terms can be derived using (Rodgers, 1990), where a i and a j are the ith and j th impact parameter values.
(Note that it is important to ensure that an 'artificial' covariance matrix constructed with this method is positive definite). A simpler approach, common in the processing of RO data (Hocke, 1997;Gorbunov and Gurvich, 1998;Steiner et al., 1999), is to assume that the background errors are uncorrelated (l = 0). The optimal solution then reduces to, for the ith value, where σ 2 b (i) and σ 2 o (i) are the background and observation variance values. An alternative expression has also been used (Hocke, 1997;Steiner et al., 1999), but formally this will not produce the most probable solution, even if the background and observation errors are uncorrelated and the standard deviation values are known accurately. Comparing Eqs. 8 and 9, it can be shown that Eq. 9 will give more weight to the measurements if σ o (i) > σ b (i) and, conversely, more weight to the background if σ o (i) < σ b (i).
In addition, Hocke approximates Consequently, this form of smoothing is not usually applied below 40 km because the solution can be biased towards the background. It should be emphasised that the assumption of uncorrelated background errors conveniently simplifies the calculation but it is questionable since any climatological estimate will have highly correlated errors. Furthermore, although Eqs. 8 and 9 produce (different) solutions that are smoother than the measurement vector, the degree of smoothing is limited because each bending angle value,α(i), is bounded by α o (i) and α b (i), the measured and background values respectively. Since the background profile is smooth by definition, a solution that is constrained to be between the observed and background profiles will always retain a degree of high frequency noise.

Comparison of smoothing methods
In this work the full vector form (Eq. 6) of the statistically optimal smoothing has been employed. The numerical solution of Eq. 6 is straightforward and can be performed efficiently using Cholesky decomposition (Press et al. 1993). In general, the observation errors are assumed uncorrelated (correlated errors are considered in 4.2.) and have a constant standard deviation value of σ o = 5 µrads. This value is consistent with the characteristic noise observed in the GPS/MET corrected bending angles around 50 km, provided by JPL. The background bending angle profile, α b , has been calculated using Eq. 2 for a refractivity profile derived from the MSIS climatology (Hedin, 1991). We are not aware of an error covariance matrix C for the MSIS model from which to calculate B. Therefore, following Hocke (1997), the standard deviations (σ b (i)) of the background bending angle values are estimated as σ b (i) = 0.2 × α b (i). The error values are in reasonable agreement with a background error covariance matrix (B = KCK T ) derived from a global climatology (C) of 1200 radiosonde and rocketsonde profiles obtained from NOAA/NESDIS. For example, the percentage errors (ie, the standard deviation values derived from the B matrix) vary from ∼ 10% at 40 km to ∼ 20% at 60 km. Since the MSIS climatology is a more sophisticated model, varying both with location and season, we believe the assumed errors are conservatively high estimates. The error covariances above 30 km are assumed to have a scale length of l = 6 × 10 3 m. Note that the covariances estimated from B = KCK T suggest a significantly larger scale length (exceeding 20 km), but we believe this primarily arises as a result of using a single global climatology for the calculation. Below 30 km the errors are assumed uncorrelated for numerical noise considerations (the B matrix is found to have negative eigenvalues with a single precision calculation) and to ensure that the retrieved profiles tends towards the measurement below this height. Since the magnitude of the observation error is assumed constant, but the background error is specified as a fraction of the simulated bending angle value, the smoothed profile tends towards the observation at lower altitudes. Typically at 30 km σ b /σ o ∼ 12. Conversely, the observation error exceeds the background error estimate (σ b < σ o ) for tangent heights above ∼45 km. This is reasonable since we are required to assume an a priori temperature value at about 50 km in order to retrieve a temperature profile. Therefore, the method produces a smooth transition between the observed and background profiles over altitudes where the error levels are comparable. Figures 1-3 are examples of the statistically optimal smoothing above 40 km for GPS/MET profiles provided by JPL. The smoothed profiles derived, assuming uncorrelated errors and using Eqs. 8 and 9 (the latter is employed assuming Hocke's estimates for the observation errors), are also shown for comparison. Note that the background profiles calculated from the MSIS climatology provide a very good first guess in each case. When the smoothed profile is evaluated with Eq. 6 it is clear that the noise level is reduced and the solution remains in good agreement, in a least squared sense, with the observations. Note that the retrieved profiles retain the structure below 30 km. In most cases the profiles derived, assuming uncorrelated errors, are noticeably noisier between 30-50 km. By definition (Eqs. 8 and 9), the smoothed bending angles are always between the measured and background values when the errors are uncorrelated, limiting the degree of smoothing. However, Eq. 9 gives less weight to the background above 45 km, where σ o (i) > σ b (i), because the ratio of standard deviations is used rather than variances.
The temperature retrievals that correspond to the bending profiles are shown in Figs. 4-6. In each calculation, the temperature is initialised with the MSIS value at 50 km. In Figs. 4 and 6, the solutions derived with the correlated background errors produce the smooth retrievals which remain a good fit to the 'measured' (i.e. unsmoothed) temperature values. The retrievals derived assuming Eqs. 8 and 9 are reasonable but, as expected, they clearly retain some of the high frequency noise. Figure 5 is a case where the smoothing has a significant effect on the temperature retrievals. The smoothed solutions are broadly in agreement, all clearly correcting a gross error in the measurement, but Fig. 5a is the smoothest.

Sensitivity to assumed errors and correlations
The sensitivity of the smoothed temperature retrievals to the assumed background errors and correlation length is shown in Figs. 7 and 8, for the '95-07-01-02:44met-gps15' measurement used in Figs. 2 and 5. We have chosen this measurement because it is the most sensitive to the assumed error levels (the largest temperature differences calculated for the two other profiles were less than 2 K). This may be related to the fact that the corrected bending angle profile is only reaches ∼52 km above the surface in this case, compared with ∼60 km for the other profiles. As a consequence, the solution is less well constrained. Figure 7 shows that the smoothed temperatures can differ by up to ∼12 K near 45 km when the fractional background errors are assumed to be 0.1 Fig. 1. Illustrating the smoothing of corrected bending angles for JPL profile '95-05-05-01:33met-gps25' using (a) Eq. 6, (b) Eq. 8 and (c) Eq. 9. The vertical axis is 'impact parameter minus radius of curvature'. and 0.3. This is a large difference, but it should be put into the context that the unsmoothed profile is 40 K cooler in this region. The general trends are as expected. As the background error is increased, the solution moves towards the 'measured' profile. At lower altitudes the solutions converge. For example, at 30 km the temperature values agree to within 2.6 K, whereas the 'measured' value is 10 K cooler.
In Fig. 8, as expected, reducing the background correlation length gradually introduces high frequency structure into the retrieval and, in the limit of l = 0, the retrieval will converge to Fig. 5b. Varying the correlation length leads to temperature differences as high as 20 K around 45 km, but once again the solutions converge to within 2.8K at 30 km. In Sect. 3 it was noted that when the observation errors are assumed to be uncorrelated, the optimal solution takes the general shape of bending angle profile from the background Fig. 3. Illustrating the smoothing of corrected bending angles for JPL profile '95-07-01-06:22met-gps16' using (a) Eq. 6, (b) Eq. 8 and (c) Eq. 9. The vertical axis is 'impact parameter minus radius of curvature'. α b but centred on observation data α o . Introducing correlated observation errors leads to a slightly more complicated picture, because the retrieval gives more weight to high frequency structure of the observation and less weight is given to the low frequencies. An advantage of solving the full vector form (Eq. 6) of the statistically optimal smoothing problem is that it is relatively simple to investigate the importance of vertically correlated observation errors. It is likely that the Fig. 4. The retrieved temperature profiles corresponding to the smoothed bending angle profiles given in Fig. 1. processing of the observations from phase measurements to corrected bending angle and impact parameter will introduce some degree of correlation. Possible sources of correlated errors include satellite position and velocity errors, clock error and incomplete calibration of the ionospheric bending (Kursinski et al. 1997); but further work is required to determine the typical magnitude for the vertical scale length. In this work we have investigated the sensitivity of the tem- perature retrievals to correlations in the observation errors by assuming that off diagonal covariance terms can be expressed in a form analogous to Eq. 7, where h is the assumed observation error scale length. The value used here is h = 3 km. We have found that the '95-05-05-01:33met-gps25' and the '95-07-01-06:22met-gps16' smoothed temperature profiles differ by only ∼ 0.6 K at 30 km when the correlations in the observations are introduced. More generally, between 30-50 km the temperature differences are less than 2.3 K. Once again the '95-07-01-02:44met-gps15' result is the most sensitive to the change in the assumed errors, as illustrated in Fig. 9. The temperature values at 30 km differ by ∼ 5.2 K and the largest difference between 30-50 km is 23 K. This illustrates that some profiles  can be sensitive to the assumed error levels and correlations. It suggests that it may be useful to perform the smoothing with a range of assumed error levels and correlations in order to gauge the sensitivity and robustness of the results.

Sensitivity to assumed climatology
The results outlined above have been repeated using the US Standard Atmosphere (1976) to derive the background bending angle profiles. Clearly this is a much cruder approach as it does not account for latitude or seasonable variations, which are included in the MSIS model (Hedin, 1991). We have found that significantly poorer results are obtained when using the standard atmosphere. Figures 10 and 11 show a comparison of the smoothed bending angle profiles and retrieved temperatures calculated with the MSIS climatology and US standard atmosphere for the '95-07-01-06:22met-gps16' measurement which is the worst case. In both calculations, the background covariance matrix is evaluated assuming that the errors are given by 0.2 × α b and the correlation length l = 6 km. The observations errors are assumed to be uncorrelated. In Fig. 9, both solutions are smooth and represent a reasonable fit to the measurements. However, the background values evaluated from the standard atmosphere are systematically higher for this particular measurement and the profile converges to significantly larger background values between 55-65 km, where the measurement error exceeds the background (σ o /σ b > 10). This results in larger refractive index values and the increased temperatures shown in Fig. 10. The retrieval based on the MSIS climatology is clearly a superior fit of the measurement. However, it could justifiably be argued that the background error levels for the standard atmosphere case should be inflated. Increasing the error level to 0.5 × α b and broadening the correlation length to 12 km removes some of the warm bias but it is not possible to produce as good a retrieval with the standard atmosphere for this measurement.

Discussion
The results presented in Sect. 4 illustrate that the temperature profiles derived from RO data are sensitive to any smoothing of the corrected bending angles and it has been demonstrated that introducing background error correlations with a simple model can produce significantly smoother solutions. We believe the approach adopted previously (Hocke, 1997;Steiner et al., (1999)) is problematic because Eq. 9 does not give the statistically optimal, most probable bending angle profile even if the background and observation error estimates are correct. In addition, in the present work we have attempted to use improved error estimates, although it is clear, given the importance of the assumed errors, that further work is required in this area. In general, the observation errors are as- sumed uncorrelated, with a magnitude of 5 µrads. This figure was chosen after analysing the noise structure above 40 km, apparent in GPS/MET corrected bending angle profiles provided by JPL. Clearly, it should be reduced when processing data from improved receivers with lower noise levels.
For example, the 'GRAS' instrument is expected have corrected bending errors of ∼ 1 µrads (GRAS-SAG, 1998) and, in general, future receivers should produce corrected bending angles with sub-micron noise levels. These should produce considerably better temperature retrievals in the stratosphere. The background error (standard deviation) values are estimated at 20% of the bending angle values forward modelled from the MSIS climatology and the correlations are derived assuming a Gaussian form, with a scale length estimate l = 6 km. We have found examples where temperature retrievals between 40-50 km are sensitive to variations in both of these values but, generally, the solutions converge at lower altitudes and agree to within 2K at 30 km. Similar results were found when correlated observation errors were considered, assuming a scale length of h = 3 km. Since the smoothing is not computationally expensive, it may be useful to process RO data with a range of background and observation error levels and correlation lengths in order to gauge the sensitivity of the retrievals and to estimate probable temperature errors. Clearly, the uncertainty in these parameters means that the solution profiles,α, will still not be statistically optimal in a formal sense and further 'tuning' is required. Whilst acknowledging limitations in the current error estimates, it has been demonstrated that the correlations in the background errors play an important role in producing a smooth solution. At altitudes where the smoothing modifies the profile, the solution takes the general shape from the background centred on the observation data. Hence, this approach produces profiles which are smoothed considerably but remain consistent with the measurement. The importance of the choice of climatology is notable as Fig. 11. The sensitivity of the temperature retrievals to the assumed climatology, using '95-07-01-06:22met-gps16' measurement.
it was not originally anticipated. This work was initially performed using the US Standard atmosphere and this appeared to give reasonable results. However, it is now clear that retrievals based upon the MSIS climatology (Hedin, 1991) are considerably better as a result of improved accuracy of the background estimates. In addition, we have found that the approach outlined by Hocke (1997) can produce discontinuities at 40 km (below which no smoothing is applied) in the smoothed bending angle profile and temperature retrievals when implemented using the standard atmosphere in this work. A possible limitation of the approach outlined in this work is that it may oversmooth some real, high frequency, wavelike stratospheric information contained in the measurements. For example, in Fig. 8 it is shown that high frequency structure above 30 km is re-introduced into the retrieved temperature profiles if the assumed background error correlation length is reduced to 3 km. Once again, this clearly illustrates the need for further work to obtain an improved estimate of the correlation lengths. In practice, the smoothing applied here has little effect on the wave-like structure below 30 km apparent in the measurement. This is because the background errors are assumed to be uncorrelated to ensure that the retrieved bending angle values tend strongly towards the measured values. Note that in many potential applications of the RO data, for example the validation of NWP models, it is preferable to sacrifice some resolution in order to improve the overall accuracy of the retrieval. Therefore, we believe this implementation of statistically optimal smoothing will prove useful in the processing of RO data, particularly for NWP users.

Conclusions
The previous implementations of smoothing, based on statistical optimization (Hocke, 1997;Steiner et al., (1999)), used in the processing of RO data have been extended to introduce correlations in the background bending angle errors. We have found smooth profiles, consistent with the measurements, can be derived assuming large background errors (σ = 0.2 × α b ) which are vertically correlated with a Gaussian form (scale length l = 6 km). Further 'tuning' of these values is desirable before this approach is used routinely but we have, nevertheless, demonstrated that the correlations have an important role in determining the degree of smoothing and should be considered carefully in any such approach. Correlations in the observation errors have been considered, assuming a vertical scale length of h = 3 km. It has been shown that in some cases the retrieved temperature profile can be sensitive to assumed observation error correlation length. We have found that retrievals based on the MSIS climatology (Hedin, 1991) are considerably better than calculations based on a globally averaged climatology.