Analytical Reproduction of Radar Images of Rainfall
Spatial (2D) Distribution as an Input to GIS Based
Flood Modelling
B. Dodov1, R. Arsov1, B. Marinov2
1 University of Architecture, Civil Engineering and Geodesy; Faculty of Hydrotechnics;
1 Chr.Smirnensky Blvd., 1421 Sofia, Bulgaria; Tel./ Fax: +359 2 668 995;
E-mail: ;
2Bulgarian Academy of Science; Central Laboratory for Parallel Processing;
Acad. G. Bonchev Str., Bl. 25 A; 1113 Sofia, Bulgaria; E-mail:
Abstract
A new approach for radar rainfall images compression, recovering and interpolation have been suggested, based on an adaptive wavelet transform with fractal basis. It combines advantages of the wavelet transform (WT) and the adaptive wavelet approximation (AWA) approaches and offers a possibility for semi-deterministic reproduction of radar rainfall images by means of restricted number of parameters on which statistical analyses could be performed. This technology promises to be a powerful tool for reliable prediction of storms detailed spatial and temporal distribution with prescribed probability.
Keywords
adaptive wavelet transform, cascade processes, fractals, multifractals, radar rainfall images, modelling
1. Objectives
Increasingly implementation of the remote sensing technologies for atmosphere observation and availability of automated radar, satellite and airborne rainfall images in particular, is a good base for improving the accuracy of the rain fields assessment as an input to the adequate floods and pollution transfer modelling.Although the information from these observations is enough detail to support directly the sophisticated hydrological models, it is huge in quantity and is organised in a manner, which makes very difficult the extraction of whatsoever general trends with respect to process development. Regardless of the fact that this information is much more detail than the conventional one, the resolution of the most popular and accessible radar images is usually not less than 2x2km per pixel, with time interval of about 4-5 minutes.
*Published in Flood Issues in Contemporart Water Management, NATO Science Series, Environmental Security – Vol.71, Kluwer Academic Publishers,1999, The Netherlands
Such a resolution is still too “coarse” for adequate flood modelling, particularly at urban territories. It is reasonable question to be asked if there is such a technique, which could
allow extraction of the most important features of any sequence of images, describing a rainfall event with minimum information resources in a manner, that would allow subsequent restoration of the process to prescribed resolution, preserving the visual resemblance and statistical properties of the original data. The importance of this question motivated the reported research and fortunately now we can answer it positively.
In order to explain better the idea, we find it reasonable to introduce some important definitions and explanations, concerning the terms and analytical procedures, used in the paper.
There is not a unique definition of fractals, but they could be described as objects of a complex structure, consisting of parts, which are at any scale similar to the whole in some way. An important feature of the most of the fractals is the scale invariance (or scaling). A fractal function (t) is considered to be scaling with factor > 0 if the homogeneity relation (t) = (t) holds for all > 0 [4].
While the fractal is a set, the multifractal could be considered as a measure or distribution (of mass, energy, probability, feature, etc.). The Hausdorf dimensions of the sets, represented in a relevance to their singularityexponents [4], define the multifractal spectrum of the relevant density (precipitation) distribution. The singularity exponent measures the “strength” of singularity of the relevant distribution at particular position x0 in the domain under consideration. The estimated spectrum is known also as Legendre’s one, because it is calculated trough Legendre transform of the so-called mass exponents – τ(q), accounting for the scaling lows of the q-th order moments of the data. For the distribution this parameter could be defined as follows [20]:
(1a,b)
where (C) is the total mass of over the (hyper)cube C with side ;
G - the set of all the non-overlapping cubes, covering the support of .
If the distribution is scale invariant, the values of S for particular q form a straight line in a log-log plot (1a), which slope is the corresponding value of τ. Then the Legendre spectrum could be estimated through the transform:
(2a,b)
where is the singularity exponent, f() - corresponding value of the singularity spectrum.
The multifractal theory implementation with respect to the atmospheric fields modelling is founded on the assumption that fluxes of mass (water) and energy in the atmosphere could be represented by the so called multiplicative cascade processes, successively transferring these quantities from larger to smaller scales [12,13,28]. The cascade processes could be considered as eddies of fluxes, broken down into smaller sub-eddies, receiving a part of the flux of its parent eddy and the amount of which is being determined by a random factor.
Applying different strategies and “levels” of randomness, as well as some flux dissipation, an imitation of very complex structures such as rainfall patterns or clouds could be simulated [28]. The artificial fields produced by cascade processes exhibit usually good multifractal properties, which could be changed at different stages of resolution [9].Regardless of this, the fact that they are stochastically based makes the chance to fit some real data insignificant. May be themost promising deterministic technique to produce “naturally looking” fields is the so-called “Fractal/Multifractal” (FM) approach, which could be considered as a deterministic cascade process “refracted” through a Fractal Interpolation Functions (FIF) [1, 6, 14-18].
Morlet [10] introduced the wavelet transform (WT) in connection with the analyses of geophysical signals with very different scales. This transformation basically leads to decomposition of the signal in terms of some elementary functions b,a,obtained from a “mother” (basis) function by dilations and translations, such as:
b,a(x) = a-1/2b,a((x-b)/a), (3)
where a and b are real numbers.
The wavelet transform applied continuously with real values of a and b is named continuous wavelet transform (CWT) since it produces continuous function in scale-space domain (called often wavelet space). Later on, the multiresolution analysis discovered by Mallat [7, 8] and describing any orthogonal wavelet basis as common multiscale structure, was combined with the compactly supported wavelet basis [2, 3] to give a fast transform, working on discrete scale-space domain. This transform was respectively called discrete wavelet transform (DWT). The inherent scaling structure of the wavelet transform makes it a perfect tool for multifractal analyses, including continuous [11] and recently - discrete wavelet transform techniques [19, 20].
The so-called Haar scaling and wavelet functions provide the simplest example of an ortho-normal wavelet basis. The Haar scaling and waveletcoeficients could be represented as:
Sj,k = 2-1/2 (Sj+1,2k+ Sj+1,2k+1),Wj,k = 2-1/2 (Sj+1,2k Sj+1,2k+1). (4)
In the notion of the multiresolution analyses the same expression could be written as
(5)
where j+1,1j+1,2 1.
Wavelet decompositions contain considerable information about the oscillations (singularities) of the recorded signal. The greater the oscillations at particular scale are, the greater are the WT coefficients at positions at which the basis “fits” best the data. As has been observed in many contexts, the most wavelet coefficients are very small and the large ones draw the so-called “maxima” lines in the wavelet space. These features are the reasons for the image compression and detecting of the singularities (edge detection in 2D) to be among the main applications of the WT.
In his recent work Riedi [20] introduces the term asymptotically wavelet self-similar with parameter Hwave for the increments X of process Y if the following limit exist:
(6)
where Wj,k are the wavelet coefficients at the j-th level of transformation.
The parameter Hwave could be used as equivalent to the H parameter of the so-called long-range dependence (LRD), governing the correlation at large lags. This means that changing the trend of the wavelet coefficients’ decay, we could control the power low behavior at prescribed scale intervals and therefore to modell a long-range dependence observed at real data.
The values of the parameters j,1andj,2 could be optimized in a way that provides the lowest possible values of the squares of the wavelet coefficients, generated at the j-th level of transformation. Such optimization, called adaptive wavelet approximation (AWA), has a unique and relatively simple solution, which is discussed elsewhere [32].
The support of the wavelet basis could be extended to work over four or eight values of the approximated function using the so-called Rademacher -matrices [26]. This extension provides possibilities to apply the Haar WT and the adaptive wavelet approximation with respect to the surfaces and volumes. More information could be found in [21-27].
The subject of this paper is a new approach for the radar rainfall images compression, recovering and interpolation, based on the Haar WT and the AWA. Since the AWA’s basis (mother) function is generically multiplicative, it is very suitable tool for representation and analyses of multiplicative (cascade) processes, such as the rainfall distributions are expected to be because of the turbulent environment of their genesis.
2. Results and Discussion
Since the wavelet transform (WT) is a linear operation, it could not be expected its product – the wavelet transform coefficients (WTC), their absolute values in particular, to have different type of distribution then the original signal. Instead, the parameters of the distribution of the WTC will change at any successive level of the WT. As the processes under consideration are expected to be multiplicative cascades, log-normal distribution of the wavelet coefficients could be observed. This means that exploiting (6), a linear correlation could be expected in regard to the parameters of the normalized distribution of the WTC, produced at any successive level of the WT. If confirmed, this would allow adequate interpolation of data, resulted from multiplicative cascade process.
To confirm the above statement we processed 20 images, representing values of the radar reflectivity (so-called dBz-reflectivity - Z), related to the physical reflectivity - z, [mm6/m3] according to the expression 10log10z = Z. As it is shown in Fig. 1, a normal statistical distribution of the relevant data Z has been obtained. Therefore, the statistical distribution of the radar reflectivity z itself is obviously a log-normal.
In the term of the multiplicative cascades, the application of the adaptive approximation means that for a deterministic binomial process the mean and the variance of the wavelet coefficients at any successive level of transformation will be practically zero. Arising the randomness of the process, the mean and the variance of the coefficients arise up to the level, where “absence of adaptivity” is observed and measured by the corresponding Haar transform. Figures 2a and 2b represent comparison between Haar WT and AWA in regards to mean and standard deviation of the first level wavelet coefficients of binomial cascades with base 0.1 and 0.4, respectively for arising tolerance of randomness from 0 to 100 % (from 0 to 0.1 and from 0 to 0.4, respectively). Obviously the AWA insures higher compression of the signal than the conventional Haar WT . The values of the ratios j1 /(j1+j2)for cascade with base 0.4are shownin Fig. 3.
As it could be seen in the Fig. 3, irrespective of the degree of randomness, the adaptive wavelet approximation detects well the base of the underlying cascade. The base at any level can be calculated for data with length 2SN through the following expression:
bj,i = j,i/(j,1+…+j, 2s), (7)
where j = 1,…,N; S = 1,2,3; i = 1,…,2S and N represents the index of the highest level of resolution. It can be seen from Fig.3, that the higher the randomness is, the higher is the variance of the -s.
Fig. 4 represents the ratios j1 /(j1+j2) = bj1 of a binomial cascade with base 0.4 and length of 64, 128, 256, 512, 1024 and 2048, respectively with 100 % tolerance of randomness. It can be seen that the values bi1 (and respectively, the base vectors) alternatively change their directions, decaying to bN1. This corresponds to the eddies’ direction and size changing at the successive scales of the turbulence, which at the finest resolution must disappear. Therefore the AWA application allows revealing of the availability of cascade behavior of given process (signal) and could be denominated also as inverse cascade approximation (ICA).
Some of our preliminary analyses of real data indicate that the values of N,i tend to one. This could be expected, since any other value means that the energy and the mass will never dissipate in the turbulent processes. In this sense the AWA is a good tool for verification of different types of turbulent processes’ models. The first model which we tested and which exhibited good properties in this sense is the FM one.
The analyses carried out so far show that the adaptive wavelet approximation exhibits high sensitivity with respect to the existence of a “cascadeness” in the analyzed signal. This sensitivity reflects in the possibility to extract the main features of the signal within minimum amount of information, or in other words – with higher compression in comparison to the corresponding Haar transform. The rate of the compression is strictly related to the base of the underlying cascade and the “level of randomness” implanted in the signal. The closer the underlying cascade is to the Lebesque measure (equal bases) and the higher the degree of randomness is, the less is the possibility for additional compression, compared to the Haar WT. In this sense the practical aspects of application of the adaptive wavelet approximation are related mainly to the possibilities for higher degree of compression of turbulence-related data and respectively, the possibilities for minimizing of the randomness, implied in the restored (modelled) ones.
An important feature of the WT is the possibility for the data interpolation based on (6), which allows more detailed representation of the relevant process development in space and time, preserving its inherent multifractal properties. Interpolation of the data must be performed by the Haar transform or by AWA if the last three -sets are relatively equal and the trend of the WTC is captured. If the last condition is not fulfilled the AWA is not applicable, because the -s distortions the WTC distribution and (6) will not be valid.
The practical conclusions from the above notions in the context of the radar data analyses are related to the possibilities for:
- storage and analyses of the -coefficients, the core wavelet coefficients marking the general maxima lines and the trends of decay (or distributions at any level of transformation) of these coefficients;
- analyses of the obtained information in the multiparametric probability space to obtain the trend in a particular process, or to find the position of this process in the probability space;
- restoration, interpolation and modelling of radar images, based on the derived coefficients, filling in the gaps of the missing coefficients with random ones (but derived from distribution, inherent for the data) and finally - permutation of the random coefficients toprovide positive values for the recovered data;
- governing the scaling behavior of the model through control over the rates of decay of the wavelet coefficients at any level of the transform and therefore to match the bifurcation points (or zones) of the natural turbulent process. This is essential for the successful interpolation of this process in space and time.
Example of 1D interpolation of data, produced by means of the FM approach is shown in Fig. 5. The corresponding τ(q) and f() are given in Fig.6 and the power spectrum is given in Fig. 7. A similar example but in 2D is represented in Fig. 8a and Fig. 8b.
Interpolation of real data of a radar image is illustrated in Fig.9. The interpolated image is not an exact copy of the natural one in visual aspect (as it is at the examples, illustrated in Fig. 5 and Fig. 8 a, b), because of missing of detailed data of the registered radar signal. The low intensity signals are cut automatically by the radar being used (as could be seen in Fig. 1), since they are currently considered to be of a little practical importance. Taking into account as many details as possible however is of crucial importance for the adequate modelling of the multiplicative cascade processes [18]. Missing of low intensity signals data changes the pattern of “cascadnes” of the turbulent process (and respectively, the distribution of the wavelet coefficients) and does not allow an adequate interpolation. This reflects the accuracy of the modelling, as well.
3. Conclusions
The adequate assessment and modelling of the rainfall processes have been recognized as a key issue in the successive flood prediction and mitigation. The rainfall/runoff processes modelling in the current practice emphasize mainly on the developments on the ground. Even in this cases the possibilities, offered by the GIS have not been used in the advanced models in the way, which allows an accurate accounting for the topology and land use on the runoff development. Higher space and time resolution can be implemented for more adequate modelling both of rainfall and runoff processes in a GIS environment, than the currently used approaches offer.
Since many of the existing representations are limited by their analytical tractability and because there has been a recognition of chaotic effects in atmospheric precipitation, a new class of models, based on Fractal/Multifractal (FM) approach, wavelet transform (WT) and adaptive wavelet approximation (AWA) have been considered in the rainfall pattern description.
The multifractal theory implementation with respect to the atmospheric fields is founded on the assumption that fluxes of water and energy in the atmosphere could be represented by multiplicative cascade processes, successively transferring these quantities from larger to smaller scales. For this considering of the smallest details of the turbulent process is of crucial importance.
The new technology discussed in this paper is based on the advantages of the WT and AWA approaches and allows analyses, modelling and detailed interpolation of turbulence-related processes, such as rainfalls. It offers a possibility for semi- deterministic reproduction of radar rainfall images by means of restricted number of parameters and allows statistical analyses to be performed on them. This technology promises to be a powerful tool for reliable prediction of storms detailed spatial and temporal distribution with prescribed probability.
4. Recommendation
It is reasonable a catalogue to be created with optimal number of the coefficients, which for prescribed size of radar images, could allow detailed restoration of the data in statistical and visual sense. Fortunately the numerical aspects of the notion “similarity in statistical and visual sense” has it’s mathematical analogue, such as: conventional statistics and multifractal analyses for statistical similarity, as well as integratedHausdorf distance [27] for visual similarity.