- Email: [email protected]
Fast monochromatic radiative transfer calculations for limb sounding
Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
www.elsevier.com/locate/jqsrt
Fast monochromatic radiative transfer calculations for limb sounding Anu Dudhia ∗ , Paul E. Morris, Robert J. Wells Atmospheric, Oceanic and Planetary Physics, University of Oxford, Parks Rd, Oxford OX1 3PU, UK Received 20 August 2001; accepted 29 November 2001
Abstract Satellite observations of atmospheric infrared spectra can be modeled accurately with line-by-line calculations, but these are too slow to be incorporated into operational retrieval schemes. However, a monochromatic calculation is still feasible if the line-by-line summation is replaced by pre-tabulated absorption coe7cients, requiring a three-way optimization of storage space, accuracy and access time. Such a scheme is used for the operational processing of data from MIPAS, a limb-viewing interferometer. The tabulated data are compressed to a manageable size using singular value decomposition, although the reconstruction adds a small overhead. The number of monochromatic radiative transfer calculations is reduced by determining suitable quadrature points in the spectral domain, which reduces both processing time and data storage requirements. An important aspect of such optimizations is the control of the associated errors in the forward model calculation. The result is an acceleration of monochromatic forward model calculations by one or two orders of magnitude compared to a line-by-line calculation without any signi:cant loss of accuracy, while obtaining data compression factors of 100 or more compared to a direct tabulation of absorption coe7cients. ? 2002 Elsevier Science Ltd. All rights reserved. Keywords: Radiative transfer; Absorption coe7cient look-up tables; Singular value decomposition; Fast forward model; Remote sensing; MIPAS
1. Introduction The inversion of infrared satellite measurements to obtain pro:les of atmospheric temperature and composition generally starts with an assumed pro:le. A ‘forward model’ is then used to calculate the expected measurements from such a pro:le, these are compared with the observed measurements and the pro:le adjusted iteratively until convergence is achieved. The forward model is usually the ∗
Corresponding author. Tel.: +44-1865-272922; fax: +44-1865-272923. E-mail address: [email protected] (A. Dudhia).
0022-4073/02/$ - see front matter ? 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 2 - 4 0 7 3 ( 0 1 ) 0 0 2 8 5 - 0
746
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
most time consuming part of the retrieval process and, while the precision of the retrieved products may be determined by the instrument noise, the accuracy is often limited by the approximations used in the forward model. To model a radiance measurement R requires calculations along the following lines (subscript denoting spectrally varying quantities and superscript i denoting molecule-dependent quantities): R = L d; (1) L = =
B
d ds; ds
i ;
(2) (3)
i
i
i i = exp − k ds ;
(4)
where L is the monochromatic radiance at wavenumber ; is the channel spectral response, B is the Planck function (local thermodynamic equilibrium assumed), the transmittance to the satellite from point s along the ray path, i the transmittance of species i with density i and ki the absorption coe7cient (we will drop the superscript i from now for simplicity). For limb-sounders, there is usually an additional convolution in the elevation-angle domain to represent the :nite :eld-of-view. The most accurate method of evaluating the absorption coe7cient k is a ‘line-by-line’ calculation [1– 4]: a summation of the contributions of all spectroscopic lines in the vicinity, each modeled using the appropriate local path conditions of pressure, temperature, and, occasionally, absorber density. However, such calculations are usually too time consuming to form part of a near realtime retrieval. To model satellite measurements in near real time, the usual approach is to use a band model [5 –8]. For these, the order of the spectral and path integrations in Eqs. (1) and (2) is reversed, so that the monochromatic transmittance can be replaced by a spectrally averaged value , J either computed for the entire channel response or smaller intervals. These band-transmittances are pre-computed for a variety of path conditions, leaving just the path integration in the forward model. While this is extremely fast, band transmittances do not follow Beer’s Law (Eq. (4)) so further assumptions have to be made when obtaining the transmittance increments d =ds J and handling multiple absorbers (J = i Ji ), and this fundamentally limits the accuracy of this approach. Given current computing speeds, a reasonable compromise might be to use the monochromatic radiative transfer (Eqs. (1–3)) but, rather than evaluate k using a line-by-line model, use look-up tables containing values which have been pre-computed for a range of path conditions (e.g. [9]). The main drawback to this direct tabulation approach is that, in order not to introduce signi:cant interpolation errors, the absorption coe7cient tables tend to be rather large and if these cannot be contained within the computer memory signi:cant time can spent accessing the data from disk. Various techniques for compressing the look-up tables have been suggested, such as parameterizing the temperature dependence by a polynomial, or using singular value decomposition (SVD) [10]. These require additional CPU time in order to reconstruct the absorption coe7cients, but this may be outweighed by the bene:ts of reducing the data to a size which can be contained within the computer memory.
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
747
So far, such techniques have only been considered for nadir-viewing instruments. Here we present the approach adopted for a limb-sounding interferometer: the Michelson interferometer for passive atmospheric sounding (MIPAS) due for launch on the Envisat satellite in 2001. The Level 2 operational processing will retrieve pro:les of temperature and composition in near real time [11] with a monochromatic forward model using SVD-compressed look-up tables (LUTs) of absorption coe7cients and further optimized by using quadrature points (‘irregular grids’) for the spectral integration. This paper describes the generation of such tables, the associated optimizations, and discusses various alternatives. 2. CPU considerations MIPAS will acquire spectra covering 685 –2410 cm−1 at 0:025 cm−1 spacing every 4:6 s. A typical limb scan sequence from 68–6 km comprises 17 such spectra and from these data pro:les of temperature and six species are retrieved (although not for all tangent altitudes). Only small subsets (‘microwindows’) of the spectra will be used in the retrievals [12,13], totalling around 10 cm−1 per spectrum per retrieved species (Table 1). With 40 measurements=cm−1 , this means that the near real time processing absorbs ∼500 measurements every second. The spectral integration (Eq. (1)) behaves monochromatically if it is performed at a su7ciently :ne spacing to resolve atmospheric lines. A resolution of 0:0005 cm−1 is usually considered adequate (since Doppler broadening limits lines to a minimum half-width of around 0:001 cm−1 ), giving 50 radiance calculations for each MIPAS measurement on a 0:025 cm−1 grid. By simultaneously modeling multiple tangent paths for the entire limb-scan, a high degree of replication is possible so that each additional spectrum, including FOV convolution, only requires about 10 additional calculations of k for the whole atmospheric path. So to compute transmittances for one gas for one MIPAS measurement requires about 50 × 10 = 500 computations of k . Allowing for 10 absorbing species per microwindow and four retrieval iterations this gives a :gure of 107 k computations every second. To compute the absorption coe7cient with a line-by-line model using some version of the HumlPQcR ek algorithm [14] requires 10 –100 Soating-point operations (FLOPs) per line so, considering 10 –100 local lines, typically 103 FLOPs per absorption coe7cient. Thus, to retrieve MIPAS data in near real time using a line-by-line forward model would require a 10 Giga-FLOP processor. By using look-up tables, where the absorption coe7cients can be obtained in ¡100 operations, Table 1 Summary of the microwindows selected for operational processing of MIPAS data (more details in Appendix B). Columns show the number of microwindows per retrieval, the average number of absorbers per microwindow and the average microwindow width (cm−1 )
Retrieval
NMW
NJ gas
J T
Retrieval
NMW
NJ gas
J T
pT CH4 H2 O HNO3
10 5 6 5
6.0 11.0 7.7 9.6
1.2 2.2 2.0 3.0
N2 O NO2 O3
6 5 3
12.2 5.6 12.0
2.0 2.6 3.0
748
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
this is reduced to ¡1 Giga-FLOPs. In practice, the radiative transfer calculations become the limiting term rather than the k computation, but both can be further reduced by a factor 3–10 if the integration is performed explicitly on a subset of quadrature points from the 0:0005 cm−1 grid. 3. Accuracy criteria The representation of absorption coe7cient by means of look-up tables is a three-way compromise between storage space, CPU time and accuracy, so some accuracy requirement also has to be speci:ed. It is simple to express the accuracy in terms of the diUerence between the interpolated absorption coe7cient and a ‘true’ value calculated with a line-by-line method. However, it is the accuracy of the forward model calculation that is of interest, and this can only be determined by comparing radiative transfer calculations using both a line-by-line model and look-up tables. For this work, the RFM [4] was used since this could use both line-by-line and LUT calculations of absorption coe7cients while keeping all other aspects of the forward calculation the same. For MIPAS, the requirement is that these forward model calculations should agree to better than 10% of the noise equivalent spectral radiance (NESR), i.e., that any errors due to the absorption coe7cient representation are at least a factor 10 smaller than the expected random noise on each measurement. When viewing a 230 K black body, MIPAS signal=noise values vary from 150 to 50 in the spectral ranges used for retrievals, so the NESR=10 criterion corresponds to around 0.1% accuracy for an opaque path. While the NESR=10 value may appear conservative, it should be noted that the noise contribution to retrieval error is random and therefore reduced by adding more measurements or averaging pro:les, but the forward model errors may be systematic and their contribution remains constant. The forward model error limit is not necessarily the best criterion that could be used. A relatively straightforward modi:cation would be to use a single-layer Jacobian analysis [15] to convert this into retrieval error, which would then weight the error of each measurement according to its impact on the retrieval. However, the simpler criterion was chosen mainly on the grounds that it does not distort the residual spectra from the retrievals. 4. Procedure It is assumed that k can be tabulated as a function of wavenumber , pressure p and temperature T only, ignoring any dependence on absorber density. Self-broadening of water vapor lines can be signi:cant in the lower troposphere, but clouds usually prevent infrared limb sounders from viewing these altitudes. The :rst step is to use a line-by-line model to create (large) look-up tables (LUTs) each representing k(; p; T ) for one absorbing species in one microwindow. The wavenumber axis has a spacing 0:0005 cm−1 and extends ±0:175 cm−1 beyond the nominal microwindow boundaries to allow for convolution with the apodized instrument line shape (AILS). The temperature axis spans 180 –310 K, the extreme range of atmospheric temperatures expected, while the pressure axis depends on the range of tangent altitudes for which the microwindow is de:ned (Appendix B). For the
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
749
Table 2 Summary of p, T -axis optimization∗
Retrieval
NLUT
NJ
NJ p
NJ T
Npts
pT CH4 H2 O HNO3 N2 O NO2 O3
60 55 46 48 73 28 36
3330. 4813. 4084. 6701. 4769. 5758. 6701.
26.5 37.1 32.0 29.9 33.8 31.2 28.3
10.3 12.6 11.9 13.8 12.6 8.1 14.0
86 × 106 188 × 106 101 × 106 237 × 106 209 × 106 72 × 106 140 × 106
∗
Columns show total number of tables, average number of pressure, temperature and wavenumber points per table, and total number of data points required for each retrieval.
major absorbers within each microwindow, a total of 2500 (p; T ) points are calculated, with typically 12 points for every factor 10 along the pressure axis and 2 K increments along the temperature axis. The optimization consists of three stages: 1. Subsample the p; T axes. 2. Compress the table using SVD. 3. Subsample the wavenumber axis. These are described in the following subsections. Each stage has associated optimization errors, so an increase of NESR=30 is allocated to each stage to maintain the overall error within NESR=10. 4.1. Pressure and temperature axes Having generated LUTs with high resolution in all three domains (; p; T ), the :rst stage is determine the coarsest p; T grids that can be used, so reducing the overall size of the tables. The original p- and T -axes for each absorber LUT are iteratively subsampled with increasing integer steps. After each iteration, the forward calculation is rerun for a variety of tangent heights and atmospheres, and the process stops when the maximum discrepancy compared to the line-by-line calculation reaches NESR=30 for any spectrum. The results are summarized in Table 2. For example, a total of 55 LUTs are required for the CH4 retrieval. The p-; T -axes in each LUT are reduced to average dimensions 37.1 and 12.6, respectively (∼500(p; T ) points compared to the original 2500), but with an average of 4813 points along the spectral axis a total of 188 million data points are still required. This approach maintains uniform axis increments in ln p and T , which has the advantage of giving fast access to any particular (p; T ) location required for interpolation. Two alternative approaches are noted, potentially yielding further size reduction for the same accuracy, at the expense of more complicated access. One would be to use a temperature axis representing the di4erence relative to some standard pro:le T (p) − T0 (p), thus avoiding redundant ‘dead zones’ within the table such as high absolute temperatures at tropopause pressures. Secondly, the axes could be arbitrarily subsampled, giving irregular axis increments in ln p and T and leaving maximum density where k(ln p; T )
750
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
Fig. 1. Tabulation of the CO2 absorption coe7cient in the ‘PT 0001’ microwindow (686.4 –689:4 cm−1 ). The index of the pT axis shown represents 10 pressures (30 –0:001 mb), repeated for 5 temperatures (180 –312 K), while the wavenumber axis represents 0:05 cm−1 increments. (The original tabulation axes have been undersampled for plotting purposes.)
varies most non-linearly. However, due to the manipulation of the large LUTs, this stage is already the most time-consuming part so procedures requiring yet more iterations are undesirable. 4.2. Compression using SVD After optimizing the p-; T - axes, Table 2 indicates that each retrieval still requires approximately 108 tabulated data points, an inconveniently large number. Strow et al. [10] considered this problem for a nadir-viewing instrument (AIRS). Their solution was to compress the tables using SVD (Appendix A). The idea is that since an absorption coe7cient tabulation contains a large amount of regular structure (Fig. 1) SVD allows the principal patterns to be identi:ed and retained, while the ‘minor’ detail is eliminated. For this procedure, the 3-axis LUT is :rst converted into two-dimensional matrix by multiplexing the p and T -axes into a single dimension, the wavenumber axis forming the other dimension (SVD is a more e7cient compression tool for matrices which are nearly square). The matrix is then decomposed using the leading 30 singular values, a number found to give negligible reconstruction errors. The limited decomposition of matrices containing ∼108 points can be achieved e7ciently with standard software packages [16]. The number of singular values is then reduced iteratively for each absorber in turn, performing forward calculations for each iteration (as in the p-; T - axis optimization) and the process stops when the maximum discrepancy with line-by-line calculations reaches NESR=15 (i.e., NESR=30 from the SVD decomposition=reconstruction, plus NESR=30 from the (p; T ) interpolation). The results are summarized in Table 3. Comparing the values of Npts in Tables 2 and 3, it can be seen that the data volume has been reduced by 2 orders of magnitude. The reconstruction=interpolation of k(; p; T ) from the com-
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
751
Table 3 Summary of SVD compression∗
Retrieval
NLUT
NJ l
NJ NJ l
NJ l NJ z
Npts
pT CH4 H2 O HNO3 N2 O NO2 O3
60 55 46 48 73 28 36
8.0 9.0 7.8 8.9 9.3 8.9 8.8
28 × 103 45 × 103 35 × 103 60 × 103 47 × 103 54 × 103 59 × 103
7 × 103 11 × 103 8 × 103 12 × 103 9 × 103 8 × 103 11 × 103
2:1 × 106 3:1 × 106 2:0 × 106 3:5 × 106 4:1 × 106 1:7 × 106 2:5 × 106
∗
Columns show total number of tables, average number of singular values per table, average number of points in the U and W matrices (Appendix A), and total number of data points required for each retrieval.
pressed representation requires around 10 operations per singular value (Appendix A). For a typical tabulation requiring 10 singular values this means that the break-even point compared to a line-by-line calculation (assuming 10 –100 operations per line) occurs for somewhere between 1 and 10 lines. Since ∼100 lines have to be considered in most regions of the mid-infrared spectrum, using SVD-compressed LUTs would still be expected to produce a signi:cant acceleration compared to line-by-line calculations. Strow et al. note that high and low values of k contribute least to the radiance calculation (opaque √ and transparent limits), so suggest a tabulation of 4 k rather than k in order that the SVD is biased towards :tting the mid-range structure. In our case we have used a ln k tabulation which de-emphasizes the high k values but not the low k, and we found little diUerence between this and a √ 4 k tabulation in terms of the number of singular values required for a given forward model accuracy. This may be because weak absorption features contribute more signi:cantly in limb-viewing mode (against a space background) than in nadir mode where there is usually a background of stronger absorption from the same line lower in the atmosphere. A further advantage of the ln k tabulation is that, having reconstructed the value from the SVD, it is the natural function for interpolation against ln p. This is due to the pressure-dependence of the Lorentz line-shape: k˙1=p at line centre, k˙p in line√wings, so that ln k˙ln p in both regions. Performing the SVD on any other function of k (e.g., 4 k) would :rst require conversion to ln k before performing the interpolation in the (ln p; T ) domain. (Strow et al. avoid this problem since their nadir forward model can use the same :xed pressure values as the tabulation, hence requires no interpolation in the p-domain). 4.3. Wavenumber axis Monochromatic limb radiance spectra (Fig. 2) exhibit structure on a range of scales varying from the narrow, Doppler-broadened line centers from high altitudes, to wide, Lorentz-broadened line wings from low altitudes. In the smoother regions, it is reasonable to expect that the spectrum could be obtained by simple interpolation from nearby points, with computationally expensive full radiative transfer calculations performed only for the remaining quadrature points. The aim of this stage is to :nd such a set of points (‘irregular grid’) for each microwindow.
752
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
Fig. 2. Reconstruction of radiances from linear interpolation. The solid line in the upper plot shows the ‘true’ radiance at 30 km tangent altitude calculated on the 0:0005 cm−1 grid, the dashed line the radiance convolved with the apodized instrument line shape (0:025 cm−1 grid) and the vertical bars the positions of the spectral grid points used for interpolation. The diUerences between the true and interpolated spectra are shown in the lower plot: the solid line represents the true–interpolated diUerences on the :ne grid (left-hand scale) while the dashed line shows the diUerence between the convolved spectra (right hand scale). Radiance units are nW=(sr cm2 cm−1 ) and NESR=10 for this spectral region is 7:0 nW=(sr cm2 cm−1 ).
Analogously with the optimization of the p-, T - axes, one method might be to remove spectral points iteratively from the original k(; p; T ) tabulation and repeat the full forward calculation to check the error. However, as spectral points are removed, the monochromatic radiances calculated for the remaining points are unchanged and diUerences in the forward calculation arise only from the application of the AILS convolution. The actual method used was to start with a set of ‘true’ monochromatic (0:0005 cm−1 resolution) spectra calculated with the line-by-line model for the variety of tangent altitudes and atmospheres. An equivalent set of ‘test’ spectra are calculated using the SVD-compressed LUTs. Then, for each point in the monochromatic spectra an interpolation ‘cost’ is determined: this is the maximum absolute forward model error that would result if the point were replaced by an interpolated value from the surrounding points in the test spectra (this involves convolving the
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
753
Table 4 Summary of -axis optimization∗
Retrieval
NLUT
NJ
NJ NJ l
NJ l NJ z
Npts
pT CH4 H2 O HNO3 N2 O NO2 O3
60 55 46 48 73 28 36
645. 873. 416. 1642. 814. 1306. 2294.
5 × 103 8 × 103 4 × 103 15 × 103 8 × 103 13 × 103 20 × 103
7 × 103 11 × 103 8 × 103 12 × 103 9 × 103 8 × 103 11 × 103
0:8 × 106 1:1 × 106 0:5 × 106 1:3 × 106 1:3 × 106 0:6 × 106 1:1 × 106
∗ Columns show total number of tables, average number of reduced spectral grid points per table, average number of points in the U and W matrices (Appendix A), and total number of data points required for each retrieval.
diUerence between the true and test spectra with the AILS function). The point with the lowest cost is then removed, and the test spectra adjusted to use the interpolated value. The interpolation costs of any aUected points are then recalculated and the process repeated until no point can be removed without generating an interpolation cost which exceeds the NESR=10 error criterion. Fig. 2 illustrates this. Note the high density of grid points around the ‘sharper’ spectral features and that, while the interpolated monochromatic spectra often show large deviation ‘spikes’, these are reduced within the NESR=10 criterion by the AILS convolution. The results are summarized in Table 4. Comparison of the number of spectral points (NJ ) in Tables 4 and 2 indicates that 10 –35% of the original spectral grid points are retained, less compression being obtainable for species with dense line spectra (O3 , NO2 and HNO3 ). The number of radiative transfer calculations, hence CPU time should be reduced in proportion. An additional advantage is that the LUTs only have to be stored for these spectral points so that the U matrix (dimension N × Nl ), which is the larger of the two matrices from the SVD, is reduced by the same factor hence a net reduction in storage requirements (Npts ) from around 3 × 106 to 1 × 106 data points per retrieval. In principle, the spectral grid optimization could be performed before the generation of LUTs. However, an important bene:t of the sequence used is the wavenumber interpolation can be used to oUset some of the errors introduced by the (p; T ) tabulation and SVD compression. This is because the interpolation procedure can produce forward model errors of either sign, and, for the :rst few points, the lowest interpolation cost will be associated any interpolation that reduces the amplitude of the maximum LUT error. Various interpolation functions were investigated (e.g., inverse quadratics to model Lorentzian line wings), producing a factor 2 spread in compression ratios. However, simple linear interpolation was ultimately selected due to its insensitivity to variations in the interpolation points, its low computing cost when reconstructing spectra and its ability to handle negative values. The latter requirement is to handle the Jacobian spectra used within the retrieval, which are derived on the same irregular grid. Single grids have been selected which are applicable for the entire tangent altitude range for each microwindow. Grids could also be selected for each altitude individually, requiring fewer points.
754
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
However, the forward model is signi:cantly more e7cient if monochromatic radiance calculations are duplicated for the full set of tangent altitudes, and it is thought that this outweighs any bene:t from using altitude-dependent grids. 5. Conclusions The MIPAS operational retrieval uses a monochromatic forward model with SVD-compressed LUTs of absorption coe7cient and optimized spectral grids. This allows processing in near real time with relatively modest computational resources and negligible loss of accuracy compared to a line-by-line forward model. The disadvantage is the lack of Sexibility: tables have to be generated in advance for prede:ned microwindows, a procedure requiring one or two days processing per microwindow. Singular value decomposition is used to compress the LUTs by a factor 100, which allows them to be stored in memory during the retrieval. However, the reconstruction with, typically, 10 singular values introduces an overhead approximately equivalent to the cost of a line-by-line calculation including 1–10 lines. The number of monochromatic radiative transfer calculations can be reduced by a factor 3–10 if suitable subsets of the 0:0005 cm−1 grid (‘irregular grids’) are determined for each microwindow. This also yields a further factor 3 reduction in storage of SVD-compressed data. Acknowledgements This work was performed under ESA contract 11886=96=NL=GS. Appendix A. Singular value decomposition Any matrix K (m × n) can be decomposed as the product of three other matrices: K = U$V;
(5)
where U (m × n) and V (n × n) are orthonormal matrices, and $ is a diagonal matrix containing n singular values ([17], for example, contains further details). Assuming that most of the information is contained in the l (n) largest singular values, the decomposition matrices can be truncated in the n dimension to give: KU $ V = U W ;
(6)
where the reduced matrices U (m × l) and W (l × n) are much smaller matrices than the original matrix K, giving a compression factor l=n (assuming mn). In this application, the matrix K represents ln k tabulated for m wavenumber points and n (p; T ) combinations. The compression factor l=n is therefore given by the ratio of the number of singular values to the number of (p; T ) points.
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
755
Table 5 MIPAS microwindows for each species, indicating wavenumber (cm−1 ) and tangent altitude (km) ranges. The microwindows for each retrieval are listed in priority order
Microwindow
Wavenumber
Altitude
Microwindow
Wavenumber
Altitude
PT PT PT PT PT PT PT PT PT PT
0001 0004 0017 0022 0021 0038 0037 0034 0026 0039
686.400 728.300 696.200 1353.325 1932.850 700.475 694.800 1357.200 1349.400 685.700
689.400 729.125 698.375 1354.825 1934.350 701.000 695.100 1358.000 1350.875 685.825
30 15 27 12 12 21 27 9 12 33
68 27 68 47 60 30 36 24 52 47
HNO30001 HNO30006 HNO30012 HNO30003 HNO30021
876.375 885.100 895.675 1324.175 1319.050
879.375 888.100 898.675 1327.175 1322.050
6 6 12 33 12
68 42 68 68 68
CH4 CH4 CH4 CH4 CH4
0001 0012 0013 0005 0022
1350.875 1227.175 1247.775 1256.675 1607.750
1353.875 1230.175 1248.650 1257.650 1610.750
12 6 6 9 15
68 60 30 30 60
N20 N20 N20 N20 N20 N20
0001 0012 0021 0004 0008 0005
1272.050 1233.275 1161.625 1256.675 1265.750 1262.350
1275.050 1236.275 1164.625 1257.975 1266.800 1263.125
12 6 6 9 15 18
60 27 52 30 27 33
H20 H20 H20 H20 H20 H20
0001 0002 0021 0011 0027 0026
1650.025 807.850 1454.525 1574.800 1374.125 1394.475
1653.025 808.450 1457.525 1577.800 1375.075 1395.775
15 9 15 15 12 12
68 18 68 68 24 24
N02 N02 N02 N02 N02
0001 0003 0010 0006 0013
1607.275 1613.725 1619.125 1624.800 1622.550
1610.275 1616.600 1622.125 1627.800 1623.475
15 15 15 47 6
68 68 60 68 30
03 0021 03 0001 03 0012
763.375 1122.800 1073.800
766.375 1125.800 1076.800
6 6 9
68 68 68
The following FORTRAN code illustrates the reconstruction procedure for the spectrum of ln k(; p; T ) (ABSLOG) assuming that U, W contain U , W , and that WII, WIJ, WJI, WJJ contain the pre-computed weights for the four points II, IJ, JI, JJ in the n dimension required for the (p; T ) interpolation DO IV = 1, NV ! loop over wavenumber ABSLOG(IV) = 0.0 ! initialise lnk(nu) = 0 DO IL = 1, NL ! loop over s.v. ABSLOG(IV) = ABSLOG(IV) + U(IV,IL)* & (WII∗ W(IL,II) + WIJ∗ W(IL,IJ)+ & WJI∗ W(IL,JI) + WJJ∗ W(IL,JJ)) END DO END DO Assuming that the wavenumber dimension is large enough for the overhead computation of the weights and interpolation points to be negligible, this corresponds to 10l operations to reconstruct ln k(; p; T ).
756
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
Appendix B. MIPAS microwindows Table 5 lists the microwindows initially selected for the operational processing. Microwindows are constrained to have a maximum width of 3 cm−1 and maximum tangent altitude range 6 –68 km. References [1] Clough SA, Iacono MJ, Moncet J-L. Line-by-line calculations of atmospheric Suxes and cooling rates: application to water vapour. J Geophys Res 1992;97:15761–85. [2] Edwards DP. GENLN2: a general line-by-line atmospheric transmittance and radiance model. Version 3.0— Description and User’s Guide. Report NCAR=TN-367+STR, NCAR, Boulder, CO, USA, 1992. [3] Gordley LL, Marshall BT, Chu DA. LINEPAK: algorithms for modeling spectral transmittance and radiance. JQSRT 1994;52:563–80. [4] Dudhia A. Reference Forward Model (RFM) Software User’s Manual. ESA Document PO-MA-OXF-GS-0003, ESTEC, Noordwijk, 1997 (www.atm.ox.ac.uk/RFM/). [5] Marshall BT, Gordley LL, Chu DA. BANDPAK: algorithms for modeling broadband transmission and radiance. JQSRT 1994;52:581–99. [6] Marks CJ, Rodgers CD. A retrieval method for atmospheric composition from limb emission measurements. J Geophys Res 1993;98:14939–53. [7] McMillin LM, Crone LJ, Kleespies TJ. Atmospheric transmittance of an absorbing gas. 5. Improvements to the OPTRAN approach. Appl Opt 1995;34:8396–9. [8] Saunders R, Matricardi M, Brunel P. An improved fast radiative transfer model for assimilation of satellite radiance observations. QJR Meteorol Soc 1999;125:1407–25. [9] Turner DS. Absorption coe7cient estimation using a two-dimensional interpolation procedure. JQSRT 1995;53:633–77. [10] Strow LL, Motteler HE, Benson RG, Hannon SE, Souza-Machado S de. Fast computation of monochromatic infrared atmospheric transmittances using compressed look-up tables. JQSRT 1998;59:481–93. [11] Ridol: M, Carli B, Carlotti M, Clarmann T von, Dinelli BM, Dudhia A, Flaud J-M, H[opfner M, Morris PE, Raspollini P, Stiller G, Wells RJ. Optimized forward model and retrieval scheme for MIPAS near-real-time data processing. Appl Opt 2000;39:1323–40. [12] Echle G, Clarmann T von, Dudhia A, Flaud J-M, Funke B, Glatthor N, Kerridge BJ, LPopez-Puertas M, MartPQn-Torres FJ, Stiller GP. Optimized spectral microwindows for MIPAS-Envisat data analysis. Appl opt 2000;39:5531–40. [13] Dudhia A, Jay VL, Rodgers CD. Microwindow selection for high spectral resolution sounders. Appl Opt 2001, submitted for publication. [14] HumlPQcR ek J. Optimized computation of the Voigt and complex probability functions. JQSRT 1982;27:437–44. [15] Clarmann T von, Echle G. Selection of optimized microwindows for atmospheric spectroscopy. Appl Opt 1998;37:7661–9. [16] ARPACK software http://www.caam.rice.edu/software/ARPACK/. [17] Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical recipes. Cambridge: Cambridge University Press, 1986.
www.elsevier.com/locate/jqsrt
Fast monochromatic radiative transfer calculations for limb sounding Anu Dudhia ∗ , Paul E. Morris, Robert J. Wells Atmospheric, Oceanic and Planetary Physics, University of Oxford, Parks Rd, Oxford OX1 3PU, UK Received 20 August 2001; accepted 29 November 2001
Abstract Satellite observations of atmospheric infrared spectra can be modeled accurately with line-by-line calculations, but these are too slow to be incorporated into operational retrieval schemes. However, a monochromatic calculation is still feasible if the line-by-line summation is replaced by pre-tabulated absorption coe7cients, requiring a three-way optimization of storage space, accuracy and access time. Such a scheme is used for the operational processing of data from MIPAS, a limb-viewing interferometer. The tabulated data are compressed to a manageable size using singular value decomposition, although the reconstruction adds a small overhead. The number of monochromatic radiative transfer calculations is reduced by determining suitable quadrature points in the spectral domain, which reduces both processing time and data storage requirements. An important aspect of such optimizations is the control of the associated errors in the forward model calculation. The result is an acceleration of monochromatic forward model calculations by one or two orders of magnitude compared to a line-by-line calculation without any signi:cant loss of accuracy, while obtaining data compression factors of 100 or more compared to a direct tabulation of absorption coe7cients. ? 2002 Elsevier Science Ltd. All rights reserved. Keywords: Radiative transfer; Absorption coe7cient look-up tables; Singular value decomposition; Fast forward model; Remote sensing; MIPAS
1. Introduction The inversion of infrared satellite measurements to obtain pro:les of atmospheric temperature and composition generally starts with an assumed pro:le. A ‘forward model’ is then used to calculate the expected measurements from such a pro:le, these are compared with the observed measurements and the pro:le adjusted iteratively until convergence is achieved. The forward model is usually the ∗
Corresponding author. Tel.: +44-1865-272922; fax: +44-1865-272923. E-mail address: [email protected] (A. Dudhia).
0022-4073/02/$ - see front matter ? 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 2 - 4 0 7 3 ( 0 1 ) 0 0 2 8 5 - 0
746
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
most time consuming part of the retrieval process and, while the precision of the retrieved products may be determined by the instrument noise, the accuracy is often limited by the approximations used in the forward model. To model a radiance measurement R requires calculations along the following lines (subscript denoting spectrally varying quantities and superscript i denoting molecule-dependent quantities): R = L d; (1) L = =
B
d ds; ds
i ;
(2) (3)
i
i
i i = exp − k ds ;
(4)
where L is the monochromatic radiance at wavenumber ; is the channel spectral response, B is the Planck function (local thermodynamic equilibrium assumed), the transmittance to the satellite from point s along the ray path, i the transmittance of species i with density i and ki the absorption coe7cient (we will drop the superscript i from now for simplicity). For limb-sounders, there is usually an additional convolution in the elevation-angle domain to represent the :nite :eld-of-view. The most accurate method of evaluating the absorption coe7cient k is a ‘line-by-line’ calculation [1– 4]: a summation of the contributions of all spectroscopic lines in the vicinity, each modeled using the appropriate local path conditions of pressure, temperature, and, occasionally, absorber density. However, such calculations are usually too time consuming to form part of a near realtime retrieval. To model satellite measurements in near real time, the usual approach is to use a band model [5 –8]. For these, the order of the spectral and path integrations in Eqs. (1) and (2) is reversed, so that the monochromatic transmittance can be replaced by a spectrally averaged value , J either computed for the entire channel response or smaller intervals. These band-transmittances are pre-computed for a variety of path conditions, leaving just the path integration in the forward model. While this is extremely fast, band transmittances do not follow Beer’s Law (Eq. (4)) so further assumptions have to be made when obtaining the transmittance increments d =ds J and handling multiple absorbers (J = i Ji ), and this fundamentally limits the accuracy of this approach. Given current computing speeds, a reasonable compromise might be to use the monochromatic radiative transfer (Eqs. (1–3)) but, rather than evaluate k using a line-by-line model, use look-up tables containing values which have been pre-computed for a range of path conditions (e.g. [9]). The main drawback to this direct tabulation approach is that, in order not to introduce signi:cant interpolation errors, the absorption coe7cient tables tend to be rather large and if these cannot be contained within the computer memory signi:cant time can spent accessing the data from disk. Various techniques for compressing the look-up tables have been suggested, such as parameterizing the temperature dependence by a polynomial, or using singular value decomposition (SVD) [10]. These require additional CPU time in order to reconstruct the absorption coe7cients, but this may be outweighed by the bene:ts of reducing the data to a size which can be contained within the computer memory.
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
747
So far, such techniques have only been considered for nadir-viewing instruments. Here we present the approach adopted for a limb-sounding interferometer: the Michelson interferometer for passive atmospheric sounding (MIPAS) due for launch on the Envisat satellite in 2001. The Level 2 operational processing will retrieve pro:les of temperature and composition in near real time [11] with a monochromatic forward model using SVD-compressed look-up tables (LUTs) of absorption coe7cients and further optimized by using quadrature points (‘irregular grids’) for the spectral integration. This paper describes the generation of such tables, the associated optimizations, and discusses various alternatives. 2. CPU considerations MIPAS will acquire spectra covering 685 –2410 cm−1 at 0:025 cm−1 spacing every 4:6 s. A typical limb scan sequence from 68–6 km comprises 17 such spectra and from these data pro:les of temperature and six species are retrieved (although not for all tangent altitudes). Only small subsets (‘microwindows’) of the spectra will be used in the retrievals [12,13], totalling around 10 cm−1 per spectrum per retrieved species (Table 1). With 40 measurements=cm−1 , this means that the near real time processing absorbs ∼500 measurements every second. The spectral integration (Eq. (1)) behaves monochromatically if it is performed at a su7ciently :ne spacing to resolve atmospheric lines. A resolution of 0:0005 cm−1 is usually considered adequate (since Doppler broadening limits lines to a minimum half-width of around 0:001 cm−1 ), giving 50 radiance calculations for each MIPAS measurement on a 0:025 cm−1 grid. By simultaneously modeling multiple tangent paths for the entire limb-scan, a high degree of replication is possible so that each additional spectrum, including FOV convolution, only requires about 10 additional calculations of k for the whole atmospheric path. So to compute transmittances for one gas for one MIPAS measurement requires about 50 × 10 = 500 computations of k . Allowing for 10 absorbing species per microwindow and four retrieval iterations this gives a :gure of 107 k computations every second. To compute the absorption coe7cient with a line-by-line model using some version of the HumlPQcR ek algorithm [14] requires 10 –100 Soating-point operations (FLOPs) per line so, considering 10 –100 local lines, typically 103 FLOPs per absorption coe7cient. Thus, to retrieve MIPAS data in near real time using a line-by-line forward model would require a 10 Giga-FLOP processor. By using look-up tables, where the absorption coe7cients can be obtained in ¡100 operations, Table 1 Summary of the microwindows selected for operational processing of MIPAS data (more details in Appendix B). Columns show the number of microwindows per retrieval, the average number of absorbers per microwindow and the average microwindow width (cm−1 )
Retrieval
NMW
NJ gas
J T
Retrieval
NMW
NJ gas
J T
pT CH4 H2 O HNO3
10 5 6 5
6.0 11.0 7.7 9.6
1.2 2.2 2.0 3.0
N2 O NO2 O3
6 5 3
12.2 5.6 12.0
2.0 2.6 3.0
748
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
this is reduced to ¡1 Giga-FLOPs. In practice, the radiative transfer calculations become the limiting term rather than the k computation, but both can be further reduced by a factor 3–10 if the integration is performed explicitly on a subset of quadrature points from the 0:0005 cm−1 grid. 3. Accuracy criteria The representation of absorption coe7cient by means of look-up tables is a three-way compromise between storage space, CPU time and accuracy, so some accuracy requirement also has to be speci:ed. It is simple to express the accuracy in terms of the diUerence between the interpolated absorption coe7cient and a ‘true’ value calculated with a line-by-line method. However, it is the accuracy of the forward model calculation that is of interest, and this can only be determined by comparing radiative transfer calculations using both a line-by-line model and look-up tables. For this work, the RFM [4] was used since this could use both line-by-line and LUT calculations of absorption coe7cients while keeping all other aspects of the forward calculation the same. For MIPAS, the requirement is that these forward model calculations should agree to better than 10% of the noise equivalent spectral radiance (NESR), i.e., that any errors due to the absorption coe7cient representation are at least a factor 10 smaller than the expected random noise on each measurement. When viewing a 230 K black body, MIPAS signal=noise values vary from 150 to 50 in the spectral ranges used for retrievals, so the NESR=10 criterion corresponds to around 0.1% accuracy for an opaque path. While the NESR=10 value may appear conservative, it should be noted that the noise contribution to retrieval error is random and therefore reduced by adding more measurements or averaging pro:les, but the forward model errors may be systematic and their contribution remains constant. The forward model error limit is not necessarily the best criterion that could be used. A relatively straightforward modi:cation would be to use a single-layer Jacobian analysis [15] to convert this into retrieval error, which would then weight the error of each measurement according to its impact on the retrieval. However, the simpler criterion was chosen mainly on the grounds that it does not distort the residual spectra from the retrievals. 4. Procedure It is assumed that k can be tabulated as a function of wavenumber , pressure p and temperature T only, ignoring any dependence on absorber density. Self-broadening of water vapor lines can be signi:cant in the lower troposphere, but clouds usually prevent infrared limb sounders from viewing these altitudes. The :rst step is to use a line-by-line model to create (large) look-up tables (LUTs) each representing k(; p; T ) for one absorbing species in one microwindow. The wavenumber axis has a spacing 0:0005 cm−1 and extends ±0:175 cm−1 beyond the nominal microwindow boundaries to allow for convolution with the apodized instrument line shape (AILS). The temperature axis spans 180 –310 K, the extreme range of atmospheric temperatures expected, while the pressure axis depends on the range of tangent altitudes for which the microwindow is de:ned (Appendix B). For the
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
749
Table 2 Summary of p, T -axis optimization∗
Retrieval
NLUT
NJ
NJ p
NJ T
Npts
pT CH4 H2 O HNO3 N2 O NO2 O3
60 55 46 48 73 28 36
3330. 4813. 4084. 6701. 4769. 5758. 6701.
26.5 37.1 32.0 29.9 33.8 31.2 28.3
10.3 12.6 11.9 13.8 12.6 8.1 14.0
86 × 106 188 × 106 101 × 106 237 × 106 209 × 106 72 × 106 140 × 106
∗
Columns show total number of tables, average number of pressure, temperature and wavenumber points per table, and total number of data points required for each retrieval.
major absorbers within each microwindow, a total of 2500 (p; T ) points are calculated, with typically 12 points for every factor 10 along the pressure axis and 2 K increments along the temperature axis. The optimization consists of three stages: 1. Subsample the p; T axes. 2. Compress the table using SVD. 3. Subsample the wavenumber axis. These are described in the following subsections. Each stage has associated optimization errors, so an increase of NESR=30 is allocated to each stage to maintain the overall error within NESR=10. 4.1. Pressure and temperature axes Having generated LUTs with high resolution in all three domains (; p; T ), the :rst stage is determine the coarsest p; T grids that can be used, so reducing the overall size of the tables. The original p- and T -axes for each absorber LUT are iteratively subsampled with increasing integer steps. After each iteration, the forward calculation is rerun for a variety of tangent heights and atmospheres, and the process stops when the maximum discrepancy compared to the line-by-line calculation reaches NESR=30 for any spectrum. The results are summarized in Table 2. For example, a total of 55 LUTs are required for the CH4 retrieval. The p-; T -axes in each LUT are reduced to average dimensions 37.1 and 12.6, respectively (∼500(p; T ) points compared to the original 2500), but with an average of 4813 points along the spectral axis a total of 188 million data points are still required. This approach maintains uniform axis increments in ln p and T , which has the advantage of giving fast access to any particular (p; T ) location required for interpolation. Two alternative approaches are noted, potentially yielding further size reduction for the same accuracy, at the expense of more complicated access. One would be to use a temperature axis representing the di4erence relative to some standard pro:le T (p) − T0 (p), thus avoiding redundant ‘dead zones’ within the table such as high absolute temperatures at tropopause pressures. Secondly, the axes could be arbitrarily subsampled, giving irregular axis increments in ln p and T and leaving maximum density where k(ln p; T )
750
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
Fig. 1. Tabulation of the CO2 absorption coe7cient in the ‘PT 0001’ microwindow (686.4 –689:4 cm−1 ). The index of the pT axis shown represents 10 pressures (30 –0:001 mb), repeated for 5 temperatures (180 –312 K), while the wavenumber axis represents 0:05 cm−1 increments. (The original tabulation axes have been undersampled for plotting purposes.)
varies most non-linearly. However, due to the manipulation of the large LUTs, this stage is already the most time-consuming part so procedures requiring yet more iterations are undesirable. 4.2. Compression using SVD After optimizing the p-; T - axes, Table 2 indicates that each retrieval still requires approximately 108 tabulated data points, an inconveniently large number. Strow et al. [10] considered this problem for a nadir-viewing instrument (AIRS). Their solution was to compress the tables using SVD (Appendix A). The idea is that since an absorption coe7cient tabulation contains a large amount of regular structure (Fig. 1) SVD allows the principal patterns to be identi:ed and retained, while the ‘minor’ detail is eliminated. For this procedure, the 3-axis LUT is :rst converted into two-dimensional matrix by multiplexing the p and T -axes into a single dimension, the wavenumber axis forming the other dimension (SVD is a more e7cient compression tool for matrices which are nearly square). The matrix is then decomposed using the leading 30 singular values, a number found to give negligible reconstruction errors. The limited decomposition of matrices containing ∼108 points can be achieved e7ciently with standard software packages [16]. The number of singular values is then reduced iteratively for each absorber in turn, performing forward calculations for each iteration (as in the p-; T - axis optimization) and the process stops when the maximum discrepancy with line-by-line calculations reaches NESR=15 (i.e., NESR=30 from the SVD decomposition=reconstruction, plus NESR=30 from the (p; T ) interpolation). The results are summarized in Table 3. Comparing the values of Npts in Tables 2 and 3, it can be seen that the data volume has been reduced by 2 orders of magnitude. The reconstruction=interpolation of k(; p; T ) from the com-
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
751
Table 3 Summary of SVD compression∗
Retrieval
NLUT
NJ l
NJ NJ l
NJ l NJ z
Npts
pT CH4 H2 O HNO3 N2 O NO2 O3
60 55 46 48 73 28 36
8.0 9.0 7.8 8.9 9.3 8.9 8.8
28 × 103 45 × 103 35 × 103 60 × 103 47 × 103 54 × 103 59 × 103
7 × 103 11 × 103 8 × 103 12 × 103 9 × 103 8 × 103 11 × 103
2:1 × 106 3:1 × 106 2:0 × 106 3:5 × 106 4:1 × 106 1:7 × 106 2:5 × 106
∗
Columns show total number of tables, average number of singular values per table, average number of points in the U and W matrices (Appendix A), and total number of data points required for each retrieval.
pressed representation requires around 10 operations per singular value (Appendix A). For a typical tabulation requiring 10 singular values this means that the break-even point compared to a line-by-line calculation (assuming 10 –100 operations per line) occurs for somewhere between 1 and 10 lines. Since ∼100 lines have to be considered in most regions of the mid-infrared spectrum, using SVD-compressed LUTs would still be expected to produce a signi:cant acceleration compared to line-by-line calculations. Strow et al. note that high and low values of k contribute least to the radiance calculation (opaque √ and transparent limits), so suggest a tabulation of 4 k rather than k in order that the SVD is biased towards :tting the mid-range structure. In our case we have used a ln k tabulation which de-emphasizes the high k values but not the low k, and we found little diUerence between this and a √ 4 k tabulation in terms of the number of singular values required for a given forward model accuracy. This may be because weak absorption features contribute more signi:cantly in limb-viewing mode (against a space background) than in nadir mode where there is usually a background of stronger absorption from the same line lower in the atmosphere. A further advantage of the ln k tabulation is that, having reconstructed the value from the SVD, it is the natural function for interpolation against ln p. This is due to the pressure-dependence of the Lorentz line-shape: k˙1=p at line centre, k˙p in line√wings, so that ln k˙ln p in both regions. Performing the SVD on any other function of k (e.g., 4 k) would :rst require conversion to ln k before performing the interpolation in the (ln p; T ) domain. (Strow et al. avoid this problem since their nadir forward model can use the same :xed pressure values as the tabulation, hence requires no interpolation in the p-domain). 4.3. Wavenumber axis Monochromatic limb radiance spectra (Fig. 2) exhibit structure on a range of scales varying from the narrow, Doppler-broadened line centers from high altitudes, to wide, Lorentz-broadened line wings from low altitudes. In the smoother regions, it is reasonable to expect that the spectrum could be obtained by simple interpolation from nearby points, with computationally expensive full radiative transfer calculations performed only for the remaining quadrature points. The aim of this stage is to :nd such a set of points (‘irregular grid’) for each microwindow.
752
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
Fig. 2. Reconstruction of radiances from linear interpolation. The solid line in the upper plot shows the ‘true’ radiance at 30 km tangent altitude calculated on the 0:0005 cm−1 grid, the dashed line the radiance convolved with the apodized instrument line shape (0:025 cm−1 grid) and the vertical bars the positions of the spectral grid points used for interpolation. The diUerences between the true and interpolated spectra are shown in the lower plot: the solid line represents the true–interpolated diUerences on the :ne grid (left-hand scale) while the dashed line shows the diUerence between the convolved spectra (right hand scale). Radiance units are nW=(sr cm2 cm−1 ) and NESR=10 for this spectral region is 7:0 nW=(sr cm2 cm−1 ).
Analogously with the optimization of the p-, T - axes, one method might be to remove spectral points iteratively from the original k(; p; T ) tabulation and repeat the full forward calculation to check the error. However, as spectral points are removed, the monochromatic radiances calculated for the remaining points are unchanged and diUerences in the forward calculation arise only from the application of the AILS convolution. The actual method used was to start with a set of ‘true’ monochromatic (0:0005 cm−1 resolution) spectra calculated with the line-by-line model for the variety of tangent altitudes and atmospheres. An equivalent set of ‘test’ spectra are calculated using the SVD-compressed LUTs. Then, for each point in the monochromatic spectra an interpolation ‘cost’ is determined: this is the maximum absolute forward model error that would result if the point were replaced by an interpolated value from the surrounding points in the test spectra (this involves convolving the
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
753
Table 4 Summary of -axis optimization∗
Retrieval
NLUT
NJ
NJ NJ l
NJ l NJ z
Npts
pT CH4 H2 O HNO3 N2 O NO2 O3
60 55 46 48 73 28 36
645. 873. 416. 1642. 814. 1306. 2294.
5 × 103 8 × 103 4 × 103 15 × 103 8 × 103 13 × 103 20 × 103
7 × 103 11 × 103 8 × 103 12 × 103 9 × 103 8 × 103 11 × 103
0:8 × 106 1:1 × 106 0:5 × 106 1:3 × 106 1:3 × 106 0:6 × 106 1:1 × 106
∗ Columns show total number of tables, average number of reduced spectral grid points per table, average number of points in the U and W matrices (Appendix A), and total number of data points required for each retrieval.
diUerence between the true and test spectra with the AILS function). The point with the lowest cost is then removed, and the test spectra adjusted to use the interpolated value. The interpolation costs of any aUected points are then recalculated and the process repeated until no point can be removed without generating an interpolation cost which exceeds the NESR=10 error criterion. Fig. 2 illustrates this. Note the high density of grid points around the ‘sharper’ spectral features and that, while the interpolated monochromatic spectra often show large deviation ‘spikes’, these are reduced within the NESR=10 criterion by the AILS convolution. The results are summarized in Table 4. Comparison of the number of spectral points (NJ ) in Tables 4 and 2 indicates that 10 –35% of the original spectral grid points are retained, less compression being obtainable for species with dense line spectra (O3 , NO2 and HNO3 ). The number of radiative transfer calculations, hence CPU time should be reduced in proportion. An additional advantage is that the LUTs only have to be stored for these spectral points so that the U matrix (dimension N × Nl ), which is the larger of the two matrices from the SVD, is reduced by the same factor hence a net reduction in storage requirements (Npts ) from around 3 × 106 to 1 × 106 data points per retrieval. In principle, the spectral grid optimization could be performed before the generation of LUTs. However, an important bene:t of the sequence used is the wavenumber interpolation can be used to oUset some of the errors introduced by the (p; T ) tabulation and SVD compression. This is because the interpolation procedure can produce forward model errors of either sign, and, for the :rst few points, the lowest interpolation cost will be associated any interpolation that reduces the amplitude of the maximum LUT error. Various interpolation functions were investigated (e.g., inverse quadratics to model Lorentzian line wings), producing a factor 2 spread in compression ratios. However, simple linear interpolation was ultimately selected due to its insensitivity to variations in the interpolation points, its low computing cost when reconstructing spectra and its ability to handle negative values. The latter requirement is to handle the Jacobian spectra used within the retrieval, which are derived on the same irregular grid. Single grids have been selected which are applicable for the entire tangent altitude range for each microwindow. Grids could also be selected for each altitude individually, requiring fewer points.
754
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
However, the forward model is signi:cantly more e7cient if monochromatic radiance calculations are duplicated for the full set of tangent altitudes, and it is thought that this outweighs any bene:t from using altitude-dependent grids. 5. Conclusions The MIPAS operational retrieval uses a monochromatic forward model with SVD-compressed LUTs of absorption coe7cient and optimized spectral grids. This allows processing in near real time with relatively modest computational resources and negligible loss of accuracy compared to a line-by-line forward model. The disadvantage is the lack of Sexibility: tables have to be generated in advance for prede:ned microwindows, a procedure requiring one or two days processing per microwindow. Singular value decomposition is used to compress the LUTs by a factor 100, which allows them to be stored in memory during the retrieval. However, the reconstruction with, typically, 10 singular values introduces an overhead approximately equivalent to the cost of a line-by-line calculation including 1–10 lines. The number of monochromatic radiative transfer calculations can be reduced by a factor 3–10 if suitable subsets of the 0:0005 cm−1 grid (‘irregular grids’) are determined for each microwindow. This also yields a further factor 3 reduction in storage of SVD-compressed data. Acknowledgements This work was performed under ESA contract 11886=96=NL=GS. Appendix A. Singular value decomposition Any matrix K (m × n) can be decomposed as the product of three other matrices: K = U$V;
(5)
where U (m × n) and V (n × n) are orthonormal matrices, and $ is a diagonal matrix containing n singular values ([17], for example, contains further details). Assuming that most of the information is contained in the l (n) largest singular values, the decomposition matrices can be truncated in the n dimension to give: KU $ V = U W ;
(6)
where the reduced matrices U (m × l) and W (l × n) are much smaller matrices than the original matrix K, giving a compression factor l=n (assuming mn). In this application, the matrix K represents ln k tabulated for m wavenumber points and n (p; T ) combinations. The compression factor l=n is therefore given by the ratio of the number of singular values to the number of (p; T ) points.
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
755
Table 5 MIPAS microwindows for each species, indicating wavenumber (cm−1 ) and tangent altitude (km) ranges. The microwindows for each retrieval are listed in priority order
Microwindow
Wavenumber
Altitude
Microwindow
Wavenumber
Altitude
PT PT PT PT PT PT PT PT PT PT
0001 0004 0017 0022 0021 0038 0037 0034 0026 0039
686.400 728.300 696.200 1353.325 1932.850 700.475 694.800 1357.200 1349.400 685.700
689.400 729.125 698.375 1354.825 1934.350 701.000 695.100 1358.000 1350.875 685.825
30 15 27 12 12 21 27 9 12 33
68 27 68 47 60 30 36 24 52 47
HNO30001 HNO30006 HNO30012 HNO30003 HNO30021
876.375 885.100 895.675 1324.175 1319.050
879.375 888.100 898.675 1327.175 1322.050
6 6 12 33 12
68 42 68 68 68
CH4 CH4 CH4 CH4 CH4
0001 0012 0013 0005 0022
1350.875 1227.175 1247.775 1256.675 1607.750
1353.875 1230.175 1248.650 1257.650 1610.750
12 6 6 9 15
68 60 30 30 60
N20 N20 N20 N20 N20 N20
0001 0012 0021 0004 0008 0005
1272.050 1233.275 1161.625 1256.675 1265.750 1262.350
1275.050 1236.275 1164.625 1257.975 1266.800 1263.125
12 6 6 9 15 18
60 27 52 30 27 33
H20 H20 H20 H20 H20 H20
0001 0002 0021 0011 0027 0026
1650.025 807.850 1454.525 1574.800 1374.125 1394.475
1653.025 808.450 1457.525 1577.800 1375.075 1395.775
15 9 15 15 12 12
68 18 68 68 24 24
N02 N02 N02 N02 N02
0001 0003 0010 0006 0013
1607.275 1613.725 1619.125 1624.800 1622.550
1610.275 1616.600 1622.125 1627.800 1623.475
15 15 15 47 6
68 68 60 68 30
03 0021 03 0001 03 0012
763.375 1122.800 1073.800
766.375 1125.800 1076.800
6 6 9
68 68 68
The following FORTRAN code illustrates the reconstruction procedure for the spectrum of ln k(; p; T ) (ABSLOG) assuming that U, W contain U , W , and that WII, WIJ, WJI, WJJ contain the pre-computed weights for the four points II, IJ, JI, JJ in the n dimension required for the (p; T ) interpolation DO IV = 1, NV ! loop over wavenumber ABSLOG(IV) = 0.0 ! initialise lnk(nu) = 0 DO IL = 1, NL ! loop over s.v. ABSLOG(IV) = ABSLOG(IV) + U(IV,IL)* & (WII∗ W(IL,II) + WIJ∗ W(IL,IJ)+ & WJI∗ W(IL,JI) + WJJ∗ W(IL,JJ)) END DO END DO Assuming that the wavenumber dimension is large enough for the overhead computation of the weights and interpolation points to be negligible, this corresponds to 10l operations to reconstruct ln k(; p; T ).
756
A. Dudhia et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 745 – 756
Appendix B. MIPAS microwindows Table 5 lists the microwindows initially selected for the operational processing. Microwindows are constrained to have a maximum width of 3 cm−1 and maximum tangent altitude range 6 –68 km. References [1] Clough SA, Iacono MJ, Moncet J-L. Line-by-line calculations of atmospheric Suxes and cooling rates: application to water vapour. J Geophys Res 1992;97:15761–85. [2] Edwards DP. GENLN2: a general line-by-line atmospheric transmittance and radiance model. Version 3.0— Description and User’s Guide. Report NCAR=TN-367+STR, NCAR, Boulder, CO, USA, 1992. [3] Gordley LL, Marshall BT, Chu DA. LINEPAK: algorithms for modeling spectral transmittance and radiance. JQSRT 1994;52:563–80. [4] Dudhia A. Reference Forward Model (RFM) Software User’s Manual. ESA Document PO-MA-OXF-GS-0003, ESTEC, Noordwijk, 1997 (www.atm.ox.ac.uk/RFM/). [5] Marshall BT, Gordley LL, Chu DA. BANDPAK: algorithms for modeling broadband transmission and radiance. JQSRT 1994;52:581–99. [6] Marks CJ, Rodgers CD. A retrieval method for atmospheric composition from limb emission measurements. J Geophys Res 1993;98:14939–53. [7] McMillin LM, Crone LJ, Kleespies TJ. Atmospheric transmittance of an absorbing gas. 5. Improvements to the OPTRAN approach. Appl Opt 1995;34:8396–9. [8] Saunders R, Matricardi M, Brunel P. An improved fast radiative transfer model for assimilation of satellite radiance observations. QJR Meteorol Soc 1999;125:1407–25. [9] Turner DS. Absorption coe7cient estimation using a two-dimensional interpolation procedure. JQSRT 1995;53:633–77. [10] Strow LL, Motteler HE, Benson RG, Hannon SE, Souza-Machado S de. Fast computation of monochromatic infrared atmospheric transmittances using compressed look-up tables. JQSRT 1998;59:481–93. [11] Ridol: M, Carli B, Carlotti M, Clarmann T von, Dinelli BM, Dudhia A, Flaud J-M, H[opfner M, Morris PE, Raspollini P, Stiller G, Wells RJ. Optimized forward model and retrieval scheme for MIPAS near-real-time data processing. Appl Opt 2000;39:1323–40. [12] Echle G, Clarmann T von, Dudhia A, Flaud J-M, Funke B, Glatthor N, Kerridge BJ, LPopez-Puertas M, MartPQn-Torres FJ, Stiller GP. Optimized spectral microwindows for MIPAS-Envisat data analysis. Appl opt 2000;39:5531–40. [13] Dudhia A, Jay VL, Rodgers CD. Microwindow selection for high spectral resolution sounders. Appl Opt 2001, submitted for publication. [14] HumlPQcR ek J. Optimized computation of the Voigt and complex probability functions. JQSRT 1982;27:437–44. [15] Clarmann T von, Echle G. Selection of optimized microwindows for atmospheric spectroscopy. Appl Opt 1998;37:7661–9. [16] ARPACK software http://www.caam.rice.edu/software/ARPACK/. [17] Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical recipes. Cambridge: Cambridge University Press, 1986.