The physical retrieval methodology for IASI: the δ-IASI code

Environmental Modelling & Software 20 (2005) 1111–1126 www.elsevier.com/locate/envsoft

The physical retrieval methodology for IASI: the d-IASI code A. Carissimoa, I. De Feisb, C. Serioa,c,* a

Dipartimento di Ingegneria e Fisica dell’Ambiente, Universita´ della Basilicata, Via dell’Ateneo Lucano, 10, 85100 Potenza, Italy b Istituto per le Applicazioni del Calcolo del CNR, Via Pietro Castellino 111, 80131 Napoli, Italy c Istituto Nazionale per la Fisica della Materia, Unita` di Napoli, Gruppo Coordinato di Potenza, C.da Macchia Romana, Via dell’ Ateneo Lucano 10, 85100 Potenza, Italy Received 17 December 2003; received in revised form 2 July 2004; accepted 9 July 2004

Abstract This paper describes a physical-based methodology for the retrieval of geophysical parameters (temperature, water vapor and ozone) from highly resolved infrared radiance, and presents the algorithm which implements the procedure. The algorithm we have implemented is mostly intended for the Infrared Atmospheric Sounding Interferometer which is planned to be flown on the first European Meteorological Operational Satellite (Metop/1) in 2006. Nevertheless, with minor modifications, the code is well suited for any nadir viewing satellite and airborne infrared sensor with a sampling rate in the range of 0.1–2 cm
1. Basically, the implementation of the inverse scheme follows Rodgers’ Statistical Regularization method. However, an additional regularization parameter is introduced in the inverse scheme which gains to the algorithm the capability of improving the retrieval accuracy and to constraint the step size of Newton updates in such a way to lead iterates toward the feasible region of the inverse solution. Although, the paper mostly focuses on documenting and discussing the mathematical details of the inverse method, retrieval exercises have been provided, which exemplify the use and potential performance of the method. These retrieval exercises have been performed for the Infrared Atmospheric Sounding Interferometer. In addition, examples of application to real observations have been discussed based on the Interferometric Monitoring Greenhouse (IMG) gases Fourier Transform Spectrometer which has flown on the Japanese Advanced Earth Observation Satellite.
2004 Elsevier Ltd. All rights reserved. Keywords: Computer modeling and simulation; Mathematical Inversion; Remote observing techniques; Infrared

Software availability d-IASI is a FORTRAN based inverse radiative transfer code designed to match the spectral range (645– 2700 cm
1) of the Infrared Atmospheric Sounding Interferometer (IASI). With minor modifications, it could also be used to retrieve geophysical parameters from other sounding instruments. The functionality of the dIASI software package has been set up to be as indepen* Corresponding author. Istituto Nazionale per la Fisica della Materia, Unita` di Napoli, Gruppo Coordinato di Potenza, C. da Macchia Romana, Via dell’ Ateneo Lucano 10, 85100 Potenza, Italy. Tel.: C39 0971 427 261; fax: C39 0971 427 271. E-mail address: [email protected] (C. Serio). 1364-8152/$ - see front matter
2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2004.07.003

dent as possible of the details of the computer environment in which it is run. The package has been developed for the ALPHA UNIX and HP UNIX workstations. A release for MS-Windows PC platforms is available. Minimal requirements: 1 free GB hard disk and 128 MB RAM. Upon EUMETSAT approval, d-IASI UNIX implementation is available from Prof. Carmine Serio (e-mail: [email protected]).

1. Introduction The Infrared Atmospheric Sounding Interferometer, IASI (Cayla, 1995), is a spectrometer interferometer

1112

A. Carissimo et al. / Environmental Modelling & Software 20 (2005) 1111–1126

(Michelson design) which is presently under development at the laboratories of the French National Space Agency (CNES). The instrument is a part of the operational meteorological payload managed by the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT). IASI will be flown in 2006 on the operational meteorological satellite Metop in the frame of the European Polar System (EPS). The spectrometer will be launched to support operational meteorology and climate monitoring. The package is being specified to possess the high spectral resolution required to achieve improvement in the sounding quality, which the community of meteorologists considers mandatory in order to improve Numerical Weather Prediction. In addition, the instrument is expected to give information on greenhouse trace gases and on surface and cloud properties: sea surface temperature, land surface temperature and emissivity, cloud transmittance and reflectance. IASI belongs to the class of new generation space borne sensors which will be put in orbit within this decade. They include AIRS (Atmospheric InfraRed Sounder) on Aqua (launched in April 2002), GIFTS (Geosynchronous Imaging Fourier Spectrometer) on EOS-3 (2005), CrIS (Cross Track Infrared Sounder) on the NPP (2006) and NPOESS (2008) satellites. All of these spectrometers are characterized by a broad band spectral coverage (3.7 to 15.5 mm) and a spectral sampling rate in the range 0.25–2 cm
1. This paper describes and fully documents the physical inversion of infrared spectral radiance for geophysical parameters, such as temperature, water vapor and ozone profiles, in a form which is suitable for the next generation infrared sensors aiming at meteorological missions, and presents its software implementation. The radiative transfer model embedded in the inversion scheme is s-IASI (Amato et al., 2002), a line-by-line fast forward model which has the capability of calculating analytical Jacobian matrix for the geophysical parameters. The present d-IASI implementation includes s-IASI in the form of a subroutine, so that the user does not have to provide additional forward calculations. Furthermore, the user could easily replace the radiative transfer module, if needed, to speed up the retrieval process. The radiative transfer module is the one which is the most time consuming because of line-by-line calculations. Although computationally efficient, the on-line radiative transfer calculations slow down d-IASI and make it not particularly suitable for applications which need real time processing. However, d-IASI is mostly oriented for research applications rather than operational end-uses. It incorporates, indeed, a full lineby-line forward module (that is s-IASI, Amato et al., 2002) and may be run in several of many modes, which include simultaneous or sequential retrieval strategies, user selection of the geophysical parameters to be

retrieved and user choice of the spectral ranges or channels to be used in the inversion process. Research uses include: retrieval and error analysis for the assessment of the expected performance of infrared sensors, off-line inversions for geophysical parameters mostly intended for climatology and quality check of non-physical inversion scheme products such as retrieval by linear Least Squares regression and neural network. Freezing some options and modes of d-IASI in order to speed up the line-by-line calculations, the code could also be used to back up operational non-physical inversions, that is it could be used to provide just one iteration for temperature and water vapor using the output of a suitable regression or neural network retrieval scheme. The use of d-IASI has been exemplified by considering retrieval exercises applied to both simulated and real observations. Simulation exercises for IASI have been performed for a wide range of atmospheric states which are representative of the global meteorological variability. Examples of application to real observations have been discussed by considering spectral infrared radiance recorded through the Japanese Fourier transform spectrometer IMG (Interferometric Monitoring of Greenhouse gases). IMG (Kobayashi et al., 1999) is a Fourier transform spectrometer which has successfully flown on board the Japanese polar platform Advanced Earth Observing Satellite (ADEOS) from August 1996 to June 1997. An early version of the scheme described in this paper may be found in Amato and Serio (2000). Applications of the scheme to infrared spectral radiance for the inversion of geophysical parameters have been mostly described and discussed (Amato et al., 1999; Lubrano et al., 2000; Masiello et al., 2002; Lubrano et al., 2002) in the framework of the IMG project. The paper is organized as follows. In Section 2, the radiative transfer embedded in d-IASI is briefly summarized and then the inversion method is discussed. Section 3 is dedicated to the retrieval exercises which exemplify the inversion procedure for clear sky. Application of the inversion approach to IMG is discussed in this section, as well. Conclusions are drawn in Section 4.

2. Methodology 2.1. The forward module The radiative transfer calculations performed within d-IASI are carried out in the forward module subroutine s-IASI (Amato et al., 2002). This is a fast line-by-line model that, in addition to spectral radiance, can compute analytical Jacobian matrixes for the major geophysical parameters (e.g. temperature, water vapor, ozone, carbon dioxide, nitrous oxide, surface emissivity).

A. Carissimo et al. / Environmental Modelling & Software 20 (2005) 1111–1126

The code computes monochromatic radiance from look-up tables of monochromatic layer optical depth generated using the line-by-line model LBLRTM (e.g. Clough et al., 1992; Clough and Iacono, 1995). Following the LBLRTM strategy, for each individual species, the monochromatic layer optical depth is generated and stored in the appropriate look-up table at a sampling rate, Ds, equal to 0.25 times the width of the narrowest gas absorption line, which gives Ds Z 0.000169 cm
1 for the highest-altitude layer and Ds Z 0.01 cm
1 for the lower troposphere. The code s-IASI borrows from LBLRTM the CO2 line mixing scheme. Cross sections of heavy molecules such as chlorofluorocarbons are not parameterized in terms of look-up tables, since absorption can be computed quickly and separately. The code s-IASI accounts for the presence of CFC species, namely CCl3F (CFC-11) and CCl2F2 (CFC-12). Water vapor continuum absorption is computed separately; we use the CKD standard (Clough et al., 1989), version 2.4. All the required spectroscopic parameters have been extracted from the database HITRAN 2000 (Rothman et al., 2003). The atmospheric layering embedded in s-IASI consists of a pressure layer grid of Nl Z 44 grid points (e.g., Amato et al., 2002). The layering of the atmosphere we adopt in s-IASI has been first developed and used by Matricardi and Saunders (1999) who showed that the error introduced by limiting the number of layers to 43 is negligible. This was done by computing the difference between brightness temperature spectra obtained by dividing the atmosphere into 43 and 98 layers, respectively. A similar exercise has been repeated by us and we have found that differences reach at most 0.1 K. If the user needs higher accuracy, a version is also available which considers a grid of 60 layers, although in this way the radiative transfer is slowed down. It should be considered that in the present version, the inverse module takes a few milliseconds to run against minutes taken by the forward module. The code s-IASI allows one to input any user defined surface emissivity. The default in the code is the Masuda model (Masuda et al., 1988) for sea surface emissivity. Finally, s-IASI computes the surface emissivity Jacobian matrix and this quantity may be either a user supplied parameter or a retrieval parameter. Further details about s-IASI may be obtained in Amato et al. (2002).

2.2. The inversion methodology Let Rt Z (R1,.,RN) be the up welling spectral radiance vector (here and afterwards the suffix t denotes transpose), where N indicates the number of measured spectral radiance data points. Let F(v) be the forward model function which relates R to the atmospheric state

1113

parameter vector, vt Z (v1,.,vM), where the M parameters may include temperature profile, gas constituent profiles and surface emissivity function. We expand the direct model equation as a Taylor series around a suitable first guess v0 of the solution: RZFðvÞZR0 CKðv
v0 ÞChigher order termsCn3

ð1Þ

with R0 Z F(v0) and K the N by M Jacobian matrix of the function F(v) computed at v Z v0; n3 is an additive noise term which accounts for observational and forward model errors. These are assumed to be Gaussian and zero mean, that is E(n3) Z 0, with E( ) denoting expectation value. Setting m Z v
v0 and d Z R
R0, Eq. (1) can be written in the more compact form dZKmChigher order termsCn3 :

ð2Þ

Eq. (1) has to be inverted for v or, equivalently m. A suitable procedure to obtain a stable estimate ^v of the true state v is to minimize the distance or norm between v and a suitable background state, vg (Rodgers, 1976; Twomey, 1977; Tarantola, 1987). In addition, we require the observations and fitted model to be close to each other in the Least Squares sense. Mathematically, the inverse problem is stated as 8
t
 < v
vg L v
vg min! subject to
 : ðR
FðvÞÞt S
1 R
FðvÞ %c2a

ð3Þ

where L is a suitable smoothing matrix and S is the covariance matrix of the observations: S Z OR C OF, with OR the radiometric noise matrix and OF the forward model error covariance matrix. Furthermore, c2a is the value of a c2 variable with N degree of freedom at the confidence interval 1
a. In practice N is the number of spectral radiance data points; if, as usual, pffiffiffiffiffiffiffiwe consider a 2s tolerance interval, then c2a ZNC2 2N. Following Rodgers’ method (Rodgers, 1976, 2000) we choose L Z B
1, where B is the a priori covariance matrix of vg. Using a covariance operator has the additional advantage of homogenizing the possible different physical units in v, through appropriate physical-unit normalization. It can be shown (Colton and Kress, 1992) that the constrained optimization problem (3) is equivalent to the quadratic optimization problem t
 
t QðvÞZðR
FðvÞÞ S
1 ðR
FðvÞÞCg v
vg B
1 v
vg ð4Þ where g R 0 and whose minimum may be sought by the Newton method

1114

A. Carissimo et al. / Environmental Modelling & Software 20 (2005) 1111–1126



 ^vnC1 Z^vn
Vv Vv Q
1 Vv Q

ð5Þ

where Vv denotes the gradient operation with respect to the variable v. If as usual, we retain only the first order term in the Taylor expansion for F(v), and use this linear approximation in the quadratic form (4), then the minimization scheme above is more correctly referred to as the Gauss– Newton method. Within this scheme, we use d  KmCn3

ð6Þ

and the quadratic form may be written as   t t QðmÞZðd
KmÞ S
1 ðd
KmÞCg ðmCq0 Þ B
1 ðmCq0 Þ ð7Þ where we have used v
vg Z (v0
vg) C m Z q0 C m. Considering for now a truly linear problem the Gauss– Newton converge in one step and gives the solution  
1  t
1 ^ K S d
gB
1 q0 mZ gB
1 CKt S
1 K


ð8Þ

and if we choose g Z 1, the usual statistically regularized solution is obtained (Rodgers, 1976, 2000). The advantage of having a further regularization parameter, g has been described in Amato et al. (1996). This parameter adds flexibility to the usual statistical regularization scheme since, in practice, the a priori constraint, B could not be the best possible a priori information for the problem at hand. It could be too tight or too loose as it may happen when B is obtained, e.g., by climatology. The additional regularization parameter, g may be tuned to trade off between data and a priori constraint. Some of the results about the use of g are here anticipated.  It will be shown in Section 2.2.2 that the introduction of the parameter g generalizes the Statistical ^ gZ1 to the class of ridge Regularization estimator, m estimators (Kendall and Stuart, 1979). Unlike g Z 1, a suitable choice of g can improve the root means square error (Kendall and Stuart, 1979), that is the ^ and the Euclidean distance between the estimator m true state vector, m.  For a non-linear problem, the solution above needs to be iterated and setting g Z 1 we cannot be sure that the Gauss–Newton algorithm will converge. Then, the introduction of the g-regularization parameter has the effect to dump the Newton updates (see Section 2.2.2) since it can diminish the a priori variances in B. In other words, g dumps the ^ nC1
m ^ n between two consecuwidth of the step m tive iterates and leads its direction toward the steepest descent of the cost function. This has ^ nC1 the action of preventing the current update m to jump beyond the linear domain around the ^ n . This approach is almost equivalent current guess m

(Tarantola, 1987) to the so-called Levenberg–Marquardt method (Marquardt, 1963). Eq. (8) may be put in terms of the difference between the solution vector ^v and vg, a form which will be extensively used in the following: 
^ gB
1 CKt S
1 K
1 Kt S
1 y xZ 0

ð9Þ

^ ^v
vg and y0 Z d C Kq0; q0 Z v0
vg. with xZ The problem of estimating the regularization parameter will be addressed in Section 2.2.3, here we focus on the covariance of the inverse solution. If we assume Gaussian errors and that the background and observations errors are independent, then the covariance matrix, Sv of the inverse solution may be easily obtained by considering the expectation value of the product

t d^v d^v . We have

1
t
1  Sv Z Kt S
1 KCgB
1 K S KCg2 B
1

1 ! Kt S
1 KCgB
1

ð10Þ

which simplifies in the usual expression for the covariance matrix derived by Rodgers (1976) for g Z 1. It is also interesting to note that for g / N, we have Sv Z B and the final solution coincides with the first guess, whereas, for g / 0, the solution tends to be that of unconstrained Least Squares and Sv Z (KtS
1K)
1. Like any regularization scheme, the present one also yields a biased estimation for the true state vector, m or v. We have, in general,
 bZE ^v
vs0 ð11Þ where, again, E( ) denotes expectation value and v denotes the true state vector. The degree of the bias b, for a given inverse solution, may be easily evaluated under the assumptions of zero mean noise affecting the observations. Using Eq. (8) and Eq. (2), we have
 ^
mZ^v
vZ gB
1 CKt S
1 K
1 m   ! Kt S
1 KmCn3
gB
1 q0
m

ð12Þ

Considering expectation values and using E(n3) Z 0, we have


 ^
mZE ^v
vZ gB
1 CKt S
1 K
1 bZE m  ! ½Kt S
1 Km
gB
1 q0
m

ð13Þ

From this last equation we easily derive

 ^ 
1
vE ^v vb vE m Z
IZ
IZ gB
1 CKt S
1 K vm vm vv ! Kt S
1 K
I

ð14Þ

and, finally, we get the Averaging Kernel

A. Carissimo et al. / Environmental Modelling & Software 20 (2005) 1111–1126



 ^ 
1 vE m vE ^v

1 Z Z gB CKt S
1 K Kt S
1 K vm vv

ð15Þ

The formalism of averaging kernel, or resolving power, was first introduced in Backus and Gilbert (1968) and then revisited by Rodgers (2000). It is now a common practice to use the above derivative matrix as a measure of the spatial vertical resolution of the retrieval. Our derivation shows that Eq. (15) gives how the bias propagates from the nearby layers or grid points to a given layer. Its use as a quantitative definition of spatial resolution is not immune to criticism since it addresses only one side of the problem, that is the bias, whereas completely ignores the covariance of the retrieval, Eq. (10). As an example, in the case g Z 0 we have that the averaging kernel is the identity matrix. This just expresses the well known unbiasedness property of unconstrained least squares (Kendall and Stuart, 1979). If we interpret this result in terms of spatial resolution, then this last one would coincide with the pressure or altitude grid used for representing the state vector, v, and we would have achieved the highest spatial resolution with the given data. However, this conclusion, which may be read in many averaging kernel papers, is inconsistent since the covariance matrix has, in general, non-null off-diagonal elements even for g Z 0. The retrieval is, hence, correlated and the state vector has not been independently resolved by the data set. Taking the derivative of the bias vector with respect to a given potential interfering parameter, say X
X0, with X0 the value assumed in the forward calculation against its true value X, we arrive at another useful bias operator. It may be shown, e.g. (Serio et al., 2000), that

 ^ 
1 vE m vE ^v

1 Z Z gB CKt S
1 K Kt S
1 KX ð16Þ vX vX where KX is the Jacobian matrix of the parameter X computed at X0. The above derivative can be used to assess the bias introduced, e.g., by a given atmospheric gas which is not included in the retrieved state vector and whose actual value may deviate from that we use in the forward calculation. It is interesting to note that the bias introduced by the inverse algorithm (Eq. (13)) tends to zero for g / 0, whereas the bias introduced by a given interfering factor (Eq. (16)) goes to zero for g / CN. These are two conflicting conditions, which, once again, show the limitation of defining the spatial resolution through the bias affecting the final solution. 2.2.1. Iterization For mildly non-linear problems, the scheme above may be easily iterated. Towards this objective, we use the basic form of the solution given by Eq. (9). We get


^nC1 Z gn B
1 CKt S
1 Kn
1 Kt S
1 yn x n n

1115

ð17Þ

where the subscript n Z 0, 1, 2,. indicates quantities at the iteration n. The procedure is initialized by choosing v0 and setting R0 ZFðv0 Þ  vFðvÞ K0 Z vv vZv0
 y0 ZðR
R0 ÞCK0 v0
vg The points obtained by Eq. (17) can be considered as successive approximate solutions of the minimum of the quadratic form given by Eq. (4), obtained by the Gauss– Newton method. At each step the non-linear function F(v) is substituted with its Taylor expansion truncated at the first order term. When implementing in practice the above iterative scheme, a suitable criterion is needed to stop the iterations. An obvious criterion is to look directly at the c2 constraint (the second of Eq. (3)). The iterative process is stopped when the computed value of c2 is below the threshold c2a . In practice the c2 constraint may not be met because of unknown systematic effects which are likely to characterize the forward model, therefore the iterations are stopped based on the following three-level constraint, implemented according to an IF structure with an OR logic c2n %c2a c2n
c2n
1 %0:1 c2n
1 c2n Rc2n
1 where c2n is the calculated value of the c2 variable at the iteration step of order n. A normal run should reach the end through the first constraint or second constraint above. The exit through the second constraint is a clue of some inconsistency between forward model and observations, then the computed c2-value gives a quantitative measure of such a misfit. When the cycle stops because of the third constraint, it is likely that the iterative Newtonian mechanism just fails. In our experience, the exit through the third constraint is a rare event. Once the iterative process has been stopped, and assuming that the forward model is linearizable at the solution point, which is always the case for clear sky, then the covariance matrix is still approximated by Eq.

1116

A. Carissimo et al. / Environmental Modelling & Software 20 (2005) 1111–1126

(10) with the Jacobian matrix computed at the solution point. It is important to stress that the parameter gn can depend on the order of the iteration and it has to be properly chosen at each iteration step. 2.2.2. Efficient computation An optimization stage for the regularization parameter involves to re-evaluate Eq. (17) for a range of values of gn. For this reason we need to rewrite it in a more efficient form. An efficient algorithm can be found by decomposing the solution through Singular Value Decomposition (SVD, see, e.g., Golub and Reinsch (1970)). In the present section we first reduce the problem (17) to a classical ridge regression one and then we use SVD to get an explicit equivalent expression for Eq. (17). Towards this objective, we have first to consider the square root of a covariance matrix. This is defined by B1/2Bt/2 Z B. In principle B1/2 could be obtained by Choleski factorization, since B is symmetric and positive definite. However, this method is not recommended here since B itself may be ill-conditioned. A more accurate computation is obtained by resorting to the SVD theorem. The details of the calculations are not shown here, the interested reader is referred to Tarantola (1987). Assuming for now that we know how to numerically compute the square root of B, formal manipulation of Eq. (17) leads us to   t 1 ^nC1 ZJt zn gn B
2 B
2 CJtn Jn x ð18Þ n where Jn ZS
1=2 Kn and zn ZS
1=2 yn . Eq. (18) can be still rearranged to yield   1 t t 1 ^nC1 ZJt zn B
2 gn ICB2 Jtn Jn B2 B
2 x ð19Þ n ^nC1 , Eq. (17) Defining Gn ZJn B1=2 and ^ unC1 ZB
1=2 x assumes the form of a ridge regression unC1 ZGtn zn : ðgn ICGtn Gn Þ^

ð20Þ

Eq. (20), which is fully dimensionless, constitutes our major results in the mathematical derivation of the inverse solution. As anticipated in Section 2.2, the above result generalizes the Rodgers’ statistical regularization estimator (Rodgers, 1976) to the class of ridge regression (e.g., Kendall and Stuart, 1979) estimators and allow us to choose the positive scalar g to minimize the root mean square difference between the solution, ^v and the true state vector, v. Furthermore, for a truly non-linear problem, it has the additional important advantage to bend Newton iterates towards the feasibility region of the inverse solution, therefore it improves the conver-

gence of the Newton–Gauss minimization algorithm which is deeply embedded in our implementation of the inversion process (see also last paragraph of this section). Solution (20) may be still manipulated to have a more efficient form for numerical computation. Considering the SVD decomposition of the symmetric matrix Gtn Gn , we have Gtn Gn ZUn D2n Vtn , with Vn Z Un (the matrix Gtn Gn , is symmetric!). In addition Vn and Un are unitary matrixes, i.e. Utn Un ZI). Finally, the matrix D2n is diagonal, and the squared notation stresses that its elements, the singular values, are positively definite. Inserting the SVD decomposition of equation GntGn back into Eq. (20), we have ^unC1 ZVn ðgn ICD2 Þ
1 Vt Gt zn n n n

ð21Þ

Recalling that D2n is diagonal, the inversion of the kernel ðgn ICD2n Þ is trivial:

1 1 2
1 ðgn ICDn Þ Zdiag ; /; ð22Þ gn Cs21 gn Cs2M where the diag operation means: build up the diagonal M-by-M matrix with the M-elements included within the diag brackets; s2i is the ith singular value, i.e., the ith element of the matrix D2n . Finally, let us come back to our original objective of deriving an efficient computational form which, for any given iteration order n, allows us to quickly compute the solution vector. Let us consider Eq. (21) in the light of a loop in which, for a given n, ^unC1 has to be computed for a range of values of gn. It is seen that Eq. (21) requires only one initial decomposition at each iteration n to be computed, after which all computations for different values of gn do not need any additional matrix inversion or a new SVD decomposition. We have only to update
1 the inverse kernel ðgn ICD2n Þ , which according to Eq. (22) just consists of inputting the current value for gn, no additional operation is required for the matrices Vn and D2n ; which need only an ab initio calculation. From a computational point of view, the method is very fast and efficient. As anticipated in Section 2.2, it is also numerically stable in the presence of possible outliers given by very small singular values, which could have the effect of drawing away the solution from the convergence region. This last property is clearly seen by looking at the formal SVD decomposition of the solution kernel given by Eq. (22). Here the damping action of gn over the smallest singular values is mathematically obvious. 2.2.3. Choice of the regularization parameter The present section deals with the problem of the optimal choice of the regularization parameter. A

1117

A. Carissimo et al. / Environmental Modelling & Software 20 (2005) 1111–1126

general method for the choice of gn which applies to any linear or linearized inverse problem is the so-called Lcurve criterion (Hansen, 1992). The method is based on the study of the curve, called L-curve, obtained by plotting the L2-norm of the regularized solution (i.e., ^nC1 ) x ^t B
1 x ^nC1 k2 Zx ^nC1 k B
1=2 x nC1 2




 00 0 gICD2 Vt ^u C2Vt ^u Z0

ð26Þ

The first and second derivatives are then easily obtained by Eq. (25) and Eq. (26), respectively. We get
 ^u0 Z
V gICD2
1 Vt ^u ð27Þ


^u Z
2V gICD2
1 Vt ^u0 00

against the L2-norm of the corresponding residual h

 
i ^nC1
yn k2 Z Kn x ^nC1
yn t S
1 Kn x ^nC1
yn k S
1=2 Kn x 2 for gn ˛ [0, CN]. This implies that the solution vector has to be computed for a series of values for gn. The optimal regularization parameter is then the value of g that corresponds to the point where the L-curve has its maximum curvature (the sharpest corner). In the rest of the present subsection we shall drop the iteration index for simplicity’s sake; it is intended that ^ u is referred at the iteration n C 1 and all other variables at the iteration n. The L-curve is more efficiently derived in terms of the dimensionless solution (Eq. (21)). Let us define aðgÞZ^ ut ^ u

where we stress that ^ u depends on g through Eq. (21). The curvature of the L-curve may be then explicitly expressed in terms of the parameter g,   0 a ðgÞb00 ðgÞ
a00 ðgÞb0 ðgÞ ð23Þ cðgÞZ h  2 i32 2 0 0 ða ðgÞÞ C b ðgÞ where the prime indicates derivative with respect to g. The L-curve criterion, then, consists in choosing the value of g that maximizes the curve c(g). To derive an analytical expression for c(g), we need the first and second derivative of the solution vector, ^u, with respect to g. In fact, we have  0 t

t  0  ^ a0 ðgÞZ ^ u u u C ^ u ^ ð24Þ

which shows that the problem is reduced to compute the first and second derivative of ^ u. To this end, let us consider Eq. (21). It can be directly differentiated with respect to g to give  0
gICD2 Vt ^ u CVt ^ uZ0 ð25Þ which, differentiating a second time, gives


t 0 t
0 b0 ðgÞZðG^u Þ G^u
z C G^u
z ðG^u Þ 00 t

0 t

0

b00 ðgÞZðG^u Þ ðG^u
zÞC2ðG^u Þ ðG^u Þ t 00 CðG^u
zÞ ðG^u Þ 0

t u
zÞ bðgÞZðG^ u
zÞ ðG^

 00 t
  0 t 0
t  00  ^ a00 ðgÞZ ^ u uC ^ u u C2 ^ u ^ u ^

which contain the solution ^u and matrices that have been already computed to calculate ^u itself. a#(g) and 0 00 a$(g) are easily obtained by substituting back ^u and ^ u into Eq. (24). The computation of b#(g) and b$(g) involves the same kind of calculation.

00

with ^u and ^u given by Eq. (27). Note that by virtue of Eq. (23) and of the explicit expression of the solution (21) made possible by SVD, an explicit formula is developed also for the curvature, that makes the search of the maximum over g feasible from a computational point of view. In practice, a global search algorithm has to be used, since several maxima are, in general, present in the L-curve. In order to implement the L-curve criterion, for each iteration order n, the solution vector ^unC1 has to be computed for a range of g-values to obtain the curve c(g) as a function of g. We define g Z 10j, with the integer index, j which takes values within a regularly spaced interval [–n, n]. For a given application n needs to be optimized by trial and error, for our application we use n Z 4. The division of the interval [–n, n] in one hundred parts is typically enough, that is j Z
n C step ! (k
1), with step Z (2n)/99 and k Z 1,.,100.

2.2.4. Choice of the first guess The choice of the first guess is important in order to speed up the iterative method for solving the inverse problem, by reducing the number of iterations. Moreover, it is a well established fact that a poor choice for the first-guess could reduce the capability of the inverse method to find the solution at all. To gain flexibility, the package d-IASI does not include a module for the initialization of the inversion process (first guess state of the atmosphere), which is left to the specific user application. Of the many possible

1118

A. Carissimo et al. / Environmental Modelling & Software 20 (2005) 1111–1126

initialization schemes (not provided within the d-IASI software), we quote here those that have gained some popularity,  Climatology  Linear Least Squares Regression  Neural Networks In the case of Climatology, the initial guess state of the atmosphere is selected on the basis of a library search method. The library is composed of a suitable set of profiles for temperature, water vapor and other atmospheric gas constituents. For each state of the atmosphere defined in the library, a synthetic spectrum is computed through a forward module and stored in the library itself. In practice, the first guess is chosen according to the lowest c2-index between observed and synthetic spectral radiance. The above approach has been taken by Chedin et al. (1985). In addition, the library provides the input profiles to build up the background covariance matrix, B. The Linear Least Squares Regression method gets the first guess using a least squares regression between the radiance vector, R and the state vector v: R Z Av, where the matrix of weights, A, is obtained on the basis of a suitable training data set. The method may attain a high level of sophistication in its mathematical formulation, e.g, it may include an Empirical Orthogonal Function decomposition of the pair (R, v) such as the one considered in Zhou et al. (2002), however the basic strategy is to seek for a linear regression between R and v. Based on the training, or better, on a independent validation data set, it is then possible to compute the performance of the scheme by computing the difference dvZ^v
v, with ^v the linear regression estimation of v. Choosing the first guess by Least Square regression forces us to consider the difference between background and first guess, that is the difference between the vectors vg and v0 (see Section 2). With the first guess obtained by the observations, v0 is no more independent of the observations themselves, a fact that violates basic assumptions needed to derive the inverse solution and its covariance matrix. For this case, vg and its covariance matrix are still obtained by climatology, and v0, that is the regression solution, is just the first guess: the first point we consider as a candidate for the retrieval. Any second guess is then found by linearizing about the first guess point. The Neural Network approach is similar in its logic to the regression method above. The difference is that the relationship between input and output may be non linear. As before, the method needs a training data set to compute the weights of the net, and a validation data set to assess its performance. A neural network approach has been considered by Luchetta et al. (2002, 2003). As for Least Squares regression, it is here important to

distinguish between background and first guess state vector.

3. Application to IASI and IMG 3.1. Simulation of a hybrid statistical-physical IASI retrieval processor To exemplify the use of d-IASI here we discuss the performance of a hybrid statistical-physical scheme to get temperature, water vapor and ozone from IASI spectra. The hybrid scheme consists in retrieving the state vector from a neural net followed just from one single d-IASI iteration initialized with the neural net retrieval. In other words, the exercise mimics a possible operational retrieval processor which obtains the solution from a statistical algorithm and uses one physical iteration to improve the statistical inversion. The neural network used in the exercise has been developed by Luchetta et al. (2002, 2003). This retrieval exercise applies to clear-sky conditions and sea surface. The Masuda model (Masuda et al., 1988) for emissivity has been used. The use of this hybrid method is presently under consideration for IASI. An alternate strategy, under present consideration, too, uses an Empirical Orthogonal Function (EOF) regression instead of the neural net method. IASI has three spectral bands which are conventionally referred to as Bands 1, 2 and 3. The definition of the three bands is (1) Band 1: 645–1210 cm
1; (2) Band 2: 1210–2000 cm
1; (3) Band 3: 2000–2760 cm
1. which, taking into account for the IASI sampling rate of 0.25 cm
1, gives 8461 spectral channels. 3.1.1. Training and test data set The weights of the network have to be properly determined on the basis of specified input/output pairs. To this end a suitable data set consisting of atmospheric temperature, water vapor and ozone profiles has been build up. The data set consists of 2311 diverse atmospheric states for the training and validation phase and 663 diverse profiles for the cross-check or test phase. The data set has been generated from two data sets. 1. The TIGR data set (Chedin et al., 1985) which is widely used by the meteorological remote sensing community as it represents the global variability of the Earth’s atmosphere in terms of temperature, humidity, and ozone.

A. Carissimo et al. / Environmental Modelling & Software 20 (2005) 1111–1126

2. A globally representative set of temperature and moisture profiles, collected by NOAA/NESDIS from radiosonde ascent in the years 1988 and 1989. This data set does not include ozone profiles. Instead, they have been taken from the TIGR data set, according to the closest geographical coordinates and month of the year. The 2311 C 663 states of the atmosphere have no overlap in the sense that temperature and water vapor profiles have no repetition in the global set. The situation is different for ozone. There are only 400 diverse profiles for ozone, because observations of such atmospheric gas are much more limited than those for temperature and water vapor. Therefore the same ozone profile may occur in various single states of the atmosphere, although with different temperature and water vapor content. Furthermore, there is overlap between the training/validation and the test data sets, as far as ozone is concerned. The training of the neural net is performed by means of an early stopping procedure which explains why the full data set is divided into three subsets: a training set, a validation set and finally a test set. The frequency of air mass types occurring in the training/ validation and test data base is shown in Table 1.

3.1.2. Details of the d-IASI inversion strategy The physical inversion has been performed for the 663 test cases shown in Table 1, and as said before, d-IASI was initialized with the neural net output. The inversion strategy we exemplify in this section consists in retrieving skin temperature, temperature and water vapor profiles simultaneously. This is done by exploiting the spectral range 660–830 cm
1, that is the core of the 15 mm CO2 absorption band, and the shortwave side of the CO2 absorption band: 2180–2250 cm
1. If needed, water vapor retrieval is refined by inverting H2O alone in the spectral range 1100–2000 cm
1. The sequential approach for water vapor saves computation time and avoids to use together weak and strong absorption water vapor lines which affect the linearization of the Radiative Transfer equation at a very different degree. In addition, to alleviate the problem of saturation in the water vapor band we use a selection of channels whose spectral radiance is best linearly Table 1 Number of profiles in the training/validation and test data sets as a function of the air mass type Air mass type

Training/validation data set

Test data set

Tropical Mid-latitude summer Mid-latitude winter High-latitude summer High-latitude winter

595 305 388 283 740

221 43 155 80 164

1119

correlated with the water vapor mixing ratio (DeFeis et al., 2003). The information for temperature and water vapor is then used to retrieve sequentially for ozone in the spectral domain 1000–1080 cm
1. The choice above for the spectral ranges could be changed based on the user specification. In principle, dIASI allows for any spectral range to be used for the inversion. Simultaneous or, alternatively, sequential inversion strategies, or a combination of both, may be chosen by the user, as well. 3.2. Results The results are here shown in terms of root means square error between retrieval and test or true state vector, computed directly as the expectation value E(^v
v2 )1/2. This quantity addresses both the effect of bias and variance over the retrieval. It is beyond the aim of this study to provide a comprehensive analysis of retrieval interdependency and vertical resolution of the inverse products. This analysis has been addressed in many studies performed by the IASI science working group and are available at the IASI web site: http:// smsc.cnes.fr/IASI/. Based on the current IASI specification, for temperature the spectrometer should achieve independent retrieval in 1 km layers in the lower troposphere, for water vapor the figure is 2 km layers. Fig. 1 shows the result for the skin temperature, Ts. The histogram of Fig. 1a shows the difference, DTZT^sfg
Tts , between the first guess value, obtained on the basis of the neural network scheme, and the test value. Fig. 1b shows the final difference reached with the d-IASI physical inversion. We see that the neural network scheme reaches very close to the solution, the root mean square difference is 0.67 K. However, the final physical inversion improves the result even more, getting to a root mean square difference of 0.45 K. The number of test cases considered for this exercise is 663 and include all the test cases shown in Table 1. For the case of temperature, water vapor, ozone profiling accuracy we limit ourselves to show the results for a tropical atmosphere and a high-latitude winter air mass. A detailed account of the results for the remaining air mass types is given in Serio et al. (2001). Fig. 2 exemplifies the retrieval performance for temperature, in the case of a tropical air mass. The exercise shown in Fig. 2 includes a total of 221 test cases (see Table 1). The figure shows the root mean square error achieved with the first guess and that reached after the final step of the physical inversion. It is possible to see that the final physical inversion is able to get an accuracy below 1 K in most of the troposphere. For the case shown in Fig. 2, the temperature profile has been retrieved over a pressure layer grid of 44 levels (see, e.g., Amato et al. (2002)).

1120

A. Carissimo et al. / Environmental Modelling & Software 20 (2005) 1111–1126

Number of events

200

a)

150 100 50 0 −4

−3

−2

−1

0

1

2

3

4

Number of events

300

b)

250 200 150 100 50 0 −4

−3

−2

−1

0

1

2

3

4

∆ T (K) Fig. 1. Retrieval performance for skin temperature. Panel a shows the difference DTZT^sfg
Tts , between the first guess value, obtained on the basis of the neural network scheme, and the test value. Panel b shows the final difference reached with the d-IASI, physical, inversion. The number of test cases considered for this exercise is 663.

The analysis for water vapor is shown in Fig. 3. Again the physical inversion improves over the first guess and reaches a root mean square error in the boundary layer that is better than 2 g/kg (water vapor mixing ratio). At the tropics the water vapor concentration in the boundary layer may be of the order of 10–20 g/kg, so that the above values correspond to a performance of 10–20% accuracy in the lower atmosphere.

0

Fig. 4, again for the case of a tropical air mass, shows the expected performance for ozone. The first guess provided by the neural network is accurate within 10% in between 10 and 50 km altitude range. However, the final physical inversion improves even more over the first guess throughout the atmosphere. In considering this result, the limitation of the ozone data set should be kept in mind. It is likely this finding is in part due to the highly redundant information used for the training/ validation and test data set (see Section 3.1.1). As opposite to a tropical air mass, we consider now the coldest and dry model of atmosphere. The results for temperature, water vapor and ozone are shown in Figs.

100

300

0

Neural Net F.G. Physical Retrieval

400 500

200

600

300

700 800 900 1000

0

0.5

1

1.5

2

2.5

Neural Net F.G. Physical Retrieval

100

Pressure (mbar)

Pressure (mbar)

200

3

RMS Temperature Error (K)

400 500 600 700 800 900

Fig. 2. The figure exemplifies the retrieval performance for temperature, in the case of a tropical air mass. The exercise includes a total of 255 test cases. The figure shows the root mean square error achieved with the first guess and that reached after the final step of the physical inversion.

1000 0

0.5

1

1.5

2

RMS H2O Error (g/kg) Fig. 3. As Fig. 2, but for water vapor.

2.5

3

1121

A. Carissimo et al. / Environmental Modelling & Software 20 (2005) 1111–1126 0

10-1

100 Neural Net F.G. Physical Retrieval

200

Pressure (mbar)

Pressure (mbar)

100 Neural Net F.G. Physical Retrieval 10% reference line

101

102

300 400 500 600 700 800 900

103

1000 0

5

10

15

20

25

30

35

40

45

0

50

0.1

0.2

0.3

rms O3 error, (%)

5–7, respectively. For temperature, one more physical inversion step improves significantly the result over that provided by the neural net. For water vapor the improvement mostly limits to the lower part of the atmosphere. For ozone again we have a nice performance around 10%. However, for this last case, as said before, the limitation of the ozone data set has to be considered. Similar results have been obtained for air masses different from the tropical ones. As a rule, the physical inversion improves over the first guess in any condition. 3.3. Application to IMG In this section we consider, for illustrative purposes, examples of temperature, water vapor and ozone retrieval based on IMG observations. The set of IMG spectra we have used in our analysis is shown in Table 2. 0

0.6

0.7

0.8

0.9

1

We have a total of 10 spectra, seven of which have been recorded in the tropical belt, whereas of the remaining three, two may be classified as mid-latitude winter, and one as sub-arctic winter. This last spectrum has been recorded close to Point Barrow, Alaska. All the spectra are for clear sky as shown by co-located Ocean Color and Temperature Scanner (OCTS) imagery. OCTS (Shimada et al., 1999) is an imager which has been a companion instrument of IMG on board the ADEOS/ 1 platform. Apart from the sub-arctic sounding, which has been taken on sea–ice surface, the soundings have been obtained on sea surface. IMG has a nominal sampling rate of 0.05 cm
1. However, our analysis is mostly intended to exemplify the potential performance of IASI. For this purpose, IMG spectra have been transformed into IASI-like spectra. For the generation of synthetic spectral radiance, this is done by convolving the monochromatic spectral radiance with a pure cardinal sine function,

Neural Net F.G. Physical retrieval

200

10-1

300 400

100

500

Pressure (mbar)

Pressure (mbar)

0.5

Fig. 6. As Fig. 5, but for water vapor.

Fig. 4. As Fig. 2, but for ozone.

100

0.4

rms H2O error, (g/kg)

600 700 800 900 1000 0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Neural Net F.G. Physical Retrieval 10% reference line

101

102

rms Temperature Error (K) Fig. 5. The figure exemplifies the retrieval performance for temperature, in the case of a high-latitude winter air mass. The exercise includes a total of 164 test cases. The figure shows the root mean square error achieved with the first guess and that reached after the final step of the physical inversion.

103

0

5

10

15

20

25

30

35

rms O3 error, (%) Fig. 7. As Fig. 5, but for ozone.

40

45

50

1122

A. Carissimo et al. / Environmental Modelling & Software 20 (2005) 1111–1126

Table 2 List of IMG spectra analyzed in this paper IMG record

Obs.

Lat., Lon. (degrees)

Date

166922 166923 166923 166923 166934 167229 167530 327906 357011 374126

5 4 5 6 5 2 1 1 1 2

(33.15 N, 136.19 E) (27.44 N, 134.4 E) (26.68 N, 134.2 E) (25.91 N, 134.01 E) (44.54 S, 116.3 E) (8.16 S, 50.45 E) (12.14 S, 26.10 W) (39.9 N, 8.9 E) (72.01 N, 204 E) (11.92 N, 52.37 E)

12 December 1996 12 December 1996 12 December 1996 12 December 1996 12 December 1996 12 December 1996 12 December 1996 3 April 1997 24 April 1997 6 May 1997

whose parameters are set up to correspond to the instrumental line shape of an ideal interferometer spectrometer of maximum optical delay of 2 cm, and, finally, apodized with a Gaussian function of 0.25 cm
1 half-width half-height. For consistency IMG spectra (real observations) have been apodized with the same function as that used for the synthetic spectral radiance. The sampling rate of the final spectra (simulated and observed) is Ds Z 0.25 cm
1. We choose Gaussian apodization because it is consistent with what for now has been prescribed for the IASI level 1C data, so that our calculations yield IASI-like spectra. For the present exercise, because of limitation of IMG, only the spectral range 667–1600 cm
1 has been used. The strategy for the inversion process of IMG spectra is:  retrieve simultaneously skin temperature, temperature and water vapor profiles by exploiting the spectral range 667–830 cm
1;  improve, if needed, the H2O retrieval by considering a sequential inversion step which exploits a set of channels in the spectral range 1100–1600 cm
1; in this step water vapor is retrieved alone;  finally use temperature and H2O retrieved above to invert for ozone in the spectral band 1000– 1080 cm
1. The inversion process will be illustrated step-by-step by considering the IMG spectrum # 166923, Obs. 4, recoded in the tropics. For all the cases, illustrated in this section, of IMG inversions, the initial state of atmosphere has been obtained by climatology. The simultaneous step for temperature and water vapor takes three iterations to converge. Convergence is here checked through the c2 constraint (the second of Eq. (3)). According to this constraint and considering the normalized c2 index, that is c2o Zc2 =N, with N the number of data points, wephave ffiffiffiffiffiffiffiffiffi that the iterations are stopped when c2o !1C2 2=N  1. Note that the number of data points, N, is in any case of a magnitude pffiffiffiffiffiffiffiffiffi as large as to make the additional term 2 2=N

negligible with respect to 1. The sequence of c2o values for the order of iteration 0, 1, 2 and 3 is 51.2, 12.3, 1.4, and 1.00, respectively. Note that the 0th iteration corresponds to the first guess. The fact that the normalized c2o index goes below its critical value (that in this case is 1.11 because N Z 653) means that our forward model is consistent with the observations within the IMG radiometric noise. Of course, this does not necessarily imply a good quality of retrieval. The quality of the inverse results may be directly inspected by Fig. 8 which compares temperature to the ECMWF analysis for the same location and date. The ECMWF analysis is of poor spatial resolution if compared to the 8 km2 IMG Field of View. In addition, our retrieval is obtained on a 43 pressure layer grid, whereas the ECMWF analysis is only available at 15 pressure levels, so that the comparison has to be taken with some caution. With this in mind, Fig. 8 shows that the retrieval is highly consistent with the ECMWF analysis. The comparison for water vapor is shown in Fig. 9 (left-hand-side panel) and, again, we see a fair agreement with the ECMWF analysis. Finally, a good agreement is found for the skin temperature, Ts. We have 296.6 K (ECMWF analysis) against the 296.2 K retrieved. triplet of thermodynamic parameters 
The Ts ; T; qH2 O , that is, skin temperature, temperature profile and water vapor profile, feeds now the sequential step for water vapor alone. The inversion here considers a set of N Z 1375 channels, selected according to (DeFeis et al., 2003) from
the spectral  range 1100– 1600 cm
1. Using the triplet Ts ; T; qH2 O obtained from the simultaneous stage, we first compute the synthetic spectrum for this new range and then we calculate the related c2, which, for the case at hand, yields the value 3.1. This is higher than the critical value (1.08). We assume here that the discrepancy may be due mostly to water vapor, therefore we try to improve the final product for H2O by considering a new cycle of iterations for H2O alone. This cycle converges in one iteration since we have c2o Z1:05 after the first iteration. The high convergence rate says that the final solution is close to that found at the previous step, as it is possible to see by looking at the results shown in Fig. 9 (right-hand-side panel). We are now at the last inversion for ozone which uses the spectral range 1000–1080 cm
1 (N Z 321). The c2 drops by more than one order of magnitude (from 25 to 0.86) in a single iteration. This demonstrates the efficiency of our scheme. The results for ozone are shown in Fig. 10. For all the IMG spectra here analyzed, we observed a rapid convergence, in terms of the c2-criterion, even though the inversion process was initialized, as said before, by climatology. The first two inversion stages (for ðTs ; T; qH2 O Þ and qH2 O alone, respectively) took 2–3

1123

A. Carissimo et al. / Environmental Modelling & Software 20 (2005) 1111–1126 0 10-1

100 200

100

Pressure (mbar)

Pressure (mbar)

300 400 500

F.G. it1 it2 it3 ecmwf

600

101

700 102

800 900 1000 180

200

220

240

260

280

103 180

300

Temperature (K)

200

220

240

260

280

300

Temperature (K)

Fig. 8. Illustrating the d-IASI inversion process for temperature. The case shown applies to the IMG spectrum # 166923, Obs. 4. The left-hand-side panel shows the first guess, the first, second and third iteration (final iteration). For comparison the ECMWF analysis (up to 10 mbar) for the same location and date is shown, as well. The right-hand-side panel shows the same drawing as the left-hand-side panel but in logarithmic scale for the pressure-axis.

Simultaneous Step

Sequential Step

0

0 F.G. it1 it2 it3 ecmwf

100 200

100 200

300

300

Pressure (mbar)

Pressure (mbar)

F.G. Final Product ecmwf

400 500 600

400 500 600

700

700

800

800

900

900

1000

1000 0

5

10

H2O mixing ratio (g/kg)

15

0

5

10

15

H2O mixing ratio (g/kg)

Fig. 9. Illustrating the d-IASI inversion process for water vapor. As for Fig. 8, the case shown applies to the IMG spectrum # 166923, Obs. 4. The left-hand-side drawing shows the three iterations which are taken simultaneously with temperature, whereas the right-hand-side shows the sequential inversion step for water vapor alone, see text for details. For comparison the ECMWF analysis (up to 10 mbar) for the same location and date is shown, as well.

1124

F.G. it1

Pressure (mbar)

10-1

100

101

102

103

0

5

10

15

O3 mixing ratio (ppm)

and 1 iterations, respectively, to converge. For ozone, in three cases we needed two iterations to get to convergence. For all the cases we have examined, the retrieval for temperature and water vapor showed a good agreement with the ECMWF analysis. Figs. 11 and 12 allow us to have an overall check of the quality of the inverse results. Fig. 11 compares the IMG retrieval for skin temperature to the ECMWF analysis (the figure is limited to the nine sea-surface IMG soundings). The root mean square difference is within 0.54 K. Fig. 12 is an interesting check for ozone. It compares the columnar ozone, derived from IMG, to that obtained by the TOMS/ADEOS instrument for the same location and date. We see that the comparison evidences no important bias. The root mean square difference is 6 Dobson, which is less than 2% on average.

Skin Temperature (IMG retrieval), (K)

305

300

295

290

285

285

290

295

450

400

350

300

250

200 200

250

300

350

400

450

Columnar Ozone (TOMS/ADEOS), (Dobson units)

Fig. 10. Illustrating the d-IASI inversion process for ozone. As for Fig. 8, the case shown applies to the IMG spectrum # 166923, Obs. 4. The figure shows the first (and final) iteration along with the first guess for the final sequential step of d-IASI which retrieves ozone.

280 280

Columnar Ozone (IMG retrieval), (Dobson units)

A. Carissimo et al. / Environmental Modelling & Software 20 (2005) 1111–1126

300

305

Skin Temperature (ECMWF Analysis), (K) Fig. 11. Checking the accuracy of the retrieved skin temperature. The figure shows the scatterplot of the skin temperature from the ECMWF analysis and the skin temperature retrieved from IMG observations.

Fig. 12. Checking the accuracy of the retrieved columnar ozone. The figure shows the scatterplot of the columnar ozone retrieved from TOMS/ADEOS and that retrieved from the IMG observations.

4. Conclusions The paper has described and discussed a physical inversion strategy to be used for the inversion of thermal up welling spectral radiance. The scheme is mostly intended for IASI, in that it matches the spectral coverage of this spectrometer, however, the mathematics is quite general and the method can be applied to the broad class of infrared spectrometers with a sampling rate between 0.1 and 2 cm
1. We have developed a novel approach in which the usual statistically regularized estimator (Rodgers’ solution to inverse problems, e.g., Rodgers (1976)) is generalized to the class of ridge regression estimators which allows us to effectively optimize the root mean square error between the solution and the true state vector and also has the advantage to control the convergence of the Gauss– Newton minimization algorithm. The implementation of the method makes full exploitation of the Singular Value Decomposition theorem, which yields a very computationally efficient algorithm. The d-IASI method is mostly intended for research applications rather than for operational end-uses. Although, for operational uses, the method could back up, e.g., statistical inversion approaches, such as least square regression and/or neural network. Under likely weather conditions, these approaches may fail to reach the goal accuracy because the training data set they need could well be lacking completeness as far as the variety of weather conditions is concerned. A sequential physical inversion, limited, to one iteration, could improve the final products. It should be considered that in the present version, the inverse module takes a few milliseconds per iteration to run against 4–6 min per iteration taken by the forward module (1.0 GHz Pentium III platform, with Compaq Visual Fortran 6

A. Carissimo et al. / Environmental Modelling & Software 20 (2005) 1111–1126

compiler). This could be improved by a factor of about 2–3 by allowing for a parallel computation of the three Jacobian matrixes for temperature, H2O and O3, respectively, which in the present version is simply sequential. The method has been illustrated with the help of two examples: (1) simulation of a hybrid statistical-physical IASI retrieval processor; (2) retrieval exercises for the IMG spectrometer. For IASI we have shown that the spectrometer easily achieves the required accuracy of 1 K in 1 km layer in the troposphere for temperature. For water vapor the final accuracy lies between 10 and 20% depending on the atmospheric load of water vapor. For cold and dry atmospheres, where the water vapor load may go below 1 g/kg, IASI may be highly inaccurate. Finally, good results are found for ozone, although our findings may be optimistically biased because of the ozone climatology we have used to define the training/validation and test data set. IMG retrievals have been considered for different air masses, encompassing tropical and sub-arctic winter situations. The results have been compared to ECMWF analysis and a fair agreement has been found. The comparison of columnar retrieved ozone to ADEOS/ TOMS products shows a fine agreement and evidences the potential capability of infrared sensors, such as IASI, to get high quality retrieval of such a gas which is of a major climatological concern.

Acknowledgments Work supported by Italian Space Agency and EUMETSAT (contract EUM/CO/99/688/DD).

References Amato, U., Masiello, G., Serio, C., Viggiano, M., 2002. The s-IASI code. Environmental Software & Modeling 17, 651–667. Amato, U., Serio, C., 2000. Remote sensing of the atmosphere from high spectral resolution infrared radiation. Recent Research Developments in Geophysics 3 (Part I), 17–36. Amato, U., Cuomo, V., DeFeis, I., Romano, F., Serio, C., Kobayashi, H., 1999. Inverting for geophysical parameters from IMG radiances. IEEE Transactions on Geoscience and Remote Sensing 37 (3), 1620–1632. Amato, U., DeFeis, I., Serio, C., 1996. Linearization pseudo-noise and its effect on the retrieval of atmospheric state from infrared spectral radiances. Geophysical Research Letters 23 (18), 2565–2568. Backus, G., Gilbert, F., 1968. The resolving power of gross Earth data. Geophysical Journal of the Royal Astronomical Society 16, 169–205. Cayla, F.-R., 1995. IASI Infrared Interferometer for Operations and Research. In: Che´din, A., Chahine, M.T., Scott, N.A. (Eds.), High Spectral Resolution Infrared Remote Sensing for Earth’s Weather and Climate Studies. In: NATO ASI Series I9. Springer Verlag, Berlin, pp. 9–19. Chedin, A., Scott, N.A., Wahiche, C., Moulnier, P., 1985. The improved initialization inversion method: a high resolution method

1125

for temperature retrievals from satellites of the TIROS-N series. Journal of Climate and Applied Meteorology 24, 128–143. Clough, S.A., Iacono, M.J., Moncet, J.-L., 1992. Line-byline calculation of atmospheric fluxes and cooling rates, 1, Application to water vapor. Journal of Geophysical Research 97, 15761–15785. Clough, S.A., Iacono, M.J., 1995. Line-by-line calculation of atmospheric fluxes and cooling rates II: application to carbondioxide, ozone, methane, nitrous-oxide and the halocarbons. Journal of Geophysical Research 100, 16519–16535. Clough, S.A., Kneizys, F.X., Davies, R.W., 1989. Line shape and the water vapor continuum. Atmospheric Research 23, 229–241. Colton, D., Kress, R., 1992. Inverse Acoustic and Electromagnetic Scattering Theory. Springer-Verlag, New York. DeFeis, I., Lubrano, A.M., Serio, C., 2003. Channel selection for water vapor retrieval. SPIE 4882, 363–374. Golub, G.H., Reinsch, C., 1970. Singular value decomposition and least squares solutions. Numerische Mathematik 14, 403–420. Hansen, C., 1992. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Review 34 (4), 561–580. Kendall, M., Stuart, A., 1979. Inference and Relationship. The Advanced Theory of Statistics, vol. 2. Charles Griffin, London & High Wycombe. Kobayashi, H., Shimota, A., Yoshigahara, C., Yoshida, I., Uehara, Y., Kondo, K., 1999. Satellite-borne high-resolution FTIR for lower atmosphere sounding and its evaluation. IEEE Transactions on Geoscience and Remote Sensing 37 (3), 1496–1507. Lubrano, A.M., Serio, C., Clough, S.A., Kobayashi, H., 2000. Simultaneous inversion for temperature and water vapor from IMG radiances. Geophysical Research Letters 27, 2533–2536. Lubrano, A.M., Masiello, G., Serio, C., Matricardi, M., Rizzi, R., 2002. IMG evidence of chlorofluorocarbon absorption in the atmospheric window region 800–900 cm
1. Journal of Quantitative Spectroscopy & Radiative Transfer 72, 623–635. Luchetta, A., Schlu¨ssel, P., Serio, C., 2002. Neural network approach to the inversion of high resolution IASI spectra in troposphere and stratosphere. Atti della Fondazione Giorgio Ronchi Anno LVII (4), 731–736. Luchetta, A., Serio, C., Viggiano, M., 2003: Neural Network to Retrieve Atmospheric Parameters from Infrared High Resolutionb Sensor Spectra, Technical Proceedings of the 2003 IEEE International Symposium on Circuits and Systems, IEEE Catalog Number: 03CH37430C, ISBN: 0-7803-7762-1 (CD-ROM), vol. V, pp. 745–748. Masiello, G., Matricardi, M., Rizzi, R., Serio, C., 2002. Homomorphism between clear-sky and cloudy spectra in the 800–900 cm
1 window region. Applied Optics 41/6, 965–973. Matricardi, M., Saunders, R., 1999. Fast radiative transfer model for simulation of infrared atmospheric sounding interfererometer radiances. Applied Optics 38 (27), 5679–5691. Masuda, K., Takashima, T., Takayama, Y., 1988. Emissivity of pure and sea waters for the model sea surface in the infrared window regions. Remote Sensing of Environment 24, 313–329. Marquardt, D.W., 1963. An algorithm for least squares estimation of nonlinear parameters. SIAM Journal 11, 431–441. Rothman, L.S., Barbe, A., Chris Benner, D., Brown, L.R., CamyPeyret, C., Carleer, M.R., Chancea, K., Clerbaux, C., Dana, V., Devi, V.M., Fayt, A., Flaud, J.-M., Gamache, R.R., Goldman, A., Jacquemart, D., Jucks, K.W., Lafferty, W.J., Mandin, J.-Y., Massie, S.T., Nemtchinov, V., Newnham, D.A., Perrin, A., Rinsland, C.P., Schroeder, J., Smith, K.M., Smith, M.A.H., Tang, K., Toth, R.A., Vander Auwera, J., Varanasi, P., Yoshino, K., 2003. The HITRAN molecular spectroscopic database: edition of 2000 including updates through 2001. JQSRT 82, 5–44. Rodgers, C.D., 1976. Retrieval of atmospheric temperature and composition from remote measurements of thermal radiation. Review of Geophysics and Space Physics 14 (4), 609–624.

1126

A. Carissimo et al. / Environmental Modelling & Software 20 (2005) 1111–1126

Rodgers, C.D., 2000. Inverse Methods for Atmospheric Sounding – Theory and Practice. Series on Atmospheric, Oceanic and Planetary Physics, vol. 2. World Scientific, Singapore. Serio, C., Lubrano, A.M., Masiello G., 2000. CO2 retrieval feasibility study for IASI, EUMETSAT technical report #3, EUM/CO/99/ 688/DD, EUMETSAT, Darmstadt, Germany. Serio, C., Amato, U., Cuomo, V., DeFeis, I., Lubrano, A.M., Masiello, G., Viggiano, M., 2001. Performance evaluation of d-IASI. EUMETSAT Technical Report # 6, EUM/CO/99/688/ DD, EUMETSAT, Darmstadt, Germany. Shimada, M., Oauku, H., Mitomi, Y., Murukami, H., Kawamua, H., 1999. Calibration of the ocean color and temperature scanner.

IEEE Transactions on Geoscience and Remote Sensing 37 (3), 1484–1495. Twomey, S., 1977. Introduction to the Mathematics of Inversion in Remote Sensing and Indirect Measurements. Elsevier, Amsterdam. Tarantola, A., 1987. Inverse Problem Theory. Elsevier, New York. Zhou, D.K., Smith, W.L., Li, J., Howell, H.B., Cantwell, G.W., Larar, A.M., Knuteson, R.O., Tobin, D.C., Revercomb, H.E., Bingham, G.E., Tsou, J.-.J., Mango, S.A., 2002. Thermodynamic product retrieval methodology for NAST-I and validation. Applied Optics 41 (33), 6957–6970.