Hyper fast radiative transfer for the physical retrieval of surface parameters from SEVIRI observations

This paper describes the theoretical aspects of a fast scheme for the physical retrieval of surface temperature and emissivity from SEVIRI data, their implementation and some sample results obtained. The scheme is based on a Kalman Filter approach, which effectively exploits the temporal continuity in the observations of the geostationary Meteosat Second Generation (MSG) platform, on which SEVIRI (Spinning Enhanced Visible and InfraRed Imager) operates. Such scheme embodies in its core a physical retrieval algorithm, which employs an hyper fast radiative transfer code highly customized for this retrieval task. Radiative transfer and its customizations are described in detail. Fastness, accuracy and stability of the code are fully documented for a variety of surface features, showing a peculiar application to the massive Greek forest fires in August 2007.


Introduction
Geostationary satellite platforms have the capability of monitoring large portions of our planet with the advantage of having temporal continuity between consecutive observations. Such continuity could enable to solve the diurnal cycle of many of geophysical parameters, and should be properly exploited when data are analyzed, in order to maximize the information content extracted by data themselves. In particular, the observational continuity can be useful to correctly constrain the physical retrieval of the geophysical parameters of interest, even with different time scales of variability. This job has been already done with success [1] on geostationary sensors data by adopting a Kalman Filter (KF) strategy [2].
Even though the observational continuity is essential to follow the evolution of the geophysical parameters of interest, the continuous data collection is very demanding from the computational point of view: an enormous amount of data must be processed in a very reduced time. This means that, besides a solid mathematical framework for inversion, one should also have a very fast and accurate way to perform radiative transfer calculations, which are the bases for physical inversion. In this paper, we will address this peculiar problem by describing the implementation of an hyper fast radiative transfer model for simulating the radiance as it is observed by the SEVIRI instrument, which is on board of the Meteosat Second Generation (MSG) platform, currently the operative core of the geostationary meteorological programme by EUMETSAT.
The aim of the KF methodology implementation, largely addressed in [1], is to develop a scheme in which the retrieval of the parameters of interest was driven exclusively by data, not combining a Numerical Weather Prediction scheme (NWP, see e.g. [3]) with this machinery. This points out why having a fast radiative transfer code is really a fundamental piece of the work: this approach relies only on data and the efficiency of Optimal Estimation in fitting observed radiance with radiative transfer calculations. Since the task we have dealt with is limited to surface temperature and emissivity retrieval, we have optimized the radiative transfer code paying particular attention to the accuracy of calculations in those spectral regions where atmospheric absorption is minimal. Section 2 will describe the model in its details; in section 3 we will show the results obtained with the whole procedure on the case study of the Greek fires occurred during August 2007. Finally, conclusions will be drawn in Section 4.

The SEVIRI imager
The SEVIRI imager, on board of Meteosat-9 geostationary platform, produces a complete scan (from now on Full Disk, FD) of the visible hemisphere once every 15 minutes. The FD image is a cube with a spatial resolution of 3 km, and 12 broad spectral channels, 8 of which are in the thermal infrared band, that of our interest. These eight SEVIRI channels cover the spectral range [3.9 µm; 13.4 µm]. The spectral response for each channel is plotted in Fig. 1 Figure 1. Spectral response of the eight thermal infrared SEVIRI channels. They are superposed to a typical IASI [4] spectrum at a resolution of 0.5 cm −1 to highlight the main spectral features that yield the radiance in each SEVIRI channel. SEVIRI spectral response functions are in arbitrary units.
It is clearly visible that four of the eight SEVIRI infrared channels (numbers 1, 4, 6 and 7 in Table 1) are centred on atmospheric windows, hence the corresponding radiance values are only marginally affected by gaseous absorption. It should also be noted that Channel 1, though it is entirely located in an atmospheric window, is also characterized by the lowest radiance value; what is more, the radiance observed in this channel is not easy to be modelled given the solar contamination during the day and other non-LTE effects (i.e. multiple scattering).
The way in which the radiative transfer code has been tailored is mainly driven by the fact that the problem of retrieving surface properties does not require to work with a very complicated (and consequently, slow) atmospheric model to obtain accurate values of radiance in the three SEVIRI channels we use for retrieval. In addition, SEVIRI spectral response function is far broader than any gas absorption line: this enables us to perform radiative transfer calculations at a relatively coarse resolution, which in turn steeply reduces computational burden and time.

Radiative transfer basic elements
The forward radiative transfer model we describe here is called σ-SEVIRI, and has been specifically developed for this instrument. The model's name is inspired by the σ-IASI [5] from which this model (and other similar, see e.g. [6]) has been derived. The radiative transfer scheme we adopt has been already extensively validated in many studies [7,8,9,10] as far as both radiative transfer itself and retrieval applications are concerned. The σ-SEVIRI model solves the complete radiative transfer equation (RTE), which expresses the total radiance observed by the sensor, and has the following form: where σ is the wave number, ǫ g surface emissivity, T s surface temperature, B the Planck function, τ the atmospheric transmittance, z the height on the surface, I s the solar irradiance. The first term of the RTE is the radiance directly emitted from the surface, scaled with the atmospheric total transmittance τ 0 ; the second term represents the radiance emitted from the atmosphere integrated over the whole observation path; the third term is the down welling radiance, that is that emitted from the atmosphere, reflected from surface and then reaching the sensor. Both Lambertian and specular reflections can be modelled, in order to use the code both on land and sea surface. Finally, the last term is the solar irradiance which is reflected from surface; this last term becomes important only at wave numbers greater than ≈2000 cm −1 .
The RTE in the form (1) cannot be directly implemented, but it requires to be a bit manipulated. The key passage consists in its discretization, which is done by dividing the atmosphere in a certain number of layers, namely N L . We choose to use a fixed pressure grid, where each layer is rigorously defined by its pressure boundaries; each layer is assumed to have an uniform temperature and composition. Equation (1) then becomes: with the same notation as above. Equation (2) reduces the computed radiance to sums/products of transmittances τ j from the bottom of the j-th layer to the atmospheric top. To accomplish calculations of transmittances, the σ-SEVIRI model uses a wide look-up table for the optical depth. This table is pre-calculated using the reference Line-By-Line Radiative Transfer Model (LBLRTM) most recent available version [11].
Before giving any detail about the two steps by which we have speeded up the σ-SEVIRI code, it is worth to explain how transmittances are expressed and how the look-up table is built. For simplicity, let us consider a single layer j in our pressure grid. At a given wave number σ, the partial transmittance of the layer is expressed through its relation with the optical depths of the S absorbing species in that layer: where χ i,j (σ) is the optical depth of the i-th absorber in the j-th layer. In σ-SEVIRI, we take into account only gases as absorbing species, and the strategy we adopt is to exploit the dependence of the gases optical depths both on their concentration and the layer temperature T j . In detail, the optical depth χ i,j (σ) is expressed as a simple second-order polynomial of the temperature for carbon dioxide, ozone and all the other gases but water vapour treated all together as mixed gases (i=2, 3,4). For water vapour (i=1), instead, we add a further coefficient in order to take into account effects depending on its abundance, such as self-broadening [12]. Hence where q i,j is the molecule concentration, and ∆T is the difference between the actual layer temperature and a reference temperature, that is the central point of the parameterization. The c ·,i,j coefficients are stored in the look-up table used by the code. In principle, calculations should be performed at infinite resolution in such a way to take into account the single spectral line contribution. For Earth atmosphere, this would correspond to ∆σ=10 −4 cm −1 , making calculations extremely slow.
In the case of SEVIRI, thanks to its broad channels, we can work with a far coarser resolution, without losing in accuracy. The equivalent optical depths for water vapour and the other gases, can be written as follows for all the other gases. The symbol σ indicates the average over wave number. The averaged coefficients are obtained via averaging transmittances and not optical depths, because of the form in which RTE is solved by the model, which deals with layers' transmittances. In the work described in [1], optical depths coefficients' database was generated at a resolution of 10 −1 cm −1 , and with an atmosphere model made by N L =60 layers. In this work, instead, the first step to speed up the code has been to enlarge the binning step until 5×10 −1 cm −1 . Radiances calculated at this resolution are then convolved to SEVIRI resolution through its Instrument Spectral Response Function (ISRF). On the one hand, this has yielded a drastic reduction of the look-up table size; on the other hand, this has cut by a factor of 5 the time needed to compute a single spectrum and all its derivatives, that are used in the retrieval framework. The second optimization step is related to the fact that we use only window channels, where atmospheric absorption is only marginal. Thus, there should be no or negligible bias in computed radiances if we reduce the number N L of atmospheric layers used to perform radiative transfer calculations. Hence, we have chosen to further optimize the code choosing N L =25. Such a choice is for convenience, since in the KF machinery we use the atmospheric profiles provided by ECMWF (European Centre for Medium-range Weather Forecast) whose analyses data are given on a 25-layer pressure grid. As a consequence, no preliminary interpolation is required to fit data into the pressure grid, further speeding calculations. The comparison between the two pressure grids (60 and 25 layers) is given in Fig. 2, while an inter-comparison between sample calculations performed with 60 and 25 layers is provided in Table 1. The tabulated radiances evidence how the coarsening of the resolution at which calculations are performed has introduced practically no bias in the results, while the reduction of N L has an effect which is at most of the same order of magnitude of the SEVIRI radiometric noise, again introducing no significant bias in the radiance calculated in the three atmospheric window channels we use in KF retrieval.
With these modifications to the σ-SEVIRI model, the inversion scheme to retrieve surface parameters has been speeded up by a factor of two with respect to the original configuration.
Because of these results, this boosted version of σ-SEVIRI suits well in an application, like that we are going to show here, in which a very huge amount of data has to be analyzed in near-real time. σ-SEVIRI has been integrated in a physical retrieval framework, whose core is an iterative, non linear Optimal Estimation Gauss-Newton scheme, that has been already  Table 1. Sample radiance calculation results for the three window channels used in the KF scheme. The comparison points out that the differences between calculations made at 10 −1 cm −1 (first row), those made at 5×10 −1 cm −1 (second row), both with a 60-layer pressure grid, and those made at a resolution of 5×10 −1 cm −1 and a 25-layer pressure grid (third row) are below SEVIRI radiometric noise (fourth row). Radiances are expressed in W m −2 sr −1 (cm −1 ) −1 and are calculated for the U.S. standard climatology profile AFGL6 [13]. Sampling-layers Radiance @12.0 µm @10.8 µm @8.7 µm 10 −1 cm −1 -60 layers implemented in many contexts. In order not to be too long-winded, here we will not recall any element neither of the inversion procedure nor of the Kalman Filter application to this case. The interested reader can refer to [1,8,14,15] for details about these aspects.

Results on a case study: the Greek fires in August 2007
The case study we are going to show here proves the double potentiality of both the radiative transfer model σ-SEVIRI to calculate rapidly and accurately the radiance, and the capability of the Kalman Filter to follow the time cycle of skin temperature and emissivity independently of data voids and/or clouds, and without wasting computational resources. We examine SEVIRI data that have been acquired on the Greek peninsula along a month, from the 23rd of August to the 21st of September, 2007. Such a choice is due to the fact that, in the last decade of August 2007, many very extended fires affected the Peloponnese region. Such events are usually associated with sharp peaks in the observed radiance of the corresponding pixel, and to some kinds of discontinuities in the emissivity. Such phenomena, together with extended smoke clouds, are not only interesting from a scientific point of view, but they also offer the ideal terrain to test the robustness and accuracy of the KF procedure. The target area we have investigated covers the whole Greek peninsula. We have selected 8619 MSG-9 pixels, 4047 over land and 4572 over sea for the period mentioned before. SEVIRI observations have undergone through a first selection, aimed to exclude from the KF analysis all the data corresponding to cloudy scenes, equalizing them to data voids. Then the KF has been run on each single SEVIRI pixel, without considering the spatial correlation between them. As said, we have adopted a simple persistence model to make the system strongly data-driven.
The parameters we want to retrieve are the three emissivity values in the SEVIRI channels 4, 6 and 7, and the skin temperature T s . Hence, the parameters column vector v at a generic time t is v=(ǫ 7 , ǫ 6 , ǫ 4 , T s ). To begin with, Fig. 3 shows the averaged maps of the retrieved skin temperature (land and sea) and land emissivities in the three SEVIRI channels @8.7, @10.8 and @12.0 µm. Such results are useful to point out the capability of the KF approach, together with  a fast radiative transfer model, to process effectively and accurately a large volume of data. It is evident how the KF captures the most salient surface features as far as emissivity is concerned, while the orography is well visible from skin temperature map. In this case, one can evidence sea currents and/or the presence of deep water fluxes towards sea surface.
As already stated, this peculiar dataset was selected because of the occurrence of large fires. In order to verify the KF scheme behaviour in presence of such events, we have selected the SEVIRI pixels interested by fires, and plotted the corresponding retrieved skin temperature and emissivities time series. Results are shown in Fig. 4. Fire events are mostly associated to sharp, short peaks in both retrieved skin temperature and observed radiance, more visible in the SEVIRI channel @3.9 µm, not used for retrieval purposes. Such peaks in the skin temperature usually associate with to similar sharp variations in the retrieved emissivity. During the fire event, indeed, the emissivity @8.7 µm goes down, reaching a value of about 0.9, typical of a grey body [16]. It is interesting that, notwithstanding such anomalous events, the KF is stable with respect to them, as the temporal behaviour of the retrieved quantities shows. Considering that, for a single SEVIRI pixel, the KF scheme requires no more than one minute for a one-month data series for a single pixel on a single modern CPU, one can perform a realtime data analysis for a region like that considered in this case, even using a reduced amount of computational resources.

Conclusions
We have presented a new version of the already existing σ-IASI model which is customized for the SEVIRI radiometer. The model, which is able to calculate the Earth radiance as it is observed by SEVIRI in a few hundredths of second, has been speeded-up preserving accuracy in those SEVIRI channels where the atmospheric absorption is negligible, making them useful for the retrieval of skin temperature and spectral emissivity.
The code has been used in a Kalman Filter framework to retrieve the temporal evolution of skin temperature and emissivity. The forward radiative transfer model, joint with the KF methodology, has proved to be able to retrieve such parameters despite possible data voids responsible for local temporal discontinuities in the retrieval. The whole procedure works efficiently both on restricted and on large spatial regions, where each pixel is processed individually. The machinery has been applied on a challenging case study, which involves SEVIRI observations acquired on the Greek peninsula during the huge fires occurred in August 2007. The KF has successfully processed such data, highlighting emissivity anomalies in the spatial pixels involved in fire events.
For its completeness and speed, the KF methodology is suitable to be applied on large regions, even on the SEVIRI Full Disk, in order to build a comprehensive skin temperature and emissivity atlas [17], with the aim to support other retrieval procedures, providing some background for surface, and to interface with Numerical Weather Prediction (NWP) models.