Identification and filtering of rainfall and ground radar echoes using textural features

This study deals with the identification of precipitation and ground echoes in radar images using their textural features. The images were collected by meteorological radars in the regions of Setif (Algeria) and Bordeaux (France). Two kinds of texture-based techniques have been considered, consisting in calculating either the histograms of grey levels or the histograms of their sum and difference. Hence, the first-order probability distributions were found to be sufficient to account for the textural features of radar images. Energy is found to be the textural parameter that clearly distinguishes between precipitation and ground echoes. With both methods, fixed ground echoes and anaprops are efficiently rejected, whereas precipitation echoes are kept almost unchanged. These methods have the advantages of effectiveness and simplicity. The threshold of discrimination is independent of the geographical and climatic conditions in the regions under study. Because the computation time needed by these techniques is small, the radar images can be processed in real-time.


Introduction
Meteorological radars are today acknowledged as efficient tools for detecting precipitation fields in a given region. These radars collect images made of echoes mostly caused by rainfall. The rainfall echoes vary in intensity and surface. A serious problem in meteorological radar systems is the intrusion of permanent and sporadic ground clutter echoes. These echoes are undesirable since they considerably reduce the radar performance and lead to non-negligible errors in the evaluation of precipitation and identification of storms. The ground echoes observed around the radar are caused by the reliefs or the neighbouring constructions. Although, they exist permanently at fixed positions and their coordinates are well known, they cannot easily be removed since they fluctuate from pulse to pulse (Ulaby and Dobson 1994). When super-refraction settles down and when the air refractive index strongly varies with altitude, electromagnetic waves arising from the radar propagate through atmospheric ducts and are backscattered by the ground obstacles. In such conditions, another type of disturbing ground echoes, called anaprops or AP echoes (anomalous propagation), sometimes appears at various distances from the radar. Anaprops depend on the state of the lower atmosphere. They are sporadic, vary slowly in position and cause many false observations of precipitation. Their suppression is then uneasy.
Computer animation can be used to distinguish between rainfall echoes, anaprops, and fixed-short-distance ground echoes in radar images. To perform AP rejection, several approaches were proposed, for example real-time analysis of the autocorrelation of the radar signal (Tatehira and Shimizu 1978, Aoyagi 1983, Hamuzu and Wakabayashi 1991, climatology of anaprops (Smith 1984), comparison of radar to satellite images (Fiore et al. 1986), neural networks (Pamment andConway 1998, Pankiewicz 1998), three-dimensional reflectivity structure (Steiner and Smith 2002). In the case of coherent or polarimetric radars, ground clutters can be eliminated by cancellation techniques such as Doppler signal processing (Doviak andZrnic 1993, Seltmann andReidl 1999), or polarization filtering (Giuli et al. 1991, Ryzhkov and Zrnic 1998, da Silveira and Holt 1999. Many radars equipping the meteorological stations are neither coherent nor polarimetric. For such radars, ground clutter removal can be performed by discriminating between the statistical properties of rainfall and ground echoes. In particular, efficient non-coherent filtering approaches can be implemented if Markov processes (Haddad et al. 2000), neural networks (Grecu and Krajewski 2000), or texture analyses are used. For example, Moscowicz et al. (1994) used the Bayesian approach to get a textural classification of the radar echoes. Gabella et al. (1999) filtered the undesirable radar echoes by removing the pixels which are spatially weakly correlated with surrounding ones. Texture is the word used to characterize the spatial distribution of the elements of a given surface in an image. In the literature, it has also been interpreted as a set of measurements of the spatial statistical distributions of grey level (Davis et al. 1979, Haralick 1979, Sun and Wee 1983, Harlow et al. 1986, Unser 1986. For a given image, textural features are then defined by the processing of these measurements. These features are, for example, fineness, granulation, smoothness, and randomness. Texture-based methods are extensively used to process satellite images. The main models used to account for the textural properties are those based on the power spectrum representation (Pratt 1976), the autoregressive processes (Haralick et al. 1973, Haralick, 1979, Huang and Gong 1988, the fractals (Mandelbrot 1982, Lovejoy andSchertzer 1986), and the mosaic approach (Ahuya and Rosenfeld 1981).
The goal of the present paper is to emphasize the ability of the textural approaches to accurately segment and classify the objects composing radar images. It is then interesting to apply such methods to the echo type identification. To get the most important textural features of radar images and use them to remove the ground clutter from precipitation echoes, we have considered two sorts of statistical approach. The first one consists of the calculation of the histograms of grey levels and the second one is based on the sum and difference of histograms of grey levels. Both methods are applied, hereafter, to filter precipitation echoes in radar images collected in the regions of Setif (Algeria) and Bordeaux (France).

Data
The number of processed images is 10 052 for Bordeaux and 5799 for Setif. These images were collected with non-coherent pulsed radars used for meteorological observations. The technical specifications of these radars are reported in table 1. All the images are radar reflectivity fields plotted following the PPI (Plan Position Indicator) mode. They consist of 512 pixels6512 pixels and cover an area of 512 km6512 km. The number of levels encoding each pixel of these images is 52 for Bordeaux and 16 for Setif (due to technical choices of the radar system manufacturers).
The radar set up near the town of Setif, at a latitude of 36 ‡11'N, a longitude of 5 ‡25'E W and an altitude of 1700 m, is one of the seven radars of the Algerian Meteorological Network. The Setif region is part of the High Tablelands of Algeria ( figure 1(a)). Several ground obstacles are distributed around the radar and the nearest ground echoes mainly arise from the constructions of the town of Setif. Beyond the horizon, the ground obstacles cause a lot of ground echoes in radar images. For example, in the southeast, 60 km away from the radar, rise the mountains of Djurdjura whose height reaches an altitude of 2300 m. In the southwest, 40 km away from the radar, are the mounts of Bibans with a highest point at an altitude of 1417 m. In the north-east, at shorter distance (about 30 km from the radar), the mounts of Babors reach 2004 m of altitude. The climate of Setif is continental. In winter, most of the atmospheric perturbations come from the north-west. In summer, masses of clouds usually move from south-west to northeast. The yearly mean cumulative rainfall height is about 400 mm falling mainly in winter. Hence, the images under consideration are those collected every fifteen minutes by the radar of Setif from January to February 2000.
The radar of Bordeaux, located at the Mérignac airport, at 44 ‡52'N of latitude, 0 ‡30'W of longitude and 70 m of altitude, is an element of French operational radar network managed by Météo-France (figure 2(a)). The region of Bordeaux is rather flat. The first obstacles are the hills of the western part of Massif Central with an altitude inferior to 400 m situated at about 100 km from the radar in the northeastern direction. Farther down, in the south, at a distance of 200 km from the radar, the Pyrenees Range rises up to 3000 m altitude. The climate of the Bordeaux area is oceanic. The climatic conditions are rather homogeneous. Most rainy events observed in this region are zonal frontal cyclonic disturbances, moving eastward from sea to land. However, the frontal systems are often reduced to a convective line in summer. The cumulative annual rainfall is about 800 mm. Rainfall is more intense in summer than in winter. This justifies the present choice to study the radar images collected from June to September 1996. These images were stored every five minutes following the CAPPI mode (Constant Altitude Plan Position Indicator), with a beam elevation angle of 1.5 ‡ for distances inferior to 50 km and 0.4 ‡ above 50 km.
On the radar images used to test the methods presented in this paper, the identification of permanent ground clutter, anaprops and precipitation echoes was performed by hand (i.e. by an operator) using visual criteria associated with the change of the elevation angle of the radar beam, the comparison with a geographical map and the animation of radar images. Such an identification by a trained radar operator is very easy and not ambiguous (Sauvageot 1992).
3. Evaluation of textural features of radar images 3.1. Statistical approaches In general, the texture of an image is properly described by the statistical features deduced from the first-and second-order probability distribution functions of grey levels. The methods most commonly used for texture analysis are based on the evaluation of the joint probability that two pixels may be separated by a distance (d ) making an angle (Q) with the horizontal direction. Therefore, the set of joint probabilities obtained for all the pairs of pixels forms a Grey Level Cooccurrence Matrix (GLCM). The GLCM-based method amounts to compute the textural features of an image by assessing a set of co-occurrences matrices for different directions and orientation angles. With the GLCM algorithm, texture is completely represented by fourteen typical statistical parameters (Maisel 1971, Haralick 1979, Welch et al. 1988, Joe 1991. However, the resulting computation is complicated since each of these parameters is obtained by performing a double summation over the pairs of pixels. Consequently, the implementation of such an algorithm needs a large computation memory. To both reduce the memory requirements and reach the same performance as for GLCM, either grey level histograms or histograms of pairs of pixels can be used to estimate the first and second-order probability distribution functions (Pratt 1976, Unser 1986). Such an estimation amounts, in the first step, to build these histograms and, in the second step, to get the textural parameters which discriminate between radar echoes. Hereafter, the histogram-based approaches have been applied to two kinds of radar images, denoted by A and B, consisting of either only ground or only precipitation echoes respectively. The images of ground echoes have been collected in clear sky conditions. The images of precipitation echoes have been obtained by selecting all the radar images where rainfall echoes are distinctly separated from ground echoes, and then by masking the latter using a retouching image software. Texture is analysed by scanning the radar images with a 5 pixel65 pixel window. The frequency distributions of textural parameters are calculated for ground and for precipitation echoes. Both kinds of frequency distributions are then compared. For typical textural features, these distributions are found to be completely distinct or to slightly overlap. Hence, such a property can be usefully used to separate precipitation from ground echoes.

Sum and difference histograms
Two pixels of a radar image lying in a textural region (D), at the respective coordinates (m, n) and (m', n') (with m, n, m' or n'~1, … , 512) are considered. (i) and ( j) are the grey level of each pixel, respectively (with i or j~0, … , L-1). The number L of grey levels is 16 for Setif and 52 for Bordeaux. The domain (D) is a window of 5 pixels65 pixels centred on the pixel to be analysed. The coordinates of the displacement (d ) from one pixel to another are (Dx, Dy). The brightness of both pixels can be described by random variables, respectively, denoted y m,n and y mzDx, nzDy . The probability of getting grey levels (i) and ( j) at the relative position (Dx, Dy) is: The set of P(i, j) probabilities forms a co-occurrence matrix from which the textural features can be evaluated. Next, two new random variables expressing the sum and difference of grey levels, are defined, respectively, as: The normalized sum histograms are given by: In these expressions, # denotes the number of elements in the (D) region under consideration.
The relationships (4) and (5) indicate that the relative frequencies arising from the sum and difference histograms are almost equal to P s (i) and P d ( j), the probabilities of having the (i) and ( j) grey levels, respectively.
The textural features used in this study are (Chen et al. 1989): Variance: v~s 2~1 2 Energy: Contrast: Correlation: Entropy: Local homogeneity: C~X Skewness: The textural features of the radar images of Setif and Bordeaux have been calculated by using the above definitions, choosing the distance between pixels so that d~1 and considering four of the principal orientations Q~0 ‡, 45 ‡, 90 ‡ and 135 ‡ (with dw1, the resulting improvements are insignificant). Histograms of each textural parameter are then obtained for ground and precipitation echoes, respectively. All of them overlap in pairs, except energy and local homogeneity histograms. Figure 3 illustrates such histograms obtained for d~1 and Q~0 ‡, when processing the images A and B of Setif associated with figure 1(a) and defined in §3.1. In the case of energy and local homogeneity, it indicates that the frequency distributions characterizing rainfall echoes are distinct from those representing ground echoes. The radar echoes are then identified by applying the following conditions: Energy: ground echoes if Wv0.05 and precipitation echoes if Wo0.05. Local homogeneity: ground echoes if Cv0.6 and precipitation echoes if Co0.6.
These conditions have been verified for all the images in the database. They are then valid for the two regions under study. If energy is taken as the discriminating parameter, the ground echoes will be removed by considering all the pixels for which Wv0.05 and by replacing their colour by white. If Wo0.05, the pixel colour will obviously remain unchanged. The procedure is similar when local homogeneity is used. After such a filtering, residual dots due to the ground reflection sometimes remain in the radar images. They are then suppressed by median filtering. To illustrate this approach, a radar image where rainfall echoes are distinctly separated from ground echoes is chosen. This image denoted by C, is given in figure 1(a). It was collected in the region of Setif on 13 January 2000 at 23:20 UTC. Figure 1(b) represents the image filtered by using energy. An almost identical result is obtained using local homogeneity as the discriminating parameter. Figure 1(b) shows that the ground echoes are efficiently suppressed and that precipitation echoes remain unchanged. To characterize the texture of the image C, the mean value of the textural parameters running from equations (6)-(14) has been calculated. The results obtained are reported in table 2. Energy and local homogeneity are found to be much greater for rainfall echoes than for ground echoes and prove that the texture of rainfall echoes is more uniform. In this table, contrast, skewness, and kurtosis are different from a direction to another. The values obtained show that the texture anisotropy is more important for ground echoes. However, the other textural parameters are nearly independent from the direction. This means that the texture of the radar images under study is sufficiently described by the first-order probability distribution and that the calculation of the parameters of the secondorder probability distribution is then unnecessary.

Histogram approach
According to the results reported in §3.2, the image processing has been limited to the calculation of the textural parameters arising from the first-order probability distribution of grey levels. Let (i) be the grey level of each pixel of the radar image (with i~0, … , L-1), n i , the number of pixels at grey level (i) and N T , the total number of pixels in the image. Its relative frequency is given by: The textural features are then defined as: Energy : Contrast: Histograms have been constructed for each of the textural features of the radar images of Bordeaux and Setif. The histograms obtained by processing the ground echoes are found to be confused with those characterizing the precipitation echoes, except for the energy histograms. The conclusion here is that energy is the best discriminating parameter of radar echoes. For example, figure 4 exhibits the histograms of mean and energy characterizing the ground and precipitation echoes observed in the radar image of figure 1(a). Nothing is distinguished in the mean histogram whereas the energy histogram consists of two different parts. The first one delimited by Wv0.35 is associated to ground echoes. The second is defined by Wo0.35 and corresponds to precipitation echoes. As in §3.2, ground echoes can be suppressed in radar images if energy is chosen as the discriminating parameter. White is assigned to the pixels for which Wv0.35, the colour of the other pixels remains the same and the residual dots are eliminated by median filtering. Such a processing of our database results in filtered images where almost all the ground echoes are suppressed whereas the precipitation echoes are preserved. The resulting filtered image is found to be very similar to figure 1(b) (that is why it is not presented). The same result has been obtained for every image of the present database. According to the first-order approximation, the approaches, respectively, based on the sum and difference histogram of pairs of pixels and histogram of grey levels enable thus similar filtering performance.

Discussion
As shown in §3, the texture-based techniques eliminate ground echoes with a very good efficiency and keep the precipitation echoes almost intact. To evaluate the performances of these techniques, four data samples have been considered. The first one gathers the pixels identified by histogram and sum and difference histogram approaches as precipitation echoes. The second one consists of the pixels classified by both methods as ground echoes. The pixels sorted by histogram approach as precipitation and ground echoes are assigned to the third data sample and those classified by sum and difference histogram approach form the fourth data sample. For the sake of clarity in the computation, the number of pixels of echoes identified by histogram and sum and difference histogram approaches are, respectively, denoted X a and Y a (where a~0 stands for precipitation echoes and a~1 for ground echoes). The following correlation coefficients are then calculated: where cov denotes the covariance of the couple of random variables and s the standard deviation of each variable. Table 3 gives the values of these correlation coefficients computed for the images of Setif and Bordeaux. It exhibits a strong correlation between the same type of echoes and a quasi independence between precipitation and ground echoes. It also shows that the sum and difference histograms of pairs of pixels and histograms of grey levels yield almost the same results. Consequently, only the latter approach is considered hereafter.
Figures 2 (a) and (b) give original and filtered radar images, respectively, collected in the region of Bordeaux in clear sky conditions on 20 October 1996 at 08:15 UTC. The ground echoes appearing in the original image, are mainly caused by an anomalous propagation. These anaprops are transmitted through a 400 m thick surface duct. On 20 October 1996, nocturne radiative cooling persisted in the region under study for about twenty hours and produced this atmospheric duct. On figure 2(b) (filtered), all the ground echoes have been eliminated, except residual echoes which occupy 7% of the primary surface of ground echoes. Figure 5(a) shows a radar image of Bordeaux on 21 November 1995 at 23:45 UTC with both precipitation and ground echoes. When this image is filtered, the precipitation echoes are almost totally preserved ( figure 5(b)).
The radar image of figure 6(a) from the region of Setif on 12 January 2000 at 12:30 UTC displays the addition of a great number of ground echoes and a field of precipitating clouds coming from the Mediterranean Sea and moving from the north-west to the south-east. Near the radar, the ground echoes mainly arise from the industrial zone of Setif. Beyond the sky-line they are probably due to the anomalous propagation. After filtering, the image obtained is that of figure 6(b) To get an estimate of the performance of this method, a recognition ratio has been calculated for each type of echo from the images where rainfall echoes are distinctly separated from ground echoes. For a given type of echo, this parameter is defined as the ratio of the number of pixels obtained after filtering to that computed in the original image. We get that, thanks to the histogram approach, 93% of ground echoes are removed and 95% of precipitation echoes are preserved in the radar images of Bordeaux. In the radar images of Setif, only 90% of the ground echoes are eliminated and 93% of the precipitation echoes are faithfully reproduced. During the above-described processing, the reflectivity of each pixel of the radar images still looks the same. The degradation of the performance of the histogram approach observed for the images of Setif is probably due to the smallness of the number of grey levels of these images. When the radar images are processed by means of the histogram approach, the time of calculation is inferior to 40 s per image. The use of the sum and difference histogram method lasts about 30 s for each radar image. These results mean that both approaches enable us to identify the radar echoes in real time.
Results similar to the ones of Bordeaux and Setif have been obtained on images of other sites, notably Dakar (Senegal) in Africa.

Conclusions
The two texture-based approaches presented in this paper, namely the sum and difference histograms of pairs of pixels and the histograms of grey levels, enable the identification without ambiguity of the type of echo in the images collected by meteorological radars. Applied to the filtering of such images, they are able to remove more than 90% of the ground echoes and to preserve more than 93% of the precipitation echoes. Rather than the usual statistical parameters such as variance or correlation, a criterion consisting in comparing the frequency distributions of their textural features has been chosen to distinguish between both types of echoes. Energy was found to be the best discriminating parameter. The combination of the energy with the other textural parameters does not improve the ability of the histogram-based approaches to identify radar echoes. The histogram approach has been implemented by devising an appropriate algorithm playing the role of a realtime automatic system of radar echo filtering. The radar images under study have been collected every five minutes or more, whereas the computation time needed by this algorithm does not exceed forty seconds per image. Such a processing could be usefully employed to get a good estimate of the precipitation intensity in the regions of Bordeaux and Setif. The discriminating conditions between ground and precipitation echoes are independent from the different climates prevailing in the regions under study, which are oceanic, continental, and Sahelian respectively. These results show that the texture-based approaches associated to the abovementioned discriminating criterion improve the performance of meteorological radars for precipitation field estimates.