Variational cloud decontamination and homogeneous component decomposition of IASI measurements

Pascal PRUNET, Bernard TOURNIER, Laure CHAUMAT, Dorothée COPPENS, Carsten STANDFUSS.
NOVELTIS, Toulouse, FRANCE.

Abstract

The aim of this communication is to present the development of a variational cloud-decontamination scheme for the IASI measurements, and its extension to the decomposition of the IASI spectrum. This work was made in the context of the EUMETSAT Polar System (EPS) Programme. We use a one-dimensional variational framework (1DVAR) for the treatment of cloudy data by cloud decontamination. The extraction of clear component spectrum is made through comparison of IASI partly cloudy adjacent pixels. The 1DVAR approach simultaneously extracts cloud-clearing parameters and information about the atmospheric and surface state from infrared measurements. This IASI cloud decontamination scheme allows for complex cloud structures including multiple cloud layers with wavelength dependent radiative properties.

The extension to the IASI measurement decomposition uses the level 1c IASI product of AVHRR radiances analysis in the IASI field of view. This product provides the classification and localisation of the different components inside the IASI footprint.

The variational IASI cloud-decontamination scheme is described. The capability of the scheme to extract cloudy components is highlighted. Its applications for IASI data cloud decontamination, and for the spectral correction of IASI measurement errors due to the heterogeneities of the radiative surfaces inside the footprint is presented.

Introduction

The cloud decontamination refers to the procedure of removing the effect of clouds from infrared measurements. This process applies to partially cloudy measured scenes, and consists in the extraction of the radiance emitted from the cloud-free portion of the instrument footprint. The estimation of the clear-sky component of the infrared measurement is of great importance for use in meteorological operational weather forecast systems. The better correction of cloud-contaminated observations for sounders such as IASI is essential for NWP, as they can represent about 75% of the measurements provided by infrared sounders.

The cloud decontamination can be generalised to the decomposition of the measured radiance into several components, each of them corresponding to a physical surface involved in the instrument pixel or field-of-view (FOV), such as clouds, sea surfaces or land surfaces. The estimation of emission spectra of cloud types or various soil types within the FOV has relevance for various research applications. The understanding of cloud scattering effects and others cloud radiation interactions, or the geographical surface emissivity variability at fine scales, are major challenges for atmospheric and environmental applications.

This paper presents scientific aspects of the developments and the validations of a scheme for the cloud decontamination and component decomposition of the IASI measurements. The work was carried out in the frame of the development of the EUMETSAT polar System (EPS) ground-segment processing chain. NOVELTIS was in charge to provide EUMETSAT with a prototype software for IASI cloud decontamination.

This communication presents first the principle of the cloud-decontamination and the component decomposition of the measured spectrum. The variational IASI decomposition scheme and the synthetic database of IASI observations in cloudy conditions, generated for validation purposes, are described. The capability of the scheme to extract the clear-sky component is analysed. The problem of spectral correction of IASI measurement errors due to the heterogeneities of the radiative surfaces inside the FOV is addressed. The potential of the variational scheme to consistently provides retrieved atmospheric and surface parameter for homogeneous surfaces is illustrated.

Principle

Eyre and Watts (1987) provide a detailed review of cloud decontamination methods. The work presented here is derived from the cloud decontamination algorithm developped by Joiner and Rokke (1997) for the TOVS data. It is an extension of the adjacent FOV approach proposed by Chahine (1974, 1977).

Scene radiance decomposition

The radiance in a partly cloudy FOV is the linear average of clear-sky and cloudy radiances, weighted by their corresponding fractional areas. More generally, the observed radiance Rik in channel i and FOV k, in a scene with J different homogeneous surfaces (scene types), can be expressed with the scene type fraction jk as

Rik = jkRi,j (1)

where Ri,j is the radiance corresponding to the homogeneous component j.

For the scene decomposition, two or more adjacent IASI FOVs are required, where each FOV produces a different realisation of (1). It is assumed that the Ri,j are constant in adjacent FOVs but that the jk vary.

The radiance of component j for channel i can be written as a linear combination of the measured radiances in K adjacent FOVs

Rji - Rip = jk(Rip - Ri,k )(2)

Equation (2) is adapted from Chahine (1977) and Joiner and Rokke (1997) by introducing the pivot point radiance at channel i, Rip, i.e., the average of the adjacent FOV measured radiances at channel i. The jk, hereafter named the decomposition parameters, are channel-independent constants since they depends only on the jk. Equation (1) can be rewritten

Rik - Rip = -jk (Rip - Ri,j)(3)

The equation system (3) consists of K equations, each of them associated to a FOV. In general, as many FOVs are required as homogeneous components are present.The associated matrix equation can be written, for each channel i.

(Rik - Rip)=  (Rip - Rij )(4)

where (Rik - Rip)is the measurement vector of dimension K, (Rip - Rij ) is the component vector of dimension J, and  is the matrix of jk. By comparison with equation system (2) of J equations, the inversion of (4) give an estimate of the jk parameters, and each set of the jk for the jth homogeneous component is the jth line of the inverse matrix -1.

The unknowns of (2) are the radiance components Rij , which are obtained through the estimation of the jk parameters. The available information is the radiance measurements Rik in the K adjacent FOVs, and for some instrumental systems such as IASI, an estimation of the -jk. The following sections propose a method to derive the jk and finally the Rij from this information.

Variational approach

The general approach to solve jk for a given j involves the following steps:

  1. provide a guess estimate Rguessij of the radiance component Rij, for a subset of infrared channels;
  2. Use this guess estimate to solve Equation (2) on this subset of channels to obtain an estimate of the jk;
  3. reconstruct the radiance component for all infrared channels using Equation (2);
  4. retrieve the atmospheric state above the homogeneous surface and associated surface properties using the reconstructed radiance component and independent information on the atmospheric and surface state;
  5. repeat steps 1-4 using a new guess radiance component Rguessij computed from the retrieved atmospheric and surface state.

Following the method of Joiner and Rokke (1997), we use a 1 dimensional variational (1DVAR) technique to simultaneously estimate the decomposition parameters and the atmospheric and surface parameters. In our scheme, we define the state vector as the vector x containing the atmospheric and surface parameters (e.g., temperature and humidity at different layers above the surface, surface pressure, temperature and emissivity) and the vector j containing the K decomposition parameters associated with component j. In practice, we start from a guess estimate xg of x, and the associated error covariance matrix B. We also suppose a guess estimate jg of j and the associated guess error vector g.

1. Rguessij is computed from xgthrough a radiative transfer model (RTM) for a subset of channels. This estimate is associated with the observation operator error standard deviation ci, which is channel dependent. No inter-channel correlation is supposed. In practice, Rguessij is useful in case of a clear sky component j. In case of cloudy component, this estimate can be strongly biased, and no realistic estimation of ci is available.

2. A least-square inversion of (2) with Rij = Rguessij is performed to find an estimate ja of j and the associated guess error vector a. This set of equation is solved by considering a diagonal observation error

(oi)2 = (di)2 + (ci)2(5)

where di is the instrument noise standard deviation at channel i.

3. In parallel with the radiance component Raij computation from the estimate ja, we compute the channel independent decomposition noise amplification factor of the measurement noise due to cloud decontamination procedure.

C = (1 + k)2 + k2(6)

4. Raij and its amplified measurement error standard deviation C×di is used, together with the guess estimate xg and the guess error covariance matrix B, to retrieve an estimate xa of the atmospheric and surface state through 1DVAR inversion.

5. The above scheme is iterated with xa as input of the RTM to compute an improved Rguessij.

The 1DVAR approach simultaneously (i.e., at the same iteration step) and consistently extracts cloud-clearing parameters and information about the atmospheric and surface state from infrared measurements.

If no guess information on the cloud decontamination parameters is available, jg is set equal to zero and its associated error g is set to 108, a very large value. In this hypothesis, the quality of the radiance decomposition is strongly dependent on the quality of the guess estimates xg and Rguessij. The approach gives good results in case of a clear sky component j (e.g., Joiner and Rokke, 1997 for TOVS) but fails for cloudy components.

Exploitation of the IASI sub-pixel information

In the case of IASI, a sub-pixel information is available which can be used to provide an accurate guess estimate of the decomposition parameters.

The AVHRR radiance analysis in the IASI field of view will be available for the level 1c IASI product (Tournier et al., 2002). This product will provide, at the full resolution AVHRR measurement scale (1 km pixel size), the classification and localisation of the different homogeneous components inside the IASI FOV. Each component is characterised by its gravity centre angular position, jk, and its relative fraction coverage jk. Furthermore, this information will be treated by an adapted version of the CMS/Météo-France classification algorithm, providing a characterisation of each component in terms of physical surface type.

The elements jk of  (Equation 4) can be given directly by the AVHRR radiance analysis.  then gives an estimate of the decomposition parameters jk for each of all selected channels. The decomposition parameter vector, channel independent, jg for each component j is provided by their mean values, and the associated error g by their dispersion.

Considering the high quality of the fine scale AVHRR analysis in the IASI FOV, the accuracy of the jk leads to an accurate guess estimate jg. This guess estimate can be introduced in the variational decomposition scheme described above. In this case, the approach uses both the FOV contrast measurement and the AVHRR radiance analysis to give the best estimate of the decomposition parameters. Due to the good accuracy of the jg, the accurate estimate of the homogeneous component radianceRi,j can be obtained without any iteration. in particular the decomposition of the IASI measurement becomes independent from the Rguessij.

Development of the decomposition scheme

The implementation of the radiance component decomposition algorithm, concerns IASI level 1c data. These data consist of calibrated and apodised spectra for each adjacent FOV, and the output parameters of the AVHRR radiance analysis, i.e., jk for each component j and each adjacent FOV k.The organisation of adjacent FOVs is adapted to the characteristics of the IASI measurements. We define a minimum array (or field-of-regard, FOR) of K=4 adjacent FOVs to run the decomposition algorithm. The FOR is a matrix of 2*2 IASI pixels corresponding to one scan position, covering 50*50 km2 at the nadir. It allows to retrieve radiance component coming from a maximum of 4 non homogeneous surfaces (clear sky and three cloud formations for example). However, the number of FOVs in the FOR is parameterised, and can be increased for example to 8, 12 or 16.

The IASI variational decomposition scheme implementation involves several steps including pre-processing and quality control. Our approach considers many aspects of the work of Joiner and Rokke (1997). The main steps of the implementation are described below.

In the preprocessing phase, the user may specify a channel selection to which the algorithm is applied. Component decomposition should only be applied to a channel selection which depends on the particular application. This is because component decomposition is an extrapolation process that amplifies the measurement error. We attempt to determine whether or not some channels are useful in the decomposition process. In this work, the IASI super channel selection is used (Serio et al., 2000) which is configured for temperature, humidity and ozone atmospheric and surface parameters retrievals. It considers 1624 lead channels instead of the 8461 individual IASI channels. The computation of the average of the FOVs radiances in the FOR, the so-called pivot point radiance Rip at channel i, is also computed in the preprocessing phase, for the given channel selection.

The whole state vector is defined at the centre of the FOR. It is composed of the atmospheric and surface state vector x and of the decomposition parameter vector j. The x vector consists of nt layer-mean values of temperature from the surface to 0.1 hPa, nh layer-mean values of the specific humidity from the surface to 0.1 hPa, and the surface skin temperature. The j vector contains the K decomposition parameters jk for the component j. nt, nh, and K can be parameterised, and it is straight-forward to modify the state vector to accommodate additional variables such as ozone profile or spectral infrared surface emissivity.

The input observation vector consists of the level 1c IASI data for each FOV constituting the FOR, i.e., the K calibrated and apodised IASI spectra, and the matrix from the AVHRR radiances analysis associated to the J components of the FOR.

The observation error covariance on the IASI spectra is given by Equation (5), and is assumed to be diagonal (i.e., no inter-channel correlations). d is provided by the specification of the IASI radiometric noise. c is given by the analysis of the differences between spectra generated by the 4A-Model (used for the generation of the synthetic data set) and by the RTIASI model (used in the decomposition scheme for direct and inverse RTM computation).

The amplification of the measurement error by radiance component decomposition is taken into account by multiplication with the amplification factor C given by Equation (6). Note that the larger the values of the decomposition parameters , the larger the amplification effect will be.

The background atmospheric and surface state vector xg is generated from the 3-D atmospheric fields and surface fields used for the simulation of the synthetic IASI data (i.e., the true state vector) by adding noise according the background error covariance matrix B.

A statistical (state independent) forecast error covariance matrix B is considerd as background error. The forecast error covariance matrix was obtained from the ECMWF for temperature and humidity profiles and for surface temperature. From the background error covariances, the background state vector is generated by random perturbation of the true state vector.

The background decomposition parametersgclr and the associated error g are derived from the outputs of the AVHRR radiance analysis if available. Otherwhise, gclris set equal to zero and its associated error g is set to 108.

Ozone profile, as well as other atmospheric and surface parameters other than temperature and humidity, are supposed to be perfectly known. For these parameters, the true state values given in the synthetic database are thus considered for the direct radiative transfer simulations.

The inputs of the 1DVAR inversion are the first guess estimate of the whole state vector, composed with the background atmospheric and surface states and the first guess estimate of the jk, and the IASI level 1c data, as well as their respective error figures. An incremented estimate of the state vector is obtained from 1DVAR inversion, which considers the Jacobians of the RTM regarding atmospheric and surface parameters, and the Jacobians of equation (4) regarding the decomposition parameters. This estimate of the state vector theoretically allows to compute a new set of Rguessij and Ri,j. Then the 1DVAR is iterated. The Newtonian iteration scheme proposed by Rodgers (1976) has been set up. A convergence check is used as follow: We compute at each iteration i the quantity  given by

(9)

where N is the number of vertical discretisation levels and xij represents the parameter vector or analysis state (temperature/humidity) at the iteration i and level j. Ajj is the element jj of the error covariance matrix for the analysis state, which is mainly constant regarding the iterations. The maximum value of the N level values j is retained in order to check the convergence is achieved for all levels. The iterative process is terminated when  is lower than 0.01. A sounding is rejected if the solution has not converged within 20 iterations.

The decomposition parameters are computed by inversion of Equation (2). If the AVHRR radiances analysis is not available, the algorithm is limited to the cloud decontamination functionality, and atmospheric and surface parameter retrievals and iterations of the algorithm are required.

The RTIASI Radiative Transfer Model (Matricardi and Saunders, 1999), is used to calculate the IASI radiances for a non-scattering atmosphere. RTIASI is based on the parameterisation of the atmospheric optical depths in order to compute quickly the transmittances of the different gases. The algorithm uses a fixed gas transmittance coefficient file to compute the transmittances of the CO2, N2O, CO, CH4, O2, N2, F11 and F12 gases. This parameterisation allows to reduce considerably the processing time. It needs as input 43 layer-mean values of temperature, water vapour and ozone mixing ratio as well as surface variables like the emissivity, surface pressure, surface air temperature, surface type index and the skin temperature. It computes level 1c IASI clear sky spectra and the associated Jacobians.

A quality control step is implemented, mainly to check if the retrieved temperature at any level or the surface skin temperature exceeds 400 K or falls below 150 K. In this case the souding Is rejected. We also reject the sounding if the specific humidity exceeds 0.06 g/g or becomes negative.