A New Data Fusion Neural Network Scheme for Rainfall Retrieval Using Passive Microwave and Visible/Infrared Satellite Data

Abstract

A new data fusion technique based on Artificial Neural Networks (ANN) for the design of a rainfall retrieval algorithm is presented. The use of both VIS/IR (VISible and InfraRed) data from GEO (Geostationary Earth Orbit) satellite and of passive microwave data from LEO (Low Earth Orbit) satellite can take advantage of both types of sensors reducing their limitations. The technique can reconstruct the surface rain field with the MSG-SEVIRI (Meteosat Second Generation–Spinning Enhanced Visible Infrared Imager) spatial and temporal resolution, which means 3 km at the sub satellite point and 5 km at mid-latitudes, every 15 min, respectively. Rainfall estimations are also compared with H-SAF (Hydrology Satellite Application Facility) PR-OBS3A operational product showing better performance both on the identification of rainy areas and on the retrieval of the amount of precipitation. In particular, in the considered test cases, results report an improvement in average of 83% in terms of probability of rainy areas detection, of 45% in terms of false alarm rate, and of 47% in terms of root mean square error in the retrieval of the amount of precipitation.

1. Introduction

The prevailing element in the Earth water cycle is the precipitation, both liquid and solid, which therefore assumes a crucial role in hydrological budget. A better knowledge of the spatial and temporal distribution of this variable has positive consequences in many different fields, like climatology, weather forecasting, agriculture, fishing, and flooding early warning systems. Despite its importance, precipitation is one of the most difficult geophysical parameter to be quantified, particularly in those regions (e.g., ocean, remote lands, and developing countries) where there are sparse or even lacking ground measurements. A significant improvement to monitor precipitation for those regions come from meteorological satellites because of their spatial coverage and temporal sampling capabilities over large areas. Areas covered by rain gauges network or weather radars can also benefit from satellites precipitation estimation to compensate the usual irregular spatial distribution of rain gauges or the orography problems that could affect weather radars.

Many algorithms for satellite-based rainfall retrieval exist. An important category exploits the attenuation, measurable by microwave radiometers on satellites, due to the presence of precipitation particles on Earth [1,2,3,4]. However, the passive microwave (PMW)-based techniques observations are currently only available from LEO satellites, typically resulting in two observations per day per satellite, which represent a significant operational constraint. In addition, the spatial resolution also may be not adequate, being of the order of 50 × 50 km for the low frequency channels that are used over the ocean and no better than 10 × 10 km for the higher-frequency channels used over the land.

Algorithms based on data acquired from GEO satellites (VIS/IR sensors) represent another approach [5,6]. In this case both spatial and temporal resolution is increased but the physical connection, based on scattering and absorption mechanisms, between the visible and infrared part of the spectrum and precipitation is not as strong as for microwave, affecting the final accuracies. A rising number of techniques are being also developed to exploit the synergy between the polar-orbiting PMW retrievals (irregular, and more direct) with the geostationary VIS/NIR/IR observations (regular, and less direct). Most of them rely on some empirical inter-calibration operations [7]. This type of approach, however, may have generalization limitations as the calibration holds only in a specific geographic area or for a specific class of events. A motion-based technique has been also proposed [8], where the IR data provide a reasonable measure of cloud movement, which is then used to advect or morph the more direct PMW data between the successive satellite overpasses. A more conceptual cloud-development model, driven by GEO satellite imagery, and locally updated using MW-based rainfall measurements from LEO platforms, is described in Bellerby et al. [9], such an approach has been refined in [10] where cloud classification is a specific step incorporated in the precipitation estimation method.

In the context of the data fusion methods, the use of the Meteosat Second Generation (MSG) satellite mission from the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) has been only sparsely addressed [11]. MSG is a very important mission in the context of GEO satellites. It carries the SEVIRI instrument, a line-by-line scanning radiometer, which provides image data in four Visible (VIS) and Near-InfraRed (VNIR) channels and eight InfraRed (IR) channels. Although “operational” precipitation products are routinely produced within concerted programs, such as the EUMETSAT Satellite Applications Facility on Support to Operational Hydrology and Water Management (H-SAF) [12], more investigation is needed to analyze the potential of a data fusion approach merging information provided by PMW and MSG. The objective of this paper is to provide a contribution along this direction. The final goal is to improve the accuracy of up-to-date precipitation products based on EUMETSAT SEVIRI data. To this purpose, we propose an innovative methodology for merging MW and VIS/NIR/IR MSG observations based on neural networks (NN).

In fact, rainfall estimation is a problem that contains many unknown processes. The inverse problem of the radiative transfer process involves an ill-posed scenario which yields many solutions for a given set of brightness temperatures. In such circumstances, NN models, which can learn input–output significant relationships directly from the training data and, if the training process is carried out properly, are able to generalize out of them, can be very useful. The relationship between the amount of radiation reaching the satellite sensor and the intensity (and the type) of surface rain is nonlinear and non-Gaussian, then it is well-matched to NNs estimation capabilities. Additionally, NNs can handle huge amount of data in an automatic and quick way without the need of new training sessions. On the other hand, physically based models often require a great quantity of computational resources, then NNs might take a considerable advantage in an operational context. Lastly, noisy data are well managed by NNs, fault tolerance is indeed another advantage of using such an approach [13].

On the other hand, the use of NN algorithms has been often considered suitable for parameter retrieval in atmospheric remote sensing [14,15] and neural techniques for rainfall monitoring involving some data fusion approach have been applied in the past [16,17]. More in detail, if we restrict the research to the use of MSG data, a first attempt of using NN and a fusion technique for rainfall retrieval can be found in [18]. In this case only the IR10.8 SEVIRI channel was considered together with SSM/I measurements. More recently, in [19] a NN approach using MSG and ground based radar measurements is presented to estimate rainfall over north of Algeria.

In this study, a novel NN architecture using as inputs MSG SEVIRI multispectral VIS and IR data and capable to generate rainfall maps at the MSG SEVIRI time resolution is presented. For the training phase, the possibility to use LEO PMW data from all in-orbit satellites, “fusing” data from six different sensors has been explored. The methodology is based on the design of two neural architectures, the first one to extract rainy pixels, the second one to provide the corresponding estimated rainfall value. The resolution issues are handled by means of a rigorous and accurate collocation procedure and by considering mean and standard deviation values over 3 × 3 and 5 × 5 pixel neighborhood. The algorithm performance is tested through comparison with a more conventional one, such as PR-OBS-3A, the EUMETSAT Hydrological Satellite Application Facility (HSAF) precipitation product. It has to be noted that both daylight and night time scenarios have been addressed. These latter characteristics together with the concurrent richness of the data used, either in terms of MSG bands considered as NN inputs or in terms of the microwave information exploited in the training phase, represent an additional significant element of innovation in this field.

2. Datasets

2.1. Passive Microwave Data-GPM Mission

The GPM Mission, is an international space network of satellites designed to provide precipitation observations every three hours anywhere around the world. As shown in

Table 1. Launch schedules and life spans of Global Precipitation Measurement constellation satellites. Satellites used in this study are in bold.

SATELLITE	SENSOR	LIFE SPAN
Global Precipitation Measurement Core	GPM Microwave Imager/Dual Frequency Precipitation Radar	2014–beyond 2021
Global Change Observation Mission-Water 1	Advance Microwave Scanning Radiometer 2	2012–2018
Megha-Tropiques	Sondeur Atmospherique du Profil d’Humidite Intertropical par Radiometrie	2011–2017
Defence Meteorological Satellite Program-Flight 17	Special Sensor Microwave Imager/Sounder	2011–2015
Defence Meteorological Satellite Program-Flight 18	Special Sensor Microwave Imager/Sounder	2011–2016
Defence Meteorological Satellite Program-Flight 19	Special Sensor Microwave Imager/Sounder	2014–beyond 2021
Defence Meteorological Satellite Program-Flight 20	Special Sensor Microwave Imager/Sounder	2020–N/A
Meteorological Operational Satellites-A	Microwave Humidity Sounder	2011–2014
Meteorological Operational Satellites-B	Microwave Humidity Sounder	2012–2017
Meteorological Operational Satellites-C	Microwave Humidity Sounder	2016–beyond 2021
National Oceanic and Atmospheric Administration-18	Microwave Humidity Sounder	2011–2017
National Oceanic and Atmospheric Administration-19	Microwave Humidity Sounder	2011–2018
National Polar-orbiting Partnership	Advanced Technology Microwave Sounder	2012–2019
Joint Polar Satellite System-1	Advanced Technology Microwave Sounder	2017–beyond 2021

, the GPM constellation consists of several operational, research, and development satellites that include (or will include) the operational US NOAA and European MetOp polar orbiting series, the US Defence Meteorological Program satellites, the US Joint Polar Satellite System, together with the Japanese GCOM-W1 (Global Change Observation Mission-Water) mission, and the French-Indian Megha-Tropiques mission. The GPM core satellite, the “observatory”, carries a dual frequency precipitation radar (DPR) using Ku band at 13 GHz and Ka band at 35 GHz and a canonical-scanning microwave imager (GMI) operating from 10 to 183 GHz with 13 channels. It builds upon the Tropical Rainfall Measuring Mission (TRMM) heritage [7] but provides retrieval accuracy and higher sensitivity to slight rainfall (smaller than 2.5 mm/h according to the definition of the Glossary of Meteorology [20], such a reference will be considered also for other rainfall categories mentioned in the paper) and snowfall detection relative to TRMM. DPR and GMI, during their life-time, will generate an observational database on the microphysical properties of precipitating particles. The GPM concept is based on the deployment of a “Core” observatory carrying active and passive microwave sensors in a non-Sun-synchronous orbit to serve as a physics observatory to gain insights into precipitation systems and as a calibration reference to unify and refine precipitation estimates from a constellation of diverse satellites. Through the inter-calibration of radiometric data obtained from this constellation of different sensors, an accurate unified precipitation retrieval algorithm, using a common hydrometeor database and working in the different environmental conditions, has been developed [21]. These data set will be used to provide rainy pixels and to train the NN performing the rainfall retrieval.

2.2. Visible/Infrared MSG SEVIRI Data

The Meteosat Second Generation satellites have an altitude of about 36,000 km above the equator at 0° W and 41.5° E. With respect to its precursor instrument MVIRI (Meteosat Visible and InfraRed Imager,), the EUMETSAT Spinning Enhanced Visible and InfraRed Imager (SEVIRI) [22] provides significantly increased information due to an imaging-repeat cycle of 15 min (30 min for MVIRI) and 12 spectral channels (3 channels for MVIRI), quantization with 10 bits per pixel (8 for MVIRI), and image sampling distances of 3 km at nadir for all channels except the high-resolution visible with 1 km (for MVIRI, 5 and 2.5 km, respectively). These enhanced time-space resolution performances allow to retrieve timely information on rapid weather development and are particularly helpful to rainfall retrieval algorithms. The high resolution visible channel has an IFOV of 1.67 km, and the oversampling factor is 1.67 that corresponds to a sampling distance of 1 km at nadir. The corresponding values for the eight thermal IR and the other three solar channels are 4.8 km IFOV, with an oversampling factor of 1.6 and a sampling distance of 3 km for nadir view. Eight of the channels are in the thermal infrared, providing, among other information, permanent data about the temperatures of clouds, land, and sea surfaces. One of the channels is called the High Resolution Visible (HRV) channel and has a sampling distance at nadir of 1 km, as opposed to the 3 km resolution of the other visible channels. Using radiation with wavelengths physically and significantly interacting with ozone, water vapor and carbon dioxide molecules, MSG satellites allow meteorologists to analyze the characteristics of atmospheric air masses in the three dimensions. The improved horizontal image resolution for the visible light spectral channel (1 km as opposed to 2.5 km) also helps weather forecasters in detecting and predicting the onset or end of severe weather. The crucial properties to detect rain (optical thickness, droplet size, cloud phase, cloud top temperature, and water vapor content) can be retrieved using the SEVIRI channels [19,23]. This aspect will be more deeply addressed in Section 3.1 dealing with the selection of the radiances measured by SEVIRI to be used as input for our NN algorithms.

2.3. MSG SEVIRI Cloud Mask

The MSG SEVIRI cloud mask (CMa) product permits the identification of cloud free pixels in a satellite scene with a high confidence [24]. Indeed, cloudy pixels can be often identified, because they appear colder (at 10.5 μm) and/or brighter (at 0.6 or 0.8 μm) than cloud free areas. A fine analysis of their respective spectral behavior is nevertheless needed to perform a full cloud detection. For example, low clouds identification at night-time relies on their low emissivity at 3.8 μm, whereas thin cirrus clouds can be identified due to their different emissivity at 10.5 μm and 12 μm. Cloud free areas covered by snow or ice are identified at daytime with their very low reflectivity at 3.7 μm and high reflectivity at 0.6 μm, whereas oceanic cloud free areas affected by sun glint are identified with their very high reflectivity at 3.8 μm. The CMa identifies pixels that are contaminated by either clouds, dust, or snow/sea ice. The problem to be solved is to automatically predict, with sufficient accuracy, brightness temperatures and reflectance of cloud free areas, so that any discrepancy between the measured and predicted values can be used to detect contaminated pixels. Therefore, the algorithm is based on multispectral threshold technique applied to each pixel of the image. A first process allows the identification of pixels contaminated by clouds or snow/ice. This is complemented by an analysis of the temporal variation (between individual image repeat cycles) of some spectral combination of channels (to detect rapidly moving clouds) and a specific treatment combining temporal coherency analysis and region growing technique (to improve the detection of low clouds). This first process is carried out to determine the cloud cover category of each pixel (cloud-free, cloud contaminated, cloud filled, snow/ice contaminated, or undefined/non-processed) and compute a quality flag on the processing itself. A second process, allowing the identification of dust clouds and volcanic ash clouds, is applied to all pixels (even already classified as cloud-free or contaminated by clouds). The result is stored in the dust cloud and volcanic ash cloud flags. Within the encoding of the product the pixels are identified as one of the following:

• Clear sky over water;
• Clear sky over land;
• Cloud;
• No data (if none of the classifications above is possible).

This product, as provided by EUMETSAT [25], will be used to select, at the SEVIRI space resolution, for both the NN detecting rainy pixels and the NN providing rainfall estimation, only cloudy pixels.

2.4. Validation Weather Radar Data

To validate the methodology proposed in this paper, weather radar over Europe are used. They are provided by the Operational Programme for the Exchange of Weather Radar Information (OPERA) in form of a composite mosaic of 200 weather radar across the Europe [26] (see Figure 1). The number of dual-polarization radars in Europe is growing steadily and the current count is 48. Most of the radars operate at C-band (168 sites), but many S-band radars are installed in the south of Europe (33 sites). The most northern radar is in Hasvik (Norway) at 70.6° N, and the most southern one is in Gran Canaria (Spanish island off the African coast, southwest of Europe; not shown on the map) at 28.0° N. The most western radar is located in Keflavik (Iceland) at 22.6° W, and the most eastern one in Berlevåg (Norway) at 29.0° E. The median distance between two neighboring radars in the network is 128 km, and the smallest distance between a pair is just 7 km (a C-band and an S-band radar in Oradea, Romania) and the largest distance is only 276 km. These mosaic products are generated on a Cartesian grid covering the whole of Europe (area: 3800 × 4400 km) with a horizontal grid spacing of 2.0 km. Every mosaic pixel is a weighted average of the lowest valid pixels of the contributing radars, weighted by the inverse of the beam altitude. Polar volume data cells within a search radius of 2.5 km of the composite pixel are considered. Each pixel of the mosaic is provided with an error-based quality index (range 0 to 1). This index is calculated by combining information on radar properties, measurement position, and meteorological conditions in a form that relates directly to rain rate error. The performance of the quality index have been demonstrated to be positive, smoothing rain rates at the boundaries between radar fields, and reinstating some rain that was missed in the operational composite. As reported in [27], the statistics show no significant reduction in compositing skill, and in fact suggest improved skill at high rain rates. These data indicate that use of a quality index for compositing increases the accuracy of surface rain rate estimates at high intensity, without significantly decreasing accuracy at lower rain rates.

2.5. PR-OBS-3A Precipitation Product

PR-OBS-3A, a EUMETSAT Hydrological Satellite Application Facility (HSAF) precipitation product, is the benchmark to evaluate the data fusion technique here proposed. In fact, this product and the one yielded by our methodology have the same spatial and temporal resolutions which is a very important aspect for a reliable comparison. A complete description of the algorithm is provided in [28]. PR-OBS-3A provides instantaneous precipitation rain rate retrieved from blending information produced by PMW and IR observations using SSMIS, MHS and SEVIRI sensors (channel 9, wavelength 10.8 µm) [29]. The product is designed to provide a better time continuity and even spatial coverage with respect to the relatively infrequent precipitation observation of LEO instruments, processing, and calibrating IR observation from GEO platforms. PR-OBS-3A is calibrated directly by means of a “rapid update technique” lookup table of precipitation based on time-space coincident GEO and LEO information.

PR-OBS-03 performs relatively well for the highest precipitation class when compared either with the radar or with rain gauge observations. Its performance starts to decrease in the intermediate precipitation class, until to become very poor for the lower precipitation regime. The analysis of monthly statistical scores shows a seasonal trend: it performs better during winter time, but an underestimation of precipitation rate is always observed [30,31].

2.6. Case Studies

To reflect the seasonal meteorological variability and the different type of precipitations, four case studies, summarized in

Table 2. Case Studies.

Date	Time Range	Type of Event	Country and Boundaty
5 June 2015	14:00–23:00	Thunderstorm	Belgium (48 to 52)° N-(2 to 8)° E
20 February 2016	00:00–22:00	Rain associated to winter storm	Belgium (48 to 52)° N-(2 to 8)° E
12 October 2015	03:00–12:00	Thunderstorm	Poland (49 to 55)° N-(14 to 25)° E
8 July 2015	13:00–17:00	Cyclone over North Sea	Hungary (45 to 50)° N-(16 to 22)° E

, have been selected to create an appropriate dataset for this study.

5 June 2015-Belgium

The weather in Belgium was influenced in the beginning by a high-pressure zone in central Europe. This system was redirecting tropical continental air masses to western Europe while a low-pressure zone from the Atlantic, accompanied by cold fronts, was slowly approaching the British Isles. Such conditions were favorable for the development of a thermal trough over Belgium, resulting in thunderstorms and intense rainfall in the afternoon and evening. The zones of highest rainfall accumulation were the central parts of the country and a narrow zone near the coast.

20 February 2016-Belgium

This scenario is characterized by a vast zone of very low pressure in the south of Iceland, moving towards Scandinavia. Several fronts and strong winds are accompanying this system. Belgium is under the influence of warm and then cold fronts. During this event, warm air from south-west caused a significant increase in temperatures while precipitation lasted for many hours over several places in Belgium.

12 October 2015-Poland

On 12 October 2015, a significant weather change took place over Poland. The sunny conditions were replaced by cloudy weather characterized by both solid and liquid precipitation in the southern regions of the country. The change was induced by low weather system directing the advection and thus the precipitating cloud movement. The occluded front located on the Polish/Czech border was the main driver for the rainfall events over Poland on that day. The massive precipitation system was travelling from SE to NW with slight northerly deviation. The stable stratiform precipitation was recorded over southern regions of the country.

8 July 2015-Hungary

On 8 July 2015, an extra-tropical cyclone developed over North Sea. The cold front of this cyclone caused lot of convective clouds in Europe. In Hungary (around Budapest) lot of thunderstorms brought hail, storm wind and high precipitation intensity during the afternoon.

3. Methodology

The proposed technology is summarized in Figure 2 where it can be seen that, for the specific event, the core processing is based on two NNs, named NN1 and NN2. Part (a) and part (b) show how NN1, providing the rain mask, and NN2, performing the precipitation rate estimation, respectively, are trained. Part 9(c) reports the whole operational processing chain. Once the two NNs are created, they can be used to reconstruct the precipitation field associated to the event with a time sample of 15 min and with SEVIRI spatial resolution, filling the temporal and geographical gaps between the LEO satellite passages. The following paragraphs will address the main technical issues involved in the design of the double-NN retrieval scheme.

3.1. Selection of the NN Input Data

A crucial issue in the design of a NN algorithm is the selection of appropriate inputs. The information considered as input for the two NN algorithms presented in this study, one for day time and the other one for night time, is now described. According to [32,33] a first piece of information directly linked to the probability of a cloud of being rainy is represented by the Cloud Water Path (CWP). The CWP is the amount of water vertically integrated in the cloud, and given some assumptions about cloud vertical structures, it is related to the cloud optical thickness and to the droplet effective radius where the effective particle radius is defined by the ratio between the third and the second power of the droplet spectrum [34]. To characterize the geometrical thickness of a cloud, the integration of the extinction coefficient integrated over the cloud geometrical thickness, called Cloud Optical Thickness (τ), is also used. To retrieve CWP using SEVIRI data, reflectance of VIS0.6 and NIR1.6 can be used [6,32]. High values of VIS0.6 reflectance correspond to high optical depth of cloud and low values of NIR1.6 reflectance indicate large particles in the cloud. This means that a large CWP is obtained when high values of VIS0.6 coincide with low values of NIR1.6. During the night, the cloud optical thickness is well defined by brightness temperature difference BTIR10.8-BTIR12.1 [35], while to infer information about microphysical and optical cloud properties, the brightness temperature difference BTIR3.9-BTIR10.8 and BTIR3.9-BTWV7.3 may be used [36,37,38]. Another distinctive aspect of the raining clouds is the thermodynamic phase. It may be assumed that if both liquid (droplets) and ice water (crystals) coexist then the cloud is more likely to bring rain [32,33,35,39,40]. Two good indicators of the thermodynamic phase of the cloud are the cloud top temperature (CTT) (provided by BTIR10.8) and the brightness temperature difference BTIR8.7-BTIR10.8, used concurrently [33,41]. To take into account the different temporal signature of convective (and orographic) processes with the weather fronts, [18,42] suggest to consider not only the current value of BT of a pixel, but also the value corresponding at the previous five acquisitions. Furthermore, to minimize geolocation and synchrony problems between SEVIRI and GPM dataset and to provide useful texture information of the clouds, average and standard deviation values over boxes of 3 × 3 and 5 × 5 surrounding the pixel were added to the input vector. On the base of the previous considerations, for each day, the input vector of the two NNs, for day-time and for night-time, are listed in

Table 3. Day time input.

Input #	Data Type	Explanation
1–6	BTIR10.8 t0 … t5	Brightness temperature of the pixel for the current IR image (t0) and the previous five images. Cloud Top Temperature, Cloud Phase, Frontal System (time).
7	Refl. VIS0.6	Cloud Optical Thickness, Cloud Water Path.
8	Refl. NIR1.6	Effective Radius, Cloud Water Path.
9	BTIR8.7-BTIR10.8	Cloud Phase.
10	BTIR10.8-BTIR12.1	Cloud Optical Thickness, Cirrus test.
11–15	MEAN3(1, 7, 8, 9, 10)	Mean value over a 3 by 3 pixel neighborhood. Minimization of geolocation and synchrony problems. Texture information about the cloud cover.
16–20	STD3(1,7,8,9,10)	Standard deviation of pixel value over the 3 by 3 pixel neighborhood. Texture information about the cloud cover.
21–25	MEAN5(1, 7, 8, 9, 10)	Mean value over a 5 by 5 pixel neighborhood. Minimization of geolocation and synchrony problems. Texture information about the cloud cover.
26–30	STD5(1, 7, 8, 9, 10)	Standard deviation of pixel value over the 5 by 5 pixel neighborhood. Texture information about the cloud cover.

and

Table 4. Night time input.

Input #	Data Type	Explanation
1–6	BTIR10.8 t0 … t5	Brightness temperature of the pixel for the current IR image (t0) and the previous five images. Cloud Top Temperature, Cloud Phase, Frontal System (time).
7	BTIR3.9-BTIR10.8	Microphysical and optical cloud properties.
8	BTIR3.9-BTWV7.3	Water Vapor Content (CWP).
9	BTIR8.7-BTIR10.8	Cloud Phase.
10	BTIR10.8-BTIR12.1	Cloud Optical Thickness, Cirrus test.
11–15	MEAN3(1, 7, 8, 9, 10)	Mean value over a 3 by 3 pixel neighborhood. Minimization of geolocation and synchrony problems. Texture information about the cloud cover.
16–20	STD3(1, 7, 8, 9, 10)	Standard deviation of pixel value over the 3 by 3 pixel neighborhood. Texture information about the cloud cover.
21–25	MEAN5(1, 7, 8, 9, 10)	Mean value over a 5 by 5 pixel neighborhood. Minimization of geolocation and synchrony problems. Texture information about the cloud cover.
26–30	STD5(1, 7, 8, 9, 10)	Standard deviation of pixel value over the 5 by 5 pixel neighborhood. Texture information about the cloud cover.

, respectively.

3.2. Neural Network Topology

Standard MLP architectures have been already proven very effective in atmospheric remote sensing retrieval problems [43,44] hence we used these types of models also in our study with tangent sigmoid nodes [45]. The selection of the best topology was based on an empirical way. The number of input neurons, hidden layers and nodes within it, and the training algorithm were chosen to obtain the best scores in statistical tests (see below) using an independent ground truth (weather radar). The two NNs (day-time and night-time) are based on a Multi-Layer Perceptron (MLP) feed-forward model with 30 input neurons, 45 units in a unique hidden layer, and one output neuron. The Matlab development environment has been considered for the design of both NNs.

3.3. Data Collocation

As far as the collocation of GPM data (providing rain rate in mm/h unit) with SEVIRI data is concerned, this is achieved in 5 steps.

• A timestamp is associated to each SEVIRI pixel. It is calculated starting with the nominal date and time provided in the SEVIRI data file-name and adding the scansion time at the latitude of the pixel (about 12 min at 45° N).
• Geographical area and the temporal window are defined. The spatial and the temporal window should be large enough to obtain a statistically significant dataset but not too wide in order to avoid climatologically different zones, for example the Mediterranean area and tropics, being processed at the same time.
• Using the corresponding EUMETSAT cloud mask, all the cloud-free SEVIRI pixels are filtered out from the dataset.
• For each GPM pixel belonging to the chosen area and temporal window, all the SEVIRI pixels at a distance less than 2.5 km (half SEVIRI pixel resolution at mid-latitude) and with a temporal distance less than 7.5 min (half SEVIRI repeat cycle) are selected.
• Among the SEVIRI pixels from step 4, the nearest in time is selected and collocated with the GPM pixel.

3.4. Training Issues

Precipitation is a statistical rare event. Analyzing many polar orbits, it turns out that just a few percent (or even less) of the pixels are rainy. Furthermore, considering cloudy pixels, only 10% to 30% of them are rainy as it is shown in Table 9. Using indiscriminately all the pixels of the images, the dataset would be totally unbalanced towards zero-rain points, and a NN trained with such a dataset will not be able to appreciate differences in rain rate during a storm or to distinguish between a short violent thunderstorm from a long stratiform rain. To avoid this problem a neural network classifier has been first created. A subset of the total data, choosing randomly from the original dataset an equal number of rainy and not-rainy pixels, is used to train a neural network classifier. GPM data are used as output of this classifier with a threshold of 0.2 mm/h (above the threshold the pixel is considered rainy; below not rainy).

Using the collocation procedure described above, with geographical areas and temporal windows provided in Table 2, the subset to create the rainy mask is obtained. The subset is randomly divided in three parts, 70% for training set, 15% for the validation set, and 15% for the test set. The training set is used for computing the network weights and biases. The error on the validation set is monitored during the training process. The validation error normally decreases during the initial phase of training, as does the training set error. However, when the network begins to overfit the data, the error on the validation set typically begins to rise. The network weights and biases are saved at the minimum of the validation set error. The test set error is not used during training, but it is used for the final evaluations. The numbers of pixels used to create the NNs are reported in

Table 5. Rainy mask subset-day-time.

Date	Total Pixels	Training Pixels	Validation Pixels	Testing Pixels
5 June 2015	61,418	42,992	9213	9213
20 February 2016	52,769	36,939	7915	7915
12 October 2015	89,709	62,797	13,456	13,456
8 July 2015	26,844	18,792	4026	4026

and

Table 6. Rainy mask subset-night-time.

Date	Total Pixels	Training Pixels	Validation Pixels	Testing Pixels
5 June 2015	15,826	11078	2374	2374
20 February 2016	150,653	105,457	22,598	22,598
12 October 2015	37,384	26170	5607	5607

for the day-time and night-time cases, respectively (in Table 6 the 8 July 2015 event is not reported due to night time data unavailability). As training algorithm, a Bayesian Regularization has been preferred. This algorithm typically requires more time, but provides good generalization for difficult, small, or noisy datasets like those used in this study.

As a final step, the dataset composed only of rainy pixel obtained at the previous step is used to train the neural network with GPM rain rate as output. Moreover, in this case, the subset is randomly divided in three parts, 70% for training set, 15% for the validation set, and 15% for the test set as reported in

Table 7. Rain rate estimation dataset-daytime.

Date	Total Pixels	Training Pixels	Validation Pixels	Testing Pixels
5 June 2015	30,709	23,031	4606	4606
20 February 2016	26,385	19,788	3957	3957
12 October 2015	44,855	33,641	6728	6728
8 July 2015	15,054	11,290	2258	2258

and

Table 8. Rain rate estimation dataset-night-time.

Date	Total Pixels	Training Pixels	Validation Pixels	Testing Pixels
5 June 2015	7913	5934	1186	1186
20 February 2016	75,322	56,491	11,298	11,298
12 October 2015	18,693	14,019	2803	2803

(as for Table 6, in Table 8 the 8 July 2015 event is not reported due to night time data unavailability). Actually, the NN was trained to estimate the base-e logarithm of the precipitation rate so that training would not be dominated by heavy precipitation rates at the expense of lighter precipitation rates which are more common.

Summarizing, we can say that the rain rate from SEVIRI data is retrieved in three steps:

• After the re-gridding, not-cloudy pixels are filtered out using the cloud mask.
• The first NN is applied only to cloudy pixels to obtain rainy pixels.
• The second NN is applied only to rainy pixels to retrieve rain rate.

4. Results and Performance Evaluation

To test and validate the data-fusion technique proposed in this study, precipitation fields for the case studies in Table 2 have been recreated and compared with radar data (only pixels with quality index ≥ 0.6) and PR-OBS-3A rainfall estimation product. To make these three products comparable, they must be brought to the same spatial grid. This has been done with the nearest neighbor algorithm, assigning an average value to each node that has one or more points within a radius of 3 km centered on the node. The distance between each of the nodes is 5 km. The average value is computed as a weighted mean of the nearest point from each pixel inside the search radius. The weighting function used is:

$$w\left(r\right)=\frac{1}{1+d^2}, d=\frac{3\ast r}{s\_r}$$

(1)

where s_r is the search radius and r is the distance from the node. Furthermore, MSG SEVIRI images, weather radar rain rate and PR-OBS-3A products have the same time sampling of 15 min.

In

Table 9. Number of pixels processed for performance evaluation.

Date	Images	Pixels Classified Day	Pixels Classified Night	Pixels Retrieved Day	Pixel Retrieved Night
5 June 2015	37	48,7547	658,634	57,097	49,732
20 February 2016	89	1,777,737	2,914,099	376,591	272,805
12 October 2015	37	739,551	409,922	34,193	19,412
8 July 2015	17	342,565	−	52,542	−

the number of images processed, the number of pixels classified as rainy or not rainy (both night and day) and the number of pixels for which rain rate has been retrieved (both night and day) are presented for each case study.

For all the case studies, two different statistical scores have been considered, continuous and dichotomous statistical scores. Dichotomous statistical scores are calculated all over the pixels classified as rainy or not. Continuous statistical scores are calculated all over the pixels classified as rainy and for which rain rate has been retrieved.

4.1. Dichotomous Statistical

The dichotomous statistical scores help to evaluate the accuracy of rain detection of the technique and then of the classifier. Six scores are used: probability of detection ($POD)$, false alarm rate ($FAR$), critical success index ($CSI$), the frequency bias index ($BIAS$), probability of false detection ($POFD$), and the percentage of corrects ($PC$). The rain/no rain threshold was set to 0.2 mm/h. Their definitions are:

$$POD=\frac{s1}{s1+s3}, FAR=\frac{s2}{s1+s2}, CSI=\frac{s1}{s1+s2+s3},$$

(2)

$$BIAS=\frac{s1+s2}{s1+s3}, POFD=\frac{s2}{s2+s4}, PC=\frac{s1+s4}{N}.$$

(3)

where s1, s2, s3, s4 are defined according to the contingency

Table 10. Contingency table.

	Radar
		Rain	No Rain	Total
Satellite	Rain	s1	s2
	No Rain	s3	s4
	Total	s1 + s3	s2 + s4	s1 + s2 + s3 + s4 = n

.

• POD measures the fraction of observed events that were correctly identified, the best value is one.
• FAR measures the fraction of estimated rainy pixels that were not raining, its best value is zero.
• CSI measures the fraction of estimated events that were correctly identified, its best value is one.
• BIAS measures the overestimation or underestimation of the method, then a BIAS greater than one indicates an overestimation, while a BIAS lower than one indicates an underestimation.
• POFD indicates the fraction of pixels incorrectly identified by the satellite and its optimal value is zero.
• PC is the percentage of correct estimations and its best value is one.

A pixel is considered “rainy” if the rain rate is greater than 0.2 mm/h. In

Table 11. Dichotomous scores for neural network technique (NN) and PR-OBS-3A (H03) algorithm.

Case Studies	POD NN/H03	FAR NN/H03	CSI NN/H03	BIAS NN/H03	POFD NN/H03	PC NN/H03
DAY-TIME
5 June 2015	0.79/0.45	0.34/0.65	0.69/0.38	1.03/1.51	0.07/0.09	0.91/0.34
20 February 2016	0.73/0.39	0.40/0.71	0.61/0.31	0.92/2.37	0.09/0.08	0.89/0.33
12 October 2015	0.74/0.37	0.35/0.69	0.61/0.29	0.81/2.21	0.08/0.07	0.91/0.35
8 July 2015	0.81/0.51	0.31/0.62	0.64/0.41	1.13/1.01	0.05/0.08	0.93/0.35
NIGHT-TIME
5 June 2015	0.72/0.41	0.39/0.64	0.62/0.33	1.05/1.67	0.06/0.09	0.88/0.39
20 February 2016	0.71/0.40	0.43/0.73	0.60/0.35	0.91/2.11	0.10/0.07	0.88/0.39
12 October 2015	0.69/0.33	0.40/0.73	0.59/0.42	0.79/1.99	0.09/0.07	0.87/0.38

the statistical results are listed.

According to the values in Table 11, globally, the fraction of raining cloud detected by the proposed technique is quite consistent with the one detected by weather radars. Indeed, high POD values go always together with very low POFD and low FAR. Bias values show a more complex situation. During spring/summer case studies the technique seems to overestimate raining area and the opposite in autumn/winter. As it appears clearly, the data fusion technique is always much better the PR-OBS3 product which has performance consistent with results of its validation report [31].

Performance during day time are usually better, showing that the use of visible channel has a measurable impact on the algorithm. As suggested in [46] it seems that information on CWP (provided by VIS0.6 and NIR1.6) helps to identify precipitation area better. This behavior is not present for PR-OBS3 product because it only uses 10.8 µm infrared channel during day and night.

4.2. Continuous Statistical Scores

The continuous statistical scores are useful to evaluate the accuracy of the rain rate estimation. To perform a better comparison between the NN technique and PR-OBS3 algorithm, the same PR-OBS3 continuous scores are used [30]: mean error (ME), mean absolute error (MAE), standard deviation (SD), the root mean squared error (RMSE), and root mean squared error percent (RMSE%). Their definitions are:

$$ME=\frac{1}{N}{\displaystyle \sum }\left(NN-WR\right),$$

(4)

$$MAE=\frac{1}{N}{\displaystyle \sum }\left|NN-WR\right|, SD=\sqrt{\frac{1}{N}{\displaystyle \sum }{\left(NN-WT-ME\right)}^2}$$

(5)

$$RMSE=\sqrt{\frac{1}{N}{\displaystyle \sum }{\left(NN-WR\right)}^2}, RMSE\%=\sqrt{\frac{\frac{1}{N}\sum {\left(NN-WR\right)}^2}{W{R}^2}}$$

(6)

where NN is the rain rate retrieved by the data fusion technique and WR is the rain rate from weather radar. The sum is all over the pixels of all 15 min images.

• ME provides information on the local behavior of the rain. If the errors compensate with reverse signs the score can be good even with large errors. Range: −∞ to ∞, perfect score: 0, unit (mm/h).
• MAE measures the average magnitude of satellite data respect to radar data; in this way it is a scalar measure of forecast accuracy. It must be used in conjunction with ME. Range: 0 to ∞, perfect score: 0, unit (mm/h).
• SD measures how local behavior move away from ME. Range: 0 to ∞, perfect score: 0, unit (mm/h).
• RMSE is more sensitive to large errors than the MAE. Range: 0 to ∞, perfect score: 0, unit (mm/h).
• RMSE% compensates large sensitivity of RMSE to large errors. Range: 0 to ∞, perfect score: 0, unit (×100%)

The product has been evaluated only for the precipitation class for which rain rate is greater than 0.25 mm/h. In

Table 12. Continuous scores for neural network technique (NN) and PR-OBS-3A (H03) algorithm.

Case Studies	ME NN/H03	MAE NN/H03	SD NN/H03	RMSE NN/H03	RMSE% NN/H03
Day-Time
5 June 2015	−1.23/−2.21	1.12/2.51	3.11/4.75	2.07/5.99	1.28/2.56
20 February 2016	−1.11/−1.98	1.48/2.41	2.54/4.59	2.39/3.56	1.34/1.99
12 October 2015	−1.09/−2.01	1.23/3.01	2.11/5.12	2.25/4.01	1.31/2.01
8 July 2015	−1.31/−1.89	1.57/2.22	3.09/2.09	2.48/6.13	1.43/2.43
Night-Time
5 June 2015	−1.28/−2.25	1.64/2.49	2.99/4.84	2.89/6.03	1.62/2.67
20 February 2016	−1.31/−2.03	1.23/2.38	2.77/4.77	2.37/4.83	1.33/2.54
12 October 2015	−1.19/−1.99	1.85/2.98	3.23/5.00	3.15/4.01	1.87/2.00

continuous statistical scores are reported.

The interpretation of the continuous statistical scores is less straightforward than the dichotomous ones even when a visual comparison is provided (see below). Globally, looking at ME score, it appears clear that the technique underestimates the amount of precipitation at ground, and this behavior is more pronounced during summer/spring case studies when the precipitation is higher. MAE score confirms this fact. Standard deviation is usually huge compared with ME, in particular for low regime raining. Comparing these results with those of PR-OBS3 it is evident that both products suffer of the same problems but in general the proposed technique provides a significant improved estimation of the rainfall.

A visual comparison is presented in Figure 3 that shows precipitation data on 5 June at 17:30 UTC over Belgium. The images confirm the interpretation of the statistical scores. Figure 3a,b represent rain rate of the data-fusion technique and rain rate from radar, respectively. The main path of the precipitation crossing France, Belgium, and North Sea over Holland is well recognized. The isolated storm in the low part of the image around 49° N 5° E is also fine detected. The underestimation of the model is evident as well. At peaks of precipitation around and beyond 20 mm/h in various pixel of the radar image correspond values of precipitation around 12 mm/h in the model. Values of FAR of about 0.4 are consistent for instance observing the precipitation area over the North Sea that is overestimated in the NN model. On the other hand, the model shows an output more similar to the nearest in time LEO passage of a PMW sensors (SSMIS Figure 3c). This is not unexpected because of the use of GPM data as truth to train the NN. It should also be considered that the quality of the radar data, according to Figure 4, is quite low (<0.5) over the North Sea.

In Figure 3d, the rain rate from PR-OBS3A product is shown. To better understand its behavior, the SEVIRI brightness temperature at 10.8 µm has been added in Figure 5. The main path of the precipitation is recognized, but the area is overestimated confirming the high values of FAR of Table 11. Comparing the two images it appears clear that PR-OBS3A is mostly driven by the cloud top temperature spreading high precipitation over a larger area. It should be noted that SSMIS product also tend to spread precipitation over a larger area then the radar.

Furthermore, during the study it has been noticed that the results significantly benefit choosing a higher threshold of the quality index of radar data (e.g., QI ≥ 0.8 instead of the actual value 0.6). Nonetheless, because of the re-gridding process, the number of radar pixels available for the validation decreases until not to have a valid statistical sample. This is clear in Figure 5 where typical quality index values for OPERA radar composite are shown (they do not change significantly between images through time). Many studies conclude that precipitation measurements from ground (e.g., radar and rain gauge) are subject to large uncertainty, especially when they are used to validate retrievals from satellite data and the results presented in this study confirm this fact [30,47].

5. Conclusions and Future Improvements

A new data fusion technique based on NNs was presented. It uses both visible/infrared data from GEO satellite and passive microwave data from LEO satellite to avoid their single limitations and to take advantage of both types of sensors. It provides, for a given geographical area and a time interval, a reconstruction of the surface rain field with the MSG-SEVIRI spatial and temporal resolution, e.g., 3 km at the sub satellite point and 5 km at mid-latitudes, every 15 min.

The present method differs from previous NNs based works because it presents a significantly greater number of input using all the information available from VIS/IR data, two different NNs for day and night scenario, and the possibility to use all LEO PMW data. Moreover, these preliminary results, even if over a limited number of cases not covering the huge variability of the European climate, has shown its capacity over different climatological scenarios.

The results that were obtained in the identification of rainy clouds in the various test areas, compared with the benchmark precipitation product, PR-OBS3A, show an improvement ranging from 59% up to 109% in terms of probability of rainy areas detection and of from 39% to 50% in terms of false alarm rate. The use of two different NNs for day and night has shown the impact of the visible channels in the estimation of the cloud water path. During the night, when the information at the visible bands is missing, the measure of such a parameter is more difficult and POD, FAR, and CSI have worse results, but always better than PR-OBS3A that uses a single infrared channel.

As far as the retrieval of the amount of rainfall is concerned, compared with the HSAF product, in the considered test cases our algorithm increases the accuracy from 21% up to 65%. Nevertheless, a global underestimation is present, during summer/spring case studies. This behavior seems to be common to precipitation products based on visible and infrared data; however, the quality of radar data maybe an additional issue which should be more investigated.

A possible future improvement might be the inclusion of a continuous updating and re-training of the NNs as new measurements are available for the area of interest. In this way, the technique might be more accurate and also considered for tracking or even prediction purposes although this would involve the availability of suitable processing capabilities. This would also permit to make the product operational and then to validate it in many different climatological scenes, possibly covering the European weather variability.